Merge branch 'SpawnerRunAnalysis' of gitlab-ssh.irstea.fr:SimAquaLife/GR3D...

Merge branch 'SpawnerRunAnalysis' of gitlab-ssh.irstea.fr:SimAquaLife/GR3D into exploration_GR3D_process

.gitignore

@@ -117,11 +117,10 @@ org.*
 # R: files
 *.Rproj
 *.Rhistory
-*.RData
+#*.RData
 *.Rproj.user
 *.~lock.*
 
 # qgis file
 *.qix
-
 .Rproj.user

data/input/northeastamerica/fishRIOBasin_Sapidissima_Rjava.xml

@@ -202,8 +202,8 @@
 		<lengthAtHatching>2.8</lengthAtHatching>
 		<linfVonBertForFemale>76.0</linfVonBertForFemale>
 		<linfVonBertForMale>76.0</linfVonBertForMale>
-		<lFirstMaturityForFemale>44.8</lFirstMaturityForFemale>
-		<lFirstMaturityForMale>40.2</lFirstMaturityForMale>
+		<lFirstMaturityForFemale>47.26</lFirstMaturityForFemale>
+		<lFirstMaturityForMale>40.6</lFirstMaturityForMale>
 		<processes>
 			<processesAtBegin>
 
@@ -237,10 +237,10 @@
 					<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
 					<tempMinGrow>1.7</tempMinGrow>
 					<tempMaxGrow>27.9</tempMaxGrow>
-					<tempOptGrow>4.5</tempOptGrow>
-					<kOptForFemale>0.57</kOptForFemale>
-					<kOptForMale>0.61</kOptForMale>
-					<sigmaDeltaLVonBert>0.4</sigmaDeltaLVonBert>
+					<tempOptGrow>6.2</tempOptGrow>
+					<kOptForFemale>0.25</kOptForFemale>
+					<kOptForMale>0.23</kOptForMale>
+					<sigmaDeltaLVonBert>0.2</sigmaDeltaLVonBert>
 				</species.Grow>
 
 				<species.MigrateFromOffshoreToInshore>
@@ -278,6 +278,12 @@
 					<mortalityRateInOffshore>0.4</mortalityRateInOffshore>
 				</species.Survive>
 
+				<species.SurviveWhenBecomeOld>
+					<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
+					<ageMax>12.0</ageMax>
+					<verbose>true</verbose>
+				</species.SurviveWhenBecomeOld>
+
 				<!--<species.WriteEffectiveAndBiomassImportFluxes> <synchronisationMode>ASYNCHRONOUS</synchronisationMode> <exportSeason>SPRING</exportSeason> 
 					<fileNameOutput>effectiveBiomassFluxesBeforeReproduction</fileNameOutput> </species.WriteEffectiveAndBiomassImportFluxes> -->
 
@@ -319,13 +325,14 @@
 					<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
 					<afterReproductionSeason>SUMMER</afterReproductionSeason>
 					<maximalSurvivalRate>0.18</maximalSurvivalRate>
-					<temperatureEffectSurvivalAfterReproduction class="temperatureEffect.LogitEffect">
-					<Tref>19.9</Tref>
-					<alpha>-0.58</alpha>
+					<temperatureEffectSurvivalAfterReproduction
+						class="temperatureEffect.LogitEffect">
+						<Tref>19.9</Tref>
+						<alpha>-0.58</alpha>
 					</temperatureEffectSurvivalAfterReproduction>
-<!-- 					<temperatureEffectSurvivalAfterReproduction -->
-<!-- 						class="temperatureEffect.NoEffect"> -->
-<!-- 					</temperatureEffectSurvivalAfterReproduction> -->
+					<!-- <temperatureEffectSurvivalAfterReproduction -->
+					<!-- class="temperatureEffect.NoEffect"> -->
+					<!-- </temperatureEffectSurvivalAfterReproduction> -->
 				</species.SurviveAfterReproduction>
 
 				<species.MigrateFromRiverToInshore>

exploration/GR3D_Rdescription/GR3Dfunction.R

