Re: Error in multifit_lmsder (?)

Re: Error in multifit_lmsder (?)

[Brian Gough]

Nicolai Hanssing writes:
 > I still think that it's probably my own routine,
 > but it follows the example ing gsl-ref closely...
 > 
 > The "lmsder" converges correctly the firste time, but when the routine is
 > called again [with the same parameters] it fails, and sends X= ['nan'
 > 'nan' 'nan'] to the solverfunctions f and df.

I think there was a problem in your code,

bjg@hppav:~$ gcc -g -Wall main.c fit_Asterix.c -lgsl -lgslcblas -lm
fit_Asterix.c: In function `print_state':
fit_Asterix.c:187: warning: implicit declaration of function `gsl_blas_dnrm2'
fit_Asterix.c:187: warning: double format, different type arg (arg 6)

This seemed to cause some random corruption of the stack, depending on
how the program was compiled.

Adding the header file

  #include <gsl/gsl_blas.h>

fixed the problem for me.

While looking for the problem I did find a floating point exception
that occurred for singular jacobian matrices, but this did not cause
the difference between the runs.

Regarding the GPL, as was previously said, there is information in the
new GPL FAQ linked off the front page at www.gnu.org.  In particular
you will want to look at the questions,

  "If a library is released under the GPL (not the LGPL), does that mean that
   any program which uses it has to be under the GPL?"

and

  "What if my school might want to make my program into its own proprietary
   software product?"

The direct url is http://www.gnu.org/copyleft/gpl-faq.html

regards
Brian Gough

Re: Error in multifit_lmsder (?)

[Brian Gough]

I've compiled the program but it will take me a while to study it.

Nicolai Hanssing writes:
 > Dear listmembers
 > 
 > I have struggeled with the "gsl_multifit_fdfsolver_lmsder", and it works
 > fine *once*. I have tried to look over the source of "gsl_multifit_", but
 > haven't found anything. 
 > 
 > I still think that it's probably my own routine,
 > but it follows the example ing gsl-ref closely...
 > 
 > The "lmsder" converges correctly the firste time, but when the routine is
 > called again [with the same parameters] it fails, and sends X= ['nan'
 > 'nan' 'nan'] to the solverfunctions f and df.
 > 
 > I do allocate space for the solver, initialize the solver, and I do free
 > the solver afterwards. And I cant find any reason for the behaivour of the
 > solver.
 > 
 > 3 calls as the following fails... [And I dont get it :-)]
 > solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha);  /* OK */ 
 > solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha);  /* FAILS */
 > solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha);  /* FAILS */
 > 
 > ANy insight will be much appreciated.
 > 
 > The problem solved is n>=3 equations: 0 = P0 - X0*sin(X1*P1+X2)
 > 
 > Regards 
 > 	Nicolai Hanssing
 > 	Denmark
 > 

Error in multifit_lmsder (?)

[Nicolai Hanssing]

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1: Type: text/plain, Size: 7138 bytes --]

Dear listmembers

I have struggeled with the "gsl_multifit_fdfsolver_lmsder", and it works
fine *once*. I have tried to look over the source of "gsl_multifit_", but
haven't found anything. 

I still think that it's probably my own routine,
but it follows the example ing gsl-ref closely...

The "lmsder" converges correctly the firste time, but when the routine is
called again [with the same parameters] it fails, and sends X= ['nan'
'nan' 'nan'] to the solverfunctions f and df.

I do allocate space for the solver, initialize the solver, and I do free
the solver afterwards. And I cant find any reason for the behaivour of the
solver.

3 calls as the following fails... [And I dont get it :-)]
solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha);  /* OK */ 
solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha);  /* FAILS */
solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha);  /* FAILS */

ANy insight will be much appreciated.

The problem solved is n>=3 equations: 0 = P0 - X0*sin(X1*P1+X2)

Regards 
	Nicolai Hanssing
	Denmark


I have included a compilable listing of the messy source:
------------

fit_Asterix.c:
--
/*
 * fit_Asterix.h for use in Asterix_main.c
 *
 * Nicolai Hanssing, nuh@kampsax.dtu.dk
 *
 * Masters Thesis 01, Foss-Electric
 *
 */

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_multifit_nlin.h>
#include <math.h>
#include "fit_Asterix.h"

/* "Private" prototyper */
int ulineval_f(const gsl_vector *x, void *params, gsl_vector *f);
int ulineval_df (const gsl_vector * x, void *params,gsl_matrix *df);
int ulineval_fdf (const gsl_vector *x, void *params, gsl_vector *f,
gsl_matrix *df);
int print_state(size_t iter, gsl_multifit_fdfsolver * s);

struct data {
		double L[5];
		gsl_vector *v_K;
	};

#define DEBUG 
#define DEBUG_DELRIN
#define DEBUG_SNS


int solve_nonlin_sys_multifit_nlin(gsl_vector *v_topp,gsl_vector
*v_alpha){
	/* Løser det overbestemte ulineære system
	 * med: 0 = L - a3sin(a4*k+a5)
	 */
	
	const gsl_multifit_fdfsolver_type * T;
	gsl_multifit_fdfsolver *s;
	gsl_multifit_function_fdf FDF;

	struct data params;
	int status, i;
  	size_t iter = 0;

	const size_t n = 3; /* number of equations >= 3 */
  	const size_t p = 3; /* number of parameters (3) */

  	double x_init[3] = { 0.1600, 0.0, 0.0 };

  	gsl_vector_view x = gsl_vector_view_array (x_init, p);
	params.v_K = gsl_vector_alloc(v_topp->size);
	gsl_vector_memcpy(params.v_K,v_topp);
	gsl_vector_scale (params.v_K, 0.01); /* To avoid badly conditioned
Jacobian */
	for (i=0;i<5;i++) params.L[i] = Delrin_lambda[i]/10000.0; /* To
avoid badly conditioned Jacobian */
	
	FDF.n = n; 	/* 3 equations 	*/
	FDF.p = p;	/* 3 variables 	*/
	FDF.params = (void *) &params;	
	FDF.f = &ulineval_f;
	FDF.df = &ulineval_df;
	FDF.fdf = &ulineval_fdf;
	
	T = gsl_multifit_fdfsolver_lmsder;
  	s = gsl_multifit_fdfsolver_alloc (T, n, p);
  	status = gsl_multifit_fdfsolver_set (s, &FDF, &x.vector);
	
	#ifdef DEBUG_SNS
	printf ("\n   status = %s\n", gsl_strerror (status));
    print_state (iter, s);
	#endif
	
	do {
      	iter++;
		status = gsl_multifit_fdfsolver_iterate (s);
    
		#ifdef DEBUG_SNS	
	  	if (!(int)fmod(iter,1)) {
			printf ("\n   status = %s\n", gsl_strerror
(status));
      		print_state (iter, s);
		}
    	#endif
	
	 	if (status)	break;
     	status = gsl_multifit_test_delta (s->dx, s->x, 0.000001,
0.000001);
    } while (status == GSL_CONTINUE && iter < 500);



	/* Scaling back */
	gsl_vector_set(v_alpha,2,gsl_vector_get(s->x,0)*10000);
	gsl_vector_set(v_alpha,3,gsl_vector_get(s->x,1)/100);
	gsl_vector_set(v_alpha,4,gsl_vector_get(s->x,2));
	
	gsl_vector_free(params.v_K);
	gsl_multifit_fdfsolver_free(s);
	
	return GSL_SUCCESS;
}


