#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<gsl/gsl_multimin.h>
#include<gsl/gsl_vector.h>
#define RND (double)rand()/RAND_MAX

int main(void){
	const int N=10;
	double x[N],y[N];
	for(int i=0;i<N;i++){
		x[i]=7.0*i/(N-1);
		y[i]=3.3*sin(1.1*x[i]-0.3)+0.5*RND;
	}

	double master(const gsl_vector *v, void* params){
		double A=gsl_vector_get(v,0);
		double K=gsl_vector_get(v,1);
		double B=gsl_vector_get(v,2);
		double s=0;
		for(int i=0;i<N;i++)
			s+=pow(A*sin(K*x[i]+B)-y[i],2);
		return s;
	}

	gsl_multimin_function F;
	F.f=&master;
	F.n=3;
	F.params=NULL;

	gsl_multimin_fminimizer *s
		= gsl_multimin_fminimizer_alloc
			(gsl_multimin_fminimizer_nmsimplex2rand,F.n);

	gsl_vector *v = gsl_vector_alloc(F.n);
	gsl_vector_set(v,0,3);
	gsl_vector_set(v,1,1);
	gsl_vector_set(v,2,1);

	gsl_vector *h = gsl_vector_alloc(F.n);
	gsl_vector_set(h,0,0.1);
	gsl_vector_set(h,1,0.1);
	gsl_vector_set(h,2,0.1);

	gsl_multimin_fminimizer_set(s,&F,v,h);

	int status,iter=0;
	do{
		iter++;
		status = gsl_multimin_fminimizer_iterate(s);
		if(status)break;	

		double size = gsl_multimin_fminimizer_size (s);
		status = gsl_multimin_test_size (size, 1e-2);

		if (status == GSL_SUCCESS)fprintf(stderr,"converged to minimum at\n");

	fprintf(stderr,"%i %g %g %g f() = %7.3f size = %.3f\n", 
              iter,
              gsl_vector_get (s->x, 0), 
              gsl_vector_get (s->x, 1), 
              gsl_vector_get (s->x, 2), 
              s->fval, size);
	}
	while(status == GSL_CONTINUE && iter < 200);

	double dt=0.1;
	for(double t=x[0];t<x[N-1]+dt;t+=dt){
              double A=gsl_vector_get (s->x, 0); 
              double K=gsl_vector_get (s->x, 1); 
              double B=gsl_vector_get (s->x, 2); 
		printf("%g %g\n",t,A*sin(K*t+B));
	}

		printf("\n\n");
		for(int i=0;i<N;i++)
		printf("%g %g\n",x[i],y[i]);
	
		

return 0;
}