@@ -39,6 +39,7 @@ vonBertalanffyInverse = function(L, L0, Linf, K){
 
 # von Bertalanffy increment
 # pas cohérent avec la temperature effet sur le coeff de croissance mais ca marche
+
 vonBertalanffyIncrement = function(nStep, L0, Linf, K, timeStepDuration, sigma, withTempEffect=FALSE, TrefAtSea = c(9.876946, 13.489854, 15.891487, 11.554104) ){
   tempEffect = temperatureEffect( TrefAtSea , 3, 17, 26)
   L = matrix(nrow = nStep + 1)
@@ -54,7 +55,7 @@ vonBertalanffyIncrement = function(nStep, L0, Linf, K, timeStepDuration, sigma,
   return(L)
 }
 
-vonBertalanffyWithNextIncrement = function(L, L0, Linf, K, timeStepDuration, sigma, tempEffect ){
+vonBertalanffyWithNextIncrement = function(L,  Linf, K, timeStepDuration, sigma, tempEffect ){
   if (sigma == 0) {
     mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
     increment = exp(mu)
@@ -67,6 +68,128 @@ vonBertalanffyWithNextIncrement = function(L, L0, Linf, K, timeStepDuration, sig
   return(L)
 }
 
+# generate a cohort of length trajectories
+computeMultipleLengthTrajectories = function(temperaturePattern, Nind = 10, growPar){
+  ages = temperaturePattern$age
+  
+  res <- expand_grid(ind = seq.int(Nind),  age = ages) %>% 
+    inner_join(temperaturePattern, by = 'age') %>% 
+    arrange(age, ind) %>% 
+    mutate(temperatureEffect = temperatureEffect(temperature, 
+                                                 growPar$tempMinGrow, 
+                                                 growPar$tempOptGrow, 
+                                                 growPar$tempMaxGrow), 
+           L = if_else(age == 0, growPar$lengthAtHatching, 0))
+  
+  for (i in 2:length(ages)) {
+    previousAge  = ages[i - 1]
+    currentAge = ages[i]
+    tempEffect = res %>% filter(age == currentAge) %>% select(temperatureEffect) %>% unlist(use.names = FALSE)
+    previousL <- res %>% filter(age == previousAge) %>% select(L) %>% unlist(use.names = FALSE)
+    
+    currentL <- vonBertalanffyWithNextIncrement(L = previousL, 
+                                                Linf = growPar$Linf, 
+                                                K = growPar$kOpt,
+                                                timeStepDuration = currentAge - previousAge, 
+                                                sigma = growPar$sigmaDeltaLVonBert, 
+                                                tempEffect = tempEffect)
+    res = res %>% mutate(L =replace(L, age == ages[i], currentL) ) 
+  }
+  return(res)
+}
+
+computeMultipleLengthTrajectoriesWithRandomSeed = function(temperaturePattern, 
+                                                           Nind = 10, 
+                                                           growPar, 
+                                                           RNGseed =1){
+  set.seed(RNGseed)
+  ages = temperaturePattern %>% 
+    distinct(age) %>% 
+    unlist(use.names = FALSE)
+  
+  res <- expand_grid(ind = seq.int(Nind),  age = ages) %>% 
+    mutate(random = rnorm(Nind * length(ages)))  %>% 
+    inner_join(temperaturePattern, by = 'age') %>% 
+    arrange(age, ind) %>% 
+    mutate(temperatureEffect = temperatureEffect(temperature, 
+                                                 growPar$tempMinGrow, 
+                                                 growPar$tempOptGrow, 
+                                                 growPar$tempMaxGrow), 
+           L = if_else(age == 0, growPar$lengthAtHatching, 0))
+  
+  for (i in 2:length(ages)) {
+    previousAge  = ages[i - 1]
+    currentAge = ages[i]
+    tempEffect <-  res %>% filter(age == currentAge) %>% 
+      select(temperatureEffect) %>% unlist(use.names = FALSE)
+    
+    previousL <- res %>% filter(age == previousAge) %>% 
+      select(L) %>% unlist(use.names = FALSE)
+    rnd <- res %>% filter(age == currentAge) %>% select(random) %>% unlist(use.names = FALSE)
+    currentL <- vonBertalanffyWithRandomVector(L = previousL, 
+                                               Linf = growPar$Linf, 
+                                               K = growPar$kOpt,
+                                               timeStepDuration = currentAge - previousAge, 
+                                               randomVector = rnd,
+                                               sigma = growPar$sigmaDeltaLVonBert, 
+                                               tempEffect = tempEffect)
+    res = res %>% mutate(L = replace(L, age == ages[i], currentL) ) 
+  }
+  return(res)
+}
+
+#computeMultipleLengthTrajectoriesWithRandomSeed(temperaturePattern, Nind = 10, 
+#growPar = growParUnisex,
+# RNGseed = 1)
+
+vonBertalanffyWithRandomVector = function(L, Linf, K, timeStepDuration, randomVector, sigma, tempEffect ){
+  if (sigma == 0) {
+    mu = if_else(L < Linf, log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))),-Inf)
+    #mu =  log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration)))
+    increment = exp(mu)
+  }
+  else {
+    mu = if_else(L < Linf, log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) - (sigma * sigma) / 2, -Inf)
+    #      mu = log((Linf - L) * (1 - exp(-K * tempEffect * timeStepDuration))) - (sigma * sigma) / 2
+    increment = exp(randomVector * sigma + mu)
+  }
+  L = pmin(Linf, L + increment) 
+  return(L)
+}
+
+computeOgive = function(lengthTrajectories, lengthAtMaturity){
+  ogive <- lengthTrajectories %>%
+    group_by(age) %>%
+    summarise(nTotal = n()) %>%
+    left_join(lengthTrajectories %>%  group_by(age) %>%
+                filter(L >= lengthAtMaturity) %>%
+                summarise(taller = n(), .groups = 'drop'),
+              by = 'age') %>%
+    replace_na(list(taller = 0)) %>%
+    mutate(mature = c(0,diff(taller)), 
+           immature = nTotal - cumsum(mature),
+           p = if_else(mature + immature > 0 , mature / (mature + immature), 1)) %>% 
+    select(age, p)
+  return(ogive)
+}
+
+# deprecated
+computeOgive3 = function(lengthTrajectories, lengthAtMaturity) {
+  ogive <-
+    lengthTrajectories %>% 
+    group_by(age) %>% 
+    summarise(nTotal=n()) %>% 
+    left_join(lengthTrajectories %>%  group_by(age) %>% 
+                filter(L > lengthAtMaturity) %>% 
+                summarise(mature = n(), .groups = 'drop'), 
+              by = 'age') %>% 
+    replace_na(list(mature = 0)) %>% 
+    mutate(p =  mature/nTotal) %>% 
+    select(age, p)
+  
+  return(ogive)
+} 
+
 # -------------------------------------------------------
 # Dispersal      
 # see see (Chapman et al., 2007; Holloway et al., 2016)

src/main/java/analysis/AnalyseSpawnerFeatures.java

@@ -20,7 +20,6 @@
 package analysis;
 
 import java.util.ArrayList;
-import java.util.Arrays;
 import java.util.HashMap;
 import java.util.List;
 import java.util.Map;
@@ -163,7 +162,7 @@ public class AnalyseSpawnerFeatures extends AquaNismsGroupProcess<DiadromousFish
 
 					if (nbAgeForPrimiparousFemale > 0.)
 						ageOfPrimiparousFemaleMemories.get(riverBasin)
-						.push(meanAgeForPrimiparousFemale / nbAgeForPrimiparousFemale);
+								.push(meanAgeForPrimiparousFemale / nbAgeForPrimiparousFemale);
 					else
 						ageOfPrimiparousFemaleMemories.get(riverBasin).push(Double.NaN);
 					if (nbAgeForPrimiparousMale > 0.)
@@ -223,13 +222,12 @@ public class AnalyseSpawnerFeatures extends AquaNismsGroupProcess<DiadromousFish
 			range[0] = ageEffective.calculateMedian();
 			range[1] = ageMin;
 			range[2] = ageMax;
-		}
-		else { 
+		} else {
 			range[0] = Double.NaN;
 			range[1] = Double.NaN;
 			range[2] = Double.NaN;
 		}
