Thread: memory error

Someone failed to read the << !! Posting Code? Read this First !! >> post at the top of the forum, as well as the Forum Guidelines.

[puke]

please use code tags and indenting.

Originally Posted by abachler

[puke]

please use code tags and indenting.

(emphasis mine)

excuse me all, i am new at this and i am currently trying too debbug by myself. if someone would please help without caring too much about editing??

If you aren't even willing to spend the little time to indent your code to make it easier for us to read, why should we spend OUR time helping you?

Originally Posted by joeata

excuse me all, i am new at this and i am currently trying too debbug by myself. if someone would please help without caring too much about editing??

Originally Posted by cyberfish

If you aren't even willing to spend the little time to indent your code to make it easier for us to read, why should we spend OUR time helping you?

And not to put too fine a point on it, how can you even read it without it being indented?

I can do Ctrl-F to find "malloc" but I don't see it in the code. Presumably it comes back from one of the gsl functions? If so, first you need to figure out which one and then figure out why. If you're passing it something incorrect, then you can fix that.

is this better now? i have attached it too

Code:

/*
* likelihood.c
* 
*/
#include <unistd.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_pow_int.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_linalg.h>
#include <ctime>
#include<gsl/gsl_multimin.h>
#include<gsl/gsl_deriv.h>

struct ML_params {
gsl_matrix *vy;
gsl_matrix *vyinit;
}parameters;


gsl_rng * r;

gsl_vector* sumr(gsl_matrix *mX )
{
size_t M=mX -> size1;
size_t N=mX -> size2;
int m, n;
gsl_vector *sum = gsl_vector_calloc(M);
double aux;

for(m=0;m<M;m++)
{
aux=0.0;
for(n=0;n<N;n++)
{
aux +=gsl_matrix_get(mX,m,n);
}
gsl_vector_set(sum,m,aux);
}
return sum;
}


gsl_matrix* compute_mD(const gsl_matrix*my,const gsl_matrix*myinit,const gsl_matrix *mbeta, const gsl_vector *vsigma, int p,int K)
/*

input: data matrix, location and scale parameters and size of data
output:f(y_t|y_t-1,theta)
*/
{
size_t T=my -> size1;
size_t cy=my -> size2;
int i,j;
int P=p;
int k; 
int t,ctemp=T+p-1;
gsl_matrix *mD = gsl_matrix_calloc(T,K);
gsl_matrix *mX = gsl_matrix_calloc(T,P);
gsl_matrix *mytemp = gsl_matrix_calloc(ctemp,cy);
gsl_vector *vx = gsl_vector_calloc(P);
gsl_vector *vbeta = gsl_vector_calloc(P);
double aux,xbeta=0;



for(i=0;i<ctemp;i++)
{
for(j=0;j<cy;j++)
{
if(i<p-1)
gsl_matrix_set(mytemp,i,j,gsl_matrix_get(myinit,i, j));
else
gsl_matrix_set(mytemp,i,j,gsl_matrix_get(my,i-p+1,j));
}
}
for(t=0;t<T;t++)
gsl_matrix_set(mX,t,0,1);
for(t=0;t<T;t++)
{
for(i=1;i<P;i++)
{
gsl_matrix_set(mX,t,i,gsl_matrix_get(mytemp,(P-i-1+t),0));
}
}

for(t=0;t<T;t++)
{
for(k=0;k<K;k++)
{
gsl_matrix_get_row(vx, mX, t);
gsl_matrix_get_col(vbeta, mbeta, k);
xbeta=0.0;
for(i=0;i<P;i++)
xbeta+=gsl_vector_get(vx,i)*gsl_vector_get(vbeta,i );
aux =(exp(-0.5*((gsl_pow_int ((gsl_matrix_get(my,t,0)-xbeta),2)*1.0/gsl_vector_get(vsigma,k)))))*pow((2.0*M_PI*gsl_vec tor_get(vsigma, k)),-0.5);	
gsl_matrix_set(mD, t, k, aux);
}
}


//	deallocate memory

gsl_matrix_free (mX);
gsl_matrix_free (mytemp);
gsl_vector_free (vx);
gsl_vector_free (vbeta);
return mD;	
}
gsl_matrix * compute_mF(const gsl_matrix *vy,
const gsl_matrix *vyinit,
const gsl_matrix *mbeta,
const gsl_vector *vsigma,
const gsl_matrix *mP,
int p, 
int K) 
{	
int t,i,k;
size_t T=vy -> size1;
gsl_matrix *mF = gsl_matrix_calloc(T,K);

gsl_matrix *mD = gsl_matrix_calloc(T,K);

gsl_matrix_set(mF, 0, 0, 1);
gsl_matrix_set(mF, T-1,K-1, 1);

mD=compute_mD(vy,vyinit,mbeta,vsigma,p,K);

gsl_vector *vtemp = gsl_vector_calloc(K);
gsl_vector *vF = gsl_vector_calloc(K);

gsl_vector *vD = gsl_vector_calloc(K);

double sum,aux;

for( t=1;t<T-1;t++)
{
gsl_matrix_get_row (vF, mF, t-1);
gsl_matrix_get_row (vD, mD, t);

for(k=0;k<K;k++)
{
aux=0.00;
for(i=0;i<K;i++)
aux+= gsl_matrix_get(mF,t-1,i) * gsl_matrix_get(mP,i,k);
gsl_vector_set(vtemp,k,aux);
}
gsl_vector_mul (vtemp, vD);//vtemp=vtemp.*mD[t][]

sum =0.0;
for(k=0;k<K;k++)
sum+= gsl_vector_get(vtemp,k);
gsl_vector_scale (vtemp, 1.0/sum);

for(k=0;k<K;k++)
gsl_matrix_set(mF, t, k,gsl_vector_get(vtemp,k));
}


//	deallocate memory

gsl_matrix_free (mD);
gsl_vector_free (vtemp);
gsl_vector_free (vF);
gsl_vector_free (vD);
return mF;
}

