From: N. Brouard Date: Sun, 21 May 2023 05:05:17 +0000 (+0000) Subject: Summary: Temporary change for imachprax X-Git-Tag: imach-099s7~47 X-Git-Url: https://henry.ined.fr/git/?a=commitdiff_plain;h=ecc0c07372b76ac46908c2dc0be1a85ccfcad5d4;p=.git Summary: Temporary change for imachprax --- diff --git a/src/imach.c b/src/imach.c index bf41793..fbbde0c 100644 --- a/src/imach.c +++ b/src/imach.c @@ -1,6 +1,9 @@ /* $Id$ $State$ $Log$ + Revision 1.353 2023/05/08 18:48:22 brouard + *** empty log message *** + Revision 1.352 2023/04/29 10:46:21 brouard *** empty log message *** @@ -1609,7 +1612,7 @@ int **nbcode, *Tvar; /**< model=V2 => Tvar[1]= 2 */ /* Tage[cptcovage]=k 5 8 10 */ /* Position in the model of ith cov*age */ /* model="V2+V3+V4+V6+V7+V6*V2+V7*V2+V6*V3+V7*V3+V6*V4+V7*V4+age*V2+age*V3+age*V4+age*V6+age*V7+age*V6*V2+age*V6*V3+age*V7*V3+age*V6*V4+age*V7*V4\r"*/ /* p Tvard[1][1]@21 = {6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0}*/ -/* p Tvard[2][1]@21 = {7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0 } +/* p Tvard[2][1]@21 = {7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0 } */ /* p Tvardk[1][1]@24 = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0}*/ /* p Tvardk[1][1]@22 = {0, 0, 0, 0, 0, 0, 0, 0, 6, 2, 7, 2, 6, 3, 7, 3, 6, 4, 7, 4, 0, 0} */ /* Tvard[1][1]@4={4,3,1,2} V4*V3 V1*V2 */ /* Position in model of the ith prod without age */ @@ -4550,7 +4553,7 @@ double funcone( double *x) * 3 ncovta=15 +age*V3*V2+age*V2+agev3+ageV4 +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 * 3 TvarAVVA[1]@15= itva 3 2 2 3 4 6 7 6 3 7 3 6 4 7 4 * 3 ncovta 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 - * TvarAVVAind[1]@15= V3 is in k=2 1 1 2 3 4 5 4,2 5,2, 4,3 5 3}TvarVVAind[] + *?TvarAVVAind[1]@15= V3 is in k=2 1 1 2 3 4 5 4,2 5,2, 4,3 5 3}TvarVVAind[] * TvarAVVAind[1]@15= V3 is in k=6 6 12 13 14 15 16 18 18 19,19, 20,20 21,21}TvarVVAind[] * 3 ncovvta=10 +age*V6 + age*V7 + age*V6*V3 +age*V7*V3 + age*V6*V4 +age*V7*V4 * 3 we want to compute =cotvar[mw[mi][i]][TvarVVA[ncovva]][i] at position TvarVVAind[ncovva] @@ -4565,6 +4568,7 @@ double funcone( double *x) * 6, 8, 9, 10, 11} * TvarFind[itv] 0 0 0 * FixedV[itv] 1 1 1 0 1 0 1 0 1 0 0 + *? FixedV[itv] 1 1 1 0 1 0 1 0 1 0 1 0 1 0 * Tvar[TvarFind[ncovf]]=[1]=2 [2]=3 [4]=4 * Tvar[TvarFind[itv]] [0]=? ?ncovv 1 à ncovvt] * Not a fixed cotvar[mw][itv][i] 6 7 6 2 7, 2, 6, 3, 7, 3, 6, 4, 7, 4} @@ -4576,10 +4580,14 @@ double funcone( double *x) ipos=TvarVVind[ncovv]; /* TvarVVind={2, 5, 5] gives the position in the model of the ncovv th varying covariate*/ /* if(TvarFind[itv]==0){ /\* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv *\/ */ if(FixedV[itv]!=0){ /* Not a fixed covariate? Could be a fixed covariate of a product with a higher than ncovcol+nqv, itv */ + /* printf("DEBUG ncovv=%d, Varying TvarVV[ncovv]=%d\n",ncovv, TvarVV[ncovv]); */ cotvarv=cotvar[mw[mi][i]][TvarVV[ncovv]][i]; /* because cotvar starts now at first ncovcol+nqv+ntv+nqtv (1 to nqtv) */ + /* printf("DEBUG Varying cov[ioffset+ipos=%d]=%g \n",ioffset+ipos,cotvarv); */ }else{ /* fixed covariate */ /* cotvarv=covar[Tvar[TvarFind[itv]]][i]; /\* Error: TvarFind gives the name, that is the true column of fixed covariates, but Tvar of the model *\/ */ + /* printf("DEBUG ncovv=%d, Fixed TvarVV[ncovv]=%d\n",ncovv, TvarVV[ncovv]); */ cotvarv=covar[itv][i]; /* Good: In V6*V3, 3 is fixed at position of the data */ + /* printf("DEBUG Fixed cov[ioffset+ipos=%d]=%g \n",ioffset+ipos,cotvarv); */ } if(ipos!=iposold){ /* Not a product or first of a product */ cotvarvold=cotvarv; @@ -4588,6 +4596,7 @@ double funcone( double *x) } iposold=ipos; cov[ioffset+ipos]=cotvarv; + /* printf("DEBUG Product cov[ioffset+ipos=%d] \n",ioffset+ipos); */ /* For products */ } /* for(itv=1; itv <= ntveff; itv++){ /\* Varying dummy covariates single *\/ */ @@ -4930,14 +4939,15 @@ void mlikeli(FILE *ficres,double p[], int npar, int ncovmodel, int nlstate, doub double fret; double fretone; /* Only one call to likelihood */ /* char filerespow[FILENAMELENGTH];*/ - + + double * p1; /* Shifted parameters from 0 instead of 1 */ #ifdef NLOPT int creturn; nlopt_opt opt; /* double lb[9] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL, -HUGE_VAL }; /\* lower bounds *\/ */ double *lb; double minf; /* the minimum objective value, upon return */ - double * p1; /* Shifted parameters from 0 instead of 1 */ + myfunc_data dinst, *d = &dinst; #endif @@ -5272,6 +5282,7 @@ double hessij( double x[], double **hess, double delti[], int thetai,int thetaj, kmax=kmax+10; } if(kmax >=10 || firstime ==1){ + /* What are the thetai and thetaj? thetai/ncovmodel thetai=(thetai-thetai%ncovmodel)/ncovmodel +thetai%ncovmodel=(line,pos) */ printf("Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol); fprintf(ficlog,"Warning: directions %d-%d, you are not estimating the Hessian at the exact maximum likelihood; you could increase ftol=%.2e\n",thetai,thetaj, ftol); printf("%d %d k=%d, k1=%.12e k2=%.12e k3=%.12e k4=%.12e delti*k=%.12e deltj*k=%.12e, xi-de*k=%.12e xj-de*k=%.12e res=%.12e k1234=%.12e,k1-2=%.12e,k3-4=%.12e\n",thetai,thetaj,k,k1,k2,k3,k4,delti[thetai]/k,delti[thetaj]/k,x[thetai]-delti[thetai]/k,x[thetaj]-delti[thetaj]/k, res,k1-k2-k3+k4,k1-k2,k3-k4); @@ -8357,8 +8368,8 @@ true period expectancies (those weighted with period prevalences are also\ /******************* Gnuplot file **************/ void printinggnuplot(char fileresu[], char optionfilefiname[], double ageminpar, double agemaxpar, double bage, double fage , int prevfcast, int prevbcast, char pathc[], double p[], int offyear, int offbyear){ - char dirfileres[132],optfileres[132]; - char gplotcondition[132], gplotlabel[132]; + char dirfileres[256],optfileres[256]; + char gplotcondition[256], gplotlabel[256]; int cpt=0,k1=0,i=0,k=0,j=0,jk=0,k2=0,k3=0,k4=0,kf=0,kvar=0,kk=0,ipos=0,iposold=0,ij=0, ijp=0, l=0; int lv=0, vlv=0, kl=0; int ng=0;