Changeset - a31c9fa6dae1
Parent rev.
Child rev.
[Not reviewed]
Merge
0 0 31
Tom Bannink - 8 years ago 2017-09-18 11:15:45
tom.bannink@cwi.nl
Merge remote-tracking branch 'cwigit/master'
9 files changed with 1045 insertions and 0 deletions:
README
6
exact/compute_linesegment.m
313
montecarlo/Makefile
10
montecarlo/data/data_reach_chain_end.m
319
montecarlo/fenwicktree.hpp
51
montecarlo/plots.m
30
montecarlo/resampler.hpp
100
montecarlo/resampler_cycle.cpp
124
montecarlo/resampler_segment.cpp
92
0 comments (0 inline, 0 general)
README
➞
Show inline comments
 
new file 100644
 
# 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.

			




	
	
	
		
 

			




        

    


    
    

    

        

                

                    exact/compute_linesegment.m
                

                

                  
                      
                        

                      ➞

                      
                        
                  

                  
                      
                  
                      
                  
                      
                  
                      
                  
                  
                

                

                    Show inline comments
                    
                

        

        

            



	
	
		
 
new file 100644

			




	
	
	
		
 
(* ::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]

			




        

    


    
    

    

        

                

                    montecarlo/Makefile
                

                

                  
                      
                        

                      ➞

                      
                        
                  

                  
                      
                  
                      
                  
                      
                  
                      
                  
                  
                

                

                    Show inline comments
                    
                

        

        

            



	
	
		
 
new file 100644

			




	
	
	
		
 
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)

			




	
	
	
		
 

			




        

    


    
    

    

        

                

                    montecarlo/data/data_reach_chain_end.m
                

                

                  
                      
                        

                      ➞

                      
                        
                  

                  
                      
                  
                      
                  
                      
                  
                      
                  
                  
                

                

                    Show inline comments
                    
                

        

        

            



	
	
		
 
new file 100644

			




	
	
	
		
 
(* 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}

			




	
	
	
		
 
}

			




        

    


    
    

    

        

                

                    montecarlo/fenwicktree.hpp
                

                

                  
                      
                        

                      ➞

                      
                        
                  

                  
                      
                  
                      
                  
                      
                  
                      
                  
                  
                

                

                    Show inline comments
                    
                

        

        

            



	
	
		
 
new file 100644

			




	
	
	
		
 
// 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 <typename T>

			




	
	
	
		
 
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;

			




	
	
	
		
 
    }

			




	
	
	
		
 
};

			




        

    


    
    

    

        

                

                    montecarlo/plots.m
                

                

                  
                      
                        

                      ➞

                      
                        
                  

                  
                      
                  
                      
                  
                      
                  
                      
                  
                  
                

                

                    Show inline comments
                    
                

        

        

            



	
	
		
 
new file 100644

			




	
	
	
		
 
(* ::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}]]]

			




        

    


    
    

    

        

                

                    montecarlo/resampler.hpp
                

                

                  
                      
                        

                      ➞

                      
                        
                  

                  
                      
                  
                      
                  
                      
                  
                      
                  
                  
                

                

                    Show inline comments
                    
                

        

        

            



	
	
		
 
new file 100644

			




	
	
	
		
 
#include "fenwicktree.hpp"

			




	
	
	
		
 
#include <cassert>

			




	
	
	
		
 
#include <random>

			




	
	
	
		
 
#include <vector>

			




	
	
	
		
 

			




	
	
	
		
 
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 <typename Derived, typename RNG>

			




	
	
	
		
 
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<VState>& 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<VState> state;

			




	
	
	
		
 
    BIT<unsigned int> fwTree; // fenwicktree

			




	
	
	
		
 
    size_t nBads;

			




	
	
	
		
 

			




	
	
	
		
 
    RNG& rng;

			




	
	
	
		
 
    std::bernoulli_distribution distr;

			




	
	
	
		
 
};

			




	
	
	
		
 

			




        

    


    
    

    

        

                

                    montecarlo/resampler_cycle.cpp
                

                

                  
                      
                        

                      ➞

                      
                        
                  

                  
                      
                  
                      
                  
                      
                  
                      
                  
                  
                

                

                    Show inline comments
                    
                

        

        

            



	
	
		
 
new file 100644

			




	
	
	
		
 
#include <iostream>

			




	
	
	
		
 
#include <random>

			




	
	
	
		
 
#include "resampler.hpp"

			




	
	
	
		
 

			




	
	
	
		
 
std::mt19937 mt(std::random_device{}());

			




	
	
	
		
 

			




	
	
	
		
 
class CycleResampler : public ResamplerBase<CycleResampler, std::mt19937> {

			




	
	
	
		
 
  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;

			




	
	
	
		
 
}

			




        

    


    
    

    

        

                

                    montecarlo/resampler_segment.cpp
                

                

                  
                      
                        

                      ➞

                      
                        
                  

                  
                      
                  
                      
                  
                      
                  
                      
                  
                  
                

                

                    Show inline comments
                    
                

        

        

            



	
	
		
 
new file 100644

			




	
	
	
		
 
#include <iostream>

			




	
	
	
		
 
#include <random>

			




	
	
	
		
 
#include "resampler.hpp"

			




	
	
	
		
 

			




	
	
	
		
 
std::mt19937 mt(std::random_device{}());

			




	
	
	
		
 

			




	
	
	
		
 
class SegmentResampler : public ResamplerBase<SegmentResampler, std::mt19937> {

			




	
	
	
		
 
  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;

			




	
	
	
		
 
}

			




        

    



    


    



    



      

      





    
    0 comments (0 inline, 0 general)
    





    


  

  







  


    




    








    
        Server instance: kalli-3714881
    
    
        This site is powered by
            Kallithea,
        which is
        © 2010–2022 by various authors & licensed under GPLv3.
            – Support