(* Package for automated analysis of oscillating signals. Its main goal is to determine equation structures that may generate given oscillatory signals *)
(* Author: Tomasz Konopka *)
(* Date: v1 September 2010 - January 2011 *)

BeginPackage["Autoosc`"]

(* Definitions of functions available in the package. For usage notes, see separate package documentation *)
AOAbout::usage="AOAbout[] - to give some information about the package";

AOBuildEqList::usage="AOBuildEqList[maxorder] - to recreate a list of model structures";

AOComputeParameters::usage="AOComputeParameters[ser1, ser2, numcyc, numharm, modeldef, depth, numer, toquiet] - to compute one or more sets of values for parameters appearing in a model equation";

AOComputeParametersAllModels::usage="AOComputeParametersAllModels[ser1, ser2, numcyc, numharm, depth, numer, toquiet] - to batch compute the parameters associated with many model structures";

AOComputeTimeDomainErrList::usage=  "AOComputeTimeDomainErrList[ser1, ser2, numcyc, numharm, modellist, tstep] - to compute the error between a reconstructed signal and its original for a list of proposed model structures";

AOCorruptSeries::usage="AOCorruptSeries[ser, level] - to add normal noise to a time series";

AOCorruptSpec::usage="AOComputeSpec[ser1] - to add normal noise to a spectrum";

AOComputeSpec::usage="AOCorruptSpec[ser, level] - to compute the spectrum of a time series. Note: the function essentially uses the built-in Fourier[· · · ]
functionality in Mathematica. However, it also rescales some values outputed by that command. The
other functions in the package assume all spectra have been computed with this function
";

AODelaySpec::usage="AODelaySpec[ser, tau] - to produce a spectrum of a time-delayed signal";

AOGetEqList::usage="AOGetEqList[] - to check what types of model structures are currently set to be tested";

AOGetLevels::usage="AOGetLevels[] - to verify what criteria are currently set for parameter interpretation";

AOGetRecreateSignalJump::usage="AOGetRecreateSignalJump[] - to retrieve the current value of the jump parameter used in signal reconstruction";

AOInterpretParameters::usage="AOInterpretParameters[prevres] - to apply the parameter interpretation criteria determined with AOSetLevels to sets of parameter
values";

AOInterpretParametersList::usage="AOInterpretParametersList[prevres] - to interpret parameters for many models at once";

AOInterpretTimeDomainErrList::usage="AOInterpretTimeDomainErrList[prevres, thresh] - to interpret and invalidate models based on the time domain reconstruction criteria";

AOModelSubSample::usage="AOModelSubSample[ser, pos1, pos2, ssz] - to produce a short time series from a long time series";

AOPlotModel::usage="AOPlotModel[ser, is, title] - to generate a plot of a time series with some standardized settings";

AOPlotModelSpec::uage="AOPlotModelSpec[spec, speclen, is, title] - to produce a plot of a signal spectrum";

AOPlotModelSpecAbs::usage="AOPlotModelSpecAbs[spec, speclen, is, tolog, title] - to produce a plot of the aboslute value of a signal spectrum";

AORecreateSignal::usage="AORecreateSignal[ser1, ic, numcyc, numharm, modeldef, params, tstep] - to build a signal φ using an initial condition, model equation, and a driving signal";

AOSetLevels::usage="AOSetLevels[dep, ranges, tole, reimtole] - to set options and tolerance levels used in parameter interpretation";

AOSetEqList::usage="AOSetEqList[newlist] - to custom set a list of model structures to be tested";

AOSetRecreateSignalJump::usage="AOSetRecreateSignalJump[jump] - to set the jump parameter used in signal reconstruction";

AOShowAnsatz::usage="AOShowAnsatz[ser1, numcyc, numharm] - to construct a time series of the ansatz used in parameter reconstruction";

AOShowTerms::usage="AOShowTerms[pp1not, pp2not, ordnot, signalX, signalXD] - to list the types of terms and their codes";

AOSuggestDepth::usage="AOSuggestDepth[ser1, numcyc, noisetype, thr] - to estimate how many parameter sets should be computed from the signal";


(* Deprecated *)
AOCleanSpec::usage="AOCleanSpec[ser, interv, noisetype, thr] - deprecated";

(* This is a list of parameters used in the package. The statements define them as global variables *)
pp1
pp2
qq1
qq2
t
X[t]
XD[t]

Begin["`Private`"]

(* Some private variables for this module. *)
(* These are for plotting *)
iss=8;
m1=Graphics[{Purple,Thick,Circle[{0,0},1]},ImageSize->iss];
m2=Graphics[{Purple,Thick,Line[{{-0.5,-0.5},{0.5,-0.5},{0.5,0.5},{-0.5,0.5},{-0.5,-0.5}}]},ImageSize->iss];
m3=Graphics[{Gray,Circle[{0,0},1]},ImageSize->iss];
m4=Graphics[{Gray,Line[{{-0.5,-0.5},{0.5,-0.5},{0.5,0.5},{-0.5,0.5},{-0.5,-0.5}}]},ImageSize->iss];

(* These are for parameter interpretation *)
rangevalues=Table[i*10^6,{j,1,4},{i,-1,1,2}];
rangevalues[[2,1]]=0;
rangevalues[[4,1]]=0;
depthlevel=3;
tolelevels=Table[1,{i,1,depthlevel}];
reimtolelevels=Table[1,{i,1,depthlevel}];

(* These are for the parameter computations *)
ansatz1
ansatz2
ExpToZ
z
EqList
EqListLen

(* This is for speeding up the time domain signal reconstruction *)
recreatesignaljump=10;

(* display some basic info about the package - name, version, author *)
AOAbout[]:=Module[{},
Print["Package Autoosc\nVersion: 0.1\nAuthor: Tomasz Konopka"];
Return[];
];


(* Takes a long series and returns a subsampled one *)
AOModelSubSample[ser_,pos1_,pos2_,ssz_]:=Module[{newser,newlen},
newlen=Floor[(pos2-pos1)*ssz]-1;
newser=Table[ser[[pos1+Floor[t/ssz]]],{t,0,newlen}];
Return[newser];
];

