/* $Id$
$State$
$Log$
+ Revision 1.185 2015/03/11 13:26:42 brouard
+ Summary: Inclusion of compile and links command line for Intel Compiler
+
Revision 1.184 2015/03/11 11:52:39 brouard
Summary: Back from Windows 8. Intel Compiler
*/
#define POWELL /* Instead of NLOPT */
-/* #define POWELLORIGINAL */ /* Don't use Directest to decide new direction but original Powell test */
-/* #define MNBRAKORIGINAL */ /* Don't use mnbrak fix */
+/* #define POWELLORIGINAL /\* Don't use Directest to decide new direction but original Powell test *\/ */
+/* #define MNBRAKORIGINAL /\* Don't use mnbrak fix *\/ */
#include <math.h>
#include <stdio.h>
/* $Id$ */
/* $State$ */
-char version[]="Imach version 0.98q0, March 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
+char version[]="Imach version 0.98q1, April 2015,INED-EUROREVES-Institut de longevite-Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research 25293121), Intel Software 2015";
char fullversion[]="$Revision$ $Date$";
char strstart[80];
char optionfilext[10], optionfilefiname[FILENAMELENGTH];
the name of the file (name), its extension only (ext) and its first part of the name (finame)
*/
char *ss; /* pointer */
- int l1, l2; /* length counters */
+ int l1=0, l2=0; /* length counters */
l1 = strlen(path ); /* length of path */
if ( l1 == 0 ) return( GLOCK_ERROR_NOPATH );
if ( l2 == 0 ) return( GLOCK_ERROR_NOPATH );
strcpy( name, ss ); /* save file name */
strncpy( dirc, path, l1 - l2 ); /* now the directory */
- dirc[l1-l2] = 0; /* add zero */
+ dirc[l1-l2] = '\0'; /* add zero */
printf(" DIRC2 = %s \n",dirc);
}
/* We add a separator at the end of dirc if not exists */
cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
/*printf("prevalim ij=%d k=%d Tvar[%d]=%d nbcode=%d cov=%lf codtab[%d][Tvar[%d]]=%d \n",ij,k, k, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k], ij, 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++) /\* Useless *\/ */
- /* cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]]; */
+ /*wrong? for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+ for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]]*cov[2];
+ for (k=1; k<=cptcovprod;k++) /* Useless */
+ cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]] * nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
/*printf("ij=%d cptcovprod=%d tvar=%d ", ij, cptcovprod, Tvar[1]);*/
/*printf("ij=%d cov[3]=%lf cov[4]=%lf \n",ij, cov[3],cov[4]);*/
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<=cptcovage;k++)
- cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
+ for (k=1; k<=cptcovage;k++) /* Should start at cptcovn+1 */
+ /* cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+ cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2];
for (k=1; k<=cptcovprod;k++) /* Useless because included in cptcovn */
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
for (k=1; k <= maxncov; k++) ncodemax[k]=0; /* Horrible constant again replaced by NCOVMAX */
/* Loop on covariates without age and products */
- for (j=1; j<=(cptcovs); j++) { /* model V1 + V2*age+ V3 + V3*V4 : V1 + V3 = 2 only */
- for (i=1; i<=imx; i++) { /* Lopp on individuals: reads the data file to get the maximum value of the
+ for (j=1; j<=(cptcovs); j++) { /* From model V1 + V2*age+ V3 + V3*V4 keeps V1 + V3 = 2 only */
+ for (i=1; i<=imx; i++) { /* Loop on individuals: reads the data file to get the maximum value of the
modality of this covariate Vj*/
ij=(int)(covar[Tvar[j]][i]); /* ij=0 or 1 or -1. Value of the covariate Tvar[j] for individual i
* If product of Vn*Vm, still boolean *:
} /* Ndum[-1] number of undefined modalities */
/* j is a covariate, n=Tvar[j] of Vn; Fills nbcode */
- /* For covariate j, modalities could be 1, 2, 3, 4. If Ndum[2]=0 ncodemax[j] is not 4 but 3 */
- /* If Ndum[3}= 635; Ndum[4]=0; Ndum[5]=0; Ndum[6]=27; Ndum[7]=125;
+ /* For covariate j, modalities could be 1, 2, 3, 4, 5, 6, 7.
+ If Ndum[1]=0, Ndum[2]=0, Ndum[3]= 635, Ndum[4]=0, Ndum[5]=0, Ndum[6]=27, Ndum[7]=125;
modmincovj=3; modmaxcovj = 7;
- There are only 3 modalities non empty (or 2 if 27 is too few) : ncodemax[j]=3;
- which will be coded 0, 1, 2 which in binary on 3-1 digits are 0=00 1=01, 2=10; defining two dummy
- variables V1_1 and V1_2.
+ There are only 3 modalities non empty 3, 6, 7 (or 2 if 27 is too few) : ncodemax[j]=3;
+ which will be coded 0, 1, 2 which in binary on 2=3-1 digits are 0=00 1=01, 2=10;
+ defining two dummy variables: variables V1_1 and V1_2.
nbcode[Tvar[j]][ij]=k;
nbcode[Tvar[j]][1]=0;
nbcode[Tvar[j]][2]=1;
for (k=0; k<= cptcode; k++) { /* k=-1 ? k=0 to 1 *//* Could be 1 to 4 */
/*recode from 0 */
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.
+ nbcode[Tvar[j]][ij]=k; /* stores the modality k 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++;
*/
/* nbcode[1][1]=0 nbcode[1][2]=1;*/
}
- for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2];
+ /* for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=cov[2+Tage[k]]*cov[2]; */
+ for (k=1; k<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2];
for (k=1; k<=cptcovprod;k++)
cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
ij=1;/* To be checked else nbcode[0][0] wrong */
for(j=3; j <=ncovmodel; j++) {
- /* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /\* Bug valgrind *\/ */
- /* /\*fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);*\/ */
- /* ij++; */
- /* } */
- /* else */
+ if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { /* Bug valgrind */
+ fprintf(ficgp,"+p%d*%d*x",i+j-1,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+ ij++;
+ }
+ else
fprintf(ficgp,"+p%d*%d",i+j-1,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
}
fprintf(ficgp,")/(1");
fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
ij=1;
for(j=3; j <=ncovmodel; j++){
- /* if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) { */
- /* fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]); */
- /* ij++; */
- /* } */
- /* else */
+ if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
+ fprintf(ficgp,"+p%d*%d*x",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][Tvar[j-2]]]);
+ ij++;
+ }
+ else
fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
}
fprintf(ficgp,")");
/* 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<=cptcovage;k++) cov[2+Tage[k]]=nbcode[Tvar[Tage[k]]][codtab[ij][Tvar[Tage[k]]]]*cov[2]; */
/*
* Treating invertedly V2+V1+V3*age+V2*V4 is as if written V2*V4 +V3*age + V1 + V2 */
for(k=cptcovt; k>=1;k--) /**< Number of covariates */
cutl(stre,strb,strd,'V'); /* strd=V3(input): stre="3" */
Tvar[k]=atoi(stre); /* V2+V1+V4+V3*age Tvar[4]=3 ; V1+V2*age Tvar[2]=2 */
cptcovage++; /* Sums the number of covariates which include age as a product */
- Tage[cptcovage]=k; /* Tage[1] = 4 */
+ Tage[cptcovage]=k; /* Tvar[4]=3, Tage[1] = 4 or V1+V1*age Tvar[2]=1, Tage[1]=2 */
/*printf("stre=%s ", stre);*/
} else if (strcmp(strd,"age")==0) { /* or age*Vn */
cptcovprod--;
/Zc:wchar_t /Zi /O2 /Fd"Release\vc120.pdb" /D "WIN32" /D "NDEBUG" /D
"_CONSOLE" /D "_LIB" /D "_USING_V110_SDK71_" /D "_UNICODE" /D
"UNICODE" /Qipo /Zc:forScope /Gd /Oi /MT /Fa"Release\" /EHsc /nologo
- /Fo"Release\" /Qprof-dir "Release\" /Fp"Release\IMaCh.pch"
- /* 64 bits */
+ /Fo"Release\" /Qprof-dir "Release\" /Fp"Release\IMaCh.pch"
+ */
+ /* 64 bits */
/*
/GS /W3 /Gy
/Zc:wchar_t /Zi /O2 /Fd"x64\Release\vc120.pdb" /D "WIN32" /D "NDEBUG"
/Zc:forScope /Oi /MD /Fa"x64\Release\" /EHsc /nologo /Qparallel
/Fo"x64\Release\" /Qprof-dir "x64\Release\" /Fp"x64\Release\IMaCh.pch"
*/
- /* Link is $/ /* /OUT:"visual studio
+ /* Link is */ /* /OUT:"visual studio
2013\Projects\IMaCh\Release\IMaCh.exe" /MANIFEST /NXCOMPAT
/PDB:"visual studio
2013\Projects\IMaCh\Release\IMaCh.pdb" /DYNAMICBASE
/*-------- arguments in the command line --------*/
- /* Log file */
+ /* Main Log file */
strcat(filelog, optionfilefiname);
strcat(filelog,".log"); /* */
if((ficlog=fopen(filelog,"w"))==NULL) {
strcat(fileres, optionfilefiname);
strcat(fileres,".txt"); /* Other files have txt extension */
- /*---------arguments file --------*/
+ /* Main ---------arguments file --------*/
if((ficpar=fopen(optionfile,"r"))==NULL) {
printf("Problem with optionfile '%s' with errno='%s'\n",optionfile,strerror(errno));
fscanf(ficpar,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%lf stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d model=%s\n",title, datafile, &lastobs, &firstpass,&lastpass,&ftol, &stepm, &ncovcol, &nlstate,&ndeath, &maxwav, &mle, &weightopt,model);
numlinepar++;
- printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model);
+ /* printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncovcol, nlstate,ndeath, maxwav, mle, weightopt,model); */
+ printf("title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\n",title, datafile, lastobs, firstpass,lastpass);
+ /*
+
+
+
+ */
+ printf("\nftol=%e \n", ftol);
+ printf("stepm=%d \n", stepm);
+ printf("ncovcol=%d nlstate=%d \n", ncovcol, nlstate);
+ printf("ndeath=%d maxwav=%d mle=%d weight=%d\n", ndeath, maxwav, mle, weightopt);
+ printf("model=%s\n",model);
fprintf(ficparo,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
fprintf(ficlog,"title=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncovcol=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol,stepm,ncovcol,nlstate,ndeath,maxwav, mle, weightopt,model);
fflush(ficlog);
goto end;
exit(0);
}
- else if(mle==-3) {
+ else if(mle==-3) { /* Main Wizard */
prwizard(ncovmodel, nlstate, ndeath, model, ficparo);
printf(" You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
fprintf(ficlog," You choose mle=-3, look at file %s for a template of covariance matrix \n",filereso);
fprintf(ficres,"#%s\n",version);
} /* End of mle != -3 */
-
+ /* Main data
+ */
n= lastobs;
num=lvector(1,n);
moisnais=vector(1,n);
Tage[1=V3*age]= 4; Tage[2=age*V4] = 3
*/
+/* Main decodemodel */
+
if(decodemodel(model, lastobs) == 1)
goto end;
Ndum =ivector(-1,NCOVMAX);
if (ncovmodel > 2)
tricode(Tvar,nbcode,imx, Ndum); /**< Fills nbcode[Tvar[j]][l]; */
+ /* Nbcode gives the value of the lth modality of jth covariate, in
+ V2+V1*age, there are 3 covariates Tvar[2]=1 (V1).*/
+ /* 1 to ncodemax[j] is the maximum value of this jth covariate */
codtab=imatrix(1,100,1,10); /* codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) */
/*printf(" codtab[1,1],codtab[100,10]=%d,%d\n", codtab[1][1],codtab[100][10]);*/
+ /* codtab gives the value 1 or 2 of the hth combination of k covariates (1 or 2).*/
h=0;
if (h>m)
h=1;
/**< codtab(h,k) k = codtab[h,k]=( (h-1) - mod(k-1,2**(k-1) )/2**(k-1) + 1
- * h 1 2 3 4
+ * For k=4 covariates, h goes from 1 to 2**k
+ * codtabm(h,k)= 1 & (h-1) >> (k-1) ;
+ * h\k 1 2 3 4
*______________________________
* 1 i=1 1 i=1 1 i=1 1 i=1 1
* 2 2 1 1 1
* 16 2 2 2 1
*/
codtab[h][k]=j;
+ /* codtab[12][3]=1; */
/*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]]);
}
- /*------------ gnuplot -------------*/
+ /* Initialisation of ----------- gnuplot -------------*/
strcpy(optionfilegnuplot,optionfilefiname);
if(mle==-3)
strcat(optionfilegnuplot,"-mort");
fprintf(ficgp,"set datafile missing 'NaNq'\n");
}
/* fclose(ficgp);*/
- /*--------- index.htm --------*/
+
+
+ /* Initialisation of --------- index.htm --------*/
strcpy(optionfilehtm,optionfilefiname); /* Main html file */
if(mle==-3)
p=param[1][1]; /* *(*(*(param +1)+1)+0) */
globpr=0; /* To get the number ipmx of contributions and the sum of weights*/
-
+ /* For mortality only */
if (mle==-3){
ximort=matrix(1,NDIM,1,NDIM);
-/* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */
+ /* ximort=gsl_matrix_alloc(1,NDIM,1,NDIM); */
cens=ivector(1,n);
ageexmed=vector(1,n);
agecens=vector(1,n);
/* Initialize method and iterate */
/* p[1]=0.0268; p[NDIM]=0.083; */
-/* gsl_vector_set(x, 0, 0.0268); */
-/* gsl_vector_set(x, 1, 0.083); */
+ /* gsl_vector_set(x, 0, 0.0268); */
+ /* gsl_vector_set(x, 1, 0.083); */
gsl_vector_set(x, 0, p[1]);
gsl_vector_set(x, 1, p[2]);
free_ivector(dcwave,1,n);
free_matrix(ximort,1,NDIM,1,NDIM);
#endif
- } /* Endof if mle==-3 */
-
+ } /* Endof if mle==-3 mortality only */
+ /* Standard maximisation */
else{ /* For mle >=1 */
globpr=0;/* debug */
+ /* Computes likelihood for initial parameters */
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]);
printf("\n");
- globpr=1; /* to print the contributions */
+ globpr=1; /* again, to print the contributions */
likelione(ficres, p, npar, nlstate, &globpr, &ipmx, &sw, &fretone, funcone); /* Prints the contributions to the likelihood */
printf("Second Likeli=%12.6f ipmx=%ld sw=%12.6f",fretone,ipmx,sw);
for (k=1; k<=npar;k++)
printf(" %d %8.5f",k,p[k]);
printf("\n");
- if(mle>=1){ /* Could be 1 or 2 */
+ if(mle>=1){ /* Could be 1 or 2, Real Maximisation */
mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
}
fprintf(ficres,"# agemin agemax for life expectancy, bage fage (if mle==0 ie no data nor Max likelihood).\n");
fprintf(ficres,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
fprintf(ficparo,"agemin=%.0f agemax=%.0f bage=%.0f fage=%.0f estepm=%d\n",ageminpar,agemaxpar,bage,fage, estepm);
-
+
+ /* Other stuffs, more or less useful */
while((c=getc(ficpar))=='#' && c!= EOF){
ungetc(c,ficpar);
fgets(line, MAXLINE, ficpar);
fclose(ficres);
+ /* Other results (useful)*/
+
+
/*--------------- Prevalence limit (period or stable prevalence) --------------*/
/*#include "prevlim.h"*/ /* Use ficrespl, ficlog */
prlim=matrix(1,nlstate,1,nlstate);
/* fprintf(ficlog,"Warning %d!! You can only forecast the prevalences if the optimization\n has been performed with stepm = 1 (month) instead of %d or model=. instead of '%s'\n", erreur, stepm, model); */
/* } */
}
-
+
+ /* ------ Other prevalence ratios------------ */
/* Computes prevalence between agemin (i.e minimal age computed) and no more ageminpar */
git://henry.ined.fr / .git / commitdiff