this self contained program

Hi I have found an implementation of an algorithm called the Levenberg Marquardt algorithm that states that it is self contained.
I am trying to understand an implementation of it and most of the web and HUGE but this one seems nice and small however I cannot get it to compile (using codeblocks). I have attached what I think is required, and below is the page I downloaded it from.
Could someone please have a look and see if they can get it to compile.
http://joachimwuttke.de/lmfit/
(download link under resources)

Thanks
Alex

What operating system are you using? The normal ./configure/make/make install commands seemed to work for me, on Linux.

Jim

Oh! Im on windows. Can I run a similar command from here.
out of interest what does it give you? .o files or something? I need to be able to compile this into an embedded program.
I am looking for a standalone C implementation of this algorithm but there seem to be few out there - there are ones built into MUCH MUCH more complex libraries, but I dont really have the C expertise to extract them
Alex

Haven't used codeblocks before but I'm guessing that you'd have to start by adding all the source files to a new project before attempting to build the application. Have you done that? If not then do so and try to compile/link and report the exact error messages you get if any.

Quote:

---

Originally Posted by a.mlw.walker

Oh! Im on windows. Can I run a similar command from here.
out of interest what does it give you? .o files or something?

---

Yes. On Windows you can install MinGW, which gives you the MSYS environment, allowing you to run scripts like configure.

The result is a bunch of .o files which are combined into an archive file (.a). In GCC, you normally name the file something like libfoo.a and then you can compile a program using -lfoo to link with that library. Compilers for embedded systems will have a similar facility, but this varies from compiler to compiler.

After a quick look everything looks fairly standard. It looks like you should be able to compile it with a Windows compiler without any problems. You should be able to add the source files to your project and compile it normally. Just be sure to put the source and include files into your project source directory.

Jim

Any chance someone can post the "lib"' s because Im not sure to compile the rest of the files, except the ones I attached.
Yeah I made a new project in codeblocks, added the files in my zip (attached above) and it wouldnt compile. I can see that there must be more code somewhere because there aren't any algorithms that do what levenberg marquardt does in the files I can see (the .c ones I attached).

Post the compiler errors! Please use code tags for posting the errors.

Tim S.

Your file lmmin.c DOES NOT have the correct contents. It is the same as lmcurve.h

Compiles for me using MinGW GCC C under Code::Blocks when I used the correct code in the file lmmin.c.

Tim S.

Bang on, I dont know how I managed that. Thanks
Alex

Hey guys,
sorry to bring this up again, I am trying to change the function that the example is fitting for my equation:

Code:

---

double myfunction(double Actual_Output[9], double *p){
double Fitted_Curve[5] = {0};
int i;
double invtrig = 0, x = 0;


  for (i = 0; i <sizeof(Actual_Output)-1; i++)
  {
        invtrig = arcsinh(sinh(p[2])*exp(p[0]*(i+1)*2*M_PI));
        Fitted_Curve[i] = (1/(p[0]*p[1]))*(p[2]-invtrig);
    x += Actual_Output[i+2]-Fitted_Curve[i];

    }
return x;
}
/* model function: a quadratic function p0 + p1 * (t-p2)^2 */

double f( double t, const double *p )
{
    return p[0] + p[1]*(t-p[2])*(t-p[2]);
}

---

the example function is f, mine is myfunction...
I thought I would have to pass the function my data somehow, but they seem to only pass a double and a pointer to a double array.
Can anyone help me to implement my algorithm instead of theirs?
Thanks
Alex

Please show where you are using your function. If you are trying to use your function like they are using their function then your function's signature must match their function's signature. You can change the interior of the function but it must return the same type of variable and the parameters must be exactly the same.

Jim

You may want to study this link: The Function Pointer Tutorial. It might help explain how that function is actually being used.

Jim

Quote:

---

Originally Posted by a.mlw.walker

I thought I would have to pass the function my data somehow, but they seem to only pass a double and a pointer to a double array.
Can anyone help me to implement my algorithm instead of theirs?

---

As I understand it after reading the example program (curve1.c) the function "f" should describe the type of function which you think will best fit your data points. For example if you think there is a linear relationship "f" could look like:

Code:

---

double f(double t, const double *p)
{
    return p[0] * t + p[1]
}

---

If you think a logarithmic relationship works best it could look like:

Code:

---