(* Functions for plotting waves and spectra. They return standard listplot objects with predefined characteristics *)
AOPlotModelSpec[spec_,speclen_,is_,title_]:=Module[{ans,plotspecRe,plotspecIm},
plotspecRe=Table[{i-1,Re[spec[[i]]]},{i,1,speclen}];
plotspecIm=Table[{i-1,Im[spec[[i]]]},{i,1,speclen}];
ans=ListPlot[{plotspecRe,plotspecIm},ImageSize->is,PlotMarkers->{m1,m2,m3,m4},PlotStyle->{Purple,Black},PlotRange->All,Frame->True,AxesOrigin->{0,0},PlotLabel->title];
Return[ans];
];
AOPlotModelSpecAbs[spec_,speclen_,is_,tolog_,title_]:=Module[{ans,plotspec,tt1,tt2,pmm},
If[Length[spec]==2,
pmm={m1,m3};
plotspec=Table[0,{i,1,2}];
tt1=Table[{i-1,If[tolog,Log[Abs[spec[[1,i]]]],Abs[spec[[1,i]]]]},{i,1,speclen}];
plotspec[[1]]=tt1;
tt1=Table[{i-1,If[tolog,Log[Abs[spec[[2,i]]]],Abs[spec[[2,i]]]]},{i,1,speclen}];
plotspec[[2]]=tt1;
,
pmm=m1;
plotspec=Table[{i-1,If[tolog,Log[Abs[spec[[i]]]],Abs[spec[[i]]]]},{i,1,speclen}];
];

ans=ListPlot[plotspec,ImageSize->is,PlotRange->All,PlotMarkers->pmm,Frame->True,AxesOrigin->{0,0},PlotLabel->title];
Return[ans];
];
AOPlotModel[ser_,is_,title_]:=Module[{ans,ps},
If[Length[ser]==2,ps={Purple,Darker[Gray]},ps=Purple];
ans=ListLinePlot[ser,ImageSize->is,PlotStyle->ps,PlotRange->All,Frame->True,Axes->None,PlotLabel->title];
Return[ans];
];

(* Adds noise to the signals in the time and frequency domains *)
AOCorruptSeries[ser_,nl_]:=Module[{newser,len,i,mean},
len=Length[ser];
mean=0;
Do[mean+=ser[[i]],{i,1,len}];
mean=mean/len;
newser=Table[Max[0,RandomReal[NormalDistribution[ser[[i]],Abs[nl*mean]]]],{i,1,len}];
Return[newser];
];
AOCorruptSpec[ser_,nl_]:=Module[{newserre,newserim,len,i},
len=Length[ser];
newserre=Table[N[Re[ser[[i]]]],{i,1,len}];
newserim=Table[N[Im[ser[[i]]]],{i,1,len}];
Do[
If[newserre[[i]]≠0,
newserre[[i]]=RandomReal[NormalDistribution[newserre[[i]],Abs[nl*newserre[[i]]]]];
];
If[newserim[[i]]≠0,
newserim[[i]]=RandomReal[NormalDistribution[newserim[[i]],Abs[nl*newserim[[i]]]]];
];
,{i,1,len}];
Return[newserre+I*newserim];
];

(* Sets some parts of a spectrum to zero if lower than a threshold about a noise *)
AOCleanSpec[ser_,interv_,noisetype_,thr_]:=Module[{ans,len,drawres,i,nlRe,nlIm,co},
ans=ser;
drawres={};
len=Floor[Length[ser]/2];

(* will store average of the real and imaginary parts of the non-fourier-peak components *)
nlRe=0;
nlIm=0;
(* number of non-fourier-peak components *)
co=0;

(* Do different actions depending on the noise type in the data. For the moment, only normal noise is implemented. *)
If[noisetype=="normal",
Do[
If[Mod[i-1,interv]≠0,
nlRe+=Abs[Re[ser[[i]]]];
nlIm+=Abs[Im[ser[[i]]]];
co++;
]
,{i,2,len}];
nlRe=nlRe/co;
nlIm=nlIm/co;
Do[
If[Mod[i-1,interv]==0 &&Abs[Re[ser[[i]]]]<thr*nlRe && Abs[Re[ser[[i]]]]<thr*nlIm,
ans[[i]]=0;
];
,{i,2,len}];

];

Return[ans];
];

(* This introduces a delay tau to the spectrum ser. Returns another spectrum *)
AODelaySpec[ser_,tau_]:=Module[{ans,len,i,omeg,rotre,rotim,oldmu,oldnu},
ans=ser;
len=Length[ser];
omeg=2*Pi/len;
Do[
rotre=Re[Exp[(i-1)*I*omeg*tau]];
rotim=Im[Exp[(i-1)*I*omeg*tau]];
oldmu=Re[ser[[i]]];
oldnu=Im[ser[[i]]];
ans[[i]]=(rotre*oldmu-rotim*oldnu)+I*(rotre*oldnu+rotim*oldmu);
,{i,1,len}];
Return[ans];
];


(* This returns a list of numbers made from the ansatz *)
AOShowAnsatz[ser1_,numcyc_,numharm_]:=Module[{ans,len,i},

(* Use the BuildAnsatz function to create ansatz1 and ansatz2 from the spectrum in ser1 *)
(* This is somewhat inefficient but avoids rewriting the BuildAnsatz code *)
AOBuildAnsatz[ser1,ser1,numcyc,numharm];

(* now use ansatz1 to create a time series *)
len=Length[ser1];
ans=Table[ReplaceAll[ansatz1,t->i],{i,0,len-1}];

Return[ans];
];

(* Computes the spectrum of a series. It is readjusted as described in the package documentation *)
AOComputeSpec[timeser_]:=Module[{ans},
ans=2*Fourier[timeser,FourierParameters->{-1,1}];
ans[[1]]/=2;
Return[ans];
];