-		System.out.println(Arrays.toString(range));
+		// System.out.println(Arrays.toString(range));
 		return range;
 	}
 

src/main/java/species/Grow.java

 package species;
 
-
-import java.util.List;
+import org.openide.util.lookup.ServiceProvider;
 
 import com.thoughtworks.xstream.XStream;
 import com.thoughtworks.xstream.io.xml.DomDriver;
 
 import environment.Basin;
-import environment.Time;
-import fr.cemagref.simaqualife.kernel.AquaNismsGroup;
 import fr.cemagref.simaqualife.kernel.processes.AquaNismsGroupProcess;
 import fr.cemagref.simaqualife.kernel.util.TransientParameters.InitTransientParameters;
 import fr.cemagref.simaqualife.pilot.Pilot;
 import miscellaneous.Miscellaneous;
-
-import org.openide.util.lookup.ServiceProvider;
-
 import species.DiadromousFish.Gender;
 import species.DiadromousFish.Stage;
 import umontreal.iro.lecuyer.probdist.NormalDist;
@@ -27,30 +21,32 @@ public class Grow extends AquaNismsGroupProcess<DiadromousFish, DiadromousFishGr
 
 	/**
 	 * temperature minimum for growth
+	 * 
 	 * @unit °C
 	 */
 	private double tempMinGrow = 3.;
 
 	/**
 	 * temperature maximum for growth
+	 * 
 	 * @unit °C
 	 */
