memory error

This is a discussion on memory error within the C++ Programming forums, part of the General Programming Boards category; hi everyone I have a code for maximum likelihood using gsl_multimin . i must perform the maximization several time in ...

  1. #1
    Registered User
    Join Date
    Aug 2009
    Posts
    8

    Unhappy memory error

    hi everyone
    I have a code for maximum likelihood using gsl_multimin . i must perform the maximization several time in a loop but i keep getting a malloc error.
    this is the code :
    Code:
    /*
     *  likelihood.c
     *  
     *
     *  Created by jeroen Rombouts on 25/08/09.
     *  Copyright 2009 HEC Montreal. All rights reserved.
     *
     */
    #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_vector_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_vector_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(vparamML,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,vsigmastar_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(vllikelihood,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_get(vytemp,i,j));
    		
    		for(i=0;i<ct;i++)
    			gsl_matrix_set(parameters.vy,i,j,gsl_matrix_get(vytemp,(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;
    }
    Last edited by Salem; 09-02-2009 at 10:55 AM. Reason: Fixed code tags, which magically fixes the indent as well

  2. #2
    a_capitalist_story
    Join Date
    Dec 2007
    Posts
    2,655
    Someone failed to read the << !! Posting Code? Read this First !! >> post at the top of the forum, as well as the Forum Guidelines.

  3. #3
    Malum in se abachler's Avatar
    Join Date
    Apr 2007
    Posts
    3,189
    [puke]

    please use code tags and indenting.
    Until you can build a working general purpose reprogrammable computer out of basic components from radio shack, you are not fit to call yourself a programmer in my presence. This is cwhizard, signing off.

  4. #4
    Registered User
    Join Date
    Aug 2009
    Posts
    8
    sorry!!!

    <<DELETED UNINDENTED CODE>>
    FYI, past code from your IDE, not copy/paste from your previous post if you want indent.
    Last edited by Salem; 09-02-2009 at 10:57 AM. Reason: Removed bad code - see original post for good indent

  5. #5
    Registered User
    Join Date
    Dec 2006
    Location
    Canada
    Posts
    3,183
    Quote Originally Posted by abachler View Post
    [puke]

    please use code tags and indenting.
    (emphasis mine)

  6. #6
    Registered User
    Join Date
    Aug 2009
    Posts
    8
    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??

  7. #7
    Registered User
    Join Date
    Dec 2006
    Location
    Canada
    Posts
    3,183
    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?

  8. #8
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,185
    Quote Originally Posted by joeata View Post
    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??
    Quote Originally Posted by cyberfish View Post
    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.

  9. #9
    Registered User
    Join Date
    Aug 2009
    Posts
    8
    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;
    }
    Attached Files Attached Files

  10. #10
    Malum in se abachler's Avatar
    Join Date
    Apr 2007
    Posts
    3,189
    what error is it giving you?
    Until you can build a working general purpose reprogrammable computer out of basic components from radio shack, you are not fit to call yourself a programmer in my presence. This is cwhizard, signing off.

  11. #11
    Registered User
    Join Date
    Sep 2004
    Location
    California
    Posts
    3,254
    What is the text of the error you are getting?
    bit∙hub [bit-huhb] n. A source and destination for information.

  12. #12
    Registered User
    Join Date
    Aug 2009
    Posts
    8
    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

  13. #13
    Registered User
    Join Date
    Sep 2004
    Location
    California
    Posts
    3,254
    It sounds like you are running out of memory. How big is the input file you are using?
    bit∙hub [bit-huhb] n. A source and destination for information.

  14. #14
    Registered User
    Join Date
    Aug 2009
    Posts
    8
    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

  15. #15
    a_capitalist_story
    Join Date
    Dec 2007
    Posts
    2,655
    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 = &par;
    actually need to be done each time through the loop in main()?

Page 1 of 2 12 LastLast
Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Quantum Random Bit Generator
    By shawnt in forum C++ Programming
    Replies: 62
    Last Post: 06-18-2008, 11:17 AM
  2. We Got _DEBUG Errors
    By Tonto in forum Windows Programming
    Replies: 5
    Last Post: 12-22-2006, 05:45 PM
  3. Dikumud
    By maxorator in forum C++ Programming
    Replies: 1
    Last Post: 10-01-2005, 07:39 AM
  4. Please Help - Problem with Compilers
    By toonlover in forum C++ Programming
    Replies: 5
    Last Post: 07-23-2005, 11:03 AM
  5. ras.h errors
    By Trent_Easton in forum Windows Programming
    Replies: 8
    Last Post: 07-15-2005, 11:52 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21