From 4b6c189e4ae468d01dc14a4316d660356f7322d7 2017-06-28 17:44:22
From: Tom Bannink <tom.bannink@cwi.nl>
Date: 2017-06-28 17:44:22
Subject: [PATCH] Add mixing time analysis using exponential fits

---

diff --git a/triangle_analysis.m b/triangle_analysis.m
index e972a7a311431b897ad816f8955efaa20e4d814c..d20b4e5fa65f4bbd6933972a84721872447805de 100644
--- a/triangle_analysis.m
+++ b/triangle_analysis.m
@@ -1,5 +1,8 @@
 (* ::Package:: *)
 
+Quit[]
+
+
 Needs["ErrorBarPlots`"]
 
 
@@ -197,6 +200,22 @@ movingAvg=Map[MovingAverage[#[[2]],2000][[1;;-1;;skipPts]]-Mean[#[[2,-20000;;-1]
 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"*)
 
@@ -234,6 +253,16 @@ getMixingTime[values_]:=Module[{avg},
     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}];
@@ -245,14 +274,23 @@ ErrorListPlot[mixingTimesBars,Joined->True,PlotMarkers->Automatic,AxesLabel->{"n
 
 (* 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}]]