#include "longnam.h"
#include "math.h"
#include <stdio.h>
#include <string.h>
#include </users/hanieh/softwares/include/gsl/gsl_multimin.h>

double my_f (const gsl_vector *v, void *params) // keyparam *map
{
   int       npixel, i; 
   double    *I_model, *I_obs;
   double    deviation, least_square, delta, sum, average;
   double    R;
   double    *p = (double *)params;

npixel = 2000;

// allocating space for Intensity vectors
I_model = malloc(npixel * sizeof(double));
I_obs = malloc(npixel * sizeof(double));

// getting values to I_obs
for (i=0; i< npixel; i++)
   {   
//     printf ("observed pixel value: %f",map->new_pixel_value_chip1[i]);
     I_obs[i] = (i+0.01); //map->new_pixel_value_chip1[i];
   }

for (i=0; i< npixel; i++)
   {   
     I_model[i] = gsl_vector_get(v, i);
   }

//calculating mean value of pixels
sum = 0.0;
for (i=0; i< npixel; i++)
   {   
     sum = I_obs[i] + sum;
   }
average = sum / npixel;
//printf ("mean of observed pixel values:%f\n",average);


//calculating least-square value of pixels
least_square = 0.0;
for (i=0; i< npixel; i++)
   {
     delta = I_model[i] - I_obs[i];
     deviation = pow(I_obs[i]-average,2)/npixel;
     least_square = pow(delta,2)/deviation + least_square;
   }

R = 0.0;
for (i=0; i< npixel; i++)
   {
     R = I_model[i] * (I_model[i+1] - I_model[i]) + R;
   }

return  least_square + p[0] * R;
}


void my_df (const gsl_vector *v, void *params, gsl_vector *df) // keyparam *map
{
   int       npixel, i; 
   double    *I_model, *I_obs;
   double    deviation, dleast_square, delta, sum, average;
   double    dR;
   double    *p = (double *)params;

npixel = 2000;

// allocating space for Intensity vectors
I_model = malloc(npixel * sizeof(double));
I_obs = malloc(npixel * sizeof(double));

// getting values to I_obs
for (i=0; i< npixel; i++)
   {   
//     printf ("observed pixel value: %f",map->new_pixel_value_chip1[i]);
     I_obs[i] = (i+0.01)/100.0; //map->new_pixel_value_chip1[i];
   }

for (i=0; i< npixel; i++)
   {   
     I_model[i] = gsl_vector_get(v, i);
   }

//calculating mean value of pixels
sum = 0.0;
for (i=0; i< npixel; i++)
   {   
     sum = I_obs[i] + sum;
   }
average = sum / npixel;
//printf ("mean of observed pixel values:%f\n",average);

//calculating df/dI_model[i]
for (i=0; i< npixel; i++)
   {
     delta = I_model[i] - I_obs[i];
     deviation = pow(I_obs[i]-average,2)/npixel;
     dleast_square = 2 * delta/deviation;
     dR = -2.0 * I_model[i] + I_model[i+1];
     gsl_vector_set (df, i, dleast_square + p[0] * dR);
   }
}

void my_fdf(const gsl_vector *x,void *params, double *f, gsl_vector *df) // void *params, keyparam *map
{
*f = my_f(x, params);
my_df(x, params, df);
}


int main (void)
{
    const    gsl_multimin_fdfminimizer_type *T;
    int      status, i;
    int      npixel;
    double   par[1]={2.0};

npixel = 2000;
size_t iter = 0;

gsl_multimin_fdfminimizer *s;
gsl_vector *x; 
gsl_multimin_function_fdf my_func;

my_func.n = npixel; 
my_func.f = my_f;
my_func.df = my_df;
my_func.fdf = my_fdf;
my_func.params = par;

x = gsl_vector_alloc (npixel);
gsl_vector_set_all (x, 0.0);

T = gsl_multimin_fdfminimizer_conjugate_fr;
s = gsl_multimin_fdfminimizer_alloc (T, npixel);

gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);

do
  {
    iter++;
    status = gsl_multimin_fdfminimizer_iterate (s);

    if (status) { printf ("\n\n  **********  we failded  **********  \n\n"); break;}
    
    status = gsl_multimin_test_gradient (s->gradient, 1e-5);

    if (status == GSL_SUCCESS) printf ("minimum found at:\n");

    for (i=0; i<npixel; i++)
      {
         printf ("%5d X_%d %.5e f = %10.5f \n", iter, i, gsl_vector_get (s->x, i), s->f); 
      }
  }

while (status == GSL_CONTINUE && iter < 1000);

gsl_vector_free(x); 
//gsl_multinim_fdfminimizer_free(s);

return 0;
}