Example of Model Reconstruction Using the Package

Author : Tomasz Konopka
Date : January 2011

The objective of this notebook is to demonstrate model reconstruction  using the package. 

The input to the analysis is a set of two time-series.  In this notebook, these are imported from two text files. 

The objective of the analysis is to determine what equations/model structures are consistent with the time-series data.

Loading data

It is nice to clear variables before starting a new computation to avoid conflicts

ClearAll["Global`*"];

Load the package. Sometimes the location of the package must be added to the path in order for Mathematica to find it

$Path=Prepend[$Path,"/home/tomasz/MyP_Bruxelles/Autoosc/v0.1"]; 
Needs["Autoosc`"];

To check if the package is loaded properly, execute the following command. If all is ok, it should display some text

AOAbout[]

Read data from disk and show it

LoadedSignals={0,0};
LoadedSignals[[1]]=ReadList["/home/tomasz/MyP_Bruxelles/Autoosc/exampledata1.dat"];LoadedSignals[[2]]=ReadList["/home/tomasz/MyP_Bruxelles/Autoosc/exampledata2.dat"];
{AOPlotModel[LoadedSignals[[1]],250,"Signal 1"],AOPlotModel[LoadedSignals[[2]],250,"Signal 2"]}

Add some noise if desired

Signals={0,0};
Signals[[1]]=AOCorruptSeries[LoadedSignals[[1]],8/100];
Signals[[2]]=AOCorruptSeries[LoadedSignals[[2]],8/100];
{AOPlotModel[Signals[[1]],250,"Signal 1"],AOPlotModel[Signals[[2]],250,"Signal 2"],AOPlotModel[{Signals[[1]],Signals[[2]]},250,"Both Signals"]}

Compute Spectra and show them

Spectra={0,0};
Spectra[[1]]=AOComputeSpec[Signals[[1]]];
Spectra[[2]]=AOComputeSpec[Signals[[2]]];
{AOPlotModelSpec[Spectra[[1]],22,250,"Spectrum 1"],AOPlotModelSpecAbs[Spectra[[1]],22,250,False,"Spectrum 1 Magnitude"],AOPlotModelSpecAbs[Spectra[[1]],22,250,True,"Spectrum 1 Magnitude (Log)"]}
{AOPlotModelSpec[Spectra[[1]],22,250,"Spectrum 2"],AOPlotModelSpecAbs[Spectra[[2]],22,250,False,"Spectrum 2 Magnitude"],AOPlotModelSpecAbs[Spectra[[2]],22,250,True,"Spectrum 2 Magnitude (Log)"]}

Preliminaries

Set a variable numcycles, setting the spacing between peaks in the spectra. In this example it is 6

numcycles=6;

Get suggestions for depths. The results here will change differ here for every time a noisy signal is generated

Do[
Print["Signal ",i," - suggested depth: ",AOSuggestDepth[Spectra[[i]],numcycles,"normal",1]];
,{i,1,2}];

The procedures only use the spectral peak heights. Here see how information in those peaks compare with the full signal.

{AOPlotModel[{AOShowAnsatz[Spectra[[1]],numcycles,4],Signals[[1]]},250,"expansion from harmonics - Signal 1"],AOPlotModel[{AOShowAnsatz[Spectra[[2]],numcycles,4],Signals[[2]]},250,"expansion from harmonics - Signal 2"]}

Based on the above, decide how many harmonics to use in the procedure below

numharmonics=4;

Set the type of models to analyze. Here chose all models with maximal order equal to 1

AOBuildEqList[1];
Print["Number of model structures considered will be ",Length[AOGetEqList[]]];

Set tolerance levels. The arrays can be of different lengths from the specified depth

AOGetLevels[]
GoDepth=2;
ToleranceArray={1,1,1};
RIToleranceArray={0.5,1,2};
parameterranges ={{pp1,-6000,6000},{pp2,0,300},{qq1,-6000,6000},{qq2,0,300}};
AOSetLevels[GoDepth,parameterranges,ToleranceArray,RIToleranceArray]
AOGetLevels[]

Model reconstruction

Using the frequency domain

Compute parameters for one model structure.

try1=AOComputeParameters[Spectra[[1]],Spectra[[2]],numcycles,numharmonics,{{3,0,1},{6,1,0}},2,"S",False]

The result shows the equation tested, the equation code, and the computed parameter values. 
Interpret the result, i.e. apply the parameter and tolerance criteria to suggest whether the model is correct or not

AOInterpretParameters[try1[[3]]]

The result shows whether the parameter values are consistent within the defined tolerance levels. If yes, one set of parameter values is also given.
Alternatively, compute parameters for all model structures. This can take some time.