(* looks at the real and imaginary parts of a number. Returns boolean depending on which part is greater. Uses a multiplicative threshold *)
MoreImThanRe[inp_,tthr_]:=Module[{},
If[tthr*Abs[Re[inp]]>Abs[Im[inp]],Return[False];,Return[True];];
Return[True];
];

(* This uses the spectrum to build two ansatz functions. The two ansatz are stored in variables private to the module *)
AOBuildAnsatz[ser1_,ser2_,numcyc_,numharm_]:=Module[{kk,omeg,Xacore,XC,XS},
ansatz1=0;
ansatz2=0;

(* define the base frequency of the signal *)
omeg=2*Pi*numcyc/Length[ser1];

(* Create the oscillating part to z substitution step *)
ExpToZ={};
Do[
ExpToZ=Append[ExpToZ,Exp[I*i*omeg*t]->z^i]
,{i,-Floor[20*numharm],Floor[20*numharm]}];

(* define a generic ansatz function with general coefficients *)
Xacore= TrigToExp[Sum[XC[n]*Cos[n*omeg*t]+1*XS[n]*Sin[n*omeg*t],{n,0,numharm}]];
ansatz1 = Xacore;
ansatz2 = Xacore;
(* use the spectra to substitute values into the generic ansatz *)
Do[
ansatz1=ReplaceAll[ansatz1,{XC[kk]->Re[ser1[[1+kk*numcyc]]],XS[kk]->Im[ser1[[1+kk*numcyc]]]}];
ansatz2=ReplaceAll[ansatz2,{XC[kk]->Re[ser2[[1+kk*numcyc]]],XS[kk]->Im[ser2[[1+kk*numcyc]]]}];
,{kk,0,numharm}];

Return[];
];



(* This build one ansatz from spectrum and leaves it in trig form *)
AOBuildTrigAnsatz[ser1_,numcyc_,numharm_]:=Module[{ans,kk,omeg,XC,XS,n},
ans=0;

(* define the base frequency of the signal *)
omeg=2*Pi*numcyc/Length[ser1];

(* define a generic ansatz function with general coefficients *)
ans= Sum[XC[n]*Cos[n*omeg*(t-1)]+1*XS[n]*Sin[n*omeg*(t-1)],{n,0,numharm}];

(* use the spectra to substitute values into the generic ansatz *)
Do[
ans=ReplaceAll[ans,{XC[kk]->Re[ser1[[1+kk*numcyc]]],XS[kk]->Im[ser1[[1+kk*numcyc]]]}];
,{kk,0,numharm}];

Return[ans];
];


(* This takes a model code and returns a list of two equations. One is an equation showing fractions, for human use. The other is one with the fractions multiplied out. *)
AOBuildModelEq[modeldef_]:=Module[{modeq,type1,type2,ord1X,ord1XD,ord2X,ord2XD,ans},
(* ans will contain two things. 1 is nicely formatted model equation. 2 is model format used in calcualtions (no fractions). *)
ans={0,0};

type1=modeldef[[1,1]];
ord1X=Floor[modeldef[[1,2]]];
ord1XD=Floor[modeldef[[1,3]]];
type2=modeldef[[2,1]];
ord2X=Floor[modeldef[[2,2]]];
ord2XD=Floor[modeldef[[2,3]]];

modeq=-D[X[t],t];
Which[type1==1,modeq+=pp1,
type1==2,modeq+=pp1*(X[t]^ord1X),
type1==3,modeq+=pp1*(XD[t]^ord1XD),
type1==4,modeq+=pp1*(X[t]^ord1X)*(XD[t]^ord1XD),
type1==5 ,modeq+=pp1/(pp2+(X[t]^ord1X)),
type1==6 ,modeq+=pp1*(X[t]^ord1X)/(pp2+(X[t]^ord1X)),
type1==7 ,modeq+=pp1*(XD[t]^ord1XD)/(pp2+(X[t]^ord1X)),
type1==8 ,modeq+=pp1*(X[t]^ord1X)*(XD[t]^ord1XD)/(pp2+(X[t]^ord1X)),
type1==9 ,modeq+=pp1/(pp2+(XD[t]^ord1XD)),
type1==10 ,modeq+=pp1*(X[t]^ord1X)/(pp2+(XD[t]^ord1XD)),
type1==11 ,modeq+=pp1*(XD[t]^ord1XD)/(pp2+(XD[t]^ord1XD)),
type1==12,modeq+=pp1*(X[t]^ord1X)*(XD[t]^ord1XD)/(pp2+(XD[t]^ord1XD))
];
Which[type2==1,modeq+=qq1,
type2==2,modeq+=qq1*(X[t]^ord2X),
type2==3,modeq+=qq1*(XD[t]^ord2XD),
type2==4,modeq+=qq1*(X[t]^ord2X)*(XD[t]^ord2XD),
type2==5 ,modeq+=qq1/(qq2+(X[t]^ord2X)),
type2==6 ,modeq+=qq1*(X[t]^ord2X)/(qq2+(X[t]^ord2X)),
type2==7 ,modeq+=qq1*(XD[t]^ord2XD)/(qq2+(X[t]^ord2X)),
type2==8 ,modeq+=qq1*(X[t]^ord2X)*(XD[t]^ord2XD)/(qq2+(X[t]^ord2X)),
type2==9 ,modeq+=qq1/(qq2+(XD[t]^ord2XD)),
type2==10 ,modeq+=qq1*(X[t]^ord2X)/(qq2+(XD[t]^ord2XD)),
type2==11 ,modeq+=qq1*(XD[t]^ord2XD)/(qq2+(XD[t]^ord2XD)),
type2==12,modeq+=qq1*(X[t]^ord2X)*(XD[t]^ord2XD)/(qq2+(XD[t]^ord2XD))
];
ans[[1]]=modeq;

(*multiply out the fractions*)
If[type1≥5 && type1≤8,
modeq *= (pp2+(X[t]^ord1X));
modeq=Simplify[modeq];
];
If[type1≥9 && type1≤12,
modeq *= (pp2+(XD[t]^ord1XD));
modeq=Simplify[modeq];
];
If[type2≥5 && type2≤8,
modeq *= (qq2+(X[t]^ord2X));
modeq=Expand[Simplify[modeq]];
];
If[type2≥9 && type2≤12,
modeq *= (qq2+(XD[t]^ord2XD));
modeq=Simplify[modeq];
];
ans[[2]]=Simplify[Expand[Simplify[modeq]]];

Return[ans];
];


