Skip to content

Commit 9f4ad4b

Browse files
gkanwargkanwar
committed
Adds Mathematica lib
1 parent 1a6a7b2 commit 9f4ad4b

File tree

13 files changed

+13622
-0
lines changed

13 files changed

+13622
-0
lines changed

Diff for: ‎mathlib/Makefile

+10
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
APPNAME=QCDLib
2+
3+
all: qcdlib.paclet
4+
5+
clean:
6+
rm qcdlib.paclet
7+
qcdlib.paclet: $(APPNAME)/PacletInfo.m $(APPNAME)/Kernel/*.m
8+
zip -r $@ $(APPNAME)
9+
deploy: qcdlib.paclet
10+
scp qcdlib.paclet [email protected]:web_scripts/mma_paclets/

Diff for: ‎mathlib/QCDLib/Kernel/init.m

+20
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
(* ::Package:: *)
2+
3+
(* Just load all sub-packages *)
4+
Get["QCDLib`visualize`"];
5+
Get["QCDLib`util`"];
6+
Get["QCDLib`plot`"];
7+
Get["QCDLib`analyze`"];
8+
Get["QCDLib`sn`"];
9+
Get["QCDLib`sun`"];
10+
Get["QCDLib`phase`"];
11+
(* Set up some useful options *)
12+
Module[{nb},
13+
nb = EvaluationNotebook[];
14+
If[!MatchQ[nb, $Failed],
15+
SetOptions[nb, ShowGroupOpener -> True]];
16+
On[Assert];
17+
];
18+
19+
BeginPackage["QCDLib`"];
20+
EndPackage[];

Diff for: ‎mathlib/QCDLib/PacletInfo.m

+20
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
Paclet[
2+
Name -> "QCDLib",
3+
Version -> "1.0.1",
4+
Description -> "Various utilities related to LQCD calculations.",
5+
Creator -> "[email protected]",
6+
MathematicaVersion -> "11+",
7+
Extensions ->
8+
{
9+
{"Kernel", "Context" -> {
10+
"QCDLib`",
11+
"QCDLib`visualize`",
12+
"QCDLib`util`",
13+
"QCDLib`plot`",
14+
"QCDLib`analyze`",
15+
"QCDLib`sn`",
16+
"QCDLib`sun`",
17+
"QCDLib`phase`"
18+
}}
19+
}
20+
]

Diff for: ‎mathlib/QCDLib/QCDLib.m

+2
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,2 @@
1+
(* :: Package :: *)
2+
(* Does anything go here? *)

Diff for: ‎mathlib/QCDLib/analyze.m

+144
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,144 @@
1+
(* ::Package:: *)
2+
(* Data analysis functionality *)
3+
BeginPackage["QCDLib`analyze`"];
4+
5+
(* Unprotect and clear everything. *)
6+
Unprotect["`*"];
7+
ClearAll["`*"];
8+
9+
(* Define everything implicitly by giving usages. *)
10+
bootstrap::usage =
11+
"bootstrap[series, Nboot, f] = bootstrap estimate of sample \
12+
error of f on the series.";
13+
autocorrs::usage =
14+
"autocorrs[series, maxt] = evaluate samples of rho(t) for t in [0,maxt].";
15+
stdDev::usage =
16+
"stdDev[l] = standard dev of numerical ensemble, or zeros if length 1.";
17+
constantFit::usage =
18+
"constantFit[means, covars] = best fit constant, error, and chi squared.";
19+
constantSysFit::usage =
20+
"constantSysFit[means, covars, ti, tf] = best fit constant, stat + syst \
21+
error and chi squared given a window [ti,tf].";
22+
getMeffs::usage =
23+
"getMeffs[corrs] = the meffs and errors of a given ensemble \
24+
of correlators with dimensions Nmeas x Lt";
25+
meff::usage =
26+
"meff[corrs] = compute effective masses given time-sequence of correlators.";
27+
meffAcosh::usage =
28+
"meffAcosh[corrs] = compute acosh effective masses given correlators.";
29+
30+
Begin["`Private`"];
31+
(* Bootstrapping error analysis *)
32+
Options[bootstrap] = {"BinSize"->1, "RawBoot"->False, "Progress"->False};
33+
bootstrap[series_List, Nboot_Integer, f_, OptionsPattern[]] :=
34+
Module[{Nmeas, binSize, rawBoot, progress, binSeries, bootInds, bi, boots},
35+
Nmeas = Length[series];
36+
binSize = OptionValue["BinSize"];
37+
rawBoot = OptionValue["RawBoot"];
38+
progress = OptionValue["Progress"];
39+
(* Assert[Mod[Nmeas, binSize] == 0]; *)
40+
(* Binning *)
41+
If[binSize != 1,
42+
binSeries = Map[Mean, Partition[series, UpTo[binSize]]],
43+
binSeries = series];
44+
Nmeas = Length[binSeries];
45+
boots = {};
46+
If[progress, PrintTemporary[
47+
ProgressIndicator[Dynamic[bi], {1, Nboot+1}]]];
48+
For[bi = 1, bi <= Nboot, bi++,
49+
bootInds = RandomInteger[{1,Nmeas}, Nmeas];
50+
AppendTo[boots, f[binSeries[[bootInds]]]];
51+
];
52+
If[rawBoot, boots,
53+
{2 f[binSeries] - Mean[boots], StandardDeviation[boots]}]
54+
];
55+
56+
(* (normalized) autocorr samples, which can be fed into bootstrap mean
57+
for errors *)
58+
autocorrs[series_List, maxt_Integer] :=
59+
Module[{mean = Mean@series, out},
60+
out = Table[(series[[;;-t]] - mean) (series[[t;;]] - mean), {t, 1, maxt}];
61+
out / Mean[out[[1]]]
62+
];
63+
64+
(* Custom StdDev that handles single element lists *)
65+
stdDev[{element_}] := ConstantArray[0, Dimensions@element];
66+
stdDev[l_List] := StandardDeviation[l];
67+
68+
(* Constant value fitter, h/t mlwagman *)
69+
constantFit[means_, covars_] :=
70+
Module[{length, M, MSolver, MInv1, H, Delta, deltam, m, chiSq, dof, chiSqDof},
71+
length = Length@means;
72+
M = 1/2 (covars + ConjugateTranspose[covars]);
73+
MSolver = LinearSolve[M, Method->"Cholesky"];
74+
MInv1 = MSolver[ConstantArray[1,length]];
75+
If[Not@MatchQ[MInv1, _List], Throw["Ill-conditioned covariance matrix!"]];
76+
H = 2 ConstantArray[1,length].MInv1;
77+
Delta = 2 / H;
78+
deltam = Sqrt[Delta];
79+
m = Delta means . MInv1;
80+
chiSq = (means - m) . MSolver[means - m];
81+
dof = length - 1;
82+
chiSqDof = chiSq/dof;
83+
{m, deltam, chiSqDof}];
84+
85+
(* Constant value fitter with syst estimate, h/t mlwagman *)
86+
constantSysFit[means_, covars_, ti_, tf_] :=
87+
Module[{
88+
m1, \[Delta]m1, \[Chi]sq1, m2, \[Delta]m2, \[Chi]sq2,
89+
m3, \[Delta]m3, \[Chi]sq3, \[Delta]msys},
90+
91+
{m1,\[Delta]m1,\[Chi]sq1} = constantFit[
92+
means[[ti;;tf]], covars[[ti;;tf,ti;;tf]]];
93+
{m2,\[Delta]m2,\[Chi]sq2} = constantFit[
94+
means[[ti+1;;tf+1]], covars[[ti+1;;tf+1,ti+1;;tf+1]]];
95+
{m3,\[Delta]m3,\[Chi]sq3} = constantFit[
96+
means[[ti+2;;tf+2]], covars[[ti+2;;tf+2,ti+2;;tf+2]]];
97+
\[Delta]msys = 1/2 (Max[m1,m2,m3] - Min[m1,m2,m3]);
98+
99+
{m2, \[Delta]m2, \[Delta]msys, \[Chi]sq2, {ti,tf}}];
100+
101+
(* Effective masses from correlator data *)
102+
Options[getMeffs] = {"Nboot" -> 8, "BinSize" -> 1, "RawBoot" -> False};
103+
getMeffs[corrs_, OptionsPattern[]] :=
104+
Module[{
105+
Nmeas, Lt, Nboot, binSize, rawBoot, binNmeas, binCorrs,
106+
bi, bootInds, meffs, bootCorrs, bootMeffs},
107+
Nboot = OptionValue["Nboot"];
108+
binSize = OptionValue["BinSize"];
109+
rawBoot = OptionValue["RawBoot"];
110+
{Nmeas, Lt} = Dimensions@corrs;
111+
(* Binning *)
112+
Assert[Mod[Nmeas, binSize] == 0];
113+
binNmeas = Nmeas / binSize;
114+
binCorrs = Map[Mean, Partition[corrs, binSize]];
115+
(* Bootstrap *)
116+
meffs = {};
117+
For[bi = 1, bi <= Nboot, bi++,
118+
bootInds = RandomInteger[{1,binNmeas}, binNmeas];
119+
bootCorrs = Re@Mean@binCorrs[[bootInds]];
120+
bootMeffs = Log@bootCorrs[[;;-2]] - Log@bootCorrs[[2;;]];
121+
AppendTo[meffs, bootMeffs];
122+
];
123+
If[rawBoot,
124+
meffs,
125+
{Mean@meffs, stdDev@meffs}]];
126+
127+
(* Raw meff calculation on time-series *)
128+
meff[corrs_] := Module[{meanCorr},
129+
meanCorr = Mean@corrs;
130+
Re[Log[meanCorr[[;; -2]]] - Log[meanCorr[[2 ;;]]]]
131+
];
132+
133+
meffAcosh[corrs_] := Module[{meanCorr},
134+
meanCorr = Mean@corrs;
135+
ArcCosh[(meanCorr[[;; -3]] + meanCorr[[3 ;;]]) / (2 meanCorr[[2;;-2]])]
136+
];
137+
138+
End[]; (* `Private` *)
139+
140+
(* Log and protect everything exported. *)
141+
Print["Loaded " <> Context[] <> " with symbols: ", Names["`*"]];
142+
Protect["`*"];
143+
144+
EndPackage[];

Diff for: ‎mathlib/QCDLib/phase.m

+25
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
(* ::Package:: *)
2+
(* Functionality related to phase unwrapping. *)
3+
BeginPackage["QCDLib`phase`", {"QCDLib`util`"}];
4+
5+
(* Unprotect and clear everything. *)
6+
Unprotect["`*"];
7+
ClearAll["`*"];
8+
9+
(* Define everything implicitly by giving usages. *)
10+
get2DResids::usage =
11+
"get2DResids[phases] = array of resid value out of {-1,0,1} for each coord.";
12+
13+
Begin["`Private`"];
14+
get2DResids[phases_] := Round[(wrap[phases - RotateLeft[phases, {1, 0}]]
15+
+ wrap[RotateLeft[phases, {1, 0}] - RotateLeft[phases, {1, 1}]]
16+
+ wrap[RotateLeft[phases, {1, 1}] - RotateLeft[phases, {0, 1}]]
17+
+ wrap[RotateLeft[phases, {0, 1}] - phases]) / (2 Pi)];
18+
19+
End[]; (* `Private` *)
20+
21+
(* Log and protect everything exported. *)
22+
Print["Loaded " <> Context[] <> " with symbols: ", Names["`*"]];
23+
Protect["`*"];
24+
25+
EndPackage[];

Diff for: ‎mathlib/QCDLib/plot.m

+135
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,135 @@
1+
(* ::Package:: *)
2+
(* Plotting functionality *)
3+
BeginPackage["QCDLib`plot`", {"QCDLib`util`"}];
4+
5+
(* Unprotect and clear everything. *)
6+
Unprotect["`*"];
7+
ClearAll["`*"];
8+
9+
(* Define everything implicitly by giving usages. *)
10+
bigErrBarFn::usage =
11+
"bigErrBarFn[halfWidth] = function to plot error bar with end width 2d.";
12+
constFitEpilog::usage =
13+
"constFitEpilog[fitMean, fitErr, window] = list of epilog actions to plot constant line fit.";
14+
frameOpts::usage =
15+
"frameOpts[] = list of standard plot sizing and frame options.";
16+
title::usage =
17+
"title[str] = str in plot title format.";
18+
errorListLogPlot::usage =
19+
"errorListLogPlot[series] = log plot with error bars. Series should either be\
20+
a single list of {mean,err} tuples, or a list of series, each of which are a list of\
21+
{mean,err} tuples.";
22+
23+
24+
Begin["`Private`"];
25+
26+
Options[bigErrBarFn] = {"Directives" -> {}};
27+
bigErrBarFn[d_, OptionsPattern[]] := Function[{coords, errs},
28+
Module[{out},
29+
out = Join[OptionValue["Directives"], {
30+
Line[{coords + {0, errs[[2, 1]]},
31+
coords + {0, errs[[2, 2]]}}]}];
32+
If[d > 0,
33+
out = Join[out, {
34+
Line[{coords + {-d, errs[[2, 1]]}, coords + {d, errs[[2, 1]]}}],
35+
Line[{coords + {-d, errs[[2, 2]]}, coords + {d, errs[[2, 2]]}}]}]];
36+
out]];
37+
38+
constFitEpilog[fitMean_, fitErr_, {ti_,tf_}] := {
39+
Pink, Opacity[0.5], Rectangle[{ti, fitMean-fitErr}, {tf, fitMean+fitErr}],
40+
Red, Thick, Opacity[1.0], Line[{{ti, fitMean}, {tf, fitMean}}]
41+
};
42+
43+
frameOpts[] = {
44+
ImageSize -> 500, Frame -> True,
45+
FrameStyle -> Thick, FrameTicksStyle -> Directive[18]
46+
};
47+
48+
title[str_] := Style[str, Bold, Black, 18];
49+
50+
Options[errorListLogPlot] = {
51+
"LogPad" -> 0.05, "LogHardPad" -> 0.10,
52+
"YRange" -> Automatic,
53+
"Colors" -> mmaColors, "MeanJoined" -> False,
54+
"ErrorMarker" -> {Graphics[{Thick, Line[{{-0.5,0},{0.5,0}}]}], 0.03},
55+
"FillingStyle" -> {},
56+
"MainMarker" -> Automatic};
57+
errorListLogPlot[series_, opts:OptionsPattern[{errorListLogPlot, ListLogPlot}]] :=
58+
Module[{allseries, logplotmax, logplotmin, loghardmax, loghardmin, logrange,
59+
logpad, loghardpad, xpts, mpts, epts, pts, xmax, xmin, xspan, yrange, colors, extraPlotOpts},
60+
colors = OptionValue["Colors"];
61+
(* allseries has dims {nseries, npts, 2} or {nseries, npts, 3} *)
62+
allseries = If[Length@Dimensions@series == 2, {series}, series];
63+
allseries = Table[
64+
Table[
65+
If[Length[allseries[[i,j]]] == 3, allseries[[i,j]],
66+
Join[{j}, allseries[[i,j]]]],
67+
{j, Length[allseries[[i]]]}],
68+
{i, Length@allseries}];
69+
{xpts, mpts, epts} = Table[
70+
Table[allseries[[i,j,piece]], {j, Length[allseries[[i]]]}],
71+
{piece, 3}, {i, Length@allseries}];
72+
(* {xpts, mpts, epts} = Transpose[allseries, {2,3,1}]; *)
73+
logpad = OptionValue["LogPad"];
74+
loghardpad = OptionValue["LogHardPad"];
75+
logplotmax = Max@Abs@mpts;
76+
logplotmin = Min@Abs@mpts;
77+
logrange = Log[logplotmax] - Log[logplotmin];
78+
logplotmax *= Exp[logpad*logrange];
79+
logplotmin /= Exp[logpad*logrange];
80+
yrange = If[MatchQ[OptionValue["YRange"], Automatic],
81+
{logplotmin, logplotmax}, OptionValue["YRange"]];
82+
loghardmax = logplotmax * Exp[loghardpad * logrange];
83+
loghardmin = logplotmin / Exp[loghardpad * logrange];
84+
{xmin, xmax} = {Min[xpts], Max[xpts]};
85+
xspan = xmax-xmin;
86+
{xmin, xmax} = {xmin-xspan/10, xmax+xspan/10};
87+
ptsMean = Table[{mpts[[i]]}, {i, Length@mpts}];
88+
ptsErr = Table[{mpts[[i]] - epts[[i]],
89+
mpts[[i]] + epts[[i]]}, {i, Length@mpts}];
90+
ptsMean = Map[
91+
If[# > loghardmax, loghardmax,
92+
If[# < loghardmin, loghardmin, #]] &, ptsMean, {3}];
93+
ptsErr = Map[
94+
If[# > loghardmax, loghardmax,
95+
If[# < loghardmin, loghardmin, #]] &, ptsErr, {3}];
96+
colorsMean = colors[[;;Length@ptsMean]];
97+
colorsErr = Join@@Map[ConstantArray[#, 2]&, colorsMean];
98+
markersMean = OptionValue["MainMarker"];
99+
If[MatchQ[markersMean, Automatic],
100+
markersMean = mmaPlotMarkers[[;;Length@ptsMean]]];
101+
markersErr = ConstantArray[OptionValue["ErrorMarker"], 2 Length@ptsErr];
102+
fillingsErr = Table[2 i -> {2 i - 1}, {i, Length@ptsErr}];
103+
extraPlotOpts = getOptsFor[ListLogPlot, opts];
104+
extraPlotOptsNoLegend = FilterRules[extraPlotOpts, Except[PlotLegends]];
105+
If[OptionValue["MeanJoined"],
106+
AppendTo[extraPlotOpts, Joined->True]];
107+
tracesErr = Join@@ptsErr;
108+
tracesMean = Join@@ptsMean;
109+
tracesMean = Table[
110+
Transpose@{xpts[[i]], tracesMean[[i]]}, {i, Length@tracesMean}];
111+
tracesErr = Table[
112+
Transpose@{xpts[[Floor[(i+1)/2]]], tracesErr[[i]]}, {i, Length@tracesErr}];
113+
Show[{
114+
ListLogPlot[tracesErr,
115+
PlotRange -> {{xmin, xmax}, yrange},
116+
Filling -> fillingsErr, FillingStyle ->
117+
Directive@Join[{Opacity[1.0]}, OptionValue["FillingStyle"]],
118+
PlotStyle -> colorsErr, PlotMarkers -> markersErr,
119+
extraPlotOptsNoLegend],
120+
ListLogPlot[tracesMean,
121+
PlotRange -> {{xmin, xmax}, yrange},
122+
PlotStyle -> colorsMean, PlotMarkers -> markersMean,
123+
extraPlotOpts
124+
]
125+
}]
126+
];
127+
128+
End[]; (* `Private` *)
129+
130+
131+
(* Log and protect everything exported. *)
132+
Print["Loaded " <> Context[] <> " with symbols: ", Names["`*"]];
133+
Protect["`*"];
134+
135+
EndPackage[];

0 commit comments

Comments
 (0)