parametersall={0,0};
parametersall[[1]]=AOComputeParametersAllModels[Spectra[[1]],Spectra[[2]],numcycles,numharmonics,2,"A",False];
parametersall[[2]]=AOComputeParametersAllModels[Spectra[[2]],Spectra[[1]],numcycles,numharmonics,2,"A",False];

Since the results from parameter calculations are large, it is usually only interesting to view them one line at a time

parametersall[[1,3]]

Interpret all the model structures in a row

interpretations={0,0};
interpretations[[1]]=AOInterpretParametersList[parametersall[[1]]];
interpretations[[2]]=AOInterpretParametersList[parametersall[[2]]];

The results after  interpreations are often small, so it is sometimes interesting to list them all at once.
Get the number of remaining models

Print["Number of remaining models for signal 2 driven by signal 1: ",Length[interpretations[[1]]]];
Print["Number of remaining models for signal 1 driven by signal 2: ",Length[interpretations[[2]]]];

Using the time domain

The time domain can be used to check that the suggested model equations do really produce signals similar to those input

Can show plots of reconstructed curves for all of the model structures shortlisted in the frequency domain 
Here, the plots show signal 2 being driven by signal 1. The items before the plots show the equation code and the equation itself.

Do[
temp=interpretations[[1,i]];
temp2=AORecreateSignal[Spectra[[1]],Signals[[2,1]],numcycles,numharmonics,temp[[3]],temp[[2]],0];
Print[temp[[3]],"   ",temp[[1]],"   ",AOPlotModel[{temp2,Signals[[2]]},250,""]];
,{i,1,Length[interpretations[[1]]]}];

Such reconstructions can be used for model selection. First compute all the errors

interpretationsWerr={0,0};
interpretationsWerr[[1]]=AOComputeTimeDomainErrList[Spectra[[1]],Signals[[2]],numcycles,numharmonics,interpretations[[1]],0];
interpretationsWerr[[2]]=AOComputeTimeDomainErrList[Spectra[[2]],Signals[[1]],numcycles,numharmonics,interpretations[[2]],0];

Then apply interpretation criteria. Here the threshold is set as a multiple of the signal offset

timeinterpretations={0,0};
timeinterpretations[[1]]=AOInterpretTimeDomainErrList[interpretationsWerr[[1]],0.04*Re[Spectra[[2,1]]]];
timeinterpretations[[2]]=AOInterpretTimeDomainErrList[interpretationsWerr[[2]],0.04*Re[Spectra[[1,1]]]];

Again print the numbers of remaining models

Print["After time domain analysis, number of remaining models for signal 2 driven by signal 1: ",Length[timeinterpretations[[1]]]];
Print["After time domain analysis, number of remaining models for signal 1 driven by signal 2: ",Length[timeinterpretations[[2]]]];

Finally, can look at the reduced set of models and their time domain reconstructions
Here, this is done for the two cases, 1 driving 2, and 2 driving 1

Print["Models where signal 1 drives signal 2"];
Do[
temp=timeinterpretations[[1,i]];
temp2=AORecreateSignal[Spectra[[1]],Signals[[2,1]],numcycles,numharmonics,temp[[3]],temp[[2]],0];
Print[temp[[3]],"   ",temp[[1]],"   ",AOPlotModel[{temp2,Signals[[2]]},250,""]];
,{i,1,Length[timeinterpretations[[1]]]}];
Print["\nModels where signal 2 drives signal 1"];
Do[
temp=timeinterpretations[[2,i]];
temp2=AORecreateSignal[Spectra[[2]],Signals[[1,1]],numcycles,numharmonics,temp[[3]],temp[[2]],0];
Print[temp[[3]],"  ",temp[[1]],"   ",AOPlotModel[{temp2,Signals[[1]]},250,""]];
,{i,1,Length[timeinterpretations[[2]]]}];

Conclusion

That's all! 

The number of model equations listed at the end is usually small. They represent educated guesses for what mechanisms might be responsible for generating the data signals.

The equations originally used to generate the curves in the imported signals were actually the Lotka - Volterra equations. The correct equation code for both the relation between signal 1 driving signal 2 and for signal 2 driving signal 1 is therefore {{2,1,0},{4,1,1}}. Are these codes among the short list? 

When noise levels are low and tolerance levels are set at reasonable values, the Lotka-Volterra equations should appear on the short list. They may even be the only solutions left after all the model selection steps. 

When noise levels are high, when tolerance levels are set too strictly, or both, the short list might be empty or show model structures that do not include the Lotka-Volterra equations.