forum.alglib.net

ALGLIB forum
It is currently Wed Aug 10, 2022 2:20 pm

All times are UTC


Forum rules


1. This forum can be used for discussion of both ALGLIB-related and general numerical analysis questions
2. This forum is English-only - postings in other languages will be removed.



Post new topic Reply to topic  [ 3 posts ] 
Author Message
 Post subject: Match GSL sample output
PostPosted: Sat Jun 18, 2016 12:17 am 
Offline

Joined: Fri Jun 17, 2016 5:47 pm
Posts: 2
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.


Top
 Profile  
 
 Post subject: Re: Match GSL sample output
PostPosted: Mon Jun 20, 2016 2:47 pm 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 878
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.


Top
 Profile  
 
 Post subject: Re: Match GSL sample output
PostPosted: Mon Jun 20, 2016 2:55 pm 
Offline

Joined: Fri Jun 17, 2016 5:47 pm
Posts: 2
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!


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 3 posts ] 

All times are UTC


Who is online

Users browsing this forum: No registered users and 1 guest


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group