diff --git a/src/main/java/species/ReproduceAndSurviveAfterReproductionWithDiagnose.java b/src/main/java/species/ReproduceAndSurviveAfterReproductionWithDiagnose.java index 9e4ca55954bb99f666783f51c6dc73a78e696cc1..4599a0483ec229bc5809d1d219796ed1396ede8c 100644 --- a/src/main/java/species/ReproduceAndSurviveAfterReproductionWithDiagnose.java +++ b/src/main/java/species/ReproduceAndSurviveAfterReproductionWithDiagnose.java @@ -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)); }