Commit b0bb1ed0 authored by Lambert Patrick's avatar Lambert Patrick
Browse files

first attemp to calibrate GR3D avec CMA-ES

parent 560966ef
<list>
<species.DiadromousFishGroup>
<name>species A</name>
<color>
<red>255</red>
<green>0</green>
<blue>0</blue>
<alpha>255</alpha>
</color>
<dMaxDisp>300.0</dMaxDisp>
<linfVonBertForFemale>70.0</linfVonBertForFemale>
<linfVonBertForMale>70.0</linfVonBertForMale>
<lFirstMaturityForFemale>55.0</lFirstMaturityForFemale>
<lFirstMaturityForMale>40.0</lFirstMaturityForMale>
<lengthAtHatching>2</lengthAtHatching>
<nutrientRoutine>
<nutrientsOfInterest>
<string>N</string>
<string>P</string>
</nutrientsOfInterest>
<residenceTime>30.0</residenceTime>
<excretionRate class="hashtable">
<entry>
<string>P</string>
<double>2.17E-6</double>
</entry>
<entry>
<string>N</string>
<double>2.471E-5</double>
</entry>
</excretionRate>
<fishFeaturesPreSpawning class="hashtable">
<entry>
<species.DiadromousFish_-Gender>MALE</species.DiadromousFish_-Gender>
<hashtable>
<entry>
<string>bLW_Gonad</string>
<double>3.3838</double>
</entry>
<entry>
<string>aLW_Gonad</string>
<double>-8.8744</double>
</entry>
<entry>
<string>bLW</string>
<double>2.1774</double>
</entry>
<entry>
<string>aLW</string>
<double>0.27144384321211973</double>
</entry>
</hashtable>
</entry>
<entry>
<species.DiadromousFish_-Gender>FEMALE</species.DiadromousFish_-Gender>
<hashtable>
<entry>
<string>bLW_Gonad</string>
<double>2.6729</double>
</entry>
<entry>
<string>aLW_Gonad</string>
<double>-5.2425</double>
</entry>
<entry>
<string>bLW</string>
<double>3.147</double>
</entry>
<entry>
<string>aLW</string>
<double>0.007388725660209693</double>
</entry>
</hashtable>
</entry>
</fishFeaturesPreSpawning>
<fishFeaturesPostSpawning class="hashtable">
<entry>
<species.DiadromousFish_-Gender>MALE</species.DiadromousFish_-Gender>
<hashtable>
<entry>
<string>bLW_Gonad</string>
<double>3.8331</double>
</entry>
<entry>
<string>aLW_Gonad</string>
<double>-11.285</double>
</entry>
<entry>
<string>bLW</string>
<double>2.9973</double>
</entry>
<entry>
<string>aLW</string>
<double>0.010383887012522573</double>
</entry>
</hashtable>
</entry>
<entry>
<species.DiadromousFish_-Gender>FEMALE</species.DiadromousFish_-Gender>
<hashtable>
<entry>
<string>bLW_Gonad</string>
<double>2.8545</double>
</entry>
<entry>
<string>aLW_Gonad</string>
<double>-6.6234</double>
</entry>
<entry>
<string>bLW</string>
<double>2.9418</double>
</entry>
<entry>
<string>aLW</string>
<double>0.013199187556948952</double>
</entry>
</hashtable>
</entry>
</fishFeaturesPostSpawning>
<juvenileFeatures class="hashtable">
<entry>
<string>bLW</string>
<double>3.0306</double>
</entry>
<entry>
<string>aLW</string>
<double>0.006986429759979109</double>
</entry>
</juvenileFeatures>
<compoCarcassPreSpawning class="hashtable">
<entry>
<species.DiadromousFish_-Gender>MALE</species.DiadromousFish_-Gender>
<hashtable>
<entry>
<string>P</string>
<double>0.00666</double>
</entry>
<entry>
<string>N</string>
<double>0.02941</double>
</entry>
</hashtable>
</entry>
<entry>
<species.DiadromousFish_-Gender>FEMALE</species.DiadromousFish_-Gender>
<hashtable>
<entry>
<string>P</string>
<double>0.006730000000000001</double>
</entry>
<entry>
<string>N</string>
<double>0.029580000000000002</double>
</entry>
</hashtable>
</entry>
</compoCarcassPreSpawning>
<compoCarcassPostSpawning class="hashtable">
<entry>
<species.DiadromousFish_-Gender>MALE</species.DiadromousFish_-Gender>
<hashtable>
<entry>
<string>P</string>
<double>0.00961</double>
</entry>
<entry>
<string>N</string>
<double>0.0279</double>
</entry>
</hashtable>
</entry>
<entry>
<species.DiadromousFish_-Gender>FEMALE</species.DiadromousFish_-Gender>
<hashtable>
<entry>
<string>P</string>
<double>0.00997</double>
</entry>
<entry>
<string>N</string>
<double>0.03216</double>
</entry>
</hashtable>
</entry>
</compoCarcassPostSpawning>
<compoGametes class="hashtable">
<entry>
<species.DiadromousFish_-Gender>MALE</species.DiadromousFish_-Gender>
<hashtable>
<entry>
<string>P</string>
<double>0.00724</double>
</entry>
<entry>
<string>N</string>
<double>0.0325</double>
</entry>
</hashtable>
</entry>
<entry>
<species.DiadromousFish_-Gender>FEMALE</species.DiadromousFish_-Gender>
<hashtable>
<entry>
<string>P</string>
<double>0.0032</double>
</entry>
<entry>
<string>N</string>
<double>0.03242</double>
</entry>
</hashtable>
</entry>
</compoGametes>
<compoJuvenile class="hashtable">
<entry>
<string>P</string>
<double>0.00887</double>
</entry>
<entry>
<string>N</string>
<double>0.02803</double>
</entry>
</compoJuvenile>
</nutrientRoutine>
<fileNameInputForInitialObservation>data/input/reality/Obs1900.csv</fileNameInputForInitialObservation>
<centileForRange>0.95</centileForRange>
<parameterSetfileName>data/input/reality/parameterSet.csv</parameterSetfileName>
<parameterSetLine>0</parameterSetLine>
<yearOfTheUpdate>0</yearOfTheUpdate>
<basinsToUpdateFile>data/input/reality/basinsToUpdate.csv</basinsToUpdateFile>
<outputPath>data/output/</outputPath>
<processes>
<processesAtBegin>
<species.PopulateBasinNetworkWithANorthLimit>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<nbSIPerBasin>200</nbSIPerBasin>
<initialLength>20.0</initialLength>
<nbFishPerSI>2500</nbFishPerSI>
<northLimitLatitude>43.54</northLimitLatitude>
</species.PopulateBasinNetworkWithANorthLimit>
</processesAtBegin>
<processesEachStep>
<environment.InformTime>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
</environment.InformTime>
<species.PlopProcess>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<temporisation>0</temporisation>
</species.PlopProcess>
<species.Age>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
</species.Age>
<!-- 3 -->
<species.Grow>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<tempMinGrow>3.0</tempMinGrow>
<tempMaxGrow>26.0</tempMaxGrow>
<tempOptGrow>17.0</tempOptGrow>
<kOptForFemale>0.5363472</kOptForFemale>
<kOptForMale>0.3900707</kOptForMale>
<sigmaDeltaLVonBert>0.2</sigmaDeltaLVonBert>
</species.Grow>
<species.DisperseAndMigrateToRiverWithMultiNomDistriAndDeathBasin>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<alpha0Rep>-2.9</alpha0Rep>
<alpha1Rep>19.7</alpha1Rep>
<alpha3Rep>0.0</alpha3Rep>
<meanBvSurface>17351</meanBvSurface>
<standardDeviationBvSurface>35594</standardDeviationBvSurface>
<meanInterDistance>300.0</meanInterDistance>
<standardDeviationInterDistance>978.0</standardDeviationInterDistance>
<pHomingForReachEquil>0.75</pHomingForReachEquil>
<pHomingAfterEquil>0.75</pHomingAfterEquil>
<NbYearForInstallPop>0</NbYearForInstallPop>
<riverMigrationSeason>SPRING</riverMigrationSeason>
<alpha2Rep>0.0</alpha2Rep>
<meanSpawnersLengthAtRepro>45.0</meanSpawnersLengthAtRepro>
<standardDeviationOfSpawnersLengthAtRepro>2.0</standardDeviationOfSpawnersLengthAtRepro>
<weightOfDeathBasin>0.4</weightOfDeathBasin>
</species.DisperseAndMigrateToRiverWithMultiNomDistriAndDeathBasin>
<species.Survive>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<tempMinMortGenInRiv>10.0</tempMinMortGenInRiv>
<tempMaxMortGenInRiv>23.0</tempMaxMortGenInRiv>
<tempOptMortGenInRiv>20.0</tempOptMortGenInRiv>
<survivalProbOptGenInRiv>1.0</survivalProbOptGenInRiv>
<mortalityRateInRiver>0.4</mortalityRateInRiver>
<mortalityRateInSea>0.4</mortalityRateInSea>
<mortalityRateInOffshore>0.4</mortalityRateInOffshore>
</species.Survive>
<!-- 6 -->
<species.ReproduceAndSurviveAfterReproductionWithDiagnose>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<reproductionSeason>SPRING</reproductionSeason>
<tempMinRep>9.846984</tempMinRep>
<tempMaxRep>26.0</tempMaxRep>
<tempOptRep>20.0</tempOptRep>
<eta>2.4</eta>
<ratioS95__S50>1.9</ratioS95__S50>
<a>135000.0</a>
<delta__t>0.33</delta__t>
<survOptRep>0.0017</survOptRep>
<lambda>4.1E-4</lambda>
<proportionOfFemaleAtBirth>0.5</proportionOfFemaleAtBirth>
<initialLength>2.0</initialLength>
<sigmaRecruitment>0.2</sigmaRecruitment>
<survivalRateAfterReproduction>0.1</survivalRateAfterReproduction>
<maxNumberOfSuperIndividualPerReproduction>50.0
</maxNumberOfSuperIndividualPerReproduction>
<withDiagnose>false</withDiagnose>
</species.ReproduceAndSurviveAfterReproductionWithDiagnose>
<species.MigrateToSea>
<seaMigrationSeason>SUMMER</seaMigrationSeason>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
</species.MigrateToSea>
<environment.updateTemperatureInRealBasin>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<offshoreTemperature>12.0</offshoreTemperature>
</environment.updateTemperatureInRealBasin>
</processesEachStep>
<processesAtEnd>
<species.IdentifyPopulation>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<consoleDisplay>false</consoleDisplay>
<fluxesSeason>SPRING</fluxesSeason>
<years>
<long>2000</long>
<long>2100</long>
</years>
<fileNameOutput>fluxes</fileNameOutput>
</species.IdentifyPopulation>
<species.TypeTrajectoryCV>
<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
<fileNameOutput>JeuParam100_2100RCP85_A_essai</fileNameOutput>
</species.TypeTrajectoryCV>
</processesAtEnd>
</processes>
<useCemetery>false</useCemetery>
</species.DiadromousFishGroup>
</list>
......@@ -19,16 +19,274 @@
*/
package species;
import java.io.IOException;
import java.lang.reflect.InvocationTargetException;
import java.util.Arrays;
import java.util.Hashtable;
import java.util.Map;
import fr.cemagref.simaqualife.extensions.pilot.BatchRunner;
import fr.cemagref.simaqualife.pilot.Pilot;
import fr.inria.optimization.cmaes.CMAEvolutionStrategy;
import fr.inria.optimization.cmaes.fitness.IObjectiveFunction;
import miscellaneous.Duo;
import miscellaneous.ReflectUtils;
/**
*
*/
public class Calibrate {
public class Calibrate {
public Calibrate() {
}
public static void main(String[] args) {
GR3DObjeciveFunction fitfun = new GR3DObjeciveFunction(10.,10.);
double par[] = {9., 0.5, 0.6};
double x[] = fitfun.par2x(par);
System.out.println(Arrays.toString(par));
System.out.println(Arrays.toString(x));
//calibrate.valueOf(x);
System.out.println(fitfun.valueOf(x));
// new a CMA-ES and set some initial values
CMAEvolutionStrategy cma = new CMAEvolutionStrategy();
cma.setDimension(fitfun.getParameterRanges().size());
cma.setInitialX(5.);
cma.setInitialStandardDeviation(0.2);
cma.options.stopTolFun=1e-4; // function value range within iteration and of past values
// from CMAEvolutionStrategy.properties, to avoid to load the file
cma.options.stopTolFunHist = 1e-13 ; // function value range of 10+30*N/lambda past values
cma.options.stopTolX = 0.0 ; // absolute x-changes
cma.options.stopTolXFactor = 1e-11; // relative to initial stddev
cma.options.stopTolUpXFactor = 1000; // relative to initial stddev
// initialize cma and get fitness array to fill in later
double[] fitness = cma.init(); // new double[cma.parameters.getPopulationSize()];
double[][] pop;
// iteration loop
while(cma.stopConditions.getNumber() == 0) {
// --- core iteration step ---
pop = cma.samplePopulation(); // get a new population of solutions
for(int i = 0; i < pop.length; ++i) { // for each candidate solution i
// a simple way to handle constraints that define a convex feasible domain
// (like box constraints, i.e. variable boundaries) via "blind re-sampling"
// assumes that the feasible domain is convex, the optimum is
while (!fitfun.isFeasible(pop[i])) // not located on (or very close to) the domain boundary,
pop[i] = cma.resampleSingle(i); // initialX is feasible and initialStandardDeviations are
// sufficiently small to prevent quasi-infinite looping here
// compute fitness/objective value
fitness[i] = fitfun.valueOf(pop[i]); // fitfun.valueOf() is to be minimized
}
cma.updateDistribution(fitness); // pass fitness array to update search distribution
// --- end core iteration step ---
// output to files and console
cma.writeToDefaultFiles(0);
int outmod = 100;
if (cma.getCountIter() % (100*outmod) == 1)
cma.printlnAnnotation(); // might write file as well
if (cma.getCountIter() % outmod == 1)
cma.println();
//
// evaluate mean value as it is the best estimator for the optimum
cma.setFitnessOfMeanX(fitfun.valueOf(cma.getMeanX())); // updates the best ever solution
}
// final output
//cma.writeToDefaultFiles(1);
cma.println();
cma.println("Terminated due to");
for (String s : cma.stopConditions.getMessages())
cma.println(" " + s);
cma.println("Best function value " + cma.getBestFunctionValue() + " at evaluation " + cma.getBestEvaluationNumber());
}
}
/**
* objective function based on GR3D
*/
class GR3DObjeciveFunction implements IObjectiveFunction {
private double a_femaleLengthPenalty=100.;
private double a_maleLengthPenalty=100.;
private Map<String, Duo<Double, Double>> parameterRanges;
private transient Pilot pilot;
public GR3DObjeciveFunction(double a_femaleLengthPenalty, double a_maleLengthPenalty) {
loadSimulation();
this.a_femaleLengthPenalty = a_femaleLengthPenalty;
this.a_maleLengthPenalty = a_maleLengthPenalty;
parameterRanges = new Hashtable<String, Duo<Double,Double>>();
parameterRanges.put("tempMinRep", new Duo<Double, Double>(2., 15.));
parameterRanges.put("KOptFemale", new Duo<Double, Double>(0.2, .9));
parameterRanges.put("KOptMale", new Duo<Double, Double>(0.2, .9));
}
/**
*
* @return the parameterRanges
*/
public Calibrate() {
// TODO Auto-generated constructor stub
public Map<String, Duo<Double, Double>> getParameterRanges() {
return parameterRanges;
}
@Override
public double valueOf(double[] x) {
// x[0] tempMinRep
// x[1] kOptFemale
// x[2] kOptMale
double[] par = x2par(x); // in natural unit
try {
pilot.load();
ReflectUtils.setFieldValueFromPath(pilot, "aquaticWorld.aquaNismsGroupsList.0.processes.processesEachStep.6.tempMinRep", par[0]);
//System.out.println("tempMinRep: " + (double) ReflectUtils.getValueFromPath(pilot, "aquaticWorld.aquaNismsGroupsList.0.processes.processesEachStep.6.getTempMinRep"));
ReflectUtils.setFieldValueFromPath(pilot, "aquaticWorld.aquaNismsGroupsList.0.processes.processesEachStep.3.kOptForFemale", par[1]);
//System.out.println("KOptFemale: " + (double) ReflectUtils.getValueFromPath(pilot, "aquaticWorld.aquaNismsGroupsList.0.processes.processesEachStep.3.getkOptForFemale"));
ReflectUtils.setFieldValueFromPath(pilot, "aquaticWorld.aquaNismsGroupsList.0.processes.processesEachStep.3.kOptForMale",par[2]);
//System.out.println("KOptMale: " + (double) ReflectUtils.getValueFromPath(pilot, "aquaticWorld.aquaNismsGroupsList.0.processes.processesEachStep.3.getkOptForMale"));
} catch (Exception e1) {
e1.printStackTrace();
}
pilot.run();
double likelihood =0;
double femaleLengthPenalty = 0.;
double maleLengthPenalty =0.;
try {
likelihood = (double) ReflectUtils.getValueFromPath(pilot, "aquaticWorld.aquaNismsGroupsList.0.computeLikelihood");
System.out.println("likelihood: "+ likelihood);
femaleLengthPenalty = (double) ReflectUtils.getValueFromPath(pilot, "aquaticWorld.aquaNismsGroupsList.0.computeFemaleSpawnerForFirstTimeSummaryStatistic");
System.out.println("femaleLengthPenalty: "+femaleLengthPenalty);
maleLengthPenalty = (double) ReflectUtils.getValueFromPath(pilot, "aquaticWorld.aquaNismsGroupsList.0.computeMaleSpawnerForFirstTimeSummaryStatistic");
System.out.println("maleLengthPenalty: "+maleLengthPenalty);
} catch (NoSuchFieldException | IllegalArgumentException | IllegalAccessException | InvocationTargetException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
return likelihood +a_femaleLengthPenalty * femaleLengthPenalty + a_maleLengthPenalty * maleLengthPenalty ;
}
@Override
public boolean isFeasible(double[] x) {
double[] par = x2par(x);
boolean test = true;
double a, b;
a = parameterRanges.get("tempMinRep").getFirst();
b = parameterRanges.get("tempMinRep").getSecond();
if (par[0] < a | par[0]>b)
test = false;
a = parameterRanges.get("KOptFemale").getFirst();
b = parameterRanges.get("KOptFemale").getSecond();
if (par[1] < a | par[1]>b)
test = false;
a = parameterRanges.get("KOptMale").getFirst();
b = parameterRanges.get("KOptMale").getSecond();
if (par[2] < a | par[2]>b)
test = false;
return test;
}
public void loadSimulation() {
String[] args= {"-simDuration", "100", "-simBegin", "1", "-RNGStatusIndex", "1",
"-groups", "data/input/fishTryRealBV_calibration.xml",
"-env", "data/input/BNtryRealBasins.xml",
"q", "true"
};
pilot = new Pilot();
try {
// iniatialize the simulation
pilot.init();
// no display on the console ????
pilot.setQuiet();
new BatchRunner(pilot).parseArgs(args, true, true, false);
} catch (IOException | IllegalArgumentException e) {
e.printStackTrace();
}
}
/**
* @param x paremters in the range [0,10]
* @return parameter in natural unit (based on parameters range)
*/
public double[] x2par(double[] x) {
double[] naturalPar = new double[x.length]; // par in natural unit
double a, b;
a = parameterRanges.get("tempMinRep").getFirst();
b = parameterRanges.get("tempMinRep").getSecond();
naturalPar[0] = a + (b-a) * x[0] /10.;
a = parameterRanges.get("KOptFemale").getFirst();
b = parameterRanges.get("KOptFemale").getSecond();
naturalPar[1] = a + (b-a) * x[1] /10.;
a = parameterRanges.get("KOptMale").getFirst();
b = parameterRanges.get("KOptMale").getSecond();
naturalPar[2] = a + (b-a) * x[2] /10.;
return naturalPar;
}
public double[] par2x(double[] naturalPar) {
double[] x = new double[naturalPar.length];
double a, b;
a = parameterRanges.get("tempMinRep").getFirst();
b = parameterRanges.get("tempMinRep").getSecond();