forum.alglib.net
http://forum.alglib.net/

Match GSL sample output
http://forum.alglib.net/viewtopic.php?f=2&t=3756
Page 1 of 1

Author:  cbellmore [ Sat Jun 18, 2016 12:17 am ]
Post subject:  Match GSL sample output

Hi, I'm trying to replicate the Levenberg-Marquardt optimization on an exponential curve from here: www .csse.uwa.edu.au/programming/gsl-1.0/gsl-ref_35.html in alglib with F#.

When I run what I have, it calculates the Jacobian and then terminates saying it met the gradient condition(4).
I figure I'm either calculating the Jacobian incorrectly or I'm using some other part of the library incorrectly. I've been stuck on this for a day now, and I figure it should be something obvious about how to setup the end conditions or a fundamental misunderstanding of the optimizer.

Here's the alglib code in F# I'm using now:
Code:
let expbF (expData : float[]) (x : float[]) (fi : float[]) myobj =
    for i in 0..fi.Length-1 do
        fi.[i] <- ((x.[0] * System.Math.Exp(-x.[1] * float i) + x.[2]) - expData.[i]) / 0.1

let expbDf (x : float[]) (fi : double[]) (jac : double[,]) myobj =
    let a = x.[0]
    let lambda = x.[1]
    for i in 0..fi.Length-1 do
        let t = float i   
        let e = System.Math.Exp(-lambda * t)
        jac.[i, 0] <- e / 0.1
        jac.[i, 1] <- -t * a * e / 0.1
        jac.[i, 2] <- 1.0 / 0.1

let testDV rawExpData =
    let mutable fitExp : float[] = [|1.0; 0.0; 0.0; |]
    let mutable state : alglib.minlmstate  =  (alglib.minlmstate())
    let mutable rep : alglib.minlmreport = alglib.minlmreport()

    alglib.minlmcreatevj(Array.length rawExpData, fitExp, &state)
    alglib.minlmsetcond(state, 0.0, 0.0, 0.0, 0)
    alglib.minlmoptimize(state, alglib.ndimensional_fvec(expbF rawExpData), expbDf, null, null )
    alglib.minlmresults(state, &fitExp, &rep)
    rep.terminationtype, alglib.ap.format(fitExp, 2)

I'm calling testDV with an array with the same values for the data that I'm generating in C++, but it doesn't matter because it never even calls expbF. I tried to follow the source, but I got really lost in the many different cases/gotos of minlmiteration.

And I'll in line my slightly modified GSL C++ code that I'm running to compare as well (VS2015):
Code:
namespace GSL {
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_fit.h>

}

struct data {
   size_t n;
   const double * y;
};

int expb_f(const GSL::gsl_vector * x, void *data, GSL::gsl_vector * f)
{
   size_t n = ((struct data *)data)->n;
   const double *y = ((struct data *)data)->y;

   double A = gsl_vector_get(x, 0);
   double lambda = gsl_vector_get(x, 1);
   double b = gsl_vector_get(x, 2);

   size_t i;

   for (i = 0; i < n; i++)
   {
      /* Model Yi = A * exp(-lambda * i) + b */
      double t = i;
      double Yi = A * exp(-lambda * t) + b;
      gsl_vector_set(f, i, (Yi - y[i]) / 0.1);
   }

   return GSL::GSL_SUCCESS;
}
int expb_df(const GSL::gsl_vector * x, void *data, GSL::gsl_matrix * J)
{
   size_t n = ((struct data *)data)->n;

   double A = GSL::gsl_vector_get(x, 0);
   double lambda = GSL::gsl_vector_get(x, 1);

   size_t i;

   for (i = 0; i < n; i++)
   {
      /* Jacobian matrix J(i,j) = dfi / dxj, */
      /* where fi = (Yi - yi)/sigma[i],      */
      /*       Yi = A * exp(-lambda * i) + b  */
      /* and the xj are the parameters (A,lambda,b) */
      double t = i;
      double s = 0.1;
      double e = exp(-lambda * t);
      gsl_matrix_set(J, i, 0, e / s);
      gsl_matrix_set(J, i, 1, -t * A * e / s);
      gsl_matrix_set(J, i, 2, 1 / s);
   }
   return GSL::GSL_SUCCESS;
}
int expb_fdf(const GSL::gsl_vector * x, void *data, GSL::gsl_vector * f, GSL::gsl_matrix * J)
{
   expb_f(x, data, f);
   expb_df(x, data, J);

   return GSL::GSL_SUCCESS;
}
double* getTestData() {
   const size_t n = 40;
   auto y = new double[n];

   GSL::gsl_rng_env_setup();

   auto type = GSL::gsl_rng_default;
   GSL::gsl_rng * r;
   r = GSL::gsl_rng_alloc(type);
   for (size_t i = 0; i < n; i++)
   {
      double t = i;
      y[i] = 1.0 + 5 * exp(-0.1 * t) + gsl_ran_gaussian(r, 0.1);
   };
   return y;
}
int main(void)
{
   double* const y = ExpFit::getTestData();

   const GSL::gsl_multifit_fdfsolver_type * T = GSL::gsl_multifit_fdfsolver_lmsder;
   GSL::gsl_multifit_fdfsolver * s = gsl_multifit_fdfsolver_alloc(T, 40, 3);
   struct data d = { 40, y };
   GSL::gsl_multifit_function_fdf f;
   f.f = &expb_f;
   f.df = &expb_df;
   f.fdf = &expb_fdf;
   f.n = 40;
   f.p = 3;
   f.params = &d;

   double x_init[3] = { 1.0, 0.0, 0.0 };
   GSL::gsl_vector_view x = GSL::gsl_vector_view_array(x_init, 3);
   gsl_multifit_fdfsolver_set(s, &f, &x.vector);

   int status;
   unsigned int iter = 0;
   do
   {
      iter++;
      status = gsl_multifit_fdfsolver_iterate(s);

      if (status)
      {
         break;
      }

      status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
   } while (status == GSL::GSL_CONTINUE && iter < 500);
   gsl_multifit_fdfsolver_free(s);
        delete[](y)
   return 0;
}


Thoughts?
At the very least this code should be a good example of this for others when it's done. I assumed this sort of thing would be more common.

Author:  Sergey.Bochkanov [ Mon Jun 20, 2016 2:47 pm ]
Post subject:  Re: Match GSL sample output

I did not tested your code yet, but it seems that you forgot to set function vector fi in expbDf. It calculates both function vector fi and Jacobian jac. Because you set only jacobian, optimizer treats function as constant and stops.

Author:  cbellmore [ Mon Jun 20, 2016 2:55 pm ]
Post subject:  Re: Match GSL sample output

Yep, that was it! It works as expected with that change.
Code:
let expbDf (expData : float[]) (x : float[]) (fi : double[]) (jac : double[,]) myobj =
    let a = x.[0]
    let lambda = x.[1]
    for i in 0..fi.Length-1 do
        fi.[i] <- ((x.[0] * System.Math.Exp(-x.[1] * float i) + x.[2]) - expData.[i]) / 0.1
    for i in 0..fi.Length-1 do
        let t = float i   
        let e = System.Math.Exp(-lambda * t)
        jac.[i, 0] <- e / 0.1
        jac.[i, 1] <- -t * a * e / 0.1
        jac.[i, 2] <- 1.0 / 0.1

I see it calculates fi in both functions in the example as well, I just didn't know how it was being used!

Thanks for that!

Page 1 of 1 All times are UTC
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group
http://www.phpbb.com/