/* * #%L * IsisFish data * %% * Copyright (C) 2006 - 2011 Ifremer, CodeLutin * %% * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as * published by the Free Software Foundation, either version 2 of the * License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public * License along with this program. If not, see * . * #L% */ package rules; import fr.ifremer.isisfish.datastore.ResultStorage; import fr.ifremer.isisfish.entities.Metier; import fr.ifremer.isisfish.entities.Population; import fr.ifremer.isisfish.entities.PopulationGroup; import fr.ifremer.isisfish.rule.AbstractRule; import fr.ifremer.isisfish.simulator.SimulationContext; import fr.ifremer.isisfish.types.Month; import fr.ifremer.isisfish.types.TimeStep; //import fr.ifremer.isisfish.util.Doc; import fr.ifremer.isisfish.annotations.Doc; import java.util.Arrays; import java.util.HashMap; import java.util.Map; import java.util.List; import java.util.ArrayList; import java.util.Collection; import org.apache.commons.logging.Log; import org.apache.commons.logging.LogFactory; import org.nuiton.math.matrix.*; import org.nuiton.math.matrix.MatrixFactory; import org.nuiton.math.matrix.MatrixIterator; import org.nuiton.math.matrix.MatrixND; import scripts.MinimisationUtil; import scripts.ObjectiveFunctionBaranov; import scripts.RuleUtil; import scripts.SiMatrix; import resultinfos.*; /** * Le but de cette regle est de permettre de simuler une harvest control rule, * et plus precisemment la transition entre une gestion par approche de precaution * et une gestion par le MSY. * * La premiere annee de transition est 2011 et la derniere 2014, pour une gestion * au MSY des 2015. * *Ici cette regle va s'appliquer a partir de la quatrieme annee de simu (annee initiale=2008) * et a chacune des 4 populations simulees, alors que jusqu'a 2010 on avait un TAC commun pour la plie. * * Created: * * @author Sigrid * @version $Revision: 22 aout 2018 * Nouveau calcul du TAC * * Last update: $Date: 131212 $ * by : $Author: Loic $ */ public class HCR_PGY_Fupper_multisp extends AbstractRule { /** to use log facility, just put in your code: log.info("..."); */ static private Log log = LogFactory.getLog(HCR_PGY_Fupper_multisp.class); public MatrixND SSBinitYearPop ; //public MatrixND MatTest ; //@Doc("Affected population") //public Population param_pop = null; @Doc("Populations considered considered in the Pretty Good Yield computation") public String param_popNames = "Lepidorhombus_whiffiagonis;Leucoraja_naevus;Lophius_piscatorius;Merluccius_merluccius;Nephrops_norvegicus;Raja_clavata;Solea_solea"; @Doc("MSY Btrigger IN TONS of populations specified in param_popNames") public String param_msyBtrigger = "9882.001;1884.756;8276.524;291828.298;0.01;1429.915;10437.234"; // Attention aux unites : kg VS t @Doc("FMSY for populations specified in param_popNames") public String param_fmsy = "0.265;0.303;0.102;0.197;0.326;0.095;0.252"; // Default values are estimates from Ricouard, pers. comm. 28/05/26 @Doc("Fupper for populations specified in param_popNames") public String param_fupper = "0.361;0.384;0.157;0.281;0.389;0.117;0.420"; // //@Doc("Fpa for populations specified in param_popNames") //public String param_fpa = "0.4;0.4"; // @Doc("Begin date") public TimeStep param_beginStep = new TimeStep(12); @Doc("End date") public TimeStep param_endStep = new TimeStep(480); //@Doc("Fpa") //public double param_fpa = 0.4; //@Doc("Fmsy") //public double param_fmsy = 0.29; //@Doc("MSY Btrigger") //public double param_msyBtrigger = 8000; // Attention aux unites : kg VS t //@Doc("Duration of transition period (year)") //public double param_transitionDuration = 5; @Doc("Variation annuelle max du TAC ") public double param_varTac = 0.20; //@Doc("TAC the year before it starts in tons") //public double param_Tac_startYear = 4219; protected double Fstart; boolean affectation = false; public double param_RelStability = 0.57; public boolean param_LandingObligation = true; @Doc("Discard rate to use if TAC is based on catches (and not landings)") public double param_DiscardsRate = 10; protected MatrixND TacPopMatrix; protected Map abaaMap; protected Map waaMap; protected Map avRMap; protected Map ssbMap; List listPopNames; List listPop = new ArrayList(); protected String[] necessaryResult = { // put here all necessary result for this rule // example: MatrixTacPerPop.NAME, MatrixBiomass.NAME, MatrixTotalFishingMortality.NAME, MatrixFishingMortalityPerGroup.NAME, MatrixAbundance.NAME }; /** * @return the necessaryResult */ @Override public String[] getNecessaryResult() { return this.necessaryResult; } /** * Permet d'afficher a l'utilisateur une aide sur la regle. * * @return L'aide ou la description de la regle */ @Override public String getDescription() { return "Harvest Control Rule simulating the transition between a precautionary approach management and a MSY management, for a given populaion"; } /** * Appele au demarrage de la simulation, cette methode permet d'initialiser * des valeurs * * @param context La simulation pour lequel on utilise cette regle */ @Override public void init(SimulationContext context) throws Exception { SiMatrix siMatrix = SiMatrix.getSiMatrix(context); TimeStep date = new TimeStep(0); String[] allPopNames = param_popNames.split(";"); List allPop = context.getPopulationDAO().findAll(); List allPopNamesList = Arrays.asList(allPopNames); String[] StrBtrigger = param_msyBtrigger.split(";"); String[] StrFmsy = param_fmsy.split(";"); String[] StrFupper = param_fupper.split(";"); //String[] StrFpa = param_fpa.split(";"); for (Population pop : allPop){ //if(pop.getName().matches(allPopNames)){ if(allPopNamesList.contains(pop.getName())){ listPop.add((Population) pop); } System.out.println("taille de list pop :" + listPop.size()); } for ( int i = 0 ; i < allPopNames.length; i++) { context.setValue("Btrigger_"+allPopNames[i], Double.parseDouble(StrBtrigger[i])); context.setValue("Fmsy_"+allPopNames[i], Double.parseDouble(StrFmsy[i])); context.setValue("Fupper_"+allPopNames[i], Double.parseDouble(StrFupper[i])); //context.setValue("Fpa_"+allPopNames[i], Double.parseDouble(StrFpa[i])); } // Add TAC values to a vector indexed by TAC year (from TAC starting year on) // String[] StrTacSeries = param_tacInKgs.split(";"); // tacSeries = new double[StrTacSeries.length]; // for ( int i = 0 ; i < StrTacSeries.length; i++) { // //for ( int i = 0 ; i < tacSeries.length; i++) { // if(param_tacInKgs == "0"){ // tacSeries[i] = 0 ; // }else { // tacSeries[i] = Double.parseDouble(StrTacSeries[i]); // } // } } /** * La condition qui doit etre vrai pour faire les actions. * * @param context la simulation pour lequel on utilise cette regle * @param step le pas de temps courant * @param metier le metier concerne * @return vrai si on souhaite que les actions soit faites */ @Override public boolean condition(SimulationContext context, TimeStep step, Metier metier) throws Exception { boolean result = false; if (step.afterOrEquals(param_beginStep) && step.before(param_endStep) && step.getMonth() == Month.JANUARY){ result = true; affectation =false; } return result; } /** * Si la condition est vrai alors cette action est executee avant le pas * de temps de la simulation. * * @param context la simulation pour lequel on utilise cette regle * @param step le pas de temps courant * @param metier le metier concerne */ @Override public void preAction(SimulationContext context, TimeStep step, Metier metier) throws Exception { if( !affectation){ //log.info("Mise en place de l'HCR"); ResultStorage matResult = context.getSimulationStorage().getResultStorage(); // MatrixND MatTest = matResult.getMatrix(step.previous(), MatrixTotalFishingMortality.NAME); // MatrixND MatTest = matResult.getMatrix(step.previous(), MatrixTotalFishingMortality.NAME); //MatTest_PGY = getSubMatrix(0, listPop.toArray(new TimeStep[0])); //getSubMatrix(0, sequenceOrhago.toArray(new TimeStep[0])); // creer un MAp de stockage de Ab@age et rec@age et W@age indicé par pop boolean Fupper_switch = true; // loop for Btrigger abaaMap = new HashMap (); waaMap = new HashMap (); // util ? maaMap = new HashMap (); avRMap = new HashMap (); ssbMap = new HashMap (); for (Population pop : listPop) { // To do : if pop.getName().equals(merlu!) // prendre direct la ssb de previous month B * mat // Import de la biomasse de l'annee en cours //Calcul SSB, Biomass, Abondance and mean weight MatrixND MatAb = matResult.getMatrix(step.previous(), pop, MatrixAbundance.NAME); MatAb = MatAb.sumOverDim(1).reduce(); abaaMap.put(pop,MatAb); List listWaa = new ArrayList(); List listMaa = new ArrayList(); // Calcul de la SSB de l'annee en cours (anticipe janvier a partir des nombres de decembre + average R...) double SSBiom = 0; List groups = pop.getPopulationGroup(); /*int size = pop.getPopulationGroup().size()+1; double[] mwa = new double[size]; double[] mat = new double[size];*/ for(PopulationGroup pg:groups){ //mwa[pg.getId()] = pg.getMeanWeight(); //mat[pg.getId()] = pg.getMaturityOgive(); listWaa.add(pg.getMeanWeight()); listMaa.add(pg.getMaturityOgive()); } listMaa.add(listMaa.get(listMaa.size() -1)); listWaa.add(listWaa.get(listWaa.size() -1)); waaMap.put(pop,listWaa); //maaMap.put(pop,listMaa); System.out.println("Pop "+pop.getName() + "waa : "+ listWaa); //mwa[size-1] = mwa[size-2]; //mat[size-1] = mat[size-2]; /* to do : average of last 3 years // double avR = context.get(pop).getAsDouble("averageR"); // entré dans l'interface de la pop */ //TimeStep tmpStep = new TimeStep(step.getStep()); int past3years = Math.min(step.getStep(),36); System.out.println("past3year " + past3years); //TimeStep ThreeyearsAgo = new TimeStep(step.getStep() - past3years); double avR = 0; int threeyearsago = step.getStep() - past3years; for (TimeStep tmpStep = new TimeStep(threeyearsago); !tmpStep.after(step.previous()); tmpStep = tmpStep.next()) { double tmpRec = matResult.getMatrix(tmpStep, pop ,MatrixRecruitment.NAME).sumAll(); // context.getPopulationMonitor().getRecruitment(tmpStep, pop).sumAll(); avR += tmpRec; System.out.println("dans loop rec, step "+tmpStep.getStep() + " tmpRec = "+tmpRec); } avR = avR/past3years*12; avRMap.put(pop,avR); System.out.println("Pop "+pop.getName() + "avR : "+avR); for (MatrixIterator i=MatAb.iterator(); i.hasNext();) { i.next(); Object [] sems = i.getSemanticsCoordinates(); PopulationGroup group = (PopulationGroup)sems[0]; double val = i.getValue()* listMaa.get(group.getId()+1) * listWaa.get(group.getId()+1); SSBiom = SSBiom + val; } SSBiom = SSBiom + avR * listMaa.get(0) * listWaa.get(0); SSBiom = SSBiom / 1000; ssbMap.put(pop,SSBiom); System.out.println("Pop "+pop.getName() + "ssb : "+SSBiom); // if(context.getValue("errSSB") != null) { // double[] errlist = (double[])context.getValue("errSSB"); // SSBiom = SSBiom * errlist[step.getYear()-3]; //} //---------------------------------------------------------------------------------------------------- /////// HCR /////////////////// // Comparaison de la SSB a MSYBtrigger et calcul de fMsyHcr double fMsyHcr = 0; // if the pop is part of a larger stock on which the HCR is based //SSBiom = SSBiom / context.get(pop).getAsDouble("PropOfTheStock") ; if ( SSBiom < (double)context.getValue("Btrigger_"+pop.getName())) { // => ftarget groups = pop.getPopulationGroup(); int size = pop.getPopulationGroup().size()+1; double[] mwa = new double[size]; double[] mat = new double[size]; for(PopulationGroup pg:groups){ mwa[pg.getId()] = pg.getMeanWeight(); mat[pg.getId()] = pg.getMaturityOgive(); } mwa[size-1] = mwa[size-2]; mat[size-1] = mat[size-2]; double avR = context.get(pop).getAsDouble("averageR"); */ // Calcul de la SSB de l'annee en cours (anticipe janvier a partir des nombres de decembre + average R...) //double SSBiom = 0; MatrixND MatAb = abaaMap.get(pop); List listWaa = waaMap.get(pop); double avR = avRMap.get(pop); // Calcul Fbar de l'anne double Fbar; // Fbar ISIS MatrixND MatFy = matResult.getMatrix(step.previous(), pop, MatrixTotalFishingMortality.NAME); // Fbar y-1 //............................................................................. // Foth bar if necessary /*int grpNbbar = (int)(pop.getFbarGroupMax()-pop.getFbarGroupMin() +1); MatrixND fothaa = pop.getFishingMortalityOtherFleetsMatrix().meanOverDim(1).reduce(); // mean over zones double FbarOTH = fothaa.getSubMatrix(0, pop.getFbarGroupMin(),grpNbbar).meanAll(); */ //........................................................................... Fbar = MatFy.getValue(0) ;//+ FbarOTH; System.out.println("Pop "+pop.getName() + "Fbar : "+Fbar); /* // Store Fstart in context if(step.getYear() == param_beginStep.getYear()){ context.setValue("Fstart_"+pop.getName(),Fbar); Fstart = Fbar; // hyp F de l'an dernier = F de cette annee } else { Fstart = (Double) context.getValue("Fstart_"+pop.getName()); }*/ double fMsyHcr; if(Fupper_switch == true){ fMsyHcr = (double)context.getValue("Fupper_" + pop.getName()); }else{ fMsyHcr = (double)context.getValue("Fref_" + pop.getName()); } /* Pas de transition // Transition (si besoin) : Calcul des multiplicateurs a appliquer a Fstart et fMsyHCR double multFstart = 1 - ((step.getYear() - (param_beginStep.getYear() - 1)) * (1/param_transitionDuration)); if (multFstart < 0){ multFstart = 0; } double multFmsy = (step.getYear() - (param_beginStep.getYear() - 1)) * (1/param_transitionDuration); if (multFmsy > 1){ multFmsy = 1; } // Calcul de fMsyHcrTransition double fMsyHcrTransition = multFstart * Fstart + multFmsy * fMsyHcr; // Si fMsyHcrTransition > fpa on applique fpa double AppliedF = 0; if (fMsyHcrTransition < (double)context.getValue("Fpa_"+pop.getName())) { AppliedF = fMsyHcrTransition; } else { AppliedF = (double)context.getValue("Fpa_"+pop.getName()); } */ double AppliedF = fMsyHcr; System.out.println("Pop "+pop.getName() + "appliedF : "+AppliedF); //--------------------------------------------------------------------------------------------------- // Calcul du niveau de captures correspondant au F choisi (hypothese Ab dec precedent connue) double TAC = 0; double multF = AppliedF/ Fbar; MatrixND MatFyPG = matResult.getMatrix(step.previous(), pop, MatrixFishingMortalityPerGroup.NAME); // + Foth for(PopulationGroup pg:pop.getPopulationGroup()){ double FPGtar = (MatFyPG.getValue(pg))*multF; double MPG = pop.getNaturalDeathRateMatrix().meanOverDim(1).reduceDims(1).getValue(pg.getId()); double NPG; //double WPG = (mwa[pg.getId()] + mwa[pg.getId()+1])/2; // otherwise weight in january double WPG = (listWaa.get(pg.getId()) + listWaa.get(pg.getId()+1))/2; // otherwise weight in january if(pg.getId() == 0){ //= recruitment, new N assumed = average R NPG = avR; } else if (pg.getId() == (pop.getPopulationGroup().size()-1) && pop.getPlusGroup() ){ // last group NPG = MatAb.getValue(pg.getId()-1) + MatAb.getValue(pg.getId()); // +group } else NPG = MatAb.getValue(pg.getId()-1); // if group change in january, ab(year-1,pg-1) = ab(year,pg) (MatAb previous step = december) TAC += (FPGtar/(FPGtar+MPG))*(1-Math.exp(-(FPGtar+MPG)))*NPG*WPG ; } System.out.println("Pop "+pop.getName() + "multF : "+multF); System.out.println("Pop "+pop.getName() + "TAC : "+TAC); //---------------------------------------------------------------------------------------------------- // Total allowable landing and TAC uplifts... /*if(context.get(pop).getAsDouble("FonCatch") == 1 && !param_LandingObligation){ // If TAC is based on a F computed on catches, remove the part of discards TAC = (1-param_DiscardsRate)*TAC; }else if(context.get(pop).getAsDouble("FonCatch") == 0 && param_LandingObligation){ // If TAC is based on a F computed on landing, add a TAC uplift TAC = (1+param_DiscardsRate)*TAC; }*/ //---------------------------------------------------------------------------------------------------- // Si niveau de captures inferieur de plus de 15 pcts au TAC de l'annee precedente limiter a 15 pcts double TacReel = TAC; /*double diffTac = 0; double TacReel = 0; double TacPrevYr = 0; if (step.getYear() == param_beginStep.getYear()){ // on s'embete pas abev la premiere annee /*diffTac = TAC / (param_Tac_startYear*1000) - 1; //log.info("diffTac = " + diffTac); if (diffTac < -1 * param_varTac){ TacReel = (param_Tac_startYear*1000) - param_varTac * (param_Tac_startYear*1000); } else if (diffTac >= param_varTac){ TacReel = (param_Tac_startYear*1000) + param_varTac * (param_Tac_startYear*1000); } else { TacReel = TAC; }*/ /*}else { // retrouver le niveau de TAC de l'annee precedente TacPrevYr = (Double) context.getValue("TAC_" + pop.getName()); //log.info("TacPrevYr = " + TacPrevYr); diffTac = TAC / TacPrevYr - 1; //log.info("diffTac = " + diffTac); if (diffTac < -1 * param_varTac){ TacReel = TacPrevYr - param_varTac * TacPrevYr; } else if (diffTac >= param_varTac){ TacReel = TacPrevYr + param_varTac * TacPrevYr; } else { TacReel = TAC; } } TacReel = TacReel;*/ //---------------------------------------------------------------------------------------------------- // Faire en sorte que le TAC soit utilise par la regle TACPoidsPop_PourHCR context.setValue("TAC_PGY_" + pop.getName() , TacReel); //TacReel context.setValue("TACfr_PGY_"+pop.getName(), TacReel * param_RelStability); //---------------------------------------------------------------------------------------------------- // FOTH // Modify natural mortality /*double newmult; if(context.getValue("multF_"+ pop.getName()) !=null){ double oldmult = (Double) context.getValue("multF_"+pop.getName()); newmult = oldmult * multF; } else newmult = multF; context.setValue("multF_"+ pop.getName(),newmult);*/ //---------------------------------------------------------------------------------------------------- // Store infos in MatrixTacPerPop int fup = 0; if(Fupper_switch)fup=1; TacPopMatrix = MatrixFactory.getInstance().create( MatrixTacPerPop.NAME, new List[] { Arrays.asList(new String[]{"TAC","TACfr","SSBiom","AppliedF","TACavtVar","Fbar","FbarOTH","multF","StopFishing","Fupper"})}); TacPopMatrix.setValue("TAC", TacReel); TacPopMatrix.setValue("TACfr", TacReel * param_RelStability); TacPopMatrix.setValue("SSBiom", ssbMap.get(pop)); TacPopMatrix.setValue("AppliedF", AppliedF); //TacPopMatrix.setValue("TACavtVar", TAC); TacPopMatrix.setValue("Fbar", Fbar); ///TacPopMatrix.setValue("FbarOTH", FbarOTH); TacPopMatrix.setValue("multF", multF); TacPopMatrix.setValue("Fupper", fup); matResult.addResult(step, pop, TacPopMatrix); } affectation = true; } } /** // * Si la condition est vrai alors cette action est executee apres le pas * de temps de la simulation. * * @param context La simulation pour lequel on utilise cette regle * @param step le pas de temps courant * @param metier le metier concerne */ @Override public void postAction(SimulationContext context, TimeStep step, Metier metier) throws Exception { // Creation de la matrice de stockage de TAC } }