lmdemo.c

levenberg marquit算法的.C源码

  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=1.0;
  jac[j++]=-0.5/sqrt(p[3]);
}

/* Problem hatfldc (box constrained), minimum at (1.0, 1.0, 1.0, 1.0)
 * constri: p[i]>=0.0; (i=1..4)
 * constri+4: p[i]<=10.0; (i=1..4)
 */
void hatfldc(double *p, double *x, int m, int n, void *data)
{
register int i;

  x[0]=p[0]-1.0;

  for(i=1; i<m-1; ++i)
     x[i]=p[i-1]-sqrt(p[i]);

  x[m-1]=p[m-1]-1.0;
}

void jachatfldc(double *p, double *jac, int m, int n, void *data)
{
register int j=0;

  jac[j++]=1.0;
  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=0.0;

  jac[j++]=1.0;
  jac[j++]=-0.5/sqrt(p[1]);
  jac[j++]=0.0;
  jac[j++]=0.0;

  jac[j++]=0.0;
  jac[j++]=1.0;
  jac[j++]=-0.5/sqrt(p[2]);
  jac[j++]=0.0;

  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=1.0;
}

/* Equilibrium combustion problem, constrained nonlinear equation from the book by Floudas et al.
 * Minimum at (0.0034, 31.3265, 0.0684, 0.8595, 0.0370)
 * constri: p[i]>=0.0001; (i=1..5)
 * constri+5: p[i]<=100.0; (i=1..5)
 */
void combust(double *p, double *x, int m, int n, void *data)
{
  double R, R5, R6, R7, R8, R9, R10;

  R=10;
  R5=0.193;
  R6=4.10622*1e-4;
  R7=5.45177*1e-4;
  R8=4.4975*1e-7;
  R9=3.40735*1e-5;
  R10=9.615*1e-7;

  x[0]=p[0]*p[1]+p[0]-3*p[4];
  x[1]=2*p[0]*p[1]+p[0]+3*R10*p[1]*p[1]+p[1]*p[2]*p[2]+R7*p[1]*p[2]+R9*p[1]*p[3]+R8*p[1]-R*p[4];
  x[2]=2*p[1]*p[2]*p[2]+R7*p[1]*p[2]+2*R5*p[2]*p[2]+R6*p[2]-8*p[4];
  x[3]=R9*p[1]*p[3]+2*p[3]*p[3]-4*R*p[4];
  x[4]=p[0]*p[1]+p[0]+R10*p[1]*p[1]+p[1]*p[2]*p[2]+R7*p[1]*p[2]+R9*p[1]*p[3]+R8*p[1]+R5*p[2]*p[2]+R6*p[2]+p[3]*p[3]-1.0;
}

void jaccombust(double *p, double *jac, int m, int n, void *data)
{
register int j=0;
  double R, R5, R6, R7, R8, R9, R10;

  R=10;
  R5=0.193;
  R6=4.10622*1e-4;
  R7=5.45177*1e-4;
  R8=4.4975*1e-7;
  R9=3.40735*1e-5;
  R10=9.615*1e-7;

  for(j=0; j<m*n; ++j) jac[j]=0.0;

  j=0;
  jac[j]=p[1]+1;
  jac[j+1]=p[0];
  jac[j+4]=-3;

  j+=m;
  jac[j]=2*p[1]+1;
  jac[j+1]=2*p[0]+6*R10*p[1]+p[2]*p[2]+R7*p[2]+R9*p[3]+R8;
  jac[j+2]=2*p[1]*p[2]+R7*p[1];
  jac[j+3]=R9*p[1];
  jac[j+4]=-R;

  j+=m;
  jac[j+1]=2*p[2]*p[2]+R7*p[2];
  jac[j+2]=4*p[1]*p[2]+R7*p[1]+4*R5*p[2]+R6;
  jac[j+4]=-8;

  j+=m;
  jac[j+1]=R9*p[3];
  jac[j+3]=R9*p[1]+4*p[3];
  jac[j+4]=-4*R;

  j+=m;
  jac[j]=p[1]+1;
  jac[j+1]=p[0]+2*R10*p[1]+p[2]*p[2]+R7*p[2]+R9*p[3]+R8;
  jac[j+2]=2*p[1]*p[2]+R7*p[1]+2*R5*p[2]+R6;
  jac[j+3]=R9*p[1]+2*p[3];
}



int main()
{
register int i, j;
int problem, ret;
double p[5], // 6 is max(2, 3, 5)
	   x[16]; // 16 is max(2, 3, 5, 6, 16)
int m, n;
double opts[LM_OPTS_SZ], info[LM_INFO_SZ];
char *probname[]={
    "Rosenbrock function",
    "modified Rosenbrock problem",
    "Powell's function",
    "Wood's function",
    "Meyer's (reformulated) problem",
    "helical valley function",
    "Boggs & Tolle's problem #3",
    "Hock - Schittkowski problem #28",
    "Hock - Schittkowski problem #48",
    "Hock - Schittkowski problem #51",
    "Hock - Schittkowski problem #01",
    "Hock - Schittkowski (modified) problem #21",
    "hatfldb problem",
    "hatfldc problem",
    "equilibrium combustion problem"
};

  opts[0]=LM_INIT_MU; opts[1]=1E-15; opts[2]=1E-15; opts[3]=1E-20;
  opts[4]=LM_DIFF_DELTA; // relevant only if the finite difference jacobian version is used 

  /* uncomment the appropriate line below to select a minimization problem */
  problem=
		  //0; // Rosenbrock function
		  //1; // modified Rosenbrock problem
		  //2; // Powell's function
      //3; // Wood's function
		  4; // Meyer's (reformulated) problem
      //5; // helical valley function
#ifdef HAVE_LAPACK
      //6; // Boggs & Tolle's problem 3
      //7; // Hock - Schittkowski problem 28
      //8; // Hock - Schittkowski problem 48
      //9; // Hock - Schittkowski problem 51
#else // no LAPACK
#ifdef _MSC_VER
#pragma message("LAPACK not available, some test problems cannot be used")
#else
#warning LAPACK not available, some test problems cannot be used
#endif // _MSC_VER

#endif /* HAVE_LAPACK */
      //10; // Hock - Schittkowski problem 01
      //11; // Hock - Schittkowski (modified) problem 21
      //12; // hatfldb problem
      //13; // hatfldc problem
      //14; // equilibrium combustion problem
				
  switch(problem){
  default: fprintf(stderr, "unknown problem specified (#%d)! Note that some minimization problems require LAPACK.\n", problem);
           exit(1);
    break;
  case 0:
  /* Rosenbrock function */
    m=2; n=2;
    p[0]=-1.2; p[1]=1.0;
    for(i=0; i<n; i++) x[i]=0.0;
    ret=dlevmar_der(ros, jacros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    //ret=dlevmar_dif(ros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL);  // no jacobian
  break;

  case 1:
  /* modified Rosenbrock problem */
    m=2; n=3;
    p[0]=-1.2; p[1]=1.0;
    for(i=0; i<n; i++) x[i]=0.0;
    ret=dlevmar_der(modros, jacmodros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    //ret=dlevmar_dif(modros, p, x, m, n, 1000, opts, info, NULL, NULL, NULL);  // no jacobian
  break;

  case 2:
  /* Powell's function */
    m=2; n=2;
    p[0]=3.0; p[1]=1.0;
    for(i=0; i<n; i++) x[i]=0.0;
    ret=dlevmar_der(powell, jacpowell, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    //ret=dlevmar_dif(powell, p, x, m, n, 1000, opts, info, NULL, NULL, NULL);		// no jacobian
  break;

  case 3:
  /* Woods's function */
    m=4; n=6;
    p[0]=-3.0; p[1]=-1.0; p[2]=-3.0; p[3]=-1.0;
    for(i=0; i<n; i++) x[i]=0.0;
    ret=dlevmar_dif(wood, p, x, m, n, 1000, opts, info, NULL, NULL, NULL);  // no jacobian
  break;

  case 4:
  /* Meyer's data fitting problem */
    m=3; n=16;
    p[0]=8.85; p[1]=4.0; p[2]=2.5;
    x[0]=34.780;	x[1]=28.610; x[2]=23.650; x[3]=19.630;
    x[4]=16.370;	x[5]=13.720; x[6]=11.540; x[7]=9.744;
    x[8]=8.261;	x[9]=7.030; x[10]=6.005; x[11]=5.147;
    x[12]=4.427;	x[13]=3.820; x[14]=3.307; x[15]=2.872;
    //ret=dlevmar_der(meyer, jacmeyer, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian

   { double *work, *covar;
    work=malloc((LM_DIF_WORKSZ(m, n)+m*m)*sizeof(double));
    if(!work){
    	fprintf(stderr, "memory allocation request failed in main()\n");
      exit(1);
    }
    covar=work+LM_DIF_WORKSZ(m, n);

    ret=dlevmar_dif(meyer, p, x, m, n, 1000, opts, info, work, covar, NULL); // no jacobian, caller allocates work memory, covariance estimated

    printf("Covariance of the fit:\n");
    for(i=0; i<m; ++i){
      for(j=0; j<m; ++j)
        printf("%g ", covar[i*m+j]);
      printf("\n");
    }
    printf("\n");

    free(work);
   }

/* uncomment the following block to verify jacobian */
/*
   {
    double err[16];
    dlevmar_chkjac(meyer, jacmeyer, p, m, n, NULL, err); 
    for(i=0; i<n; ++i) printf("gradient %d, err %g\n", i, err[i]);
   }
*/

  break;

  case 5:
  /* helical valley function */
    m=3; n=3;
    p[0]=-1.0; p[1]=0.0; p[2]=0.0;
    for(i=0; i<n; i++) x[i]=0.0;
    ret=dlevmar_der(helval, jachelval, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    //ret=dlevmar_dif(helval, p, x, m, n, 1000, opts, info, NULL, NULL, NULL);  // no jacobian
  break;

#ifdef HAVE_LAPACK
  case 6:
  /* Boggs-Tolle problem 3 */
    m=5; n=5;
    p[0]=2.0; p[1]=2.0; p[2]=2.0;
    p[3]=2.0; p[4]=2.0;
    for(i=0; i<n; i++) x[i]=0.0;

    {
      double A[3*5]={1.0, 3.0, 0.0, 0.0, 0.0,  0.0, 0.0, 1.0, 1.0, -2.0,  0.0, 1.0, 0.0, 0.0, -1.0},
             b[3]={0.0, 0.0, 0.0};

    ret=dlevmar_lec_der(bt3, jacbt3, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, analytic jacobian
    //ret=dlevmar_lec_dif(bt3, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no jacobian
    }
  break;
  case 7:
  /* Hock - Schittkowski problem 28 */
    m=3; n=3;
    p[0]=-4.0; p[1]=1.0; p[2]=1.0;
    for(i=0; i<n; i++) x[i]=0.0;

    {
      double A[1*3]={1.0, 2.0, 3.0},
             b[1]={1.0};

    ret=dlevmar_lec_der(hs28, jachs28, p, x, m, n, A, b, 1, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, analytic jacobian
    //ret=dlevmar_lec_dif(hs28, p, x, m, n, A, b, 1, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no jacobian
    }
  break;
  case 8:
  /* Hock - Schittkowski problem 48 */
    m=5; n=5;
    p[0]=3.0; p[1]=5.0; p[2]=-3.0;
    p[3]=2.0; p[4]=-2.0;
    for(i=0; i<n; i++) x[i]=0.0;

    {
      double A[2*5]={1.0, 1.0, 1.0, 1.0, 1.0,  0.0, 0.0, 1.0, -2.0, -2.0},
             b[2]={5.0, -3.0};

    ret=dlevmar_lec_der(hs48, jachs48, p, x, m, n, A, b, 2, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, analytic jacobian
    //ret=dlevmar_lec_dif(hs48, p, x, m, n, A, b, 2, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no jacobian
    }
  break;
  case 9:
  /* Hock - Schittkowski problem 51 */
    m=5; n=5;
    p[0]=2.5; p[1]=0.5; p[2]=2.0;
    p[3]=-1.0; p[4]=0.5;
    for(i=0; i<n; i++) x[i]=0.0;

    {
      double A[3*5]={1.0, 3.0, 0.0, 0.0, 0.0,  0.0, 0.0, 1.0, 1.0, -2.0,  0.0, 1.0, 0.0, 0.0, -1.0},
             b[3]={4.0, 0.0, 0.0};

    ret=dlevmar_lec_der(hs51, jachs51, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, analytic jacobian
    //ret=dlevmar_lec_dif(hs51, p, x, m, n, A, b, 3, 1000, opts, info, NULL, NULL, NULL); // lin. constraints, no jacobian
    }
  break;
#endif /* HAVE_LAPACK */

  case 10:
  /* Hock - Schittkowski problem 01 */
    m=2; n=2;
    p[0]=-2.0; p[1]=1.0;
    for(i=0; i<n; i++) x[i]=0.0;
    //ret=dlevmar_der(hs01, jachs01, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    {
      double lb[2], ub[2];

      lb[0]=-DBL_MAX; lb[1]=-1.5;
      ub[0]=ub[1]=DBL_MAX;

      ret=dlevmar_bc_der(hs01, jachs01, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    }
    break;
  case 11:
  /* Hock - Schittkowski (modified) problem 21 */
    m=2; n=2;
    p[0]=-1.0; p[1]=-1.0;
    for(i=0; i<n; i++) x[i]=0.0;
    //ret=dlevmar_der(hs21, jachs21, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    {
      double lb[2], ub[2];

      lb[0]=2.0; lb[1]=-50.0;
      ub[0]=50.0; ub[1]=50.0;

      ret=dlevmar_bc_der(hs21, jachs21, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    }
    break;
  case 12:
  /* hatfldb problem */
    m=4; n=4;
    p[0]=p[1]=p[2]=p[3]=0.1;
    for(i=0; i<n; i++) x[i]=0.0;
    //ret=dlevmar_der(hatfldb, jachatfldb, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    {
      double lb[4], ub[4];

      lb[0]=lb[1]=lb[2]=lb[3]=0.0;

      ub[0]=ub[2]=ub[3]=DBL_MAX;
      ub[1]=0.8;

      ret=dlevmar_bc_der(hatfldb, jachatfldb, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    }
    break;
  case 13:
  /* hatfldc problem */
    m=4; n=4;
    p[0]=p[1]=p[2]=p[3]=0.9;
    for(i=0; i<n; i++) x[i]=0.0;
    //ret=dlevmar_der(hatfldc, jachatfldc, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    {
      double lb[4], ub[4];

      lb[0]=lb[1]=lb[2]=lb[3]=0.0;

      ub[0]=ub[1]=ub[2]=ub[3]=10.0;

      ret=dlevmar_bc_der(hatfldc, jachatfldc, p, x, m, n, lb, ub, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    }
    break;
  case 14:
  /* equilibrium combustion problem */
    m=5; n=5;
    p[0]=p[1]=p[2]=p[3]=p[4]=0.0001;
    for(i=0; i<n; i++) x[i]=0.0;
    //ret=dlevmar_der(combust, jaccombust, p, x, m, n, 1000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    {
      double lb[5], ub[5];

      lb[0]=lb[1]=lb[2]=lb[3]=lb[4]=0.0001;

      ub[0]=ub[1]=ub[2]=ub[3]=ub[4]=100.0;

      ret=dlevmar_bc_der(combust, jaccombust, p, x, m, n, lb, ub, 5000, opts, info, NULL, NULL, NULL); // with analytic jacobian
    }
    break;
  } /* switch */
  
  printf("Results for %s:\n", probname[problem]);
  printf("Levenberg-Marquardt returned %d in %g iter, reason %g\nSolution: ", ret, info[5], info[6]);
  for(i=0; i<m; ++i)
    printf("%.7g ", p[i]);
  printf("\n\nMinimization info:\n");
  for(i=0; i<LM_INFO_SZ; ++i)
    printf("%g ", info[i]);
  printf("\n");

  return 0;
}