* {*********************************************************************

File: opt/inv.cm

  This file shows a very simple example of solution of an inverse problem
with the shell. We have an experiment in which we measure a set of data and
we have a model with which we describe this experiment mathematically. Some
parameters of this model are unknown and we would like to estimate them.
  Imagine a particle which moves in one dimension. We for example know that
we can describe its motion by x(t)=a+b*t+c*sin(d*t), but we don't know exactly
which are the values of the coefficients a, b, c and d. We simply measure
something (in our case we decide to measre positions of the particle at
specific times) and then try to fit the unknown coefficients so that the
deviation between the measurements and the appropriate quantities calculated
with the model will be as small as possible.
  We don't have real measurements for this example. Instead, we will generate
imaginary measurements by adding stohastic errors to the values which would
be exact for our model with pre-assumed values of model parameters.

Author: Igor Gresovnik, April 1998

************************************************************************}

setfile{outfile inv.ct}

*{ We define an expression evaluator's function which describes the assumed
model. It will be useful later on: }

${position[t,a,b,c,d]:a+b*t+c*sin(d*t)}

*{ We generate an imaginary set of measurements. First we assume some values
of model parameters: }

={a:4} ={b:1} ={c:4} ={d:0.5}

write{"$$$ before definition of times"\n}

*{ We define vector of times at which the positions will be measured... }

setvector{times 20}
={i:1}
while{ (i<=getvector["times",0])
[
  setvector
  { times 
    { ${i}:${0.5*(i-1)} }
  }
  ={i:i+1}
]}
printvector{times}
fprintvector{times}
write{\n} fwrite{\n}

*{ ... and then form a vector of the exact measurements; dimension is the same
as for vector "times": }

setvector{measexact ${getvector["times",0]} } 

*{ Components are calculated according to our model: }

={i:1}
while{ (i<=getvector["measexact",0])
[
  setvector
  { measexact
    { ${i}:${position[getvector["times",i],a,b,c,d]} }
  }
  ={i:i+1}
]}
printvector{measexact}
fprintvector{measexact}
write{\n} fwrite{\n}

*{ Next we define a vector of standard deviations of measurements "sigma": }

={error:0.05}
setvector{sigma ${getvector["times",0]} } 
={i:1}
while{ (i<=getvector["sigma",0])
[
  setvector
  { sigma 
    { ${i}:$ error }
  }
  ={i:i+1}
]}
printvector{sigma}
fprintvector{sigma}
write{\n} fwrite{\n}

simmeas{}
printvector{meas}
fprintvector{meas}
write{\n} fwrite{\n}

setoption{chi2spec} *{ Automatic evaluation of chi square }
setvector{measmom ${getvector["times",0]}}
setvector{parammom 4}

*{ Analysis block. This block  is interpreted each time chi square is
evaluated. It contains the actual definition of chi square: }

analysis
{
  *{ Within the analysis block we must perform sumulation and calculate
  simulated measurements. Normally, a program for simulation of the experiment
  would be run and its results would be picked to obtain simulated measurements
  at a specific set of the unknown parameters. In ous case, measurements at
  a specific set of parameters are simply evaluated using function "position",
  which we have defined to replace the simulation. Values of parameters are
  passed to this bloch through vector "parammom" from functions which require
  evaluation of chi square, e.g. optimization algoritms or tabulating
  functions. Function "meas" which will be called in a loop assignes values to
  components of the "measmom" vector: }
  ={pa:getvector["parammom",1]}
  ={pb:getvector["parammom",2]}
  ={pc:getvector["parammom",3]}
  ={pd:getvector["parammom",4]}
  ={i:1}
  while{ (i<=getvector["times",0])
  [
    meas{ ${i} ${position[getvector["times",i],pa,pb,pc,pd]} }
    ={i:i+1}
  ]}
  
  ={chi2buf:0}  *{ Temporary storage for the value of chi square will  }
  *{ Evaluation procedure for chi square: }
  ={i:1}
  while{ (i<=getvector["times",0])
  [
    ={chi2buf:chi2buf+
      sqr[getvector["measmom",i]-getvector["meas",i]]/sqr[getvector["sigma",i]]}
    *{ We have added a square of the difference between a component of the
    vector of simulated measurements and the appropriate component of the
    vector of measurements, divided by the appropriate component of the vector
    of measurement errors, to variable chi2buf. }
    ={i:i+1}
  ]}

  *{ With function chi2 we specify the value of chi square which will be
  returned to optimization algorithm. If option chi2spec would not be set,
  this function would be ignored and chi square would be automatically
  calculated vrom vectors "measmom", "meas" and "sigma" according to the
  least squares formula. }
  chi2{chi2buf}
}


write{"After analysis."\n}


*{ To test basic parts of the command file, we run analysis at specific
parameters. Function "analyse" runs the analysis at the parameters which are
in the vector "parammom". }

setvector{parammom 4 { $a  $b  $c  $d }}
analyse{}

*{ When everything looks O.K., we can run inverse analysis. In this case we
run the simplex algorithm which will minimize the value of chi square. The
algorithm successively changes the parameters, runs the argument block of the
"analysis" function and decides which values to take for the next guess
according to the obtained value of chi square and its values obtaines obtained
a few steps before. The first two arguments of function "inverse" are
specifications which specify the algorithm which will be used, the third
algorithm is a tolerance for the value of chi square in the minimum, the fourth
argument is maximum allowed number of iterations and the last parameter is
a matrix of five starting guesses: }

inverse{nd simplex 0.0001 500
  5 4
  { 1 1 : 6 }
  { 1 2 : 0 }
  { 1 3 : 7 }
  { 1 4 : 1 }
  
  { 2 1 : 2.1 }
  { 2 2 : 0 }
  { 2 3 : 7 }
  { 2 4 : 1 }
  
  { 3 1 : 6 }
  { 3 2 : 3.3 }
  { 3 3 : 7 }
  { 3 4 : 1 }
  
  { 4 1 : 6 }
  { 4 2 : 0 }
  { 4 3 : 0.2 }
  { 4 4 : 1 }
  
  { 5 1 : 6 }
  { 5 2 : 0 }
  { 5 3 : 7 }
  { 5 4 : 0.24 }
}

*{ We can evaluate different types of tebles, e.g. a table on a line between
two points; In the example below we make a table between the calculated
optimum ("param0") and a point which was used as a starting guess for inverse
analysis: }

printvector{paramopt}
fprintvector{paramopt}

tabline{lin 40 1 1
  { 4 
    {   ${getvector["paramopt",1]} ${getvector["paramopt",2]}
        ${getvector["paramopt",3]} ${getvector["paramopt",4]}
    }
  }
  { 4 { 6 0 7 1 } }
}