/* $Id$
$State$
$Log$
+ Revision 1.130 2009/05/26 06:44:34 brouard
+ (Module): Max Covariate is now set to 20 instead of 8. A
+ lot of cleaning with variables initialized to 0. Trying to make
+ V2+V3*age+V1+V4 strb=V3*age+V1+V4 working better.
+
Revision 1.129 2007/08/31 13:49:27 lievre
Modification of the way of exiting when the covariate is not binary in order to see on the window the error message before exiting
#define GLOCK_ERROR_NOPATH -1 /* empty path */
#define GLOCK_ERROR_GETCWD -2 /* cannot get cwd */
-#define MAXPARM 30 /* Maximum number of parameters for the optimization */
+#define MAXPARM 128 /* Maximum number of parameters for the optimization */
#define NPARMAX 64 /* (nlstate+ndeath-1)*nlstate*ncovmodel */
#define NINTERVMAX 8
/* $Id$ */
/* $State$ */
-char version[]="Imach version 0.98i, June 2006, INED-EUROREVES-Institut de longevite ";
+char version[]="Imach version 0.98k, June 2006, INED-EUROREVES-Institut de longevite ";
char fullversion[]="$Revision$ $Date$";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
double *weight;
int **s; /* Status */
double *agedc, **covar, idx;
-int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;
+int **nbcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;
double *lsurv, *lpop, *tpop;
double ftol=FTOL; /* Tolerance for computing Max Likelihood */
int i, ii,j,k;
double min, max, maxmin, maxmax,sumnew=0.;
double **matprod2();
- double **out, cov[NCOVMAX], **pmij();
+ double **out, cov[NCOVMAX+1], **pmij();
double **newm;
double agefin, delaymax=50 ; /* Max number of years to converge */
for(i=1; i<= nlstate; i++){
s1=0;
- for(j=1; j<i; j++)
+ for(j=1; j<i; j++){
s1+=exp(ps[i][j]);
- for(j=i+1; j<=nlstate+ndeath; j++)
+ /*printf("debug1 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
+ for(j=i+1; j<=nlstate+ndeath; j++){
s1+=exp(ps[i][j]);
+ /*printf("debug2 %d %d ps=%lf exp(ps)=%lf s1+=%lf\n",i,j,ps[i][j],exp(ps[i][j]),s1); */
+ }
ps[i][i]=1./(s1+1.);
for(j=1; j<i; j++)
ps[i][j]= exp(ps[i][j])*ps[i][i];
*/
int i, j, d, h, k;
- double **out, cov[NCOVMAX];
+ double **out, cov[NCOVMAX+1];
double **newm;
/* Hstepm could be zero and should return the unit matrix */
/* Covariates have to be included here again */
cov[1]=1.;
cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
- for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+ for (k=1; k<=cptcovn;k++)
+ cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
for (k=1; k<=cptcovage;k++)
cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
for (k=1; k<=cptcovprod;k++)
double func( double *x)
{
int i, ii, j, k, mi, d, kk;
- double l, ll[NLSTATEMAX], cov[NCOVMAX];
+ double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
double **out;
double sw; /* Sum of weights */
double lli; /* Individual log likelihood */
{
/* Same as likeli but slower because of a lot of printf and if */
int i, ii, j, k, mi, d, kk;
- double l, ll[NLSTATEMAX], cov[NCOVMAX];
+ double l, ll[NLSTATEMAX+1], cov[NCOVMAX+1];
double **out;
double lli; /* Individual log likelihood */
double llt;
ipmx +=1;
sw += weight[i];
ll[s[mw[mi][i]][i]] += 2*weight[i]*lli;
-/* printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]); */
+ printf("i=%6d s1=%1d s2=%1d mi=%1d mw=%1d dh=%3d prob=%10.6f w=%6.4f out=%10.6f sav=%10.6f\n",i,s1,s2,mi,mw[mi][i],dh[mi][i],exp(lli),weight[i],out[s1][s2],savm[s1][s2]);
if(globpr){
fprintf(ficresilk,"%9d %6d %2d %2d %1d %1d %3d %11.6f %8.4f\
%11.6f %11.6f %11.6f ", \
/* Tvar[i]=atoi(stre); /* find 'n' in Vn and stores in Tvar. If model=V2+V1 Tvar[1]=2 and Tvar[2]=1 */
- int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=19;
+ int Ndum[20],ij=1, k=0, j=0, i=0, maxncov=NCOVMAX;
int cptcode=0;
cptcoveff=0;
for (k=0; k<maxncov; k++) Ndum[k]=0;
- for (k=1; k<=7; k++) ncodemax[k]=0;
+ for (k=1; k<=7; k++) ncodemax[k]=0; /* Horrible constant again */
- for (j=1; j<=(cptcovn+2*cptcovprod); j++) {
+ for (j=1; j<=(cptcovn+2*cptcovprod); j++) { /* For each covariate */
for (i=1; i<=imx; i++) { /*reads the data file to get the maximum
modality*/
- ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual*/
- Ndum[ij]++; /*store the modality */
+ ij=(int)(covar[Tvar[j]][i]); /* ij is the modality of this individual, might be -1*/
+ Ndum[ij]++; /*counts the occurence of this modality */
/*printf("i=%d ij=%d Ndum[ij]=%d imx=%d",i,ij,Ndum[ij],imx);*/
- if (ij > cptcode) cptcode=ij; /* getting the maximum of covariable
+ if (ij > cptcode) cptcode=ij; /* getting the maximum value of the modality of the covariate (should be 0 or 1 now)
Tvar[j]. If V=sex and male is 0 and
female is 1, then cptcode=1.*/
}
- for (i=0; i<=cptcode; i++) {
- if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariates. In fact ncodemax[j]=2 (dichotom. variables) but it can be more */
- }
+ for (i=0; i<=cptcode; i++) { /* i=-1 ?*/
+ if(Ndum[i]!=0) ncodemax[j]++; /* Nomber of modalities of the j th covariate. In fact ncodemax[j]=2 (dichotom. variables only) but it can be more */
+ } /* Ndum[-1] number of undefined modalities */
ij=1;
- for (i=1; i<=ncodemax[j]; i++) {
- for (k=0; k<= maxncov; k++) {
- if (Ndum[k] != 0) {
- nbcode[Tvar[j]][ij]=k;
- /* store the modality in an array. k is a modality. If we have model=V1+V1*sex then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
-
+ for (i=1; i<=ncodemax[j]; i++) { /* i= 1 to 2 */
+ for (k=0; k<= maxncov; k++) { /* k=-1 ?*/
+ if (Ndum[k] != 0) { /* If at least one individual responded to this modality k */
+ nbcode[Tvar[j]][ij]=k; /* stores the modality in an array nbcode.
+ k is a modality. If we have model=V1+V1*sex
+ then: nbcode[1][1]=0 ; nbcode[1][2]=1; nbcode[2][1]=0 ; nbcode[2][2]=1; */
ij++;
}
if (ij > ncodemax[j]) break;
for (k=0; k< maxncov; k++) Ndum[k]=0;
- for (i=1; i<=ncovmodel-2; i++) {
+ for (i=1; i<=ncovmodel-2; i++) { /* -2, cste and age */
/* Listing of all covariables in statement model to see if some covariates appear twice. For example, V1 appears twice in V1+V1*V2.*/
- ij=Tvar[i];
+ ij=Tvar[i]; /* Tvar might be -1 if status was unknown */
Ndum[ij]++;
}
ij++;
}
}
-
- cptcoveff=ij-1; /*Number of simple covariates*/
+ ij--;
+ cptcoveff=ij; /*Number of simple covariates*/
}
/*********** Health Expectancies ****************/
free_vector(gmp,nlstate+1,nlstate+ndeath);
free_matrix(gradgp,1,npar,nlstate+1,nlstate+ndeath);
free_matrix(trgradgp,nlstate+1,nlstate+ndeath,1,npar); /* mu or p point j*/
- fprintf(ficgp,"\nset noparametric;set nolabel; set ter png small;set size 0.65, 0.65");
+ fprintf(ficgp,"\nunset parametric;unset label; set ter png small;set size 0.65, 0.65");
/* for(j=nlstate+1; j<= nlstate+ndeath; j++){ *//* Only the first actually */
- fprintf(ficgp,"\n set log y; set nolog x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
+ fprintf(ficgp,"\n set log y; unset log x;set xlabel \"Age\"; set ylabel \"Force of mortality (year-1)\";");
/* fprintf(ficgp,"\n plot \"%s\" u 1:($3*%6.3f) not w l 1 ",fileresprobmorprev,YEARM/estepm); */
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3+1.96*$4)*%6.3f) t \"95\%% interval\" w l 2 ",fileresprobmorprev,YEARM/estepm); */
/* fprintf(ficgp,"\n replot \"%s\" u 1:(($3-1.96*$4)*%6.3f) not w l 2 ",fileresprobmorprev,YEARM/estepm); */
fprintf(ficresprobcov,"\n");
fprintf(ficresprobcor,"\n");
*/
- xp=vector(1,npar);
+ xp=vector(1,npar);
dnewm=matrix(1,(nlstate)*(nlstate+ndeath),1,npar);
doldm=matrix(1,(nlstate)*(nlstate+ndeath),1,(nlstate)*(nlstate+ndeath));
mu=matrix(1,(nlstate)*(nlstate+ndeath), (int) bage, (int)fage);
/* Confidence intervalle of pij */
/*
- fprintf(ficgp,"\nset noparametric;unset label");
+ fprintf(ficgp,"\nunset parametric;unset label");
fprintf(ficgp,"\nset log y;unset log x; set xlabel \"Age\";set ylabel \"probability (year-1)\"");
fprintf(ficgp,"\nset ter png small\nset size 0.65,0.65");
fprintf(fichtm,"\n<br>Probability with confidence intervals expressed in year<sup>-1</sup> :<a href=\"pijgr%s.png\">pijgr%s.png</A>, ",optionfilefiname,optionfilefiname);
for (i=1; i<= nlstate ; i ++) {
if (i==cpt) fprintf(ficgp," \%%lf (\%%lf)");
- else fprintf(ficgp," \%%*lf (\%%*lf)");
+ else fprintf(ficgp," \%%*lf (\%%*lf)");
}
fprintf(ficgp,"\" t\"Period (stable) prevalence\" w l 0,\"%s\" every :::%d::%d u 1:($2+1.96*$3) \"\%%lf",subdirf2(fileres,"vpl"),k1-1,k1-1);
for (i=1; i<= nlstate ; i ++) {
ncovmodel=2+cptcovn; /*Number of variables = cptcovn + intercept + age */
nvar=ncovmodel-1; /* Suppressing age as a basic covariate */
npar= (nlstate+ndeath-1)*nlstate*ncovmodel; /* Number of parameters*/
-
+ if(npar >MAXPARM || nlstate >NLSTATEMAX || ndeath >NDEATHMAX || ncovmodel>NCOVMAX){
+ printf("Too complex model for current IMaCh: npar=(nlstate+ndeath-1)*nlstate*ncovmodel=%d >= %d(MAXPARM) or nlstate=%d >= %d(NLSTATEMAX) or ndeath=%d >= %d(NDEATHMAX) or ncovmodel=(k+age+#of+signs)=%d(NCOVMAX) >= %d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX);
+ fprintf(ficlog,"Too complex model for current IMaCh: %d >=%d(MAXPARM) or %d >=%d(NLSTATEMAX) or %d >=%d(NDEATHMAX) or %d(NCOVMAX) >=%d\n",npar, MAXPARM, nlstate, NLSTATEMAX, ndeath, NDEATHMAX, ncovmodel, NCOVMAX);
+ fflush(stdout);
+ fclose (ficlog);
+ goto end;
+ }
delti3= ma3x(1,nlstate,1,nlstate+ndeath-1,1,ncovmodel);
delti=delti3[1][1];
/*delti=vector(1,npar); *//* Scale of each paramater (output from hesscov)*/
ungetc(c,ficpar);
matcov=matrix(1,npar,1,npar);
+ for(i=1; i <=npar; i++)
+ for(j=1; j <=npar; j++) matcov[i][j]=0.;
+
for(i=1; i <=npar; i++){
fscanf(ficpar,"%s",&str);
if(mle==1)
for(i=1;i<=n;i++) weight[i]=1.0; /* Equal weights, 1 by default */
mint=matrix(1,maxwav,1,n);
anint=matrix(1,maxwav,1,n);
- s=imatrix(1,maxwav+1,1,n);
+ s=imatrix(1,maxwav+1,1,n); /* s[i][j] health state for wave i and individual j */
tab=ivector(1,NCOVMAX);
ncodemax=ivector(1,8);
for (j=maxwav;j>=1;j--){
cutv(stra, strb,line,' ');
- errno=0;
- lval=strtol(strb,&endptr,10);
+ if(strb[0]=='.') { /* Missing status */
+ lval=-1;
+ }else{
+ errno=0;
+ lval=strtol(strb,&endptr,10);
/* if (errno == ERANGE && (lval == LONG_MAX || lval == LONG_MIN))*/
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %d %s for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);
- exit(1);
+ if( strb[0]=='\0' || (*endptr != '\0')){
+ printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a status of wave %d. Setting maxwav=%d might be wrong. Exiting.\n", strb, linei,i,line,j,maxwav);fflush(ficlog);
+ goto end;
+ }
}
s[j][i]=lval;
month=99;
year=9999;
}else{
- printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);
- exit(1);
+ printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of interview (mm/yyyy or .) at wave %d. Exiting.\n",strb, linei,i, line,j);fflush(ficlog);
+ goto end;
}
anint[j][i]= (double) year;
mint[j][i]= (double)month;
month=99;
year=9999;
}else{
- printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);
- exit(1);
+ printf("Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld for individual %d, '%s'\nShould be a date of death (mm/yyyy or .). Exiting.\n",strb, linei,i,line);fflush(ficlog);
+ goto end;
}
andc[i]=(double) year;
moisdc[i]=(double) month;
year=9999;
}else{
printf("Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j);
- exit(1);
+ fprintf(ficlog,"Error reading data around '%s' at line number %ld %s for individual %d, '%s'\nShould be a date of birth (mm/yyyy or .). Exiting.\n",strb, linei,i,line,j);fflush(ficlog);
+ goto end;
}
annais[i]=(double)(year);
moisnais[i]=(double)(month);
dval=strtod(strb,&endptr);
if( strb[0]=='\0' || (*endptr != '\0')){
printf("Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei);
- exit(1);
+ fprintf(ficlog,"Error reading data around '%f' at line number %ld, \"%s\" for individual %d\nShould be a weight. Exiting.\n",dval, i,line,linei);
+ fflush(ficlog);
+ goto end;
}
weight[i]=dval;
strcpy(line,stra);
for (j=ncovcol;j>=1;j--){
cutv(stra, strb,line,' ');
- errno=0;
- lval=strtol(strb,&endptr,10);
- if( strb[0]=='\0' || (*endptr != '\0')){
- printf("Error reading data around '%d' at line number %ld %s for individual %d, '%s'\nShould be a covar (meaning 0 for the reference or 1). Exiting.\n",lval, linei,i, line);
- exit(1);
+ if(strb[0]=='.') { /* Missing status */
+ lval=-1;
+ }else{
+ errno=0;
+ lval=strtol(strb,&endptr,10);
+ if( strb[0]=='\0' || (*endptr != '\0')){
+ printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);
+ fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\nShould be a covariate value (=0 for the reference or 1 for alternative). Exiting.\n",lval, linei,i, line);fflush(ficlog);
+ goto end;
+ }
}
if(lval <-1 || lval >1){
printf("Error reading data around '%d' at line number %ld for individual %d, '%s'\n \
and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
output of IMaCh is often meaningless.\n \
Exiting.\n",lval,linei, i,line,j);
+ fprintf(ficlog,"Error reading data around '%d' at line number %ld for individual %d, '%s'\n \
+ Should be a value of %d(nth) covariate (0 should be the value for the reference and 1\n \
+ for the alternative. IMaCh does not build design variables automatically, do it yourself.\n \
+ For example, for multinomial values like 1, 2 and 3,\n \
+ build V1=0 V2=0 for the reference value (1),\n \
+ V1=1 V2=0 for (2) \n \
+ and V1=0 V2=1 for (3). V1=1 V2=1 should not exist and the corresponding\n \
+ output of IMaCh is often meaningless.\n \
+ Exiting.\n",lval,linei, i,line,j);fflush(ficlog);
goto end;
}
covar[j][i]=(double)(lval);
else weight[i]=1;*/
/* Calculation of the number of parameters from char model */
- Tvar=ivector(1,15); /* stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
+ Tvar=ivector(1,NCOVMAX); /* Was 15 changed to NCOVMAX. Stores the number n of the covariates in Vm+Vn at 1 and m at 2 */
Tprod=ivector(1,15);
Tvaraff=ivector(1,15);
Tvard=imatrix(1,15,1,2);
strcpy(modelsav,model);
if ((strcmp(model,"age")==0) || (strcmp(model,"age*age")==0)){
printf("Error. Non available option model=%s ",model);
- fprintf(ficlog,"Error. Non available option model=%s ",model);
+ fprintf(ficlog,"Error. Non available option model=%s ",model);fflush(ficlog);
goto end;
}
/* Routine tricode is to calculate cptcoveff (real number of unique covariates) and to associate covariable number and modality */
- Tcode=ivector(1,100);
nbcode=imatrix(0,NCOVMAX,0,NCOVMAX);
ncodemax[1]=1;
if (cptcovn > 0) tricode(Tvar,nbcode,imx);
h=0;
m=pow(2,cptcoveff);
- for(k=1;k<=cptcoveff; k++){
- for(i=1; i <=(m/pow(2,k));i++){
- for(j=1; j <= ncodemax[k]; j++){
- for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
+ for(k=1;k<=cptcoveff; k++){ /* scans any effective covariate */
+ for(i=1; i <=(m/pow(2,k));i++){ /* i=1 to 8/1=8; i=1 to 8/2=4; i=1 to 8/8=1 */
+ for(j=1; j <= ncodemax[k]; j++){ /* For each modality of this covariate */
+ for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){ /* cpt=1 to 8/2**(3+1-1 or 3+1-3) =1 or 4 */
h++;
- if (h>m) h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;
+ if (h>m)
+ h=1;codtab[h][k]=j;codtab[h][Tvar[k]]=j;
printf("h=%d k=%d j=%d codtab[h][k]=%d Tvar[k]=%d codtab[h][Tvar[k]]=%d \n",h, k,j,codtab[h][k],Tvar[k],codtab[h][Tvar[k]]);
}
}
codtab[1][2]=1;codtab[2][2]=2; */
/* for(i=1; i <=m ;i++){
for(k=1; k <=cptcovn; k++){
- printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);
+ printf("i=%d k=%d %d %d ",i,k,codtab[i][k], cptcoveff);
}
printf("\n");
}
strcat(optionfilehtm,"-mort");
strcat(optionfilehtm,".htm");
if((fichtm=fopen(optionfilehtm,"w"))==NULL) {
- printf("Problem with %s \n",optionfilehtm), exit(0);
+ printf("Problem with %s \n",optionfilehtm);
+ exit(0);
}
strcpy(optionfilehtmcov,optionfilefiname); /* Only for matrix of covariance */
} /* Endof if mle==-3 */
else{ /* For mle >=1 */
-
- likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
+ globpr=1;/* debug */
+ /* likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone);*/ /* Prints the contributions to the likelihood */
printf("First Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
for (k=1; k<=npar;k++)
printf(" %d %8.5f",k,p[k]);
for(cptcov=1,k=0;cptcov<=i1;cptcov++){
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;
- /*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
+ /* to clean */
+ printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,codtab[cptcod][cptcov],nbcode);
fprintf(ficrespl,"\n#******");
printf("\n#******");
fprintf(ficlog,"\n#******");
free_ma3x(probs,1,AGESUP,1,NCOVMAX, 1,NCOVMAX);
} /* mle==-3 arrives here for freeing */
+ endfree:
free_matrix(prlim,1,nlstate,1,nlstate);
free_matrix(pmmij,1,nlstate+ndeath,1,nlstate+ndeath);
free_matrix(oldms, 1,nlstate+ndeath,1,nlstate+ndeath);
free_ivector(Tprod,1,15);
free_ivector(Tvaraff,1,15);
free_ivector(Tage,1,15);
- free_ivector(Tcode,1,100);
free_imatrix(nbcode,0,NCOVMAX,0,NCOVMAX);
free_imatrix(codtab,1,100,1,10);