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