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.
Best Regards
Hanieh Mahmoudian
#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;
}
_______________________________________________
Help-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-gsl
