gsl_vector_memcpy and dynamic memory allocation
<p>I have an array of pointers to gsl_vectors which is to constantly change size. Each vector stores an ID number and parameters for a mathematical model; specifically a Gaussian i.e. amplitude, mean, variance. Two arrays are in use, vector_ptr_array_current and vector_ptr_array_proposed which store pointers to gsl_vectors. Throughout a for loop the number of gsl_vectors used in the model, each representing an individual Gaussian, changes as the code is trying to determine how many are needed to best fit the data. The array vector_ptr_array_current stores the vectors currently being considered while vector_ptr_array_proposed stores vectors under consideration. If the ones (number of vectors in proposed array may be of a different length) under consideration are accepted I want the code to copy and store them into vector_ptr_array_current. As of right now the code I post below stores the first pointer form vector_ptr_array_propsed into vector_ptr_array_current successfully through the use of gsl_vector_memcpy(gsl_vector * dest, gsl_vector * src) in a for loop, but fails to do so for a second Gaussian. The gsl_vector_memcpy command successfully works, but if I try to print the vector contents of the second Gaussian I get a Segmentation fault: 11.</p>

<p>Now I will run through the code. The headers included are</p>

<pre><code>/*--------Include Files-----------------------------------------------------------------*/
#include &lt;stdio.h&gt;
#include &lt;gsl/gsl_rng.h&gt;
#include &lt;gsl/gsl_randist.h&gt;
#include &lt;gsl/gsl_math.h&gt;
#include &lt;math.h&gt;
</code></pre>

<p>The relevant globals are </p>

<pre><code>int no_comp_current;             // counter for number of components in use currently
int no_comp_proposed;            // number of components proposed
gsl_vector **vector_ptr_array_current;
gsl_vector **vector_ptr_array_proposed;
</code></pre>

<p>where a "component" is a Gaussian. Inside main(int arg, char* arg[]) two gsl_vectors are initialized </p>

<pre><code>gsl_vector *x,*y;            // current state gsl_vector, and proposed
x = gsl_vector_alloc(4);     // allocate memory for gsl_vector
y = gsl_vector_alloc(4);     // for proposal components
</code></pre>

<p>and memory is allocated for the gsl_vector pointer arrays</p>

<pre><code>vector_ptr_array_current = (gsl_vector **)malloc(sizeof(gsl_vector)*1);
vector_ptr_array_proposed = (gsl_vector **)malloc(sizeof(gsl_vector)*1);
</code></pre>

<p>no_comp_current and no_comp_proposed are set to one and a gsl_vector is initialized and stored inside vector_ptr_array_current. The function which prints the components stored in a pointer array is</p>

<pre><code>void print_components(gsl_vector *vector_ptr_array[], int length) {
    int i; // loop counter
    gsl_vector *current_comp_ptr; // gsl_vector ptr for current component under consideration

    // variable initialization
    int id;
    double amp,mean,sigma;

    for (i=0;i&lt;length;i++) {
        current_comp_ptr = vector_ptr_array[i];
        id = gsl_vector_get(current_comp_ptr,0); // get the id from current vector
        amp = gsl_vector_get(current_comp_ptr,1);
        mean = gsl_vector_get(current_comp_ptr,2); 
        sigma = gsl_vector_get(current_comp_ptr,3);  
        printf("ID: %i, amplitude: %f, mean: %f, sigma %f \n",id,amp,mean,sigma);
    }

    return;
}
</code></pre>

<p>which is what ultimately is running when the Segmentation fault: 11 is returned. Now for the sake of testing while developing I force the code to generate a new proposal Gaussian components with a function random_birth_move. First I updated the number of proposed components and reallocate memory for vector_ptr_array_proposed</p>

<pre><code>no_comp_proposed = no_comp_current+1;
vector_ptr_array_proposed = realloc(vector_ptr_array_proposed, sizeof(gsl_vector)*(no_comp_proposed));
</code></pre>

<p>Then call random_birth_move( arg )</p>

<pre><code>random_birth_move(r,vector_ptr_array_current,vector_ptr_array_proposed);
</code></pre>

<p>where r and a gsl_rng pointer for random number generation. Below is the random_birth_move function</p>