int ulineval_f(const gsl_vector * x, void *params, gsl_vector * f){

	size_t i;
	double f_val;
	long double Alpha3=gsl_vector_get(x,0);
	long double Alpha4=gsl_vector_get(x,1);
	long double Alpha5=gsl_vector_get(x,2);
	long double Ki,L;

	#ifdef DEBUG_SNS	
	printf("  X = %3.6e %3.6e %3.6e\n  ",
						gsl_vector_get(x,0),
						gsl_vector_get(x,1),
						gsl_vector_get(x,2));
	#endif

		
	for (i=0;i<3;i++){
		Ki = gsl_vector_get(((struct data *)params)->v_K,i);
		L = ((struct data *)params)->L[i];
		f_val= L - Alpha3 * sin(Alpha4 * Ki + Alpha5);
		gsl_vector_set(f,i,f_val);

		#ifdef DEBUG_SNS	
		printf("%e  ",f_val);
		#endif

	}

	#ifdef DEBUG_SNS	
	printf("\n");
	#endif

  	
	return GSL_SUCCESS;
}


int ulineval_df (const gsl_vector * x, void *params,gsl_matrix *df){

  	size_t i,h;
	double Alpha3=gsl_vector_get(x,0);
	double Alpha4=gsl_vector_get(x,1);
	double Alpha5=gsl_vector_get(x,2);
	double j[3];
	double Ki;

	for (i=0;i<3;i++){
		Ki=gsl_vector_get(((struct data *)params)->v_K,i);
		j[0] = - sin(Alpha4 * Ki + Alpha5);
		j[1] = - Alpha3 * cos(Alpha4 * Ki + Alpha5)*Ki;
		j[2] = - Alpha3 * cos(Alpha4 * Ki + Alpha5);
		for (h=0;h<3;h++) gsl_matrix_set (df, i, h, j[h]);
		#ifdef DEBUG_SNS	
		printf("\n |  %e %e %e | ",j[0],j[1],j[2]);
		#endif

	}
	#ifdef DEBUG_SNS	
	printf("\n ");
	#endif

	return GSL_SUCCESS;
}

int ulineval_fdf (const gsl_vector *x, void *params, gsl_vector *f,
gsl_matrix *df)
{
  ulineval_f(x, params, f);
  ulineval_df(x, params, df);

  return GSL_SUCCESS;
}

int print_state(size_t iter, gsl_multifit_fdfsolver * s){
  printf ("\niter: %3u x = % 15.8f % 15.8f % 15.8f |f(x)| = %g\n",
          iter,
          gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1),
          gsl_vector_get (s->x, 2), gsl_blas_dnrm2 (s->f));
}


--


main.c
--
/* Testprogram for debugging lmsder */

#include "fit_Asterix.h"
#include <gsl/gsl_vector.h>


int main(void){
	double topp[] = {525.80226278891245783598,
					620.25431012228921190399,
					725.58581627818944070896,
				  	854.57557708963383902301,
					960.14100685986863936705};

	gsl_vector_view v_topp;
	gsl_vector *v_alpha;
	
	v_alpha=gsl_vector_alloc(5);
	v_topp = gsl_vector_view_array(topp, 5);
	
	printf("\n\n  Calling solver\n");
	solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha);

/* The program fails on the following calls */	
	printf("\n\n\n  Calling solver, again...\n");
	solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha);
	printf("\n\n\n  Calling solver, again...\n");
	solve_nonlin_sys_multifit_nlin(&v_topp.vector,v_alpha);

	gsl_vector_free(v_alpha);
	
	return 0;
}
--

fit_Asterix.h:
--
#define Topp_samples	7	/* Antal målinger der skal 
							 * bruges ved fit
af toppunkt
							 * OBS: Skal være
ulige
							 */

#define Delrin_MinMax	5	/* Antal kendte lokale
							 * min/max
punkter:
							 */

static const double Delrin_lambda[]={903.408416839855,

960.556012816831,

1022.49708053353,
									0,

0}; /* i nm */
									
static const double Delrin_Arb[]={	1.61981074275959,

1.5127518464565,

1.59433851771899,
									0,

0};

/* OBS - Fundet via splinefit.... [wb04.m] */
/* 1.61981074275959           1.5127518464565          1.59433851771899
Arb */
/* 903.408416839855          960.556012816831          1022.49708053353
lambda */


#include <gsl/gsl_vector.h>

int solve_nonlin_sys_multifit_nlin(gsl_vector *v_topp,gsl_vector
*v_alpha);
--

fit_Asterix.c