-	private double tempMaxGrow = 26.; 
+	private double tempMaxGrow = 26.;
 
 	/**
 	 * temperature optimal for growth
+	 * 
 	 * @unit °C
 	 */
 	private double tempOptGrow = 17.;
 
 	/**
-	 * K, Brody growth rate at optimal temperature
-	 * L = Linf *(1-exp(-K*(t-t0))
+	 * K, Brody growth rate at optimal temperature L = Linf *(1-exp(-K*(t-t0))
+	 * 
 	 * @unit year -1
 	 */
-	private double kOptForFemale= 0.3;
+	private double kOptForFemale = 0.3;
 
-	
 	/**
 	 * @return the kOptForFemale
 	 */
@@ -58,13 +54,16 @@ public class Grow extends AquaNismsGroupProcess<DiadromousFish, DiadromousFishGr
 		return kOptForFemale;
 	}
 
+
 	/**
-	 * @param kOptForFemale the kOptForFemale to set
+	 * @param kOptForFemale
+	 *            the kOptForFemale to set
 	 */
 	public void setkOptForFemale(double kOptForFemale) {
 		this.kOptForFemale = kOptForFemale;
 	}
 
+
 	/**
 	 * @return the kOptForMale
 	 */
@@ -72,95 +71,112 @@ public class Grow extends AquaNismsGroupProcess<DiadromousFish, DiadromousFishGr
 		return kOptForMale;
 	}
 
+
 	/**
-	 * @param kOptForMale the kOptForMale to set
+	 * @param kOptForMale
+	 *            the kOptForMale to set
 	 */
 	public void setkOptForMale(double kOptForMale) {
 		this.kOptForMale = kOptForMale;
 	}
 
 	/**
-	 * K, Brody growth rate at optimal temperature
-	 * L = Linf *(1-exp(-K*(t-t0))
+	 * K, Brody growth rate at optimal temperature L = Linf *(1-exp(-K*(t-t0))
+	 * 
 	 * @unit year -1
 	 */
-	private double kOptForMale= 0.3;
-	
+	private double kOptForMale = 0.3;
+
 	/**
 	 * standart deviation for the lognormal random draw of growth increment
+	 * 
 	 * @unit cm
 	 */
 	private double sigmaDeltaLVonBert = 0.2;
 
 	private transient NormalGen genNormal;
 
-	public static void main(String[] args) { System.out.println((new
-			XStream(new DomDriver())) .toXML(new Grow())); }
+	public static void main(String[] args) {
+		System.out.println((new XStream(new DomDriver())).toXML(new Grow()));
+	}
+
 
 	@Override
 	@InitTransientParameters
 	public void initTransientParameters(Pilot pilot) {
 		super.initTransientParameters(pilot);
-		genNormal = new NormalACRGen( pilot.getRandomStream(),
-				new NormalDist(0., 1.));		
+		genNormal = new NormalACRGen(pilot.getRandomStream(), new NormalDist(0., 1.));
 	}
 
