Select Git revision
-
Xisco Jimenez Forteza authoredXisco Jimenez Forteza authored
DataFits.m 111.47 KiB
(* ::Package:: *)
(************************************************************************)
(* This file was generated automatically by the Mathematica front end. *)
(* It contains Initialization cells from a Notebook file, which *)
(* typically will have the same name as this file except ending in *)
(* ".nb" instead of ".m". *)
(* *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs. Doing so is equivalent *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end. *)
(* *)
(* DO NOT EDIT THIS FILE. This entire file is regenerated *)
(* automatically each time the parent Notebook file is saved in the *)
(* Mathematica front end. Any changes you make to this file will be *)
(* overwritten. *)
(************************************************************************)
(* ::Code::Initialization:: *)
BeginPackage["DataFits`",{"NRLists`","ErrorBarPlots`"}];
(* ::Code::Initialization:: *)
\[Eta]::usage="\[Eta] for local usage";
S::usage="S for local usage";
\[Chi]1::usage="\[Chi]1 for local usage";
\[Chi]2::usage="\[Chi]2 for local usage";
\[Delta]\[Chi]::usage="\[Delta]\[Chi] for local usage";
a0::usage="a0 for local usage";
a1::usage="a1 for local usage";
a2::usage="a2 for local usage";
sTot::usage="";
sTot3::usage="";
\[Chi]Tot::usage="";
\[Chi]diffplus::usage="";
\[Chi]diffstan::usage="";
DataFitFunction::usage="DataFitFunction[data_?ListQ,ansatzList_?ListQ,OptionsPattern[]]";
DataFitFunctionAll::usage="DataFitFunctionAll[dataraw_?ListQ,OptionsPattern[]]";
DataFitFunctionAllNoWeights::usage="DataFitFunctionAllWeights[dataraw_?ListQ,OptionsPattern[]]";
AnsatzRestrictions::usage="AnsatzRestrictions[nsfit_, q1fit_, ansatzGen_]";
Plot2DFits::usage="Plot2DFits[data_?ListQ,fitlist_?ListQ,fitvars_?ListQ, OptionsPattern[]]";
myListPlot3D::usage="myListPlot3D[data_?ListQ,options]: Plot 3d data + interpolated function";
ColorGradient::usage="ColorGradient[data,colors,options]. Create a linear gradient of colors for the data";
CreateColors::usage="CreateColors[codes,colors]. Set different colors for each NR code";
CleanAnsatzParams::usage="cleanAnsatzParams[paramsGuess_, vars_]";
Generate1DPolynomialAnsatz::usage="Generate1DPolynomialAnsatz[CoefficientPrefixString_?StringQ,variable_,MinOrder_?IntegerQ,MaxOrder_?IntegerQ] creates a polynomial function ansatz.";
Residuals::usage="Residuals[data_,fit_,vars_] computes fit residuals.";
AtomsList::usage="Take the coefficients out";
GeneralizeFunction::usage="GeneralizeFunction[expr_,x_] insert free coefficient at all (non-exact) real numbers in a function";
GeneralizeFunction::usage="GeneralizeFunction[expr_,x_,coord_,orderNum,orderDenom_] insert polynomial with free coefficients at all (non-exact) real numbers in a function";
ExactPade::usage="ExactPade[ansatz_,coeffRules_,padeOptions_] get Pade approximant while keeping exact coefficients fixed.";
Generate1DPadeAnsatzList::usage="Generate1DPadeAnsatzList[ansatz_,\!\(\*
StyleBox[\"coeffRules_\", \"Code\"]\)\!\(\*
StyleBox[\",\", \"Code\"]\)\!\(\*
StyleBox[\"variable_\", \"Code\"]\)\!\(\*StyleBox[\",\", \"Code\"]\)\!\(\*
StyleBox[\"paramName_\", \"Code\"]\)\!\(\*
StyleBox[\",\", \"Code\"]\)\!\(\*
StyleBox[\"minPadeTypeSum_\", \"Code\"]\)\!\(\*
StyleBox[\",\", \"Code\"]\)\!\(\*
StyleBox[\"maxPadeTypeSum_\", \"Code\"]\)] make a list of all Pade ansaetze with type (num+denom) summing up between [min,max].";
InvertRules::usage="InvertRules[rules_] invert a list of rules";
Get1dInverseRules::usage="Get1dInverseRules[ansatz_,coeffRules_] invert a set of 1d rules in eta or S, taking care of the Pade form";
Get1dInverseRulesS::usage="Get1dInverseRulesS[ansatz_,coeffRules_,productAnsatz_] invert the 1d rules in S, taking care of normalization and Pade form";
FitPredictionIntervalFunctionFinalOnly::usage="FitPredictionIntervalFunctionFinalOnly[fit_,ansatzRules_] return function for the fit uncertainty (prediction interval), final fit statistics only";
FitPredictionIntervalFunctionFinalOnlyq1::usage="FitPredictionIntervalFunctionFinalOnlyq1[fit_,ansatzRules_] return function for the fit uncertainty (prediction interval), in the limit of q=1, final fit statistics only";
FitPredictionIntervalFinalOnly::usage="FitPredictionIntervalFinalOnly[fit_,ansatzRules_,etain_,chi1in_,chi2in_] estimate the fit uncertainty (prediction interval) at a given point, final fit statistics only";
FitPredictionIntervalStderrSq::usage="FitPredictionIntervalStderrSq[\!\(\*
StyleBox[\"finalFit_\", \"Code\"]\)\!\(\*
StyleBox[\",\", \"Code\"]\)\!\(\*
StyleBox[\"finalAnsatzRules_\", \"Code\"]\)\!\(\*
StyleBox[\",\", \"Code\"]\)\!\(\*
StyleBox[\"fit2dParts_\", \"Code\"]\)\!\(\*
StyleBox[\",\", \"Code\"]\)\[Eta]0stuff_,extraCoeffRules_,productAnsatz_,q1_] evaluate the various stderrsq contributions for fit uncertainty intervals";
FitPredictionIntervalFunction::usage="FitPredictionIntervalFunction[finalFit_,finalAnsatzRules_,fit2dParts_,\[Eta]0stuff_,extraCoeffRules_,productAnsatz_] return function for the fit uncertainty (prediction interval)";
FitPredictionIntervalFunctionq1::usage="FitPredictionIntervalFunctionq1[finalFit_,finalAnsatzRules_,fit2dParts_,\[Eta]0stuff_,extraCoeffRules_,productAnsatz_] return function for the fit uncertainty (prediction interval), in the limit of q=1";
FitPredictionInterval::usage="FitPredictionInterval[finalFit_,finalAnsatzRules_,fit2dParts_,\[Eta]0fit_,extraCoeffRules_,productAnsatz_,etain_,chi1in_,chi2in_] estimate the fit uncertainty (prediction interval) at a given point";
TeXExportCoeffTable::usage="TeXExportCoeffTable[fit_,filename_] export a fit coefficient table into a nice TeX file.";
TeXExportCovarMatrixTable::usage="TeXExportCovarMatrixTable[fitCovar_,fitCoeffRules_,filename_] export fit covariance matrix into a nice TeX file.";
SubscriptRules::usage="SubscriptRules[rules_] make a subscript rule from a numeric coefficient rule.";
GetRawTwoDAnsatz::usage="GetRawTwoDAnsatz[finalfit_,fit2dParts_,productAnsatz_,constrained_] put together raw 2d ansatz equation.";
GetAllRules::usage="GetAllRules[finalfit_,fit2dParts_,twodrules_] collect all 1d, 2d and 3d coefficient rules";
TeXFormatAnsatz::usage="TeXFormatAnsatz[ansatz_,formattingRules_,coeffRules_] format an ansatz nicer in TeX";
TeXExportAnsatz::usage="TeXExportAnsatz[ansatz_,formattingRules_,coeffRules_,filename_] export a nicely formatted ansatz to TeX.";
TeXExportTwoDAnsatz::usage="TeXExportTwoDAnsatz[finalfit_,fit2dParts_,twodrules_,formattingRules_,productAnsatz_,constrained_,filename_] format the 2d part of the final ansatz.";
TeXExportChiDiffTerms::usage="TeXExportChiDiffTerms[chidiffAnsaetze_,formattingRules_,coeffRules_,filename_] concatenate and export to TeX the various chidiff ansatz terms.";
TeXExportFinalAnsatz::usage="TeXExportFinalAnsatz[finalfit_,fit2dParts_,twodrules_,chidiffAnsaetze_,formattingRules_,productAnsatz_,twodConstrained_,filename_] format the final ansatz nicely.";
TexExportCoeffTable::usage="TeXExportCoeffTable[coeffRules_,covar_,filename_] format a fit coefficient table without t-stat nor p-value.";
TeXExportTabularTable::usage="TeXExportTabularTable[datatable_,filename_,rowHeadings_,colHeadings_,padZeroes_] export a table with array->tabular and formatting fixes.";
PyExportFinalFit::usage="PyExportFinalFit[fit_,extraFormattingRules_,filename_] export the final fit (with numerical coefficients) for python LAL uise.";
PyExportFinalAnsatz::usage="PyExportFinalAnsatz[finalfit_,fit2dParts_,chidiffAnsaetze_,productAnsatz_,extraFormattingRules_,filename_] export the final ansatz (with symbolic coefficients) for python LAL uise.";
PyExportFinalFitCoeffs::usage="PyExportFinalFitCoeffs[finalfit_,fit2dParts_,\!\(\*
StyleBox[\"all2dconstraints_\", \"Code\"]\)\!\(\*
StyleBox[\",\", \"Code\"]\)filename_] export the numerical coefficients to be used together with the ansatz.";
SupplExportAllFitCoeffs::usage="SupplExportAllFitCoeffs[finalfit_,fit2dParts_,\!\(\*
StyleBox[\"all2dconstraints_\", \"Code\"]\)\!\(\*
StyleBox[\",\", \"Code\"]\)\!\(\*
StyleBox[\"eta0covar_\", \"Code\"]\)\!\(\*
StyleBox[\",\", \"Code\"]\)filename_] export the numerical coefficients to a plain ASCII table file.";
SupplExportCovarMatrix::usage="SupplExportCovarMatrix[covar_,coeffRules_,fileName_] export covariance matrices as simple ASCII table files.";
AIC::usage="AIC[data,fitname,variables,number of parameters]: Explicit computation of the AIC value given by NonlinearModelFit. In data [[#,-2]]: values, [[#,-1]]: weights. Fit is the fit name e.g FinalSpin0815. Does not work for RITFinalSpinNonPrec2014";
AICc::usage="AIC[data,fitname,variables,number of parameters]: Explicit computation of the AICc value given by NonlinearModelFit. In data [[#,-2]]: values, [[#,-1]]: weights. Fit is the fit name e.g FinalSpin0815. Does not work for RITFinalSpinNonPrec2014";
BIC::usage="AIC[data,fitname,variables,number of parameters]: Explicit computation of the BIC value given by NonlinearModelFit. In data [[#,-2]]: values, [[#,-1]]: weights. Fit is the fit name e.g FinalSpin0815. Does not work for RITFinalSpinNonPrec2014";
KullbagLeiblerDiv::usage="KullbagLeiblerDiv[p1_?ListQ,p2_?ListQ]: Computation of the KL criterion";
JensenShanonDiv::usage="JensenShanonDiv[p1_?ListQ,p2_?ListQ]: Computation of the JS criterion";
KL::usage="KL[p1_?ListQ,p2_?ListQ],lim: Computation of the KL criterion";
JS::usage="JS[p1_?ListQ,p2_?ListQ,lim]: Computation of the JS criterion";
CredibleRegion::usage="CredibleRegion[data_,level_]: Computation of the credible region";
ComputeEdges::usage="ComputeEdges[pts]";
CredibleInterval::usage="CredibleInterval[data_,level_]: Computation of the credible intervals";
(* ::Code::Initialization:: *)
Begin["`Private`"];
(* ::Code::Initialization:: *)AtomsList[expr_]:=Union@Select[Level[expr,{0,Infinity}],AtomQ];
InterpolationDomain[fun_]:=Module[{min,max},{min,max}={fun[[1,1,1]],fun[[1,1,2]]}];
q\[Eta][\[Eta]_]:=-(1.-1/(2\[Eta]))+Sqrt[(1-1/(2\[Eta]))^2 -1]
\[Eta]q[q_]:=q/(1.+q)^2.
\[Chi]diffplus[\[Eta]_,xx_,yy_]:=Module[{s},s=(1/2(1+Sqrt[1-4 \[Eta]])xx - 1/2(1-Sqrt[1-4 \[Eta]])yy)]
\[Chi]diffstan[\[Eta]_,xx_,yy_]:=Module[{s},s=xx-yy]
myListPlot3D[list_,opt___]:=Module[{p1,p2},
p1=ListPlot3D[If[(Length@Dimensions@list)>1,list[[1]],list],InterpolationOrder-> 3];
p2=ListPointPlot3D[list,PlotStyle-> PointSize-> Large,opt];
Show[p1,p2]
];
myListPlot3D[list_,opt___]:=Module[{p1,p2},
p1=ListPlot3D[list,InterpolationOrder-> 3,opt];
p2=ListPointPlot3D[list,PlotStyle-> PointSize-> Large,opt];
Show[p1,p2]
];
\[Chi]Tot[\[Eta]_,\[Chi]1_,\[Chi]2_]:=Module[{m1,m2,s},
m1=1/2 (1+Sqrt[1-4 \[Eta]]);
m2=1/2 (1-Sqrt[1-4 \[Eta]]);
s=(m1 \[Chi]1 + m2 \[Chi]2)/(m1+ m2)
]
sTot[\[Eta]_,\[Chi]1_,\[Chi]2_]:=Module[{m1,m2},
m1=1/2 (1+Sqrt[1-4 \[Eta]]);
m2=1/2 (1-Sqrt[1-4 \[Eta]]);
(m1^2 \[Chi]1 + m2^2 \[Chi]2)
]
sTotR[\[Eta]_,\[Chi]1_,\[Chi]2_]:=Module[{m1,m2},
m1=1/2 (1+Sqrt[1-4 \[Eta]]);
m2=1/2 (1-Sqrt[1-4 \[Eta]]);
(m1^2 \[Chi]1 + m2^2 \[Chi]2)/(m1^2 + m2^2)
]
sTot3[\[Eta]_,\[Chi]1_,\[Chi]2_]:=Module[{m1,m2},
m1=1/2 (1+Sqrt[1-4 \[Eta]]);
m2=1/2 (1-Sqrt[1-4 \[Eta]]);
(m1^2 \[Chi]1 + m2^2 \[Chi]2)/(m1^2 + m2^2)
]
countSummands[expr_]:=If[Head[expr]===Plus,Length[expr],If[expr===0,0,1]]
GeneralizeFunction[expr_,x_]:=Module[{num,numTermsNum,numSymbolsFirstTerm,i,termsNum,termsDenom,iStart,rule0,rulesNum,rulesDenom,reals},
num=Numerator[expr];
numTermsNum = countSummands[num];
If[numTermsNum==1,
termsNum = {num};
numSymbolsFirstTerm=1;
,
termsNum = Level[num,{1,1}];
numSymbolsFirstTerm=Length[Cases[Level[termsNum[[1]],{1,Infinity}],y_Symbol]];
];
If[numSymbolsFirstTerm==0,
iStart=2;
reals=Cases[Tally[Select[Level[termsNum[[1]],{0,Infinity}],AtomQ]][[All,1]],y_Real];
If[Length[reals]>0,
rule0={termsNum[[1]]->termsNum[[1]] ToExpression[ToString@x<>"0"]};
,
rule0={};
];
,
iStart=1;
rule0={};
];
rulesNum=rule0;
For[i=iStart,i<=Length[termsNum],i++,
reals=Cases[Tally[Select[Level[termsNum[[i]],{0,Infinity}],AtomQ]][[All,1]],y_Real];
If[Length[reals]>0,
rulesNum=Join[rulesNum,{termsNum[[i]] -> termsNum[[i]] ToExpression[ToString@x<>ToString[i-iStart+1]]}];
];
];
(* currently assumes that denominator always has a -1 in front, can be easily generalized with countSummands as for numerator *)
termsDenom = Level[Denominator[expr],{1,1}];
rulesDenom = {};
For[i=1,i<=Length[termsDenom],i++,
reals=Cases[Tally[Select[Level[termsDenom[[i]],{0,Infinity}],AtomQ]][[All,1]],y_Real];
If[Length[reals]>0,
rulesDenom=Join[rulesDenom,{termsDenom[[i]] -> termsDenom[[i]] ToExpression[ToString@x<>ToString[Length[termsNum]-iStart+1+i]]}];
];
];
Return[(Numerator[expr]/.rulesNum)/(Denominator[expr]/.rulesDenom)]
]
GeneralizeFunction[expr_,x_,coord_,orderNum_,orderDenom_]:=Module[{num,numTermsNum,numSymbolsFirstTerm,i,termsNum,termsDenom,iStart,rule0,rulesNum,rulesDenom,reals},
num=Numerator[expr];
numTermsNum = countSummands[num];
If[numTermsNum==1,
termsNum = {num};
numSymbolsFirstTerm=1;
,
termsNum = Level[num,{1,1}];
numSymbolsFirstTerm=Length[Cases[Level[termsNum[[1]],{1,Infinity}],y_Symbol]];
];
If[numSymbolsFirstTerm==0,
iStart=2;
reals=Cases[Tally[Select[Level[termsNum[[1]],{0,Infinity}],AtomQ]][[All,1]],y_Real];
If[Length[reals]>0,
rule0={termsNum[[1]]->termsNum[[1]](Sum[coord^j ToExpression/@ToExpression@(ToString@x<>"0"<>ToString@j),{j,0,orderNum}])};
,
rule0={};
];
,
iStart=1;
rule0={};
];rulesNum=rule0;
For[i=iStart,i<=Length[termsNum],i++,
reals=Cases[Tally[Select[Level[termsNum[[i]],{0,Infinity}],AtomQ]][[All,1]],y_Real];
If[Length[reals]>0,
rulesNum=Join[rulesNum,{termsNum[[i]] -> termsNum[[i]](Sum[coord^j ToExpression@(ToString@x<>ToString[i-iStart+1]<>ToString@j),{j,0,orderNum}])}];
];
];
(* currently assumes that denominator always has a -1 in front, can be easily generalized with countSummands as for numerator *)
termsDenom = Level[Denominator[expr],{1,1}];
rulesDenom = {};
For[i=1,i<=Length[termsDenom],i++,
reals=Cases[Tally[Select[Level[termsDenom[[i]],{0,Infinity}],AtomQ]][[All,1]],y_Real];
If[Length[reals]>0,
rulesDenom=Join[rulesDenom,{termsDenom[[i]] -> termsDenom[[i]](Sum[coord^j ToExpression/@ToExpression@(ToString@x<>ToString[Length[termsNum]-iStart+1+i]<>ToString@j),{j,0,orderDenom}])}];
];
];
(Numerator[expr]/.rulesNum)/(Denominator[expr]/.rulesDenom)
]
(*wrapper to Pade approximant that keeps exact coefficients like Sqrt[2] exact, so that GeneralizeFunction won't assign extra coefficients to them.*)
ExactPade[ansatz_,coeffRules_,padeOptions_]:=PadeApproximant[ansatz,padeOptions]/.coeffRules;
myRound[x_,d_] := Module[{mantExp,mant,exp},
mantExp = MantissaExponent[x];
mant = mantExp[[1]];
exp = mantExp[[2]];
If[mant<1.0,
mant = 10*mant;
exp = exp-1;
];
mant = Round[mant,0.1^d];
(*
If[exp<0,
mant = Round[mant,0.1^(d-1)];
,
mant = Round[mant,0.1^d];
];
*)
Return[mant*10^exp]
]
Generate1DPadeAnsatzList[ansatz_,coeffRules_,variable_,paramName_,minPadeTypeSum_,maxPadeTypeSum_]:=Module[{i,j,ansatzDegree,padeAnsatzList,padeType,exactPade,genPade,genPadeList,exactPadeReals,exactPadeRounding},
For[i=1,i<maxPadeTypeSum,i++,
j=Max[1,minPadeTypeSum-i];
While[(j<=i)&&(i+j<=maxPadeTypeSum),
ansatzDegree = Exponent[ansatz,variable];
If[(j==1)&&(i==ansatzDegree),
Print["Not generating Pade approximant for type ("<>ToString[i]<>","<>ToString[j]<>") with polynomial ansatz of order "<>ToString[ansatzDegree]<>", as it would just duplicate the polynomial."];
,
padeType = {i,j};
exactPade = ExactPade[ansatz,coeffRules,{variable,0,padeType}];
exactPadeReals=Cases[AtomsList[exactPade],y_Real];
exactPadeRounding=Table[exactPadeReals[[i]]->myRound[exactPadeReals[[i]],2],{i,1,Length@exactPadeReals}];
exactPade = exactPade /. exactPadeRounding;
genPade = GeneralizeFunction[exactPade,paramName];
If[!ValueQ[genPadeList],
genPadeList={{genPade,{variable}}};
,
AppendTo[genPadeList,{genPade,{variable}}];
]
];
j++;
];
];Print["Generated these Pade approximants:"];
Print[genPadeList];
Return[genPadeList];
]
InvertRules[rules_]:=Module[{},
Table[rules[[i,2]]->rules[[i,1]],{i,Length@rules}]
]
Get1dInverseRules[ansatz_,coeffRules_]:=Module[{origCoeffs,inverseRules},
(* in case this is a Pade rational ansatz, disentangle it *)
origCoeffs = Cases[Select[Join[Level[Numerator[ansatz],{1,Infinity}],Level[Denominator[ansatz],{1,Infinity}]],AtomQ],y_Real];
If[Length[origCoeffs]>0,
inverseRules = Table[origCoeffs[[i]]*coeffRules[[i,2]] -> origCoeffs[[i]]*coeffRules[[i,1]], {i,Length@coeffRules}];
,
inverseRules = InvertRules[coeffRules];
];
Return[inverseRules];
]
Get1dInverseRulesS[ansatz_,coeffRules_,productAnsatz_]:=Module[{ansatzWithoutS0,inverseRules},
If[productAnsatz,
ansatzWithoutS0 = ansatz/(ansatz/.S->0)/.{1.->1};
,
ansatzWithoutS0 = ansatz-(ansatz/.S->0)/.{0.->0};
];
inverseRules = Get1dInverseRules[ansatzWithoutS0,coeffRules];
Return[inverseRules];
]
ColorGradient[data_,ColorList_,OptionsPattern[{"Weights"->"","Verbose"->False}]]:=Module[{min,max,weights,wc,blendtab,\[Delta]n,verbose},
weights=OptionValue["Weights"];
verbose=OptionValue["Verbose"];
min=Min[data];
max=Max[data];
\[Delta]n=(max-min)1./(Length@ColorList-1);
If[Not@ListQ@weights,weights=Table[min + \[Delta]n(i-1),{i,Length@ColorList}];,weights=weights Table[min + \[Delta]n(i-1),{i,Length@ColorList}];];
blendtab=Table[{weights[[i]],ColorList[[i]]},{i,Length@ColorList}];
If[verbose,Print["Weights: ",weights]];
Table[Blend[blendtab,data[[i]]],{i,Length@data}]
]
CreateColors[code_,colors_]:=Module[{},
Which[code=="BAM",colors[[1]],code=="SXS",colors[[2]],code=="GaTech",colors[[3]],code=="RIT",colors[[4]],True,colors[[5]]]]
FitPredictionIntervalFunctionFinalOnly[fit_,ansatzRules_]:=Module[{ansatz,coeffNames,coeffRules,Ncoeff,Ndata,EstVar,quant95,covarMatrix,coeffGrad,stderrsq},
ansatz = fit[[15]];
coeffRules = fit[[2]];
coeffNames = coeffRules[[All,1]];
Ncoeff = Length@coeffNames;
Ndata = Length@Last@fit; (* length of residuals vector *)EstVar = fit[[16]]; (* would be = Total[resid^2]/(Ndata-Ncoeff) without weights *)
(* take gradient vector in the coefficients *)
coeffGrad = Table[D[ansatz/.ansatzRules,coeffNames[[i]]],{i,Ncoeff}] /. coeffRules;
(* estimate of fit error: multiply coefficient gradient with covariance matrix *)
covarMatrix = fit[[14]];
stderrsq = (coeffGrad.covarMatrix.coeffGrad);
(* report back the 95% student-t quantile (applied on both sides, this gives a 90% interval) *)
quant95 = Quantile[StudentTDistribution[Ndata-Ncoeff],0.95];
Return[quant95 * Sqrt[ EstVar + stderrsq ]];
]
FitPredictionIntervalFunctionFinalOnlyq1[fit_,ansatzRules_]:=Module[{ansatz,coeffNames,coeffRules,Ncoeff,Ndata,EstVar,quant95,covarMatrix,coeffGrad,stderrsq,chidiffcoeff,chidiff2coeff},
ansatz = fit[[15]];
coeffRules = fit[[2]];
coeffNames = coeffRules[[All,1]];
Ncoeff = Length@coeffNames;
Ndata = Length@Last@fit; (* length of residuals vector *)
EstVar = fit[[16]]; (* would be = Total[resid^2]/(Ndata-Ncoeff) without weights *)
(* take gradient vector in the coefficients *)
coeffGrad = Table[D[ansatz/.ansatzRules,coeffNames[[i]]],{i,Ncoeff}] /. coeffRules;
(* take care of issues in the q=1, eta=0.25 limit of the spin-diff terms *)
chidiffcoeff = Coefficient[coeffGrad,(\[Chi]1-\[Chi]2)];
chidiff2coeff = Coefficient[coeffGrad,(\[Chi]1-\[Chi]2)^2];
coeffGrad = coeffGrad - chidiffcoeff(\[Chi]1-\[Chi]2) + Limit[chidiffcoeff,\[Eta]->0.25](\[Chi]1-\[Chi]2) - chidiff2coeff (\[Chi]1-\[Chi]2)^2 + Limit[chidiff2coeff,\[Eta]->0.25](\[Chi]1-\[Chi]2)^2;
(* estimate of fit error: multiply coefficient gradient with covariance matrix *)
covarMatrix = fit[[14]];
stderrsq = (coeffGrad.covarMatrix.coeffGrad);
(* report back the 95% student-t quantile (applied on both sides, this gives a 90% interval) *)
quant95 = Quantile[StudentTDistribution[Ndata-Ncoeff],0.95];
Return[quant95 * Sqrt[ EstVar + stderrsq ]];
]
FitPredictionIntervalFinalOnly[fit_,ansatzRules_,etain_,chi1in_,chi2in_]:=Module[{fiterrFunc,fiterrFuncq1,fiterrs,i,eta,chi1,chi2},
(* so we can deal with both scalars and lists *)
If[Length[etain]==0,
eta = {etain};
chi1 = {chi1in};
chi2 = {chi2in};
,
eta = etain;
chi1 = chi1in;
chi2 = chi2in;
];
(* evaluate the general error estimate *)
fiterrFunc = FitPredictionIntervalFunctionFinalOnly[fit,ansatzRules];
(* avoid indeterminates at q=1, eta=0.25 *)
If[MemberQ[eta,0.25],
fiterrFuncq1 = FitPredictionIntervalFunctionFinalOnlyq1[fit,ansatzRules];
];
fiterrs = Table[fiterrFunc,{i,1,Length@eta}];
For[i=1,i<=Length@eta,i++,
If[eta[[i]]==0.25,
fiterrs[[i]]=fiterrFuncq1;
];
];
(* switch spin parametrization and insert user values *)fiterrs = fiterrs/.{S->sTot3[\[Eta],\[Chi]1,\[Chi]2]};
Return[Table[fiterrs[[i]]/.{\[Eta]->eta[[i]],\[Chi]1->chi1[[i]],\[Chi]2->chi2[[i]]},{i,1,Length@eta}]];
]
FitPredictionIntervalStderrSq[finalFit_,finalAnsatzRules_,fit2dParts_,\[Eta]0stuff_,extraCoeffRules_,productAnsatz_,q1_]:=Module[{finalAnsatz,finalCoeffNames,finalCoeffRules,NfinalCoeffs,finalCovarMatrix,finalCoeffGrad,chidiffcoeff,chidiff2coeff,
etaCoeffRules,etaCoeffNames,NetaCoeffs,etaAnsatz,etaInverseRules,etaCoeffGrad,etaCovarMatrix,
SCoeffRules,SCoeffNames,NSCoeffs,SAnsatz,SInverseRules,SCoeffGrad,SCovarMatrix,
twodCoeffRules,ansatz2dRaw,zeroRules,
\[Eta]0fit,\[Eta]0Constraints,\[Eta]0CoeffRules,\[Eta]0CoeffNames,N\[Eta]0Coeffs,\[Eta]0CovarMatrix,\[Eta]0Coeffgrad,
finalStderrSq,etaStderrSq,SStderrSq,twodStderrSq,\[Eta]0StderrSq,
allStderrSq},
finalAnsatz = finalFit[[15]];
finalCoeffRules = finalFit[[2]];
finalCoeffNames = finalCoeffRules[[All,1]];
NfinalCoeffs = Length@finalCoeffNames;
(* take gradient vector in the coefficients *)
finalCoeffGrad = Table[D[finalAnsatz/.finalAnsatzRules,finalCoeffNames[[i]]],{i,NfinalCoeffs}] /. finalCoeffRules;
(* take care of issues in the q=1, eta=0.25 limit of the spin-diff terms *)
If[q1,
chidiffcoeff = Coefficient[finalCoeffGrad,(\[Chi]1-\[Chi]2)];
chidiff2coeff = Coefficient[finalCoeffGrad,(\[Chi]1-\[Chi]2)^2];
finalCoeffGrad = finalCoeffGrad - chidiffcoeff(\[Chi]1-\[Chi]2) + Limit[chidiffcoeff,\[Eta]->0.25](\[Chi]1-\[Chi]2) - chidiff2coeff (\[Chi]1-\[Chi]2)^2 + Limit[chidiff2coeff,\[Eta]->0.25](\[Chi]1-\[Chi]2)^2;
finalCoeffGrad = Limit[finalCoeffGrad,\[Eta]->0.25];
];
(* estimate of fit error: multiply coefficient gradient with covariance matrix *)
finalCovarMatrix = finalFit[[14]];
finalStderrSq = (finalCoeffGrad.finalCovarMatrix.finalCoeffGrad);
(* also get contributions from 1D eta and S fits *)
etaAnsatz = fit2dParts[[2]];
etaCoeffRules = fit2dParts[[3]];
etaCoeffNames = etaCoeffRules[[All,1]];
NetaCoeffs = Length@etaCoeffNames;
etaCovarMatrix = fit2dParts[[4]];
SAnsatz = fit2dParts[[5]];
SCoeffRules = fit2dParts[[6]];
SCoeffNames = SCoeffRules[[All,1]];
NSCoeffs = Length@SCoeffNames;
SCovarMatrix = fit2dParts[[7]];
(* construct raw 2D ansatz *)
etaInverseRules = Get1dInverseRules[etaAnsatz,etaCoeffRules];
SInverseRules = Get1dInverseRulesS[SAnsatz,SCoeffRules,productAnsatz];
twodCoeffRules = Join[fit2dParts[[8]],fit2dParts[[9]]];
If[Not@productAnsatz,
zeroRules = Table[ToExpression["f"<>ToString@i<>"0"]->0,{i,0,Exponent[Numerator[fit2dParts[[1]]],S]}];
twodCoeffRules = Join[twodCoeffRules/.zeroRules,zeroRules];
];
ansatz2dRaw = fit2dParts[[1]]/.etaInverseRules/.SInverseRules;
(* take gradient vectors of the final ansatz (only 2D part needed) in the previously-determined coefficients *)
etaCoeffGrad = Table[D[ansatz2dRaw, etaCoeffNames[[i]]], {i,NetaCoeffs}] /. Join[etaCoeffRules,SCoeffRules,twodCoeffRules];
SCoeffGrad = Table[D[ansatz2dRaw, SCoeffNames[[i]]], {i,NSCoeffs}] /. Join[etaCoeffRules,SCoeffRules,twodCoeffRules];
(* also get contribution from extreme-mass ratio limit fit *)
If[ListQ[\[Eta]0stuff]&&(Length[\[Eta]0stuff]==0),
\[Eta]0StderrSq=0;
etaCoeffGrad = etaCoeffGrad /. Join[finalCoeffRules,extraCoeffRules];
SCoeffGrad = SCoeffGrad /. Join[finalCoeffRules,extraCoeffRules];
,
If[ListQ[\[Eta]0stuff],
\[Eta]0fit = \[Eta]0stuff[[1]];
\[Eta]0Constraints = \[Eta]0stuff[[2]]; ,
\[Eta]0fit = \[Eta]0stuff;
\[Eta]0Constraints = {};
];
\[Eta]0CoeffRules = \[Eta]0fit["BestFitParameters"];
\[Eta]0CoeffNames = \[Eta]0CoeffRules[[All,1]];
N\[Eta]0Coeffs = Length@\[Eta]0CoeffNames;
\[Eta]0CovarMatrix = \[Eta]0fit["CovarianceMatrix"];
\[Eta]0Coeffgrad = Table[D[ansatz2dRaw/.twodCoeffRules/.\[Eta]0Constraints, \[Eta]0CoeffNames[[i]]], {i,N\[Eta]0Coeffs}] /. Join[etaCoeffRules,SCoeffRules] /. Join[\[Eta]0CoeffRules,finalCoeffRules,extraCoeffRules];
If[q1,
\[Eta]0Coeffgrad = Limit[\[Eta]0Coeffgrad, \[Eta]->0.25];
];
\[Eta]0StderrSq = (\[Eta]0Coeffgrad.\[Eta]0CovarMatrix.\[Eta]0Coeffgrad);
etaCoeffGrad = etaCoeffGrad /. \[Eta]0Constraints /. Join[\[Eta]0CoeffRules,finalCoeffRules,extraCoeffRules];
SCoeffGrad = SCoeffGrad /. \[Eta]0Constraints /. Join[\[Eta]0CoeffRules,finalCoeffRules,extraCoeffRules];
];
(* dot gradient vectors together with the corresponding covar matrices *)
etaStderrSq = (etaCoeffGrad.etaCovarMatrix.etaCoeffGrad);
SStderrSq = (SCoeffGrad.SCovarMatrix.SCoeffGrad);
If[q1,
etaCoeffGrad = Limit[etaCoeffGrad, \[Eta]->0.25];
SCoeffGrad = Limit[SCoeffGrad, \[Eta]->0.25];
];
allStderrSq = {finalStderrSq,etaStderrSq,SStderrSq,\[Eta]0StderrSq};
Return[allStderrSq];
]
FitPredictionIntervalFunction[finalFit_,finalAnsatzRules_,fit2dParts_,\[Eta]0fit_,extraCoeffRules_,productAnsatz_]:=Module[{finalCoeffRules,finalCoeffNames,NfinalCoeffs,Ndata,EstVar,quant95,stderrsq,stderrSqContribs,q1},
finalCoeffRules = finalFit[[2]];
finalCoeffNames = finalCoeffRules[[All,1]];
NfinalCoeffs = Length@finalCoeffNames;
Ndata = Length@Last@finalFit; (* length of residuals vector *)
EstVar = finalFit[[16]]; (* would be = Total[resid^2]/(Ndata-Ncoeff) without weights *)
q1=False;
stderrSqContribs = FitPredictionIntervalStderrSq[finalFit,finalAnsatzRules,fit2dParts,\[Eta]0fit,extraCoeffRules,productAnsatz,q1];
(* add up all error contributions in quadrature *)
stderrsq = Total[stderrSqContribs];
(* report back the 95% student-t quantile (applied on both sides, this gives a 90% interval), adding up also the variance contribution for PREDICTION interval (not mean confidence) *)
quant95 = Quantile[StudentTDistribution[Ndata-NfinalCoeffs],0.95];
Return[quant95 * Sqrt[ EstVar + stderrsq ]];
]
FitPredictionIntervalFunctionq1[finalFit_,finalAnsatzRules_,fit2dParts_,\[Eta]0stuff_,extraCoeffRules_,productAnsatz_]:=Module[{finalCoeffRules,finalCoeffNames,NfinalCoeffs,Ndata,EstVar,quant95,stderrsq,stderrSqContribs,q1},
finalCoeffRules = finalFit[[2]];
finalCoeffNames = finalCoeffRules[[All,1]];
NfinalCoeffs = Length@finalCoeffNames;
Ndata = Length@Last@finalFit; (* length of residuals vector *)
EstVar = finalFit[[16]]; (* would be = Total[resid^2]/(Ndata-Ncoeff) without weights *)
q1=True;
stderrSqContribs = FitPredictionIntervalStderrSq[finalFit,finalAnsatzRules,fit2dParts,\[Eta]0stuff,extraCoeffRules,productAnsatz,q1];
(* add up all error contributions in quadrature *)
stderrsq = Total[stderrSqContribs];
(* report back the 95% student-t quantile (applied on both sides, this gives a 90% interval), adding up also the variance contribution for PREDICTION interval (not mean confidence) *)
quant95 = Quantile[StudentTDistribution[Ndata-NfinalCoeffs],0.95];
Return[quant95 * Sqrt[ EstVar + stderrsq ]];
]
FitPredictionInterval[finalFit_,finalAnsatzRules_,fit2dParts_,\[Eta]0stuff_,extraCoeffRules_,productAnsatz_,etain_,chi1in_,chi2in_]:=Module[{fiterrFunc,fiterrFuncq1,fiterrs,i,eta,chi1,chi2},
(* so we can deal with both scalars and lists *)
If[Length[etain]==0,
eta = {etain};
chi1 = {chi1in};
chi2 = {chi2in};
,
eta = etain;
chi1 = chi1in;
chi2 = chi2in;
];
(* evaluate the general error estimate *)
fiterrFunc = FitPredictionIntervalFunction[finalFit,finalAnsatzRules,fit2dParts,\[Eta]0stuff,extraCoeffRules,productAnsatz];
(* avoid indeterminates at q=1, eta=0.25 *)
If[MemberQ[eta,0.25],
fiterrFuncq1 = FitPredictionIntervalFunctionq1[finalFit,finalAnsatzRules,fit2dParts,\[Eta]0stuff,extraCoeffRules,productAnsatz];
];
fiterrs = Table[fiterrFunc,{i,1,Length@eta}];
For[i=1,i<=Length@eta,i++,
If[eta[[i]]==0.25,
fiterrs[[i]]=fiterrFuncq1;
];
];
(* switch spin parametrization and insert user values *)
fiterrs = fiterrs/.{S->sTot3[\[Eta],\[Chi]1,\[Chi]2]};
Return[Table[fiterrs[[i]]/.{\[Eta]->eta[[i]],\[Chi]1->chi1[[i]],\[Chi]2->chi2[[i]]},{i,1,Length@eta}]];
]
TeXExportCoeffTable[coeffTable_,filename_]:=Module[{texTable},
texTable = StringReplace[ToString[TeXForm[NumberForm[coeffTable,3]]],{"$\\grave{ }$*${}^{\\wedge}$"->"}\\cdot10^{"}];
Export[filename, texTable, "Text"];
]
TeXExportCovarMatrixTable[fitCovar_,fitCoeffRules_,filename_]:=Module[{covarMatrix,paramNames,fullMatrix,texTable},
covarMatrix = fitCovar;
paramNames = fitCoeffRules[[All,1]];
fullMatrix = Join[List/@Join[{""},paramNames],Join[{paramNames},covarMatrix,1],2];
texTable = StringReplace[ToString[TeXForm[NumberForm[fullMatrix,3]]],{"$\\grave{ }$*${}^{\\wedge}$"->"}\\cdot10^{"}];
Export[filename,texTable,"Text"];
(* then run on the resulting file: sed -i 's/{/[/g;s/}/],/g;s/*^/e/g' PeakLuminosityUIBCovMatrix_S5.txt *)
]
SubscriptRules[rules_]:=Module[{},
Return[Table[rules[[i,1]]->Subscript[ToExpression[StringTake[ToString@rules[[i,1]],1]],ToExpression[StringTake[ToString@rules[[i,1]],2;;]]],{i,1,Length@rules}]];
]
TeXFormatAnsatz[ansatz_,formattingRules_,coeffRules_]:=Module[{subscripts,formattedAnsatz,TeXansatz},
subscripts = SubscriptRules[coeffRules];
formattedAnsatz = ansatz/.formattingRules/.subscripts;
(*
ansatzReals=Cases[AtomsList[ansatzRaw],y_Real];
ansatzRounding=Table[ansatzReals[[i]]\[Rule]NumberForm[ansatzReals[[i]],{2,2}],{i,1,Length@ansatzReals}];
formattedAnsatz = formattedAnsatz ./ .ansatzRounding;
*)
TeXansatz = TeXForm[formattedAnsatz];
TeXansatz = StringReplace[ToString@TeXansatz," _"->"_"];TeXansatz = StringReplace[TeXansatz,"{}"->""];
TeXansatz = StringReplace[TeXansatz,"\\overset{\\land }"->"\\widehat"];
TeXansatz = StringReplace[TeXansatz,"\\hat"->"\\widehat"];
Return[TeXansatz];
]
TeXExportAnsatz[ansatz_,formattingRules_,coeffRules_,filename_]:=Module[{TeXansatz},
TeXansatz = TeXFormatAnsatz[ansatz,formattingRules,coeffRules];
Export[filename, TeXansatz, "Text"];
]
GetRawTwoDAnsatz[finalfit_,fit2dParts_,productAnsatz_,constrained_]:=Module[{etaAnsatz,etaCoeffRules,etaInverseRules,SAnsatz,SCoeffRules,SInverseRules,ansatzRaw},
etaAnsatz = fit2dParts[[2]];
etaCoeffRules = fit2dParts[[3]];
etaInverseRules = Get1dInverseRules[etaAnsatz,etaCoeffRules];
SAnsatz = fit2dParts[[5]];
SCoeffRules = fit2dParts[[6]];
SInverseRules = Get1dInverseRulesS[SAnsatz,SCoeffRules,productAnsatz];
ansatzRaw = (fit2dParts[[1]]/.etaInverseRules/.SInverseRules);
If[constrained,
ansatzRaw = ansatzRaw /. fit2dParts[[8]] /. fit2dParts[[9]] /. {4.->4, 16.->16, 64.->64, 256.->256} /. {-4.->-4, -16.->-16, -64.->-64, -256.->-256};
];
Return[ansatzRaw];
]
GetAllRules[finalfit_,fit2dParts_,twodrules_]:=Module[{allRules},
allRules=Join[fit2dParts[[3]],fit2dParts[[6]],fit2dParts[[8]],fit2dParts[[9]],finalfit[[2]]]; (* etaCoeffRules, SCoeffRules, 2dconstraintsNumerator, 2dconstraintsDenominator, finalfitCoeffRules *)
If[ListQ@twodrules&&(Length@twodrules>0),
allRules=Join[allRules,twodrules];
];
Return[allRules];
]
TeXExportTwoDAnsatz[finalfit_,fit2dParts_,twodrules_,formattingRules_,productAnsatz_,constrained_,filename_]:=Module[{ansatzRaw,allRules,zeroRules},
ansatzRaw = GetRawTwoDAnsatz[finalfit,fit2dParts,productAnsatz,constrained];
allRules = GetAllRules[finalfit,fit2dParts,twodrules];
If[Not@productAnsatz,
zeroRules=Table[ToExpression["f"<>ToString@i<>"0"]->0,{i,0,Exponent[Numerator[fit2dParts[[1]]],S]}];
allRules=Join[allRules,zeroRules];
];
TeXExportAnsatz[ansatzRaw,formattingRules,allRules,filename];
]
TeXExportChiDiffTerms[chidiffAnsaetze_,formattingRules_,coeffRules_,filename_]:=Module[{fullTeX,numChiDiffTerms,thisTerm},
fullTeX = "";
numChiDiffTerms = Length@chidiffAnsaetze;
If[numChiDiffTerms>=1,
thisTerm = TeXFormatAnsatz[chidiffAnsaetze[[1]],formattingRules,coeffRules];
fullTeX = "A_1(\\eta)&="<>thisTerm;
];
If[numChiDiffTerms>=2,
thisTerm = TeXFormatAnsatz[chidiffAnsaetze[[2]],formattingRules,coeffRules];
fullTeX = fullTeX<>" \\\\\nA_2(\\eta)&="<>thisTerm;
];
If[numChiDiffTerms>=3,
thisTerm = TeXFormatAnsatz[chidiffAnsaetze[[3]],formattingRules,coeffRules];
fullTeX = fullTeX<>" \\\\\nA_3(\\eta)&="<>thisTerm;];
Export[filename, fullTeX, "Text"];
]
TeXExportFinalAnsatz[finalfit_,fit2dParts_,twodrules_,chidiffAnsaetze_,formattingRules_,productAnsatz_,twodConstrained_,filename_]:=Module[{ansatzRaw,numChiDiffTerms,allRules},
ansatzRaw = GetRawTwoDAnsatz[finalfit,fit2dParts,productAnsatz,twodConstrained];
numChiDiffTerms = Length@chidiffAnsaetze;
If[numChiDiffTerms>=1,
ansatzRaw = ansatzRaw + chidiffAnsaetze[[1]] (\[Chi]1-\[Chi]2);
];
If[numChiDiffTerms>=2,
ansatzRaw = ansatzRaw + chidiffAnsaetze[[2]] (\[Chi]1-\[Chi]2)^2;
];
If[numChiDiffTerms>=3,
ansatzRaw = ansatzRaw + chidiffAnsaetze[[3]] S (\[Chi]1-\[Chi]2);
];
allRules = GetAllRules[finalfit,fit2dParts,twodrules];
TeXExportAnsatz[ansatzRaw,formattingRules,allRules,filename];
]
TeXExportCoeffTable[coeffRules_,covar_,filename_]:=Module[{coeffNames,coeffErrs,headerLine,coeffTable,texTable},
coeffNames = coeffRules[[All,1]]/.SubscriptRules[coeffRules];
coeffErrs = Sqrt[Diagonal[covar]];
headerLine = {"estimate","std.err.","rel.err.[%]"};
(*
coeffTable = TableForm[Table[{coeffRules[[i,2]],coeffErrs[[i]],Abs[100*coeffErrs[[i]]/coeffRules[[i,2]]]},{i,Length@coeffRules}],TableHeadings\[Rule]{coeffNames,headerLine}];
texTable = StringReplace[ToString[TeXForm[NumberForm[coeffTable,3]]],{"$\\grave{ }$*${}^{\\wedge}$"\[Rule]"}\\cdot10^{"}];
Export[filename,texTable, "Text"];
*)
(*
TeXExportTabularTable[Map[NumberForm[#,3]&,Table[{coeffRules[[i,2]],coeffErrs[[i]],Abs[100*coeffErrs[[i]]/coeffRules[[i,2]]]},{i,Length@coeffRules}],{2}], filename, coeffNames, headerLine]
*)
TeXExportTabularTable[Table[{NumberForm[coeffRules[[i,2]],3],NumberForm[coeffErrs[[i]],3],NumberForm[Abs[100*coeffErrs[[i]]/coeffRules[[i,2]]],{3,1}]},{i,Length@coeffRules}],
filename, coeffNames, headerLine, 0]
]
TeXExportTabularTable[datatable_,filename_,rowLabels_,colHeadings_,padZeroes_]:=Module[{nRows,tableFormatted,texTable,tableAtoms,maxDecDigits},
nRows = Length@datatable;
tableFormatted = TableForm[datatable, TableHeadings->{rowLabels,colHeadings}];
texTable = ToString[TeXForm[tableFormatted]];
texTable = StringReplace[texTable,{"array"->"tabular"}];
texTable = StringReplace[texTable,{"&"->"$&$"}]; (* close and open math mode at each field separator *)
texTable = StringReplace[texTable,{"$&$"->"&$"},1]; (* remove the extra open$ at first separator *)
texTable = StringReplace[texTable,{"$\$&$$"->"\&"},1];(* remove the extra elements when exporting & *)
texTable = StringReplace[texTable,{"\\\\"->"$\\\\"}]; (* close math mode at each line break *)
texTable = StringReplace[texTable,{"\\\\"->"\\\\$"},nRows]; (* reopen math mode after each line break, but only for numrows (not at last one) *)
texTable = StringReplace[texTable,{"\\\\"->"\\\\\\hline"},1]; (* add a horizontal line after the first (headers) line *)
texTable = StringReplace[texTable,{"{tabular}{c"->"{tabular}{l"}]; (* make row-headers column left-centered *)
texTable = StringReplace[texTable,{"cc}"->"cc}\\hline\\hline"}]; (* add initial double line *)
texTable = StringReplace[texTable,{"\\end{tabular}"->"\\hline\\hline\\end{tabular}"}]; (* final double lines at end of table *)
If[padZeroes==1, (* pad trailing zeroes *)
texTable = StringReplace[texTable,{".\\times"->".0\\times"}];
,
If[padZeroes==2,
texTable = StringReplace[texTable,{".\\times"->".00\\times"}];
texTable = StringReplace[texTable,Table["."<>ToString@i<>"\\times"->"."<>ToString@i<>"0\\times",{i,0,9}]];
,
If[padZeroes==3,texTable = StringReplace[texTable,{".\\times"->".000\\times"}];
texTable = StringReplace[texTable,Table["."<>ToString@i<>"\\times"->"."<>ToString@i<>"00\\times",{i,0,9}]];
texTable = StringReplace[texTable,Flatten[Table[Table["."<>ToString@i<>ToString@j<>"\\times"->"."<>ToString@i<>ToString@j<>"0\\times",{i,0,9}],{j,0,9}]]];
];
];
];
Export[filename, texTable, "Text"];
]
PyExportFinalFit[fit_,extraFormattingRules_,filename_]:=Module[{fitCform,pyRules1,pyRules2,pyRules3,fitPyform},
fitCform = ToString@CForm[fit];
(* formatting hacks: greek letters, powers, Shat (FIXME: not yet variable wrt to effective-spin parameter choice!) *)
pyRules1 = {"\[Eta]"->"eta","\[Chi]"->"chi","Sqrt(2)"->"sqrt2","Sqrt(3)"->"sqrt3"};
pyRules2 = Join[{"(chi1 - chi2)"->"chidiff","Power(chi1 - chi2,2)"->"chidiff2"},
Table["Power(eta,"<>ToString@i<>")"->"eta"<>ToString@i,{i,2,9}],
Table["Power(S,"<>ToString@i<>")"->"S"<>ToString@i,{i,2,9}],
{"Power(1 - 4*eta,0.5)"->"sqrt1m4eta"}];
pyRules3 = {"S"->"Shat","1 "->"1. ","0 "->"0. "};
fitPyform = StringReplace[StringReplace[StringReplace[StringReplace[fitCform,pyRules1],pyRules2],pyRules3],extraFormattingRules];
fitPyform = fitPyform<>"\n";
Export[filename, fitPyform, "Text"];
]
PyExportFinalAnsatz[finalfit_,fit2dParts_,chidiffAnsaetze_,keepfi0_,extraFormattingRules_,filename_]:=Module[{ansatzRaw,numChiDiffTerms,zeroRules,ansatzCform,
pyRules1,pyRules2,pyRules3,ansatzPyform,sqrti,sqrtRules,fracs,fracStrings,fracReals,fracRules,twodConstrained},
twodConstrained = True;
ansatzRaw = GetRawTwoDAnsatz[finalfit,fit2dParts,productAnsatz,twodConstrained];
numChiDiffTerms = Length@chidiffAnsaetze;
If[numChiDiffTerms>=1,
ansatzRaw = ansatzRaw + chidiffAnsaetze[[1]] (\[Chi]1-\[Chi]2);
];
If[numChiDiffTerms>=2,
ansatzRaw = ansatzRaw + chidiffAnsaetze[[2]] (\[Chi]1-\[Chi]2)^2;
];
If[numChiDiffTerms>=3,
ansatzRaw = ansatzRaw + chidiffAnsaetze[[3]] S (\[Chi]1-\[Chi]2);
];
If[Not@keepfi0,
zeroRules = Table[ToExpression["f"<>ToString@i<>"0"]->0,{i,0,Exponent[Numerator[fit2dParts[[1]]],S]}];
ansatzRaw = ansatzRaw/.zeroRules;
];
(* hack to avoid rounding of square-roots *)
sqrti = {2,3,5,6,7};
sqrtRules = Table[Sqrt[sqrti[[i]]]->ToExpression["sqrt"<>ToString@sqrti[[i]]<>""],{i,Length@sqrti}];
(* hack to avoid rounding of fractions, first part: replace by named string before ToString@CForm[SetPrecision[...]] *)
fracs = {1/3,2/3,4/3,5/3};
fracs = Join[fracs,-fracs];
fracStrings = Table[ToExpression["frac"<>StringReplace[ToString@Numerator@fracs[[i]],{"-"->"m"}]<>ToString@Denominator@fracs[[i]]],{i,Length@fracs}];
fracRules = Table[fracs[[i]]->fracStrings[[i]],{i,Length@fracs}];
ansatzCform = ToString@CForm[SetPrecision[ansatzRaw/.sqrtRules/.fracRules,3]];
(* hack to avoid rounding of fractions, second part: replace named string back to real number fractions *)
fracReals = Table[ToString@SetPrecision[Numerator@fracs[[i]],2]<>"/"<>ToString@SetPrecision[Denominator@fracs[[i]],2],{i,Length@fracs}];
fracRules = Table[ToString@fracStrings[[i]]->fracReals[[i]],{i,Length@fracs}];
ansatzCform = StringReplace[ansatzCform,fracRules];
(* more formatting hacks: greek letters, powers, Shat (FIXME: not yet variable wrt to effective-spin parameter choice!) *)
pyRules1 = {"\[Eta]"->"eta","\[Chi]"->"chi"};pyRules2 = Join[{"(chi1 - 1.*chi2)"->"chidiff","Power(chi1 - 1.*chi2,2)"->"chidiff2"},
Table["Power(eta,"<>ToString@i<>")"->"eta"<>ToString@i,{i,2,9}],
Table["Power(S,"<>ToString@i<>")"->"S"<>ToString@i,{i,2,9}],
{"Power(1. - 4.*eta,0.5)"->"sqrt1m4eta"},{"Sqrt(1. - 4.*eta)"->"sqrt1m4eta"}];
pyRules3 = {"S"->"Shat"};
ansatzPyform = StringReplace[StringReplace[StringReplace[StringReplace[ansatzCform,pyRules1],pyRules2],pyRules3],extraFormattingRules];
ansatzPyform = ansatzPyform<>"\n";
Export[filename, ansatzPyform, "Text"];
]
PyExportFinalFitCoeffs[finalfit_,fit2dParts_,all2dconstraints_,filename_]:=Module[{pyCoeffsList,pyCoeffsString},
pyCoeffsList = Join[Table[ToString@fit2dParts[[3,i,1]]<>" = "<>ToString@CForm@fit2dParts[[3,i,2]],{i,Length@fit2dParts[[3]]}],
Table[ToString@fit2dParts[[6,i,1]]<>" = "<>ToString@CForm@fit2dParts[[6,i,2]],{i,Length@fit2dParts[[6]]}],
Table[ToString@all2dconstraints[[i,1]]<>" = "<>ToString@CForm[all2dconstraints[[i,2]]/.{1->1.,0->0.}],{i,Length@all2dconstraints}],
(*{ToString@\[Eta]0derconstv3[[1,1]]<>" = "<>ToString@CForm[\[Eta]0derconstv3[[1,2]]]}, *)
Table[ToString@finalfit[[2,i,1]]<>" = "<>ToString@CForm@finalfit[[2,i,2]],{i,Length@finalfit[[2]]}]
];
pyCoeffsString = "";
For[i=1, i <= Length@pyCoeffsList, i++,
pyCoeffsString = pyCoeffsString<>" "<>pyCoeffsList[[i]]<>"\n"];
Export[filename, pyCoeffsString, "Text"];
]
SupplExportAllFitCoeffs[finalfit_,fit2dParts_,all2dconstraints_,eta0covar_,filename_]:=Module[{coeffsList,coeffsString,stdErrsEta,stdErrsS,stdErrs2D,stdErrsFinal},
stdErrsEta = Sqrt[Diagonal[fit2dParts[[4]]]];
stdErrsS = Sqrt[Diagonal[fit2dParts[[7]]]];
stdErrs2D = Join[Sqrt[Diagonal[eta0covar]],ConstantArray[0.0,Length[all2dconstraints]-Length[eta0covar]]];
stdErrsFinal = Sqrt[Diagonal[finalfit[[14]]]];
coeffsList = Join[Table[ToString@fit2dParts[[3,i,1]]<>" "<>ToString@CForm@fit2dParts[[3,i,2]]<>" "<>ToString@CForm@stdErrsEta[[i]],{i,Length@fit2dParts[[3]]}],
Table[ToString@fit2dParts[[6,i,1]]<>" "<>ToString@CForm@fit2dParts[[6,i,2]]<>" "<>ToString@CForm@stdErrsS[[i]],{i,Length@fit2dParts[[6]]}],
Table[ToString@all2dconstraints[[i,1]]<>" "<>ToString@CForm[all2dconstraints[[i,2]]/.{1->1.,0->0.}]<>" "<>ToString@CForm@stdErrs2D[[i]],{i,Length@all2dconstraints}],
(*{ToString@\[Eta]0derconstv3[[1,1]]<>" "<>ToString@CForm[\[Eta]0derconstv3[[1,2]]]}, *)
Table[ToString@finalfit[[2,i,1]]<>" "<>ToString@CForm@finalfit[[2,i,2]]<>" "<>ToString@CForm@stdErrsFinal[[i]],{i,Length@finalfit[[2]]}]
];
coeffsString = "# coeff estimate stderr\n";
For[i=1, i <= Length@coeffsList, i++,
coeffsString = coeffsString<>coeffsList[[i]]<>"\n"];
coeffsString = StringReplace[coeffsString,{". "->".0 ",".\n"->".0\n"}];
Export[filename, coeffsString, "Text"];
]
SupplExportCovarMatrix[covar_,coeffRules_,fileName_]:=Module[{coeffNames,fullMatrix,covarString},
coeffNames = coeffRules[[All,1]];
covarString = "# coeffs:";
For[i=1, i <= Length@coeffRules, i++,
covarString = covarString<>" "<>ToString@coeffRules[[i,1]];
];
covarString = covarString<>"\n";
For[i=1, i <= Length@covar, i++,
For[j=1, j <= Length@covar, j++,
covarString = covarString<>" "<>ToString@NumberForm[covar[[i,j]], {16, 16}, ExponentFunction -> (If[-10 < # < 10, Null, #] &)];
];
covarString = covarString<>"\n"; ];
Export[fileName, covarString, "Text"];
]
Options[AIC]={"Weights"->False};
AIC[data_,fit_,fitvars_,OptionsPattern[]]:=Module[{bracketedvars,coeff,datapnts,err,n,res,ress,weigths},
n=Length@data;
coeff=Length@fitvars;
weigths=OptionValue["Weights"];
If[weigths,err=data[[All,-1]];datapnts=data[[All,-2]];,err=1;datapnts=data[[All,-1]];];
res=Table[data[[i,-2]]-(fit/.Table[fitvars[[j]]->data[[i,j]],{j,Length@fitvars}]),{i,Length@data}];
ress=Total[err (res)^2];
n + n Log[2\[Pi]]+n Log[ress/n]+2(coeff+1)
]
Options[AICc]=Options[AIC];
AICc[data_,fit_,fitvars_,OptionsPattern[]]:=Module[{coeff,n,res,newfit,bracketedvars,ress,weigths},
weigths=OptionValue["Weights"];
n=Length@data;
coeff=Length@fitvars;
AIC[data,fit,fitvars,"Weights"->weigths]+(2*coeff(coeff+1) )/(n - coeff -1)
]
BIC[data_,fit_,fitvars_,coeff_]:=Module[{n,res,newfit,bracketedvars,ress},
n=Length@data;
bracketedvars=StringReplace[StringReplace[ToString@fitvars,"{"->"["],"}"->"]"];
newfit=ToExpression[ToString@fit<>bracketedvars];
res=Table[data[[i,-2]]-(newfit/.Table[fitvars[[j]]->data[[i,j]],{j,Length@fitvars}]),{i,Length@data}];
ress=Total[data[[All,-1]](res)^2];
n + n Log[2\[Pi]]+n Log[ress/n]+Log[n](coeff+1)
]
Residuals[data_,fit_,vars_,OptionsPattern[{"Verbose"->False,"Relative"->False}]]:=Module[{datvals,relative,fitvals,res,verbose},
verbose=OptionValue["Verbose"];
relative=OptionValue["Relative"];
If[verbose,Print["variables -> ",vars]];
datvals=TakeColumn[data,-1];
fitvals=fit/.Table[(vars[[j]]->data[[All,j]]),{j,Length@vars}];
If[relative,res=1-fitvals/datvals,res=datvals-fitvals]
]
KullbagLeiblerDiv[p1_?ListQ,p2_?ListQ]:=Module[{fun},
fun=1.Function[{p,q},Limit[p*Log[(p+\[Epsilon])/(q+\[Epsilon])],\[Epsilon]->0]][p2,p1];
Total[fun]
]
KL[data1_,data2_,lim_]:=Module[{n,res},
n=Length@data1;
NIntegrate[PDF[data1,x] Log[PDF[data1,x]/PDF[data2,x]], {x,-lim,lim}]
]
JensenShanonDiv[p1_?ListQ,p2_?ListQ]:=Module[{m,fun},
m=0.5 (p1+p2);
fun=Function[{p,q,m},Limit[p*Log[(p+\[Epsilon])/(m+\[Epsilon])]+ q*Log[(q+ \[Epsilon])/(m+\[Epsilon])],\[Epsilon]->0]][p1,p2,m];
Total[fun]
]
JS[data1_,data2_,lim_]:=Module[{n,res},
n=Length@data1;
1/2 NIntegrate[PDF[data1,x] Log[PDF[data1,x]/(PDF[data2,x]+ PDF[data1,x])], {x,-lim,lim}] + 1/2 NIntegrate[PDF[data2,x] Log[PDF[data2,x]/(PDF[data2,x]+PDF[data1,x])] , {x,-lim,lim}]
]
ComputeEdges[pts_]:=Module[{ptsx,auxvar,auxvar2,nears,i},
ptsx=SortBy[pts,First];
auxvar={};
i=1;
AppendTo[auxvar,{ptsx[[i]]}];
While[i<= Length@ptsx-1,
If[ptsx[[i+1,1]]==ptsx[[i,1]],i=i+1,AppendTo[auxvar,{ptsx[[i]]}];i=i+1]
];
AppendTo[auxvar,{ptsx[[i]]}];
auxvar=Flatten[auxvar,1];
i=Length@ptsx-1;
While[i> 1,
If[ptsx[[i+1,1]]==ptsx[[i,1]],i=i-1,AppendTo[auxvar,ptsx[[i+1]]];i=i-1]
];
AppendTo[auxvar,ptsx[[i]]];
Do[auxvar=AppendTo[auxvar,0.5auxvar[[1]]+0.5auxvar[[-1]]],{i,3}];
auxvar
]
CredibleRegion[data_,level_]:=Module[{datasrt,prob,cumprob,pbound},
(* Last column must be the PDF *)
datasrt=SortBy[data,Last];
prob=TakeColumn[datasrt,-1];
cumprob=Accumulate[prob]/Total[prob];
pbound=Quiet@Position[cumprob,_?(#>= (1-level) cumprob[[-1]]& ),1][[1,1]];
ComputeEdges[datasrt[[pbound-1;;-1]]][[All,1;;-2]]
]
CredibleInterval[data_,level_]:=Module[{datasrt,prob,cumprob,pbound},
(* Last column must be the PDF *)
datasrt=SortBy[data,Last];
prob=TakeColumn[datasrt,-1];
cumprob=Accumulate[prob]/Total[prob];
pbound=Quiet@Position[cumprob,_?(#>= (1-level) cumprob[[-1]]& ),1][[1,1]];
datasrt[[pbound;;-1]]
]
Generate1DPolynomialAnsatz[CoefficientPrefixString_?StringQ,variable_,MinOrder_?IntegerQ,MaxOrder_?IntegerQ]:=Module[{ansatz,i,j},
ansatz = Total/@Table[ToExpression[CoefficientPrefixString<>ToString@i] * variable^i,{j,MinOrder,MaxOrder},{i,0,j}];
Transpose[{ansatz,Table[{variable},{i,Length@ansatz}]}]
]
CleanAnsatzParams[paramsGuess_, vars_]:=Module[{pos,params,varsStr,paramsStr,tmp},
params = DeleteDuplicates@Flatten@paramsGuess;
tmp=Select[params, Not@NumberQ@#& ];
varsStr = ToString/@vars;
paramsStr= ToString/@tmp;
(* Clean the variables from the parameters and identify variables *)
pos = Flatten[Position[paramsStr,#]&/@varsStr,1];
params=Delete[tmp,pos];
params
];
(* ::Code::Initialization:: *)
Options[DataFitFunction]={
"FitAll" -> True,
"AxesTag" -> "Amplitude",
"Verbose" -> 1,
"StatisticalTest" -> "AIC",
"Sorted" -> True,
"Weights" -> {},
"PlotRange" -> All,
"ToolTipTags" -> "",
"Domain" -> {0,1},
"GetIntervals" -> False
};
DataFitFunction[dataRAW_?ListQ,ansatzList_?ListQ,OptionsPattern[]]:=Module[{ansatz,ansatzaparams,ansatzvars,cleartrue,fit,stats,verbose,
ansatzparamsStr,ansatzvarsStr,outparams,fittab,plot1,plot2,datasorted,plotdomain,myvar,plotdomain1,axestag,aic,aicc,bic,rsquared,rmse,
statisticaltest,stattoplot,sortfield,residuals,minexp,maxexp,minexp2,maxexp2,sorted,paramerrors,residualsplot,weights,plotrange,print,
ansatzaparamsGuess,paramTStat,overview,data,dimWeights,tooltiptags,inner,confidenceinttable,userPlotDomain,vcov,
estvar,confbands,predbands,getintervals,loglikelihood},
verbose = OptionValue["Verbose"];
axestag = OptionValue["AxesTag"];
statisticaltest = OptionValue["StatisticalTest"];
sorted = OptionValue["Sorted"];
weights = OptionValue["Weights"];
plotrange = OptionValue["PlotRange"];
tooltiptags = OptionValue["ToolTipTags"];
userPlotDomain = OptionValue["Domain"];
getintervals = OptionValue["GetIntervals"];
If[TrueQ@weights,
dimWeights = Last@Dimensions@dataRAW;
weights = TakeColumn[dataRAW,dimWeights];
data = TakeColumn[dataRAW, Range[dimWeights-1]],
data = dataRAW[[All,1;;Length@dataRAW[[1]]]];
If[verbose,Print["Raw data with no weights on the lists!"]];
];
If[verbose>=2,
print[x___]:=Print[x],
print[x___]:={}
];
ansatz = Chop/@ (ansatzList[[All,1]]);
ansatzvars = ansatzList[[All,2]];
ansatzaparamsGuess = N/@AtomsList/@ansatz;
ansatzaparams = Table[CleanAnsatzParams[ansatzaparamsGuess[[i]], ansatzvars[[i]]],{i,Length@ansatzvars}];
(*print["Ans\[ADoubleDot]tze -> ", ansatz//TableForm];
print["#\[NonBreakingSpace]data points -> ", Length@data];
print["#\[NonBreakingSpace]data points -> ", Length@Union[data, SameTest->(Abs[#1[[1]]-#2[[1]]] < 0.01 &)], " identifying x-values deviating < 0.01"];
*)
(*Print["Weights: ", If[Length@weights>0,weights,Automatic]];*)
fittab=Table[
If[Length@ansatzaparams[[i]]>=1,
fit = NonlinearModelFit[data,ansatz[[i]],ansatzaparams[[i]],ansatzvars[[i]],MaxIterations->1000,Weights->If[Length@weights>0,weights,Automatic]];
outparams = fit["BestFitParameters"];
stats = fit["ParameterTable"];
paramerrors = fit["ParameterErrors"];
paramTStat = fit["ParameterTStatistics"];
residuals = fit["FitResiduals"];
aic = fit["AIC"];
aicc = aic + (2*Length@outparams(Length@outparams+1) )/(Length@data - Length@outparams -1);
bic = fit["BIC"];
rsquared = fit["RSquared"];
rmse = Sqrt@Mean[residuals^2];
confidenceinttable = fit["ParameterConfidenceIntervalTable"];
vcov = fit["CovarianceMatrix"];
estvar = fit["EstimatedVariance"];
If[getintervals,
confbands = fit["MeanPredictionBands",ConfidenceLevel->0.9];
predbands = fit["SinglePredictionBands",ConfidenceLevel->0.9];
,
confbands = {};
predbands = {};
];
(* 1 2 3 4 5 6 7 8 9 *)
{Normal@fit, outparams, stats, aic, aicc, bic, rsquared, ScientificForm@rmse, paramerrors,
(* 10 11 12 13 14 15 *)
Sort[CombineColumns[ansatzaparams[[i]],paramTStat],Abs@#1[[2]]<Abs@#2[[2]]&], Length@ansatzaparams[[i]], rmse, confidenceinttable, vcov, ansatz[[i]],
(* 16 17 18 19 *)
estvar, confbands, predbands, residuals },
{ansatz[[i]],"dummy","dummy","dummy","dummy","dummy","dummy","dummy","dummy","dummy"}]
,{i,Length@ansatzList}];
Which[statisticaltest=="AIC",
stattoplot=aic;
sortfield=4;
,
statisticaltest=="AICc",
stattoplot=aicc;
sortfield=5;
,
statisticaltest=="BIC",
stattoplot=bic;
sortfield=6;
,
statisticaltest=="RSquared",
stattoplot=rsquared;
sortfield=7;
,
statisticaltest=="RMSE",
stattoplot=rmse;
sortfield=12;
];
If[sorted, fittab = SortBy[fittab,#[[sortfield]]& ]];
overview = TakeColumn[fittab,{1,11,7,8,5,6,3}]; (* fit eqn, numParams, Rsquared, RSME, AICc, BIC, coefficient table *)
residuals = fittab[[All,19]];
(*
resVSaic = TakeColumn[fittab, {4,7}];
resVSaicc = TakeColumn[fittab, {5,7}];
resVSbic = TakeColumn[fittab, {6,7}];Print[ListLogPlot[{resVSaic,resVSaicc,resVSbic},PlotRange\[Rule]All,PlotStyle\[Rule]{Red,Blue,Green},PlotLegends\[Rule] {"AIC","AICc","BIC"},
Frame\[Rule] True,FrameLabel\[Rule] {"*IC*","RMSE"}]];
*)
If[verbose>=1,
Print[overview//TableForm];
];
If[verbose>=2,
Which[Length@ansatzvars[[1]]==1,
(* {minexp,maxexp}={Min[Exponent[#,ansatzvars]&/@ansatz],Max[Exponent[#,ansatzvars]&/@ansatz]}; *)
{minexp,maxexp}={Min[Table[Length[fittab[[i,2]]],{i,Length@ansatzList}]],Max[Table[Length[fittab[[i,2]]],{i,Length@ansatzList}]]};
datasorted=SortBy[data,First];
plotdomain={Min[First@First@datasorted,First@userPlotDomain],Max[First@Last@datasorted,Last@userPlotDomain]};
(*plot1=Plot[Evaluate[Table[fittab[[i,1]]/.ansatzvars[[i,1]]\[Rule]x,{i,Length@ansatzList}]],{x,First@plotdomain,Last@plotdomain},
Epilog\[Rule]{Point[Table[Tooltip[#,If[ListQ@tooltiptags,tooltiptags[[i]],ToString[data[[i]]]]],{i,Length@data}]&/@data]},
PlotRange\[Rule]plotrange,Frame\[Rule]True,PlotStyle\[Rule]Red,FrameLabel\[Rule]{Style[ToString@ansatzvars[[1,1]],14],Style["f("<>ToString@ansatzvars[[1,1]]<>")",14]}];*)
plot1=Plot[Evaluate[Table[fittab[[i,1]]/.ansatzvars[[i,1]]->x,{i,Length@ansatzList}]],{x,First@plotdomain,Last@plotdomain},PlotLegends->Evaluate[Table[fittab[[i,1]],{i,Length@ansatzList}]]];
plot2=ListPlot[Table[Tooltip[data[[i]],If[ListQ@tooltiptags,tooltiptags[[i]],ToString[data[[i]]]]],{i,Length@data}], PlotStyle->Red];
stattoplot=Transpose[{Table[Length[fittab[[i,2]]],{i,Length@ansatzList}],TakeColumn[fittab,sortfield]}];
inner=(Inner[List,data[[All,1]],#,List]&/@residuals);
Print@{Show[plot1, plot2, ImageSize->450, Frame->True, FrameLabel->{Style[ToString@ansatzvars[[1,1]],14],Style["f("<>ToString@ansatzvars[[1,1]]<>")",14]}],
ListPlot[Partition[stattoplot,1],Frame->True,FrameLabel->{Style["\!\(\*SubscriptBox[\(N\), \(coeff\)]\)",14],Style[statisticaltest,14]},ImageSize->450,PlotRange->{{minexp-1,maxexp+1},All},
PlotLegends->Table[fittab[[i,1]],{i,Length@ansatzList}]],
ListPlot[Table[Tooltip[#[[i]],If[ListQ@tooltiptags,tooltiptags[[i]],ToString[data[[i]]]]],{i,Length@data}]&/@inner,PlotRange->{plotdomain,All},
ImageSize->450,Frame->True,FrameLabel->{Style[ToString@ansatzvars[[1,1]],14],Style["data-fit",14]}]
};
,
Length@ansatzvars[[1]]==2,
(*
{minexp,maxexp} ={Min[Exponent[#,ansatzvars[[1,1]]]&/@ansatz],Max[Exponent[#,ansatzvars[[1,1]]]&/@ansatz]};
{minexp2,maxexp2}={Min[Exponent[#,ansatzvars[[1,2]]]&/@ansatz],Max[Exponent[#,ansatzvars[[1,2]]]&/@ansatz]};
*)
plot2=ListPointPlot3D[data,PlotRange->plotrange,PlotStyle->PointSize[0.02],PlotLegends->{"data"}];
plot1=Plot3D[Evaluate[Table[fittab[[i,1]]/.ansatzvars[[i,1]]->x/.ansatzvars[[1,2]]->y,{i,Length@ansatzList}]],
{x,First@First@userPlotDomain,Last@First@userPlotDomain}, {y,First@Last@userPlotDomain,Last@Last@userPlotDomain}, PlotRange->plotrange,
AxesLabel->{Style[ToString@ansatzvars[[1,1]],14],Style[ToString@ansatzvars[[1,2]],14],Style["f("<>ToString@ansatzvars[[1,1]]<>","<>ToString@ansatzvars[[1,2]]<>")",14]},
PlotLegends->(*Evaluate[Table[fittab[[i,1]],{i,Length@ansatzList}]]*)Table["Ansatz "<>ToString[i],{i,Length@ansatzList}]];
stattoplot=Table[{Exponent[fittab[[i,1]], ansatzvars[[i,1]]],Exponent[fittab[[i,1]], ansatzvars[[i,2]]],fittab[[i,4]]},{i,Length@ansatzList}];
residualsplot=Table[Transpose[{data[[All,1]],data[[All,2]],residuals[[i]]}],{i,Length@ansatzList}];
Print@{Show[plot1,plot2,ImageSize->450],
(* ListPointPlot3D[stattoplot,PlotStyle\[Rule]Red,AxesLabel\[Rule]{"order "<>ToString@ansatzvars[[1,1]],"order "<>ToString@ansatzvars[[1,2]],statisticaltest},ImageSize\[Rule]450,PlotRange\[Rule]{{minexp-0.2,maxexp+0.2},{minexp2-0.2,maxexp2+0.2},Automatic}],*)
(*myListPlot3D[residualsplot,AxesLabel\[Rule]{Style[ToString@ansatzvars[[1,1]],14],Style[ToString@ansatzvars[[1,2]],14],""},PlotLabel->Style["Residuals",14],ImageSize\[Rule]450],*)
myListPlot3D[residualsplot,AxesLabel->{Style[ToString@ansatzvars[[1,1]],14],Style[ToString@ansatzvars[[1,2]],14],Style["data-fit",14]},ImageSize->450,PlotRange->plotrange]};
Print[];
,True,
Print["More than 2 variables can not be plotted :) "];
];
];
Return[fittab];
]
Options[DataFitFunctionAll]={
"NSFitindex" -> 1, (* which fit to pick *)
"q1Fitindex" -> 1,
"FitCase" -> 1, (* 1: 1D fits, 2: 1D+2D fits, return 2D only: 3: 1D+2D fits, return both; 4: 1D+2D+3D fits, return only 2D+3D; 5: 1D+2D+3D fits, return all, but with reduced verbosity *)
"MassRatioUnequalFits" -> {1,1.5,2.,3.,4.,8.,18.},
"NonSpinningAnsatzList" -> {{a1 \[Eta]+a2 \[Eta]^2+a3 \[Eta]^3,{\[Eta]}},{a1 \[Eta]+a2 \[Eta]^2+a3 \[Eta]^3+a4 \[Eta]^4,{\[Eta]}},{a1 \[Eta]+a2 \[Eta]^2+a3 \[Eta]^3+a4 \[Eta]^4+a5 \[Eta]^5,{\[Eta]}}},"EqualBHAnsatzList" -> {{b1 S+b2 S^2+b3 S^3,{S}},{b1 S+b2 S^2+b3 S^3+b4 S^4,{S}},{b1 S+b2 S^2+b3 S^3+b4 S^4+b5 S^5,{S}},{b1 S+b2 S^2+b3 S^3+b4 S^4+b5 S^5+b6 S^6,{S}},{b1 S+b2 S^2+b3 S^3+b4 S^4+b5 S^5+b6 S^6+b7 S^7,{S}},
{b1 S+b2 S^2+b3 S^3+b4 S^4+b5 S^5+b6 S^6+b7 S^7+b8 S^8,{S}}},
"Addq1Ansatz" -> {},
"AddNSAnsatz" -> {},
"AddGenAnsatz" -> {},
"StatisticalTest" -> "BIC",
"AnsatzTestCombinations" -> {{1},{1}},
"PlotRange" -> All,
"RationalFunctions" -> True,
"ProductAnsatz" -> False,
"GeneralAnsatzRules" -> {},
"ToolTipTags" -> "",
"ProductGeneralizeOrder" -> 2,
"ProductGeneralizeOrderDenom" -> -1,
"ProductGeneralizeMinPowerNum" -> -1,
"TwoDSpinAnsatzOrderNum" -> -1,
"TwoDSpinAnsatzOrderDenom" -> -1,
"SpinDiffWeights" -> True,
"Weights" -> True,
"AddSpinDiffLinAnsatz" -> {},
"AddSpinDiffQuadAnsatz" -> {},
"AddSpinDiffMixAnsatz" -> {},
"SpinDifferenceRules" -> {},
"\[Eta]LimitEqualSpin" -> 0,
"GetIntervals" -> False,
"SpinParameter" -> sTot3,
"SpinDiffParameter" -> \[Chi]diffstan(*(\[Chi]1-\[Chi]2)*),
"FastAnalysis" -> False,
"Addf00ToAnsatz" -> False
};
DataFitFunctionAll[data_?ListQ,OptionsPattern[]]:=Module[{datans,nsfits,nsansatz,ansatzvars,dataq1,dataq1sorted,ansatzq1,nsfitindex,q1fit,ansatz,nsfit,sys,solNS,coeffs,ansatzS,solq1,ansatzFinal,
ansatzFinal\[Delta]\[Chi],allData,fitEq,allDataEqS,dataeqs,my\[Eta],myd\[Eta],dataq1Sym,pos,massratioplots,myfit,tabfits,uneqfits,massratiounequalfits,
uneqfit2,uneqfit,uneqfit3,q1fitopt,q1fitindex,expans,expns,mydq,thisq,oneDverbosity,perqVerbosity,
ansatzFinal\[Delta]\[Chi]1,ansatzFinal\[Delta]\[Chi]2,lplot,lqplotl,lqplotq,lmplotl,lmplotm,lmqplotl,lmqplotm,lmqplotq,
addnsansatz,addq1ansatz,\[Eta]sfit,statisticaltest,nsdegree,q1degree,addgenansatz,genansatzindex,ansatzchidifflinear,ansatzchidiffquadratic,plot1,plot2,plot3,plot4,plot4a,plot4b,plot4mirror,
testansatzcomb,myfit2,q1fits,ansatzAll,fitcase,alldataeqfit,alldataeqfitv2,eqdatafit,linearpar,linearparerr,linearpar2,linearpar2err,quadpar,quadparerr,
linearfit1,quadfit,ansatz\[Eta]1,ansatz\[Eta]2,linearfit2,plotrange,uneqredunfit,uneqredunfit2,ansatzFinalS\[Delta]\[Chi],spindiffterm,
spindiffterm2,ansatz\[Eta]3,ansatz\[Eta]4,spindiffterm21,spindiffterm22,generalAnsatzRules,weights,num,denom,zeroRules,spinAnsatz,q1fitS0,fitabc,tooltiptags,dataq1tag,mysel,dataq1pos,tooltiptagsaux,
\[Eta]0Ansatz,spinpart,spindifferencerules,spindiffweights,datauneqs,dataeqs2,addspindifflinansatz,addspindiffquadansatz,addspindiffmixansatz,ansatzchidiffmix,
tooltiptagsns,tooltiptagseq,dataq1mixterm,myfitmixterm,linearpar3,linearpar3err,mixpar,mixparerr,linearfit3,mixfit,spindiffterm31,spindiffterm32,ansatzS\[Delta]\[Chi]raw,ansatzFinalS\[Delta]\[Chi]\[Chi]2,
ansatzS\[Delta]\[Chi]\[Chi]2raw,myfitmixquadterm,dataq1mixquadterm,linearpar4err,linearpar4,linearfit4,mixparerr2,mixpar2,quadparerr2,quadpar2,mixfit2,quadfit2,spindiffterm41,spindiffterm42,
spindiffterm43,unphysical,numRestrictions,denomRestrictions,extraStuff,ansatzRaw,spinAnsatzRaw,etaAnsatz,etaCoeffRules,etaCovar,SAnsatz,SCoeffRules,SCovar,getintervals,data\[Eta]S\[Chi]diff,
spinparameter,avgerrperq,avgwgthperq,extrawght,dataeqsunmass,allDataEqSv2,eqdatafitv2,dataunconstr,allDataUnconstr,spindiffparameter,
productAnsatz,productGeneralizeOrderNum,productGeneralizeOrderDenom,productGeneralizeMinPowerNum,twoDSpinAnsatzOrderNum,twoDSpinAnsatzOrderDenom,
q1fitS0orderNum,q1fitS0orderDenom,spinAnsatzToGeneralize,spinAnsatzNum,spinAnsatzDenom,q1AnsatzRaw,q1AnsatzRawS0,spinAnsatzRawNum,spinAnsatzRawDenom,spinAnsatzRawToGeneralize,
dfitspindiffterm,dfitspindiffterm21,dfitspindiffterm22,dfitspindiffterm31,dfitspindiffterm32,dfitspindiffterm41,dfitspindiffterm42,dfitspindiffterm43,fastanalysis,plot3a,dataq1a},
(* Input field "data" is expected in the form {\[Eta], \[Chi]1, \[Chi]2, value}, convention: m1 \[GreaterEqual] m2 !!!!CHECK!!!! *)
nsfitindex = OptionValue["NSFitindex"];
q1fitindex = OptionValue["q1Fitindex"];
addnsansatz = OptionValue["AddNSAnsatz"];
addq1ansatz = OptionValue["Addq1Ansatz"];
fitcase = OptionValue["FitCase"];
statisticaltest= OptionValue["StatisticalTest"];
addgenansatz = OptionValue["AddGenAnsatz"];
testansatzcomb = OptionValue["AnsatzTestCombinations"];
plotrange = OptionValue["PlotRange"];
generalAnsatzRules = OptionValue["GeneralAnsatzRules"];
tooltiptags=OptionValue["ToolTipTags"];
spindifferencerules=OptionValue["SpinDifferenceRules"];
spindiffweights=OptionValue["SpinDiffWeights"];
weights=OptionValue["Weights"];
massratiounequalfits=\[Eta]q[OptionValue["MassRatioUnequalFits"]];
addf00toansatz=OptionValue["Addf00ToAnsatz"];
addspindifflinansatz=OptionValue["AddSpinDiffLinAnsatz"];
addspindiffquadansatz=OptionValue["AddSpinDiffQuadAnsatz"];
addspindiffmixansatz=OptionValue["AddSpinDiffMixAnsatz"];getintervals=OptionValue["GetIntervals"];
spinparameter=OptionValue["SpinParameter"];
spindiffparameter=OptionValue["SpinDiffParameter"];
productAnsatz = OptionValue["ProductAnsatz"];
productGeneralizeOrderNum = OptionValue["ProductGeneralizeOrder"];
productGeneralizeOrderDenom = OptionValue["ProductGeneralizeOrderDenom"];
productGeneralizeMinPowerNum = OptionValue["ProductGeneralizeMinPowerNum"];
If[productGeneralizeOrderDenom<0, productGeneralizeOrderDenom = productGeneralizeOrderNum;];
If[productGeneralizeMinPowerNum<0, productGeneralizeMinPowerNum = If[productAnsatz,0,1];]; (* by default, set the fi0 terms to zero in a sum ansatz, but keep them in a product ansatz; can be user-overriden *)
twoDSpinAnsatzOrderNum = OptionValue["TwoDSpinAnsatzOrderNum"];
twoDSpinAnsatzOrderDenom = OptionValue["TwoDSpinAnsatzOrderDenom"];
fastanalysis = OptionValue["FastAnalysis"];
nsansatz = OptionValue["NonSpinningAnsatzList"];
If[Length@addnsansatz>0,Do[nsansatz=Append[nsansatz,addnsansatz[[i]]];,{i,Length@addnsansatz}]];
unphysical = Select[data,#1[[1]] > 0.25 &];
Print["data points with \[Eta] > 0.25: ", unphysical];
datans = Select[data, #1[[2]]^2+#1[[3]]^2 < 0.01 &]; (* zero spin, for 1d fit *)
dataeqs = Select[data, (Abs[#1[[2]]-#1[[3]]] < 0.01)||(#1[[1]]< OptionValue["\[Eta]LimitEqualSpin"]) &]; (* all equal spins, including zero *)
datauneqs = Complement[data,dataeqs]; (* all unequal spins *)
dataq1 = Select[dataeqs, (Abs[#1[[1]]-0.25] < 0.0001) && (#1[[2]]^2+#1[[3]]^2 >= 0.01) &]; (* equal nonzero spins, equal mass, for 1d fit *)
dataeqsunmass = Select[dataeqs, (Abs[#1[[1]]-0.25] > 0.0001) && (#1[[2]]^2+#1[[3]]^2 >= 0.01) &]; (* equal nonzero spins, unequal mass, for 2d fit *)
dataunconstr = Complement[data,Join[datans,dataq1]]; (* not 1D constrained: unequal-spin equal-mass or nonzero-spin unequal-mass, for 3d fit *)
If[ListQ@tooltiptags,
pos = Flatten[Position[data,#]&/@datans];
tooltiptagsns = tooltiptags[[pos]];
];
Print["#\[NonBreakingSpace]of data points: ", Length@data];
Print["#\[NonBreakingSpace]of non-spinning data points: ", Length@datans];
Print["#\[NonBreakingSpace]of equal-spin equal-mass data points: ", Length@dataq1];
Print["#\[NonBreakingSpace]of equal-spin unequal-mass data points: ", Length@dataeqsunmass];
Print["#\[NonBreakingSpace]of unequal-spin data points: ", Length@datauneqs];
Print["# of non-zero-spin unequal-mass data points: ", Length@dataunconstr];
Print["Sums of subspaces: S=0 + q=1 + equal-spin + unequal-spin cases: ", Length@datans + Length@dataq1 + Length@dataeqsunmass + Length@datauneqs];
Print[" S=0 + q=1 + others: ", Length@datans + Length@dataq1 + Length@dataunconstr];
Print["Eff. spin parameter: ",spinparameter];
Print["Spin-diff parameter: ",spindiffparameter];
(*
Print["Rules for General Anzatze : ", generalAnsatzRules];
Print["Rules for Spin Difference Anzatze : ", spindifferencerules];
*)
allData = TakeColumn[data, {1,2,3,4,5}]/. {\[Eta]\[Eta]_,c1_,c2_,res_,ww_}->{\[Eta]\[Eta],spinparameter[\[Eta]\[Eta],c1,c2],res,ww};
allDataEqS = TakeColumn[dataeqs, {1,2,3,4,5}]/. {\[Eta]\[Eta]_,c1_,c2_,res_,ww_}->{\[Eta]\[Eta],spinparameter[\[Eta]\[Eta],c1,c2],res,ww};
allDataEqSv2 = TakeColumn[dataeqsunmass, {1,2,3,4,5}]/. {\[Eta]\[Eta]_,c1_,c2_,res_,ww_}->{\[Eta]\[Eta],spinparameter[\[Eta]\[Eta],c1,c2],res,ww};
allDataUnconstr = TakeColumn[dataunconstr, {1,2,3,4,5}]/. {\[Eta]\[Eta]_,c1_,c2_,res_,ww_}->{\[Eta]\[Eta],spinparameter[\[Eta]\[Eta],c1,c2],res,ww};
If[(fitcase==3)||(fitcase>=5),
oneDverbosity = 0;
,
oneDverbosity = 2;
];
(* non-spinning fit *)
Print[Style["Non-Spinning Fit",Blue]];
nsfits = DataFitFunction[TakeColumn[datans,{1,4,5}], nsansatz, "Domain"->{0,0.25}, "PlotRange"->plotrange, "StatisticalTest"->statisticaltest, "Verbose"->oneDverbosity, "Weights"->weights,
"ToolTipTags"->tooltiptagsns, "GetIntervals"->getintervals&&(fitcase==1)];
Print["---------------------"];
Print["Select \[Chi]1=\[Chi]2=0 Fit with nsfitindex = ", nsfitindex, " -> ", nsfit = nsfits[[nsfitindex,1]] ];
Print["---------------------"];
(* equal mass/equal spin fit *)
Print[Style["q=1, \[Chi]1=\[Chi]2 Fit",Blue]];
dataq1sorted = Sort[TakeColumn[dataq1,{2,3,4,5}]]/. {xx_?NumberQ,yy_,zz_,ww_}->{spinparameter[0.25,xx,yy],zz,ww};
If[ListQ@tooltiptags,
pos = Flatten[Position[data,#]&/@dataq1sorted];
tooltiptagseq = tooltiptags[[pos]];
];
Clear[b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,S];
(* ansatzq1 = Total/@Table[ToExpression["b"<>ToString@i] * S^i,{j,3,8},{i,1,j}];
ansatzq1 = Transpose[{ansatzq1,Table[{S},{i,Length@ansatzq1}]}]; *)
ansatzq1 = OptionValue["EqualBHAnsatzList"];
If[Length@addq1ansatz>0,Do[ansatzq1=Append[ansatzq1,addq1ansatz[[i]]];,{i,Length@addq1ansatz}]];
If[productAnsatz,
ansatzq1[[All,1]]=(nsfit/. \[Eta]->0.25) * ansatzq1[[All,1]];
,
ansatzq1[[All,1]]=(nsfit/. \[Eta]->0.25) + ansatzq1[[All,1]];
];
If[ListQ@tooltiptags,
pos = Flatten[Position[TakeColumn[data,{2,3,4,5}]/. {xx_?NumberQ,yy_,zz_,ww_}->{spinparameter[0.25,xx,yy],zz,ww},#]&/@dataq1sorted];
tooltiptagseq = tooltiptags[[pos]];
];
q1fits = DataFitFunction[dataq1sorted, ansatzq1, "Domain"-> {-1,1}, "PlotRange"->plotrange, "StatisticalTest"->statisticaltest, "Verbose"->oneDverbosity, "Weights"->weights,
"ToolTipTags"->tooltiptagseq, "GetIntervals"->getintervals&&(fitcase==1)];
If[fitcase==1, Return[{nsfits,q1fits}]];
etaAnsatz = nsfits[[nsfitindex,15]];
etaCoeffRules = nsfits[[nsfitindex,2]];
etaCovar = nsfits[[nsfitindex,14]];
SAnsatz = q1fits[[q1fitindex,15]];
SCoeffRules = q1fits[[q1fitindex,2]];
SCovar = q1fits[[q1fitindex,14]];
q1fit = q1fits[[q1fitindex,1]];
Print["---------------------"];
Print["Select q=1 Fit with q1fitindex = ", q1fitindex, " -> ", q1fit];
Print["---------------------"];
Print[Style["{\[Eta],S} Fit",Blue]];
nsdegree = Exponent[nsfit,\[Eta],List];
q1degree = Exponent[q1fit,S,List];
ansatz={Flatten[Total@Table[Total@Table[ToExpression["f"<>ToString@i<>ToString@j] *S^j,{j,q1degree}]*\[Eta]^i,{i,nsdegree}]],{\[Eta],S}};
If[productAnsatz,
ansatz=Collect[nsfit * Flatten[Total@Table[Total@Table[ToExpression["f"<>ToString@i<>ToString@j] * S^j, {j,4}]*\[Eta]^i, {i,4}]], S];
,
ansatz=Collect[nsfit + Flatten[Total@Table[Total@Table[ToExpression["f"<>ToString@i<>ToString@j] * S^j, {j,4}]*\[Eta]^i, {i,4}]], S];
];
(*If[Length@addgenansatz>0,ansatz=addgenansatz[[1]];];*)
(*ansatzFinal=Table[AnsatzRestrictions[ansatzAll[[i,1]],ansatzAll[[i,2]],ansatzAll[[i,3]]],{i,Length@ansatzAll}];*)
If[fitcase<=4,
Print["\[InvisibleSpace]using nsfit: ", nsfit];
Print["\[InvisibleSpace]using q1fit: ", q1fit];
];
If[Length@addgenansatz>0
,
ansatzFinal=addgenansatz[[1]];
Print[Style["Selecting ansatz from AddGenAnsatz",Blue]]; ,
If[OptionValue["RationalFunctions"],
Print[];
If[productAnsatz,
q1fitS0 = q1fit/(q1fit/. S-> 0);
,
q1fitS0 = q1fit-(q1fit/. S-> 0);
];
q1fitS0 = Chop[q1fitS0]/.{1.->1};
q1fitS0orderNum = Exponent[Numerator[q1fitS0],S];
q1fitS0orderDenom = Exponent[Denominator[q1fitS0],S];
If[twoDSpinAnsatzOrderNum<0, twoDSpinAnsatzOrderNum = q1fitS0orderNum];
If[twoDSpinAnsatzOrderDenom<0, twoDSpinAnsatzOrderDenom = q1fitS0orderDenom];
If[(twoDSpinAnsatzOrderNum>q1fitS0orderNum)||(twoDSpinAnsatzOrderDenom>q1fitS0orderDenom),
spinAnsatzNum = Numerator[q1fitS0] + Sum[1.0*S^i,{i,q1fitS0orderNum+1,twoDSpinAnsatzOrderNum}];
spinAnsatzDenom = Denominator[q1fitS0] + Sum[1.0*S^i,{i,q1fitS0orderDenom+1,twoDSpinAnsatzOrderDenom}];
spinAnsatzToGeneralize = spinAnsatzNum/spinAnsatzDenom;
,
spinAnsatzToGeneralize = q1fitS0;
];
spinAnsatz = GeneralizeFunction[spinAnsatzToGeneralize, "f", \[Eta], productGeneralizeOrderNum, productGeneralizeOrderDenom];
Print["spinAnsatz before boundary conditions = ", spinAnsatz];
(* spinAnsatzRaw = spinAnsatz /. InvertRules[SCoeffRules]; *)
q1AnsatzRaw = q1fits[[q1fitindex,15]];
If[productAnsatz,
q1AnsatzRawS0 = q1AnsatzRaw/(q1AnsatzRaw/.S->0);
,
q1AnsatzRawS0 = q1AnsatzRaw-(q1AnsatzRaw/.S->0);
];
q1AnsatzRawS0 = Chop[q1AnsatzRawS0]/.{1.->1};
If[(twoDSpinAnsatzOrderNum>q1fitS0orderNum)||(twoDSpinAnsatzOrderDenom>q1fitS0orderDenom),
spinAnsatzRawNum = Numerator[q1AnsatzRawS0] + Sum[1.0*S^i,{i,q1fitS0orderNum+1,twoDSpinAnsatzOrderNum}];
spinAnsatzRawDenom = Denominator[q1AnsatzRawS0] + Sum[1.0*S^i,{i,q1fitS0orderDenom+1,twoDSpinAnsatzOrderDenom}];
spinAnsatzRawToGeneralize = spinAnsatzRawNum/spinAnsatzRawDenom;
,
spinAnsatzRawToGeneralize = q1AnsatzRawS0;
];
spinAnsatzRaw = GeneralizeFunction[spinAnsatzRawToGeneralize, "f", \[Eta], productGeneralizeOrderNum, productGeneralizeOrderDenom];
numRestrictions = AnsatzRestrictions[Numerator[q1fitS0],Numerator[spinAnsatz],"Parameters"->True];
denomRestrictions = AnsatzRestrictions[Denominator[q1fitS0],Denominator[spinAnsatz],"Parameters"->True];
num = Collect[Numerator[spinAnsatz]/.numRestrictions,{S,\[Eta]}];
If[productGeneralizeMinPowerNum>0,
zeroRules = Flatten@Table[ToExpression["f"<>ToString@i<>ToString@j]->0,{i,0,Exponent[num,S]},{j,0,productGeneralizeMinPowerNum-1}];
Print["setting the following coefficients to zero:"];
Print[zeroRules];
num = num/.zeroRules;
];
denom = Collect[Denominator[spinAnsatz]/.denomRestrictions,{S,\[Eta]}];
spinpart = num/denom;
Print["spinAnsatz after equal mass boundary conditions = ",spinpart];
Print["applying user-set generalAnsatzRules:"];
Print[generalAnsatzRules];
If[productAnsatz,
ansatzFinal = (nsfit * spinpart) /. generalAnsatzRules;
(* \[Eta]0Ansatz =Limit[(spinAnsatz/\[Eta]), \[Eta]\[Rule] 0]//Chop//Simplify; *)
ansatzRaw = nsfit * spinAnsatzRaw;
,
ansatzFinal = (nsfit + spinpart) /. generalAnsatzRules;
(* \[Eta]0Ansatz =Limit[(spinpart/\[Eta]), \[Eta]\[Rule] 0]//Chop//Simplify; *)
ansatzRaw = nsfit + spinAnsatzRaw;];
Print[];
(* 1 2 3 4 5 6 7 8 9 *)
extraStuff = {ansatzRaw, etaAnsatz, etaCoeffRules, etaCovar, SAnsatz, SCoeffRules, SCovar, numRestrictions, denomRestrictions};
,
(* Constrain the ansatz to non spinning and q=1 cases *)
(*ansatzFinal=AnsatzRestrictions[nsfit,q1fit,ansatz] /. generalAnsatzRules;*)
ansatzFinal=AnsatzRestrictions[ nsfit,q1fit,ansatz];
];
];
Print["\[InvisibleSpace]ansatzFinal: ", ansatzFinal];
Print[];
(*
(* this is redundant and should only differ from 'v2' in uncertainty interval widths *)
If[fitcase\[LessEqual]4,
Print[Style["eqS fit to all "<>ToString@Length@allData<>" data points:",Blue]];
alldataeqfit = DataFitFunction[allData, {{ansatzFinal,{\[Eta],S}}}, "PlotRange"\[Rule]{{0,0.25},{-1,1},All}, "Domain"\[Rule]{{0,0.25},{-1,1}},
"Verbose"\[Rule]2, "Weights"\[Rule]weights,"StatisticalTest"\[Rule]statisticaltest, "GetIntervals"\[Rule]getintervals];
];
*)
Print[Style["eqS fit to all "<>ToString@Length@allDataUnconstr<>" data points without 1D regions:",Blue]];
alldataeqfitv2 = DataFitFunction[allDataUnconstr, {{ansatzFinal,{\[Eta],S}}}, "PlotRange"->{{0,0.25},{-1,1},All}, "Domain"->{{0,0.25},{-1,1}},
"Verbose"->2, "Weights"->weights,"StatisticalTest"->statisticaltest, "GetIntervals"->getintervals];
(*
(* this is redundant and should only differ from 'v2' in uncertainty interval widths *)
If[fitcase\[LessEqual]4,
Print[Style["eqS fit to "<>ToString@Length@allDataEqS<>" equal-spin data points:",Blue]];
eqdatafit = DataFitFunction[allDataEqS, {{ansatzFinal,{\[Eta],S}}}, "PlotRange"\[Rule]{{0,0.25},{-1,1},All}, "Domain"\[Rule]{{0,0.25},{-1,1}},
"Verbose"\[Rule]2, "Weights"\[Rule]weights,"StatisticalTest"\[Rule]statisticaltest, "GetIntervals"\[Rule]getintervals];
];
*)
Print[Style["eqS fit to "<>ToString@Length@allDataEqSv2<>" equal-spin data points without 1D regions:",Blue]];
eqdatafitv2 = DataFitFunction[allDataEqSv2, {{ansatzFinal,{\[Eta],S}}}, "PlotRange"->{{0,0.25},{-1,1},All}, "Domain"->{{0,0.25},{-1,1}},
"Verbose"->2, "Weights"->weights,"StatisticalTest"->statisticaltest, "GetIntervals"->getintervals];
(* 2DeqS 2Dall 3D placeholders 1D2Dparts 1Deta 1DS 2DSansatz *)
If[fitcase<=3, Return[{eqdatafitv2, alldataeqfitv2, {}, {}, {}, {}, extraStuff, nsfits, q1fits, spinpart}];];
Print[Style["(fit to allData) - (fit to allDataEqS):",Blue]];
Print[Plot3D[alldataeqfitv2[[1,1]]-eqdatafitv2[[1,1]],{\[Eta],0,0.25`},{S,-1,1},AxesLabel->{Style["\[Eta]",Black,14],Style["S",Black,14],Style["diff",14]}]];
Clear[a0,a1,a2,a3,\[Chi]diff,q];
(* Fit the spin difference terms *)
Print["----------------"];
Print[Style["{\[Eta],S,\[Chi]diff}",Blue]];
mydq=0.005;
tabfits=Table[
Print[Style["MassRatio -> ",Blue],q\[Eta][i]];
my\[Eta]=i;
thisq=q\[Eta][i];
If[ thisq < 2,
myd\[Eta] = Abs[0.0001+mydq (-2 thisq/(1+thisq)^3 + 1./(1+thisq)^2) ];,
myd\[Eta] = 0.001;
];
Print["myd\[Eta] = ", myd\[Eta] ];
mysel = Select[data,my\[Eta]-myd\[Eta]<#1[[1]]<my\[Eta]+myd\[Eta]&];
dataq1 = TakeColumn[mysel,{2,3,4,5}]/. {xx_?NumberQ,yy_,zz_,ww_}->{xx,yy,eqdatafitv2[[1,1]]/. \[Eta]->my\[Eta]/. S->spinparameter[my\[Eta],xx,yy],ww};
avgwgthperq = Mean[dataq1[[All,-1]]];
If[fitcase>=5,
perqVerbosity = 0; (* for higher fit cases, still do the per-q analysis, but suppress the plots *)
,
perqVerbosity = 1;
If[ListQ@tooltiptags, dataq1pos = Flatten[Position[data,#]&/@mysel];
tooltiptagsaux = tooltiptags[[dataq1pos]];
];
If[0.25 -myd\[Eta]<my\[Eta]<0.25 +myd\[Eta],
dataq1Sym = dataq1/. {xx_?NumberQ,yy_,zz_,ww_}->{yy,xx,zz,ww};
dataq1 = Join[dataq1,dataq1Sym];
];
plot1 = myListPlot3D[dataq1[[All,1;;3]], AxesLabel->{Style["\[Chi]1",14],Style["\[Chi]2",14],""}, PlotRange->{{-1,1},{-1,1},All}, PlotLabel->Style["data(\[Chi]1,\[Chi]2)",14], ImageSize->250];
dataq1 = TakeColumn[mysel,{2,3,4,5}]/. {xx_?NumberQ,yy_,zz_,ww_}->{xx,yy,zz-(eqdatafitv2[[1,1]]/. \[Eta]->my\[Eta]/. S->spinparameter[my\[Eta],xx,yy]),ww};
If[0.25 -myd\[Eta]<my\[Eta]<0.25 +myd\[Eta],
dataq1Sym=dataq1/. {xx_?NumberQ,yy_,zz_,ww_}->{yy,xx,zz,ww};
dataq1=Join[dataq1,dataq1Sym];
];
plot2 = myListPlot3D[dataq1[[All,1;;3]], AxesLabel->{Style["\[Chi]1",14],Style["\[Chi]2",14],""}, PlotRange->{{-1,1},{-1,1},All}, PlotLabel->Style["data(\[Chi]1,\[Chi]2) - fit(S)",14], ImageSize->250];
dataq1 = TakeColumn[mysel,{2,3,4,5}]/. {xx_?NumberQ,yy_,zz_,ww_}->{spinparameter[my\[Eta],xx,yy],xx-yy,zz-(eqdatafitv2[[1,1]]/. \[Eta]->my\[Eta]/. S->spinparameter[my\[Eta],xx,yy]),ww};
If[0.25 -myd\[Eta]<my\[Eta]<0.25+myd\[Eta],
dataq1Sym=dataq1/. {xx_?NumberQ,yy_,zz_,ww_}->{xx,-yy,zz,ww};
dataq1=Join[dataq1,dataq1Sym];
];
plot3 = myListPlot3D[dataq1[[All,1;;3]], PlotStyle->Directive[RGBColor[0.880722, 0.611041, 0.142051],Opacity[0.7]], Lighting->{{"Ambient",RGBColor[0.880722, 0.611041, 0.142051]}},
AxesLabel->{Style["S",14],Style["\!\(\*SubscriptBox[\(\[Chi]\), \(diff\)]\)",14],""}, PlotRange->{{-1,1},{-2,2},All}, PlotLabel->Style["data(S,\!\(\*SubscriptBox[\(\[Chi]\), \(diff\)]\)) - fit(S)",14], ImageSize->250];
];
Clear[a,b,c,\[Chi]1,\[Chi]2,\[Delta]\[Chi]];
dataq1a= TakeColumn[dataq1,{1,2,3}];
dataq1 = TakeColumn[mysel,{2,3,4,5}]/. {xx_?NumberQ,yy_,zz_,ww_}->{spindiffparameter[my\[Eta],xx,yy],zz-(eqdatafitv2[[1,1]]/. \[Eta]->my\[Eta]/. S->spinparameter[my\[Eta],xx,yy]),ww};
(*dataq1 = TakeColumn[mysel,{2,3,4,5}]/.\[VeryThinSpace]{xx_?NumberQ,yy_,zz_,ww_}\[Rule]{(xx-yy),zz-(eqdatafitv2[[1,1]]/.\[VeryThinSpace]\[Eta]\[Rule]my\[Eta]/.\[VeryThinSpace]S\[Rule]spinparameter[my\[Eta],xx,yy]),ww}*)
Print["data dimensionality: ",Length@dataq1];
If[0.25 -myd\[Eta]<my\[Eta]<0.25+myd\[Eta],
dataq1Sym = dataq1/. {xx_?NumberQ,zz_,ww_}->{-xx,zz,ww};
dataq1 = Join[dataq1,dataq1Sym];
If[perqVerbosity>0, tooltiptagsaux = Join[tooltiptagsaux,tooltiptagsaux];];
Print["data dimensionality (including symmetric duplicates): ",Length@dataq1];
];
(*dataq1mixterm = TakeColumn[mysel,{2,3,4,5}]/.\[VeryThinSpace]{xx_?NumberQ,yy_,zz_,ww_}\[Rule]{xx-yy,spinparameter[my\[Eta],xx,yy]*(xx-yy),zz-(eqdatafitv2[[1,1]]/.\[VeryThinSpace]\[Eta]\[Rule]my\[Eta]/.\[VeryThinSpace]S\[Rule]spinparameter[my\[Eta],xx,yy]),ww};*)
dataq1mixterm = TakeColumn[mysel,{2,3,4,5}]/. {xx_?NumberQ,yy_,zz_,ww_}->{spindiffparameter[my\[Eta],xx,yy],spinparameter[my\[Eta],xx,yy]*spindiffparameter[my\[Eta],xx,yy],zz-(eqdatafitv2[[1,1]]/. \[Eta]->my\[Eta]/. S->spinparameter[my\[Eta],xx,yy]),ww};
myfit = DataFitFunction[dataq1, {{a \[Delta]\[Chi], {\[Delta]\[Chi]}}, {a \[Delta]\[Chi] + b \[Delta]\[Chi]^2,{\[Delta]\[Chi]}}}, "Verbose"->perqVerbosity, "AxesTag"->"Amplitude",
"StatisticalTest"->statisticaltest, "Sorted"->False, "Weights"->weights, "GetIntervals"->False];
myfitmixterm = DataFitFunction[dataq1mixterm, {{a \[Delta]\[Chi] + c y, {\[Delta]\[Chi],y}}, {a \[Delta]\[Chi] + c y +b \[Delta]\[Chi]^2,{\[Delta]\[Chi],y}}}, "Verbose"->perqVerbosity, "AxesTag"->"Amplitude",
"StatisticalTest"->statisticaltest, "Sorted"->False, "Weights"->weights, "GetIntervals"->False];
myfit = Join[myfit,myfitmixterm,{avgwgthperq}];
(* do NOT show mix-term fits here, as they tend to blow the scale, and are only projected down anyway *)
(* plot4=Show[Plot[{myfit[[1,1]],myfit[[2,1]]},{\[Delta]\[Chi],-0.1,0.1},PlotRange\[Rule]All,ImageSize\[Rule]250,
Frame\[Rule]True,FrameLabel\[Rule]{"Subscript[\[Chi], diff]","\[CapitalDelta](Subscript[\[Chi], diff])"},PlotLegends\[Rule]{"linear","lin+quad"}],
ListPlot[Table[Tooltip[dataq1[[i,1;;2]],If[ListQ@tooltiptags,tooltiptagsaux[[i]],ToString[{dataq1[[i,1]],dataq1[[i,2]]}]]],{i,Length@dataq1}],PlotStyle\[Rule]PointSize[Medium],PlotRange\[Rule]All,Frame->True]];
*)
dataq1a= Transpose[{dataq1a[[All,1]],dataq1a[[All,2]],dataq1a[[All,3]]-(myfit[[1,1]]/.\[Delta]\[Chi]->dataq1[[All,1]])}];
plot3a = myListPlot3D[dataq1a, PlotStyle->Directive[RGBColor[0.880722, 0.611041, 0.142051],Opacity[0.7]], Lighting->{{"Ambient",RGBColor[0.880722, 0.611041, 0.142051]}},
AxesLabel->{Style["S",14],Style["\!\(\*SubscriptBox[\(\[Chi]\), \(diff\)]\)",14],""}, PlotRange->{{-1,1},{-2,2},All}, PlotLabel->Style["data(S,\!\(\*SubscriptBox[\(\[Chi]\), \(diff\)]\)) - (fit(S) + f(\[Eta])\!\(\*SubscriptBox[\(\[Chi]\), \(diff\)]\))",14], ImageSize->250];
If[fitcase==4,
plot4a = ListPlot[Table[Tooltip[dataq1[[i,1;;2]],If[ListQ@tooltiptags,tooltiptagsaux[[i]],ToString[{dataq1[[i,1]],dataq1[[i,2]]}]]],{i,Length@dataq1}],
PlotStyle->PointSize[Medium], PlotRange->{{-2.,2.},All}, PlotLegends->{"NR data"}];
plot4b = Plot[{myfit[[1,1]],myfit[[2,1]]},{\[Delta]\[Chi],-2.,2.}, PlotLegends->{"linear","lin+quad"}, PlotRange->{{-2.,2.},All}];
If[0.25 -myd\[Eta]<my\[Eta]<0.25+myd\[Eta],
plot4mirror = ListPlot[Table[{dataq1[[i,1]],dataq1[[i,2]]},{i,Length@dataq1}], PlotStyle->{Gray,PointSize[Medium]}, PlotMarkers->\!\(\*
StyleBox["\"\<\[FilledSquare]\>\"",
StripOnInput->False,
FontSize->8]\), PlotLegends->{"mirror NR data"}];
plot4 = Show[{plot4b,plot4a}, Frame->True, FrameLabel->{"\!\(\*SubscriptBox[\(\[Chi]\), \(diff\)]\)","\[CapitalDelta](\!\(\*SubscriptBox[\(\[Chi]\), \(diff\)]\))"}, ImageSize->250]; ,
plot4 = Show[{plot4b,plot4a}, Frame->True, FrameLabel->{"\!\(\*SubscriptBox[\(\[Chi]\), \(diff\)]\)","\[CapitalDelta](\!\(\*SubscriptBox[\(\[Chi]\), \(diff\)]\))"}, ImageSize->250];
];
Print[{plot1,plot2,plot3,plot3a,plot4}];
];
myfit, {i,massratiounequalfits}]; (* this line finishes the tabfits list *)
(*Print[tabfits[[All,-1]]];*)
Print["-----------"];
Print[Style["\[Eta] dependence of \!\(\*SubscriptBox[\(\[Chi]\), \(diff\)]\) terms",Blue]];
extrawght=Join[tabfits[[All,-1]],{1.}];
If[fastanalysis,
If[Length@addspindifflinansatz>0,
ansatzchidifflinear = addspindifflinansatz[[1,1]];
,
ansatzchidifflinear = (ToExpression["a10"] \[Eta]^Abs[ToExpression["a11"]] Sqrt[1-4 \[Eta]]^Abs[ToExpression["a12"]] ) /. spindifferencerules;
];
If[Length@addspindiffquadansatz>0,
ansatzchidiffquadratic = addspindiffquadansatz[[1,1]];
,
ansatzchidiffquadratic = ((ToExpression["a20"] \[Eta]^Abs[ToExpression["a21"]] Sqrt[1-4 \[Eta]]^Abs[ToExpression["a22"]] )) /. spindifferencerules;
];
ansatzFinal\[Delta]\[Chi]2 = (ansatzFinal /. S->spinparameter[\[Eta],\[Chi]1,\[Chi]2]) + ansatzchidifflinear spindiffparameter[\[Eta],\[Chi]1,\[Chi]2] + ansatzchidiffquadratic spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]^2;
If[Length@addspindiffmixansatz>0,
ansatzchidiffmix = addspindiffmixansatz[[1,1]];
,
ansatzchidiffmix = ((ToExpression["a30"] \[Eta]^Abs[ToExpression["a31"]] Sqrt[1-4 \[Eta]]^Abs[ToExpression["a32"]] )) /. spindifferencerules;
];
ansatzS\[Delta]\[Chi]\[Chi]2raw = ansatzchidifflinear spindiffparameter[\[Eta],\[Chi]1,\[Chi]2] + ansatzchidiffmix S spindiffparameter[\[Eta],\[Chi]1,\[Chi]2] + ansatzchidiffquadratic spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]^2; (* needed to extract the coefficients later on *)
ansatzFinalS\[Delta]\[Chi]\[Chi]2 = ( ansatzFinal + ansatzS\[Delta]\[Chi]\[Chi]2raw )/.S->spinparameter[\[Eta],\[Chi]1,\[Chi]2]; (* this is the form passed to DataFitFunction[] *)
Print[Style["full 3D fit to all "<>ToString@Length@dataunconstr<>" data points without 1D regions... (see results at the end after plots)",Blue]];
uneqfit2 = DataFitFunction[dataunconstr, {{ansatzFinalS\[Delta]\[Chi]\[Chi]2,{\[Eta],\[Chi]1,\[Chi]2}}},
"Verbose"->0, "AxesTag"->"Amplitude", "StatisticalTest"->statisticaltest, "Sorted"->False, "Weights"->weights, "GetIntervals"->getintervals];
Print["Done"];
Return[{uneqfit2 ,extraStuff, nsfits, q1fits, spinpart}];
,
linearparerr=Flatten@TakeColumn[tabfits,1][[All,9,1]];
linearparerr=Join[linearparerr,{10^-9}];
linearpar=Transpose[{Join[massratiounequalfits,{0}],Join[a/.TakeColumn[tabfits,1][[All,2]],{0}],extrawght/linearparerr^2}];
linearpar[[1]]={0.25,0,extrawght[[1]]/(10^(-9))^2};
linearparerr[[1]]={10^-9};
linearpar2err=Flatten@TakeColumn[tabfits,2][[All,9,1]];
linearpar2err=Join[linearpar2err,{10^-9}];
linearpar2=Transpose[{Join[massratiounequalfits,{0}],Join[a/.TakeColumn[tabfits,2][[All,2]],{0}],extrawght/linearpar2err^2}];
linearpar2[[1]]={0.25,0,extrawght[[1]]/(10^(-9))^2};
linearpar2err[[1]]={10^-9};
quadparerr=Flatten@TakeColumn[tabfits,2][[All,9,2]];
quadparerr=Join[quadparerr,{10^-9}];
quadpar=Transpose[{Join[massratiounequalfits,{0}],Join[b/.TakeColumn[tabfits,2][[All,2]],{0}],extrawght/quadparerr^2}];
linearpar3err=Flatten@TakeColumn[tabfits,3][[All,9,1]];
linearpar3err=Join[linearpar3err,{10^-9}];
linearpar3=Transpose[{Join[massratiounequalfits,{0}],Join[a/.TakeColumn[tabfits,3][[All,2]],{0}],extrawght/linearpar2err^2}];
linearpar3[[1]]={0.25,0,extrawght[[1]]/(10^(-9))^2};linearpar3err[[1]]={10^-9};
mixparerr=Flatten@TakeColumn[tabfits,3][[All,9,2]];
mixparerr=Join[mixparerr,{10^-9}];
mixpar=Transpose[{Join[massratiounequalfits,{0}],Join[c/.TakeColumn[tabfits,3][[All,2]],{0}],extrawght/mixparerr^2}];
mixpar[[1]]={0.25,0,extrawght[[1]]/(10^(-9))^2};
mixparerr[[1]]={10^-9};
linearpar4err=Flatten@TakeColumn[tabfits,4][[All,9,1]];
linearpar4err=Join[linearpar4err,{10^-9}];
linearpar4=Transpose[{Join[massratiounequalfits,{0}],Join[a/.TakeColumn[tabfits,4][[All,2]],{0}],extrawght/linearpar4err^2}];
linearpar4[[1]]={0.25,0,extrawght[[1]]/(10^(-9))^2};
linearpar4err[[1]]={10^-9};
mixparerr2=Flatten@TakeColumn[tabfits,4][[All,9,2]];
mixparerr2=Join[mixparerr2,{10^-9}];
mixpar2=Transpose[{Join[massratiounequalfits,{0}],Join[c/.TakeColumn[tabfits,4][[All,2]],{0}],extrawght/mixparerr2^2}];
mixpar2[[1]]={0.25,0,extrawght[[1]]/(10^(-9))^2};
mixparerr2[[1]]={10^-9};
quadparerr2=Flatten@TakeColumn[tabfits,4][[All,9,2]];
quadparerr2=Join[quadparerr2,{10^-9}];
quadpar2=Transpose[{Join[massratiounequalfits,{0}],Join[b/.TakeColumn[tabfits,4][[All,2]],{0}],extrawght/quadparerr2^2}];
(*ansatz\[Eta]1=a0 \[Eta]^Abs[a1] N[1-4 \[Eta]]^0.7;*)
(*ansatz\[Eta]1=a0 \[Eta]^Abs[a1] N[1-4 \[Eta]]^Abs[a2];*)
(*ansatz\[Eta]1=a0 \[Eta] N[1-4 \[Eta]]^Abs[a2];*)
If[Length@addspindifflinansatz>0,
ansatzchidifflinear = addspindifflinansatz[[1,1]];
,
ansatzchidifflinear = (ToExpression["a10"] \[Eta]^Abs[ToExpression["a11"]] Sqrt[1-4 \[Eta]]^Abs[ToExpression["a12"]] ) /. spindifferencerules;
];
(*ansatzFinal\[Delta]\[Chi]=(ansatzFinal/. S\[Rule] spinparameter[\[Eta],\[Chi]1,\[Chi]2])+ansatzchidifflinear (\[Chi]1-\[Chi]2);*)
ansatzFinal\[Delta]\[Chi] = (ansatzFinal /. S->spinparameter[\[Eta],\[Chi]1,\[Chi]2]) + ansatzchidifflinear spindiffparameter[\[Eta],\[Chi]1,\[Chi]2];
If[Length@addspindiffquadansatz>0,
ansatzchidiffquadratic = addspindiffquadansatz[[1,1]];
,
ansatzchidiffquadratic = ((ToExpression["a20"] \[Eta]^Abs[ToExpression["a21"]] Sqrt[1-4 \[Eta]]^Abs[ToExpression["a22"]] )) /. spindifferencerules;
];
ansatzFinal\[Delta]\[Chi]2 = (ansatzFinal /. S->spinparameter[\[Eta],\[Chi]1,\[Chi]2]) + ansatzchidifflinear spindiffparameter[\[Eta],\[Chi]1,\[Chi]2] + ansatzchidiffquadratic spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]^2;
If[Length@addspindiffmixansatz>0,
ansatzchidiffmix = addspindiffmixansatz[[1,1]];
,
ansatzchidiffmix = ((ToExpression["a30"] \[Eta]^Abs[ToExpression["a31"]] Sqrt[1-4 \[Eta]]^Abs[ToExpression["a32"]] )) /. spindifferencerules;
];
ansatzS\[Delta]\[Chi]raw = ansatzchidifflinear spindiffparameter[\[Eta],\[Chi]1,\[Chi]2] + ansatzchidiffmix S spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]; (* needed to extract the coefficients later on *)
ansatzFinalS\[Delta]\[Chi] = ( ansatzFinal + ansatzS\[Delta]\[Chi]raw ) /. S->spinparameter[\[Eta],\[Chi]1,\[Chi]2]; (* this is the form passed to DataFitFunction[] *)
ansatzS\[Delta]\[Chi]\[Chi]2raw = ansatzchidifflinear spindiffparameter[\[Eta],\[Chi]1,\[Chi]2] + ansatzchidiffmix S spindiffparameter[\[Eta],\[Chi]1,\[Chi]2] + ansatzchidiffquadratic spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]^2; (* needed to extract the coefficients later on *)
ansatzFinalS\[Delta]\[Chi]\[Chi]2 = ( ansatzFinal + ansatzS\[Delta]\[Chi]\[Chi]2raw ) /. S->spinparameter[\[Eta],\[Chi]1,\[Chi]2]; (* this is the form passed to DataFitFunction[] *)
Print["Fitting ansaetze to per-q slopes, linear only:"];
linearfit1 = DataFitFunction[linearpar, {{ansatzchidifflinear,{\[Eta]}}}, "Verbose"->1, "Weights"->weights];
Print["linear+quadratic:"];
linearfit2 = DataFitFunction[linearpar2, {{ansatzchidifflinear,{\[Eta]}}}, "Verbose"->1, "Weights"->weights];
quadfit = DataFitFunction[quadpar, {{ansatzchidiffquadratic,{\[Eta]}}}, "Verbose"->1, "Weights"->weights];
Print["linear+mixture:"];
linearfit3 = DataFitFunction[linearpar3, {{ansatzchidifflinear,{\[Eta]}}}, "Verbose"->1, "Weights"->weights];
mixfit = DataFitFunction[mixpar, {{ansatzchidiffmix,{\[Eta]}}}, "Verbose"->1, "Weights"->weights];
Print["linear+mixture+quadratic:"];
linearfit4 = DataFitFunction[linearpar4, {{ansatzchidifflinear,{\[Eta]}}}, "Verbose"->1, "Weights"->weights];
mixfit2 = DataFitFunction[mixpar2, {{ansatzchidiffmix,{\[Eta]}}}, "Verbose"->1, "Weights"->weights];
quadfit2 = DataFitFunction[quadpar2, {{ansatzchidiffquadratic,{\[Eta]}}}, "Verbose"->1, "Weights"->weights];
Clear[a10,a11,a12,a20,a21,a22,a30,a31,a32,a1,a2,a0];
(*If[Length@addgenansatz>0,ansatz=addgenansatz[[1]];];*)
(*ansatzFinal=Table[AnsatzRestrictions[ansatzAll[[i,1]],ansatzAll[[i,2]],ansatzAll[[i,3]]],{i,Length@ansatzAll}];*)
(* this will only differ from uneqfit2 in the width of the error bars, not in the results, and is not needed anymore, but let's keep it here commented-out to have it handy for sanity checks *)
(*
Print[Style["full 3D fit to all "<>ToString@Length@data<>" data points...",Blue]];
uneqfit = DataFitFunction[data, {{ansatzFinal\[Delta]\[Chi],{\[Eta],\[Chi]1,\[Chi]2}},{ansatzFinal\[Delta]\[Chi]2,{\[Eta],\[Chi]1,\[Chi]2}},{ansatzFinalS\[Delta]\[Chi],{\[Eta],\[Chi]1,\[Chi]2}},{ansatzFinalS\[Delta]\[Chi]\[Chi]2,{\[Eta],\[Chi]1,\[Chi]2}}},
"Verbose"\[Rule]0, "AxesTag"\[Rule]"Amplitude", "StatisticalTest"\[Rule]statisticaltest, "Sorted"\[Rule]False, "Weights"\[Rule] weights, "GetIntervals"\[Rule]getintervals];
*)
Print[Style["full 3D fit to all "<>ToString@Length@dataunconstr<>" data points without 1D regions... (see results at the end after plots)",Blue]];
uneqfit2 = DataFitFunction[dataunconstr, {{ansatzFinal\[Delta]\[Chi],{\[Eta],\[Chi]1,\[Chi]2}},{ansatzFinal\[Delta]\[Chi]2,{\[Eta],\[Chi]1,\[Chi]2}},{ansatzFinalS\[Delta]\[Chi],{\[Eta],\[Chi]1,\[Chi]2}},{ansatzFinalS\[Delta]\[Chi]\[Chi]2,{\[Eta],\[Chi]1,\[Chi]2}}},
"Verbose"->0, "AxesTag"->"Amplitude", "StatisticalTest"->statisticaltest, "Sorted"->False, "Weights"->weights, "GetIntervals"->getintervals];
Print[Style["spin-diff fit to "<>ToString@Length@datauneqs<>" unequal-spin data points...",Blue]];
uneqfit3 = DataFitFunction[TakeColumn[datauneqs,{1,2,3,4,5}]/. {\[Eta]\[Eta]_,c1_,c2_,res_,ww_}->{\[Eta]\[Eta],c1,c2,res-(eqdatafitv2[[1,1]]/. \[Eta]->\[Eta]\[Eta]/. S->spinparameter[\[Eta]\[Eta],c1,c2]),ww},
{{ansatzchidifflinear spindiffparameter[\[Eta],\[Chi]1,\[Chi]2], {\[Eta],\[Chi]1,\[Chi]2}},
{ansatzchidifflinear spindiffparameter[\[Eta],\[Chi]1,\[Chi]2] + ansatzchidiffquadratic spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]^2, {\[Eta],\[Chi]1,\[Chi]2}},
{ansatzchidifflinear spindiffparameter[\[Eta],\[Chi]1,\[Chi]2] + ansatzchidiffmix spinparameter[\[Eta],\[Chi]1,\[Chi]2] spindiffparameter[\[Eta],\[Chi]1,\[Chi]2], {\[Eta],\[Chi]1,\[Chi]2}},
{ansatzchidifflinear spindiffparameter[\[Eta],\[Chi]1,\[Chi]2] + ansatzchidiffmix spinparameter[\[Eta],\[Chi]1,\[Chi]2] spindiffparameter[\[Eta],\[Chi]1,\[Chi]2] + ansatzchidiffquadratic spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]^2, {\[Eta],\[Chi]1,\[Chi]2}}},
"Verbose"->1, "AxesTag"->"Amplitude", "StatisticalTest"->statisticaltest, "Sorted"->False, "Weights"->weights, "GetIntervals"->getintervals];
(*
Print[Style["full 3D fit to "<>ToString@Length@datauneqs<>" unequal-spin data points...",Blue]];
uneqfit3 = DataFitFunction[datauneqs, {{ansatzFinal\[Delta]\[Chi],{\[Eta],\[Chi]1,\[Chi]2}},{ansatzFinal\[Delta]\[Chi]2,{\[Eta],\[Chi]1,\[Chi]2}},{ansatzFinalS\[Delta]\[Chi],{\[Eta],\[Chi]1,\[Chi]2}},{ansatzFinalS\[Delta]\[Chi]\[Chi]2,{\[Eta],\[Chi]1,\[Chi]2}}},
"Verbose"\[Rule]0, "AxesTag"\[Rule]"Amplitude", "StatisticalTest"\[Rule]statisticaltest, "Sorted"\[Rule]False, "Weights"\[Rule] weights, "GetIntervals"\[Rule]getintervals];
*)
(* eta-dependence for direct 3D fit over all non-1D data *)
spindiffterm = CoefficientList[uneqfit2[[1,1]]/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[2]];
spindiffterm21 = CoefficientList[uneqfit2[[2,1]]/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[2]];
spindiffterm22 = CoefficientList[uneqfit2[[2,1]]/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[3]];
spindiffterm31 = CoefficientList[ansatzS\[Delta]\[Chi]raw/.uneqfit2[[3,2]]/.S->0/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[2]];
spindiffterm32 = CoefficientList[ansatzS\[Delta]\[Chi]raw/.uneqfit2[[3,2]]/.S*spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->y,y][[2]];
spindiffterm41 = CoefficientList[ansatzS\[Delta]\[Chi]\[Chi]2raw/.uneqfit2[[4,2]]/.S->0/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[2]];
spindiffterm42 = CoefficientList[ansatzS\[Delta]\[Chi]\[Chi]2raw/.uneqfit2[[4,2]]/.S spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->y,y][[2]];
spindiffterm43 = CoefficientList[ansatzS\[Delta]\[Chi]\[Chi]2raw/.uneqfit2[[4,2]]/.S->0/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[3]];
(* eta-dependence for residuals-only spin-diff fit over all uneqS data *)
dfitspindiffterm = CoefficientList[uneqfit3[[1,1]]/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[2]];
dfitspindiffterm21 = CoefficientList[uneqfit3[[2,1]]/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[2]];
dfitspindiffterm22 = CoefficientList[uneqfit3[[2,1]]/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[3]];
dfitspindiffterm31 = CoefficientList[ansatzS\[Delta]\[Chi]raw/.uneqfit3[[3,2]]/.S->0/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[2]];
dfitspindiffterm32 = CoefficientList[ansatzS\[Delta]\[Chi]raw/.uneqfit3[[3,2]]/.S*spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->y,y][[2]];
dfitspindiffterm41 = CoefficientList[ansatzS\[Delta]\[Chi]\[Chi]2raw/.uneqfit3[[4,2]]/.S->0/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[2]];
dfitspindiffterm42 = CoefficientList[ansatzS\[Delta]\[Chi]\[Chi]2raw/.uneqfit3[[4,2]]/.S spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->y,y][[2]];
dfitspindiffterm43 = CoefficientList[ansatzS\[Delta]\[Chi]\[Chi]2raw/.uneqfit3[[4,2]]/.S->0/.spindiffparameter[\[Eta],\[Chi]1,\[Chi]2]->\[Chi]diff,\[Chi]diff][[3]];
lplot=Show[{Plot[{spindiffterm,dfitspindiffterm,linearfit1[[1,1]]}, {\[Eta],0,0.25}, PlotLegends->{"direct 3D fit","fit to residuals","per-q res. fits"}, PlotStyle->{{Orange},{Magenta,Dotted},{Blue,Dashed}}],
ErrorListPlot[Table[{linearpar[[i,1;;2]], ErrorBar[1/Sqrt[linearpar[[i,-1]]]]}, {i,Length@linearparerr}], PlotStyle->Red],
ErrorListPlot[Table[{linearpar[[i,1;;2]], ErrorBar[linearparerr[[i]]]}, {i,Length@linearparerr}], PlotLegends->{"per-q fit results"}]},
PlotRange->{{0,0.25},{Min@linearpar[[All,2]]-0.001,Max@linearpar[[All,2]]+0.001}}, ImageSize->330, Frame->True, FrameLabel->{Style["\[Eta]",14],Style["\!\(\*SubscriptBox[\(f\), \(a\)]\)(\[Eta])",14]}, PlotLabel->{"linear fit, \!\(\*SubscriptBox[\(f\), \(a\)]\)(\[Eta])="{ansatzchidifflinear}}];
lqplotl=Show[{Plot[{spindiffterm21,dfitspindiffterm21,linearfit2[[1,1]]}, {\[Eta],0,0.25}, PlotLegends->{"direct 3D fit","fit to residuals","per-q res. fits"}, PlotStyle->{{Orange},{Magenta,Dotted},{Blue,Dashed}}],
ErrorListPlot[Table[{linearpar2[[i,1;;2]], ErrorBar[1/Sqrt[linearpar2[[i,-1]]]]}, {i,Length@linearpar2err}], PlotStyle->Red],
ErrorListPlot[Table[{linearpar2[[i,1;;2]], ErrorBar[linearpar2err[[i]]]}, {i,Length@linearpar2err}], PlotLegends->{"per-q fit results"}]},
PlotRange->{{0,0.25},{Min@linearpar2[[All,2]]-0.001,Max@linearpar2[[All,2]]+0.001}}, ImageSize->330, Frame->True, FrameLabel->{Style["\[Eta]",14],Style["\!\(\*SubscriptBox[\(f\), \(a\)]\)(\[Eta])",14]}, PlotLabel->{"lin+quad fit, lin term, \!\(\*SubscriptBox[\(f\), \(a\)]\)(\[Eta])="{ansatzchidifflinear}}];
lqplotq=Show[{Plot[{spindiffterm22,dfitspindiffterm22,quadfit[[1,1]]}, {\[Eta],0,0.25}, PlotLegends->{"direct 3D fit","fit to residuals","per-q res. fits"}, PlotStyle->{{Orange},{Magenta,Dotted},{Blue,Dashed}}],
ErrorListPlot[Table[{quadpar[[i,1;;2]],ErrorBar[1/Sqrt[quadpar[[i,-1]]]]},{i,Length@quadparerr}], PlotStyle->Red],
ErrorListPlot[Table[{quadpar[[i,1;;2]],ErrorBar[quadparerr[[i]]]},{i,Length@quadparerr}], PlotLegends->{"per-q fit results"}]},
PlotRange->{{0,0.25},{Min@quadpar[[All,2]]-0.001,Max@quadpar[[All,2]]+0.001}}, ImageSize->330, Frame->True, FrameLabel->{Style["\[Eta]",14],Style["\!\(\*SubscriptBox[\(f\), \(b\)]\)(\[Eta])",14]}, PlotLabel->{"lin+quad fit, quad term, \!\(\*SubscriptBox[\(f\), \(b\)]\)(\[Eta])="{ansatzchidiffquadratic}}];
lmplotl=Show[{Plot[{spindiffterm31,dfitspindiffterm31,linearfit3[[1,1]]},{\[Eta],0,0.25},PlotLegends->{"direct 3D fit","fit to residuals","per-q res. fits"}, PlotStyle->{{Orange},{Magenta,Dotted},{Blue,Dashed}}],
ErrorListPlot[Table[{linearpar3[[i,1;;2]],ErrorBar[1/Sqrt[linearpar3[[i,-1]]]]},{i,Length@linearpar3err}], PlotStyle->Red],
ErrorListPlot[Table[{linearpar3[[i,1;;2]],ErrorBar[linearpar3err[[i]]]},{i,Length@linearpar3err}], PlotLegends->{"per-q fit results"}]},
PlotRange->{{0,0.25},{Min@linearpar3[[All,2]]-0.001,Max@linearpar3[[All,2]]+0.001}}, ImageSize->330, Frame->True, FrameLabel->{Style["\[Eta]",14],Style["\!\(\*SubscriptBox[\(f\), \(a\)]\)(\[Eta])",14]}, PlotLabel->{"lin+mix fit, lin term, \!\(\*SubscriptBox[\(f\), \(a\)]\)(\[Eta])="{ansatzchidifflinear}}];
lmplotm=Show[{Plot[{spindiffterm32,dfitspindiffterm32,mixfit[[1,1]]},{\[Eta],0,0.25},PlotLegends->{"direct 3D fit","fit to residuals","per-q res. fits"}, PlotStyle->{{Orange},{Magenta,Dotted},{Blue,Dashed}}],
ErrorListPlot[Table[{mixpar[[i,1;;2]],ErrorBar[1/Sqrt[mixpar[[i,-1]]]]},{i,Length@mixparerr}], PlotStyle->Red],
ErrorListPlot[Table[{mixpar[[i,1;;2]],ErrorBar[mixparerr[[i]]]},{i,Length@mixparerr}], PlotLegends->{"per-q fit results"}]},
PlotRange->{{0,0.25},{Min@mixpar[[All,2]]-0.001,Max@mixpar[[All,2]]+0.001}}, ImageSize->330, Frame->True, FrameLabel->{Style["\[Eta]",14],Style["\!\(\*SubscriptBox[\(f\), \(c\)]\)(\[Eta])",14]}, PlotLabel->{"lin+mix fit, mix term, \!\(\*SubscriptBox[\(f\), \(c\)]\)(\[Eta])="{ansatzchidiffmix}}];
lmqplotl=Show[{Plot[{spindiffterm41,dfitspindiffterm41,linearfit4[[1,1]]},{\[Eta],0,0.25},PlotLegends->{"direct 3D fit","fit to residuals","per-q res. fits"}, PlotStyle->{{Orange},{Magenta,Dotted},{Blue,Dashed}}],
ErrorListPlot[Table[{linearpar4[[i,1;;2]],ErrorBar[1/Sqrt[linearpar4[[i,-1]]]]},{i,Length@linearpar4err}], PlotStyle->Red],
ErrorListPlot[Table[{linearpar4[[i,1;;2]],ErrorBar[linearpar4err[[i]]]},{i,Length@linearpar4err}], PlotLegends->{"per-q fit results"}]},
PlotRange->{{0,0.25},{Min@linearpar4[[All,2]]-0.001,Max@linearpar4[[All,2]]+0.001}}, ImageSize->330, Frame->True, FrameLabel->{Style["\[Eta]",14],Style["\!\(\*SubscriptBox[\(f\), \(a\)]\)(\[Eta])",14]}, PlotLabel->{"lin+mix+quad fit, lin term, \!\(\*SubscriptBox[\(f\), \(a\)]\)(\[Eta])="{ansatzchidifflinear}}];
lmqplotq=Show[{Plot[{spindiffterm43,dfitspindiffterm43,quadfit2[[1,1]]},{\[Eta],0,0.25},PlotLegends->{"direct 3D fit","fit to residuals","per-q res. fits"}, PlotStyle->{{Orange},{Magenta,Dotted},{Blue,Dashed}}],
ErrorListPlot[Table[{quadpar2[[i,1;;2]],ErrorBar[1/Sqrt[quadpar2[[i,-1]]]]},{i,Length@quadparerr2}], PlotStyle->Red],
ErrorListPlot[Table[{quadpar2[[i,1;;2]],ErrorBar[quadparerr2[[i]]]},{i,Length@quadparerr2}], PlotLegends->{"per-q fit results"}]},
PlotRange->{{0,0.25},{Min@quadpar2[[All,2]]-0.001,Max@quadpar2[[All,2]]+0.001}}, ImageSize->330, Frame->True, FrameLabel->{Style["\[Eta]",14],Style["\!\(\*SubscriptBox[\(f\), \(b\)]\)(\[Eta])",14]}, PlotLabel->{"lin+mix+quad fit, quad term, \!\(\*SubscriptBox[\(f\), \(b\)]\)(\[Eta])="{ansatzchidiffquadratic}}];
lmqplotm=Show[{Plot[{spindiffterm42,dfitspindiffterm42,mixfit2[[1,1]]},{\[Eta],0,0.25},PlotLegends->{"direct 3D fit","fit to residuals","per-q res. fits"}, PlotStyle->{{Orange},{Magenta,Dotted},{Blue,Dashed}}],
ErrorListPlot[Table[{mixpar2[[i,1;;2]],ErrorBar[1/Sqrt[mixpar2[[i,-1]]]]},{i,Length@mixparerr2}], PlotStyle->Red],
ErrorListPlot[Table[{mixpar2[[i,1;;2]],ErrorBar[mixparerr2[[i]]]},{i,Length@mixparerr2}], PlotLegends->{"per-q fit results"}]},
PlotRange->{{0,0.25},{Min@mixpar2[[All,2]]-0.001,Max@mixpar2[[All,2]]+0.001}}, ImageSize->330, Frame->True, FrameLabel->{Style["\[Eta]",14],Style["\!\(\*SubscriptBox[\(f\), \(c\)]\)(\[Eta])",14]}, PlotLabel->{"lin+mix+quad fit, mix term, \!\(\*SubscriptBox[\(f\), \(c\)]\)(\[Eta])="{ansatzchidifflinear}}];
Print[{lplot,lqplotl}]
Print[{lmplotl,lmqplotl}];
Print[{lqplotq,lmqplotq,Show[lmqplotq,PlotRange->{{0.0,0.25},Automatic}]}];
Print[{lmplotm,lmqplotm,Show[lmqplotm,PlotRange->{{0.0,0.25},Automatic}]}];
(* rewrite fit and ansatz output in terms of (\[Eta],S.\[Chi]1,\[Chi]2) and not as the expanded (\[Eta],\[Chi]1,\[Chi]2) version processed by DataFitFunction *)
uneqfit2[[1,1]] = ( ansatzFinal + ansatzchidifflinear (\[Chi]1-\[Chi]2) )/.uneqfit2[[1,2]];
uneqfit2[[1,15]] = ansatzFinal + ansatzchidifflinear (\[Chi]1-\[Chi]2);
uneqfit2[[2,1]] = ( ansatzFinal + ansatzchidifflinear (\[Chi]1-\[Chi]2) + ansatzchidiffquadratic (\[Chi]1-\[Chi]2)^2 )/.uneqfit2[[2,2]];
uneqfit2[[2,15]] = ansatzFinal + ansatzchidifflinear (\[Chi]1-\[Chi]2) + ansatzchidiffquadratic (\[Chi]1-\[Chi]2)^2;
uneqfit2[[3,1]] = ( ansatzFinal + ansatzchidifflinear (\[Chi]1-\[Chi]2) + ansatzchidiffmix S (\[Chi]1-\[Chi]2) )/.uneqfit2[[3,2]];
uneqfit2[[3,15]] = ansatzFinal + ansatzchidifflinear (\[Chi]1-\[Chi]2) + ansatzchidiffmix S (\[Chi]1-\[Chi]2);
uneqfit2[[4,1]] = ( ansatzFinal + ansatzchidifflinear (\[Chi]1-\[Chi]2) + ansatzchidiffmix S (\[Chi]1-\[Chi]2) + ansatzchidiffquadratic (\[Chi]1-\[Chi]2)^2 )/.uneqfit2[[4,2]];
uneqfit2[[4,15]] = ansatzFinal + ansatzchidifflinear (\[Chi]1-\[Chi]2) + ansatzchidiffmix S (\[Chi]1-\[Chi]2) + ansatzchidiffquadratic (\[Chi]1-\[Chi]2)^2;
Print[Style["total fit (linear in spindiff) to all data:",Blue]];
Print[uneqfit2[[1,1;;8]]];
Print[Style["total fit (linear+ quadratic in spindiff) to all data:",Blue]];
Print[uneqfit2[[2,1;;8]]];
Print[Style["total fit (linear in spindiff + spineff*spindiff) to all data:",Blue]];
Print[uneqfit2[[3,1;;8]]];
Print[Style["total fit (linear+quadratic in spindiff + spineff*spindiff) to all data:",Blue]];
Print[uneqfit2[[4,1;;8]]];
Print[Style["----------------", Blue]];
(* 2DeqS 2Dall 3Dlin 3Dlinquad 3Dlinmix 3Dlinquadmix 1D2Dparts 1Deta 1DS 2DSansatz *)
Return[{eqdatafitv2, alldataeqfitv2, uneqfit2[[1;;1]], uneqfit2[[2;;2]], uneqfit2[[3;;3]], uneqfit2[[4;;4]], extraStuff, nsfits, q1fits, spinpart}];];
]
AnsatzRestrictions[q1fit_,ansatzGen_,var_,varval_,OptionsPattern[{"Parameters"->False}]]:=Module[{ansatz,nsdegree,q1degree,expans,coeffs1,coeffs2,q1var,coeff\[Eta],expns,coeffS,sys,
coeffs,solveVars,solNS,ansatzS,solq1,ansatzFinal,parameters},
parameters=OptionValue["Parameters"];
ansatz=ansatzGen;
q1var=Variables[q1fit][[1]];
expans= Exponent[ansatz,q1var];
expns = Exponent[q1fit, q1var];
coeffs1= CoefficientList[Collect[ansatz,S],q1var];
solveVars= CoefficientList[#,var]&/@ coeffs1;
solveVars = Last/@\[NonBreakingSpace]Select[solveVars, Length@# > 0 &];
solveVars =Select[AtomsList[solveVars],!NumberQ[#1]&];
Print["solveVars: ", solveVars];
Print["q1fit: ", q1fit];coeffs1=Select[CoefficientList[Collect[ansatz,S],q1var]/. var->varval,Not@NumberQ@#&];
coeffs2 = Drop[CoefficientList[q1fit,q1var],1]; (* this won't always work *)
(*
Print["coeffs1: ", coeffs1];
Print["coeffs2: ", coeffs2];
*)
sys=Table[coeffs1[[i]]==coeffs2[[i]],{i,Min[Length@coeffs1,Length@coeffs2]}];
(*
Print["sys: ",sys];
*)
coeffs=Drop[Select[AtomsList[ansatz/.var-> varval],!NumberQ[#1]&],-1];
(*
Print["coeffs:", coeffs[[1;;Min[expans+1,expns+1]]]];
*)
solq1=Simplify@First[Solve[sys,solveVars]];
ansatzFinal=ansatz/. solq1;
Print["Applying q=1 restriction : ",solq1 ];
(*Print["Applying q=1 restriction : ",solq1//TableForm ];*)
If[parameters,solq1,ansatzFinal]
];
AnsatzRestrictions[q1fit_,ansatzGen_,OptionsPattern[{"Parameters"->False}]]:=Module[{ansatz,nsdegree,q1degree,expans,coeffs1,coeffs2,q1var,coeff\[Eta],expns,coeffS,sys,coeffs,solveVars,
solNS,ansatzS,solq1,ansatzFinal,parameters},
parameters=OptionValue["Parameters"];
ansatz=ansatzGen;
q1var=Variables[q1fit][[1]];
expans= Exponent[ansatz,q1var];
expns = Exponent[q1fit, q1var];
coeffs1= CoefficientList[Collect[ansatz,S],q1var];
solveVars= CoefficientList[#,\[Eta]]&/@ coeffs1;
solveVars = Last/@\[NonBreakingSpace]Select[solveVars, Length@# > 0 &];
solveVars =Select[AtomsList[solveVars],!NumberQ[#1]&];
Print["solveVars: ", solveVars];
Print["q1fit: ", q1fit];
coeffs1 = Select[CoefficientList[Collect[ansatz,S],q1var],Not@NumberQ@#&]/. \[Eta]->0.25;
(* coeffs2 = Take[CoefficientList[q1fit,q1var],-Length@coeffs1] (* won't always work *); *)
coeffs2 = Coefficient[q1fit,q1var^Sort@Select[Exponent[MonomialList[ansatz,q1var],q1var],#>0& ] ]; (* this should work more generally, unless maybe for non-exact S-independent constant coefficients...? *)
Print["coeffs1: ", coeffs1];
Print["coeffs2: ", coeffs2];
sys=Table[coeffs1[[i]]==coeffs2[[i]],{i,Min[Length@coeffs1,Length@coeffs2]}];
Print["sys: ",sys];
coeffs=Drop[Select[AtomsList[ansatz/. \[Eta]->0.25],!NumberQ[#1]&],-1];
(*
Print["coeffs:", coeffs[[1;;Min[expans+1,expns+1]]]];
*)
solq1=Simplify@First[Solve[sys,solveVars]];
ansatzFinal=ansatz/. solq1;
Print["Applying q=1 restriction : ",solq1//TableForm ];
If[parameters,solq1,ansatzFinal]
];
AnsatzRestrictions[nsfit_,q1fit_,ansatzGen_,OptionsPattern[{"Parameters"->False}]]:=Module[{ansatz,nsdegree,q1degree,expans,expns,nsvar,q1var,coeff\[Eta],coeffS,sys,coeffs,solNS,ansatzS,
solq1,ansatzFinal,parameters},
parameters=OptionValue["Parameters"];
ansatz=ansatzGen;
nsvar=Variables[nsfit][[1]];
q1var=Variables[q1fit][[1]];
Print["nsvar: ", nsvar];
Print["q1var: ", q1var];
expans= Exponent[ansatz,nsvar];
expns = Exponent[nsfit, nsvar];
Print["expans: ", expans];
Print["expns: ", expns];
sys=Table[CoefficientList[Collect[ansatz/. q1var->0,nsvar],nsvar][[i]]==CoefficientList[nsfit,nsvar][[i]],{i,Min[expans+1,expns+1]}];
Print["sys: ",sys];
If[Length@Select[sys,Not@TrueQ]==0,
ansatzS=ansatz,
Print[Style["Something is wrong with the restrictions",Red]]
];
Print["Applying non-spinning restriction : ", ansatz//TableForm ];
expans= Exponent[ansatz,q1var];
expns = Exponent[q1fit, q1var];
Print["expans: ", expans];
Print["expns: ", expns];
sys=Table[CoefficientList[Collect[ansatzS/. nsvar->0.25,S],q1var][[i]]==CoefficientList[q1fit,q1var][[i]],{i,Min[expans+1,expns+1]}];
Print["sys: ",sys];
coeffs=Drop[Select[AtomsList[ansatzS/. nsvar->0.25],!NumberQ[#1]&],-1];
Print["coeffs: ",coeffs];
solq1=First[Solve[sys,coeffs[[1;;Min[expans+1,expns+1]]]]];
ansatzFinal=ansatzS/. solq1;
Print["Applying q=1 restriction : ",solq1//TableForm ];
If[parameters,solq1,ansatzFinal]
];
Plot2DFits[data_?ListQ,fitlist_?ListQ,fitvars_?ListQ,OptionsPattern[{"PlotRange"->Automatic,"ShowLegend"->True,"zLabel"->"","ToolTipTags"->"","Outliers"->0.005}]]:=Module[{myfitlist,fitdata,myresidual,styles,plotrange,pos,myvar1,myvar2,
myminvar1,mymaxvar2,mymaxvar1,myminvar2,myfitlistfunc,plot1,plot2,plot3,plottab,showlegend,zlabel,tooltiptags,myresidualhist,posoutliers,outliers,mylegend},
plotrange=OptionValue["PlotRange"];
showlegend=OptionValue["ShowLegend"];
zlabel=OptionValue["zLabel"];
tooltiptags=OptionValue["ToolTipTags"];
outliers=OptionValue["Outliers"];
myfitlist=Flatten[TakeColumn[#,1]&/@fitlist,1];
myvar1=Evaluate@Symbol@ToString@fitvars[[1]];
myvar2=Evaluate@Symbol@ToString@fitvars[[2]];
If[ListQ@plotrange,myminvar1=plotrange[[1,1]];
mymaxvar1=plotrange[[1,2]];
myminvar2=plotrange[[2,1]];
mymaxvar2=plotrange[[2,2]];
,
myminvar1=Min[data[[All,1]]];
mymaxvar1=Max[data[[All,1]]];
myminvar2=Min[data[[All,2]]];
mymaxvar2=Max[data[[All,2]]];
];
myfitlist=(#/.myvar1->x/.myvar2->y)&/@myfitlist;
Print[Length[data]," data points, variables plotted -> ",{myvar1, myvar2}];
If[showlegend,
mylegend=(myfitlist/.x->myvar1/.y->myvar2);
Print[plot1=Show[{Plot3D[myfitlist,{x,myminvar1,mymaxvar1},{y,myminvar2,mymaxvar2},PlotStyle->Drop[ColorData[3,"ColorList"],1],PlotLabel->Style[ToString[zlabel],14,Black],PlotLegends->mylegend,Lighting->Automatic,PlotRange->plotrange,
AxesLabel->{Style[ToString@myvar1,14,Black],Style[ToString@myvar2,14,Black],""}],ListPointPlot3D[data,PlotStyle->PointSize[0.025]]}]],
Print[plot1=Show[{Plot3D[myfitlist,{x,myminvar1,mymaxvar1},{y,myminvar2,mymaxvar2},PlotStyle->Drop[ColorData[3,"ColorList"],1],PlotLabel->Style[ToString[zlabel],14,Black],Lighting->Automatic,PlotRange->plotrange,
AxesLabel->{Style[ToString@myvar1,14,Black],Style[ToString@myvar2,14,Black],""}],ListPointPlot3D[data,PlotStyle->PointSize[0.025]]}]]
];
plottab=Table[
Print[myfitlist[[j]]/.x->myvar1/.y->myvar2];
fitdata=Transpose[{data[[All,1]],data[[All,2]],Table[myfitlist[[j]]/.x->data[[i,1]]/.y->data[[i,2]],{i,Length@data}]}];
myresidual=Transpose[{data[[All,1]],data[[All,2]],fitdata[[All,3]]-data[[All,3]]}];
myresidualhist=myresidual[[All,3]]/fitdata[[All,3]];
Print[{plot2=Show[ListPointPlot3D[data,PlotStyle->PointSize[0.025],PlotRange->All],
Plot3D[myfitlist[[j]],{x,myminvar1,mymaxvar1},{y,myminvar2,mymaxvar2},PlotLabel->Style[ToString[zlabel],14,Black],
AxesLabel->{Style[ToString@myvar1,14,Black],Style[ToString@myvar2,14,Black],""},ImageSize->350,
PlotStyle->Drop[ColorData[3,"ColorList"],1][[j]],Lighting->Automatic,PlotRange->plotrange],AxesLabel->{Style[ToString@myvar1,14,Black],
Style[ToString@myvar2,14,Black],""}],
plot3=myListPlot3D[myresidual,PlotLabel->Style[ToString["residual"],14,Black],PlotRange->plotrange,ImageSize->350,AxesLabel->{Style[ToString@myvar1,14,Black],
Style[ToString@myvar2,14,Black],""},PlotStyle->Drop[ColorData[3,"ColorList"],1][[j]]],
If[Length@fitlist[[j,1]]==1,"No stats. available",fitlist[[j,1,3]]],Show[Histogram[myresidualhist,100],Frame->True,FrameLabel -> {"1-data/fit", "#"}],
ListPlot[Table[Tooltip[myresidualhist[[i]],If[ListQ@tooltiptags,tooltiptags[[i]],ToString[{data[[i,1]],data[[i,2]]}]]],{i,Length@myresidualhist}],PlotRange->All,Frame->True,FrameLabel -> {"#","1-data/fit"}],
Sqrt@Mean[myresidual[[All,3]]^2],
Select[(Reverse@SortBy[Transpose[{Abs@myresidualhist,If[ListQ@tooltiptags,tooltiptags,data[[All,1;;2]]]}],First]),#[[1]]>outliers&][[All,2]]//TableForm
}];
{plot2,plot3}
,{j,Length@myfitlist}];
Return[{plot1,plottab}]
]
(* ::Code::Initialization:: *)
End[];
EndPackage[];