/home/hgs/tap/mo

expfitbm


USAGE:
expfitbm @ <ins-file> [:<target>]

PURPOSE:
Fits a drift model to a tide residual with restored drift according to a previous (or ad-hoc) model. Involves non-linear fit of exponentials. The model equation is
d(t) = p1 + p2 t + p3 exp(-p4 (t - t1))
      + p5 + p6 (t - t5) ] H(t - t5)  [ + p7 + p8 (t - t7) ] H(t - t7) ...  + p15 * exp(-p16 (t - t15))

where H(t) is the Heaviside function. The second line is optional, the terms are invoked if the year of t5 etc. is > 2000

Solution demands  RMS(y - d) = min

Protocol output:
Where the output says MAPPED, the exponent parameter is written in units of
sampling intervals. Internally, the parameter is arctan(-p) + π/2

The slope parameters shown are in units of nm/s2 per sampling interval
Internally they are (g(end)-g(start))/n  p2 (i - 1)/n
                                         p5(i - i(t5))/n

INPUT INS-FILE:

1. FILE OPEN BLOCK
Input:    y = time series on unit 21
Output:   d = drift model on unit 31
        y-d = residual    on unit 32
        Exploration of the uncertainty margins (source code: "XPLORE")
              mc-file on unit IUMCV, default: no output
        Eigenvalue solution for uncertainty, (source code: "EIGEN")
              mc-file on unit IUMCE , default: no output


2. NAMELIST &param:
--------------------- ITERATION
iter            - integer - max. number of iterations DFPMIN (Numerical Recipes)
gtol            - real*8  - convergence control
p(NMP)          - real*8  - a priori solution (p(1) <- DClevel by program!)
pvar_start(NMP) - real*8  - for computing the p-uncertainties, the size of the
                            initial variation.
eps             - real*8  - iteration for uncertainties: stop when
                            update is less than eps*|p(j)|
qitermon        - logical - Monitor iteration of uncertainties

--------------------- INPUT
cal             - real*8  - Gravimeter's calibration value, needed to output
                            tsf-edit commands for drift synthesis
input_measure
                - char*16 - cal is the ratio nm/s² per Input_Measure; is put as a
                            comment in the tsf-edit lines. Output of tsf-edit
                            lines for Volt-units is only done if the string
                            contains 'V' somewhere.
nend            - integer - Far limit of the input time series. Max 1,000,000
                   
--------------------- MODEL

idate_exp, ihms_exp
                - integer - Y,M,D and H,M,S,F for the start of the exponent
                            (t1 in the equation above)
idate_second_step ...
                            (t2 in the equation above)

idate_second_exp, ihms_second_exp

idate_<nth>_step, ihms_<ntrh>_step
                           - where nth = third fourth fifth sixth

 

--------------------- VARIANCES AND XPLORE
del_p           - real*8  - For calculating the derivatives of prediciton
                            with respect to a parameter, dp(i) = p(i)*(1+del_p)
                            Default = 0.01

sig_d           - real*8  - Standard deviation of data-model fit.
                            Parameter is provided to expfit_jacobian.
                            Default -1.d0 (value<0) means that the
                            misfit**2 at every data point is used.

varf            - real*8  - Factor to remedy non-positive definiteness
                            calculation of the parameter covariances.
                            The smallest eigenvalue, multiplied with varf
                            is added to the diagonal before inverting it.

iumcv           - integer - output mc-file unit for XPLORE. Default: no output
                            To be declared in open-file block with option '<'

iumce           - integer - output mc-file unit for EIGEN. Default: no output
                            To be declared in open-file block with option '<'

qrms_raw       - logical - .false.: Predictions with PEF using eigenvector
                            is based on the PEF'd RMS.
                            .true.: ... based on the nonwhitened RMS.
                            Default: .false.
n_pef           - integer - -1: don't use a Prediction Error-Filter to
                                whiten the noise of the residual
                             0: use a trivial PEF (a diagnostic test)
                            >0: PEF order. A file with a PEF-bank must be
                                opened on unit 51.

AFTERMATH:
~/TD/plot-varying-expfit

Presently without options. Harmonizes with ~/TD/expfitm.ins :BISECT
Example result: http://holt.oso.chalmers.se/hgs/4me/ag-superc/expfitvar.png
PROBLEMS:
"Not a good start for j= n information"