(* Given one model code and its model equation, this computes sets of parameter values. Number of sets is determine dy depth *)
AOComputeAllSols[modeldef_,modeq_,depth_,numer_,toquiet_]:=Module[
{allsols,i,znow,lam0,lampos,lamneg,sss,type1,type2,npara1,npara2,ntype,opttype,ptosolve,eqstosolve},

allsols={};

(* make a list of the unknowns in the model *)
ptosolve={};
ptosolve=Append[ptosolve,pp1];
ptosolve=Append[ptosolve,qq1];

type1=modeldef[[1,1]];
type2=modeldef[[2,1]];
If[type1≥5,npara1=2,npara1=1];
If[type2≥5,npara2=2,npara2=1];
If[npara1==2,ptosolve=Append[ptosolve,pp2];];
If[npara2==2,ptosolve=Append[ptosolve,qq2];];

(* Determine whether computations will be done with Solve or NSolve function *)
If[numer=="S" || numer=="N",
ntype=numer;
,
If[type1>4,ntype="N",ntype="S"];
];

(* precompute all the coefficients that will be needed *)
lam0=Expand[Coefficient[modeq,z,0]];
lampos=Table[Expand[Coefficient[modeq,z,znow]],{znow,1,depth+4}];
lamneg=Table[Expand[Coefficient[modeq,z,-znow]],{znow,1,depth+4}];

(* Loop, for each depth value, compute a set of parameters *)
Do[
eqstosolve={};

(* set up the equations to be solved  *)
Which[npara1==1 && npara2==1,eqstosolve={lam0==0,lampos[[znow]]==0};,
npara1==1 && npara2==2 && (type2==5 || type2==9),eqstosolve={lampos[[znow]]==0,lamneg[[znow]]==0,lam0==0};,
npara1==1 && npara2==2 && (type2!=5 && type2!=9),eqstosolve={lampos[[znow]]==0,lamneg[[znow]]==0,lampos[[znow+1]]==0};,
npara1==2 && npara2==1,(* will not happen *) eqstosolve={lampos[[znow]]==0,lamneg[[znow]]==0,lampos[[znow+1]]==0};,
npara1==2 && npara2==2 && (type1==5 || type1==9 || type2==5 || type2==9),eqstosolve={lam0==0,lampos[[znow]]==0,lamneg[[znow]]==0,lampos[[znow+1]]==0};,
npara1==2 && npara2==2 && (type1!=5 && type1!=9 && type2!=5 && type2!=9),eqstosolve={lampos[[znow]]==0,lamneg[[znow]]==0,lampos[[znow+1]]==0,lamneg[[znow+1]]==0};
];

(* use conventional or numerical solver *)
If[ntype=="N",
If[toquiet,
sss=Quiet[NSolve[N[eqstosolve],ptosolve]];
,
sss=NSolve[N[eqstosolve],ptosolve];
];
,
If[toquiet,
sss=Quiet[Solve[eqstosolve,ptosolve]];
,
sss=Solve[eqstosolve,ptosolve];
];
];

allsols=Append[allsols,{znow,sss}];
,{znow,1,depth}];

Return[allsols];
];


(* This prepares the output of the parameter estimation computation. It assumes the ansatz for the primary and driving are already built *)
AOComputeParametersWReadyAnsatz[numcyc_,numharm_,modeldef_,depth_,numer_,toquiet_]:=Module[
{ans,npara1,npara2,type1,type2, ord1X, ord1XD, ord2X, ord2XD,modeq,allsols},

ans={};

type1=modeldef[[1,1]];
ord1X=Floor[modeldef[[1,2]]];
ord1XD=Floor[modeldef[[1,3]]];
type2=modeldef[[2,1]];
ord2X=Floor[modeldef[[2,2]]];
ord2XD=Floor[modeldef[[2,3]]];

(*determine number of parameters *)
If[type1≥5,npara1=2,npara1=1];
If[type2≥5,npara2=2,npara2=1];

(* Build the model equation *)
modeq=AOBuildModelEq[modeldef];
ans=Append[ans,modeq[[1]]];
ans=Append[ans,modeldef];
modeq=modeq[[2]];

(* Exit if types are not strictly ordered *)
If[type1≥type2, Return[Append[ans,{}]]];
(* Reject if model does not involve the driving function at all *)
If[(type1==1 || type1==2 || type1==5)&& (type2==2 || type2==5 || type2==6), Return[Append[ans,{}]]];
(* reject the case where pp1 and qq1 cannot be independently determined *)
If[type1==5 && type2==9,Return[Append[ans,{}]]];

(* sub into the equation and convert to z *)
modeq=ReplaceAll[modeq,{XD[t]->ansatz1,D[XD[t],t]->D[ansatz1,t],X[t]->ansatz2,D[X[t],t]->D[ansatz2,t]}];
modeq=Expand[ReplaceAll[modeq,ExpToZ]];

(* solve for the parameters at each depth *)
allsols=AOComputeAllSols[modeldef,modeq,depth,numer,toquiet];
ans=Append[ans,allsols];

Return[ans];
];

(* This is the procedure that is called from outside to compute parameters given a model code *)
AOComputeParameters[ser1_,ser2_,numcyc_,numharm_,modeldef_,depth_,numer_,toquiet_]:=Module[{},
(* Build the ansatz *)
AOBuildAnsatz[ser1,ser2,numcyc,numharm];
(* Compute parameters *)
Return[AOComputeParametersWReadyAnsatz[numcyc,numharm,modeldef,depth,numer,toquiet]];
];

