Hi,
I am implementing a function as seen in the attachment. the parameters which are changing are a,b and c0. I have the x and y data from experimentation but sadly I am not getting any correct values of a,b and c. I am still not implementing the example with the hessian and gradient included. Should I implement that example? Is there anything wrong with my code? Please help :(
Attachment:
equation.png [ 3.19 KiB | Viewed 9631 times ]
Code:
#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#define PI 3.1415926535897
#include "interpolation.h"
using namespace alglib;
void function_cx_1_func(const real_1d_array &c, const real_1d_array &x, double &func, void *ptr)
{
// this callback calculates f(c,x)=exp(-c0*sqr(x0))
// where x is a position on X-axis and c is adjustable parameter
func = (1 / (c[0] * c[1])) * (c[2] - c[0] * 1 * sinh(sinh(c[2]) * exp(c[0] * x[0] * 2 *PI)));
//func= 1/ab*(c-a*sinh(sinh(c))*exp(a*x*2*pi);
//func = exp(-c[0]*pow(x[0],2));
}
void function_cx_1_grad(const real_1d_array &c, const real_1d_array &x, double &func, real_1d_array &grad, void *ptr)
{
// this callback calculates f(c,x)=exp(-c0*sqr(x0)) and gradient G={df/dc[i]}
// where x is a position on X-axis and c is adjustable parameter.
// IMPORTANT: gradient is calculated with respect to C, not to X
double temp;
grad[0] = -(c[2]);
grad[0] = grad[0] / (pow(c[0], 2) * c[1]);
temp=(2*PI*x[0]*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*sinh(c[2]));
temp=temp/c[1];
grad[0]=grad[0]-temp;
grad[1] = sinh( exp(2 * PI * c[0] * x[0]) * sinh(c[2])) / pow(c[1], 2) - c[2] / (c[0] * pow(c[1], 2));
grad[2] = 1 / (c[0] * c[1]) - ( exp(2 * PI * c[0] * x[0]) * cosh( exp(2 * PI * c[0] * x[0]) * sinh(c[2])) * cosh(c[2])) / c[1];
}
void function_cx_1_hess(const real_1d_array &c, const real_1d_array &x, double &func, real_1d_array &grad, real_2d_array &hess, void *ptr)
{
// this callback calculates f(c,x)=exp(-c0*sqr(x0)), grad G={df/dc[i]} and Hessian H={d2f/(dc[i]*dc[j])}
// where x is a position on X-axis and c is adjustable parameter.
// IMPORTANT: gradient/Hessian are calculated with respect to C, not to X
func = exp(-c[0]*pow(x[0],2));
hess[0][0]=(2*c[2])/(pow(c[0],3)*c[1])-(4*pow(PI,2)*pow(x[0],2)*exp(4*PI*c[0]*x[0])*sinh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*pow(sinh(c[2]),2))/c[1] - (4*pow(PI,2)*pow(x[0],2)*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*sinh(c[2]))/c[1];
hess[0][1]=c[2]/(pow(c[0],2)*pow(c[1],2)) + (2*PI*x[0]*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*sinh(c[2]))/pow(c[1],2);
hess[0][2]= - 1/(pow(c[0],2)*c[1]) - (2*PI*x[0]*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2]))/c[1] - (2*PI*x[0]*exp(4*PI*c[0]*x[0])*sinh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2])*sinh(c[2]))/c[1];
hess[1][0]= c[2]/(pow(c[0],2)*pow(c[1],2)) + (2*PI*x[0]*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*sinh(c[2]))/pow(c[1],2);
hess[1][1]= (2*c[2])/(c[0]*pow(c[1],3)) - (2*sinh(exp(2*PI*c[0]*x[0])*sinh(c[2])))/pow(c[1],3);
hess[1][2]=(exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2]))/pow(c[1],2) - 1/(c[0]*pow(c[1],2));
hess[2][0]= - 1/(pow(c[0],2)*c[1]) - (2*PI*x[0]*exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2]))/c[1] - (2*PI*x[0]*exp(4*PI*c[0]*x[0])*sinh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2])*sinh(c[2]))/c[1];
hess[2][1]= (exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*cosh(c[2]))/pow(c[1],2) - 1/(c[0]*pow(c[1],2));
hess[2][2]= - (exp(2*PI*c[0]*x[0])*cosh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*sinh(c[2]))/c[1] - (exp(4*PI*c[0]*x[0])*sinh(exp(2*PI*c[0]*x[0])*sinh(c[2]))*pow(cosh(c[2]),2))/c[1];
}
int main(int argc, char **argv)
{
//
// In this example we demonstrate exponential fitting
// by f(x) = exp(-c*x^2)
// using function value, gradient and Hessian (with respect to c)
//
real_2d_array x = "[[1],[2],[3],[4],[5],[6],[7]]";
real_1d_array y = "[438,469,563,797,1196,1410,1721]";
real_1d_array c = "[0.0,0.22,0.0]";
double epsf = 0;
double epsx = 0.000001;
ae_int_t maxits = 0;
ae_int_t info;
lsfitstate state;
lsfitreport rep;
double diffstep = 0.0001;
//
// Fitting without weights
//
lsfitcreatef(x, y, c, diffstep, state);
lsfitsetcond(state, epsf, epsx, maxits);
alglib::lsfitfit(state, function_cx_1_func);
lsfitresults(state, info, c, rep);
printf("%d\n", int(info)); // EXPECTED: 2
printf("%s\n", c.tostring(1).c_str()); // EXPECTED: [1.5]
//
// Fitting with weights
// (you can change weights and see how it changes result)
//
real_1d_array w = "[1,1,1,1,1,1,1]";
lsfitcreatewf(x, y, w, c, diffstep, state);
lsfitsetcond(state, epsf, epsx, maxits);
alglib::lsfitfit(state, function_cx_1_func);
lsfitresults(state, info, c, rep);
printf("%s\n", c.tostring(1).c_str()); // EXPECTED: [1.5](v
(void)getchar();
return 0;
}