[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Help-gsl] gsl minimizer
From: |
hanieh mahmoudian |
Subject: |
[Help-gsl] gsl minimizer |
Date: |
Mon, 29 Mar 2010 08:29:37 +0200 |
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;
> }
>