<pre><code>void random_birth_move(gsl_rng *r, gsl_vector *vector_ptr_array_current[], gsl_vector *vector_ptr_array_proposed[]) {
    int i; // loop counter
    // initialize and allocate memory for a new component vector
    gsl_vector *x;           
    x = gsl_vector_alloc(4);

    gsl_vector_set(x,0,no_comp_proposed); // set the vector ID

    // set the model parameters of the proposed vector
    for (i=1;i&lt;NO_PARAMS;i++) {
        // randomly set each parameter ~ U[0.0,5.0]
        gsl_vector_set(x,i,gsl_ran_flat(r,0.0,5.0)); 
    }

    // store the old vectors in the proposed pointer array
    for (i=0;i&lt;no_comp_current;i++) {
        // copy current pointers to proposed array (destination, source) 
         gsl_vector_memcpy(vector_ptr_array_proposed[i],vector_ptr_array_current[i]);
        //vector_ptr_array_proposed[i] = vector_ptr_array_current[i];
    }
    // store the new vector component
    vector_ptr_array_proposed[no_comp_current] = x;
    return;
}
</code></pre>

<p>Now if I call the print_components function to print vector_ptr_array_proposed all is well. Now I just assume the proposal is accepted at this moment in development. This acceptance is done by the following function:</p>

<pre><code>void accept_move_prop_to_curr(gsl_vector *vector_ptr_array_proposed[]) {
    int i; // loop counter
    no_comp_current = no_comp_proposed; // update current number of components

    // reallocate current vector ptr array to size of new one
    vector_ptr_array_current = realloc(vector_ptr_array_current,(sizeof(gsl_vector)*(no_comp_proposed)));


    for (i=0;i&lt;no_comp_current;i++) {
         gsl_vector_memcpy(vector_ptr_array_current[i],vector_ptr_array_proposed[i]);
    }

    return;
}
</code></pre>

<p>which is where I believe the problem lives. This function can be called and performed without errors. The seg fault error comes when I print during testing. Printing vector_ptr_array_proposed still results in the correct answer without issues. If I pass (&amp;vector_ptr_array_current[0],1) to print_components it works, but if I pass the next pointer, again a Segmentation fault: 11 occurs ( (&amp;vector_ptr_array_current[1],1) that is). I appreciate your help.</p>

<p>Edit:</p>

<p>This is the code as it stands</p>

<pre><code>/*--------Include Files-----------------------------------------------------------------*/
#include &lt;stdio.h&gt;
#include &lt;gsl/gsl_rng.h&gt;
#include &lt;gsl/gsl_randist.h&gt;
#include &lt;gsl/gsl_math.h&gt;
#include &lt;math.h&gt;

/*--------Defines-----------------------------------------------------------------------*/
#define NO_PARAMS 4            // This is the number of parameters in the model
#define PROB_BIRTH_DEATH 0.0   // probability of choosing a birth/death move
#define PROPOSAL_SIGMA   0.1   // sigma for proposal density

// upper and lower limits on uniform priors for gaussian component parameters
#define amp_lower 0.0
#define amp_upper 10.0
#define avg_lower -10.0
#define avg_upper 10.0
#define sigma_lower 0.00000001 // to prevent any division by zero in model evaluations
#define sigma_upper 10.0


/*--------Globals-----------------------------------------------------------------------*/
int no_rows, no_cols;            // number of rows and columns in data file
int no_comp_current;             // counter for number of components in use currently
int no_comp_proposed;            // number of components proposed
gsl_vector **vector_ptr_array_current;
gsl_vector **vector_ptr_array_proposed;

double amp_prior;
double avg_prior;
double sigma_prior;


/*--------Function Prototypes-----------------------------------------------------------*/

// Initialize parameters for initial gsl_vector before MCMC commences.
void initialize_parameters(gsl_vector *x, gsl_rng *r);  


void random_birth_move(gsl_rng *r, gsl_vector *vector_ptr_array_current[], gsl_vector *vector_ptr_array_proposed[]);
void random_death_move(gsl_rng *r, gsl_vector *vector_ptr_array_current[], gsl_vector *vector_ptr_array_proposed[]);

// If a move is accepted we need to reallocate vector_ptr_array_current to hold the
// stuff in vector_ptr_array_proposed
void accept_move_prop_to_curr(gsl_vector *vector_ptr_array_proposed[]);