-	
-	
+
 	@Override
 	public void doProcess(DiadromousFishGroup group) {
-		for(Basin basin : group.getEnvironment().getBasins()){
-			if (basin.getFishs(group)!=null) 
-				for(DiadromousFish fish : basin.getFishs(group)){
+		for (Basin basin : group.getEnvironment().getBasins()) {
+			double currentTemperature = basin.getCurrentTemperature(group.getPilot());
+			if (basin.getFishs(group) != null)
+				for (DiadromousFish fish : basin.getFishs(group)) {
 					double muDeltaLVonBert = 0.;
 					double kVonBert = 0.;
 					double growthIncrement = 0.;
-					
-					// 1) calculate the kVonBert 
-					//System.out.println(this.getKOpt(fish, group) );
-						kVonBert = this.getKOpt(fish, group) *
-								Miscellaneous.temperatureEffect(fish.getPosition().getCurrentTemperature(group.getPilot()), tempMinGrow, tempOptGrow, tempMaxGrow);
 
-					// 2) Update the fish length with a lognormal normal draw  of increment
+					// 1) calculate the kVonBert
+					// System.out.println(this.getKOpt(fish, group) );
+					kVonBert = this.getKOpt(fish, group)
+							* Miscellaneous.temperatureEffect(currentTemperature, tempMinGrow, tempOptGrow, tempMaxGrow);
+
+					// 2) Update the fish length with a lognormal normal draw of increment
 					// limit the fish length to Linf
-					if (fish.getLength() < group.getLinfVonBert(fish)){
-						muDeltaLVonBert = Math.log((group.getLinfVonBert(fish) - fish.getLength()) * (1 - Math.exp(-kVonBert * group.getEnvironment().getTime().getSeasonDuration()))) - (sigmaDeltaLVonBert*sigmaDeltaLVonBert)/2;
-						growthIncrement = Math.exp(genNormal.nextDouble()*sigmaDeltaLVonBert + muDeltaLVonBert);
-						
-						fish.setLength(Math.min(group.getLinfVonBert(fish), fish.getLength() + growthIncrement));											
-					}
-					else {
+					double lengthBefore = fish.getLength();
+					double alea = Double.NaN;
+					if (fish.getLength() < group.getLinfVonBert(fish)) {
+						muDeltaLVonBert = Math
+								.log((group.getLinfVonBert(fish) - fish.getLength())
+										* (1 - Math.exp(-kVonBert * group.getEnvironment().getTime().getSeasonDuration())))
+								- (sigmaDeltaLVonBert * sigmaDeltaLVonBert) / 2;
+						alea = genNormal.nextDouble();
+						growthIncrement = Math.exp(alea * sigmaDeltaLVonBert + muDeltaLVonBert);
+
+						fish.setLength(Math.min(group.getLinfVonBert(fish), fish.getLength() + growthIncrement));
+					} else {
 						fish.setLength(group.getLinfVonBert(fish));
 					}
-					//System.out.println(fish.getAge() + " -> "+ fish.getLength() + " ("+fish.getStage()+"): "+ growthIncrement);
+
+					// System.out.println(basin.getName() + " :" + fish.getAge() + " " + fish.getGender() + " : " +
+					// lengthBefore
+					// + " -> " + fish.getLength() + " (" + fish.getStage() + "): " + growthIncrement + " with kVonBert
+					// = "
+					// + kVonBert + " (" + currentTemperature + " / " + alea + ")");
+
 					// test if fish become mature
-					if (fish.getStage() == Stage.IMMATURE && fish.getLength() > group.getlFirstMaturity(fish)){
-							fish.setStage(Stage.MATURE); 
+					if (fish.getStage() == Stage.IMMATURE && fish.getLength() > group.getlFirstMaturity(fish)) {
+						fish.setStage(Stage.MATURE);
 					}
-					//System.out.println("la temp�rature du lieu de vie du poisson est :" + fish.getPosition().getCurrentTemperature() + ", la saison est :" + Time.getSeason() + " et sa nouvelle taille est :" + fish.getLength());
+					// System.out.println("la temp�rature du lieu de vie du poisson est :" +
+					// fish.getPosition().getCurrentTemperature() + ", la saison est :" + Time.getSeason() + " et sa
+					// nouvelle taille est :" + fish.getLength());
 				}
 		}
 	}
-	
+
+
 	/**
 	 * @param fish
 	 * @param group
-	 * @return the Brody coeff   from Diadromousgroup if exists or from this grow process
-	 * depends of the fish gender .In case of undifferentiaced fish, the mean for male and female is considered
+	 * @return the Brody coeff from Diadromousgroup if exists or from this grow process depends of the fish gender .In
+	 *         case of undifferentiaced fish, the mean for male and female is considered
 	 */
-	public  double getKOpt(DiadromousFish fish, DiadromousFishGroup group) {
-		 double kOpt = group.getKOpt(fish);
-		 if (Double.isNaN(kOpt)){ // no definition for the group
-				if (fish.getGender() == Gender.FEMALE)
-					kOpt = kOptForFemale;
-				else if (fish.getGender() == Gender.MALE)
-					kOpt = kOptForMale;
-				else
-					kOpt=  (kOptForFemale + kOptForMale) / 2.;
-		 }	 
-		 
-		 return kOpt;
-		
-		
+	public double getKOpt(DiadromousFish fish, DiadromousFishGroup group) {
+		double kOpt = group.getKOpt(fish);
+		if (Double.isNaN(kOpt)) { // no definition for the group
+			if (fish.getGender() == Gender.FEMALE)
+				kOpt = kOptForFemale;
+			else if (fish.getGender() == Gender.MALE)
+				kOpt = kOptForMale;
+			else
+				kOpt = (kOptForFemale + kOptForMale) / 2.;
+		}
+
+		return kOpt;
+
 	}
 }

src/main/java/species/SurviveWhenBecomeOld.java

+package species;
+
+import java.util.ListIterator;
+
+import org.openide.util.lookup.ServiceProvider;
+
+import com.thoughtworks.xstream.XStream;
+import com.thoughtworks.xstream.io.xml.DomDriver;
+
+import environment.Basin;
+import environment.Time.Season;
+import fr.cemagref.simaqualife.kernel.processes.AquaNismsGroupProcess;
+import fr.cemagref.simaqualife.pilot.Pilot;
+
+@ServiceProvider(service = AquaNismsGroupProcess.class)
+public class SurviveWhenBecomeOld extends AquaNismsGroupProcess<DiadromousFish, DiadromousFishGroup> {
+
+	/**
+	 * fish older than or equal to ageMax are killed
+	 * 
+	 * @unit: year
+	 */
+	private double ageMax = 12;
+	private boolean verbose = false;
+
+	/**
+	 * the random numbers generator for binomial draws
+	 * 
+	 * @unit --
+	 */
+	// private transient BinomialForSuperIndividualGen aleaGen;
+
+	public static void main(String[] args) {
+		System.out.println((new XStream(new DomDriver())).toXML(new SurviveWhenBecomeOld()));
+	}
+
+
+	@Override
+	public void initTransientParameters(Pilot pilot) {
+		super.initTransientParameters(pilot);
+		// aleaGen = new BinomialForSuperIndividualGen(pilot.getRandomStream());
+	}
+
+
+	@Override
+	public void doProcess(DiadromousFishGroup group) {
+		if (group.getEnvironment().getTime().getSeason(group.getPilot()) == Season.SPRING) {
+			Long nbKill = 0L;
+			for (Basin basin : group.getEnvironment().getBasins()) {
+				// System.out.print(basin.getName());
+
+				if (basin.getFishs(group) != null) {
+					// System.out.println(" y a des poissons");
+					ListIterator<DiadromousFish> fishIterator = basin.getFishs(group).listIterator();
+					while (fishIterator.hasNext()) {
+						DiadromousFish fish = fishIterator.next();
+
+						if (fish.getAge() >= ageMax) {
+							nbKill += fish.getAmount();
+							fishIterator.remove();
+						}
+					} // end on loop of fish
+				}
+			} // end loop on basin
+
+			if (verbose & nbKill > 0L)
+				System.out.println("  number of old fish killed: " + nbKill);
+		}
+	}// end of doprocess
+}// end of class