(* Here, the user can change setting for the model selection/interpretation steps *)
AOSetLevels[depth_,ranges_,tole_,reimtole_]:=Module[{quickcheck,temp,kk},
quickcheck=True;

If[depth<1,Print["Depth level must be greater or equal than 1"];quickcheck=False;];

If[Length[tole]<depth,Print["Tolerance list is too short. Must be of length ",depth]; quickcheck=False;];
If[Length[reimtole]<depth,Print["ReIm tolerance list is too short. Must be of length ",depth]; quickcheck=False;];

If[quickcheck==False,Return[]];

depthlevel=depth;
tolelevels=tole;
reimtolelevels=reimtole;

(* Extract range values from the input ranges *)
(* give each parameter an initial value *)
kk=Length[ranges];
Do[
If[Length[ranges[[i]]]==3,
temp=-1;
If[ranges[[i,1]]==pp1,temp=1;];
If[ranges[[i,1]]==pp2,temp=2;];
If[ranges[[i,1]]==qq1,temp=3;];
If[ranges[[i,1]]==qq2,temp=4;];
If[temp>0,
rangevalues[[temp,1]]=ranges[[i,2]]; rangevalues[[temp,2]]=ranges[[i,3]];
];
,
Print["Cannot interpret ranges: ",ranges[[i]]];
];
,{i,1,kk}];

Return[];
];


(* Display the currently set parameter interpretation criteria *)
AOGetLevels[]:=Module[{},
Print["depth: ",depthlevel,"\nTolerances: ",tolelevels,"\nReImTolerances: ",reimtolelevels,"\nParameter ranges: ",rangevalues];
Return[];
];


(* This takes as input the results of ComputeParameters and some other options. It then checks whether the parameters obey the conditions *)
AOInterpretParameters[prevres_]:=Module[
{ans,possiblesols,tempok,temptempok,depth,alltole,allsolslen,solok,verdict,allsols,kk,pp1sol,pp2sol,qq1sol,qq2sol,i,j,znow,quickcheck,temp},

ans={False,{}};
possiblesols={};

(* Now look through the solutions and see if the results are consistent with each other *)
allsols=prevres;
depth=Length[allsols];
depth=Min[depth,depthlevel];
If[depth==0,Return[ans];];

(* Here loop through all the solutions and subsolutions *)
allsolslen=Table[Length[allsols[[znow,2]]],{znow,1,depth}];
kk=allsolslen[[1]];
solok=Table[True,{i,1,kk}];
alltole=Table[True,{i,1,kk},{znow,1,depth},{j,1,Max[allsolslen]}];
Do[
pp1sol=pp1/.allsols[[1,2,i]];
pp2sol=pp2/.allsols[[1,2,i]];
qq1sol=qq1/.allsols[[1,2,i]];
qq2sol=qq2/.allsols[[1,2,i]];

(* Reject if pp1sol is more imaginary than real *)
If[pp1sol==pp1,(* means not evaluated pp1 *),(* never happens*),
(* here pp1 sol has a value *)
If[MoreImThanRe[pp1sol,reimtolelevels[[1]]] ||Re[pp1sol]>rangevalues[[1,2]] || Re[pp1sol]<rangevalues[[1,1]],solok[[i]]=False;];
];
If[pp2sol==pp2,(* means not evaluated pp2 *),(* never happens*),
(* here pp2 sol has a value *)
If[MoreImThanRe[pp2sol,reimtolelevels[[1]]]||Re[pp2sol]>rangevalues[[2,2]] || Re[pp2sol]<rangevalues[[2,1]],solok[[i]]=False;];
];
If[qq1sol==qq1,(* means not evaluated qq1 *),(* never happens*),
(* here qq1 sol has a value *)
If[MoreImThanRe[qq1sol,reimtolelevels[[1]]] ||Re[qq1sol]>rangevalues[[3,2]] || Re[qq1sol]<rangevalues[[3,1]],solok[[i]]=False;];
];
If[qq2sol==qq2,(* means not evaluated qq2 *),(* never happens*),
(* here qq2 sol has a value *)
If[MoreImThanRe[qq2sol,reimtolelevels[[1]]]||Re[qq2sol]>rangevalues[[4,2]] || Re[qq2sol]<rangevalues[[4,1]],solok[[i]]=False;];
];

(* proceed with checks at further depth only if the first item is acceptable *)
If[solok[[i]],
Do[
Do[
(*check current solution against first order using tolerances *)
temp=pp1/.allsols[[znow,2,j]];
If[temp==pp1,(* *),(* *),
If[Abs[(Re[pp1sol]-Re[temp])/Re[pp1sol]]>tolelevels[[znow]]  || MoreImThanRe[temp,reimtolelevels[[znow]]],alltole[[i,znow,j]]=False;];
];
temp=pp2/.allsols[[znow,2,j]];
If[temp==pp2,(* *),(* *),
If[Abs[(Re[pp2sol]-Re[temp])/Re[pp2sol]]>tolelevels[[znow]]|| MoreImThanRe[temp,reimtolelevels[[znow]]] ,alltole[[i,znow,j]]=False;];
];
temp=qq1/.allsols[[znow,2,j]];
If[temp==qq1,(* *),(* *),
If[Abs[(Re[qq1sol]-Re[temp])/Re[qq1sol]]>tolelevels[[znow]] || MoreImThanRe[temp,reimtolelevels[[znow]]],alltole[[i,znow,j]]=False;];
];
temp=qq2/.allsols[[znow,2,j]];
If[temp==qq2,(* *),(* *),
If[Abs[(Re[qq2sol]-Re[temp])/Re[qq2sol]]>tolelevels[[znow]]|| MoreImThanRe[temp,reimtolelevels[[znow]]] ,alltole[[i,znow,j]]=False;];
];

,{j,1,allsolslen[[znow]]}];
,{znow,2,depthlevel}];
];

,{i,1,kk}];


(* search the alltole table and make a final verdict *)
Do[
tempok=True;
Do[
temptempok=False;
Do[
temptempok=temptempok || alltole[[i,znow,j]];
,{j,1,allsolslen[[znow]]}];
tempok=tempok && temptempok;
,{znow,2,depth}];
solok[[i]]=solok[[i]] &&tempok;
,{i,1,kk}];

verdict=False;
Do[
If[solok[[i]], possiblesols=Append[possiblesols,allsols[[1,2,i]]]; verdict=True;];
,{i,1,kk}];

ans[[1]]=verdict;
ans[[2]]=possiblesols;
Return[ans];

];

