AENC/switchchain Changeset - 4b6c189e4ae4 · Centrum Wiskunde & Informatica (CWI)

Changeset - 4b6c189e4ae4

Parent rev.

Child rev.

[Not reviewed]

0 1 0

Tom Bannink - 8 years ago 2017-06-28 17:44:22
tom.bannink@cwi.nl

Add mixing time analysis using exponential fits

1 file changed with 41 insertions and 3 deletions:

triangle_analysis.m

0 comments (0 inline, 0 general)

triangle_analysis.m

➞

Show inline comments

 (* ::Package:: *)
 Quit[]
 Needs["ErrorBarPlots`"]
 (* ::Section:: *)
 (*TODO*)
 (* ::Text:: *)
 (*- Triangle law exponent: gather more data*)
 (**)
 (*- Why does GCM-2 start with very low #triangles*)
 (*  Do not only consider number of standard deviations but also relative number of triangles.*)
@@ @@ -188,24 +191,40 @@ coarseData=Map[#[[2,1;;maxTime;;skipPts]]&,selectedData]; @@
 labels=Map["{n,tau} = "<>ToString[#[[1]]]&,selectedData];
 plotTimeEvol=ListPlot[coarseData,Joined->True,PlotRange->{0*minCount,maxCount},DataRange->{0,measureSkip*maxTime},Frame->True,FrameLabel->{"timesteps","number of triangles"},ImageSize->300]
 (* Map[ListPlot[#,Joined->True,PlotRange\[Rule]{minCount,maxCount},DataRange\[Rule]{0,maxTime}]&,coarseData] *)
 Export[NotebookDirectory[]<>"plots/timeevol.pdf",plotTimeEvol]
 movingAvg=Map[MovingAverage[#[[2]],2000][[1;;-1;;skipPts]]-Mean[#[[2,-20000;;-1]]]&,selectedData[[1;;-1;;5]]];
 plotMovingAvg=ListPlot[movingAvg,Joined->True,PlotRange->All,DataRange->{0,measureSkip*maxTime},Frame->True,FrameLabel->{"timesteps","number of triangles"}]
 (* ::Subsection:: *)
 (*Fit exponential to triangles time evolution*)
 fitList=Map[NonlinearModelFit[#[[2]],Exp[-(t-t0)/tmix]+c,{{tmix,1000},{t0,10000},{c,2000}},t]&,selectedData];
 (* tmix*Log[1/epsilon] gives the time it takes to get a factor epsilon close to the average *)
 (* t0 gives the time it takes to be exactly 1 triangle (in absolute value) away from the average *)
 (* Use fit["BestFitParameters"] to get parameters *)
 (* Use fit[t] to get fit value *)
 fitFuncsT=Map[#[t]&,fitList];
 timeplot1=ListPlot[coarseData,Joined->True,PlotRange->{0*minCount,maxCount},DataRange->{0,measureSkip*maxTime},PlotStyle->Opacity[0.5]];
 Show[timeplot1,Plot[fitFuncsT,{t,1,maxTime},PlotRange->All]]
 (* ::Subsection:: *)
 (*Plot success rate over "time"*)
 numPlots=10;
 selectedData=gdata[[1,-1]][[-numPlots;;-1]];
 measureSkip=100;
 maxTime=Max[Map[Length[#[[4]]]&,selectedData]];
 maxTime=10000;
 coarseData=Map[#[[4,1;;maxTime/measureSkip]]&,selectedData];
 labels=Map["{n,tau} = "<>ToString[#[[1]]]&,selectedData];
 ListPlot[coarseData,Joined->True,PlotRange->{0,100},DataRange->{0,maxTime},PlotLegends->labels]
 (* ::Subsection:: *)
 (*Compute 'mixing time'*)
 (* Compute average of last part and check when the value drops below that for the first time *)
 getMixingTime[values_]:=Module[{avg},
     (* average over the last 20 percent *)
     avg=Mean[values[[-Round[Length[values]/5];;-1]]];
     FirstPosition[values,_?(#<=avg&)][[1]]
+]
 (* Get fit of Exp[-t/tmix] *)
 getMixingTime2[values_]:=Module[{avg,etmt,fitVals},
     (* average over the last 20 percent *)
     avg=Mean[values[[-Round[Length[values]/5];;-1]]];
     etmt=FirstPosition[values,_?(#<=avg&)][[1]];
     fitVals=FindFit[values,Exp[-(t-t0)/tmix]+tavg,{{tmix,etmt/4},{t0,2*etmt},{tavg,avg}},t];
     tmix/.fitVals
 (* tmix*Log[1/epsilon] gives the time it takes to get a factor epsilon close to the average *)
 (* t0 gives the time it takes to be exactly 1 triangle (in absolute value) away from the average *)
+]
 (* gdata[[ tau index, n index, run index , {ntau, #tris, ds} ]] *)
 measureSkip=1;
 mixingTimes=Map[{#[[1,1]],(1/#[[1,1]])measureSkip * getMixingTime[#[[2]]]}&,gdata,{3}];
 mixingTimesBars=Map[
     {{#[[1,1]],Mean[#[[All,2]]]},ErrorBar[StandardDeviation[#[[All,2]]](*/Sqrt[Length[#]]*)]}&
 ,mixingTimes,{2}];
 ErrorListPlot[mixingTimesBars,Joined->True,PlotMarkers->Automatic,AxesLabel->{"n","~~mixing time divided by n"},PlotLegends->taulabels]
 (* For n fixed, look at function of tau *)
 measureSkip=1;
 mixingTimes=Map[{#[[1,2]],(1/#[[1,1]])measureSkip * getMixingTime[#[[2]]]}&,gdata,{3}];
 mixingTimes=Map[(PrintTemporary[#[[1]]];
 {#[[1,2]],(1/#[[1,1]])measureSkip * getMixingTime[#[[2]]],(1/#[[1,1]])measureSkip * getMixingTime2[#[[2]]]}
 )&,gdata,{3}];
 mixingTimesBars=Map[
     {{#[[1,1]],Mean[#[[All,2]]]},ErrorBar[StandardDeviation[#[[All,2]]]]}&
 ,mixingTimes[[All,-1]],{1}];
 ,mixingTimes[[All,-1,All,{1,2}]],{1}];
 mixingTimesBars2=Map[
     {{#[[1,1]],Mean[#[[All,2]]]},ErrorBar[StandardDeviation[#[[All,2]]]]}&
 ,mixingTimes[[All,-1,All,{1,3}]],{1}];
 Show[
 ErrorListPlot[mixingTimesBars,Joined->True,PlotMarkers->Automatic,AxesLabel->{"tau","~~mixing time divided by n, for n = 1000"},PlotRange->{{2,3},{0,30}}]
 ErrorListPlot[{mixingTimesBars,mixingTimesBars2},Joined->True,PlotMarkers->Automatic,
 AxesLabel->{"tau","~~mixing time divided by n, for n = 1000"},
 PlotRange->{{2,3},{0,30}}]
 ,Plot[(32-20(tau-2)),{tau,2,3}]]
 (* ::Subsection:: *)
 (*Plot #triangles distribution for specific (n,tau)*)
 plotRangeByTau[tau_]:=Piecewise[{{{0,30000},tau<2.3},{{0,4000},2.3<tau<2.7},{{0,800},tau>2.7}},Automatic]
 histograms=Map[Histogram[#[[All,2]],PlotRange->{plotRangeByTau[#[[1,1,2]]],Automatic}]&,averagesGrouped,{2}];
 (* TableForm[histograms,TableHeadings->{taulabels,nlabels}] *)

0 comments (0 inline, 0 general)