From 354969098a1efe8af7f10f14438782b883f8ea44 Mon Sep 17 00:00:00 2001 From: Philippe Verley <philippe.verley@ird.fr> Date: Mon, 1 Aug 2016 09:23:40 +0000 Subject: [PATCH] (Zlibs.c) Updated all the mimetic_* functions to handle properly spatial point patterns with negative coordinates. --- src/Zlibs.c | 594 +++++++++++++++++++++++++++++----------------------- 1 file changed, 333 insertions(+), 261 deletions(-) diff --git a/src/Zlibs.c b/src/Zlibs.c index bbc86e3..208818a 100755 --- a/src/Zlibs.c +++ b/src/Zlibs.c @@ -1168,7 +1168,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval Rprintf("ERREUR Ripley\n"); } else { - /*comptage du nombre de |ļobs|<=|ļsimu| pour test local*/ + /*comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local*/ double gictmp,kictmp,lictmp,nictmp; for(j=0;j<*t2;j++) { gictmp=gic1[j]/(*densite*(Pi()*(j+1)*(j+1)*(*dt)*(*dt)-Pi()*j*j*(*dt)*(*dt))); @@ -1261,7 +1261,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval Rprintf("ERREUR Ripley\n"); } else - { /*comptage du nombre de |ļobs|<=|ļsimu| pour test local*/ + { /*comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local*/ double gictmp,kictmp,lictmp,nictmp; for(j=0;j<*t2;j++) { gictmp=gic1[j]/(*densite*(Pi()*(j+1)*(j+1)*(*dt)*(*dt)-Pi()*j*j*(*dt)*(*dt))); @@ -1317,7 +1317,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval return -1; } - /* definition de i0 : indice ou sera stocké l'estimation des bornes de l'IC*/ + /* definition de i0 : indice ou sera stockīŋŊ l'estimation des bornes de l'IC*/ i0=*lev/2*(*nbSimu+1); if (i0<1) i0=1; @@ -1355,7 +1355,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval Rprintf("ERREUR Ripley\n"); } else - { /*comptage du nombre de |ļobs|<=|ļsimu| pour test local*/ + { /*comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local*/ double gictmp,kictmp,lictmp,nictmp; for(j=0;j<*t2;j++) { gictmp=gic1[j]/(*densite*(Pi()*(j+1)*(j+1)*(*dt)*(*dt)-Pi()*j*j*(*dt)*(*dt))); @@ -1451,7 +1451,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval Rprintf("ERREUR Ripley\n"); } else { - /*comptage du nombre de |ļobs|<=|ļsimu| pour test local*/ + /*comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local*/ double gictmp,kictmp,lictmp,nictmp; for(j=0;j<*t2;j++) { gictmp=gic1[j]/(*densite*(Pi()*(j+1)*(j+1)*(*dt)*(*dt)-Pi()*j*j*(*dt)*(*dt))); @@ -1464,7 +1464,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval if ((float)fabs(ll[j])<=(float)fabs(lictmp)) {lval[j]+=1;} } - /*Traitement des résultats*/ + /*Traitement des rīŋŊsultats*/ ic(i,i0,gic,kic,gic1,kic1,*t2); } R_FlushConsole(); @@ -1474,7 +1474,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval i1=i0+2; i2=i0; - /*Copies des valeurs dans les tableaux résultats*/ + /*Copies des valeurs dans les tableaux rīŋŊsultats*/ for(i=0;i<*t2;i++) { gic1[i]=gic[i+1][i1]; gic2[i]=gic[i+1][i2]; @@ -2215,7 +2215,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval lt[j]=sqrt(kt[j]/Pi())-(j+1)*(*dt); } } - if (*h0==2) { //Option 2 : initialisation coordonnes points dcals + if (*h0==2) { //Option 2 : initialisation coordonnīŋŊes points dīŋŊcalīŋŊs vecalloc(&x,*point_nb1); vecalloc(&y,*point_nb1); } @@ -2262,9 +2262,9 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval if (erreur==0) { if (*h0==1)//etiquetage aleatoire erreur=intertype_rect(point_nb1,x1,y1,point_nb2,x2,y2,xmi,xma,ymi,yma,t2,dt,gic1,kic1); - if (*h0==2) // décallage avec rectangle + if (*h0==2) // dīŋŊcallage avec rectangle erreur=intertype_rect(&point_nb,x,y,point_nb2,x2,y2,xmi,xma,ymi,yma,t2,dt,gic1,kic1); - if (*h0==3) // mimtique + if (*h0==3) // mimīŋŊtique erreur=intertype_rect(point_nb1,x,y,point_nb2,x2,y2,xmi,xma,ymi,yma,t2,dt,gic1,kic1); } // si il y a une erreur on recommence une simulation @@ -2273,7 +2273,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval Rprintf("ERREUR Intertype\n"); } else { - //comptage du nombre de |ļobs|<=|ļsimu| pour test local + //comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local double gictmp,kictmp,lictmp,nictmp; for(j=0;j<*t2;j++) { gictmp=gic1[j]/(densite_2*(Pi()*(j+1)*(j+1)*(*dt)*(*dt)-Pi()*j*j*(*dt)*(*dt))); @@ -2294,7 +2294,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval } } - //Traitement des résultats + //Traitement des rīŋŊsultats ic(i,i0,gic,kic,gic1,kic1,*t2); } R_FlushConsole(); @@ -2304,7 +2304,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval i1=i0+2; i2=i0; - //Copies des valeurs dans les tableaux résultats + //Copies des valeurs dans les tableaux rīŋŊsultats for(i=0;i<*t2;i++) { gic1[i]=gic[i+1][i1]; gic2[i]=gic[i+1][i2]; @@ -2359,7 +2359,7 @@ int *nsimax, int *conv, int *rep, double *lev,double *g,double *k,double *gic1,d erreur=intertype_disq(point_nb1,x1,y1,point_nb2,x2,y2,x0,y0,r0,t2,dt,g,k); if (erreur!=0) return -1; - /*Définition de i0 : indice oų sera stocké l'estimation des bornes de l'IC*/ + /*DīŋŊfinition de i0 : indice oīŋŊ sera stockīŋŊ l'estimation des bornes de l'IC*/ i0=*lev/2*(*nbSimu+1); if (i0<1) i0=1; @@ -2434,7 +2434,7 @@ int *nsimax, int *conv, int *rep, double *lev,double *g,double *k,double *gic1,d Rprintf("Monte Carlo simulation\n"); for(i=1;i<=*nbSimu;i++) { - /*On simule les hypothčses nulles*/ + /*On simule les hypothīŋŊses nulles*/ if(*h0==1) erreur=randlabelling(x,y,*point_nb1,x1,y1,*point_nb2,x2,y2,type); if(*h0==2) @@ -2460,7 +2460,7 @@ int *nsimax, int *conv, int *rep, double *lev,double *g,double *k,double *gic1,d erreur=intertype_disq(&point_nb,x,y,point_nb2,x2,y2,x0,y0,r0,t2,dt,gic1,kic1); } if (*h0==3) { - // mimtique + // mimīŋŊtique erreur=intertype_disq(point_nb1,x,y,point_nb2,x2,y2,x0,y0,r0,t2,dt,gic1,kic1); } } @@ -2470,7 +2470,7 @@ int *nsimax, int *conv, int *rep, double *lev,double *g,double *k,double *gic1,d Rprintf("ERREUR Intertype\n"); } else { - /*comptage du nombre de |ļobs|<=|ļsimu| pour test local*/ + /*comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local*/ double gictmp,kictmp,lictmp,nictmp; for(j=0;j<*t2;j++) { gictmp=gic1[j]/(densite_2*(Pi()*(j+1)*(j+1)*(*dt)*(*dt)-Pi()*j*j*(*dt)*(*dt))); @@ -2490,7 +2490,7 @@ int *nsimax, int *conv, int *rep, double *lev,double *g,double *k,double *gic1,d if ((float)fabs(ll[j])<=(float)fabs(lictmp)) {lval[j]+=1;} } } - /*Traitement des résultats*/ + /*Traitement des rīŋŊsultats*/ ic(i,i0,gic,kic,gic1,kic1,*t2); } R_FlushConsole(); @@ -2498,7 +2498,7 @@ int *nsimax, int *conv, int *rep, double *lev,double *g,double *k,double *gic1,d } i1=i0+2; i2=i0; - /*Copies des valeurs dans les tableaux résultats*/ + /*Copies des valeurs dans les tableaux rīŋŊsultats*/ for(i=0;i<*t2;i++) { gic1[i]=gic[i+1][i1]; gic2[i]=gic[i+1][i2]; @@ -2553,7 +2553,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval erreur=intertype_tr_rect(point_nb1,x1,y1,point_nb2,x2,y2,xmi,xma,ymi,yma,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,g,k); if (erreur!=0) return -1; - /*Définition de i0 : indice oų sera stocké l'estimation des bornes de l'IC*/ + /*DīŋŊfinition de i0 : indice oīŋŊ sera stockīŋŊ l'estimation des bornes de l'IC*/ i0=*lev/2*(*nbSimu+1); if (i0<1) i0=1; @@ -2631,7 +2631,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval Rprintf("Monte Carlo simulation\n"); for(i=1;i<=*nbSimu;i++) { - /*On simule les hypothčses nulles*/ + /*On simule les hypothīŋŊses nulles*/ if(*h0==1) erreur=randlabelling(x,y,*point_nb1,x1,y1,*point_nb2,x2,y2,type); if(*h0==2) @@ -2651,11 +2651,11 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval } if (erreur==0) { - if (*h0==1) //étiquetage aléatoire + if (*h0==1) //īŋŊtiquetage alīŋŊatoire erreur=intertype_tr_rect(point_nb1,x1,y1,point_nb2,x2,y2,xmi,xma,ymi,yma,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,gic1,kic1); - if (*h0==2) //décallage avec rectangle + if (*h0==2) //dīŋŊcallage avec rectangle erreur=intertype_tr_rect(&point_nb,x,y,point_nb2,x2,y2,xmi,xma,ymi,yma,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,gic1,kic1); - if (*h0==3) // mimtique + if (*h0==3) // mimīŋŊtique erreur=intertype_tr_rect(point_nb1,x,y,point_nb2,x2,y2,xmi,xma,ymi,yma,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,gic1,kic1); } /* si il y a une erreur on recommence une simulation*/ @@ -2664,7 +2664,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval Rprintf("ERREUR Intertype\n"); } else { - /*comptage du nombre de |ļobs|<=|ļsimu| pour test local*/ + /*comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local*/ double gictmp,kictmp,lictmp,nictmp; for(j=0;j<*t2;j++) { gictmp=gic1[j]/(densite_2*(Pi()*(j+1)*(j+1)*(*dt)*(*dt)-Pi()*j*j*(*dt)*(*dt))); @@ -2685,7 +2685,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval } } - /*Traitement des résultats*/ + /*Traitement des rīŋŊsultats*/ ic(i,i0,gic,kic,gic1,kic1,*t2); } R_FlushConsole(); @@ -2695,7 +2695,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval i1=i0+2; i2=i0; - /*Copies des valeurs dans les tableaux résultats*/ + /*Copies des valeurs dans les tableaux rīŋŊsultats*/ for(i=0;i<*t2;i++) { gic1[i]=gic[i+1][i1]; gic2[i]=gic[i+1][i2]; @@ -2750,7 +2750,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval erreur=intertype_tr_disq(point_nb1,x1,y1,point_nb2,x2,y2,x0,y0,r0,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,g,k); if (erreur!=0) return -1; - /*Définition de i0 : indice oų sera stocké l'estimation des bornes de l'IC*/ + /*DīŋŊfinition de i0 : indice oīŋŊ sera stockīŋŊ l'estimation des bornes de l'IC*/ i0=*lev/2*(*nbSimu+1); if (i0<1) i0=1; @@ -2825,7 +2825,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval Rprintf("Monte Carlo simulation\n"); for(i=1;i<=*nbSimu;i++) { - /*On simule les hypothčses nulles*/ + /*On simule les hypothīŋŊses nulles*/ if(*h0==1) erreur=randlabelling(x,y,*point_nb1,x1,y1,*point_nb2,x2,y2,type); if(*h0==2) @@ -2845,11 +2845,11 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval } if (erreur==0) { - if (*h0==1) //étiquetage aléatoire + if (*h0==1) //īŋŊtiquetage alīŋŊatoire erreur=intertype_tr_disq(point_nb1,x1,y1,point_nb2,x2,y2,x0,y0,r0,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,gic1,kic1); - if (*h0==2) //décallage avec rectangle + if (*h0==2) //dīŋŊcallage avec rectangle erreur=intertype_tr_disq(&point_nb,x,y,point_nb2,x2,y2,x0,y0,r0,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,gic1,kic1); - if (*h0==3) // mimtique + if (*h0==3) // mimīŋŊtique erreur=intertype_tr_disq(point_nb1,x,y,point_nb2,x2,y2,x0,y0,r0,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,gic1,kic1); } /* si il y a une erreur on recommence une simulation*/ @@ -2858,7 +2858,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval Rprintf("ERREUR Intertype\n"); } else { - /*comptage du nombre de |ļobs|<=|ļsimu| pour test local*/ + /*comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local*/ double gictmp,kictmp,lictmp,nictmp; for(j=0;j<*t2;j++) { gictmp=gic1[j]/(densite_2*(Pi()*(j+1)*(j+1)*(*dt)*(*dt)-Pi()*j*j*(*dt)*(*dt))); @@ -2879,7 +2879,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval } } - /*Traitement des résultats*/ + /*Traitement des rīŋŊsultats*/ ic(i,i0,gic,kic,gic1,kic1,*t2); } R_FlushConsole(); @@ -2889,7 +2889,7 @@ double *gic1,double *gic2, double *kic1,double *kic2, double *gval, double *kval i1=i0+2; i2=i0; - /*Copies des valeurs dans les tableaux résultats*/ + /*Copies des valeurs dans les tableaux rīŋŊsultats*/ for(i=0;i<*t2;i++) { gic1[i]=gic[i+1][i1]; gic2[i]=gic[i+1][i2]; @@ -3406,8 +3406,8 @@ double *cx, double *cy,double prec) { } /******************************************************************************/ -/* Calcule la fonction de corrlation des marques Km(r) pour un semis (x,y) */ -/* affect d'une marque quantitative c(x,y) en parametres */ +/* Calcule la fonction de corrīŋŊlation des marques Km(r) pour un semis (x,y) */ +/* affectīŋŊ d'une marque quantitative c(x,y) en parametres */ /* dans une zone de forme rectangulaire de bornes xmi xma ymi yma */ /* Les corrections des effets de bords sont fait par la methode de Ripley, */ /* i.e. l'inverse de la proportion d'arc de cercle inclu dans la fenetre. */ @@ -3494,7 +3494,7 @@ int *t2,double *dt,double *gm,double *km) return 0; } -/*function de corrlation dans forme circulaire*/ +/*function de corrīŋŊlation dans forme circulaire*/ int corr_disq(int *point_nb,double *x,double *y,double *c, double *x0,double *y0,double *r0, int *t2,double *dt,double *gm,double *km) { int i,j,tt; @@ -3773,7 +3773,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double return -1; } - /*Définition de i0 : indice oų sera stocké l'estimation des bornes de l'IC*/ + /*DīŋŊfinition de i0 : indice oīŋŊ sera stockīŋŊ l'estimation des bornes de l'IC*/ i0=*lev/2*(*nbSimu+1); if (i0<1) i0=1; @@ -3805,7 +3805,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double Rprintf("ERREUR mark correlation\n"); } else { - /*comptage du nombre de |ļobs|<=|ļsimu| pour test local*/ + /*comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local*/ double gmictmp,kmictmp; for(j=0;j<*t2;j++) { gmictmp=gmic1[j]; @@ -3814,7 +3814,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double if ((float)fabs(kkm[j])<=(float)fabs(kmictmp)) {kmval[j]+=1;} } - /*Traitement des résultats*/ + /*Traitement des rīŋŊsultats*/ ic(i,i0,gmic,kmic,gmic1,kmic1,*t2); } R_FlushConsole(); @@ -3824,7 +3824,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double i1=i0+2; i2=i0; - /*Copies des valeurs dans les tableaux résultats*/ + /*Copies des valeurs dans les tableaux rīŋŊsultats*/ for(i=0;i<*t2;i++) { gmic1[i]=gmic[i+1][i1]; gmic2[i]=gmic[i+1][i2]; @@ -3857,7 +3857,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double return -1; } - /*Définition de i0 : indice oų sera stocké l'estimation des bornes de l'IC*/ + /*DīŋŊfinition de i0 : indice oīŋŊ sera stockīŋŊ l'estimation des bornes de l'IC*/ i0=*lev/2*(*nbSimu+1); if (i0<1) i0=1; @@ -3889,7 +3889,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double Rprintf("ERREUR mark correlation\n"); } else { - /*comptage du nombre de |ļobs|<=|ļsimu| pour test local*/ + /*comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local*/ double gmictmp,kmictmp; for(j=0;j<*t2;j++) { gmictmp=gmic1[j]; @@ -3898,7 +3898,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double if ((float)fabs(kkm[j])<=(float)fabs(kmictmp)) {kmval[j]+=1;} } - /*Traitement des résultats*/ + /*Traitement des rīŋŊsultats*/ ic(i,i0,gmic,kmic,gmic1,kmic1,*t2); } R_FlushConsole(); @@ -3908,7 +3908,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double i1=i0+2; i2=i0; - /*Copies des valeurs dans les tableaux résultats*/ + /*Copies des valeurs dans les tableaux rīŋŊsultats*/ for(i=0;i<*t2;i++) { gmic1[i]=gmic[i+1][i1]; gmic2[i]=gmic[i+1][i2]; @@ -3943,7 +3943,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double return -1; } - /*Définition de i0 : indice oų sera stocké l'estimation des bornes de l'IC*/ + /*DīŋŊfinition de i0 : indice oīŋŊ sera stockīŋŊ l'estimation des bornes de l'IC*/ i0=*lev/2*(*nbSimu+1); if (i0<1) i0=1; @@ -3975,7 +3975,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double Rprintf("ERREUR Intertype\n"); } else { - /*comptage du nombre de |ļobs|<=|ļsimu| pour test local*/ + /*comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local*/ double gmictmp,kmictmp; for(j=0;j<*t2;j++) { gmictmp=gmic1[j]; @@ -3984,7 +3984,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double if ((float)fabs(kkm[j])<=(float)fabs(kmictmp)) {kmval[j]+=1;} } - /*Traitement des résultats*/ + /*Traitement des rīŋŊsultats*/ ic(i,i0,gmic,kmic,gmic1,kmic1,*t2); } R_FlushConsole(); @@ -3994,7 +3994,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double i1=i0+2; i2=i0; - /*Copies des valeurs dans les tableaux résultats*/ + /*Copies des valeurs dans les tableaux rīŋŊsultats*/ for(i=0;i<*t2;i++) { gmic1[i]=gmic[i+1][i1]; gmic2[i]=gmic[i+1][i2]; @@ -4028,7 +4028,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double return -1; } - /*Définition de i0 : indice oų sera stocké l'estimation des bornes de l'IC*/ + /*DīŋŊfinition de i0 : indice oīŋŊ sera stockīŋŊ l'estimation des bornes de l'IC*/ i0=*lev/2*(*nbSimu+1); if (i0<1) i0=1; @@ -4060,7 +4060,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double Rprintf("ERREUR Intertype\n"); } else { - /*comptage du nombre de |ļobs|<=|ļsimu| pour test local*/ + /*comptage du nombre de |īŋŊobs|<=|īŋŊsimu| pour test local*/ double gmictmp,kmictmp; for(j=0;j<*t2;j++) { gmictmp=gmic1[j]; @@ -4069,7 +4069,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double if ((float)fabs(kkm[j])<=(float)fabs(kmictmp)) {kmval[j]+=1;} } - /*Traitement des résultats*/ + /*Traitement des rīŋŊsultats*/ ic(i,i0,gmic,kmic,gmic1,kmic1,*t2); } R_FlushConsole(); @@ -4079,7 +4079,7 @@ double *gmic1,double *gmic2, double *kmic1,double *kmic2, double *gmval, double i1=i0+2; i2=i0; - /*Copies des valeurs dans les tableaux résultats*/ + /*Copies des valeurs dans les tableaux rīŋŊsultats*/ for(i=0;i<*t2;i++) { gmic1[i]=gmic[i+1][i1]; gmic2[i]=gmic[i+1][i2]; @@ -4139,14 +4139,14 @@ void randmark(int point_nb,double *c,double *c2) /*********************************************************************************************/ -/* Fonctions de Shimatani pour les semis de points multivaris + fonctions utilitaires */ +/* Fonctions de Shimatani pour les semis de points multivariīŋŊs + fonctions utilitaires */ /********************************************************************************************/ /// fonction de shimatani pour une zone rectangulaire int shimatani_rect(int *point_nb,double *x,double *y, double *xmi,double *xma,double *ymi, double *yma,int *t2,double *dt,int *nbtype,int *type,double *surface,double *gs, double *ks,int *error) { int i,j,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds; @@ -4254,7 +4254,7 @@ int shimatani_disq(int *point_nb,double *x,double *y, double *x0,double *y0,doub int *t2,double *dt,int *nbtype,int *type,double *surface,double *gs, double *ks,int *error) { int i,j,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds; @@ -4351,7 +4351,7 @@ double *yma,int *triangle_nb, double *ax, double *ay, double *bx, double *by, do int *t2,double *dt,int *nbtype,int *type,double *surface,double *gs, double *ks,int *error) { int i,j,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds; @@ -4455,7 +4455,7 @@ int *triangle_nb, double *ax, double *ay, double *bx, double *by, double *cx, do int *t2,double *dt,int *nbtype,int *type,double *surface,double *gs, double *ks,int *error) { int i,j,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds; @@ -4554,7 +4554,7 @@ int shimatani_rect_ic(int *point_nb,double *x,double *y, double *xmi,double *xma double *gval,double *kval,int *error) { int i,j,b,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds; @@ -4613,7 +4613,7 @@ int shimatani_rect_ic(int *point_nb,double *x,double *y, double *xmi,double *xma } if(z==0) { //boucle initiale - //pour chaque espce + //pour chaque espīŋŊce for(i=0;i<*nbtype;i++) { erreur=ripley_rect(&l[i+1],xx[i],yy[i],xmi,xma,ymi,yma,t2,dt,g,k); if (erreur!=0) @@ -4730,7 +4730,7 @@ int shimatani_disq_ic(int *point_nb,double *x,double *y, double *x0,double *y0,d double *gval,double *kval,int *error) { int i,j,b,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds; @@ -4787,7 +4787,7 @@ int shimatani_disq_ic(int *point_nb,double *x,double *y, double *x0,double *y0,d if(z==0) { //boucle initiale - //pour chaque espce + //pour chaque espīŋŊce for(i=0;i<*nbtype;i++) { erreur=ripley_disq(&l[i+1],xx[i],yy[i],x0,y0,r0,t2,dt,g,k); if (erreur!=0) @@ -4899,7 +4899,7 @@ int shimatani_tr_rect_ic(int *point_nb,double *x,double *y, double *xmi,double * double *gval,double *kval,int *error) { int i,j,b,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds; @@ -4963,7 +4963,7 @@ int shimatani_tr_rect_ic(int *point_nb,double *x,double *y, double *xmi,double * } if(z==0) { //boucle initiale - //pour chaque espce + //pour chaque espīŋŊce for(i=0;i<*nbtype;i++) { erreur=ripley_tr_rect(&l[i+1],xx[i],yy[i],xmi,xma,ymi,yma,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,g,k); if (erreur!=0) @@ -5081,7 +5081,7 @@ int shimatani_tr_disq_ic(int *point_nb,double *x,double *y, double *x0,double *y double *gval,double *kval,int *error) { int i,j,b,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds; @@ -5145,7 +5145,7 @@ int shimatani_tr_disq_ic(int *point_nb,double *x,double *y, double *x0,double *y } if(z==0) { //boucle initiale - //pour chaque espce + //pour chaque espīŋŊce for(i=0;i<*nbtype;i++) { erreur=ripley_tr_disq(&l[i+1],xx[i],yy[i],x0,y0,r0,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,g,k); if (erreur!=0) @@ -5304,14 +5304,14 @@ int randomlab(double *x,double *y,int total_nb,int *type,int nb_type, double **x } /**********************************************************/ -/* Fonctions de Rao pour les semis de points multivaris */ +/* Fonctions de Rao pour les semis de points multivariīŋŊs */ /**********************************************************/ //V2 as wrapper of K12fun - 08/2013 int rao_rect(int *point_nb,double *x,double *y, double *xmi,double *xma,double *ymi,double *yma, int *t2,double *dt,int *h0,int *nbtype,int *type,double *mat,double *surface,double *HD,double *gr, double *kr,double *gs,double *ks,int *error) { int i,j,p,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds,dis; @@ -5444,7 +5444,7 @@ int rao_rect_ic(int *point_nb,double *x,double *y, double *xmi,double *xma,doubl vecalloc(&g,*t2); vecalloc(&k,*t2); - // l contient le nombre d'arbres par espce + // l contient le nombre d'arbres par espīŋŊce for(i=1;i<*nbtype+1;i++){ l[i]=0; compt[i]=0; @@ -5454,7 +5454,7 @@ int rao_rect_ic(int *point_nb,double *x,double *y, double *xmi,double *xma,doubl } } - // cration des tableaux xx et yy par espce + // crīŋŊation des tableaux xx et yy par espīŋŊce double **xx,**yy; xx=taballoca(*nbtype,l); yy=taballoca(*nbtype,l); @@ -5551,7 +5551,7 @@ int rao_rect_ic(int *point_nb,double *x,double *y, double *xmi,double *xma,doubl { gval[i]=1; kval[i]=1; } - //prparation des tableaux pour randomisation de dis (H0=2) + //prīŋŊparation des tableaux pour randomisation de dis (H0=2) int lp=0; double HDsim; int *vec, mat_size; @@ -5669,7 +5669,7 @@ int rao_disq(int *point_nb,double *x,double *y, double *x0,double *y0,double *r0 int *t2,double *dt,int *h0,int *nbtype,int *type,double *mat,double *surface,double *HD,double *gr, double *kr,double *gs,double *ks,int *error) { int i,j,p,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds,dis; @@ -5795,7 +5795,7 @@ int rao_disq_ic(int *point_nb,double *x,double *y, double *x0,double *y0,double double *gval,double *kval,int *error) { int i,j,p,b,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds,dis; @@ -6025,7 +6025,7 @@ int rao_tr_rect(int *point_nb,double *x,double *y, double *xmi,double *xma,doubl int *t2,double *dt,int *h0,int *nbtype,int *type,double *mat,double *surface,double *HD,double *gr, double *kr,double *gs,double *ks,int *error) { int i,j,p,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds,dis; @@ -6150,7 +6150,7 @@ int rao_tr_disq(int *point_nb,double *x,double *y, double *x0,double *y0, int *t2,double *dt,int *h0,int *nbtype,int *type,double *mat,double *surface,double *HD,double *gr, double *kr,double *gs,double *ks,int *error) { int i,j,p,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds,dis; @@ -6272,7 +6272,7 @@ int rao_tr_rect_ic(int *point_nb,double *x,double *y, double *xmi,double *xma,do double *gval,double *kval,int *error) { int i,j,p,b,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds,dis; @@ -6506,7 +6506,7 @@ int rao_tr_disq_ic(int *point_nb,double *x,double *y, double *x0,double *y0,doub double *gval,double *kval,int *error) { int i,j,p,b,z=0; - int *l;// contient le nombre d'arbres par espčce + int *l;// contient le nombre d'arbres par espīŋŊce int compt[*nbtype+1];// tableau de compteurs pour xx et yy vecintalloc(&l,*nbtype+1); double *ds,dis; @@ -6775,51 +6775,66 @@ int mimetic_rect(int *point_nb,double *x,double *y, double *surface,double *xmi, int compteur=0; double *l; double cout,cout_c; - double intensity=(*point_nb)/(*surface); - vecalloc(&l,*t2); + double intensity=(*point_nb)/(*surface); + vecalloc(&l,*t2); + + // Shift windows to ensure positive coordinates + double offsetx=((*xmi)<0)?(*xmi):0.; + double offsety=((*ymi)<0)?(*ymi):0.; + decalRect(*point_nb, x, y, xmi, xma, ymi, yma); //creation of a initial point pattern and cost s_alea_rect(*point_nb,x,y,*xmi,*xma,*ymi,*yma,*prec); - erreur=ripley_rect(point_nb,x,y,xmi,xma,ymi,yma,t2,dt,g,k); - if (erreur!=0) return -1; - cout=0; - for(i=0;i<*t2;i++) - { //l[i]=sqrt(k[i]/(intensity*Pi())); - l[i]=sqrt(k[i]/(intensity*Pi()))-(i+1)*(*dt); - cout+=(lobs[i]-l[i])*(lobs[i]-l[i]); - } - cost[0]=cout; - int lp=0; - if(mess!=0) - Rprintf("Simulated annealing\n"); - cout_c=0; - while(compteur<*nsimax) - { cout_c=echange_point_rect(*point_nb,x,y,*xmi,*xma,*ymi,*yma,intensity,*prec,cout,lobs,t2,dt,g,k); - if(cout==cout_c) - compteur_c++; - else - compteur_c=0; - cout=cout_c; - //Rprintf(" coût calculé : %f\n", cout); - compteur++; - cost[compteur]=cout; - if(compteur_c==*conv) - break; - if(mess!=0) { - R_FlushConsole(); - progress(compteur,&lp,*nsimax); - } - } - if(compteur==*nsimax) { - if(mess!=0) - Rprintf("Warning: failed to converge after nsimax=%d simulations",*nsimax); - r=1; - } - for(i=0;i<(*point_nb);i++) - { xx[i]=x[i]; - yy[i]=y[i]; - } - free(l); + erreur=ripley_rect(point_nb,x,y,xmi,xma,ymi,yma,t2,dt,g,k); + if (erreur!=0) return -1; + cout=0; + for(i=0;i<*t2;i++) + { //l[i]=sqrt(k[i]/(intensity*Pi())); + l[i]=sqrt(k[i]/(intensity*Pi()))-(i+1)*(*dt); + cout+=(lobs[i]-l[i])*(lobs[i]-l[i]); + } + cost[0]=cout; + int lp=0; + if(mess!=0) + Rprintf("Simulated annealing\n"); + cout_c=0; + while(compteur<*nsimax) + { + cout_c=echange_point_rect(*point_nb,x,y,*xmi,*xma,*ymi,*yma,intensity,*prec,cout,lobs,t2,dt,g,k); + if(cout==cout_c) + compteur_c++; + else + compteur_c=0; + cout=cout_c; + //Rprintf(" coīŋŊt calculīŋŊ : %f\n", cout); + compteur++; + cost[compteur]=cout; + if(compteur_c==*conv) + break; + if(mess!=0) { + R_FlushConsole(); + progress(compteur,&lp,*nsimax); + } + } + if(compteur==*nsimax) { + if(mess!=0) + Rprintf("Warning: failed to converge after nsimax=%d simulations",*nsimax); + r=1; + } + for(i=0;i<(*point_nb);i++) + { + // Shift x and y back into original windows + x[i]+=offsetx; + y[i]+=offsety; + xx[i]=x[i]; + yy[i]=y[i]; + } + // Shift windows back into original location + *xmi+=offsetx; + *xma+=offsetx; + *ymi+=offsety; + *yma+=offsety; + free(l); return r; } @@ -6843,7 +6858,7 @@ double echange_point_rect(int point_nb,double *x,double *y,double xmi,double xma for(j=0;j<4;j++) { - xcent[j]=xmi+(unif_rand()*(xr/p))*p;// coordonnées (x,y) du point tiré + xcent[j]=xmi+(unif_rand()*(xr/p))*p;// coordonnīŋŊes (x,y) du point tirīŋŊ ycent[j]=ymi+(unif_rand()*(yr/p))*p; x[num]=xcent[j]; y[num]=ycent[j]; @@ -6866,7 +6881,7 @@ double echange_point_rect(int point_nb,double *x,double *y,double xmi,double xma if(n_cout[i]<n_cout[max]) max=i; } - if(n_cout[max]<cout) // on prend le nouveau point qui minimise le coût + if(n_cout[max]<cout) // on prend le nouveau point qui minimise le coīŋŊt { x[num]=xcent[max]; y[num]=ycent[max]; @@ -6893,56 +6908,68 @@ int mimetic_disq(int *point_nb,double *x,double *y, double *surface,double *x0,d int compteur=0; double *l; double cout,cout_c; - double intensity=(*point_nb)/(*surface); + double intensity=(*point_nb)/(*surface); vecalloc(&l,*t2); + + // Shift windows to ensure positive coordinates + double offsetx=((*x0-*r0)<0)?(*x0-*r0):0.; + double offsety=((*y0-*r0)<0)?(*y0-*r0):0.; + decalCirc(*point_nb, x, y, x0, y0, *r0); //creation of a initial point pattern and cost //r=s_RegularProcess_disq(*point_nb,3,x,y,*x0,*y0,*r0,*prec); //r=s_NeymanScott_disq(5,*point_nb,3,x,y,*x0,*y0,*r0,*prec); s_alea_disq(*point_nb,x,y,*x0,*y0,*r0,*prec); erreur=ripley_disq(point_nb,x,y,x0,y0,r0,t2,dt,g,k); - if (erreur!=0) - { - return -1; - } - cout=0; - for(i=0;i<*t2;i++) - { l[i]=sqrt(k[i]/(intensity*Pi()))-(i+1)*(*dt); - cout+=(lobs[i]-l[i])*(lobs[i]-l[i]); - } - cost[0]=cout; - - int lp=0; - if(mess!=0) - Rprintf("Simulated annealing\n"); - while(compteur<*nsimax) - { cout_c=echange_point_disq(*point_nb,x,y,*x0,*y0,*r0,intensity,*prec,cout,lobs,t2,dt,g,k); - if(cout==cout_c) - { - compteur_c++; - } - else compteur_c=0; - cout=cout_c; - //Rprintf(" coût calculé : %f\n", cout); - compteur++; - cost[compteur]=cout; - if(compteur_c==*conv) - break; - if(mess!=0) { - R_FlushConsole(); - progress(compteur,&lp,*nsimax); - } - } - if(compteur==(*nsimax)) { - if(mess!=0) - Rprintf("Warning: failed to converge after nsimax=%d simulations",*nsimax); - r=1; - } - for(i=0;i<(*point_nb);i++) - { xx[i]=x[i]; - yy[i]=y[i]; - } - free(l); + if (erreur!=0) + return -1; + cout=0; + for(i=0;i<*t2;i++) + { + l[i]=sqrt(k[i]/(intensity*Pi()))-(i+1)*(*dt); + cout+=(lobs[i]-l[i])*(lobs[i]-l[i]); + } + cost[0]=cout; + + int lp=0; + if(mess!=0) + Rprintf("Simulated annealing\n"); + while(compteur<*nsimax) + { + cout_c=echange_point_disq(*point_nb,x,y,*x0,*y0,*r0,intensity,*prec,cout,lobs,t2,dt,g,k); + if(cout==cout_c) + { + compteur_c++; + } + else compteur_c=0; + cout=cout_c; + //Rprintf(" coīŋŊt calculīŋŊ : %f\n", cout); + compteur++; + cost[compteur]=cout; + if(compteur_c==*conv) + break; + if(mess!=0) { + R_FlushConsole(); + progress(compteur,&lp,*nsimax); + } + } + if(compteur==(*nsimax)) { + if(mess!=0) + Rprintf("Warning: failed to converge after nsimax=%d simulations",*nsimax); + r=1; + } + for(i=0;i<(*point_nb);i++) + { + // Shift x and y back into original windows + x[i]+=offsetx; + y[i]+=offsety; + xx[i]=x[i]; + yy[i]=y[i]; + } + // Shift windows back to original location + *x0+=offsetx; + *y0+=offsety; + free(l); return r; } @@ -6958,11 +6985,11 @@ double echange_point_disq(int point_nb,double *x,double *y,double x0,double y0,d rr=2*r0; for(j=0;j<4;j++) { - xcent[j]=x0-r0+(unif_rand()*(rr/p))*p;// coordonnées (x,y) du point tiré + xcent[j]=x0-r0+(unif_rand()*(rr/p))*p;// coordonnīŋŊes (x,y) du point tirīŋŊ ycent[j]=y0-r0+(unif_rand()*(rr/p))*p; while ((xcent[j]-x0)*(xcent[j]-x0)+(ycent[j]-y0)*(ycent[j]-y0)>=r0*r0) { - xcent[j]=x0-r0+(unif_rand()*(rr/p))*p;// coordonnées (x,y) du point tiré + xcent[j]=x0-r0+(unif_rand()*(rr/p))*p;// coordonnīŋŊes (x,y) du point tirīŋŊ ycent[j]=y0-r0+(unif_rand()*(rr/p))*p; } x[num]=xcent[j]; @@ -6992,7 +7019,7 @@ double echange_point_disq(int point_nb,double *x,double *y,double x0,double y0,d if(n_cout[i]<n_cout[max]) max=i; } - if(n_cout[max]<cout) // on prend le nouveau point qui minimise le coût + if(n_cout[max]<cout) // on prend le nouveau point qui minimise le coīŋŊt { x[num]=xcent[max]; y[num]=ycent[max]; @@ -7022,50 +7049,74 @@ int mimetic_tr_rect(int *point_nb,double *x,double *y, double *surface,double *x double cout,cout_c; double intensity=(*point_nb)/(*surface); vecalloc(&l,*t2); + + // Shift windows to ensure positive coordinates + double offsetx=((*xmi)<0)?(*xmi):0.; + double offsety=((*ymi)<0)?(*ymi):0.; + decalRectTri(*point_nb,x,y,xmi,xma,ymi,yma,*triangle_nb,ax,ay,bx,by,cx,cy); //creation of a initial point pattern and cost s_alea_tr_rect(*point_nb,x,y,*xmi,*xma,*ymi,*yma,*triangle_nb,ax,ay,bx,by,cx,cy,*prec); erreur=ripley_tr_rect(point_nb,x,y,xmi,xma,ymi,yma,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,g,k); - if (erreur!=0) - { - return -1; - } - cout=0; - for(i=0;i<*t2;i++) - { l[i]=sqrt(k[i]/(intensity*Pi()))-(i+1)*(*dt); - cout+=(lobs[i]-l[i])*(lobs[i]-l[i]); - } - cost[0]=cout; - - int lp=0; - if(mess!=0) - Rprintf("Simulated annealing\n"); - while(compteur<*nsimax) - { cout_c=echange_point_tr_rect(*point_nb,x,y,*xmi,*xma,*ymi,*yma,triangle_nb,ax,ay,bx,by,cx,cy,intensity,*prec,cout,lobs,t2,dt,g,k); - if(cout==cout_c) - compteur_c++; - else compteur_c=0; - cout=cout_c; - //Rprintf(" coût calculé : %f\n", cout); - compteur++; - cost[compteur]=cout; - if(compteur_c==*conv) - break; - if(mess!=0) { - R_FlushConsole(); - progress(compteur,&lp,*nsimax); - } - } - if(compteur==(*nsimax)) { - if(mess!=0) - Rprintf("Warning: failed to converge after nsimax=%d simulations",*nsimax); - r=1; - } - for(i=0;i<(*point_nb);i++) - { xx[i]=x[i]; - yy[i]=y[i]; - } - free(l); + if (erreur!=0) + return -1; + cout=0; + for(i=0;i<*t2;i++) + { + l[i]=sqrt(k[i]/(intensity*Pi()))-(i+1)*(*dt); + cout+=(lobs[i]-l[i])*(lobs[i]-l[i]); + } + cost[0]=cout; + + int lp=0; + if(mess!=0) + Rprintf("Simulated annealing\n"); + while(compteur<*nsimax) + { + cout_c=echange_point_tr_rect(*point_nb,x,y,*xmi,*xma,*ymi,*yma,triangle_nb,ax,ay,bx,by,cx,cy,intensity,*prec,cout,lobs,t2,dt,g,k); + if(cout==cout_c) + compteur_c++; + else compteur_c=0; + cout=cout_c; + //Rprintf(" coīŋŊt calculīŋŊ : %f\n", cout); + compteur++; + cost[compteur]=cout; + if(compteur_c==*conv) + break; + if(mess!=0) { + R_FlushConsole(); + progress(compteur,&lp,*nsimax); + } + } + if(compteur==(*nsimax)) { + if(mess!=0) + Rprintf("Warning: failed to converge after nsimax=%d simulations",*nsimax); + r=1; + } + for(i=0;i<(*point_nb);i++) + { + // Shift coordinates back into original windows + x[i]+=offsetx; + y[i]+=offsety; + xx[i]=x[i]; + yy[i]=y[i]; + } + // Shift windows back to original position + *xmi+=offsetx; + *xma+=offsetx; + *ymi+=offsety; + *yma+=offsety; + // Shift triangles back to original position + for(i=0;i<(*triangle_nb);i++) + { + ax[i]+=offsetx; + bx[i]+=offsetx; + cx[i]+=offsetx; + ay[i]+=offsety; + by[i]+=offsety; + cy[i]+=offsety; + } + free(l); return r; } @@ -7086,7 +7137,7 @@ double echange_point_tr_rect(int point_nb,double *x,double *y,double xmi,double { do { erreur_tr=0; - xcent[j]=xmi+(unif_rand()*(xr/p))*p;// coordonnées (x,y) du point tiré + xcent[j]=xmi+(unif_rand()*(xr/p))*p;// coordonnīŋŊes (x,y) du point tirīŋŊ ycent[j]=ymi+(unif_rand()*(yr/p))*p; x[num]=xcent[j]; y[num]=ycent[j]; @@ -7112,7 +7163,7 @@ double echange_point_tr_rect(int point_nb,double *x,double *y,double xmi,double if(n_cout[i]<n_cout[max]) max=i; } - if(n_cout[max]<cout) // on prend le nouveau point qui minimise le coût + if(n_cout[max]<cout) // on prend le nouveau point qui minimise le coīŋŊt { x[num]=xcent[max]; y[num]=ycent[max]; @@ -7140,56 +7191,77 @@ int mimetic_tr_disq(int *point_nb,double *x,double *y, double *surface,double *x int compteur=0; double *l; double cout,cout_c; - double intensity=(*point_nb)/(*surface); + double intensity=(*point_nb)/(*surface); vecalloc(&l,*t2); + + // Shift windows to ensure positive coordinates + double offsetx=((*x0-*r0)<0)?(*x0-*r0):0.; + double offsety=((*y0-*r0)<0)?(*y0-*r0):0.; + decalCircTri(*point_nb,x,y,x0,y0,*r0,*triangle_nb,ax,ay,bx,by,cx,cy); //creation of a initial point pattern and cost //r=s_RegularProcess_tr_disq(*point_nb,3,x,y,*x0,*y0,*r0,triangle_nb,ax,ay,bx,by,cx,cy,*prec); //r=s_NeymanScott_tr_disq(4,*point_nb,2,x,y,*x0,*y0,*r0,triangle_nb,ax,ay,bx,by,cx,cy,*prec); s_alea_tr_disq(*point_nb,x,y,*x0,*y0,*r0,*triangle_nb,ax,ay,bx,by,cx,cy,*prec); erreur=ripley_tr_disq(point_nb,x,y,x0,y0,r0,triangle_nb,ax,ay,bx,by,cx,cy,t2,dt,g,k); - if (erreur!=0) - { - return -1; - } - cout=0; - for(i=0;i<*t2;i++) - { l[i]=sqrt(k[i]/(intensity*Pi()))-(i+1)*(*dt); - cout+=(lobs[i]-l[i])*(lobs[i]-l[i]); - } - cost[0]=cout; - - int lp=0; - if(mess!=0) - Rprintf("Simulated annealing\n"); - while(compteur<*nsimax) - { cout_c=echange_point_tr_disq(*point_nb,x,y,*x0,*y0,*r0,triangle_nb,ax,ay,bx,by,cx,cy,intensity,*prec,cout,lobs,t2,dt,g,k); - if(cout==cout_c) - { - compteur_c++; - } - else compteur_c=0; - cout=cout_c; - //Rprintf(" coût calculé : %f\n", cout); - compteur++; - cost[compteur]=cout; - if(compteur_c==*conv) - break; - if(mess!=0) { - R_FlushConsole(); - progress(compteur,&lp,*nsimax); - } - } - if(compteur==(*nsimax)) { - if(mess!=0) - Rprintf("Warning: failed to converge after nsimax=%d simulations",*nsimax); - r=1; - } - for(i=0;i<(*point_nb);i++) - { xx[i]=x[i]; - yy[i]=y[i]; - } - free(l); + if (erreur!=0) + return -1; + cout=0; + for(i=0;i<*t2;i++) + { + l[i]=sqrt(k[i]/(intensity*Pi()))-(i+1)*(*dt); + cout+=(lobs[i]-l[i])*(lobs[i]-l[i]); + } + cost[0]=cout; + + int lp=0; + if(mess!=0) + Rprintf("Simulated annealing\n"); + while(compteur<*nsimax) + { + cout_c=echange_point_tr_disq(*point_nb,x,y,*x0,*y0,*r0,triangle_nb,ax,ay,bx,by,cx,cy,intensity,*prec,cout,lobs,t2,dt,g,k); + if(cout==cout_c) + compteur_c++; + else + compteur_c=0; + cout=cout_c; + //Rprintf(" coīŋŊt calculīŋŊ : %f\n", cout); + compteur++; + cost[compteur]=cout; + if(compteur_c==*conv) + break; + if(mess!=0) { + R_FlushConsole(); + progress(compteur,&lp,*nsimax); + } + } + if(compteur==(*nsimax)) { + if(mess!=0) + Rprintf("Warning: failed to converge after nsimax=%d simulations",*nsimax); + r=1; + } + for(i=0;i<(*point_nb);i++) + { + // Shift x and y back into original windows + x[i]+=offsetx; + y[i]+=offsety; + xx[i]=x[i]; + yy[i]=y[i]; + } + // Shift windows back to original location + *x0+=offsetx; + *y0+=offsety; + // Shift triangles back to original position + for(i=0;i<(*triangle_nb);i++) + { + ax[i]+=offsetx; + bx[i]+=offsetx; + cx[i]+=offsetx; + ay[i]+=offsety; + by[i]+=offsety; + cy[i]+=offsety; + } + free(l); return r; } @@ -7210,11 +7282,11 @@ double echange_point_tr_disq(int point_nb,double *x,double *y,double x0,double y do { erreur_tr=0; - xcent[j]=x0-r0+(unif_rand()*(rr/p))*p;// coordonnées (x,y) du point tiré + xcent[j]=x0-r0+(unif_rand()*(rr/p))*p;// coordonnīŋŊes (x,y) du point tirīŋŊ ycent[j]=y0-r0+(unif_rand()*(rr/p))*p; while ((xcent[j]-x0)*(xcent[j]-x0)+(ycent[j]-y0)*(ycent[j]-y0)>=r0*r0) { - xcent[j]=x0-r0+(unif_rand()*(rr/p))*p;// coordonnées (x,y) du point tiré + xcent[j]=x0-r0+(unif_rand()*(rr/p))*p;// coordonnīŋŊes (x,y) du point tirīŋŊ ycent[j]=y0-r0+(unif_rand()*(rr/p))*p; } x[num]=xcent[j]; @@ -7257,7 +7329,7 @@ double echange_point_tr_disq(int point_nb,double *x,double *y,double x0,double y if(n_cout[i]<n_cout[max]) max=i; } - if(n_cout[max]<cout) // on prend le nouveau point qui minimise le coût + if(n_cout[max]<cout) // on prend le nouveau point qui minimise le coīŋŊt { x[num]=xcent[max]; y[num]=ycent[max]; @@ -7291,7 +7363,7 @@ int shen(int *point_nb,double *x,double *y, vecalloc(&g,*t2); vecalloc(&k,*t2); - // l contient le nombre d'arbres par espce + // l contient le nombre d'arbres par espīŋŊce for(i=1;i<*nbtype+1;i++){ l[i]=0; compt[i]=0; @@ -7301,7 +7373,7 @@ int shen(int *point_nb,double *x,double *y, } } - // cration des tableaux xx et yy par espce + // crīŋŊation des tableaux xx et yy par espīŋŊce double **xx,**yy; xx=taballoca(*nbtype,l); yy=taballoca(*nbtype,l); @@ -7404,7 +7476,7 @@ int shen_ic(int *point_nb,double *x,double *y, vecalloc(&g,*t2); vecalloc(&k,*t2); - // l contient le nombre d'arbres par espce + // l contient le nombre d'arbres par espīŋŊce for(i=1;i<*nbtype+1;i++){ l[i]=0; compt[i]=0; @@ -7414,7 +7486,7 @@ int shen_ic(int *point_nb,double *x,double *y, } } - // cration des tableaux xx et yy par espce + // crīŋŊation des tableaux xx et yy par espīŋŊce double **xx,**yy; xx=taballoca(*nbtype,l); yy=taballoca(*nbtype,l); @@ -7494,7 +7566,7 @@ int shen_ic(int *point_nb,double *x,double *y, { gval[i]=1; kval[i]=1; } - //prparation des tableaux pour randomisation de dis (H0=2) + //prīŋŊparation des tableaux pour randomisation de dis (H0=2) int b,lp=0; double HDsim; int *vec, mat_size; -- GitLab