gsl_matrix* compute_mF_likelihood(const gsl_matrix *vy,const gsl_matrix *vyinit,const gsl_matrix *mbeta, const gsl_vector *vsigma, const gsl_matrix *mP,int p, int K) 
{
size_t T=vy -> size1;
gsl_matrix *mF=gsl_matrix_calloc(T,K);

mF=compute_mF(vy,vyinit,mbeta, vsigma, mP, p, K);

double aux;
gsl_vector *vtemp=gsl_vector_alloc(K);
gsl_matrix *QF=gsl_matrix_calloc(T,K);
gsl_matrix_set(QF,0,0,1);
gsl_matrix_set(QF,T-1,K-1,1);

for(int t=1;t<T-1;++t)
{
for(int k=0;k<K;k++)
{
aux=0.00;
for(int i=0;i<K;i++)
aux+= gsl_matrix_get(mF,t-1,i) * gsl_matrix_get(mP,i,k);
gsl_vector_set(vtemp,k,aux);
}

for(int k=0;k<K;k++)
gsl_matrix_set(QF, t, k,gsl_vector_get(vtemp,k));
}	
gsl_vector_free(vtemp);
gsl_matrix_free(mF);
return QF;
}

double llikelihood_ML(const gsl_vector *vparamML,const gsl_matrix *vy,const gsl_matrix * vyinit,int p,int K)
{

size_t T=vy -> size1;	
gsl_vector *vbetastar=gsl_vector_alloc(p*K);
for(int i=0;i<p*K;i++)gsl_vector_set(vbetastar,i,gsl_vecto r_get(vparamML,i));;
gsl_vector *vsigmastar_inv=gsl_vector_alloc(K);	
for(int k=0;k<K;++k)
gsl_vector_set(vsigmastar_inv,k,pow(gsl_vector_get (vparamML,k+p*K),-1));
gsl_vector *vprobstar=gsl_vector_alloc(K-1);
for(int i=0;i<K-1;i++)gsl_vector_set(vprobstar,i,gsl_vector_get(vp aramML,i+p*K+K));;


gsl_matrix *mbetastar=gsl_matrix_calloc(p,K);
gsl_vector *vtemp=gsl_vector_calloc(p);
for(int k=0;k<K;++k)
{
for(int i=0;i<p;i++) gsl_vector_set(vtemp,i,gsl_vector_get(vbetastar,i+ k*p));
gsl_matrix_set_col (mbetastar,k,vtemp);
}


gsl_matrix *mP=gsl_matrix_calloc (K,K);
gsl_matrix_set(mP,K-1,K-1,1);
for(int k=0;k<K-1;++k)
{
for(int i=0; i<K-1;++i)
{
if(i==k)
{
gsl_matrix_set(mP,k,i,gsl_vector_get(vprobstar,k)) ;
gsl_matrix_set(mP,k,i+1,1-gsl_vector_get(vprobstar,k));
}
}
}

gsl_matrix *mlik=gsl_matrix_calloc(T,K);
mlik=compute_mD(vy,vyinit,mbetastar, vsigmastar_inv, p, K);


gsl_matrix *mprob=gsl_matrix_calloc(T,K); 
mprob=compute_mF_likelihood(vy,vyinit,mbetastar,vs igmastar_inv, mP, p,K);


gsl_matrix *mtemp=gsl_matrix_calloc(T,K);
for(int t=0;t<T;t++)
{
for(int k=0;k<K;k++)
gsl_matrix_set(mtemp,t,k,gsl_matrix_get(mlik,t,k)* gsl_matrix_get(mprob,t,k));
}

gsl_vector *vllikelihood= gsl_vector_alloc(T);
vllikelihood= sumr(mtemp);
for(int t=0;t<T;t++)
gsl_vector_set(vllikelihood,t,log(gsl_vector_get(v llikelihood,t)));

double llik=0.0;
for(int t=0;t<T;++t)
llik+=gsl_vector_get(vllikelihood,t);

gsl_vector_free(vbetastar);
gsl_vector_free(vsigmastar_inv);
gsl_vector_free(vprobstar);
gsl_matrix_free(mbetastar);
gsl_vector_free(vtemp);
gsl_matrix_free(mP);
gsl_matrix_free(mlik);
gsl_matrix_free(mprob);
gsl_matrix_free(mtemp);
gsl_vector_free(vllikelihood);

return llik;
}


