# forum.alglib.net

ALGLIB forum
 It is currently Thu Aug 15, 2024 9:40 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.

 Page 1 of 1 [ 3 posts ]
 Print view Previous topic | Next topic
Author Message
 Post subject: Match GSL sample outputPosted: Sat Jun 18, 2016 12:17 am

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

 Post subject: Re: Match GSL sample outputPosted: Mon Jun 20, 2016 2:47 pm

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

 Post subject: Re: Match GSL sample outputPosted: Mon Jun 20, 2016 2:55 pm

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

 Display posts from previous: All posts1 day7 days2 weeks1 month3 months6 months1 year Sort by AuthorPost timeSubject AscendingDescending
 Page 1 of 1 [ 3 posts ]

 All times are UTC

#### Who is online

Users browsing this forum: No registered users and 9 guests

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

Search for: