Commit 8f6565b8 authored by Lambert Patrick's avatar Lambert Patrick
Browse files

cleaning commentin and bug fixed

Showing with 216 additions and 135 deletions
+216 -135
......@@ -56,6 +56,7 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
private double sigmaRecruitment = 0.3;
private double survivalRateAfterReproduction = 0.1;
private double maxNumberOfSuperIndividualPerReproduction = 50.;
private boolean withDiagnose = true;
private transient NormalGen genNormal;
......@@ -83,7 +84,7 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
List<DiadromousFish> deadFish = new ArrayList<DiadromousFish>();
for(RiverBasin riverBasin : group.getEnvironment().getRiverBasins()){
double b, c, alpha, beta, amountPerSuperIndividual , Setoile, S95, S50 ;
double b, c, alpha, beta, amountPerSuperIndividual , S95, S50 ;
double numberOfGenitors = 0.;
double numberOfAutochtones = 0.;
double numberOfSpawnerForFirstTime = 0.;
......@@ -95,50 +96,51 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
// origins of spawner during this reproduction
Map<String, Long> spawnerOriginsDuringReproduction = new HashMap<String, Long>(group.getEnvironment().getRiverBasinNames().length);
for (String basinName : group.getEnvironment().getRiverBasinNames()){
spawnerOriginsDuringReproduction.put(basinName, (long) 0);
spawnerOriginsDuringReproduction.put(basinName, 0L);
}
List<DiadromousFish> fishInBasin = riverBasin.getFishs(group);
if (fishInBasin != null){
//Calcul of b and c in stock-recruitment relationships
// --------------------------------------------------------------------------------------------------
// definition of the stock recruiment relationship
// --------------------------------------------------------------------------------------------------
// effective temperature for reproduction (priority to the ANG value)
double tempEffectRep;
if (group.getTempMinRep() == Double.NaN){
if (Double.isNaN(group.getTempMinRep())){
tempEffectRep = Miscellaneous.temperatureEffect(riverBasin.getCurrentTemperature(group.getPilot()), tempMinRep, tempOptRep, tempMaxRep);
}
else {
tempEffectRep = Miscellaneous.temperatureEffect(riverBasin.getCurrentTemperature(group.getPilot()), group.getTempMinRep(), tempOptRep, tempMaxRep);
}
// Calcul of alpha and beta of the basin
// Compute the prelimenary parameters b and c for the stock-recruitment relationship
b = (tempEffectRep == 0.) ? 0. : - Math.log(survOptRep * tempEffectRep) / delta_t;
c = lambda/riverBasin.getAccessibleSurface();
if (tempEffectRep == 0.){
b=0;
alpha=0.;
beta=0.;
}
else {
b = - Math.log(survOptRep * tempEffectRep) / delta_t;
alpha = (b * Math.exp(- b * delta_t))/(c * (1 - Math.exp(- b * delta_t)));
beta = b / (a * c * (1 - Math.exp(- b * delta_t)));
}
// Compute alpha and beta parameters of the the stock-recruitment relationship
alpha = (b * Math.exp(- b * delta_t))/(c * (1 - Math.exp(- b * delta_t)));
beta = b / (a * c * (1 - Math.exp(- b * delta_t)));
//System.out.println(a+ ", " +b + ", " + c + ", " + delta_t + "= "+ alpha);
// keep the last value of alpha (productive capacities)
riverBasin.getLastProdCapacities().push(alpha);
// Calcul of the amount per superIndividual
// Compute the amount per superIndividual
amountPerSuperIndividual = alpha / maxNumberOfSuperIndividualPerReproduction;
// Calcul of Setoile, S95 and S50
Setoile = eta * riverBasin.getAccessibleSurface();
S95 = Setoile;
// Compute the Allee effect parameters S95 and S50
S95 = eta * riverBasin.getAccessibleSurface(); // corresponds to S* in the rougier publication
S50 = S95 / ratioS95_S50;
// initilisation de la relation stock recruitment
// initilisation of the stock recruitment relationship
stockRecruitmentRelationship.init(alpha, beta, S50, S95);
// --------------------------------------------------------------------------------------------------
// calulation of the spawner number
// --------------------------------------------------------------------------------------------------
// calcul de Zcrash
//Double Zcrash = Math.log(alpha*(d-1)/(d*beta*Math.pow(d-1, 1/d)));
// age of autochnonous spawnser
Map<Integer, Long> ageOfNativeSpawners = new TreeMap<Integer, Long>();
......@@ -155,18 +157,17 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
String basinName = fish.getBirthBasin().getName();
spawnerOriginsDuringReproduction.put(basinName, spawnerOriginsDuringReproduction.get(basinName) + fish.getAmount() );
// number of autochtone and age of autochnone
// number of autochtonous fish per age
if (riverBasin == fish.getBirthBasin()){
numberOfAutochtones += fish.getAmount();
Integer age = (int) Math.floor(fish.getAge());
Integer age = (int) Math.floor(fish.getAge()); //ASK floor() or ceil()
if (ageOfNativeSpawners.containsKey(age))
ageOfNativeSpawners.put(age, ageOfNativeSpawners.get(age)+fish.getAmount());
else
ageOfNativeSpawners.put(age, fish.getAmount());
}
//System.out.println("l'�ge du poisson est :" + fish.getAge() + " et la saison est :" + Time.getSeason());
// Survive After Reproduction
// increment number of reproduction (for possible iteroparty)
fish.incNumberOfReproduction();
// survival after reproduction (semelparity or iteroparity) of SI (change the amount of the SI)
......@@ -178,76 +179,93 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
}
}
// keep the spawner number
riverBasin.setLastSpawnerNumber(numberOfGenitors);
// prepare les donn�es pour le calcul de la mortalit� associ�e aux g�niteurs autochtones
if (numberOfGenitors > 0.) {
List<Trio<Integer, Long, Long>> mortalityData= new ArrayList<Trio<Integer, Long, Long>>();
// first age
// second effective of native spwaner
// third effective of corresponding recruitment
for (Integer age : ageOfNativeSpawners.keySet()){
if (riverBasin.getLastRecruitments().getItemFromLast(age) != null){
mortalityData.add(new Trio<Integer, Long, Long>(age,
ageOfNativeSpawners.get(age),
riverBasin.getLastRecruitments().getItemFromLast(age)));
} else{
mortalityData.add(new Trio<Integer, Long, Long>(age, 0L, 0L));
// --------------------------------------------------------------------------------------------------
// Diagnose of the population dynamics in the basin
// --------------------------------------------------------------------------------------------------
if (withDiagnose==true) {
// calculate and keep the features of the stock recruitment relationships
riverBasin.setMortalityCrash(stockRecruitmentRelationship.getSigmaZcrash());
// initialise the mortality function for the autochnous spawners
// use to approximate the mortality of all the spawners to give a proxy of the Allee trap
if (numberOfGenitors > 0.) {
List<Trio<Integer, Long, Long>> mortalityData= new ArrayList<Trio<Integer, Long, Long>>();
// first age
// second effective of native spwaner
// third effective of corresponding recruitment
for (Integer age : ageOfNativeSpawners.keySet()){
if (riverBasin.getLastRecruitments().getItemFromLast(age) != null){
mortalityData.add(new Trio<Integer, Long, Long>(age,
ageOfNativeSpawners.get(age),
riverBasin.getLastRecruitments().getItemFromLast(age)));
}
else{
mortalityData.add(new Trio<Integer, Long, Long>(age, 0L, 0L));
}
}
mortalityFunction.init(mortalityData);
riverBasin.setNativeSpawnerMortality(mortalityFunction.getSigmaZ());
}
else {
riverBasin.setNativeSpawnerMortality(Double.NaN);
}
mortalityFunction.init(mortalityData);
riverBasin.setNativeSpawnerMortality(mortalityFunction.getSigmaZ());
}
else{
riverBasin.setNativeSpawnerMortality(Double.NaN);
}
riverBasin.setMortalityCrash(stockRecruitmentRelationship.getSigmaZcrash());
riverBasin.setStockTrap(stockRecruitmentRelationship.getStockTrap(riverBasin.getNativeSpawnerMortality()));
riverBasin.setLastSpawnerNumber(numberOfGenitors);
// AFFICHAGE DES RESULTATS
/*System.out.println(riverBasin.getName().toUpperCase());
//System.out.println("alpha="+alpha+ "\tbeta="+beta+"\tS50="+S50+ "\tS95="+S95);
System.out.println("\tScrash="+stockRecruitmentRelationship.getStockAtZcrash()+
"\tZcrash="+ stockRecruitmentRelationship.getSigmaZcrash() +
"\tZ="+ riverBasin.getNativeSpawnerMortality());
System.out.println("\tStrap="+stockRecruitmentRelationship.getStockTrap(riverBasin.getNativeSpawnerMortality())+
"\tStotal="+numberOfGenitors+"\tSautochthonous="+
spawnerOriginsDuringReproduction.get(riverBasin.getName()));
String diagnose;
if (Double.isNaN(riverBasin.getNativeSpawnerMortality()))
diagnose="noSense";
else {
double stockTrap=stockRecruitmentRelationship.getStockTrap(riverBasin.getNativeSpawnerMortality());
if (riverBasin.getNativeSpawnerMortality()>stockRecruitmentRelationship.getSigmaZcrash())
diagnose="overZcrash";
riverBasin.setStockTrap(stockRecruitmentRelationship.getStockTrap(riverBasin.getNativeSpawnerMortality()));
System.out.println(riverBasin.getName().toUpperCase());
//System.out.println("alpha="+alpha+ "\tbeta="+beta+"\tS50="+S50+ "\tS95="+S95);
System.out.println("\tScrash="+stockRecruitmentRelationship.getStockAtZcrash()+
"\tZcrash="+ stockRecruitmentRelationship.getSigmaZcrash() +
"\tZ="+ riverBasin.getNativeSpawnerMortality());
System.out.println("\tStrap="+stockRecruitmentRelationship.getStockTrap(riverBasin.getNativeSpawnerMortality())+
"\tStotal="+numberOfGenitors+"\tSautochthonous="+
spawnerOriginsDuringReproduction.get(riverBasin.getName()));
String message;
if (Double.isNaN(riverBasin.getNativeSpawnerMortality()))
message="noSense";
else {
if (numberOfGenitors < stockTrap)
diagnose = "inTrapWithStrayers";
double stockTrap=stockRecruitmentRelationship.getStockTrap(riverBasin.getNativeSpawnerMortality());
if (riverBasin.getNativeSpawnerMortality()>stockRecruitmentRelationship.getSigmaZcrash())
message="overZcrash";
else {
if (spawnerOriginsDuringReproduction.get(riverBasin.getName()) < stockTrap)
diagnose = "inTrapWithOnlyNatives";
else
diagnose = "sustain";
if (numberOfGenitors < stockTrap)
message = "inTrapWithStrayers";
else {
if (spawnerOriginsDuringReproduction.get(riverBasin.getName()) < stockTrap)
message = "inTrapWithOnlyNatives";
else
message = "sustain";
}
}
}
System.out.println("\t"+message);
}
System.out.println("\t"+diagnose);*/
// --------------------------------------------------------------------------------------------------
// Reproduction process (number of recruits)
// --------------------------------------------------------------------------------------------------
if (numberOfGenitors > 0.) {
//BH Stock-Recruitment relationship with logistic depensation
double meanNumberOfRecruit = stockRecruitmentRelationship.getRecruitment(numberOfGenitors);
muRecruitment = Math.log(meanNumberOfRecruit) - (Math.pow(sigmaRecruitment,2))/2;
// lognormal random draw
muRecruitment = Math.log(meanNumberOfRecruit) - (Math.pow(sigmaRecruitment,2))/2;
long numberOfRecruit = Math.round(Math.exp(genNormal.nextDouble()*sigmaRecruitment + muRecruitment));
//System.out.println(group.getPilot().getCurrentTime()+" "+Time.getSeason(group.getPilot())+" "+ riverBasin.getName()+": " + numberOfGenitors + " spwaners \tgive "+ numberOfRecruit + " recruits");
// keep last % of autochtone
riverBasin.getLastPercentagesOfAutochtones().push(numberOfAutochtones * 100 / numberOfGenitors);
// keep the number of spawners for the firt time in the basin
if (numberOfSpawnerForFirstTime>0){
riverBasin.getSpawnersForFirstTimeMeanAges().push(spawnersForFirstTimeAgesSum/numberOfSpawnerForFirstTime);
}else{
......@@ -258,31 +276,34 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
//System.out.println("nb recruit in basin " + riverBasin.getName() + " : " + numberOfRecruit);
// Creation of new superFish
if(numberOfRecruit > 0){
// stock the first year when recruitment is non nul
if(riverBasin.getYearOfFirstNonNulRep() == 0){
riverBasin.setYearOfFirstNonNulRep(Time.getYear(group.getPilot()));
}
if (numberOfRecruit > 0){
// features of the super individuals
int numberOfsuperIndividual = Math.max(1,
(int) Math.round(numberOfRecruit / amountPerSuperIndividual));
long effectiveAmount = numberOfRecruit / numberOfsuperIndividual;
// System.out.println(numberOfRecruit + " / " + amountPerSuperIndividual +" = " +numberOfsuperIndividual);
//System.out.println(numberOfRecruit + " / " + numberOfsuperIndividual +" = " +effectiveAmount);
for (int i=0; i<numberOfsuperIndividual; i++){
for (int i = 0; i < numberOfsuperIndividual; i++){
group.addAquaNism(new DiadromousFish(group.getPilot(), riverBasin, initialLength, effectiveAmount, Gender.UNDIFFERENCIED));
}
// stock the first year when recruitment is non nul
if (riverBasin.getYearOfFirstNonNulRep() == 0){
riverBasin.setYearOfFirstNonNulRep(Time.getYear(group.getPilot()));
}
// keep the last recruitments in the stack
riverBasin.getLastRecruitmentExpectations().push(Math.round(meanNumberOfRecruit));
riverBasin.getLastRecruitments().push(numberOfsuperIndividual * effectiveAmount); // on remplit la pile qui permet de stocker un nombre fix� de derniers recrutement
riverBasin.getLastRecruitments().push(numberOfsuperIndividual * effectiveAmount);
riverBasin.getLastRecsOverProdCaps().push(((double) riverBasin.getLastRecruitments().getLastItem())/riverBasin.getLastProdCapacities().getLastItem());
// keep the no null recruitment
if (numberOfAutochtones > 0){
riverBasin.getNumberOfNonNulRecruitmentForFinalProbOfPres().push(1.0);
riverBasin.getNumberOfNonNulRecruitmentDuringLastYears().push(1.0);
riverBasin.getNumberOfNonNulRecruitmentForFinalProbOfPres().push(1.);
riverBasin.getNumberOfNonNulRecruitmentDuringLastYears().push(1.);
}else{
riverBasin.getNumberOfNonNulRecruitmentForFinalProbOfPres().push(0.0);
riverBasin.getNumberOfNonNulRecruitmentDuringLastYears().push(0.0);
riverBasin.getNumberOfNonNulRecruitmentForFinalProbOfPres().push(0.);
riverBasin.getNumberOfNonNulRecruitmentDuringLastYears().push(0.);
}
}
......@@ -293,7 +314,7 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
riverBasin.getLastRecruitments().push((long) 0);
riverBasin.getLastRecsOverProdCaps().push(0.);
riverBasin.getNumberOfNonNulRecruitmentDuringLastYears().push(0.);
riverBasin.getNumberOfNonNulRecruitmentForFinalProbOfPres().push(0.0);
riverBasin.getNumberOfNonNulRecruitmentForFinalProbOfPres().push(0.);
}
}
else {
......@@ -304,10 +325,11 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
riverBasin.getLastRecsOverProdCaps().push(0.);
riverBasin.getLastPercentagesOfAutochtones().push(0.);
riverBasin.getNumberOfNonNulRecruitmentDuringLastYears().push(0.);
riverBasin.getNumberOfNonNulRecruitmentForFinalProbOfPres().push(0.0);
riverBasin.getNumberOfNonNulRecruitmentForFinalProbOfPres().push(0.);
}
// --------------------------------------------------------------------------------------------------
// Remove deadfish
// --------------------------------------------------------------------------------------------------
for (DiadromousFish fish : deadFish){
group.removeAquaNism(fish);
}
......@@ -322,45 +344,101 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
riverBasin.getSpawnerOrigins().push(spawnerOriginsDuringReproduction);
//System.out.println(" AFTER " +riverBasin.getSpawnerOrigins().keySet());
}
// on met � jour les observeurs
// --------------------------------------------------------------------------------------------------
// update the observers
// --------------------------------------------------------------------------------------------------
for (RiverBasin riverBasin : group.getEnvironment().getRiverBasins()){
riverBasin.getCobservable().fireChanges(riverBasin, pilot.getCurrentTime());
}
}
}
/**
* Berverton and Holt stock-recruitment relationship with an Allee effect simulated with a logistic function
*/
class StockRecruitmentRelationship implements UnivariateFunction{
/**
* alpha of the stock-recruitment relationship
* R = beta Seff / (beta + Seff) with Seff the stock that effectivly participates to the reproduction
*
* @unit # of individuals
*/
double alpha;
/**
* * beta of the stock-recruitment relationship
* R = alpha * Seff / (beta + Seff) with Seff the stock that effectivly participates to the reproduction
* @unit
*/
double beta;
/**
* the value of the stock for which 50% partipate to the reproduction
* @unit # of individuals
*/
double S50;
/**
* the value of the stock for which 95% partipate to the reproduction
* @unit # of individuals
*/
double S95;
double sigmaZ; // mortality
/**
* to simplify the calculation
* @unit
*/
transient double log19;
/**
* the value of the stock for which 50% partipate to the reproduction
* @unit # of individuals
*/
double sigmaZ; //
public void init(double alpha, double beta, double S50, double S95) {
this.alpha = alpha;
this.beta = beta;
this.S50 = S50;
this.S95 = S95;
log19 = Math.log(19) ;
}
public double getEffectiveStock (double stock) {
return stock / (1 + Math.exp(- log19 * (stock - S50) / (S95 - S50)));
}
public double getRecruitment(double stock){
//BH Stock-Recruitment relationship with logistic depensation
double meanNumberOfRecruit = 0.;
double effectiveStock = getEffectiveStock(stock);
if (stock >0)
meanNumberOfRecruit= Math.round((alpha * stock * (1 / (1 + Math.exp(- Math.log(19)*((stock - S50) / (S95 - S50)))))) /
(beta + stock * (1 / (1 + Math.exp(- Math.log(19)*((stock - S50) / (S95 - S50)))))));
meanNumberOfRecruit = Math.round(alpha * effectiveStock) /(beta + effectiveStock );
return meanNumberOfRecruit;
}
/**
* the stock that corresponds to the intersection between SR relationship and tahe tangent that pass through the origin
* @return the stock at
*/
public double getStockAtZcrash(){
if (beta !=0)
return(S50 + (S95 - S50) * Math.log(beta * Math.log(19) / (S95-S50)) / Math.log(19));
return(S50 + (S95 - S50) * Math.log(beta * log19 / (S95-S50)) / log19);
else
return Double.NaN;
}
/**
* the crash mortality
* (corresponds the slope of the tangent that pass through the origin)
* @return the crash mortality ( year-1)
*/
public double getSigmaZcrash(){
double stockAtZcrash= getStockAtZcrash();
if (!Double.isNaN(stockAtZcrash))
......@@ -370,13 +448,23 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
}
/**
* the objective function uses to calculate the depensation trap (Allee trap)
*/
@Override
public double value(double S) {
double res=getRecruitment(S)-S*Math.exp(sigmaZ);
public double value(double stock) {
double res=getRecruitment(stock) - stock * Math.exp(sigmaZ);
return res*res;
}
private double getPoint(double sigmaZ){
/**
* calculate the stock correspondinf to the depensation trap
* corresponds to intersection between mortality rate and the stock-recruitment relationship
* @param sigmaZ the total mortality coefficient
* @return
*/
private double getStockTrap(double sigmaZ){
if (!Double.isNaN(sigmaZ)){
this.sigmaZ=sigmaZ;
BrentOptimizer optimizer = new BrentOptimizer(1e-6, 1e-12);
......@@ -385,25 +473,23 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
new MaxEval(100),
GoalType.MINIMIZE,
new SearchInterval(0, getStockAtZcrash()));
this.sigmaZ = Double.NaN;
this.sigmaZ = Double.NaN; //WHY
return solution.getPoint();
}
else
return Double.NaN;
}
public double getStockTrap(double sigmaZ) {
return getPoint(sigmaZ);
}
}
/**
* mortatiity function for stock with dirrenet ages
*/
class MortalityFunction implements UnivariateFunction {
// first age
// second effective of native spwaner
// third effective of corresponding recruitment
List<Trio<Integer, Long, Long>> data;
// first: age
// second: effective of native spwaner
// third: effective of corresponding recruitment
List<Trio<Integer, Long, Long>> data; //WHY age as integer
public void init(List<Trio<Integer, Long, Long>> data){
this.data = data;
......@@ -419,50 +505,45 @@ public class ReproduceAndSurviveAfterReproductionWithDiagnose extends AquaNismsG
return res/effTotal;
}
@Override
public double value(double Z) {
double res=0.;
for(Trio<Integer, Long, Long> trio : data){
res += (((double) trio.getSecond())/((double) trio.getThird())) * Math.exp(Z * (double) trio.getFirst());
}
return Math.pow(res-1., 2.);
return (res-1) * (res-1); //WHY -1
}
private double getPoint(){
if (data.isEmpty()){
return Double.NaN;
}
else {
/**
* calculate by optimsation of the total mortality coefficient over the livespan
* @return
*/
private double getSigmaZ2(){
double Z = Double.NaN;
if (!data.isEmpty()){
BrentOptimizer optimizer = new BrentOptimizer(1e-6, 1e-12);
UnivariatePointValuePair solution =
optimizer.optimize(new UnivariateObjectiveFunction(this),
new MaxEval(100),
GoalType.MINIMIZE,
new SearchInterval(0, 10));
return solution.getPoint();
Z= solution.getPoint();
}
return Z * getMeanAge();
}
// return sigmaZ = Z*meanAge
public double getSigmaZ2(){
/*double Z = getPoint();
double alpha; // % of maturation for a given age
double sum=0;
double sum1=0;
for(Trio<Integer, Long, Long> trio : data){
alpha = ((double)trio.getSecond()) / (((double)trio.getThird())*Math.exp(-((double) trio.getFirst())*Z));
sum += alpha * Math.exp(-((double) trio.getFirst())*Z);
sum1 += ((double)trio.getSecond()) / ((double)trio.getThird());
}
System.out.println(getPoint()*getMeanAge()+" <> "+ Math.log(sum)+ " <> " + Math.log(sum1));*/
return getPoint()*getMeanAge();
}
/**
* simple approximation of total mortality coefficient over the lifespan
* @return
*/
public double getSigmaZ(){
double sum=0;
for(Trio<Integer, Long, Long> trio : data){
sum += ((double)trio.getSecond()) / ((double)trio.getThird());
sum += ((double)trio.getSecond()) / ((double)trio.getThird());
}
return (- Math.log(sum));
}
......
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