(* These allow the user to change the quality of the signal recreation steps *)
AOSetRecreateSignalJump[jump_]:=Module[{tj},
tj=Floor[Re[jump]];
(* Avoid cases with jump less than 1 *)
If[tj<1, recreatesignaljump=1;,recreatesignaljump=tj;];
Return[];
];
AOGetRecreateSignalJump[]:=Module[{},
Return[recreatesignaljump];
];


(* This recreates a signal, uses a manual ODE integration *)
(* The function also uses a speedup mechanism when recreatesignaljump is more than one. The speedup comes from assuming that the driving signal can be approximated by a straight line on small intervals. This avoids computing many evaluations of sine and cosine functions *)
AORecreateSignalManual[ser1_,ic_,numcyc_,numharm_,modeldef_,params_,tstep_]:=Module[
{ans,DF,len,kk,tt,omeg,toff,pp1sol,pp2sol,qq1sol,qq2sol,ok,driveansatz,DFprev,DFstart,DFend,ansadd,ansprev,
jumplen,jumpsections,jumpleftover,xc,xs,type1,type2,ord1X,ord2X,ord1XD,ord2XD},

type1=modeldef[[1,1]];
ord1X=Floor[modeldef[[1,2]]];
ord1XD=Floor[modeldef[[1,3]]];
type2=modeldef[[2,1]];
ord2X=Floor[modeldef[[2,2]]];
ord2XD=Floor[modeldef[[2,3]]];

(* Parameters are assumed to be ok. They were checked in the RecreateSignal part *)
pp1sol=pp1/.params[[1]];
pp2sol=pp2/.params[[1]];
qq1sol=qq1/.params[[1]];
qq2sol=qq2/.params[[1]];

pp1sol=Re[pp1sol];
pp2sol=Re[pp2sol];
qq1sol=Re[qq1sol];
qq2sol=Re[qq2sol];

len=Floor[(Length[ser1]+1)/tstep];
(* jumplen is the time interval over which the linear approximation to the driving signal will be made *)
jumplen=Min[Floor[1/tstep],recreatesignaljump];
jumpsections=Floor[len/jumplen];
jumpleftover=len-jumplen*jumpsections;

(* use the DF variable to store reconstructed driving signal *)
toff=Floor[1/tstep];
omeg=2*Pi*numcyc/len;
DF=Table[0,{t,1,len}];

If[jumplen==0,
(* Compute the driving signal the normal way*)

Do[
xc=Re[ser1[[1+kk*numcyc]]];
xs=Im[ser1[[1+kk*numcyc]]];
Do[
DF[[tt]]+=xc*Cos[kk*omeg*(tt-1)]+xs*Sin[kk*omeg*(tt-1)];
,{tt,1,len}];

,{kk,0,numharm}];


,
(* For other jumplen, compute the driving signal with linear approximations *)

(* Compute the value of the first point *)
DFend=0;
Do[
xc=Re[ser1[[1+kk*numcyc]]];
xs=Im[ser1[[1+kk*numcyc]]];
DFend+=xc*Cos[kk*omeg*0]+xs*Sin[kk*omeg*0];
,{kk,0,numharm}];

(* Loop over the sections of the jump *)
Do[
(* Obtain boundary values for DF within the section *)
(* DFstart is the same as DFend on previous round. The new DFend should be recomputed *)
DFstart=DFend;
DFend=0;
tt=0+js*jumplen;
Do[
xc=Re[ser1[[1+kk*numcyc]]];
xs=Im[ser1[[1+kk*numcyc]]];
DFend+=xc*Cos[kk*omeg*tt]+xs*Sin[kk*omeg*tt];
,{kk,0,numharm}];

(* loop over the time instances within the sections and apply the linear transformation *)
Do[
DF[[jt+(js-1)*jumplen]]=DFstart+ ((jt-1)*(DFend-DFstart)/jumplen);
,{jt,1,jumplen}];

,{js,1,jumpsections}];

(* loop over the time instances leftover at the end *)
DFstart=DFend;
DFend=0;
tt=len-1;
Do[
xc=Re[ser1[[1+kk*numcyc]]];
xs=Im[ser1[[1+kk*numcyc]]];
DFend+=xc*Cos[kk*omeg*tt]+xs*Sin[kk*omeg*tt];
,{kk,0,numharm}];

Do[
DF[[jt+jumpsections*jumplen]]=DFstart+ ((jt-1)*(DFend-DFstart)/jumpleftover);
,{jt,1,jumpleftover}]

];

(* now reconstruct the time series using the initial condition and the reconstructed driving signal *)
ans=Table[0,{i,1,len}]; 
ans[[1]]=ic; (* Sets the initial condition *)
Do[
ans[[t]]=ans[[t-1]];
ansprev=ans[[t-1]];
DFprev=DF[[t-1]];
ansadd=0;
Which[type1==1 ,ansadd+=pp1sol,
type1==2,ansadd+=pp1sol*(ansprev^ord1X),
type1==3,ansadd+=pp1sol*(DFprev^ord1XD),
type1==4,ansadd+=pp1sol*(ansprev^ord1X)*(DFprev^ord1XD),
type1==5 ,ansadd+=pp1sol/(pp2sol+ansprev^ord1X),
type1==6 ,ansadd+=(pp1sol*(ansprev^ord1X))/(pp2sol+(ansprev^ord1X)),
type1==7 ,ansadd+=(pp1sol*(DFprev^ord1XD))/(pp2sol+(ansprev^ord1X)),
type1==8 ,ansadd+=(pp1sol*(ansprev^ord1X)*(DFprev^ord1XD))/(pp2sol+(ansprev^ord1X)),
type1==9 ,ansadd+=pp1sol/(pp2sol+(DFprev^ord1XD)),
type1==10 ,ansadd+=(pp1sol*(ansprev^ord1X))/(pp2sol+(DFprev^ord1XD)),
type1==11 ,ansadd+=(pp1sol*DFprev)/(pp2sol+(DFprev^ord1XD)),
type1==12,ansadd+=(pp1sol*(ansprev^ord1X)*(DFprev^ord1XD))/(pp2sol+(DFprev^ord1XD))
];
Which[type2==1 ,ansadd+=qq1sol,
type2==2,ansadd+=qq1sol*(ansprev^ord2X),
type2==3,ansadd+=qq1sol*(DFprev^ord2XD),
type2==4,ansadd+=qq1sol*(ansprev^ord2X)*(DFprev^ord2XD),
type2==5 ,ansadd+=qq1sol/(qq2sol+(ansprev^ord2X)),
type2==6 ,ansadd+=(qq1sol*(ansprev^ord2X))/(qq2sol+(ansprev^ord2X)),
type2==7 ,ansadd+=(qq1sol*(DFprev^ord2XD))/(qq2sol+(ansprev^ord2X)),
type2==8 ,ansadd+=(qq1sol*(ansprev^ord2X)*(DFprev^ord2XD))/(qq2sol+(ansprev^ord2X)),
type2==9 ,ansadd+=qq1sol/(qq2sol+(DFprev^ord2XD)),
type2==10 ,ansadd+=(qq1sol*(ansprev^ord2X))/(qq2sol+(DFprev^ord2XD)),
type2==11 ,ansadd+=(qq1sol*(DFprev^ord2XD))/(qq2sol+(DFprev^ord2XD)),
type2==12,ansadd+=(qq1sol*(ansprev^ord2X)*(DFprev^ord2XD))/(qq2sol+(DFprev^ord2XD))
];
ans[[t]]=ansprev+tstep*ansadd;

,{t,2,len}];

(* Finally, subsample the signal and return it *)
ans=AOModelSubSample[ans,1,Length[ans],tstep];

Return[ans];
];