// print the components' ID and parameters  (mostly for testing)
void print_components(gsl_vector *vector_ptr_array[], int length);


// Calculates the log of the prior probability density function (pdf)
double log_prior_pdf(gsl_vector *vector_ptr_array[], int no_comps);                   

//  Calculates the log likelihood given the current and proposed parameters.
double log_lklhood(gsl_vector *vector_ptr_array[], int no_comp, gsl_matrix *data_matrix);

// evaluate the model 
double evaluate_model(gsl_vector *x, double t);       

// propose a random variate (rvs) from the proposal density
void propose_parameters(gsl_vector *x, gsl_vector *y, gsl_rng *r);    



/*----------------------------------------------------------------------------------------
 *----------------------------------------------------------------------------------------
 *----------------------------------------------------------------------------------------
 *                                      Main Program
 *----------------------------------------------------------------------------------------
 *----------------------------------------------------------------------------------------
 *----------------------------------------------------------------------------------------
 */

 int main(int argc, char* argv[]) {

    printf("---------------------------------------------------------------------------");
    printf("\n\n");



    /*
     *------------------------------------------------------------------------------------
     *--------Initialization--------------------------------------------------------------
     *------------------------------------------------------------------------------------
     */
    //////////
    int i,j,itr;                  // MCMC iteration counter
    char ch;                  // for counting of rows and columns lated of input file
    char filename[100];       // for in_file and out_file specification
    FILE *in_file;            // input data file
    FILE *out_file;           // output sample file

    double amp,mean,sigma;    // variables to describe a gaussian component
    double u,v;               // random variates ~U[0,1] for making choices
    int acceptance_id;        // identifier for acceptance or rejections.

    double lg_lklhood_y;      // log likelihood evaluated for proposed parameters
    double lg_lklhood_x;      // log likelihood evaluated for current parameters
    double lg_prior_x;        // log prior evaluated at current parameters
    double lg_prior_y;        // log prior evaluated at proposed parameters

    amp_prior = 1.0/(amp_upper-amp_lower); // for prior calculations later
    avg_prior = 1.0/(avg_upper-avg_lower);
    sigma_prior = 1.0/(sigma_upper-sigma_lower);


    double MH_ratio_num;      // metropolis hastings ratio numerator
    double MH_ratio_denom;    // the denominator
    double accept_prob;       // acceptance probability for proposed parameters


    gsl_matrix *data_matrix;  // matrix to read in stored data

    //////////


    //////////
    // Initialize random generators
    const gsl_rng_type * T;
    gsl_rng * r;

    // Allocate generate and specify type
    r = gsl_rng_alloc(gsl_rng_mt19937); 
    gsl_rng_set(r,1);       // seed is set to 1 
    //////////


    //////////
    gsl_vector *x,*y;            // current state gsl_vector, and proposed
    x = gsl_vector_alloc(4);     // allocate memory for gsl_vector
    y = gsl_vector_alloc(4);     // for proposal components

    // gsl_vector *vector_ptr_array_current[1]; // array of pointers to gsl_vector structs
    // gsl_vector ** //(in front of line below)
    vector_ptr_array_current = (gsl_vector **)malloc(sizeof(gsl_vector)*1);
    vector_ptr_array_proposed = (gsl_vector **)malloc(sizeof(gsl_vector)*1);

    // Initialize first component ID and parameters
    initialize_parameters(x,r); // Starting with one component for now
    no_comp_current = 1;
    no_comp_proposed = 1;

    vector_ptr_array_current[0] = x; // store gsl_vector pointer in an array which will change
                                     // size as code runs

    //////////


    //////// open input file and generate data matrix
    sprintf(filename,"%s",argv[1]); // convert input file to string
    in_file = fopen(filename,"rb");  // open input file for reading
    if (in_file == NULL) {
        printf("ERROR opening input file\n\n");
        exit(1);
    }

    // how many rows and columns
    no_rows = 0;
    no_cols = 0;

    while(!feof(in_file)) { // count rows
        ch = fgetc(in_file);
        if(ch == '\n'){
            no_rows++;
        }
    }
    rewind(in_file);
    ch = 'a';
    while(ch != '\n') { // count columns
        ch = fgetc(in_file);
        if (ch == ' ') {
            no_cols++;
        }
    }

    no_cols +=1;               // convert the M could to number of columns
    rewind(in_file);           // so that it can be scanned to generate a matrix
    printf("Input data file\n        no_rows: %i\n        no_cols: %i\n\n",no_rows,no_cols);

    data_matrix = gsl_matrix_alloc(no_rows,no_cols); // allocate memory for data matrix
    gsl_matrix_fscanf(in_file,data_matrix);          // scan in file to matrix
    fclose(in_file);
    ////////

    // set number of iterations from input
    itr = atoi(argv[2]);

    // create output file
    sprintf(filename,"%s",argv[3]); // convert input file to string
    out_file = fopen(filename,"wb");  // open input file for reading
    printf("Output file created: %s\n",filename);

    /*
     *------------------------------------------------------------------------------------
     *--------Initialization Completed----------------------------------------------------
     *------------------------------------------------------------------------------------
     */


    // BELOW IS SOME TESTING OF VARIOUS FUNCTIONS

    printf("\nRandom Birth move:\n");

    // Proposing a random birth move
    print_components(vector_ptr_array_current,no_comp_current);
    printf("\n");

    no_comp_proposed = no_comp_current+1;
    vector_ptr_array_proposed = realloc(vector_ptr_array_proposed, sizeof(gsl_vector)*(no_comp_proposed));

    random_birth_move(r,vector_ptr_array_current,vector_ptr_array_proposed);

    // lets assume previous stuff was accepted
    accept_move_prop_to_curr(vector_ptr_array_proposed);

    print_components(vector_ptr_array_current,no_comp_current);
    printf("You have made it!!!!!\n");
    /*
    printf("\n\nRandom Death move:\n");

    // propose a random death move
    no_comp_proposed = no_comp_current-1;
    vector_ptr_array_proposed = realloc(vector_ptr_array_proposed,(sizeof(gsl_vector)*(no_comp_proposed)));
    random_death_move(r,vector_ptr_array_current,vector_ptr_array_proposed);
    print_components(vector_ptr_array_proposed, no_comp_proposed);

    accept_move_prop_to_curr(vector_ptr_array_proposed);
    */


    /*
     *------------------------------------------------------------------------------------
     *--------MCMC Initiated--------------------------------------------------------------
     *------------------------------------------------------------------------------------
     */
    printf("\nMCMC Initialized\n");

    lg_prior_x = log_prior_pdf(vector_ptr_array_current, no_comp_current);  // calculate the log prior
    lg_lklhood_x = log_lklhood(vector_ptr_array_current, no_comp_current, data_matrix); // calculate the log likelihood

    int model_idx;                // index to specify which model is to be updated
    gsl_vector *current_comp_ptr; // currently viewed gsl_vector


    for (i=0;i&lt;itr;i++) {
        // Birth/Death move or normal?
        v = gsl_rng_uniform(r); 
        if (v &lt; PROB_BIRTH_DEATH) { // perform birth/death move
            continue;
        } else { 
            // randomly choose one component to potentially update
            model_idx = gsl_rng_uniform_int(r,no_comp_current);
            current_comp_ptr = vector_ptr_array_current[model_idx];

            // propose a new set of parameters
            propose_parameters(current_comp_ptr,y,r);

            // copy all of the current ptrs into this proposed ptr array
            vector_ptr_array_proposed = (gsl_vector **)malloc(sizeof(gsl_vector)*no_comp_current);
            for (j=0;j&lt;no_comp_current;j++) {
                if (j==model_idx) {
                    vector_ptr_array_proposed[model_idx] = y;
                } else {
                    gsl_vector_memcpy(vector_ptr_array_proposed[j],vector_ptr_array_current[j]);    
                }
            }

            // calculate the log likelihood of proposed parameters
            lg_lklhood_y = log_lklhood(vector_ptr_array_proposed, no_comp_current, data_matrix);
            // calculate priors on the proposed parameters
            lg_prior_y = log_prior_pdf(vector_ptr_array_proposed, no_comp_current);
            // calculate the acceptance probability
            MH_ratio_num = lg_prior_y + lg_lklhood_y;
            MH_ratio_denom = lg_prior_x + lg_lklhood_x;


            accept_prob = GSL_MIN(1.0, exp(MH_ratio_num-MH_ratio_denom));

            // Decide whether to accept or not
            u = gsl_rng_uniform(r);

            if (u &lt; accept_prob) { // accept the proposed moved
                // store the likelihoods and priors
                lg_lklhood_x = lg_lklhood_y;
                lg_prior_x = lg_prior_y;

                // store the proposed ptr array into the current one
                no_comp_proposed = no_comp_current; // prevents some bug, a bandaid atm
                accept_move_prop_to_curr(vector_ptr_array_proposed);
                acceptance_id = 1;


            } else {
                acceptance_id = 0;
            }


            // gotta print stuff to output file

            // print the iteration, log likelihood, number of components
            fprintf(out_file,"%i %f %i ", i, lg_lklhood_x, no_comp_current);

            // print component parameters
            for (j=0;j&lt;no_comp_current;j++) {
                current_comp_ptr = vector_ptr_array_current[j];
                amp = gsl_vector_get(current_comp_ptr,1);
                mean = gsl_vector_get(current_comp_ptr,2);
                sigma = gsl_vector_get(current_comp_ptr,3);
                fprintf(out_file,"%f %f %f ", amp, mean, sigma);
            }
            // print acceptance id
            fprintf(out_file,"%i\n", acceptance_id);

        }


    }

    /*
     *------------------------------------------------------------------------------------
     *--------MCMC Completed--------------------------------------------------------------
     *------------------------------------------------------------------------------------
     */

    // close output file
    fclose(out_file);

    // Free memory
    gsl_vector_free(x);
    gsl_vector_free(y);
    gsl_rng_free(r);

    printf("MCMC Completed\n");
    printf("\n\n");
    printf("---------------------------------------------------------------------------");
    printf("\n");

    return 0;
}

