--- imach096d/src/imach.c 2001/05/02 17:47:10 1.6
+++ imach096d/src/imach.c 2001/05/02 17:50:24 1.7
@@ -69,7 +69,7 @@
int nvar;
static int cptcov;
-int cptcovn, cptcovage=0;
+int cptcovn, cptcovage=0, cptcoveff=0;
int npar=NPARMAX;
int nlstate=2; /* Number of live states */
int ndeath=1; /* Number of dead states */
@@ -131,12 +131,12 @@ double **pmmij;
double *weight;
int **s; /* Status */
double *agedc, **covar, idx;
-int **nbcode, *Tcode, *Tvar, **codtab;
+int **nbcode, *Tcode, *Tvar, **codtab, **Tvard, *Tprod, cptcovprod, *Tvaraff;
double ftol=FTOL; /* Tolerance for computing Max Likelihood */
double ftolhess; /* Tolerance for computing hessian */
-
+/**************** split *************************/
static int split( char *path, char *dirc, char *name )
{
char *s; /* pointer */
@@ -657,13 +657,18 @@ double **prevalim(double **prlim, int nl
cov[2]=agefin;
for (k=1; k<=cptcovn;k++) {
- cov[2+k]=nbcode[Tvar[k]][codtab[ij][k]];
-
-/*printf("Tcode[ij]=%d nbcode=%d\n",Tcode[ij],nbcode[k][Tcode[ij]]);*/
+ cov[2+k]=nbcode[Tvar[k]][codtab[ij][Tvar[k]]];
+ /*printf("ij=%d Tvar[k]=%d nbcode=%d cov=%lf\n",ij, Tvar[k],nbcode[Tvar[k]][codtab[ij][Tvar[k]]],cov[2+k]);*/
}
for (k=1; k<=cptcovage;k++)
cov[2+Tage[k]]=cov[2+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]]];
+
+ /*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]);*/
+
out=matprod2(newm, pmij(pmmij,cov,ncovmodel,x,nlstate),1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath, oldm);
savm=oldm;
@@ -798,9 +803,13 @@ double ***hpxij(double ***po, int nhstep
/* Covariates have to be included here again */
cov[1]=1.;
cov[2]=age+((h-1)*hstepm + (d-1))*stepm/YEARM;
- if (cptcovn>0){
- for (k=1; k<=cptcovn;k++) cov[2+k]=nbcode[Tvar[k]][codtab[ij][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++)
+ cov[2+Tprod[k]]=nbcode[Tvard[k][1]][codtab[ij][Tvard[k][1]]]*nbcode[Tvard[k][2]][codtab[ij][Tvard[k][2]]];
+
+
/*printf("hxi cptcov=%d cptcode=%d\n",cptcov,cptcode);*/
/*printf("h=%d d=%d age=%f cov=%f\n",h,d,age,cov[2]);*/
out=matprod2(newm,oldm,1,nlstate+ndeath,1,nlstate+ndeath,1,nlstate+ndeath,
@@ -1160,7 +1169,7 @@ void freqsummary(char fileres[], int ag
freq= ma3x(-1,nlstate+ndeath,-1,nlstate+ndeath,agemin,agemax+3);
j1=0;
- j=cptcovn;
+ j=cptcoveff;
if (cptcovn<1) {j=1;ncodemax[1]=1;}
for(k1=1; k1<=j;k1++){
@@ -1175,8 +1184,8 @@ void freqsummary(char fileres[], int ag
for (i=1; i<=imx; i++) {
bool=1;
if (cptcovn>0) {
- for (z1=1; z1<=cptcovn; z1++)
- if (covar[Tvar[z1]][i]!= nbcode[Tvar[z1]][codtab[j1][z1]]) bool=0;
+ for (z1=1; z1<=cptcoveff; z1++)
+ if (covar[Tvaraff[z1]][i]!= nbcode[Tvaraff[z1]][codtab[j1][z1]]) bool=0;
}
if (bool==1) {
for(m=firstpass; m<=lastpass-1; m++){
@@ -1188,10 +1197,10 @@ void freqsummary(char fileres[], int ag
}
}
if (cptcovn>0) {
- fprintf(ficresp, "\n#Variable");
- for (z1=1; z1<=cptcovn; z1++) fprintf(ficresp, " V%d=%d",Tvar[z1],nbcode[Tvar[z1]][codtab[j1][z1]]);
+ fprintf(ficresp, "\n#********** Variable ");
+ for (z1=1; z1<=cptcoveff; z1++) fprintf(ficresp, "V%d=%d ",Tvaraff[z1],nbcode[Tvaraff[z1]][codtab[j1][z1]]);
}
- fprintf(ficresp, "\n#");
+ fprintf(ficresp, "**********\n#");
for(i=1; i<=nlstate;i++)
fprintf(ficresp, " Age Prev(%d) N(%d) N",i,i);
fprintf(ficresp, "\n");
@@ -1318,34 +1327,56 @@ float sum=0.;
/*********** Tricode ****************************/
void tricode(int *Tvar, int **nbcode, int imx)
{
- int Ndum[80],ij=1, k, j, i, Ntvar[20];
+ int Ndum[20],ij=1, k, j, i;
int cptcode=0;
- for (k=0; k<79; k++) Ndum[k]=0;
+ cptcoveff=0;
+
+ for (k=0; k<19; k++) Ndum[k]=0;
for (k=1; k<=7; k++) ncodemax[k]=0;
- for (j=1; j<=cptcovn; j++) {
+ for (j=1; j<=(cptcovn+2*cptcovprod); j++) {
for (i=1; i<=imx; i++) {
ij=(int)(covar[Tvar[j]][i]);
Ndum[ij]++;
if (ij > cptcode) cptcode=ij;
}
+
/*printf("cptcode=%d cptcovn=%d ",cptcode,cptcovn);*/
for (i=0; i<=cptcode; i++) {
if(Ndum[i]!=0) ncodemax[j]++;
}
-
ij=1;
+
for (i=1; i<=ncodemax[j]; i++) {
- for (k=0; k<=79; k++) {
+ for (k=0; k<=19; k++) {
if (Ndum[k] != 0) {
nbcode[Tvar[j]][ij]=k;
+ /* printf("ij=%d ",nbcode[Tvar[2]][1]);*/
ij++;
}
if (ij > ncodemax[j]) break;
}
}
- }
-
+ }
+ for (i=1; i<=10; i++) {
+ ij=Tvar[i];
+ Ndum[ij]++;
+ }
+ ij=1;
+ for (i=1; i<=cptcovn; i++) {
+ if((Ndum[i]!=0) && (i<=ncov)){
+ Tvaraff[i]=ij;
+ ij++;
+ }
+ }
+
+ for (j=1; j<=(cptcovn+2*cptcovprod); j++) {
+ if ((Tvar[j]>= cptcoveff) && (Tvar[j] <=ncov)) cptcoveff=Tvar[j];
+ /*printf("j=%d %d\n",j,Tvar[j]);*/
+ }
+
+ /* printf("cptcoveff=%d Tvaraff=%d %d\n",cptcoveff, Tvaraff[1],Tvaraff[2]);
+ scanf("%d",i);*/
}
/*********** Health Expectancies ****************/
@@ -1620,7 +1651,7 @@ int main()
int sdeb, sfin; /* Status at beginning and end */
int c, h , cpt,l;
int ju,jl, mi;
- int i1,j1, k1,k2,k3,jk,aa,bb, stepsize;
+ int i1,j1, k1,k2,k3,jk,aa,bb, stepsize, ij;
int jnais,jdc,jint4,jint1,jint2,jint3,**outcome,**adl,*tab;
int hstepm, nhstepm;
@@ -1858,43 +1889,59 @@ split(pathtot, path,optionfile);
imx=i-1; /* Number of individuals */
/* Calculation of the number of parameter from char model*/
- Tvar=ivector(1,15);
+ Tvar=ivector(1,15);
+ Tprod=ivector(1,15);
+ Tvaraff=ivector(1,15);
+ Tvard=imatrix(1,15,1,2);
Tage=ivector(1,15);
if (strlen(model) >1){
- j=0, j1=0;
+ j=0, j1=0, k1=1, k2=1;
j=nbocc(model,'+');
j1=nbocc(model,'*');
cptcovn=j+1;
+ cptcovprod=j1;
strcpy(modelsav,model);
- if (j==0) {
+ if (j==0) {
if (j1==0){
- cutv(stra,strb,modelsav,'V');
- Tvar[1]=atoi(strb);
+ cutv(stra,strb,modelsav,'V');
+ Tvar[1]=atoi(strb);
}
else if (j1==1) {
- cutv(stra,strb,modelsav,'*');
- /* printf("stra=%s strb=%s modelsav=%s ",stra,strb,modelsav);*/
- Tage[1]=1; cptcovage++;
- if (strcmp(stra,"age")==0) {
- cutv(strd,strc,strb,'V');
- Tvar[1]=atoi(strc);
- }
- else if (strcmp(strb,"age")==0) {
- cutv(strd,strc,stra,'V');
- Tvar[1]=atoi(strc);
- }
- else {printf("Error"); exit(0);}
+ cutv(stra,strb,modelsav,'*');
+ Tage[1]=1; cptcovage++;
+ if (strcmp(stra,"age")==0) {
+ cptcovprod--;
+ cutv(strd,strc,strb,'V');
+ Tvar[1]=atoi(strc);
+ }
+ else if (strcmp(strb,"age")==0) {
+ cptcovprod--;
+ cutv(strd,strc,stra,'V');
+ Tvar[1]=atoi(strc);
+ }
+ else {
+ cutv(strd,strc,strb,'V');
+ cutv(stre,strd,stra,'V');
+ Tvar[1]=ncov+1;
+ for (k=1; k<=lastobs;k++)
+ covar[ncov+1][k]=covar[atoi(strc)][k]*covar[atoi(strd)][k];
+ }
+ /*printf("%s %s %s\n", stra,strb,modelsav);
+printf("%d ",Tvar[1]);
+scanf("%d",i);*/
}
}
- else {
+ else {
for(i=j; i>=1;i--){
cutv(stra,strb,modelsav,'+');
- /*printf("%s %s %s\n", stra,strb,modelsav);*/
+ /*printf("%s %s %s\n", stra,strb,modelsav);
+ scanf("%d",i);*/
if (strchr(strb,'*')) {
cutv(strd,strc,strb,'*');
if (strcmp(strc,"age")==0) {
+ cptcovprod--;
cutv(strb,stre,strd,'V');
Tvar[i+1]=atoi(stre);
cptcovage++;
@@ -1902,6 +1949,7 @@ split(pathtot, path,optionfile);
printf("stre=%s ", stre);
}
else if (strcmp(strd,"age")==0) {
+ cptcovprod--;
cutv(strb,stre,strc,'V');
Tvar[i+1]=atoi(stre);
cptcovage++;
@@ -1909,17 +1957,23 @@ split(pathtot, path,optionfile);
}
else {
cutv(strb,stre,strc,'V');
- Tvar[i+1]=ncov+1;
+ Tvar[i+1]=ncov+k1;
cutv(strb,strc,strd,'V');
+ Tprod[k1]=i+1;
+ Tvard[k1][1]=atoi(strc);
+ Tvard[k1][2]=atoi(stre);
+ Tvar[cptcovn+k2]=Tvard[k1][1];
+ Tvar[cptcovn+k2+1]=Tvard[k1][2];
for (k=1; k<=lastobs;k++)
- covar[ncov+1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
+ covar[ncov+k1][k]=covar[atoi(stre)][k]*covar[atoi(strc)][k];
+ k1++;
+ k2=k2+2;
}
}
else {
cutv(strd,strc,strb,'V');
/* printf("%s %s %s", strd,strc,strb);*/
-
- Tvar[i+1]=atoi(strc);
+ Tvar[i+1]=atoi(strc);
}
strcpy(modelsav,stra);
}
@@ -1927,12 +1981,14 @@ split(pathtot, path,optionfile);
Tvar[1]=atoi(strc);
}
}
-
- /* printf("tvar=%d %d cptcovage=%d %d",Tvar[1],Tvar[2],cptcovage,Tage[1]);
- scanf("%d ",i);*/
+ /* for (i=1; i<=5; i++)
+ printf("i=%d %d ",i,Tvar[i]);*/
+ /* printf("tvar=%d %d cptcovage=%d %d",Tvar[1],Tvar[2],cptcovage,Tage[1]);*/
+ /*printf("cptcovprod=%d ", cptcovprod);*/
+ /* scanf("%d ",i);*/
fclose(fic);
- if(mle==1){
+ /* if(mle==1){*/
if (weightopt != 1) { /* Maximisation without weights*/
for(i=1;i<=n;i++) weight[i]=1.0;
}
@@ -2007,18 +2063,18 @@ printf("Total number of individuals= %d,
Tcode=ivector(1,100);
- nbcode=imatrix(1,nvar,1,8);
- ncodemax[1]=1;
- if (cptcovn > 0) tricode(Tvar,nbcode,imx);
-
+ nbcode=imatrix(1,nvar,1,8);
+ ncodemax[1]=1;
+ if (cptcovn > 0) tricode(Tvar,nbcode,imx);
+
codtab=imatrix(1,100,1,10);
h=0;
- m=pow(2,cptcovn);
+ m=pow(2,cptcoveff);
- for(k=1;k<=cptcovn; k++){
+ 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,cptcovn+1-k)); cpt++){
+ for(cpt=1; cpt <=(m/pow(2,cptcoveff+1-k)); cpt++){
h++;
if (h>m) h=1;codtab[h][k]=j;
}
@@ -2026,9 +2082,10 @@ printf("Total number of individuals= %d,
}
}
- /* for(i=1; i <=m ;i++){
+
+ /*for(i=1; i <=m ;i++){
for(k=1; k <=cptcovn; k++){
- printf("i=%d k=%d %d ",i,k,codtab[i][k]);
+ printf("i=%d k=%d %d %d",i,k,codtab[i][k], cptcoveff);
}
printf("\n");
}
@@ -2036,7 +2093,7 @@ printf("Total number of individuals= %d,
/* Calculates basic frequencies. Computes observed prevalence at single age
and prints on file fileres'p'. */
- freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax);
+ freqsummary(fileres, agemin, agemax, s, agev, nlstate, imx,Tvar,nbcode, ncodemax);
pmmij= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
oldms= matrix(1,nlstate+ndeath,1,nlstate+ndeath); /* creation */
@@ -2047,9 +2104,10 @@ printf("Total number of individuals= %d,
/* For Powell, parameters are in a vector p[] starting at p[1]
so we point p on param[1][1] so that p[1] maps on param[1][1][1] */
p=param[1][1]; /* *(*(*(param +1)+1)+0) */
-
- mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
+ if(mle==1){
+ mlikeli(ficres,p, npar, ncovmodel, nlstate, ftol, func);
+ }
/*--------- results files --------------*/
fprintf(ficres,"\ntitle=%s datafile=%s lastobs=%d firstpass=%d lastpass=%d\nftol=%e stepm=%d ncov=%d nlstate=%d ndeath=%d maxwav=%d mle=%d weight=%d\nmodel=%s\n", title, datafile, lastobs, firstpass,lastpass,ftol, stepm, ncov, nlstate, ndeath, maxwav, mle,weightopt,model);
@@ -2073,10 +2131,11 @@ printf("Total number of individuals= %d,
}
}
}
-
+ if(mle==1){
/* Computing hessian and covariance matrix */
ftolhess=ftol; /* Usually correct */
hesscov(matcov, p, npar, delti, ftolhess, func);
+ }
fprintf(ficres,"# Scales\n");
printf("# Scales\n");
for(i=1,jk=1; i <=nlstate; i++){
@@ -2131,6 +2190,8 @@ printf("Total number of individuals= %d,
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\n",agemin,agemax,bage,fage);
+
+
/*------------ gnuplot -------------*/
chdir(pathcd);
if((ficgp=fopen("graph.plt","w"))==NULL) {
@@ -2139,7 +2200,7 @@ chdir(pathcd);
#ifdef windows
fprintf(ficgp,"cd \"%s\" \n",pathc);
#endif
-m=pow(2,cptcovn);
+m=pow(2,cptcoveff);
/* 1eme*/
for (cpt=1; cpt<= nlstate ; cpt ++) {
@@ -2260,15 +2321,28 @@ fprintf(ficgp,"\nset out \"v%s%d%d.gif\"
for(k=1; k<=(nlstate+ndeath); k++) {
if (k != k2){
fprintf(ficgp," exp(p%d+p%d*x",i,i+1);
-
- for(j=3; j <=ncovmodel; j++)
+ij=1;
+ for(j=3; j <=ncovmodel; j++) {
+ if(((j-2)==Tage[ij]) &&(ij <=cptcovage)) {
+ 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,")/(1");
for(k1=1; k1 <=nlstate; k1++){
fprintf(ficgp,"+exp(p%d+p%d*x",k3+(k1-1)*ncovmodel,k3+(k1-1)*ncovmodel+1);
- for(j=3; j <=ncovmodel; j++)
+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
fprintf(ficgp,"+p%d*%d",k3+(k1-1)*ncovmodel+1+j-2,nbcode[Tvar[j-2]][codtab[jk][j-2]]);
+ }
fprintf(ficgp,")");
}
fprintf(ficgp,") t \"p%d%d\" ", k2,k);
@@ -2297,7 +2371,7 @@ chdir(path);
/*free_matrix(covar,1,NCOVMAX,1,n);*/
fclose(ficparo);
fclose(ficres);
- }
+ /* }*/
/*________fin mle=1_________*/
@@ -2335,7 +2409,7 @@ chdir(path);
fprintf(fichtm," Graphs\n");
- m=cptcovn;
+ m=cptcoveff;
if (cptcovn < 1) {m=1;ncodemax[1]=1;}
j1=0;
@@ -2344,8 +2418,8 @@ chdir(path);
j1++;
if (cptcovn > 0) {
fprintf(fichtm,"
************ Results for covariates");
- for (cpt=1; cpt<=cptcovn;cpt++)
- fprintf(fichtm," V%d=%d ",Tvar[cpt],nbcode[Tvar[cpt]][codtab[j1][cpt]]);
+ for (cpt=1; cpt<=cptcoveff;cpt++)
+ fprintf(fichtm," V%d=%d ",Tvaraff[cpt],nbcode[Tvaraff[cpt]][codtab[j1][cpt]]);
fprintf(fichtm," ************\n
");
}
fprintf(fichtm,"
- Probabilities: pe%s%d.gif
@@ -2394,7 +2468,7 @@ fclose(fichtm);
agebase=agemin;
agelim=agemax;
ftolpl=1.e-10;
- i1=cptcovn;
+ i1=cptcoveff;
if (cptcovn < 1){i1=1;}
for(cptcov=1;cptcov<=i1;cptcov++){
@@ -2402,8 +2476,8 @@ fclose(fichtm);
k=k+1;
/*printf("cptcov=%d cptcod=%d codtab=%d nbcode=%d\n",cptcov, cptcod,Tcode[cptcode],codtab[cptcod][cptcov]);*/
fprintf(ficrespl,"\n#******");
- for(j=1;j<=cptcovn;j++)
- fprintf(ficrespl," V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficrespl," V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficrespl,"******\n");
for (age=agebase; age<=agelim; age++){
@@ -2436,8 +2510,8 @@ fclose(fichtm);
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;
fprintf(ficrespij,"\n#****** ");
- for(j=1;j<=cptcovn;j++)
- fprintf(ficrespij,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficrespij,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficrespij,"******\n");
for (agedeb=fage; agedeb>=bage; agedeb--){ /* If stepm=6 months */
@@ -2495,17 +2569,17 @@ fclose(fichtm);
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;
fprintf(ficrest,"\n#****** ");
- for(j=1;j<=cptcovn;j++)
- fprintf(ficrest,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficrest,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficrest,"******\n");
fprintf(ficreseij,"\n#****** ");
- for(j=1;j<=cptcovn;j++)
+ for(j=1;j<=cptcoveff;j++)
fprintf(ficreseij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);
fprintf(ficreseij,"******\n");
fprintf(ficresvij,"\n#****** ");
- for(j=1;j<=cptcovn;j++)
+ for(j=1;j<=cptcoveff;j++)
fprintf(ficresvij,"V%d=%d ",j,nbcode[j][codtab[k][j]]);
fprintf(ficresvij,"******\n");
@@ -2566,8 +2640,8 @@ strcpy(fileresvpl,"vpl");
for(cptcod=1;cptcod<=ncodemax[cptcov];cptcod++){
k=k+1;
fprintf(ficresvpl,"\n#****** ");
- for(j=1;j<=cptcovn;j++)
- fprintf(ficresvpl,"V%d=%d ",Tvar[j],nbcode[Tvar[j]][codtab[k][j]]);
+ for(j=1;j<=cptcoveff;j++)
+ fprintf(ficresvpl,"V%d=%d ",Tvaraff[j],nbcode[Tvaraff[j]][codtab[k][j]]);
fprintf(ficresvpl,"******\n");
varpl=matrix(1,nlstate,(int) bage, (int) fage);