(* This takes a model equation, parameters, driving signal, and an initial condition. Recreates the primary signal curve *)
(* When tstep is 0, the function uses NDSolve. Otherwise, it calls another function for manual integration of the ODE *)
AORecreateSignal[ser1_,ic_,numcyc_,numharm_,modeldef_,params_,tstep_]:=Module[
{ans,len,kk,omeg,type1,type2,
ord1X,ord1XD,ord2X,ord2XD,pp1sol,pp2sol,qq1sol,qq2sol,ok,thiseq,Xsol,driveansatz},

type1=modeldef[[1,1]];
ord1X=Floor[modeldef[[1,2]]];
ord1XD=Floor[modeldef[[1,3]]];
type2=modeldef[[2,1]];
ord2X=Floor[modeldef[[2,2]]];
ord2XD=Floor[modeldef[[2,3]]];

pp1sol=pp1/.params[[1]];
pp2sol=pp2/.params[[1]];
qq1sol=qq1/.params[[1]];
qq2sol=qq2/.params[[1]];

ok=True;
(* check if model type is consistent with parameters *)
If[pp1sol==pp1,ok=False];
If[qq1sol==qq1,ok=False];
If[type1≥5 && pp2sol==pp2,ok=False];
If[type2≥5 && qq2sol==qq2,ok=False];
If[ok==False, Print["Submitted parameters do not match the model type"]; Return[];];

(* Check if needs to be computed using NDSolve or using the manual procedure *)
If[tstep>0,Return[AORecreateSignalManual[ser1,ic,numcyc,numharm,modeldef,params,tstep]]];

(* Continuing here iwth the NDSolve procedure *)
pp1sol=Re[pp1sol];
pp2sol=Re[pp2sol];
qq1sol=Re[qq1sol];
qq2sol=Re[qq2sol];

len = Length[ser1];

thiseq=AOBuildModelEq[modeldef];
thiseq=thiseq[[1]];

(* Build the ansatz for the driving signal (in trig form, without any 'i' terms) *)
driveansatz=AOBuildTrigAnsatz[ser1,numcyc,numharm];

thiseq=ReplaceAll[thiseq,XD[t]->driveansatz];
thiseq=ReplaceAll[thiseq,{pp1->pp1sol,pp2->pp2sol,qq1->qq1sol,qq2->qq2sol}];

(* Solve the differential equation. MaxSteps is set arbitrarily to a large number *)
Xsol=NDSolve[{thiseq==0,X[1]==ic},X,{t,1,len},MaxSteps->len*10^6];

(* Create a table for the signal from the Xsol solution *)
ans=Table[(X[kk]/.Xsol)[[1]],{kk,1,len}];

Return[ans];
];


(* Takes in a list of models for given primary and driving signals. Recreates the primary signal for each model in the list and computes the standard error between the primary and recreated signals. Returns a new list of models with the standard error appended for each model *)
AOComputeTimeDomainErrList[ser1_,ser2_,numcyc_,numharm_,modellist_,tstep_]:=Module[
{len,ans,i,tempser,tempmodel,lse,tslen},

len=Length[modellist];
tslen=Length[ser1];
ans={};

(* look at each model in the model list*)
Do[
tempmodel=modellist[[i]];
tempser=AORecreateSignal[ser1,ser2[[1]],numcyc,numharm,tempmodel[[3]],tempmodel[[2]],tstep];
lse=0;
Do[
lse+=(tempser[[i]]-ser2[[i]])^2;
,{i,1,tslen}];
lse=Sqrt[lse]/tslen;
tempmodel=Append[tempmodel,lse];
ans=Append[ans,tempmodel];

,{i,1,len}];

Return[ans];
];


(* Takes in a list like that outputed from the AOComputeTimeDomainErrList function, and a threshold value. Returns a shorter model list with only those models that have an error smaller than the specified threshold. *)
AOInterpretTimeDomainErrList[prevres_,thresh_]:=Module[{len,i,ans},
len=Length[prevres];
ans={};
Do[
If[prevres[[i,5]]<thresh,
ans=Append[ans,prevres[[i]]];
];
,{i,1,len}];

Return[ans];
];

(* This just returns a list of terms so that users can know how to construct the modeldef lists  *)
AOShowTerms[pp1not_,pp2not_,ordnot_,signalX_,signalXD_]:=Module[{ans,i,j},
ans=Table[i,{i,1,12},{j,1,2}];
ans[[1,2]]=pp1not;
ans[[2,2]]=pp1not*signalX;
ans[[3,2]]=pp1not*signalXD;
ans[[4,2]]=pp1not*signalX*signalXD;
ans[[5,2]]=pp1not/(pp2not+signalX^ordnot);
ans[[6,2]]=pp1not*signalX/(pp2not+signalX^ordnot);
ans[[7,2]]=pp1not*signalXD/(pp2not+signalX^ordnot);
ans[[8,2]]=pp1not*signalX*signalXD/(pp2not+signalX^ordnot);
ans[[9,2]]=pp1not/(pp2not+signalXD^ordnot);
ans[[10,2]]=pp1not*signalX/(pp2not+signalXD^ordnot);
ans[[11,2]]=pp1not*signalXD/(pp2not+signalXD^ordnot);
ans[[12,2]]=pp1not*signalX*signalXD/(pp2not+signalXD^ordnot);
Return[ans];
];

(* This does the parameter computation for all models specified in the equation list *)
AOComputeParametersAllModels[ser1_,ser2_,numcyc_,numharm_,depth_,numer_,toquiet_]:=Module[
{autoparams,ans,ordertype,m1,m2,m1min,m1max,m2min,m2max,aaaa,order1,order2},

ans={};
m1max=12;
m2max=12;

(* Build the ansatz once for all the equation structures *)
AOBuildAnsatz[ser1,ser2,numcyc,numharm];

(* Compute parameters for all the models *)
Do[
Do[aaaa=AOComputeParametersWReadyAnsatz[numcyc,numharm,EqList[[i]],depth,numer,toquiet]];
Do[ans=Append[ans,aaaa]];
,{i,1,EqListLen}];

Return[ans];
];


(* This one takes a list of computed parameters. It outputs the consistent sets only for those models that seem to be consistent *)
AOInterpretParametersList[prevres_]:=Module[{ans,i,tlen,temp,tok},
tlen=Length[prevres];
ans={};

Do[
temp=prevres[[i]];
tok = AOInterpretParameters[temp[[3]]];
If[tok[[1]]==True,
ans=Append[ans,{temp[[1]],tok[[2]],temp[[2]],temp[[3]]}];
];
,{i,1,tlen}];

Return[ans];
];


(* This looks at a spectrum, and suggests a depth *)
AOSuggestDepth[ser1_,numcyc_,noisetype_,thr_]:=Module[{ans,nlRe,nlIm,nlAbs,co,len,i},
ans=0;
len=Floor[Length[ser1]/2];

(* Suggested depth should depend on the type of noise present in the spectrum. For the moment, only normal noise is supported *)
If[noisetype=="normal",

(* compute average real and imaginary parts of non-fourier-peak components *)
nlRe=0; nlIm=0;nlAbs=0;
co=0;

Do[
If[Mod[i-1,numcyc]≠0,
nlRe+=Abs[Re[ser1[[i]]]];
nlIm+=Abs[Im[ser1[[i]]]];
nlAbs+=Abs[ser1[[i]]];
co++;
];
,{i,2,len}];

nlRe=nlRe/co;
nlIm=nlIm/co;
nlAbs=nlAbs/co;

(* now go through the fourier frequencies and find lowest number of peaks above noise *)
i=1+numcyc;
ans=0;
While[i<len &&Abs[Re[ser1[[i]]]]>(1+thr)*nlRe && Abs[Im[ser1[[i]]]]>(1+thr)*nlIm  && Abs[ser1[[i]]]>(1+thr)*nlAbs,
ans++;
i+=numcyc;
];

];

(* return computed computed level minus one because for many models depth=x calculations involved depth=x+1 coefficients *)
Return[Max[0,ans-1]];
];

(* These define the list of equation to test *)
AOBuildEqList[maxorder_]:=Module[{termlist,i,j,TList,ord1max,ord2max,ter,TListLen},
TList={{1,0,0}};
ord1max=1;
ord2max=2;

(* First build a list of terms with the maximal order specified *)
Do[
If[ter==1||ter==3||ter==9||ter==11,ord1max=0,ord1max=maxorder];
If[ter==1||ter==2||ter==5||ter==6,ord2max=0,ord2max=maxorder];

If[ord1max==0||ord2max==0,
(* if one of the ordmax is zero, one field is missing in the term, just loop through the other orders*)
If[ord1max==0&&ord2max≠0,
Do[
TList=Append[TList,{ter,ord1max,ord2}];
,{ord2,1,ord2max}];
];
If[ord2max==0&&ord1max≠0,
Do[TList=Append[TList,{ter,ord1,ord2max}];
,{ord1,1,ord1max}];
];
,
(* both fields are present and need to be included in the loop through orders *)
Do[
Do[TList=Append[TList,{ter,ord1,ord2}];
,{ord1,1,ord1max}];
,{ord2,1,ord2max}];
];
,{ter,1,12}];
TListLen=Length[TList];

(* Now use the term list to define a list of all equations *)
EqList={};
Do[
Do[
If[TList[[i,3]]==0&&TList[[j,3]]==0,
(*Skip this case because it doesn't involve a driving field*)
,
EqList=Append[EqList,{TList[[i]],TList[[j]]}]];,{j,i+1,TListLen}];
,{i,1,TListLen}];
EqListLen=Length[EqList];

];


(* Here the user can specify a list of models to test *)
AOSetEqList[newlist_]:=Module[{},
EqList=newlist;
EqListLen=Length[EqList];
];


(* Returns the list of models currently set to be tested *)
AOGetEqList[]:=Module[{},
Return[EqList];
];


(* When the module is loaded, iby default t automatically builds a list of equations. Here of order 1 *)
AOBuildEqList[1]

End[]
EndPackage[]