Commit 4356bf84 authored by Julien Papaïx's avatar Julien Papaïx
Browse files

Imrpove demo

parent 7dcca25a
......@@ -83,7 +83,7 @@ results <- st_as_sf(landscape)
st_write(habitat, paste0(filename,".gpkg"), "habitat",layer_options="OVERWRITE=yes")
st_write(results, paste0(filename,".gpkg"), "results", update = TRUE,layer_options="OVERWRITE=yes")
return(list(shapefilename=paste0(filename,".gpkg"),layername_hab="habitat",layername_res="results",rotation=rotation$habitat,Cmax0=Cmax0,Cmax1=Cmax1,propRR=propRR,strat=strat))
return(list(shapefilename=paste0(filename,".gpkg"),layername_hab="habitat",layername_res="results",rotation=as.integer(rotation$habitat),Cmax0=Cmax0,Cmax1=Cmax1,propRR=propRR,strat=strat))
}
......
......@@ -40,7 +40,9 @@ landscape <- AgriLand(landscapeTEST,filename="paysageTEST",propSR,isolSR,propRR,
#plot(resultat)
nTSpY = 120
fichierP <- "/home/julien/Documents/Documents/modelisation/CSIRO_INRA/landsepi/data/TEST_dispP.txt"
fichierH <- "/home/julien/Documents/Documents/modelisation/CSIRO_INRA/landsepi/data/TEST_dispH.txt"
modeleLandsLPI(times = list(nYears=nYears,nTSpY=nTSpY) , landscape, dispersion=list(dispP="fichier"), inits=list(C_0=0.1, PI0=5e-4), val_seed=12345, host=list(Nhote=3))
modeleLandsLPI(times = list(nYears=nYears,nTSpY=nTSpY) , landscape, dispersion=list(dispP=fichierP,dispH=fichierH), inits=list(C_0=0.1, PI0=5e-4), val_seed=12345, host=list(Nhote=3))
}
......@@ -511,7 +511,6 @@ void dynepi(const gsl_rng *gen, int nYears, int nTSpY, int Npoly, int Nhote, int
L2I = init_i4(Npoly,Npatho,Nhote,nTSpY);
I2R = init_i4(Npoly,Npatho,Nhote,nTSpY);
eqIsurv = init_i3(Npoly, Npatho, Nhote);
/* Initialisation (t=0) */
init_HHjuvSLIR(Npoly, Nhote, Npatho, H, Hjuv, S, L, I, R);
init_L2I2R(Npoly, Npatho, Nhote, nTSpY, L2I, I2R);
......@@ -611,13 +610,23 @@ void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersio
CharacterVector shapefilename = Rcpp::as<CharacterVector>(landscape["shapefilename"]);
CharacterVector layername_hab = Rcpp::as<CharacterVector>(landscape["layername_hab"]);
char * strat = (char*)Rcpp::as<std::string>(Rcpp::as<CharacterVector>(landscape["strat"])).c_str();
int * rotation = Rcpp::as<std::vector<int>>(landscape["rotation"]).data();
IntegerVector rotation_tmp = Rcpp::as<IntegerVector>(landscape["rotation"]);
double propRR = Rcpp::as<double>(landscape["propRR"]);
double CMAX0 = Rcpp::as<double>(landscape["Cmax0"]);
double CMAX1 = Rcpp::as<double>(landscape["Cmax1"]);
//fprintf(stderr,"\ntype rotation: %s\n",typeid(rotation[0]).name());
int * rotation=init_i1(nYears+1);
IntegerVector::iterator it;
int roti;
for (roti=0, it = rotation_tmp.begin() ; it != rotation_tmp.end(); ++it, roti++){
rotation[roti]=*it;
}
CharacterVector nomfile_disppatho = Rcpp::as<CharacterVector>(dispersion["dispP"]);
CharacterVector nomfile_disphote = Rcpp::as<CharacterVector>(landscape["dispH"]);
CharacterVector nomfile_disphote = Rcpp::as<CharacterVector>(dispersion["dispH"]);
double C_0 = Rcpp::as<double>(inits["C_0"]);
double PI0 = Rcpp::as<double>(inits["PI0"]);
......@@ -652,7 +661,7 @@ void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersio
{
OGRGeometry *poGeometry;
poGeometry = poFeature->GetGeometryRef();
if( poGeometry != NULL && wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon ) NPoly++;
if( poGeometry != NULL && (wkbFlatten(poGeometry->getGeometryType()) == wkbMultiPolygon || wkbFlatten(poGeometry->getGeometryType()) == wkbPolygon)) NPoly++;
else cerr<<"WARNING unknow geometry type : "<<poGeometry->getGeometryType()<<endl;
}
......@@ -663,7 +672,7 @@ void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersio
poLayer->ResetReading();
while( (poFeature = poLayer->GetNextFeature()) != NULL ){
habitat[0][ipoly] = poFeature->GetFieldAsInteger(0);
habitat[1][ipoly] = poFeature->GetFieldAsInteger(0);
habitat[1][ipoly] = poFeature->GetFieldAsInteger(1);
area[ipoly] = poFeature->GetFieldAsInteger(2);
ipoly = ipoly+1;
}
......@@ -849,10 +858,10 @@ void modeleLandsLPI(Rcpp::List times, Rcpp::List landscape, Rcpp::List dispersio
/*--------------------------------------*/
FILE *fP;
fP = fopen("parameters.txt","w");
print_param(fP, seed, (char*)Rcpp::as<std::string>(nomfile_disppatho).c_str(), (char*)Rcpp::as<std::string>(nomfile_disphote).c_str(), "nomfile_habitat1", "nomfile_habitat2", nYears, nTSpY, pSurv, eff, Tlat, repro, Tspo, NPoly, NPoly, Nhote, Npatho, Naggr, C0, pI0, Cmax, delta, reproH, area, habitat, rotation, mort, strat, propRR, resistance, adaptation, MGeff, QReff, taumut, costInfect, costAggr, beta);
print_param(fP, seed, (char*)Rcpp::as<std::string>(nomfile_disppatho).c_str(), (char*)Rcpp::as<std::string>(nomfile_disphote).c_str(), "nomfile_habitat1", "nomfile_habitat2", nYears, nTSpY, pSurv, eff, Tlat, repro, Tspo, NPoly, NPoly, Nhote, Npatho, Naggr, C0, pI0, Cmax, delta, reproH, area, habitat, rotation, mort, strat, propRR, resistance, adaptation, MGeff, QReff, taumut, costInfect, costAggr, beta);
print_param(stdout, seed, (char*)Rcpp::as<std::string>(nomfile_disppatho).c_str(), (char*)Rcpp::as<std::string>(nomfile_disphote).c_str(), "nomfile_habitat1", "nomfile_habitat2", nYears, nTSpY, pSurv, eff, Tlat, repro, Tspo, NPoly, NPoly, Nhote, Npatho, Naggr, C0, pI0, Cmax, delta, reproH, area, habitat, rotation, mort, strat, propRR, resistance, adaptation, MGeff, QReff, taumut, costInfect, costAggr, beta);
fclose(fP);
printf("Pathogen dispersal matrix [poly][poly]\n");
for (i=0 ; i<NPoly ; i++){
for (j=0 ; j<NPoly ; j++)
......
......@@ -37,8 +37,8 @@
#define PRINTON 0
#define NAGGR 6
#define STRAT "RO"
#define PROPRR 4/6
//#define STRAT "RO"
//#define PROPRR 4/6
#define COSTINFECT 0.75
#define COSTAGGR 0.5
#define TAUMUT 1e-7
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment