[Help-gsl] gsl minimizer

[hanieh mahmoudian]

Hi

  Thanks for your answer. What I want to do is that this program is part of
a
  program that I'm writing to do some calculation on ACS CCD of Hubble
  telescope. In previous part of my program I had a vector of "I_obs" as an
  out put which its length is npixel=2048*4096~8600000, The problem is that
  this vector is a structure parameter but when I want to insert it to I_obs
  it did not work so I decided to simplify it as what it was in the program
to
  make sure that it works! This npixels should reach the above number, And I
  prefer to give it as a parameter not a number (to generalize my program
for
  other telescopes!)

 For this program I want to find "I_model" which gives me an accurate model
  of brightness distribution of sky on CCD. For that I have to minimize this
  function:

  { sum_(i=0,..,npixel-1) { (I_model[i] - I_obs[i])  / deviation(I_obs[i]) }
}
  + { p[0] * sum_(i,j=0,..,npixel-1) { I_model[i] * (I_model[i+1]  -
  I_model[i]) } }

 The first part is least_square in my program and the second one is R and
  p[0] is a parameter that I will give to the program. About the "free" part
  your right, I forgot to do that! But still the program breaks!
  Also when I add "gsl_multinim_fdfminimizer_free(s)" in the "int main" I
get
  this "dfminimizer.c:(.text+0x75c): undefined reference to
 `gsl_multinim_fdfminimizer_free'
  collect2: ld returned 1 exit status"

  With this mail I attached another version of my program with more
comments!
  Again thanks for your kind help!

 Cheers,
  Hanieh


  On Sat, Mar 27, 2010 at 10:09 PM, Thomas Weber
  <address@hidden>wrote:

> > On Fri, Mar 26, 2010 at 01:59:59PM +0100, hanieh mahmoudian wrote:
> > > Hi
> > >
> > > I wrote a program in C with use of GSL but the the problem is that
during
> > > the process it breaks! can you tell me what can be the problem? I
> > attached
> > > my program with this mail.
> >
> > The error code in status is '27', which seems to stand for (see
> > gsl_errno.h)
> >
> > GSL_ENOPROG  = 27,  /* iteration is not making progress towards solution
*/
> >
> > I'm not sure what you are trying to accomplish, could you give some
> > comments about the aim of you program?
> >
> > Some remarks (although I don't know whether your code is a simplified
> > example or your real program):
> >
> > 1) You shouldn't set npixel in each function to 2000. Instead, get the
> > lenght of the vector v and use that for npixel. Otherwise, your code
> > will be a maintenance nightmare once the vector length changes.
> >
> > 2) It _might_ be helpful to use gsl_vector for I_model and I_obs as
> > well. I don't think you will gain a lot, but it helps in seeing that
> > they are used as a vector.
> >
> > 3) I see quite some malloc()'s, but not nearly as many free()'s.
> >
> > 4) Your indentation is inconsistent. This makes reading your code
> > unnecessarily difficult.
> >
> >        Thomas
> >
> >

> Hi
On Sun, Mar 28, 2010 at 10:42:24AM +0200, hanieh mahmoudian wrote:

> #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 (which I
will use structure parameter for I_obs)
> {
>    int       npixel, i;
>    double    *I_model, *I_obs;
>    double    deviation, least_square, delta, sum, average;
>    double    R;
>    double    *p = (double *)params;
>
> npixel = 2000; // -> number of pixels in CCD 2048*4096
>
> // allocating space for Intensity vectors
> I_model = malloc(npixel * sizeof(double)); // I want to find this by
minimizind function below
> 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*2.0; //"map->new_pixel_value_chip1[i]"; -> the real
value for I_obs that I want to use for my program
>    }
>
> for (i=0; i< npixel; i++)
>    {
>      I_model[i] = gsl_vector_get(v, i);
>    }
>
> //calculating mean value of pixel values in CCD array
> 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 pixel values
> 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++)
>    {
>      if (i == npixel-1) R = -1.0 * I_model[npixel-1] * I_model[npixel-1] +
R;
>      R = I_model[i] * (I_model[i+1] - I_model[i]) + R;
>    }
> //R = I_model[0] * I_model[1] - I_model[0] * I_model[0]- I_model[1] *
I_model[1];
> free (I_model);
> free (I_obs);
>
> return  least_square + p[0] * R; //-> this is the function that I
discribed in my mail!
> }
>
>
> 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*2.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.0 * delta/deviation;
>      if (i == npixel-1) dR = -2.0 * I_model[npixel-1] * I_model[npixel-1]
+ I_model[npixel-2];
>      dR = -2.0 * I_model[i] + I_model[i+1];
>      gsl_vector_set (df, i, dleast_square + p[0] * dR);
>    }
> free (I_model);
> free (I_obs);
>
> }
>
> 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]={1.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, -1.0);
> gsl_vector_set (x, 0, 1.5);
> gsl_vector_set (x, 1, 1.5);
>
> 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-3);
>
>     if (status == GSL_SUCCESS) printf ("minimum found at:\n");
>
>     for (i=0; i<npixel; i++)
>       {
>          printf ("%2d X_%d %13e ||  f = %e \n", iter, i, gsl_vector_get
(s->x, i), s->f);
>       }
>   }
>
> while (status == GSL_CONTINUE && iter < 100);
>
> gsl_vector_free(x);
> gsl_multinim_fdfminimizer_free(s);
>
> return 0;
> }
>