From a31c9fa6dae100b79b392f3115bb4a56b5818e3e 2017-09-18 11:15:45 From: Tom Bannink Date: 2017-09-18 11:15:45 Subject: [PATCH] Merge remote-tracking branch 'cwigit/master' --- diff --git a/README b/README new file mode 100644 index 0000000000000000000000000000000000000000..4ededab8761dd2ee769cbadac6ee85aa4fe3930c --- /dev/null +++ b/README @@ -0,0 +1,6 @@ +# Numerical analysis of resample markov chain + +The montecarlo directory contains cpp programs to run the markov chain and sample from it. + +The exact directory contains Mathemaitca notebooks to symbolically compute exact results. + diff --git a/exact/compute_linesegment.m b/exact/compute_linesegment.m new file mode 100644 index 0000000000000000000000000000000000000000..2fae35af144edf581f40754e6f8adb6b059beb2a --- /dev/null +++ b/exact/compute_linesegment.m @@ -0,0 +1,313 @@ +(* ::Package:: *) + +(* ::Section:: *) +(*Choice of symbols*) + + +(* ::Text:: *) +(*Sample a 0=BAD with probability p*) +(*Sample a 1=GOOD with probability 1-p = q*) +(*Pick a random 0 and resample.*) +(*All-1 state is the final good state.*) + + +(* ::Text:: *) +(*Goal: compute P(Z^n | start in 011..11) on the finite chain {1,...,n}*) +(*Method: sum over (k>=1) of P(vertex n gets zero for the first time at step k)*) +(*Let M be the Markov Chain on n-1 vertices such that if site n becomes zero, then it goes to a special state with no outgoing arrows*) +(*Furthermore the all-good state has no self-loop (but has incoming paths*) +(*This last thing seems like it is not required but removes the 1 eigenvalue*) +(*so I guess it is good for inverting*) +(*Let si be the initial state, let sf be the final special state*) +(*Goal is then to compute sum over (k>=1) of sf . M^k . si*) +(*This is equal to sf . (1 - M)^(-1) . si - sf . si*) +(*which is equal to sf . (1 - M)^(-1)*) + + +(* ::Chapter:: *) +(*Markov Chain Matrix*) + + +(* Matrix as described above! *) +generateLeakingMatrix[n_,p_]:=Module[{mcM,bits,numBads,j1,j2,j3,newbits,prob}, +n = n-1; (* See description above *) +mcM=SparseArray[{},{2^n+1,2^n+1},0]; + +Do[ + (* Compute neighbours of i *) + bits=IntegerDigits[i,2,n]; + numBads=Count[bits,0]; + (* Special cases: j==1 and j==n *) + If[bits[[1]]==0, (* j==1 *) + newbits=bits; + Do[ + newbits[[{1,2}]]={b1,b2}; + prob=((1-p)^(b1+b2)*(p)^(2-(b1+b2)))/numBads; + (* Set the matrix entry *) + mcM[[i+1,FromDigits[newbits,2]+1]]+=prob; + ,{b1,0,1},{b2,0,1}]; + ]; + If[bits[[n]]==0, (* j==n *) + (* This 'n' is actually the n-1 vertex. So we check if it sets the real n vertex to zero. *) + (* Potentially go into special state *) + newbits=bits; + Do[ + prob=((1-p)^(b1+b2+b3)*(p)^(3-(b1+b2+b3)))/numBads; + If[b3==0, + (* Special state *) + mcM[[i+1,-1]]+=prob; + , + (* Normal state *) + newbits[[{n-1,n}]]={b1,b2}; + (* Set the matrix entry *) + mcM[[i+1,FromDigits[newbits,2]+1]]+=prob; + ]; + ,{b1,0,1},{b2,0,1},{b3,0,1}]; + ]; + Do[ (* 2 <= j <= n-1 *) + (* For all sites, check if it is BAD *) + If[bits[[j]]==0, + newbits=bits; + Do[ + newbits[[{j-1,j,j+1}]]={b1,b2,b3}; + prob=((1-p)^(b1+b2+b3)*(p)^(3-(b1+b2+b3)))/numBads; + (* Set the matrix entry *) + mcM[[i+1,FromDigits[newbits,2]+1]]+=prob; + ,{b1,0,1},{b2,0,1},{b3,0,1}]; + ]; (* endif bits[[j]] *) + ,{j,2,n-1}]; +,{i,0,2^n-1}]; + +mcM +] + + +(* ::Section:: *) +(*Symbolic evaluation*) + + +(* ::Text:: *) +(*Let M1 be the process without the self-loop at the all-1 state but with paths towards that state.*) +(*Let M2 be the process with the all-1 state and all its incoming arrows completely removed.*) +(*Let n.R be the expected number of resamples times n and \[Rho] the initial state. Let J be a vector with all entries equal to 1 with size matching the equation.*) +(*Then we have the following.*) +(*Note that in case 1 we have \[Rho]\[CenterDot]J=1 but in case 2 we have \[Rho]\[CenterDot]J=1-(1-p)^n*) + + +(* ::DisplayFormula:: *) +(*nR=Subscript[sum, k>=1] P[k more or resamples]*) +(*nR=Subscript[sum, k>=1] \[Rho]\[CenterDot]*) +(*\!\(\*SubsuperscriptBox[\(M\), \(1\), \(k\)]\)\[CenterDot]J*) +(*nR=Subscript[sum, k>=0] \[Rho]\[CenterDot]*) +(*\!\(\*SubsuperscriptBox[\(M\), \(2\), \(k\)]\)\[CenterDot]J*) +(*nR=\[Rho]\[CenterDot]1/(Id-Subscript[M, 1])\[CenterDot]J-\[Rho].J*) +(*nR=\[Rho]\[CenterDot]1/(Id-Subscript[M, 2])\[CenterDot]J*) + + +(* ::Text:: *) +(*So we need to solve a . (Id-M)^-1 . b*) +(*We do this by solving (Id-M).x = b for x, and then computing a.x*) + + +Mfull=generateMatrix[nSites,p]; + +Mleak1=Mfull; +Mleak1[[-1,-1]]=0; (* Remove the self loop from the all-1 state *) + +Mleak2=Mfull[[;;-1,;;-1]]; (* Completely remove all-1 state *) + + +(* ::Subsection:: *) +(*Contributions from specific starting states*) + + +slotPositions={3}; +configuration=ConstantArray[1,nSites]; +configuration[[slotPositions]]=0; + +initialState2=ConstantArray[0,2^nSites]; +initialState2[[FromDigits[configuration,2]+1]]=1; + + +(* ::Subsection:: *) +(*Expected number of resamples using inverse method*) + + +(* This function definition causes automatic memoization of results *) +getIdMinusM1[n_]:=getIdMinusM1[n]=Module[{nSites,Mfull,Mleak1,IdMinusM1}, +nSites=n; +PrintTemporary["Generating matrix (Id-M) for n = ",n," and general p."]; +Mfull=generateMatrix[nSites,p]; +Mleak1=Mfull; +Mleak1[[-1,-1]]=0; (* Remove the self loop from the all-1 state *) +IdMinusM1=IdentityMatrix[2^nSites,SparseArray]-Mleak1; +IdMinusM1 +] + +(* Replacement rules do not work on sparsearrays. Instead we use this, taken shamelessly from stackoverflow *) +ReplaceElement[s_SparseArray, rule_]:=SparseArray[ArrayRules[s]/.rule,Dimensions[s]] + +computeFullExpression[n_,pValue_]:=Module[{nSites,initialState,IdMinusM1,inv1,J}, +(*PrintTemporary["Generating initial state for n = ",n,", p = ",pValue];*) +nSites=n; +initialState=Table[numBads=Count[IntegerDigits[i,2,nSites],0]; ((pValue)^(numBads)*(1-pValue)^(nSites-numBads)) , {i,0,2^nSites-1}]; +(*PrintTemporary["Getting (Id-M) matrix."];*) +IdMinusM1=getIdMinusM1[n]; +(*PrintTemporary["Entering p into matrix"];*) +IdMinusM1=Simplify[ReplaceElement[IdMinusM1,p->pValue]]; + +J=ConstantArray[1,2^nSites]; +PrintTemporary["Solving (Id-M).x = b for n = ",n," , p = ",pValue]; +inv1=LinearSolve[IdMinusM1,J]; +(*PrintTemporary["Computing final inner product"];*) +Simplify[(initialState.inv1-1)/n] +] + + +(expected=computeFullExpression[9,2/10])//AbsoluteTiming + + +Simplify[expected] +FullSimplify[expected] +Simplify[expected/nSites] +FullSimplify[expected/nSites] + + +(* ::Subsection:: *) +(*Expected number of resamples using power series method*) + + +(* ::Text:: *) +(*Using the expression for M1 as above.*) + + +computeSeries[n_,maxPower_]:=Module[{nSites,initialState,Mfull,Mleak1,nRsum,curState}, +PrintTemporary["Generating matrix M for n = ",n]; +nSites=n; +initialState=Table[numBads=Count[IntegerDigits[i,2,nSites],0]; ((p)^(numBads)*(1-p)^(nSites-numBads)) , {i,0,2^nSites-1}]; +Mfull=generateMatrix[nSites,p]; +Mleak1=Mfull; +Mleak1[[-1,-1]]=0; (* Remove the self loop from the all-1 state *) +PrintTemporary["Computing M^k, for k ="]; +nRsum=0; +curState=initialState; +Monitor[Do[ + curState=Simplify[curState.Mleak1]; + (* nRsum+=curState.ConstantArray[1,2^nSites]; (* curState . J *) *) + nRsum += Total[curState]; +,{power,maxPower}],power]; + +Series[nRsum/n,{p,0,maxPower}] +] + + +allSeries=Table[computeSeries[n,20],{n,3,7}]; + + +computeSeries[8,20] + + +(* ::Subsection:: *) +(*Compute and save coefficients*) + + +getRandom[]:=(Floor[RandomReal[]*10000000]/10000000); +pvalues=Table[getRandom[],{i,1,100}]; +values1=Map[{#,computeFullExpression[10,#]}&,pvalues]; + + +For[i=1,i<500,i++, +pValue=0/1000+i/500; +Rvalue=computeFullExpression[12,pValue]; +PutAppend[{12,pValue,Rvalue},"functionValuesN12.m"] +PutAppend[",","functionValuesN12.m"] +] + + +For[i=1,i<800,i++, +pValue=0/1000+i/800; +Rvalue=computeFullExpression[13,pValue]; +PutAppend[{13,pValue,Rvalue},"functionValuesN13.m"] +] + + +(* ::Section:: *) +(*Saved results - Expected number of resamples*) + + +(* expectx is for n=x. In these, p is the probability of sampling a GOOD instead of a BAD. *) +expect3=-1+1/p^3; +expect4=2/3 (-3+1/p^4+2/p); +expect5=(6-p+p^2+15 p^3+25 p^4-46 p^5+6 p^6-6 p^7)/(15 p^5+3 p^7); +expect6=-((-36-4 p-9 p^2-106 p^3-212 p^4-472 p^5+494 p^6+56 p^7+222 p^8+60 p^9+43 p^10-38 p^11+2 p^12)/(p^6 (162+74 p+97 p^2+33 p^3+7 p^4-13 p^5))); +expectArray1={expect3,expect4,expect5,expect6}; +expectArray2={expect3/3,expect4/4,expect5/5,expect6/6}; + + +Plot[expectArray2,{p,0,1},PlotRange->{0,5}] + + +(* n is already divided out!! and p is probability of BAD *) +expect3new=1/3 (-1+1/(1-p)^3); +expect4new=1/6 (-3+1/(1-p)^4+2/(1-p)); +expect5new=(-((p (90-300 p+435 p^2-325 p^3+136 p^4-36 p^5+6 p^6))/(15 (-1+p)^5 (6-2 p+p^2)))); +expect6new=(-((p (-2160+10620 p-22968 p^2+27378 p^3-18524 p^4+5464 p^5+1804 p^6-2583 p^7+1160 p^8-243 p^9+14 p^10+2 p^11))/(6 (-1+p)^6 (360-330 p+108 p^2+69 p^3-58 p^4+13 p^5)))); +expect7new=(-((p (6804000-51256800 p+181666800 p^2-401090760 p^3+619134390 p^4-711799410 p^5+634209054 p^6-449931963 p^7+258856333 p^8-121948365 p^9+47013236 p^10-14674393 p^11+3649155 p^12-711665 p^13+106776 p^14-11508 p^15+672 p^16))/(21 (-1+p)^7 (324000-820800 p+1090800 p^2-930360 p^3+582390 p^4-286260 p^5+110674 p^6-32233 p^7+6767 p^8-995 p^9+81 p^10)))); +expectArray={expect3new,expect4new,expect5new,expect6new,expect7new}; + + +Plot[expectArray,{p,0,1},PlotRange->{0,100}] + + +expect14partial=0 + 1 p^1 + 2 p^2 + 11/3 p^3 + 58/9 p^4 + 997/90 p^5 + 151957/8100 p^6 + 53492437/1701000 p^7 + 4671301969/89302500 p^8 + 5190549604343/60011280000 p^9 + 5979070536041027/42007896000000 p^10 + 135837147267480716807/582229438560000000 p^11 + 10253109683219103482669/26899000061472000000 p^12 + 150493768190667835592067474011/242333091553801248000000000 p^13 + 88145539086070409290770156317495909/87327152872327817729280000000000 p^14; + + +SeriesCoefficient[expect3new,{p,0,k},Assumptions->k>0] +SeriesCoefficient[expect4new,{p,0,k},Assumptions->k>0] +SeriesCoefficient[expect5new,{p,0,k},Assumptions->k>0] + + +a=Table[SeriesCoefficient[expect4new-expect3new,{p,0,k},Assumptions->k>0],{k,1,20}] +b=Table[SeriesCoefficient[expect5new-expect4new,{p,0,k},Assumptions->k>0],{k,1,20}] +c=Table[SeriesCoefficient[expect6new-expect5new,{p,0,k},Assumptions->k>0],{k,1,20}] +d=Table[SeriesCoefficient[expect7new-expect6new,{p,0,k},Assumptions->k>0],{k,1,20}] + +TableForm[{a,b,c,d}//N] + + +(* ::Subsection:: *) +(*Verify Andras' rational functions*) + + +getRandom[]:=(Floor[RandomReal[]*10000000]/10000000) +pvalues=Table[getRandom[],{i,1,100}] + + +values1=Map[{#,computeFullExpression[10,#]}&,pvalues]; +values2=Map[{#,(1/10)*expect10test/.{p->(1-#)}}&,pvalues]; +values2-values1 + + +(* ::Subsection:: *) +(*Investigate if all coefficients are non-negative*) + + +(* These are from Andras' file. Not divided by n yet and they are p\[Rule](1-p) *) +expect7test=-((-1 + p)*(15552 + p*(23760 + p*(34812 + p*(94908 + p*(222281 + p*(529516 + p*(1244215 + p*(1181257 + p*(1303681 + p*(1078303 + p*(608581 + p*(301063 + p*(101026 + p*(48821 + 12*p*(1233 + 7*p*(9 + 8*p)))))))))))))))))/(3*p^7*(44064 + p*(47616 + p*(71978 + p*(72983 + p*(48375 + p*(25115 + p*(7949 + p*(4197 + p*(1457 + p*(185 + 81*p))))))))))); +expect8test=(2*(-1 + p)*(44789760 + p*(139968000 + p*(303730560 + p*(679126512 + p*(1564082280 + p*(3726406168 + p*(8981205734 + p*(21485730797 + p*(38688270011 + p*(59923239869 + p*(79101090358 + p*(90734320109 + p*(89577453433 + p*(76026071407 + p*(55701775103 + p*(34983845549 + p*(18320307915 + p*(7343568305 + p*(1757401087 + p*(-223933111 + p*(-502060877 + p*(-315218371 + p*(-127155145 + p*(-39218742 + p*(-7828666 + p*(-1259858 + p*(-108337 + 150*p))))))))))))))))))))))))))))/(3*p^8*(-496419840 + p*(-1439176320 + p*(-3149570400 + p*(-5275650076 + p*(-7302622798 + p*(-8374860033 + p*(-8003980588 + p*(-6458041992 + p*(-4429656977 + p*(-2571384803 + p*(-1196664656 + p*(-388547559 + p*(-34974176 + p*(54650947 + p*(46874430 + p*(21894016 + p*(7362717 + p*(1657193 + p*(282806 + p*(27971 + 138*p))))))))))))))))))))); +expect9test=(14860167413760000 + p*(69215357553868800 + p*(200045951259770880 + p*(501753501908140032 + p*(1195146667726110720 + p*(2837237806479826944 + p*(6835141854725815296 + p*(16613543019874046208 + p*(40322016987813542144 + p*(87140086246363335040 + p*(166271656987524472832 + p*(279677853621781882688 + p*(415243711415104491744 + p*(542763239180836655456 + p*(618270674284253279912 + p*(599031597140865754981 + p*(463712082941731806804 + p*(226341699024787274328 + p*(-64228354607723832763 + p*(-339718231913772554127 + p*(-536556042700326483047 + p*(-619440820795896125668 + p*(-591727241700344933041 + p*(-487669723043225797092 + p*(-353188323854056822391 + p*(-226925248060543911149 + p*(-130068538893362461174 + p*(-66753544274542526574 + p*(-30769595187485936760 + p*(-12788017674210780563 + p*(-4816959007389203482 + p*(-1655060595318356059 + p*(-522904161346608018 + p*(-153811377319192586 + p*(-42983832966413783 + p*(-11758274800903707 + p*(-3227327741061583 + p*(-882100803393578 + p*(-228440163416365 + p*(-52518671538358 + p*(-10110092242971 + p*(-1527435215729 + p*(-166388566824 + p*(-11673329816 + p*(-406660450 + 4432849*p)))))))))))))))))))))))))))))))))))))))))))))/(2*p^9*(247979043717120000 + p*(1404556345973145600 + p*(4881679163642511360 + p*(12906262082477297664 + p*(28115590828074872832 + p*(52489318062390939648 + p*(85854608498957791872 + p*(124717889602448020720 + p*(162427702886973575968 + p*(190958305285669777432 + p*(203625593714258894940 + p*(197495414430604412153 + p*(174434489712187178104 + p*(140292913909359211099 + p*(102646528703039262092 + p*(68213672288407176716 + p*(41099902894845750809 + p*(22414213674226885613 + p*(11050572923278461790 + p*(4923881668005366819 + p*(1985990817617932439 + p*(728106834571938844 + p*(244474935800587786 + p*(76102417054683974 + p*(22376867005851128 + p*(6369114377239928 + p*(1798946431418834 + p*(509761338004573 + p*(142442762616972 + p*(37575586464977 + p*(8899420095010 + p*(1807947361702 + p*(299533377719 + p*(38240308283 + p*(3554076220 + p*(214908319 + 5480061*p))))))))))))))))))))))))))))))))))))); +expect10test=(29951249940995899392000000000000 + p*(321277074367082680811520000000000 + p*(1876043724828201407899238400000000 + p*(7957803700453574327899521024000000 + p*(27600849993152557107086371061760000 + p*(83708933060723217270945022476288000 + p*(232218200751572343369207785088614400 + p*(608803390669902070228281738846535680 + p*(1544273076828249644749847550461214720 + p*(3849000601054049367818531928453808128 + p*(9457009167395061143735316314450116608 + p*(22710137800601382492606555648833703936 + p*(52551319769714710101688154622944308224 + p*(115641157984597308322202359977324043776 + p*(239788039117985529108056897877275139072 + p*(466100810588291182793159375682571746048 + p*(847251806016175584466472956114635493248 + p*(1438845086263694938869869282245238450304 + p*(2281977933502502912556098358195689039488 + p*(3378207556119970655213352274466629112656 + p*(4662886865029498628962899390901557046004 + p*(5987475795830163704672450027075629130526 + p*(7123210965387283810827362410322037722354 + p*(7794836704260181065123102779965205340710 + p*(7743305046919138937831260351903602893220 + p*(6803120645273157342699268581459532854919 + p*(4969160722522289163522693119071239386755 + p*(2425833191587010847817187362952182398138 + p*(-478352352834109570300633185236348222377 + p*(-3308097287042200048938114586893251370936 + p*(-5647142522259365647779107129972197384766 + p*(-7197434992181972261186988755990921900697 + p*(-7839379187332261846531272309074275077586 + p*(-7637042760904826586923204492351454075152 + p*(-6793489561246748633179699757769394119944 + p*(-5578046078032320494940208847629640514464 + p*(-4252848140299219686443654842994289539291 + p*(-3020200540705779498941875503912519350604 + p*(-1999991026871906523995508917013136497134 + p*(-1234172886386747806789083278358669408871 + p*(-707890014436094626019080261007720401010 + p*(-375446170592945176587368390943415271029 + p*(-182392895811372280963423962179213451343 + p*(-79725884110136311249882273936589466620 + p*(-30193886820070827791743490097080535842 + p*(-8941903247106160500055796555591352006 + p*(-1194181997123529342246730283203054756 + p*(896587678060109594732754595405917406 + p*(1022852992210376792009315073725222268 + p*(684400659563903078679594347449757512 + p*(370212414902460927291178517581455088 + p*(174699951385056058300967356493852520 + p*(74011471448409088915716718895577422 + p*(28491927116198785339903578231282326 + p*(10006029890177055926113350394773918 + p*(3200260476270785946496244273521456 + p*(925386205163002054605998560468560 + p*(238206904285052006364586403532317 + p*(52866298651863763741544121057741 + p*(9348468391347117877049424284698 + p*(958858054559840471545883772371 + p*(-136926411301609161234965497492 + p*(-115008021854642384362351679374 + p*(-41244982952337852495390955333 + p*(-11070559287492334189423637514 + p*(-2447709113603523067317445528 + p*(-458350828193163804833495576 + p*(-72994618560338171901466600 + p*(-9788368655259968045182199 + p*(-1078274571766986784089992 + p*(-92853720471534248219698 + p*(-5575557597562066201327 + p*(-147991530451600663714 + 3*p*(2939473042278627387 + p*(332111127808852965 + 2*p*(3888241226376788 + 3*p*(-41445657273564 + 117378618257*p)))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))/(3*p^10*(682638904905198206976000000000000 + p*(8173019706412830881218560000000000 + p*(53323699054157660941398835200000000 + p*(248891573008669764067233103872000000 + p*(924687407612579930813801724641280000 + p*(2891944600003682820203778769158144000 + p*(7875333215307845589347486063080243200 + p*(19092145163533430334925059167734333440 + p*(41842060723846599378217956818057134080 + p*(83819625249296263520678599475268452352 + p*(154744761966786179578080182564885237760 + p*(264929077641518034992054266706731634688 + p*(422646641464611331099548275766631199616 + p*(630658407565185428744866514944008740608 + p*(882809399516424423060949095297359647184 + p*(1162011757403374712933484162759274746352 + p*(1440861354989945223972829409569972470756 + p*(1685474601631641160870079255922808661076 + p*(1862019471896150106118299899011256684852 + p*(1944275887100004592293260532580143761700 + p*(1919928618321355019268022528887019266326 + p*(1793535536249536699497366468967249191017 + p*(1585197435822223572162425675827591318758 + p*(1325440982570820498930344210935822181206 + p*(1048072002217543081106112398597412481982 + p*(783251663870790984476513634219923281927 + p*(552662728415568107148138472304112723190 + p*(367644054985132215092587551798222270155 + p*(230074222242873458531674653777240431720 + p*(135022211219245045125800091088296752850 + p*(73954257805233177084811932744909118208 + p*(37520406174308360962398062351812812137 + p*(17409390768768132523512955359881317010 + p*(7212631289630161215602424043408525090 + p*(2527911717628399131459446318200577066 + p*(630521501315295907688039566779316660 + p*(-3264167000026769674250037142301412 + p*(-139693374141558418039783207724595386 + p*(-119984125642330341935345402175013164 + p*(-73648953735741925869248816811906408 + p*(-38131255307772761702405407482815932 + p*(-17517515974055060595862322054303126 + p*(-7290643414588233704648436825349684 + p*(-2773291843079741728329148625684886 + p*(-966453715014487415372765637669408 + p*(-307780848263943355417434290190728 + p*(-88890463102686920044975809291432 + p*(-22925956093782388618998098510906 + p*(-5118160244307518827820812609760 + p*(-917260979612371876699699161420 + p*(-98655372948554285790307210984 + p*(11501501864542733975686221040 + p*(10761940854651464305687453022 + p*(3945723859789135370872243241 + p*(1073726996470422065262668662 + p*(240103317502273227850391290 + p*(45421861294678062478428318 + p*(7302286255304887036723447 + p*(987936432020600117567806 + p*(109789200467587590966019 + p*(9555322669271838491048 + p*(585926357618901832414 + p*(17271126997972246952 + p*(-709836412858689543 + 2*p*(-46765727964771227 + 3*p*(-442866761590467 + 6322265536535*p))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))))); + +expectArray={expect3new,expect4new,expect5new,expect6new,expect7new}; +expectArray=expectArray~Join~({expect7test/7,expect8test/8,expect9test/9,expect10test/10}/.{p->(1-p)}); + + +Plot[expectArray,{p,0,1},PlotRange->{0,100},PlotLabels->Automatic] + + +denominators=Denominator/@expectArray; +sols=Map[N[Simplify[Solve[#==0,p]]]&,denominators]; +sols2=sols[[All,All,1,2]]; +norms=Map[Norm,sols2,{2}]; +norms=Map[Sort,norms]; +TableForm[norms] diff --git a/montecarlo/Makefile b/montecarlo/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..d0d660a1548a472de516bb04bf48602c046fb41c --- /dev/null +++ b/montecarlo/Makefile @@ -0,0 +1,10 @@ +CXXFLAGS += -std=c++14 -O3 +CXXFLAGS += -Wall -Wextra -Wfatal-errors -Werror -pedantic -Wno-deprecated-declarations + +TARGETS += resampler_cycle resampler_segment + +all: $(TARGETS) + +clean: + rm -f $(TARGETS) + diff --git a/montecarlo/data/data_reach_chain_end.m b/montecarlo/data/data_reach_chain_end.m new file mode 100644 index 0000000000000000000000000000000000000000..a04e21b1ac3e2850f18e2160603409fd5d76782f --- /dev/null +++ b/montecarlo/data/data_reach_chain_end.m @@ -0,0 +1,319 @@ +(* Probability of hitting site n when starting from single BAD at 0 on segment [0,n]. +Data-format: {p, n, probability, runs within timelimit, total run attempts} *) +{ +{0.1,600,0/1000,1000,1000}, +{0.6,600,0/1000,1000,1000}, +{0.601,600,0/1000,1000,1000}, +{0.602,600,0/1000,1000,1000}, +{0.603,600,0/1000,1000,1000}, +{0.604,600,0/1000,1000,1000}, +{0.605,600,0/1000,1000,1000}, +{0.606,600,0/1000,1000,1000}, +{0.607,600,0/1000,1000,1000}, +{0.608,600,0/1000,1000,1000}, +{0.609,600,0/1000,1000,1000}, +{0.61,600,0/1000,1000,1000}, +{0.611,600,0/1000,1000,1000}, +{0.612,600,0/1000,1000,1000}, +{0.613,600,0/1000,1000,1000}, +{0.614,600,0/1000,1000,1000}, +{0.615,600,0/1000,1000,1000}, +{0.616,600,0/1000,1000,1000}, +{0.617,600,0/1000,1000,1000}, +{0.618,600,0/1000,1000,1000}, +{0.619,600,0/1000,1000,1000}, +{0.62,600,0/1000,1000,1000}, +{0.621,600,0/1000,1000,1000}, +{0.622,600,0/1000,1000,1000}, +{0.623,600,0/1000,1000,1000}, +{0.624,600,0/1000,1000,1000}, +{0.625,600,0/1000,1000,1000}, +{0.626,600,0/1000,1000,1000}, +{0.627,600,0/1000,1000,1000}, +{0.628,600,0/1000,1000,1000}, +{0.629,600,0/1000,1000,1000}, +{0.63,600,0/1000,1000,1000}, +{0.631,600,0/1000,1000,1000}, +{0.632,600,1/1000,1000,1000}, +{0.633,600,3/993,993,1000}, +{0.634,600,3/997,997,1000}, +{0.635,600,12/987,987,1000}, +{0.636,600,14/981,981,1000}, +{0.637,600,22/961,961,1000}, +{0.638,600,37/960,960,1000}, +{0.639,600,49/952,952,1000}, +{0.64,600,68/948,948,1000}, +{0.641,600,68/954,954,1000}, +{0.642,600,75/947,947,1000}, +{0.643,600,104/951,951,1000}, +{0.644,600,113/954,954,1000}, +{0.645,600,138/960,960,1000}, +{0.646,600,155/960,960,1000}, +{0.647,600,140/971,971,1000}, +{0.648,600,179/970,970,1000}, +{0.649,600,186/977,977,1000}, +{0.65,600,198/979,979,1000}, +{0.651,600,192/981,981,1000}, +{0.652,600,216/985,985,1000}, +{0.653,600,268/994,994,1000}, +{0.654,600,239/989,989,1000}, +{0.655,600,259/991,991,1000}, +{0.656,600,267/996,996,1000}, +{0.657,600,287/996,996,1000}, +{0.658,600,289/998,998,1000}, +{0.659,600,283/994,994,1000}, +{0.66,600,299/998,998,1000}, +{0.661,600,310/999,999,1000}, +{0.662,600,352/997,997,1000}, +{0.663,600,326/1000,1000,1000}, +{0.664,600,342/999,999,1000}, +{0.665,600,348/999,999,1000}, +{0.666,600,368/1000,1000,1000}, +{0.667,600,384/1000,1000,1000}, +{0.668,600,373/999,999,1000}, +{0.669,600,380/1000,1000,1000}, +{0.67,600,363/1000,1000,1000}, +{0.671,600,406/1000,1000,1000}, +{0.672,600,395/1000,1000,1000}, +{0.673,600,431/1000,1000,1000}, +{0.674,600,442/1000,1000,1000}, +{0.675,600,428/1000,1000,1000}, +{0.676,600,435/1000,1000,1000}, +{0.677,600,434/1000,1000,1000}, +{0.678,600,434/1000,1000,1000}, +{0.679,600,449/1000,1000,1000}, +{0.68,600,476/1000,1000,1000}, +{0.681,600,470/1000,1000,1000}, +{0.682,600,460/1000,1000,1000}, +{0.683,600,463/1000,1000,1000}, +{0.684,600,480/1000,1000,1000}, +{0.685,600,505/1000,1000,1000}, +{0.686,600,491/1000,1000,1000}, +{0.687,600,486/1000,1000,1000}, +{0.688,600,502/1000,1000,1000}, +{0.689,600,521/1000,1000,1000}, +{0.69,600,517/1000,1000,1000}, +{0.691,600,542/1000,1000,1000}, +{0.692,600,553/1000,1000,1000}, +{0.693,600,553/1000,1000,1000}, +{0.694,600,537/1000,1000,1000}, +{0.695,600,579/1000,1000,1000}, +{0.696,600,541/1000,1000,1000}, +{0.697,600,547/1000,1000,1000}, +{0.698,600,552/1000,1000,1000}, +{0.699,600,563/1000,1000,1000}, +{0.7,600,556/1000,1000,1000}, +{0.75,600,754/1000,1000,1000}, +{0.8,600,871/1000,1000,1000}, +{0.9,600,977/1000,1000,1000}, +{0.1,800,0/1000,1000,1000}, +{0.6,800,0/1000,1000,1000}, +{0.601,800,0/1000,1000,1000}, +{0.602,800,0/1000,1000,1000}, +{0.603,800,0/1000,1000,1000}, +{0.604,800,0/1000,1000,1000}, +{0.605,800,0/1000,1000,1000}, +{0.606,800,0/1000,1000,1000}, +{0.607,800,0/1000,1000,1000}, +{0.608,800,0/1000,1000,1000}, +{0.609,800,0/1000,1000,1000}, +{0.61,800,0/1000,1000,1000}, +{0.611,800,0/1000,1000,1000}, +{0.612,800,0/1000,1000,1000}, +{0.613,800,0/1000,1000,1000}, +{0.614,800,0/1000,1000,1000}, +{0.615,800,0/1000,1000,1000}, +{0.616,800,0/1000,1000,1000}, +{0.617,800,0/1000,1000,1000}, +{0.618,800,0/1000,1000,1000}, +{0.619,800,0/1000,1000,1000}, +{0.62,800,0/1000,1000,1000}, +{0.621,800,0/1000,1000,1000}, +{0.622,800,0/1000,1000,1000}, +{0.623,800,0/1000,1000,1000}, +{0.624,800,0/1000,1000,1000}, +{0.625,800,0/1000,1000,1000}, +{0.626,800,0/1000,1000,1000}, +{0.627,800,0/1000,1000,1000}, +{0.628,800,0/1000,1000,1000}, +{0.629,800,0/1000,1000,1000}, +{0.63,800,0/1000,1000,1000}, +{0.631,800,0/1000,1000,1000}, +{0.632,800,1/1000,1000,1000}, +{0.633,800,0/997,997,1000}, +{0.634,800,3/990,990,1000}, +{0.635,800,0/978,978,1000}, +{0.636,800,7/974,974,1000}, +{0.637,800,4/967,967,1000}, +{0.638,800,14/953,953,1000}, +{0.639,800,18/912,912,1000}, +{0.64,800,22/912,912,1000}, +{0.641,800,30/922,922,1000}, +{0.642,800,30/902,902,1000}, +{0.643,800,44/887,887,1000}, +{0.644,800,51/883,883,1000}, +{0.645,800,64/887,887,1000}, +{0.646,800,65/882,882,1000}, +{0.647,800,85/885,885,1000}, +{0.648,800,84/887,887,1000}, +{0.649,800,85/886,886,1000}, +{0.65,800,118/897,897,1000}, +{0.651,800,132/890,890,1000}, +{0.652,800,151/921,921,1000}, +{0.653,800,187/912,912,1000}, +{0.654,800,150/914,914,1000}, +{0.655,800,208/931,931,1000}, +{0.656,800,201/923,923,1000}, +{0.657,800,213/927,927,1000}, +{0.658,800,245/941,941,1000}, +{0.659,800,259/964,964,1000}, +{0.66,800,271/948,948,1000}, +{0.661,800,296/956,956,1000}, +{0.662,800,306/961,961,1000}, +{0.663,800,310/970,970,1000}, +{0.664,800,337/977,977,1000}, +{0.665,800,364/979,979,1000}, +{0.666,800,346/989,989,1000}, +{0.667,800,346/981,981,1000}, +{0.668,800,358/989,989,1000}, +{0.669,800,397/990,990,1000}, +{0.67,800,384/988,988,1000}, +{0.671,800,394/997,997,1000}, +{0.672,800,433/994,994,1000}, +{0.673,800,401/997,997,1000}, +{0.674,800,433/997,997,1000}, +{0.675,800,435/996,996,1000}, +{0.676,800,442/996,996,1000}, +{0.677,800,438/999,999,1000}, +{0.678,800,440/999,999,1000}, +{0.679,800,462/999,999,1000}, +{0.68,800,445/999,999,1000}, +{0.681,800,456/1000,1000,1000}, +{0.682,800,465/1000,1000,1000}, +{0.683,800,468/1000,1000,1000}, +{0.684,800,488/1000,1000,1000}, +{0.685,800,474/1000,1000,1000}, +{0.686,800,495/1000,1000,1000}, +{0.687,800,494/1000,1000,1000}, +{0.688,800,496/1000,1000,1000}, +{0.689,800,532/1000,1000,1000}, +{0.69,800,493/1000,1000,1000}, +{0.691,800,519/1000,1000,1000}, +{0.692,800,522/1000,1000,1000}, +{0.693,800,515/1000,1000,1000}, +{0.694,800,545/1000,1000,1000}, +{0.695,800,544/1000,1000,1000}, +{0.696,800,521/1000,1000,1000}, +{0.697,800,583/1000,1000,1000}, +{0.698,800,519/1000,1000,1000}, +{0.699,800,579/1000,1000,1000}, +{0.7,800,577/1000,1000,1000}, +{0.75,800,771/1000,1000,1000}, +{0.8,800,863/1000,1000,1000}, +{0.9,800,978/1000,1000,1000}, +{0.1,1000,0/1000,1000,1000}, +{0.6,1000,0/1000,1000,1000}, +{0.601,1000,0/1000,1000,1000}, +{0.602,1000,0/1000,1000,1000}, +{0.603,1000,0/1000,1000,1000}, +{0.604,1000,0/1000,1000,1000}, +{0.605,1000,0/1000,1000,1000}, +{0.606,1000,0/1000,1000,1000}, +{0.607,1000,0/1000,1000,1000}, +{0.608,1000,0/1000,1000,1000}, +{0.609,1000,0/1000,1000,1000}, +{0.61,1000,0/1000,1000,1000}, +{0.611,1000,0/1000,1000,1000}, +{0.612,1000,0/1000,1000,1000}, +{0.613,1000,0/1000,1000,1000}, +{0.614,1000,0/1000,1000,1000}, +{0.615,1000,0/1000,1000,1000}, +{0.616,1000,0/1000,1000,1000}, +{0.617,1000,0/1000,1000,1000}, +{0.618,1000,0/1000,1000,1000}, +{0.619,1000,0/1000,1000,1000}, +{0.62,1000,0/1000,1000,1000}, +{0.621,1000,0/1000,1000,1000}, +{0.622,1000,0/1000,1000,1000}, +{0.623,1000,0/1000,1000,1000}, +{0.624,1000,0/1000,1000,1000}, +{0.625,1000,0/1000,1000,1000}, +{0.626,1000,0/1000,1000,1000}, +{0.627,1000,0/1000,1000,1000}, +{0.628,1000,0/1000,1000,1000}, +{0.629,1000,0/1000,1000,1000}, +{0.63,1000,0/1000,1000,1000}, +{0.631,1000,0/1000,1000,1000}, +{0.632,1000,0/1000,1000,1000}, +{0.633,1000,1/997,997,1000}, +{0.634,1000,0/993,993,1000}, +{0.635,1000,1/984,984,1000}, +{0.636,1000,1/962,962,1000}, +{0.637,1000,3/944,944,1000}, +{0.638,1000,2/925,925,1000}, +{0.639,1000,3/928,928,1000}, +{0.64,1000,2/917,917,1000}, +{0.641,1000,4/891,891,1000}, +{0.642,1000,6/871,871,1000}, +{0.643,1000,5/866,866,1000}, +{0.644,1000,5/861,861,1000}, +{0.645,1000,12/840,840,1000}, +{0.646,1000,12/837,837,1000}, +{0.647,1000,8/824,824,1000}, +{0.648,1000,10/817,817,1000}, +{0.649,1000,16/811,811,1000}, +{0.65,1000,22/797,797,1000}, +{0.651,1000,22/801,801,1000}, +{0.652,1000,31/810,810,1000}, +{0.653,1000,29/787,787,1000}, +{0.654,1000,35/777,777,1000}, +{0.655,1000,65/784,784,1000}, +{0.656,1000,63/799,799,1000}, +{0.657,1000,74/777,777,1000}, +{0.658,1000,81/780,780,1000}, +{0.659,1000,96/795,795,1000}, +{0.66,1000,109/812,812,1000}, +{0.661,1000,118/799,799,1000}, +{0.662,1000,122/773,773,1000}, +{0.663,1000,144/799,799,1000}, +{0.664,1000,154/805,805,1000}, +{0.665,1000,166/848,848,1000}, +{0.666,1000,191/828,828,1000}, +{0.667,1000,218/846,846,1000}, +{0.668,1000,241/844,844,1000}, +{0.669,1000,203/861,861,1000}, +{0.67,1000,278/861,861,1000}, +{0.671,1000,270/881,881,1000}, +{0.672,1000,303/872,872,1000}, +{0.673,1000,319/897,897,1000}, +{0.674,1000,327/891,891,1000}, +{0.675,1000,358/920,920,1000}, +{0.676,1000,342/919,919,1000}, +{0.677,1000,384/941,941,1000}, +{0.678,1000,375/941,941,1000}, +{0.679,1000,404/938,938,1000}, +{0.68,1000,424/950,950,1000}, +{0.681,1000,413/953,953,1000}, +{0.682,1000,443/965,965,1000}, +{0.683,1000,422/958,958,1000}, +{0.684,1000,457/975,975,1000}, +{0.685,1000,477/982,982,1000}, +{0.686,1000,464/984,984,1000}, +{0.687,1000,466/974,974,1000}, +{0.688,1000,475/989,989,1000}, +{0.689,1000,499/991,991,1000}, +{0.69,1000,525/984,984,1000}, +{0.691,1000,523/992,992,1000}, +{0.692,1000,479/997,997,1000}, +{0.693,1000,532/992,992,1000}, +{0.694,1000,569/998,998,1000}, +{0.695,1000,533/993,993,1000}, +{0.696,1000,549/995,995,1000}, +{0.697,1000,565/997,997,1000}, +{0.698,1000,557/999,999,1000}, +{0.699,1000,529/995,995,1000}, +{0.7,1000,576/999,999,1000}, +{0.75,1000,748/1000,1000,1000}, +{0.8,1000,862/1000,1000,1000}, +{0.9,1000,973/1000,1000,1000} +} diff --git a/montecarlo/fenwicktree.hpp b/montecarlo/fenwicktree.hpp new file mode 100644 index 0000000000000000000000000000000000000000..76ef445c4897b5475c5d1185ff1bd0b98925f4e6 --- /dev/null +++ b/montecarlo/fenwicktree.hpp @@ -0,0 +1,51 @@ +// Taken from +// https://kartikkukreja.wordpress.com/2013/05/11/bit-fenwick-tree-data-structure-c-implementation/ +// Has MIT licence +// Slightly modified +#pragma once + +#define LSOne(S) (S & (-S)) + +template +class BIT { + T* ft; + size_t size; +public: + // initialize a BIT of n elements to zero + BIT(size_t n) { + size = n; + ft = new T[n+1](); // () sets it to zero + } + + ~BIT() { + delete [] ft; + } + + // returns sum of the range [1...b] + T sum(size_t b) { + T sum = 0; + for (; b; b -= LSOne(b)) sum += ft[b]; + return sum; + } + + // returns sum of the range [a...b] + T sum(size_t a, size_t b) { + return sum(b) - (a == 1 ? 0 : sum(a - 1)); + } + + // update value of the k-th element by v (v can be +ve/inc or -ve/dec) + // i.e., increment or decrement kth element by a value v + void update(size_t k, T v) { + for (; k <= size; k += LSOne(k)) ft[k] += v; + } + + // divide each original frequency by c + void scaleDown(T c){ + for (int i=1 ; i<=size ; i++) ft[i] /= c; + } + + // multiply each original frequency by c + void scaleUp(T c){ + for (int i=1 ; i<=size ; i++) ft[i] *= c; + } +}; diff --git a/montecarlo/plots.m b/montecarlo/plots.m new file mode 100644 index 0000000000000000000000000000000000000000..5fa7f225fa86604d061957f8b1d5a0e61818d19a --- /dev/null +++ b/montecarlo/plots.m @@ -0,0 +1,30 @@ +(* ::Package:: *) + +(* ::Section:: *) +(*P(Z^n | start 011..11) on the chain {1,...,n}*) + + +rawdata=Import[NotebookDirectory[]<>"data/data_reach_chain_end.m"]; +data=GatherBy[rawdata,{#[[2]]}&]; +(* data[[n index, prob index, p / n / prob / runs within timelimit / total runs]] *) +nlabels=Map["n = "<>ToString[#[[1,2]]]&,data]; + + +plotZoomed=ListPlot[data[[All,All,{1,3}]], +Joined->True, +PlotMarkers->Automatic, +PlotRange->{{0.6,0.7},{0,0.5}}, +Frame->True, +ImageSize->250, +Epilog->{Dashed,Line[{{1-1/E,0},{1-1/E,1}}]} +]; + +plotBoth=ListPlot[data[[All,All,{1,3}]], +Joined->True, +FrameLabel->{"p (probability of BAD)","P(\!\(\*SuperscriptBox[\(Z\), \(n\)]\) | start 011..11)"}, +PlotMarkers->Automatic, +PlotLegends->nlabels, +PlotRange->{{0,1},{0,1}}, +Frame->True, +ImageSize->500, +Epilog->Inset[plotZoomed,Scaled[{0.32,0.6}]]] diff --git a/montecarlo/resampler.hpp b/montecarlo/resampler.hpp new file mode 100644 index 0000000000000000000000000000000000000000..8497dbbce98e1f4c275792248d651601b2e78c99 --- /dev/null +++ b/montecarlo/resampler.hpp @@ -0,0 +1,100 @@ +#include "fenwicktree.hpp" +#include +#include +#include + +enum VState { BAD = 0, GOOD = 1 }; + +// Base class for resampler +// Contains code for storing n vertices and efficiently picking +// a random BAD using a fenwick tree data structure. +// The graph structure should be specified by the subclass, +// by implementing a function +// ``void resampleNeighbors(int i) { +// resampleVertex(i); +// resampleVertex(j1); +// // ... +// resampleVertex(jk); +// }`` +template +class ResamplerBase { + public: + // p - probability of sampling a BAD + // n - total number of vertices + // rng - random number generator, usually std::mt19937 + // Default initial state is all-GOOD + ResamplerBase(float p_, size_t n, RNG& rng_) + : p(p_), state(n, GOOD), fwTree(n), nBads(0), rng(rng_), distr(p) { + assert(p >= 0.0f && p <= 1.0f); + } + + // Returns the chosen site or -1 if in all-good state + int doMove() { + if (nBads == 0) + return -1; + assert(fwTree.sum(state.size()) == nBads); + std::uniform_int_distribution<> intDist(1, nBads); + unsigned int a = intDist(rng); + // There are 'nBads' 1's in state + // Select the a'th Bad using a binary search on the fenwicktree + // Where a is 1-based, as well as the fenwick tree + int i = 0; + { + // The a-th Bad is in the range + // [x_lower, ..., x_upper-1] + // where lower and upper are 0-based + int lower = 0; + int upper = state.size(); + + while (upper - lower > 1) { + int cur = (upper + lower) / 2; + + // Check the number of 1's in [x0,...,x_cur-1] + // fwTree is 1-based so no -1 needed + if (fwTree.sum(cur) >= a) { + // The a-th 1 is in [x_lower,...,x_cur-1] + upper = cur; + } else { + // The a-th 1 is in the range [x_cur,...,x_upper-1] + lower = cur; + } + } + i = lower; + } + + // state[i] is now the a'th Bad + ((Derived*)this)->resampleNeighbors(i); + return i; + } + + size_t getN() const { return state.size(); }; + size_t numBads() const { return nBads; } + + const std::vector& getState() const { return state; } + + protected: + void setVertex(int i, VState value) { + if (state[i] == BAD) { + fwTree.update(i + 1, -1); + nBads--; + } + state[i] = value; + if (state[i] == BAD) { + fwTree.update(i + 1, 1); + nBads++; + } + } + + // Used by Derived::resampleNeighbors + void resampleVertex(int i) { setVertex(i, (distr(rng) ? BAD : GOOD)); } + + private: + float p; + std::vector state; + BIT fwTree; // fenwicktree + size_t nBads; + + RNG& rng; + std::bernoulli_distribution distr; +}; + diff --git a/montecarlo/resampler_cycle.cpp b/montecarlo/resampler_cycle.cpp new file mode 100644 index 0000000000000000000000000000000000000000..bf25d72d052f1e622eed54392a8200595d3de762 --- /dev/null +++ b/montecarlo/resampler_cycle.cpp @@ -0,0 +1,124 @@ +#include +#include +#include "resampler.hpp" + +std::mt19937 mt(std::random_device{}()); + +class CycleResampler : public ResamplerBase { + public: + CycleResampler(float p_, size_t n, bool allBads = true) + : ResamplerBase(p_, n, mt) { + assert(n >= 3); + if (allBads) { + for (auto i = 0u; i < n; ++i) { + setVertex(i, BAD); + } + } else { + setVertex(0, BAD); + } + } + + // Cycle has periodic boundary conditions + void resampleNeighbors(int i) { + resampleVertex(i); + + if (i == 0) + resampleVertex(getN() - 1); + else + resampleVertex(i - 1); + + if (i + 1 == (int)getN()) + resampleVertex(0); + else + resampleVertex(i + 1); + } + + // Returns number of steps before all good + int run(int limit = -1) { + int t; + for (t = 0; t != limit && numBads() != 0; t++) + doMove(); + return t; + } +}; + +int main() { + //if (0) { + CycleResampler sampler(0.6f, 40); + int loc; + for (int i = 0; i < 2000; ++i) { + loc = sampler.doMove(); + if (loc == -1) + break; + std::cout << '|'; + int j = 0; + for (auto x : sampler.getState()) { + if (j++ == loc) { + std::cout << 'X'; + } else { + std::cout << (x == BAD ? '*' : ' '); + } + } + std::cout << '|' << std::endl; + } + //} + + // Measure mixing time for various n,p + // starting from all-1 state + if (0) { + std::cout << "{{0,0,0}"; + for (int i = 0; i <= 8; ++i) { + float p = 0.6f + 0.01f * i; + for (int n = 5; n <= 25; ++n) { + std::cerr << "Running p = " << p << " and n = " << n + << std::endl; + for (int j = 0; j < 1000; ++j) { + CycleResampler sampler(p, n); + int time = sampler.run(500000); + std::cout << ',' << '{' << p << ',' << n << ',' << time + << '}'; + } + std::cout << std::flush; + } + } + std::cout << '}' << std::endl; + } + + // Measure probability of hitting all-good + // starting from single-bad and stopping if some condition is met + // Interesting range is p in [0.60,0.75] + std::cout << "{{0,0,0}"; + for (int i = 0; i <= 100; ++i) { + float p = 0.60f + 0.001f * i; + for (size_t n = 1000; n <= 1000; n += 200) { + std::cerr << "Running p = " << p << " and n = " << n << std::endl; + constexpr int samples = 500; + int hitCount = 0; + int otherHitCount = 0; + for (int j = 0; j < samples; ++j) { + CycleResampler sampler(p, n, false); + for (int t = 0; t != 500000; t++) { + sampler.doMove(); + if (sampler.numBads() == 0) { + hitCount++; + break; + } + //if (sampler.numBads() >= (5*n)/10) { + if (sampler.getState()[n/2] == BAD) { + otherHitCount++; + break; + } + } + } + std::cout << ',' << '{' << p << ',' << n << ',' + << float(hitCount) / float(samples) << '}'; + std::cout << std::flush; + + std::cerr << (samples - otherHitCount - hitCount) << '/' << samples + << " runs with timelimit.\n"; + } + } + std::cout << '}' << std::endl; + + return 0; +} diff --git a/montecarlo/resampler_segment.cpp b/montecarlo/resampler_segment.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6ffd1b74b21564cecdb1ceafa091c03566092e67 --- /dev/null +++ b/montecarlo/resampler_segment.cpp @@ -0,0 +1,92 @@ +#include +#include +#include "resampler.hpp" + +std::mt19937 mt(std::random_device{}()); + +class SegmentResampler : public ResamplerBase { + public: + SegmentResampler(float p_, size_t n) : ResamplerBase(p_, n, mt) { + setVertex(0, BAD); + } + + // No periodic boundary + void resampleNeighbors(int i) { + resampleVertex(i); + + if (i != 0) + resampleVertex(i - 1); + if (i + 1 != (int)getN()) + resampleVertex(i + 1); + } +}; + +void doRun(float p, size_t n) { + std::cerr << "Running p = " << p << " and n = " << n << std::endl; + + constexpr int samples = 1000; + + int terminateCount = 0; + int hitnCount = 0; + for (int j = 0; j < samples; ++j) { + SegmentResampler sampler(p, n); + for (auto t = 0u; t != 1000 * n; t++) { + sampler.doMove(); + if (sampler.numBads() == 0) { + terminateCount++; + break; + } + if (sampler.getState()[n - 1] == BAD) { + hitnCount++; + break; + } + } + } + if (terminateCount + hitnCount != samples) { + std::cerr << (samples - terminateCount - hitnCount) << '/' << samples + << " runs with timelimit!\n"; + } + + std::cout << '{' << p << ',' << n; + std::cout << ',' << hitnCount << '/' << hitnCount + terminateCount; + std::cout << ',' << hitnCount + terminateCount; + std::cout << ',' << samples; + std::cout << '}'; + std::cout << std::flush; +} + +int main() { + // Measure probability of hitting vertex N starting from single-bad + // Interesting range is p in [0.60,0.75] + std::cout << "(* Probability of hitting site n when starting from single " + "BAD at 0 on segment [0,n].\n"; + std::cout << "Data-format: {p, n, probability, runs within timelimit, " + "total run attempts} *)\n"; + std::cout << "{\n"; + for (size_t n = 600;; n += 200) { + doRun(0.10f, n); + std::cout << ',' << std::endl; + + for (int i = 0; i <= 100; ++i) { + float p = 0.60f + 0.001f * i; + doRun(p, n); + std::cout << ',' << std::endl; + } + doRun(0.75f, n); + std::cout << ',' << std::endl; + ; + doRun(0.80f, n); + std::cout << ',' << std::endl; + ; + doRun(0.90f, n); + + if (n == 1000) + break; + else + std::cout << ',' << std::endl; + ; + } + std::cout << '\n' << '}' << std::endl; + + return 0; +}