/* Compute Gradient by central differences Maheu*/
/* For use with GSL minimization routines */
void gsl_central_diff (const gsl_multimin_function * f, const gsl_vector * x, gsl_vector * g)
{
size_t i, n = f->n;
double h = GSL_SQRT_DBL_EPSILON;

//	double h = pow(GSL_DBL_EPSILON,1.0/3.0);
double tmp,dx;

gsl_vector * x1 = gsl_vector_alloc (n);

gsl_vector_memcpy (x1, x);

for (i = 0; i < n; i++)
{
double fl, fh;

double xi = gsl_vector_get (x, i);
// double dx = fabs(xi) * h;
tmp = GSL_MAX_DBL ( fabs(xi)*h , h ) + xi;
dx = tmp - xi;

if (dx == 0.0) dx = h;

gsl_vector_set (x1, i, xi + dx);
fh = GSL_MULTIMIN_FN_EVAL(f, x1);
gsl_vector_set (x1, i, xi);

gsl_vector_set (x1, i, xi - dx);
fl = GSL_MULTIMIN_FN_EVAL(f, x1);
gsl_vector_set (x1, i, xi);

gsl_vector_set (g, i, (fh - fl) / (2.0 * dx));
}

gsl_vector_free (x1);

//	return GSL_SUCCESS;
}



double
my_f (const gsl_vector *v, void *params)
{
int *p = (int *)params;
return -(llikelihood_ML(v,parameters.vy,parameters.vyinit, p[0], p[1]));
}

/* The gradient of f, df = (df/dx, df/dy). */
void 
my_df (const gsl_vector *v, void *params, 
gsl_vector *df)
{
gsl_multimin_function F ;
F.f = my_f;
F.params = params;
F.n = v->size;
gsl_central_diff (&F,v,df);
}

/* Compute both f and df together. */
void 
my_fdf (const gsl_vector *x, void *params, 
double *f, gsl_vector *df) 
{
*f = my_f(x, params); 
my_df(x, params, df);
}