Example:

     <Main-->>> Iterating p( 6), p, p0, pvar= -5.5818D+02 -5.0236D+02 -1.0000D-01 fdf=  2.6267D+02
     <Main-->T> Not a good start for j= 6 EXP2 DELAY -5.5818E+02 -5.0236E+02 -1.35831838E+01 -1.35806259E+01
where
    p(6)          = -1.070419D+03
    pvar_start(6) = -0.1

The message is issued when the product of the values underlined is greater than zero.

    p(6) far away?   pvar(6) too small or wrong sign?

You simply have to try with other starting values. Keep a little notebook.
Write down p(i), pvar(i), and the underlined values in the output. Find out the systematics.
Ë.g. further lowering the starting value p(i) abs.below 1.0E2 didn't make any impact.
Then try with an abs.greater pvar_start(i) .
With this choice:
    p(6)          =  -1.070419D+02    
    pvar_start(6) = -100.
the result became
     <Main-->>> Congratulations  for j= 6 EXP2 DELAY -5.5804E+02 -5.6362E+04 -1.35831838E+01  2.03760801E+01

EXAMPLE:


 PREPARE:
cd ~/Ttide/SCG
setenv FNSUB div/g120707-OPNEND-1h-atm.ph04.ts
setenv FACSUB 1.0 
tslist o/g090615-OPNEND-1h-atm.ra.ts -I -Bc2012,07,07, -o div/g120707-OPNEND-1h-atm.ra.ts
tslist o/g090615-OPNEND-1h-atm.ph04.ts -I -Bc2012,07,07, -o div/g120707-OPNEND-1h-atm.ph04.ts
tslist div/g120707-OPNEND-1h-atm.ra.ts -E ra+ph-deci.tse,C -I -o div/ra+ph03-120707-1d.ts
# instead of INIT below,
cd ~/TD
source expfitm-3rd-exp.ins


  INIT:
cd ~/TD
source expfit.setenv
setenv DT -1d

#which is:
setenv EXPFIT_INPUT div/ra+ph03-090615${DT}.ts
setenv EXPFIT_OUT -pef-ra+ph03-090615${DT}
setenv EXPFIT_TSE div/expfit-pef-090615${DT}.tse

  ~/TD/expfitbm.ins  
2XPEF>
21 ^ ${EXPFIT_INPUT}
31 < div/rm${EXPFIT_OUT}.ts
32 < div/rr${EXPFIT_OUT}.ts
33 B ${EXPFIT_TSE}
42 < div/reig_${EXPFIT_OUT}.mc
51 < div/rt${DT}.pef
52 < div/rrf${EXPFIT_OUT}.ts
   q
for expfitbm (bisection) (scldp was 24.)    (just a comment)
 &param
  n_pef=3
  gtol=1.d-6
  cal=-774.421

  p= 6.131548D+02 -1.355156D+02 -7.238925D+01 -1.034555D+03,
     9.754596D+01 -1.070419D+03 -8.153017D+01  1.053482D+02,
     3.030403D-02  4.216255D+00  1.842617D+00  2.577364D+01

  pvar_start = -100d0, -10.d0, 100.d0, 20.d0, 100.d0, -1.0, -2.0, -1.,
               -1000.0, -10.0, -100.0, -10.0, -100.0, -10.0
  eps=1.d-6
  idate_exp=        2009,06,15, ihms_exp=        0,0,0,0
  idate_second_exp= 2011,02,24, ihms_second_exp= 0,0,0,0
  idate_second_step=2011,02,24, ihms_second_step=0,0,0,0
  idate_third_step= 2012,02,06, ihms_third_step= 0,0,0,0
  idate_fourth_step=2013,10,24, ihms_fourth_step=0,0,0,0
  nend=99999
  iter=1000, qitermon=.true.
  input_measure='nm/s2'
  del_p=0.01d0
  varf=1.1d0, sig_d=-1.0d0, exp_sc=1.001d0
  iumcv=41, iumce=42
 &end
FURTHER READING:
http://holt.oso.chalmers.se/hgs/hgs.man/HOW.TO-model-drift.html

.bye