/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/

void initialize_parameters(gsl_vector *x, gsl_rng *r) {
    int i; // loop iterator

    gsl_vector_set(x,0,1); // set an identification number starting at 1

    for (i=1;i&lt;NO_PARAMS;i++) {
        // randomly set each parameter ~ U[0.0,5.0]
        gsl_vector_set(x,i,gsl_ran_flat(r,0.0,5.0)); 
    }

    return;
}

/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/

void random_birth_move(gsl_rng *r, gsl_vector *vector_ptr_array_current[], gsl_vector *vector_ptr_array_proposed[]) {
    int i; // loop counter
    // initialize and allocate memory for a new component vector
    gsl_vector *x;           
    x = gsl_vector_alloc(4);

    gsl_vector_set(x,0,no_comp_proposed); // set the vector ID

    // set the model parameters of the proposed vector
    for (i=1;i&lt;NO_PARAMS;i++) {
        // randomly set each parameter ~ U[0.0,5.0]
        gsl_vector_set(x,i,gsl_ran_flat(r,0.0,5.0)); 
    }

    // store the old vectors in the proposed pointer array
    for (i=0;i&lt;no_comp_current;i++) {
        // copy current pointers to proposed array (destination, source) 
        // gsl_vector_memcpy(vector_ptr_array_proposed[i],vector_ptr_array_current[i]);
        vector_ptr_array_proposed[i] = vector_ptr_array_current[i];
    }
    // store the new vector component
    vector_ptr_array_proposed[no_comp_current] = x;

    return;
}

