Memo to use lmfit-2.4

Thanks for the author of lmfit. http://sourceforge.net/projects/lmfit/

We will have following code in lmfit, 

• lm_test.c - main( ), giving simple input data pair, fitting function and solve it
  
• lmmin.c - lm_initialize_control, lm_minimize, lm_lmdif, lm_lmpar, lm_qrfac, lm_qrsolv, lm_enorm
  
• lmmin.h -structure definition for lm_control_type
  
• lm_eval.c - lm_evaluate_default, lm_print_default 
  
• lm_eval.h - stucture definition for input data, &data in lm_minimize( )
  

Running lm_test.exe, then you will see a simple demo.

How to make the execution file?

Using linux, just follow the instruction which in lm_test.c

--------------------------------------------------

$ gcc -o lmtest -lm lmmin.c lm_eval.c lm_test.c

--------------------------------------------------

Using windows, include the lmmin.c lm_eval.c to your lm_test.c like, 

#include "lmmin.c"

#include "lm_eval.c"

The result will be

==============================================

C:\>lm_test.exe

modify or replace lm_print_default for less verbous fitting

starting minimization

  par:             1            1            1 => norm:     0.533515

determining gradient (iteration 1)

  par:             1            1            1 => norm:     0.533515

determining gradient (iteration 1)

  par:             1            1            1 => norm:     0.533515

determining gradient (iteration 1)

  par:             1            1            1 => norm:     0.533515

trying step in gradient direction

  par:             1            1            1 => norm:      1.29168

trying step in gradient direction

  par:             1            1            1 => norm:     0.247969

.....

.....

determining gradient (iteration 10)

  par:       6.25914      16.1219      6.75155 => norm:    0.0152438

trying step in gradient direction

  par:       6.25914      16.1219      6.75155 => norm:    0.0152438

terminated after 42 evaluations

  par:       6.25914      16.1219      6.75155 => norm:    0.0152438

  fitting data as follows:

    t[ 0]=        0.07 y=        0.24 fit=     0.24262 residue= -0.00261956

    t[ 1]=        0.13 y=        0.35 fit=    0.346227 residue=  0.00377306

    t[ 2]=        0.19 y=        0.43 fit=    0.423766 residue=  0.00623411

    t[ 3]=        0.26 y=        0.49 fit=    0.498947 residue= -0.00894747

    t[ 4]=        0.32 y=        0.55 fit=    0.555683 residue= -0.00568282

    t[ 5]=        0.38 y=        0.61 fit=    0.607558 residue=  0.00244152

    t[ 6]=        0.44 y=        0.66 fit=     0.65571 residue=  0.00429038

    t[ 7]=        0.51 y=        0.71 fit=    0.708095 residue=  0.00190477

    t[ 8]=        0.57 y=        0.75 fit=    0.750267 residue=-0.000266809

    t[ 9]=        0.63 y=        0.79 fit=    0.790257 residue=-0.000256849

    t[10]=        0.69 y=        0.83 fit=    0.828305 residue=  0.00169472

    t[11]=        0.76 y=        0.87 fit=    0.870492 residue=-0.000491909

    t[12]=        0.82 y=         0.9 fit=    0.904938 residue= -0.00493756

    t[13]=        0.88 y=        0.94 fit=    0.937935 residue=  0.00206541

    t[14]=        0.94 y=        0.97 fit=    0.969591 residue= 0.000409208

status: success (f) after 42 evaluations

==================================================

The input data is a pair data stored in t[  ] and y[  ], t and y are known.

p[  ] is unknown, a parameter vector for given function, called my_fit_function(  ),

In lm_test.c, the fiitting function g is

[p_0 * t_i + (1 - p_0 + p_1 +p_2) * t_i^2] / (1 + p_1 * t_i + p_2 * t_i^2)

we have to find a parameter vector that means p_0, p_1, and p_2 such that g(t_i)=y.

So unknowns are parameter vector, p[0], p[1], p[2], function g, and (t, y) are knowns. 

For example,

we have following data pair (t, y)

t                  y 

0.07           0.24 

0.13           0.35 

0.19           0.42 

0.26           0.49 

0.32           0.55 

0.38          0.61 

0.44           0.66 

0.51           0.71 

0.57           0.75 

0.63           0.79 

0.69           0.83 

0.76           0.87 

0.82           0.9 

0.88           0.94 

0.94           0.97

So we execute the procedue lm_minimize(m_dat, n_p, p, lm_evaluate_default, lm_print_default, &data, &control) and put appropriate parameters.

• m_dat : number of input, that means size of t[  ] here. Because t and y are data pair. They are in the same size.
  
• n_p : number of parameter, that means size of p[  ] here.
• lm_evaluate_default : 
• lm_print_default : show the processing messages and results
• &data : is lm_data_type, defining in lm_eval.h
• &control : a sturcture to control LM, for example, threshold of loop termination.

The initial value of control is defined by lm_initialize_control(  ) in lmmin.c

• maxcall = 100
  
• epsilon = 6.661338e-015 (Actually is 30*2.220446e-016 by definition statement)
  
• stepbound = 100
  
• ftol = 6.661338e-015
• xtol = 6.661338e-015
• gtol = 6.661338e-015

They are all defined by following definition statements,

==========================================================

/* machine-dependent constants from float.h */

#define LM_MACHEP     DBL_EPSILON   /* resolution of arithmetic */

#define LM_DWARF      DBL_MIN       /* smallest nonzero number */

#define LM_SQRT_DWARF sqrt(DBL_MIN) /* square should not underflow */

#define LM_SQRT_GIANT sqrt(DBL_MAX) /* square should not overflow */

#define LM_USERTOL    30*LM_MACHEP  /* users are recommened to require this */

==========================================================

DBL_EPSILON is 2.220446e-016 in my computer.

DBL_MIN is 2.225074e-308

DBL_MAX is 1.797693e+308

continue...