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.
Now I will run through the code. The headers included are
/*--------Include Files-----------------------------------------------------------------*/
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <math.h>
The relevant globals are
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;
where a "component" is a Gaussian. Inside main(int arg, char* arg[]) two gsl_vectors are initialized
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
and memory is allocated for the gsl_vector pointer arrays
vector_ptr_array_current = (gsl_vector **)malloc(sizeof(gsl_vector)*1);
vector_ptr_array_proposed = (gsl_vector **)malloc(sizeof(gsl_vector)*1);
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
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<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;
}
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
no_comp_proposed = no_comp_current+1;
vector_ptr_array_proposed = realloc(vector_ptr_array_proposed, sizeof(gsl_vector)*(no_comp_proposed));
Then call random_birth_move( arg )
random_birth_move(r,vector_ptr_array_current,vector_ptr_array_proposed);
where r and a gsl_rng pointer for random number generation. Below is the random_birth_move function
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<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<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;
}
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:
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<no_comp_current;i++) {
gsl_vector_memcpy(vector_ptr_array_current[i],vector_ptr_array_proposed[i]);
}
return;
}
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 (&vector_ptr_array_current[0],1) to print_components it works, but if I pass the next pointer, again a Segmentation fault: 11 occurs ( (&vector_ptr_array_current[1],1) that is). I appreciate your help.
Edit:
This is the code as it stands
/*--------Include Files-----------------------------------------------------------------*/
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>
#include <math.h>
/*--------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<itr;i++) {
// Birth/Death move or normal?
v = gsl_rng_uniform(r);
if (v < 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<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 < 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<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<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<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<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<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<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(&vector_ptr_array_proposed[0],1);
/*
for (i=0;i<no_comp_current;i++) {
printf("proposed: ");
print_components(&vector_ptr_array_proposed[i],1);
gsl_vector_memcpy(vector_ptr_array_current[i],vector_ptr_array_proposed[i]);
printf("current: ");
print_components(&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<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 < amp_lower || amplitude > amp_upper) {
return -1.0e16;
} else if (mean < avg_lower || mean > avg_upper) {
return -1.0e16;
} else if (sigma < sigma_lower || sigma > 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<no_rows;i++) {
t = gsl_matrix_get(data_matrix,i,0); // extract current time
z_t = 0.0;
for (j=0;j<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<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;
}