Cellular Automata in Matlab.docx

CellularAutomatainMatlabIntroductionCellularAutomata(CA)areaschemeforcomputingusinglocalrulesandlocalcommunication.TypicallyaCAisdefinedonagrid,witheachpointonthegridrepresentingacellwithafinitenumberofstates.Atransitionruleisappliedtoeachcellsimultaneously.Typicaltransitionrulesdependonthestateofthecellandits(4or8)nearestneighbors,althoughotherneighborhoodsareused.CAshaveapplicationsinparallelcomputingresearch,physicalsimulations,andbiologicalsimulations.ThispagewillconsiderhowtowriteefficientmatlabcodetoimplementCAsandlookatsomeinterestingrules.codemakesuseofMatlab,sveryflexibleindexingtospecifythenearestneighbors. x=2:n-l; y-2:n-l; sum(x,y)=cells(x,y-l)+cells(x,y+l)+. cells(x-l,y)+cells(x+l,y)+. cells(x-l,y-l)+cells(-l,y+l)+ cells(x+l,y-l)+cells(x+l,y+l); cells=(sum=3)(sum=2&cells); AddingasimpleGUIiseasy.Inthisexamplethreebuttonsandatextfieldwereimplmented.ThethreebuttonsallowausertoRun,Stop,orQuit.Thetextfielddisplaysthenumberofsimulationstepsexecuted. %buildtheGUI %definetheplotbutton plotbutton=uicontrol(,style,pushbutton,. ,string,Run, ,fontsize,12, ,position,100,400,50,20,. ,callback,run=l;,); %definethestopbutton erasebutton-uicontrol('style,pushbutton1,. ,string,Stop, ,fontsize,12, ,position,200,400,50,20,. ,callback,freeze=l;,); %definetheQuitbutton quitbutton=uicontrol('style','pushbutton,. ,string,Quit, ,fontsize,12, ,position,300,400,50,20, ,callback*,stop=l;close/); number=uicontrol(,style,text, ,string*,1, ,fontsize,12, ,position,20,400,50,20);Afterthecontrols(andCA)areinitialized,theprogramdropsintoaloopwhichteststhestateofthevariableswhicharesetinthecallbackfunctionsofeachbutton.Forthemoment,justlookatthenestedstructureofthewhileloopandifstatements.TheprogramloopsuntiltheQuitbuttonispushed.Theothertwobuttonscauseanifstatementtoexecutewhenpushed.stop=0;%waitforaquitbuttonpushrun=0;%waitforadrawfreeze=0;%waitforafreezewhile(StoP=O)if(run=l)%nearestneighborsumsum(x,y)二cells(x,y-l)+cells(x,y+l)+.cells(-l,y)+cells(x+l,y)+.cells(-l,y-l)+cells(-l,y+l)+.cells(3:n,y-l)+cells(x+l,y+l);%TheCArulecells=(SUnI=3)(SUm=2&cells);%drawthenewimageset(imh,cdata,cat(3,cells,z,z)%updatethestepnumberdiaplaystepnumber=1÷str2num(get(number,string,);set(number,string,num2str(stepnumber)endif(freezel)run=0;freeze=0;enddrawnow%needthisintheloopforcontrolstoworkendExamples1. Conway'Slife.Theruleis:oSumthe8nearestneighborsoIfthesum=2thenthestatedoesnotchangeoIfthesum=3thenthestate=loOtherwisethestate=OThecode:x=2:n-l;y-2:n-l;%nearestneighborsumsum(x,y)=cells(x,y-l)+cells(x,y+l)+.cells(x-l,y)+cells(x+l,y)+.cells(x-l,y-l)+cells(x-l,y+l)+.cells(3:n,y-l)+cells(x+l,y+l);%TheCArulecells=(sum=3)I(sum=2&cells);2. SUrfaCeTenSionTheruleis:oSumthe8nearestneighborsandthecellitselfoIfthesum<4orsum=5thenthestate=OoOtherwisethestate-1Thecode:x=2:n-l;y=2:n-l;sum(x,y)=cells(x,y-l)+cells(x,y+l)+.cells(-l,y)+cells(x+l,y)+.cells(-l,y-l)+cells(-l,y+l)+cells(3:n,y-l)+cells(x÷l,y+l)+.cells(x,y);%TheCArulecells=(sum<4)(sum=5);3. PercolationCIUSterTherule:oSumthe8nearestneighborsofeachcell(cellsarebinaryvalued,0/1).Cellsalsohaveaseparatestatevariable(called,visited,)whichrecordsiftheyhaveeverhadanonzeroneighborbefore.oComputearandomnumberrbetween0and1.oIfthesum>0(atleastoneneighbor)andther>threshold,andthecellhasneverhadaneighborthencell=l.oIfthesum>0setthe"visited"flagforthecel1torecordthatthecellhasanonzeroneighbor.Theupdatecode:sum(2:a-l,2:b-l)=cells(2:a-l,1:b-2)+cells(2:a-l,3:b)+.cells(1:a-2,2:bT)+cells(3:a,2:b-l)+cells(1:a-2,1:b-2)+cells(1:a-2,3:b)+cells(3:a,1:b-2)+cells(3:a,3:b);pick=rand(a,b);cells=cells(sum>=l)&(pick>=threshold)&(visit=0)visit-(sum>=l);Thevariablesaandbarethesizeoftheimage.Theinitalimageisdeterminedbygraphicsoperations.Thefollowingstatementssetupaxesofafixedsize,writetextintotheaxes,thengrabthecontentsoftheaxesandputthembackintoanarrayusingthegetframefunction.ax-axes(,units,pixels,position,11500400,color,k,);text(,units,'pixels,'position,50,255,0,string,BioNB,color,w,fontname,helvetica,fontsize,100)text(,units,pixels,position,120,120,0,.,string,44,color,w,fontname,helvetica,fontsize,100)initial-getframe(gca);a,b,c-size(initial.cdata);z=zeros(a,b);cells=double(initial,cdata(:,:,1)255);visit=Z;sum=z;Afterafewdozentimesteps(startingwiththeBioNB441image)wegetthefollowingimage.Clickonittoseeafullsizeimage.4. EXCitabIeInedia(BZreactionorheart)Therule:oCellscanbein10differentstates.State0isresting.States1-5areactive,andstates6-9arerefractory.oCountthe8nearestneighborsofeachcellwhichareinoneoftheactivestates.oIfthesumisgreaterorequalto3(atleastthreeactiveneighbors)thencell=l.oStates1to9occurstepwisewithnomoreinput.Ifstate=lthenthenextstate=2.Ifstate=2thenthenextstate=3,andsimilarlyofallthestatesupto9.Ifstate=9thenthenextstate=0andthecellisbackatrest.Theupdatecode:sum(x,y)=(cells(x,y-l)>0)(cells(x,y-l)<t)÷(cells(x,y÷l)>0)&(cells(x,y+l)<t)+.(cells(-l,y)>0)(cells(-l,y)<t)+(cells(x+l,y)>0)&(cells(x+l,y)<t)+.(cells(-l,y-l)>0)&(cells(-l,y-l)<t)+(cells(-l,y+l)>0)(cells(-l,y+l)<t)+.(cells(x+l,y-l)>0)&(cells(x+l,y-l)<t)+(cells(x+l,y+l)>0)&(cells(x+l,y+l)<t);cells-(cells-0)&(sum>=tl)+.2*(CeIIS=I)+.3*(CelIS=2)+4*(CeIIS=3)+.5*(CeIIS=4)+.6*(cells=5)+.7*(CeIIS=6)+.8*(CeIIS=7)+.9*(CeIIS=8)+.0*(CelIS=9);AnimagefromthisCAafterspiralwavesdevelopisshownbelow.5. FOreStFireTherule:oCellscanbein3differentstates.State=Oisempty,state=lisburningandstate=2isforest.oIfoneormoreofthe4neighborsifacellisburninganditisforest(state-2)thenthenewstateisburning(state-1).oThereisalowprobablity(say0.000005)ofaforestcell(state=2)startingtoburnonitsown(fromlightning).ocellwhichisburning(state=l)becomesempty(state=0).oThereisalowprobabi1ity(say,0.01)ofanemptycellbecomingforesttosimulategrowth.oThearrayisconsideredtobetoroidlyconnected,sothatfirewhichburnstoleftsidewillstartfiresontheright.Thetopandbottomaresimilarlyconnected.Theupdatecode:sum=(veg(l:n,n1:n-l)=l)+(veg(l:n,2:n1)=1)+(veg(nl:n-l,l:n)=l)+(veg(2:n1,l:n)=l);veg=.2*(veg=2)-(veg=2)&(sum>0(rand(n,n)<Plightning)+.2*(veg=0)&rand(n,n)<Pgrowth);Notethatthetoroidalconnectionisimplementedbytheorderingofsubscripts.6. GaSdynamicsThisCA(andthenexttwo)areusedtocomputemotionofparticles.Thisapplicationrequiresadifferentkindofneighborhood.Theneighborhoodofacellvarysoneachtimestepsothatmotioninacertaindirectionwillcontinueinthesamedirection.Inotherwords,theruleconservesmomentum,whichisthebasisofdynamicalcalculations.TheneighborhoodusediscalledaMargolisneighborhoodandconsistsofoverlapping2x2blocksofcells.Inthefollowingdiagram,theupper-left4cellsaretheneighborhoodduringaneventime-step,thebottom-right4cellsaretheneighborhoodduringanoddtime-step.Agivencellhas3neighborsateachtime-step,butthespecificcellswhichconstitutetheneighborsflipsbackandforth.evenevenevencelloddIIOddIoddI7. Therule:oThisiscalledaHPP-GASrule.Seereference1Chapter12.oCellshave2states.State=Oisempty,state=lrepresentsaparticleinmotion.oOnanystep,aparticleisassumedtohavejustenteredthe2x2blockonadiagonalpath.Itmustthereforecrosstheblockthroughitscenter.Soonanytime-step,swapthecontentsofeachcellwiththatofthecelldiagonalIyacrossfromit.Pairsofblocksareshownbelow.Theleftoneisbeforeandtherightafteratime-step.Sixcasesareshown,butallrotationalversionsofeachofthemaretreatedthesameway.Twospecialcases,partide-partic1ecol1isionandparticle-wallareconsideredbelow.0-goes-00tooo10goesoo00too110goes10011o011-goes011oto0111goes0110to1111goes11111o11oTomakeaparticlecollision(whichconservesmomentumandenergy),treatthecaseofexactlytwoparticlesonadiagonalasiftheyhitanddeflectedeachother90degrees.Thisisdonebyconvertingonediagonaltotheotheronthetime-step.Youcanimplementthisbyrotatingthe4cellscounterclockwisebyonecell.Thethirdruleabovebecomes:10goes0101to10oTomakeaparticlecollidewithawall,simplyleaveitsstateunchanged.Thiscausesareflection.Theupdatecode:p=mod(i,2);%margolisneighborhood,whereiisthetimestep%upperleftcellupdatexind=l+p:2:nx-2+p;yind=l+p:2:ny-2+p;%Seeifexactlyonediagonalisones%only(atmost)oneofthefollowingcanbetrue!diagl(xind,yind)=(sand(xind,yind)-l)&(sand(xind+l,yind+l)-l)&.(sand(xind+l,yind)-0)&(sand(xind,yind÷l)-0);diag2(xind,yind)=(sand(xind+l,yind)=-l)&(sand(xind,yind+l)=l)&.(sand(xind,yind)-0)&(sand(xind+l,yind+l)-0);%Thediagonalsbothnotoccupiedbytwoparticlesandl2(xind,yind)=(diagl(xind,yind)=0)&(diag2(xind,yind)=0);%0nediagonalisoccupiedbytwoparticlesorl2(xind,yind)二diagl(xind,yind)jdiag2(xind,yind);%foreverygasparticleseeifitneartheboundarysums(xind,yind)-gnd(xind,yind)gnd(xind+l,yind)gnd(xind,yind+l)Ignd(xind+l,yind+l);%celllayout:%x,yx+l,y%x,y÷lx+l,y+1%If(nowalls)and(diagonalsarebothnotoccupied)%thenthereisnocollision,somoveoppositecelltocurrentcell%If(nowalls)and(onlyonediagonalisoccupied)%thenthereisacollisionsomoveccwcelltothecurrentcell%If(awall)cell(causesareflection)%thendon,tchangethesandNew(xind,yind)=(andl2(xind,yind)(orl2(xind,yind)(sums(xind,yind)sandNew(xind+1,yind)(andl2(xind,yind)(orl2(xind,yind)(sums(xind,yind)&sums(xind,yind)&sand(xind+l,yind+l)+sums(xind,yind)&sand(xind,yind+l)+sand(xind,yind);&sums(xind,yind)&sand(xind,yind+l)+sums(xind,yind)&sand(xind,yind)+sand(xind+l,yind);sandNew(xind,yind+l)=(andl2(xind,yind)&sums(xind,yind)&sand(xind+l,yind)+(orl2(xind,yind)&sums(xind,yind)&sand(xind+l,yind+l)+(sums(xind,yind)&sand(xind,yind+l);sandNew(xind+1,yind+l)=(andl2(xind,yind)&sums(xind,yind)&sand(xind,yind)+(orl2(xind,yind)&sums(xind,yind)&sand(xind+l,yind)+(sums(xind,yind)&sand(xind+l,yind+l);sand=sandNew;8. DiffUSiOnIimitedaggregationThissystemsimulatesstickyparticlesaggregatingtoformafractalstructure,particlemotiontakesplacewitharulesimiIartotheHPP-GASruleinexample6.Thedifferenceisthatparticlesareassumedtobebouncingaroundinsomedense(butinvisible)liquid.Theeffectistorandomizethedirectionofmotionofeveryparticleateverytimestep.Putdifferently,everytime-stepisacolIisionstep.Thesimulationisalsoseededwithonefixedparticleinthecenterofthearray.Anydiffusingparticlewhichtouchesitstickstoit,anditselfbecomesanon-moving,StiCkyparticle.Therule:oUseaMargolusneighborhood.Ateverytimestep,rotatethe4cellseitherclockwiseorcounterclockwisebyonecellwithequalprobability.Therotationrandomizesthevelocities.oAfterthemove,ifoneormoreoftheeightnearestneighboorsisafixed,stickyparticle,thenfreezetheparticleandmakeitsticky.Theupdatecode:p=mod(i,2);%margolisneighborhood%upperleftcellupdatexind=l+p:2:nx-2+p;yind=l+p:2:ny-2+p;%randomvelocitychoicevary=rand(nx,ny)<.5;varyl-l-vary;%diffusionrule-margolusneighborhood%rotatethe4cellstorandomizevelocitysandNew(xind,yind)=vary(xind,yind).*sand(xind+l,yind)+%cwvaryl(xind,yind).*sand(xind,yind÷l);%ccwsandNew(xind+1,yind)=vary(xind,yind).*sand(xind+l,yind+l)+varyl(xind,yind).*sand(xind,yind);sandNew(xind,yind+l)=vary(xind,yind).*sand(xind,yind)+varyl(xind,yind).*sand(xind+l,yind+l);sandNew(xind+1,yind+l)=vary(xind,yind).*sand(xind,yind+l)+varyl(xind,yind).*sand(xind+l,yind);sand=sandNew;%foreverysandgrainseeifitnearthefixed,stickyclustersum(2:nx-l,2:ny-l)=gnd(2:nx-l,1:ny-2)+gnd(2:nx-l,3:ny)+.gnd(l:nx-2,2:ny-l)+gnd(3:nx,2:ny-l)+.gnd(l:nx-2,1:ny-2)+gnd(l:nx-2,3:ny)+.gnd(3:nx,1:ny-2)+gnd(3:nx,3:ny);%addtotheclustergnd=(sum>0)&(sand=l)gnd;%andeliminatethemovingparticlesand(find(gnd=l)=0;Thefollowingimageshowsthefixedclusteraftermanytimesteps.9. SandDiIeThecross-sectionofapileofsandcanbemodeledusingaMargolusneighborhoodtopropagatecells,butwithadifferentmotionruleTherule:oSeereference2Chapter2.2.6.oCellshave2states.State=Oisempty,state=lrepresentsagrainofsand.oOnanystep,aparticlecanfalltowardthebottomofthethe2x2block.Thepossibletransitionsareshownbelow.Wallsandfloorsstopmotion.oTomakethemotionslightlyrandom,Iaddedarulewhichsometimesreversesthestateofthetwolowercellsafterallthemovesarecomplete.00goes0000to001-goes-00to1o01goes0o001。0110goes0010to111-goes-0t（）1101goes0o101。1101goes0001to1111goes1-10to1111goes0101to11oTheupdatecode:op-mod(i,2);%margolisneighborhoodosand(nx2,ny2)=1;%addagrainatthetopo%upperleftcellupdateoxind=l+p:2:nx-2+p;oyind=l+p:2:ny-2+p;o%randomizetheflow-10%ofthetimeovary=rand(nx,ny)<.9;ovaryl=1-vary;oosandNew(xind,yind)二ognd(xind,yind).*sand(xind,yind)+.o(l-gnd(xind,yind).*sand(xind,yind).*sand(xind,yind+l).o(sand(xind+l,yind÷l)+(l-sand(xind÷l,yind+l).*sand(xind+l,yind);oosandNew(xind+1,yind)二ognd(xind+l,yind).*sand(xind+l,yind)+o(l-gnd(xind÷l,yind).*sand(xind+l,yind).*sand(xind÷l,yind÷l).*.o(sand(xind,yind÷l)+(l-sand(xind,yind+l).*sand(xind,yind);osandNew(xind,yind+l)=osand(xind,yind+l)+o(l-sand(xind,yind+l).*o(sand(xind,yind).*(l-gnd(xind,yind)+.o(l-sand(xind,yind).*sand(xind+l,yind).*(l-gnd(xind+l,yind).*sand(xind+l,yind+l);OosandNew(xind+1,yind+l)=osand(xind+l,yind+l)+.o(l-sand(xind+l,yind+l).*.o(sand(xind+l,yind).*(l-gnd(xind+l,yind)+.(l-sand(xind+l,yind).*sand(xind,yind).(l-gnd(xind,yind).*sand(xind,yind+l);%scramblethesitestomakeitlookbettertempi=sandNew(xind,yind+l).*vary(xind,yind+l)+.sandNew(xind+l,yind+l).*varyl(xind,yind+l);temp2=sandNew(xind+l,yind+l).*vary(xind,yind+l)+sandNew(xind,yind+l).*varyl(xind,yind+l);sandNew(xind,yind+l)=tempi;sandNew(xind+l,yind+l)=temp2;sand=sandNew;References1 CellularAutomataMachinesbyTommasoToffoliandNormanMargolus,MITPress,1987.2 CellularAutomataModelingofPhysicalSystemsbyBastienChopardandMichelDroz,CambridgeUniversityPress,1998.MatlabcodeconsiderationsThefollowingconsiderationswillbeillustratedusingaMatlabPrOgranlwhichcomputesConway,slife.Partsofthecodeareexplainedbelow.Matrixesandimagescanbeconvertedtoone-another,sodisplayisStraighforward.IfthearraycellscontainsthebinarystateofeachcellandthearrayZcontainszeros,thenthecatcommandbuildsanRGBimage,whichisdisplayedbytheimagecommand.Theimagecommandalsoreturnsahandle.imh=image(cat(3,cells,z,z);set(imh,erasemode,none,)axisequalaxistightMatrixesandimagescanbeconvertedtoone-another,soinitialconditionsmaybecomputedbymatrixorgraphicscomman