/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/

void print_components(gsl_vector *vector_ptr_array[], int length) {
    int i; // loop counter
    gsl_vector *current_comp_ptr; // gsl_vector ptr for current component under consideration

    // variable initialization
    int id;
    double amp,mean,sigma;

    for (i=0;i&lt;length;i++) {
        current_comp_ptr = vector_ptr_array[i];
        id = gsl_vector_get(current_comp_ptr,0); // get the id from current vector
        amp = gsl_vector_get(current_comp_ptr,1);
        mean = gsl_vector_get(current_comp_ptr,2); 
        sigma = gsl_vector_get(current_comp_ptr,3);  
        printf("ID: %i, amplitude: %f, mean: %f, sigma %f \n",id,amp,mean,sigma);
    }

    return;
}

/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/

void random_death_move(gsl_rng *r, gsl_vector *vector_ptr_array_current[], gsl_vector *vector_ptr_array_proposed[]) {
    int i; // loop counter
    int j=0; // counter of index store in propose array

    gsl_vector *current_comp_ptr; // currently viewed gsl_vector

    // randomly select an ID to be proposed for death
    long unsigned int id;
    int id_current;
    id = gsl_rng_uniform_int(r,no_comp_current)+1;

    for (i=0;i&lt;no_comp_current;i++) {
        current_comp_ptr = vector_ptr_array_current[i];
        id_current = (int)gsl_vector_get(current_comp_ptr, 0); // get the current ID value
        if (id_current != id) { // then store it in proposed array of ptrs
            gsl_vector_set(current_comp_ptr,0,j+1); // reassign ID number
            vector_ptr_array_proposed[j] = current_comp_ptr;
            j++;
        }
    }

    return;
}

