#include <math.h>
     #include <stdlib.h>
     #include <string.h>
     #include <gsl/gsl_siman.h>
#include <gsl/gsl_vector.h>
     
     /* set up parameters for this simulated annealing run */
     
     /* how many points do we try before stepping */
     #define N_TRIES 200             
     
     /* how many iterations for each T? */
     #define ITERS_FIXED_T 10
     
     /* max step size in random walk */
     #define STEP_SIZE 1.0            
     
     /* Boltzmann constant */
     #define K 1.0                   
     
     /* initial temperature */
     #define T_INITIAL 0.008         
     
     /* damping factor for temperature */
     #define MU_T 1.003              
     #define T_MIN 2.0e-10
     
     gsl_siman_params_t params 
       = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
          K, T_INITIAL, MU_T, T_MIN};
     
     /* now some functions to test in one dimension */

     double yvals[10]={70,65,90,95,110,115,120,140,155,150};
     double xvals[10]={80,100,120,140,160,180,200,220,240,260};

gsl_vector * initialconfiguration; 

double xavg = 170.0;
double yavg = 111.0;
double E1(void *xp){
       double beta0 =  gsl_vector_get(xp,0);
       double beta1 =  gsl_vector_get(xp,1);
double yscore=0.0;
	int j=0;
    for (j=0; j<10; j++){
   yscore = pow((yvals[j]-(beta0)-(beta1*xvals[j])),2.0)+yscore;
}
return yscore;
    
}
     
     double M1(void *xp, void *yp)
     {
       double x0 = gsl_vector_get(xp,0);
       double x1 = gsl_vector_get(xp,1);
       double y0 = gsl_vector_get(yp,0);
       double y1 = gsl_vector_get(yp,1);
       
     
       return sqrt(pow((y0-x0),2.0) + pow((y1-x1),2.0));
     }
     
     void S1(const gsl_rng * r, void *xp, double step_size)
     {
       double old_x0 = gsl_vector_get(xp,0);
       double new_x0;
       double u = gsl_rng_uniform(r);
       new_x0 = u * 2 * step_size - step_size + old_x0;
       u = gsl_rng_uniform(r);
       double old_x1 =  gsl_vector_get(xp,1);
       double new_x1;
       new_x1 = u * 2 * step_size - step_size + old_x1;
     
       /*memcpy(gsl_vector_ptr(xp,0), &new_x0, sizeof(new_x0));
       memcpy(gsl_vector_ptr(xp,1), &new_x1, sizeof(new_x1));*/
       gsl_vector_set(xp,0,new_x0);
       gsl_vector_set(xp,1,new_x1);
     }
     
     void P1(void *xp)
     {
       printf ("%12g %12g", gsl_vector_get(xp,0), gsl_vector_get(xp,1));
     }
     
     int
     main(int argc, char *argv[])
     {
       
     initialconfiguration  = gsl_vector_alloc(2);

       const gsl_rng_type * T;
       gsl_rng * r;
     
       double x_initial = 25.0;
       gsl_vector_set_all(initialconfiguration, 25.0);
       gsl_rng_env_setup();
     
       T = gsl_rng_default;
       r = gsl_rng_alloc(T);
     
       gsl_siman_solve(r, &initialconfiguration, E1, S1, M1, P1,
                       NULL, NULL, NULL, 
                       sizeof(initialconfiguration), params);
     
       gsl_rng_free (r);
       return 0;

     }