double f(double t, const double *p)
{
    return log(t) + p[0]    // needs <math.h>
}

---

and so on. So it's up to you what function you need.

You have to pass your data points to lmcurve_fit() as two arrays (4th and 5th parameter). Read again the example program and experiment with it (change the parameters) to see how it works.

Code:

---

double myfunction(double Actual_Output[9], double *p){
...
for (i = 0; i <sizeof(Actual_Output)-1; i++)

---

This is definitely wrong. Your first parameter is an array, hence the function only gets a pointer to the first element. Then sizeof returns the size of this pointer which is system dependent (on my system it's 4 bytes). You can't determine the size of the array (i.e. the number of its element) inside the function. You would have to pass it as an additional parameter.

Bye, Andreas

Hi guys,
so after some studying, I think that for my problem I actually want surface1.c as the basis to manipulate for my algorithm, not sure though, just think that.
I have adjusted the program so that it now compiles and runs with my function however I dont think I have quite understood my way around all the function pointers and the struct thats in there too.
The data I have that I am trying to fit is:

Code:

---

  double tx[5] = {1, 2, 3, 4, 5};
    double tz[5] = {1.2, 2.693, 4.325, 6.131, 8.125};

---

On a graph tx would be on the x axis, and tz on the y. Im not sure what the 'tx' array in the example is.
So I am trying to fit the above data to the algorithm:

Code:

---

double f(double tz, const double *p){
double Fitted_Curve[5] = {0};
int i;
double invtrig = 0, x = 0;


  for (i = 0; i <sizeof(tz)-1; i++)
  {
        invtrig = arcsinh(sinh(p[2])*exp(p[0]*(i+1)*2*3.145));
        Fitted_Curve[i] = (1/(p[0]*p[1]))*(p[2]-invtrig);
    x += pow((tz[i]-Fitted_Curve[i]), 2);

    }
return x;
}

---

However it wont compile with the above, because I am asking for tz[i]. The following algorithm it will compile and run with, however it produces incorrect results, as its not "fitting" the correct curve (I dont think)

Code:

---

double f(double tz, const double *p){
double Fitted_Curve = 0;//[5] = {0};
int i;
double invtrig = 0, x = 0;


 // for (i = 0; i <sizeof(tz)-1; i++)
 // {
        invtrig = arcsinh(sinh(p[2])*exp(p[0]*(i+1)*2*3.145));
        Fitted_Curve = (1/(p[0]*p[1]))*(p[2]-invtrig);
    x = pow((tz-Fitted_Curve), 2);

  //  }
return x;
}

---

I dont really need tx because its just 1,2,3,4,5 which I planned on using i+1 for.

Please can someone have a look and see if they can simplify the algorithm, or at least help me to simplify it so that I can "fit this equation"?
Thanks

Quote:

---

Originally Posted by a.mlw.walker

so after some studying, I think that for my problem I actually want surface1.c as the basis to manipulate for my algorithm, not sure though, just think that.

---

From the header of surface1.c

Quote:

---

Contents: Example for fitting data y_i(t_i) by a function f(t;p), where each t_i is a vector of dimension k=2.

---

If your dataset is a vector (tx, tz), what are the y values?
And if you really have a 2D-vector, your callback function (f) needs three parameters (see the model function in surface1.c)

Quote:

---

Originally Posted by a.mlw.walker

Please can someone have a look and see if they can simplify the algorithm, or at least help me to simplify it so that I can "fit this equation"?

---

I don't understand what your algorithm should produce. Does your algorithm generate different functions and you want to find the best?

Curve fitting is about finding a function which fits a set of data points best. The library uses the Levenberg-Marquardt algorithm for that. You have to provide a dataset and a general function (with some modifiable parameters). Then the library finds the parameter combination which fits best.

So I don't see where your algorithm fits here (especially the variable "Fitted_Curve" suggests that you try to calculate the curve yourself).

Bye, Andreas

Hi.
Thanks for responding. I have attached the equation I am trying to fit. Now you mention it, I wonder whether I am supposed to supply the equation within the curly braces, and the rest of the algoithm does the sum of the values.
The reason I thought surface may be closer to my problem is because '1D' (curve1) didnt seem right, I can only visualise a graph in 2D (x and y) so I thought that made more sense - also it was passing more variables/parameters to the function so I thought it would make it more obvious how to change the current code.
The values I posted before:

Code:

---

double tx[5] = {1, 2, 3, 4, 5};
    double tz[5] = {1.2, 2.693, 4.325, 6.131, 8.125};

---

are k, and tk respectively in the equation. I am trying to fit for a, b and c0 - p[0], p[1] and p[2].
Is that clearer?

Now I think we come closer to the solution.

Your image of the equation looks like the formula for the least square method (finding the minimum of the sum of the squared residuals).
Then x is your k, y is your tk, beta is (a,b,c0) and finally your model function is the part inside the braces after the first minus sign.

Are you now able to write the callback function? You should use curve1.c as an example and adjust it to your actual values.

Bye, Andreas

So what Im not sure about is how much of the equation I am supposed to put in the call back function - i.e do I only put the bit inside the curly braces or do I need to put everything
i/e do I put
y - f(y, Beta)
or do I put
SUM(y-f(y,beta))^2
?

So I have had a go, it compiles but I just want to check whether you think I have done the right thing.
i have adjusted f to this:

Code:

---

double f( double t, const double *p )
{
double invtrig = 0;
    invtrig = arcsinh(sinh(p[2])*exp(p[0]*(t)*2*M_PI));
    return (1/(p[0]*p[1]))*(p[2]-invtrig);
}

---

which is calculating just the model function. It does not include the tk (y) value as I assume then that the algorithm deals with that.
So I am just returning the value of the model function to the levenberg marquardt algorithm.
In main, I have set:

Code:

---

    double t[5] = {1, 2, 3, 4, 5};
    double y[5] = {1.2, 2.693, 4.325, 6.131, 8.125};

---

Is this right?

EDIT: It compiles, and also thinks it converged, however the parameters it calculates are way out from the expected values for the test data.
Does the LM algorithm know that it is meant to be a least squares fit?
ie based the general LM equation for least squares, I am only feeding it f(x, beta) I am kind of assuming it knows to subtract that from the other data set, and square it and sum the values....

Code:

---

invtrig = arcsinh(sinh(p[2])*exp(p[0]*(t)*2*M_PI));

---

Where is arcsinh() defined? The standard function (C99) from math.h is called asinh().

Quote:

---

EDIT: It compiles, and also thinks it converged, however the parameters it calculates are way out from the expected values for the test data.

---

I get the following result:

Code:

---

...
Results:
status after 68 function evaluations:
  success (the relative error in the sum of squares is at most tol)
obtained parameters:
  par[0] =      -1.0458
  par[1] =      4.07939
  par[2] =    -149.809
obtained norm:
      0.730656
fitting data as follows:
  t[ 0]=          1 y=        1.2 fit=    1.54018 residue=  -0.340182
  t[ 1]=          2 y=      2.693 fit=    3.08036 residue=  -0.387364
  t[ 2]=          3 y=      4.325 fit=    4.62055 residue=  -0.295545
  t[ 3]=          4 y=      6.131 fit=    6.16073 residue=  -0.0297273
  t[ 4]=          5 y=      8.125 fit=    7.70091 residue=    0.424091

---

Bye, Andreas

Yeah I think I get the same (at other computer). The results I get in Matlab, and the results I am supposed to get are
a: 0.02
b: 0.08
c0:-0.39

Im not sure if Ive adjusted the function correcty and hence why Im not getting the right results.

Quote:

---

Originally Posted by a.mlw.walker

Yeah I think I get the same (at other computer). The results I get in Matlab, and the results I am supposed to get are
a: 0.02
b: 0.08
c0:-0.39

---

So what model function did you use in Matlab?

Bye, Andreas

This is the matlab function. If you dont know Matlab, the dot is a bit like a "for" loop for matrices...

Quote:

---

a=params(1)
b=params(2)
c=params(3)

% c = c'
Fitted_Curve = (1/(a*b))*(c-asinh(sinh(c)*exp(a*Input*2*pi)));
Error_Vector=Actual_Output-Fitted_Curve ;
sse=sum(Error_Vector.^2);

---

Basically the same thing - the point is it just depends on how the method that is using it, wants it?

Quote:

---

Originally Posted by a.mlw.walker

Fitted_Curve = (1/(a*b))*(c-asinh(sinh(c)*exp(a*Input*2*pi)));

---

The weird thing is, if I use the parameters you've got and calculate the values for your inputs, then I get the following result:

Code:

---

#include <stdio.h>
#include <math.h>

int main(void)
{
    int i;
    double result;
    const double a = 0.02, b = 0.08, c = -0.39;
    for (i = 1; i <= 5; i++)
    {
        result = (1 / (a * b)) * (c - asinh(sinh(c) * exp(a * i * 2 * 3.1415)));
        printf("%d - %g\n", i, result);
    }
    return 0;
}

---

Code:

---

$ gcc -o bar bar.c -lm
$ ./bar
1 - 30.7834
2 - 64.9449
3 - 102.672
4 - 144.103
5 - 189.316

---

So they are far away from your test data.

I think it would really help if you post the exact assignment you are doing.

Bye, Andreas

Do you have Matlab?
If so, I have attached the necessary Matlab files to get the Matlab results, however if you dont I have also pasted the output of Matlab.
I think the problem is in the way that the function has been defined, but Im not sure. In Matlab the function is defined as the entire least squares equation, not just the right hand side of the levenberg equation.
Matlab function:

Code:

---

function sse=myfit(params,Input,Actual_Output)

a=params(1)
b=params(2)
c=params(3)

% c = c'
Fitted_Curve = (1/(a*b))*(c-asinh(sinh(c)*exp(a*Input*2*pi)));
Error_Vector=Actual_Output-Fitted_Curve ;
sse=sum(Error_Vector.^2);

---

The function above is then passed to Matlab's Levenberg Marquardt algorithm here:

Code:

---

Data = [1.2 2.693 4.325 6.131 8.125]%paper data
t = [1 2 3 4 5];
%%
Starting=[1, 1, -1];
options=optimset('Display','iter');
Estimates=fminsearch(@myfit,Starting,options,t,Data)

---

Estimates being the vector of the parameters returned which come out as

Quote:

---

Estimates =

0.0250 1.8944 -0.3690

---

The Matlab code then plots the graph that would be the function with those parameters so the reader can see the quality of the fit to the original data (attached).

Matlab code:

Code:

---

Data = [1.2 2.693 4.325 6.131 8.125]%paper data
t = [1 2 3 4 5];
%%
Starting=[1, 1, -1];
options=optimset('Display','iter');
Estimates=fminsearch(@myfit,Starting,options,t,Data)
%%
a = Estimates(1)
b = Estimates(2)
c = Estimates(3)
tt = 0:0.001:5;

% To check the fit
% subplot(2,2,1)
plot(t,Data,'*')
hold on
plot(tt,(1/(a*b))*(c-asinh(sinh(c)*exp(a*tt*2*pi))),'r');
xlabel('goodness of fit graph');

---

and the code in the file "myfit.m"

Code:

---

function sse=myfit(params,Input,Actual_Output)

a=params(1)
b=params(2)
c=params(3)

% c = c'
Fitted_Curve = (1/(a*b))*(c-asinh(sinh(c)*exp(a*Input*2*pi)));
Error_Vector=Actual_Output-Fitted_Curve ;
sse=sum(Error_Vector.^2);

---

So I am confused as to why this algorihtm (the c one) produces such wildly different results...

Edit, the second graph is using the results that you got - the one that has a red "step change".

So, here are probably my final thoughts after some testing and researching:

Quote:

---

Originally Posted by a.mlw.walker

The function above is then passed to Matlab's Levenberg Marquardt algorithm here:

---

You use Matlab's fminsearch function which isn't using the Levenberg-Marquardt algorithm but the Nelder-Mead simplex algorithm.

This could explain the differences.
It also seems that the Levenberg-Marquardt algorithm is more sensitive to the initial starting parameters.
If for example I started with the initial parameters 0.1 / 1 / -0.1 I get as end result

Code:

---

a = 0.0396
b = 2.756
c = -0.537

---

which will produce a similar plot as your values:Attachment 12502

As you can see, in the range [0:5] the plots are very similar.

Bye, Andreas

Still quite interesting, because you can request in Matlab for fMinSearch to use Levenberg - I should have pointed that

Quote:

---

options = optimset('MaxFunEvals', 10000, 'MaxIter', 10000);
options = optimset(options, 'LevenbergMarquardt', 'on');
Estimates=fminsearch(@myfit,Starting,options,t,Dat a);

---

However I think you are right. Its to do with the sensitivity of the algorithm. I tried with different data, and it got exactly the same as Maltab did for the parameters.
Thanks for your help I think I can take it from here.