(* ::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.                                                         *)
(************************************************************************)



BeginPackage[ "HG1G2tool`"]


formBase::usage = "Initialize the basis functions for HG1G2 and HG12 systems. Call once before using other HG1G2tools."

fitHG1G2::usage = "Fit HG1G2-system for given data. Use:
	{\!\(\*
StyleBox[\"fit\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"covmat\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\",\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"rms\",\nFontSlant->\"Italic\"]\)} = fitHG1G2[\!\(\*
StyleBox[\"data\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"weighted\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"opts\",\nFontSlant->\"Italic\"]\)], where inputs are
\!\(\*
StyleBox[\"data\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"-\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)List of {angle,V(angle)}-pairs, i.e. the observations of the phase function (in magnitude scale).
\!\(\*
StyleBox[\"weighted\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"-\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)Optional. Either one common value for the standard error of the observations, or list of standard errors for each observation (in magnitude scale). Default value is 0.03 mag.
Options:
\"\!\(\*
StyleBox[\"Degrees\",\nFontSlant->\"Italic\"]\)\" - Set to True (default) if angles are given in degrees and to False otherwise.
Output:
\!\(\*
StyleBox[\"fit\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"-\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)Best fit linear parameters, converted to H,G1,G2-values. Note: use simulateParameters-function for more reliable estimate of the best fit H,G1,G2-parameters, togehter with error estimate.
\!\(\*
StyleBox[\"covmat\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"-\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)Covariance matrix of the parameters.
\!\(\*
StyleBox[\"rms\",\nFontSlant->\"Italic\"]\) - root-mean-square value of the fit, in magnitude scale."

fitHG12::usage = "Fit HG12-system for given data. Use:
	{\!\(\*
StyleBox[\"fit\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"covmat\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\",\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"rms\",\nFontSlant->\"Italic\"]\)} = fitHG12[\!\(\*
StyleBox[\"data\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"weighted\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"opts\",\nFontSlant->\"Italic\"]\)], where inputs are
\!\(\*
StyleBox[\"data\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"-\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)List of {angle,V(angle)}-pairs, i.e. the observations of the phase function (in magnitude scale).
\!\(\*
StyleBox[\"weighted\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"-\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)Optional. Either one common value for the standard error of the observations, or list of standard errors for each observation (in magnitude scale). Default value is 0.03 mag.
Options:
\"\!\(\*
StyleBox[\"Degrees\",\nFontSlant->\"Italic\"]\)\" - Set to True (default) if angles are given in degrees and to False otherwise.
Output:
\!\(\*
StyleBox[\"fit\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"-\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)Best fit linear parameters, converted to H,G12-values. Note: use simulateParameters-function for more reliable estimate of the best fit H,G12-parameters, togehter with error estimate.
\!\(\*
StyleBox[\"covmat\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\"-\",\nFontSlant->\"Italic\"]\)\!\(\*
StyleBox[\" \",\nFontSlant->\"Italic\"]\)Covariance matrix of the parameters.
\!\(\*
StyleBox[\"rms\",\nFontSlant->\"Italic\"]\) - root-mean-square value of the fit, in magnitude scale."

simulateErrors::usage = "Simulates the best estimates and errors for HG1G2 or HG12 system parameters. Use: \!\(\*
StyleBox[\"estimates\",\nFontSlant->\"Italic\"]\) = simulateErrors[ \!\(\*
StyleBox[\"fit\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"covmat\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"simN\",\nFontSlant->\"Italic\"]\)], where \!\(\*
StyleBox[\"fit\",\nFontSlant->\"Italic\"]\) and \!\(\*
StyleBox[\"covmat\",\nFontSlant->\"Italic\"]\) are given by the fitHG1G2 or fitHG12-functions. The third parameter, \!\(\*
StyleBox[\"simN\",\nFontSlant->\"Italic\"]\), is optional and sets the number of simulated values in the estimation. By default the value is 100 000. Result \!\(\*
StyleBox[\"estimates\",\nFontSlant->\"Italic\"]\) is a table with four columns for \!\(\*
StyleBox[\"H\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"G1\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"G2\",\nFontSlant->\"Italic\"]\) and \!\(\*
StyleBox[\"k\",\nFontSlant->\"Italic\"]\) (i.e. slope) parameters, and with 8 rows. The two first rows are the estimated mean and median of the parameter. We recommend using median as best fit parameter, but mean can also be used. Next 6 rows give the estimate error limits for -3\[Sigma],-2\[Sigma],-\[Sigma],\[Sigma],2\[Sigma],3\[Sigma] (i.e. confidence interval limits 99.73% lower, 95.45% lower, 68.3% lower, 68.3% upper, 95.45% upper and 99.73% upper interval limit)."

systemVersion::usage = "Gives the version number of the basis functions of the system in use. Use: systemVersion[\!\(\*
StyleBox[\"opts\",\nFontSlant->\"Italic\"]\)]. If called with option \"Range\"->True (default is False), will also report the angle range (in radians) where the basis function system is defined."

giveBaseValues::usage = "Mainly for debuggin purposes, you can check the base functions. Use: {baseVal1, baseVal2, baseVal3} = giveBaseValues[ \!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"opts\",\nFontSlant->\"Italic\"]\)], where \!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\) is value of angle. Function will return the value of the three base functions for the given angle. Option \"Degrees\" can be given. If True (default), value \!\(\*
StyleBox[\"x\",\nFontSlant->\"Italic\"]\) is assumed to be in degrees."


HG1G2tools::notImplemented = "Initialization of base functions from file is not currently implemented."

HG1G2tools::outOfRange = "Value `1` is `2` the range `3` of the basis function."

HG1G2tools::notInited = "Base functions are not initialized, call formBase[] first."


Begin["`Private`"]


defaultKnots1 = {{7.5 \[Degree],0.75`},{30.0 \[Degree],0.33486016`},{60.0 \[Degree],0.13410560`},{90.0 \[Degree],0.051104756`},{120.0 \[Degree],0.021465687`},{150.0 \[Degree],0.0036396989`}};
defaultDerivs1 = {-1.909859317102744`,-0.091328612`};
defaultKnots2 = {{7.5 \[Degree],0.925`},{30.0 \[Degree],0.62884169`},{60.0 \[Degree],0.31755495`},{90.0 \[Degree],0.12716367`},{120.0 \[Degree],0.022373903`},{150.0 \[Degree],0.00016505689`}};
defaultDerivs2 = {-0.5729577951308232`,-8.6573138` 10^-8};
defaultKnots3 = {{0.0 \[Degree],1},{0.3 \[Degree],0.83381185`},{1.0 \[Degree],0.57735424`},{2.0 \[Degree],0.42144772`},{4.0 \[Degree],0.23174230`},{8.0 \[Degree],0.10348178`},{12.0 \[Degree],0.061733473`},{20.0 \[Degree],0.016107006`},{30.0 \[Degree],0}};
defaultDerivs3 = {-0.10630097,0};
defaultSystemVersion = 20101000;


coefG1small=0.7527; constG1small=0.06164;
coefG1large=0.9529; constG1large=0.02162;
coefG2small=-0.9612; constG2small=0.6270;
coefG2large=-0.6125; constG2large=0.5572;


funcTypeConstant = 1;
funcTypeLinear = 2;
funcTypeSpline = 3;


base1 = {};
base2 = {};
base3 = {};
systemVersionNum = 0;


Options[formBase] = {}

formBase[ ] := Module[{baseSpline1, baseSpline2, baseSpline3},
(* Init spline system *)
baseSpline1 = formSpline[defaultKnots1, defaultDerivs1];
baseSpline2 = formSpline[defaultKnots2, defaultDerivs2];
baseSpline3 = formSpline[defaultKnots3, defaultDerivs3];
(* init 1st base function *)
base1 = {{0.0, 150.0 \[Degree]},
{funcTypeLinear, 0.0 \[Degree], 7.5 \[Degree], 1.0, -1.90985931710274},
{funcTypeSpline, 7.5 \[Degree], 150.0 \[Degree], baseSpline1}
};
(* init 2nd base function *)
base2 = {{0.0, 150.0 \[Degree]},
{funcTypeLinear, 0.0 \[Degree], 7.5 \[Degree], 1.0,-0.572957795130823},
{funcTypeSpline, 7.5 \[Degree], 150.0 \[Degree], baseSpline2}
};
(* init erd base function *)
base3 = {{0.0, 150.0 \[Degree]},
{funcTypeSpline, 0.0\[Degree], 30.0 \[Degree], baseSpline3},
{funcTypeConstant, 30.0 \[Degree], 150.0 \[Degree], 0.0}
};
systemVersionNum = defaultSystemVersion
]

formBase[ fileName_] := Module[ {},
Message[ HG1G2tools::notImplemented]
]


Options[ fitHG1G2] = {"Degrees"->True}

fitHG1G2[data_, weighted_:0.03, OptionsPattern[]] := Module[ {Amat, xvec, yvec, wvec, covmat, lm, fit, rms},
If[ systemVersionNum == 0, Message[HG1G2tools::notInited];Return[]];
xvec = data[[All,1]];
If[OptionValue[ "Degrees"]===True, xvec *= \[Degree]];
yvec = 10^(-0.4 data[[All,2]]);
Amat = Map[ {getValue[ base1, #], getValue[ base2, #], getValue[ base3, #]} &, xvec];
wvec = 1/(yvec (10^(0.4 weighted)-1.0))^2;
lm = LinearModelFit[ {Amat,yvec}, Weights->wvec, VarianceEstimatorFunction->(1 &)];
fit = linearToHsystem[ lm["BestFitParameters"]];
covmat = lm["CovarianceMatrix"];
rms = Sqrt[6.25 Total[ Log[10,lm["Response"]/lm["PredictedResponse"]]^2] / (Length[data]-3)];
{fit, covmat,rms}
]


Options[ fitHG12] = {"Degrees"->True}

fitHG12[data_, weighted_:0.03, OptionsPattern[]] := Module[ {Bmat, AmatSmall, AmatLarge, xvec, yvec, wvec, covmatSmall, covmatLarge, lmSmall, lmLarge, fitSmall, fitLarge, rmsSmall, rmsLarge},
If[ systemVersionNum == 0, Message[HG1G2tools::notInited];Return[]];
xvec = data[[All,1]];
If[OptionValue[ "Degrees"]===True, xvec *= \[Degree]];
yvec = 10^(-0.4 data[[All,2]]);
Bmat = Map[ {getValue[ base1, #], getValue[ base2, #], getValue[ base3, #]} &, xvec];
AmatSmall = Map[ {constG1small #[[1]]+constG2small #[[2]]+(1-constG1small-constG2small) #[[3]], coefG1small #[[1]]+coefG2small #[[2]]-(coefG1small+coefG2small) #[[3]]} &, Bmat];
AmatLarge = Map[ {constG1large #[[1]]+constG2large #[[2]]+(1-constG1large-constG2large) #[[3]], coefG1large #[[1]]+coefG2large #[[2]]-(coefG1large+coefG2large) #[[3]]} &, Bmat];
wvec = 1/(yvec (10^(0.4 weighted)-1.0))^2;
lmSmall = LinearModelFit[ {AmatSmall,yvec}, Weights->wvec, VarianceEstimatorFunction->(1 &)];
lmLarge = LinearModelFit[ {AmatLarge,yvec}, Weights->wvec, VarianceEstimatorFunction->(1 &)];
fitSmall = linearToHsystem[ lmSmall["BestFitParameters"]];
covmatSmall = lmSmall["CovarianceMatrix"];
fitLarge = linearToHsystem[ lmLarge["BestFitParameters"]];
covmatLarge = lmLarge["CovarianceMatrix"];
rmsSmall = Sqrt[6.25 Total[ Log[10,lmSmall["Response"]/lmSmall["PredictedResponse"]]^2] / (Length[data]-2)];
rmsLarge = Sqrt[6.25 Total[ Log[10,lmLarge["Response"]/lmLarge["PredictedResponse"]]^2] / (Length[data]-2)];
If[ rmsSmall < rmsLarge,
{fitSmall, covmatSmall, rmsSmall},
{fitLarge, covmatLarge, rmsLarge}
]
]


simulateErrors[ fit_, covmat_, simN_:100000] := Module[ {linFit, simLin, simHsystem, simk, Cmat, U, \[CapitalLambda], is, perc = {0.5,0.0013499,0.0227501,0.158655,0.841345,0.97725,0.99865}},
{\[CapitalLambda], U} = Eigensystem[ covmat];
Cmat =  U\[Transpose] . DiagonalMatrix[ Sqrt[ Chop[\[CapitalLambda]]]];
is = RandomReal[ NormalDistribution[], {simN, Length[fit]}];
linFit = HsystemToLinear[ fit];
simLin = (linFit + Cmat . is\[Transpose])\[Transpose];
simHsystem = Map[ linearToHsystem[#, 3] &, simLin];
simk = Map[ -(30.0 #[[2]]+9.0 #[[3]])/(5.0 \[Pi] (#[[2]]+#[[3]])) &, simHsystem];
{Join[ {Mean[simHsystem[[All,1]]]},Quantile[ simHsystem[[All,1]],perc]], Join[ {Mean[simHsystem[[All,2]]]},Quantile[ simHsystem[[All,2]],perc]], Join[ {Mean[simHsystem[[All,3]]]},Quantile[ simHsystem[[All,3]],perc]], Join[ {Mean[simk]},Quantile[ simk,perc]]}\[Transpose]
]


Options[ systemVersion] = {"Range"->False}

systemVersion[OptionsPattern[]] := Module[ {},
If[ systemVersionNum == 0, Message[HG1G2tools::notInited];Return[]];
If[ OptionValue[ "Range"]===True,
{systemVersionNum, Max[{base1[[1,1]],base2[[1,1]],base3[[1,1]]}], Min[{base1[[1,2]],base2[[1,2]],base3[[1,2]]}]},
systemVersionNum
]
]


Options[ giveBaseValues] = {"Degrees"->True}

giveBaseValues[ x_, OptionsPattern[]] := Module[ {xx},
If[ systemVersion == 0, Message[HG1G2tools::notInited];Return[]];
xx = If[ OptionValue["Degrees"]===True, x \[Degree], x];
{getValue[ base1, xx], getValue[ base2, xx], getValue[ base3, xx]}
]


linearToHsystem[ {a1_,a2_,a3_}] := linearToHsystem[ {a1,a2,a3}, 3]
linearToHsystem[ {a1_,a2_}] := linearToHsystem[ {a1,a2}, 2]

linearToHsystem[ {a1_,a2_,a3_}, 3] := {-2.5 Log[10,a1+a2+a3], a1/(a1+a2+a3), a2/(a1+a2+a3)}

linearToHsystem[ {a1_,a2_}, 2] := {-2.5 Log[10,a1], a2/a1}

linearToHsystem[ {a1_,a2_}, 3] := Module[ {H,G1,G2},
H = -2.5*Log[10,a1];
If[a2<0.2,
G1=coefG1small a2+constG1small;
G2=coefG2small a2+constG2small,
G1=coefG1large a2+constG1large;
G2=coefG2large a2+constG2large
];
{H,G1,G2}
]


HsystemToLinear[ {H_,G1_,G2_}] := HsystemToLinear[ {H,G1,G2}, 3]
HsystemToLinear[ {H_,G12_}] := HsystemToLinear[ {H,G12}, 2]

HsystemToLinear[ {H_,G1_,G2_}, 3] := {Exp[-0.9210340371976184 H] G1, Exp[-0.9210340371976184 H] G2, Exp[-0.9210340371976184 H] (1-G1-G2)}

HsystemToLinear[ {H_,G12_}, 2] := {Exp[-0.9210340371976184 H], Exp[-0.9210340371976184 H] G12}


getValue[ basef_, x_] := Module[{funci = 2},
If[ x < basef[[1,1]],
Message[ HG1G2tools::outOfRange, x, "below", basef[[1,1]]];
Return[]
];
If[ x > basef[[1,2]],
Message[ HG1G2tools::outOfRange, x, "above", basef[[1,2]]];
Return[]
];
While[ !(basef[[funci,2]] <= x <= basef[[funci,3]]), funci++];
Switch[ basef[[funci,1]],
funcTypeConstant, basef[[funci,4]],
funcTypeLinear, basef[[funci,4]] + x basef[[funci,5]],
funcTypeSpline, splineValue[ x, basef[[funci,4]]]
]
]


splineValue[ x_, splineCoefs_] := Module[ {xd,yd,a,b,t, knoti = 1},
knoti= If[ x == splineCoefs[[-1,1]], Length[splineCoefs]-1,
Position[ splineCoefs[[All,1]], xx_ /; x < xx, 1, 1][[1,1]]-1
];
xd = splineCoefs[[knoti+1,1]]-splineCoefs[[knoti,1]];
yd = splineCoefs[[knoti+1,2]]-splineCoefs[[knoti,2]];
a = splineCoefs[[knoti,3]] xd - yd;
b = -splineCoefs[[knoti+1,3]] xd + yd;
t = (x-splineCoefs[[knoti,1]])/xd;
(1-t) splineCoefs[[knoti,2]] + t splineCoefs[[knoti+1,2]] + t (1-t) (a (1-t) + b t)
]


formSpline[ knots_, derivs_] := Module[ {Atri, A, B, n, x, y, Ds},
n = Length[knots];
{x,y} = knots\[Transpose];
Atri = Table[ {1/(x[[j]]-x[[j-1]]), 2/(x[[j]]-x[[j-1]])+2/(x[[j+1]]-x[[j]]), 1/(x[[j+1]]-x[[j]])}, {j, 2, n-1}]\[Transpose];
A = DiagonalMatrix[ Join[ {1}, Atri[[2]], {1}]] + 
DiagonalMatrix[ Join[ Atri[[1]], {0}],-1] + 
DiagonalMatrix[ Join[ {0}, Atri[[3]]],1];
B = Join[ {derivs[[1]]}, Table[ 3 ((y[[j]]-y[[j-1]])/(x[[j]]-x[[j-1]])^2+(y[[j+1]]-y[[j]])/(x[[j+1]]-x[[j]])^2), {j, 2, n-1}], {derivs[[2]]}];
Ds = LinearSolve[ A, B];
{x, y, Ds}\[Transpose]
]


End[]


EndPackage[]