/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/

void accept_move_prop_to_curr(gsl_vector *vector_ptr_array_proposed[]) {
    int i; // loop counter
    no_comp_current = no_comp_proposed; // update current number of components

    // reallocate current vector ptr array to size of new one
    vector_ptr_array_current = realloc(vector_ptr_array_current,(sizeof(gsl_vector)*(no_comp_proposed)));

    printf("\n\n");
    //gsl_vector_memcpy(vector_ptr_array_current[0],vector_ptr_array_proposed[0]);
    //printf("%i\n",no_comp_proposed);
    memcpy(vector_ptr_array_current[0],vector_ptr_array_proposed[0],sizeof(vector_ptr_array_proposed[0]));
    print_components(&amp;vector_ptr_array_proposed[0],1);
    /*
    for (i=0;i&lt;no_comp_current;i++) {
        printf("proposed: ");
        print_components(&amp;vector_ptr_array_proposed[i],1);
        gsl_vector_memcpy(vector_ptr_array_current[i],vector_ptr_array_proposed[i]);
        printf("current: ");
        print_components(&amp;vector_ptr_array_current[i],1);
        printf("\n");
    }
    */
    return;
}

/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/

double log_prior_pdf(gsl_vector *vector_ptr_array[], int no_comp) {
    double amplitude,mean,sigma;     // to store component parameter values
    double val = 0.0;                // store current value
    int i;                           // loop counter


    for (i=0;i&lt;no_comp;i++){
        amplitude = gsl_vector_get(vector_ptr_array[i],1);
        mean = gsl_vector_get(vector_ptr_array[i],2);
        sigma = gsl_vector_get(vector_ptr_array[i],3);

        if (amplitude &lt; amp_lower || amplitude &gt; amp_upper) { 
            return -1.0e16;
        } else if (mean &lt; avg_lower || mean &gt; avg_upper) {
            return -1.0e16;
        } else if (sigma &lt; sigma_lower || sigma &gt; sigma_upper) {
            return -1.0e16;
        } else {
            val += log(amp_prior*avg_prior*sigma_prior);
        }
    }
    return val;
}

/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/

double log_lklhood(gsl_vector *vector_ptr_array[], int no_comp, gsl_matrix *data_matrix) {

    int i,j;
    double t,z_t;
    double chi_square;

    chi_square = 0.0;

    for (i=0;i&lt;no_rows;i++) {
        t = gsl_matrix_get(data_matrix,i,0); // extract current time
        z_t = 0.0;
        for (j=0;j&lt;no_comp;j++) { // add up contribution from each model at this time
            z_t += evaluate_model(vector_ptr_array[j],t);
        }
        chi_square += pow((gsl_matrix_get(data_matrix,i,1) - z_t),2.0);
    }

    //return -chi_square;
    return 1.0; // test prior recovery
}

/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/

double evaluate_model(gsl_vector *x, double t) {

    double arg; // argument of exponential for gaussian model
    double amplitude,mean,sigma;

    // set the parameters
    amplitude = gsl_vector_get(x,1);
    mean = gsl_vector_get(x,2);
    sigma = gsl_vector_get(x,3);

    arg = -pow((t-mean),2.0)/2.0/pow(sigma,2.0);
    return amplitude*exp(arg);
}

/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------------------*/

void propose_parameters(gsl_vector *x, gsl_vector *y, gsl_rng *r) {
    int i;      // loop iterator
    double val; // value to set 

    gsl_vector_set(y,0,gsl_vector_get(x,0));

    for (i=1;i&lt;NO_PARAMS;i++) { // choose rvs from gaussian
        val = gsl_ran_gaussian(r,PROPOSAL_SIGMA)+gsl_vector_get(x,i);
        gsl_vector_set(y,i,val);
    }

    return;
}
</code></pre>