int main (void)
{
gsl_rng_env_setup();
//	Taus = gsl_rng_default;
r = gsl_rng_alloc (gsl_rng_taus);
gsl_rng_set (r,(unsigned long int)894526599);


int p =2;
int K=9;
int i, j;

//	lecture des données
int ctemp=735,cy=1;
gsl_matrix * vytemp = gsl_matrix_calloc(ctemp,cy); 
FILE *f = fopen("yext.txt", "rb");
if(!f) { // file couldn't be opened
printf( "Error: file could not be opened" );
}	
gsl_matrix_fscanf (f, vytemp);
fclose (f);

int ct=(ctemp-p+1);

parameters.vyinit = gsl_matrix_calloc(ctemp-ct,cy); 
parameters.vy = gsl_matrix_calloc(ct,cy);
for(j=0;j<cy;j++)
{
for(i=0;i<p-1;i++)
gsl_matrix_set(parameters.vyinit,i,j,gsl_matrix_ge t(vytemp,i,j));
for(i=0;i<ct;i++)
gsl_matrix_set(parameters.vy,i,j,gsl_matrix_get(vy temp,(i+p-1),j));
}

int R=70;
for(int rep=0;rep<R;rep++)
{
size_t iter = 0;
int status;

const gsl_multimin_fdfminimizer_type *T;
gsl_multimin_fdfminimizer *s;

/* Position of the minimum (1,2), scale factors 
10,20, height 30. */
int par[2] = { p, K };

gsl_multimin_function_fdf my_func;

my_func.n = 35;
my_func.f = &my_f;
my_func.df = &my_df;
my_func.fdf = &my_fdf;
my_func.params = &par;

/* Starting point, x = (5,7) */
gsl_vector *m_vPar=gsl_vector_calloc(35);



gsl_vector_set(m_vPar,0,0.021);
gsl_vector_set(m_vPar,1,1.001);
gsl_vector_set(m_vPar,2,0.207);
gsl_vector_set(m_vPar,3,0.908);
gsl_vector_set(m_vPar,4,0.020);
gsl_vector_set(m_vPar,5,1.006);
gsl_vector_set(m_vPar,6,0.201);
gsl_vector_set(m_vPar,7,0.972);
gsl_vector_set(m_vPar,8, 0.215);
gsl_vector_set(m_vPar,9,0.973);
gsl_vector_set(m_vPar,10,0.206);
gsl_vector_set(m_vPar,11,0.973);
gsl_vector_set(m_vPar,12,0.007);
gsl_vector_set(m_vPar,13,0.990);
gsl_vector_set(m_vPar,14,0.167);
gsl_vector_set(m_vPar,15,0.977);
gsl_vector_set(m_vPar,16,0.098);
gsl_vector_set(m_vPar,17, 0.853);
gsl_vector_set(m_vPar,18,44.066);
gsl_vector_set(m_vPar,19,4.193);
gsl_vector_set(m_vPar,20,65.428);
gsl_vector_set(m_vPar,21, 3.918);
gsl_vector_set(m_vPar,22,0.420);
gsl_vector_set(m_vPar,23, 6.440);
gsl_vector_set(m_vPar,24,21.681);
gsl_vector_set(m_vPar,25,58.375);
gsl_vector_set(m_vPar,26,4.695);
gsl_vector_set(m_vPar,27,0.988);
gsl_vector_set(m_vPar,28,0.959);
gsl_vector_set(m_vPar,29,0.980);
gsl_vector_set(m_vPar,30,0.990);
gsl_vector_set(m_vPar,31,0.961);
gsl_vector_set(m_vPar,32, 0.982);
gsl_vector_set(m_vPar,33,0.991);
gsl_vector_set(m_vPar,34,0.968);


T = gsl_multimin_fdfminimizer_vector_bfgs2;
s = gsl_multimin_fdfminimizer_alloc (T, 35);

gsl_multimin_fdfminimizer_set (s, &my_func, m_vPar, 0.005, 0.1);

do
{
iter++;
status = gsl_multimin_fdfminimizer_iterate (s);

if (status)
break;

status = gsl_multimin_test_gradient (s->gradient, 10);

if (status == GSL_SUCCESS)
printf ("Minimum found at:\n");

printf ("%5d %10.5f\n", iter,	s->f);
}
while (status == GSL_CONTINUE && iter <35);


printf("Replication %d\n ", rep+1);
if( iter==35)	
rep--;
gsl_multimin_fdfminimizer_free (s);
gsl_vector_free (m_vPar);
}

return 0;
}

what error is it giving you?

i get this message

c(223) malloc: *** mmap(size=53248) failed (error code=12)
*** error: can't allocate region
*** set a breakpoint in malloc_error_break to debug
gsl: init_source.c:46: ERROR: failed to allocate space for block data
Default GSL error handler invoked.
Abort trap

the size of the input file : 8KB , it is a vector of 734 elements, i tried with another input file which
size is 4KB ( a vector of 250 elements ) but it still does not work

Does this:

Code:

int par[2] = { p, K };	
gsl_multimin_function_fdf my_func;
my_func.n = 35;
my_func.f = &my_f;
my_func.df = &my_df;
my_func.fdf = &my_fdf;
my_func.params = ∥

actually need to be done each time through the loop in main()?