double function_temp(float p1, float p2, float p3, float p4, float var1)
{
   //float gkk = 0.95f;
   //float b = 3.f;         // Model Param
   float rho = 0.3f;         // Model Param
   //float sigma0 = 0.02f;      // Model Param
   //float m0 = 1.4f;          // Model Param
   int k = 6;             // Fixed Param  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   float lambda = -0.1f;      // Option Param
   float initPrice = 50.f;      // Option Param
//   float strikePrice   = 65.f; // Option Param
   float interest = 0.00018f;       // Option Param
   int noOfSum   = 300;          // Sym Param
   int width = 128;           // Sym Param
   int height = 256;// 256;           // Sym Param
   float temp;
return (double)clMonteCarloAMSM(            
             p3,//gkk,         // Model Param
             p1,//b,           // Model Param
             rho,         // Model Param
             p4,//sigma0,      // Model Param
             p2,          // Model Param
             k,             // Fixed Param
             lambda,      // Option Param
             initPrice,      // Option Param
             var1,//strikePrice, // Option Param
             interest,       // Option Param
             noOfSum,//noOfSum,          // Sym Param
             width,            // Sym Param
             height,           // Sym Param
             &temp);
//   return exp(-p1*pow(var1, 2));
  std::cout<<" temp = " << (double)temp <<std::endl;
}

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 = pow((float)x[0] - (float)c[0],2);

//   func = function_temp2((float)c[0],(float)x[0]);
   func = function_temp((float)c[0],(float)c[1],(float)c[2], (float)c[3], (float)x[0]);
   if (x[0]==50) {std::cout<<" c[0] = " << c[0]<<" c[1] = " << c[1] << " c[2] = " << c[2] <<" c[3] = " << c[3] << " x[0] = " << x[0] << " func = " << (double)func <<std::endl;}
}

int
main(int argc, char * argv[])
{

real_2d_array x = "[[30],[32],[34],[36],[38],[40],[42],[44],[46],[48],[50],[52],[54],[56],[58],[60],[62],[64],[66],[68],[70]]";
real_1d_array y = "[21.7538, 19.9783, 18.254, 16.5923, 15.004, 13.4978, 12.0822, 10.7645, 9.54836, 8.43455, 7.42138, 6.50722, 5.68772, 4.95832, 4.31183, 3.74198, 3.24245, 2.80582, 2.42521, 2.09441, 1.80822]";


   real_1d_array c = "[3.1, 1.38, 0.97, 0.03]";
   //real_1d_array c = "[3.0739, 1.4057, 0.9557]";

    real_1d_array bndl = "[2.0, 1.2, 0.9, 0.001]";
    real_1d_array bndu = "[10.0, 1.6, 0.9999, 0.06]";
   //real_1d_array s = "[0.8e+1, 1.0, 1.5]";
    real_1d_array s = "[8.0e-7, 0.4e-7, 0.25e-7, 0.05e-7]";

    double epsf = 0; 
    double epsx = 0.00001;  // condition on step size
    ae_int_t maxits = 0;
    ae_int_t info;
    lsfitstate state;
    lsfitreport rep;

    double diffstep =  1.0e-2;

    //
    // Fitting without weights
    //
    lsfitcreatef(x, y, c, diffstep, state);
    lsfitsetcond(state, epsf, epsx, maxits);
    lsfitsetbc(state, bndl, bndu);
   lsfitsetscale(state, s);
   //alglib::lsfitsetxrep(state, true);
    alglib::lsfitfit(state, function_cx_1_func);
    lsfitresults(state, info, c, rep);
   //printf("%d\n", int(state));
    std::cout<< "C = " << c.tostring(4).c_str() << std::endl;