/* * #%L * IsisFish data * %% * Copyright (C) 2007 - 2015 Ifremer, CodeLutin, Sigrid Lehuta * %% * 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 3 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 simulationplans; import java.io.File; import java.util.ArrayList; import java.util.Arrays; import java.util.Collections; import java.util.List; import org.apache.commons.logging.Log; import org.apache.commons.logging.LogFactory; import org.nuiton.math.matrix.MatrixIterator; import org.nuiton.math.matrix.MatrixND; import org.nuiton.topia.TopiaContext; import org.nuiton.util.StringUtil; import fr.ifremer.isisfish.datastore.ResultStorage; import fr.ifremer.isisfish.datastore.SimulationStorage; import fr.ifremer.isisfish.entities.Population; import fr.ifremer.isisfish.entities.PopulationGroup; import fr.ifremer.isisfish.simulator.SimulationPlan; import fr.ifremer.isisfish.simulator.SimulationPlanContext; import fr.ifremer.isisfish.annotations.Doc; import fr.ifremer.isisfish.util.ScriptUtil; import resultinfos.MatrixCatchWeightPerStrategyMetPerZonePop; /** * CalibrationSimplexePasVariable2Capturabilite. * * Created: 17 septembre 2007 * * @author sigrid lehuta * @version $Revision: 1.27 $ * * Last update: $Date: 2007/05/24 09:29:18 $ * by : $Author: bpoussin $ */ // //////////////////////////////////////////////////////////////////////// // USER GUIDE //------------------------------------------------------------------------- // Script must be adapted to the case study (calibration data, catchability assumptions) // Script won't compile as it stands // Comments preceeded by /////*** explain where and how to adapt the script // Access to the APIs is free // //////////////////////////////////////////////////////////////////////// // ***You can modify class name if you want // ***BUT attention : file name and class name must be the same (without the extention ".java"), // ie here : "CalibrationSimplexePasVariable2Capturabilite" public class CalibrationSimplexePasVariable2Capturabilite implements SimulationPlan { /** to use log facility, just put in your code: log.info("..."); */ static private Log log = LogFactory .getLog(CalibrationSimplexePasVariable2Capturabilite.class); enum State { STATE_INIT, STATE_0, STATE_1, STATE_2, STATE_3, STATE_4 } /////***here must appear the path to export the historic file ("Historic.csv") //in which q1, q2 and objective function computed at each simulation are stored /////***Attention : before beginning a new calibration rename any potential //old Historic.csv files or they will be lost protected File exportHistoric = new File("Historic.csv"); protected String exportHisto = ""; @Doc("Population which parameters are calibrated") public Population param_Population = null; @Doc("First initial point of the simplex: de la forme(\"xx;yy\")") public String param_M1 = "2.42e-5;2.11e-6";// devient un parametre du plan d analyse @Doc("Second initial point of the simplex") public String param_M2 = "2.34e-5;2.59e-6";// devient un parametre du plan d analyse @Doc("Third initial point of the simplex") public String param_M3 = "2.59e-5;2.41e-6";// devient un parametre du plan d analyse // ***put here the path and name of the file containing the data used to calibrate // your fishery ( here observed landings per season and age groups) @Doc(value = "file name and path of observed landings") public String param_nomfichier_debarquements = "Observedlandings2001-2003.csv"; protected File debarquementsObserves; protected MatrixND matrixDebarquement; protected State state = State.STATE_INIT; public Experiences experiences = new Experiences(); //*** write the name of the simulated matrix that contains the data corresponding // to your observations (here MATRIX_CATCH_WEIGHT_PER_STRATEGY_MET_PER_ZONE_POP) protected String[] necessaryResult = { MatrixCatchWeightPerStrategyMetPerZonePop.NAME }; public String[] getNecessaryResult() { return this.necessaryResult; } /** * Permet d'afficher a l'utilisateur une aide sur le plan. * @return L'aide ou la description du plan */ public String getDescription() throws Exception { return "Calibration using variable step Simplex method (Walters): user" + "gives a file of observations (here catches)(.csv), simulated output" + "will try to approach oservations by changing the values of catchability"; } /** * Appele au demarrage de la simulation, cette methode permet d'initialiser * des valeurs * @param context La simulation pour lequel on utilise cette regle */ public void init(SimulationPlanContext context) throws Exception { if (param_nomfichier_debarquements == null || "".equals(param_nomfichier_debarquements)) { debarquementsObserves = ScriptUtil.getFile(".*.csv", "fichier csv séparateur ';'"); } else { debarquementsObserves = new File(param_nomfichier_debarquements); } // ***Create the matrix named matrixDebarquement that will contain your observed landings // method: MatrixFactory.getInstance().create() // ***Then import your file in it // method : matrixDebarquement.importCSV(); log.info("MatrixDebarquement : " + matrixDebarquement); } double g1; double g2; double worst1; double worst2; /** * Call before each simulation. * * @param context plan context * @param nextSimulation storage used for next simulation * @return true if we must do next simulation, false to stop plan * @throws Exception */ public boolean beforeSimulation(SimulationPlanContext context, SimulationStorage nextSimulation) throws Exception { // the whole beforeSimulation part codes for the variable step simplex // algorithm (Walters et al.1999), no change to bring boolean doNext = true; boolean doBoucle = true; log.info("before simulation"); int number = nextSimulation.getParameter().getSimulationPlanNumber(); if (number < 3) { log.info("number<3"); String[] M1 = param_M1.split(";"); String[] M2 = param_M2.split(";"); String[] M3 = param_M3.split(";"); double[] q1 = StringUtil.toArrayDouble(M1[0], M2[0], M3[0]); double[] q2 = StringUtil.toArrayDouble(M1[1], M2[1], M3[1]); experiences.getExperience(number).q1 = q1[number]; experiences.getExperience(number).q2 = q2[number]; changeDB(experiences.getExperience(number), nextSimulation); } else { double q1 = 1000; double q2 = 1000; double lastCritere = experiences.getExperience(number - 1).criteria; while (doBoucle) { doBoucle = false; if (state == State.STATE_INIT) { doBoucle = false; log.info("state init"); //ordonne les 3 premieres experiences selon leur critere Collections.sort(experiences.current); //log.info("SIMPLEXE : current 0 = " + experiences.current.get(0).criteria + // "current 1 = " + experiences.current.get(1).criteria + "current 2 = " + // experiences.current.get(2).criteria); log.info("SIMPLEXE : current 0 = " + experiences.current.get(0).criteria + "current 1 = " + experiences.current.get(1).criteria + "current 2 = " + experiences.current.get(2).criteria); log.info("SIMPLEXE : Best q1 = " + experiences.current.get(0).q1 + " q2 = " + experiences.current.get(0).q2); log.info("SIMPLEXE : NextBest q1 = " + experiences.current.get(1).q1 + " q2 = " + experiences.current.get(1).q2); log.info("SIMPLEXE : Worst q1 = " + experiences.current.get(2).q1 + " q2 = " + experiences.current.get(2).q2); //Calcul et evaluation de R double g1 = (experiences.current.get(0).q1 + experiences.current .get(1).q1) / 2.0; double g2 = (experiences.current.get(0).q2 + experiences.current .get(1).q2) / 2.0; double worst1 = experiences.current.get(2).q1; double worst2 = experiences.current.get(2).q2; state = State.STATE_0; q1 = 2 * g1 - worst1; q2 = 2 * g2 - worst2; log.info("R : q1 = " + q1 + " q2 = " + q2); } else if (state == State.STATE_0) { doBoucle = false; log.info("state 0"); // on fait la 5eme avec des q qui dependent de la 4eme dans le dernier cas //log.info("g1 = " + g1 + " " + "g2 = " + g2); //log.info("worst1 = " + worst1 + " " + "worst2 = " + worst2); if (lastCritere > experiences.current.get(2).criteria) { log.info("State 0 : R : lastCtritere > current2 : R pire de W"); state = State.STATE_1; //calcul de Cw q1 = ((experiences.current.get(0).q1 + experiences.current .get(1).q1) / 2.0) - (((experiences.current.get(0).q1 + experiences.current .get(1).q1) / 2.0) - experiences.current .get(2).q1) / 2.0; q2 = ((experiences.current.get(0).q2 + experiences.current .get(1).q2) / 2.0) - (((experiences.current.get(0).q2 + experiences.current .get(1).q2) / 2.0) - experiences.current .get(2).q2) / 2.0; log.info("Cw : q1 = " + q1 + " q2 = " + q2); } else if (lastCritere > experiences.current.get(1).criteria) { log.info("State 0 :R : lastCritere > current 1 : R meilleur que W et moins bon que N"); state = State.STATE_2; // calcul de Cr q1 = ((experiences.current.get(0).q1 + experiences.current .get(1).q1) / 2.0) + (((experiences.current.get(0).q1 + experiences.current .get(1).q1) / 2.0) - experiences.current .get(2).q1) / 2.0; q2 = ((experiences.current.get(0).q2 + experiences.current .get(1).q2) / 2.0) + (((experiences.current.get(0).q2 + experiences.current .get(1).q2) / 2.0) - experiences.current .get(2).q2) / 2.0; log.info("Cr : q1 = " + q1 + " q2 = " + q2); } else if (lastCritere > experiences.current.get(0).criteria) { log.info("State 0 :R : lastCritere > current0 : R meilleur que N et moins bon que B"); state = State.STATE_INIT; experiences.current.remove(2);//remove(3)avant doBoucle = true; log.info("remove W, simplex BNR"); } else { // dernier cas possible: if (lastCritere < experiences.current.get(0).critere) { log.info("State 0 :R : lastCritere < current 0 : R meilleur que B, calcul de E"); state = State.STATE_4; q1 = experiences.getExperience(number - 1).q1 + (experiences.current.get(0).q1 + experiences.current .get(1).q1) / 2.0 - experiences.current.get(2).q1; q2 = experiences.getExperience(number - 1).q2 + (experiences.current.get(0).q2 + experiences.current .get(1).q2) / 2.0 - experiences.current.get(2).q2; //q1 = experiences.current.get(3).q1 + (experiences.current.get(0).q1 + experiences.current.get(1).q1) / 2.0 - experiences.current.get(2).q1; //q2 = experiences.current.get(3).q2 + (experiences.current.get(0).q2 + experiences.current.get(1).q2) / 2.0 - experiences.current.get(2).q2; log.info("E : q1 = " + q1 + " q2 = " + q2); } } else if (state == State.STATE_1) { log.info("state 1, simplex BNCw"); experiences.current.remove(3); experiences.current.remove(2); state = State.STATE_INIT; doBoucle = true; } else if (state == State.STATE_2) { log.info("state 2, simplex BNCr"); experiences.current.remove(3); experiences.current.remove(2); state = State.STATE_INIT; doBoucle = true; } else if (state == State.STATE_4) { log.info("state 4 :comparaison de E a B"); doBoucle = true; if (lastCritere < experiences.current.get(0).criteria) { log.info("E meilleur que B, remove 2 et 3 : simplex BNE"); experiences.current.remove(3); experiences.current.remove(2); } else { log.info("E moins bon que B, remove 2 et 4, simplex BNR"); experiences.current.remove(4); experiences.current.remove(2); } state = State.STATE_INIT; } }//fin du while //on remplit la table experiences experiences.getExperience(number).q1 = q1; experiences.getExperience(number).q2 = q2; log.info("on change Q dans la DB avec : q1 = " + q1 + " " + "q2 = " + q2); // on change la valeur de q dans la DB changeDB(experiences.getExperience(number), nextSimulation); }// fin du else (number > 3) return doNext; }// fin du before simulation /** * Call after each simulation, compute criteria for last simulation. * * @param context plan context * @param lastSimulation storage used for next simulation * @return true if we must do next simulation, false to stop plan * @throws Exception */ public boolean afterSimulation(SimulationPlanContext context, SimulationStorage lastSimulation) throws Exception { boolean doNext = true; log.info("after simulation"); int number = lastSimulation.getParameter().getSimulationPlanNumber(); ResultStorage result = lastSimulation.getResultStorage(); /////***import the matrix of simulated data (here landings) from the simulation result MatrixND L = result.getMatrix(param_Population, MatrixCatchWeightPerStrategyMetPerZonePop.NAME); /////*** extract, sum, etc to obtain the same format/data as your observation matrix // useful methods to work on matrix : sumOverDim(), getSubMatrix(), reduce() ///////////////////Calcul du critere////////////////// log.info("calcul de la fonction objectif"); log.info("dim de L" + " " + Arrays.toString(L.getDim())); log.info("dim de obs" + " " + Arrays.toString(matrixDebarquement.getDim())); double obj = 0; // *** using the matrixIterator assumes that the observation and simulated // matrix are organised in the exact same way ( columns and rows corresponding in each of them) for (MatrixIterator g = L.iterator(); g.hasNext();) { g.next(); int[] dim = g.getCoordinates(); double obs = matrixDebarquement.getValue(dim); double simules = g.getValue(); obj += Math.pow(obs - simules, 2); } //ajoute le critere dans la table experiences experiences.getExperience(number).criteria = obj; //ecriture de la table historic exportHisto += experiences.getExperience(number).q1 + ";" + experiences.getExperience(number).q2 + ";" + experiences.getExperience(number).criteria + "\n"; ScriptUtil.writeString(exportHistoric, exportHisto); return doNext; }// fin du after simulation /** * Modify nextSimulation database with q1 and q2 in exp. * * @param exp * @param nextSimulation * @throws Exception */ protected void changeDB(Experience exp, SimulationStorage nextSimulation) throws Exception { // methode appelee dans before simualtion TopiaContext db = nextSimulation.getStorage().beginTransaction();//ouvrir un context pour modifier les donnees Population pop = (Population) db.findByTopiaId(param_Population.getTopiaId()); //reccupere la pop ciblee MatrixND c = pop.getCapturability(); // reccupere la matrice de capturabilité // *** that is where you explain how to fill the catchability matrix with q1 and q2 for (MatrixIterator i = c.iterator(); i.hasNext();) { i.next(); Object[] sem = i.getSemanticsCoordinates(); PopulationGroup group = (PopulationGroup) sem[0]; // *** exemple when q2 corresponds to the 12 first groups (groups 0 to 11) if (group.getId() < 12) { i.setValue(exp.q2); } else { i.setValue(exp.q1); } // *** exemple when it depends on seasons and groups /*PopulationSeasonInfo season = (PopulationSeasonInfo) sem[1]; if (season.getFirstMonth().after(Month.JULY) && group.getId() >=18){ //month >= aout && groupID >= 18 i.setValue(exp.q2); } else { i.setValue(exp.q1); } */ }//fin du for db.commitTransaction(); // effectue la modification db.closeContext(); // ferme le context } static public class Experiences { // cree la liste experiences ou sont stoqués q1,q2 et critere pour chaque simulation /** contains last simplex and potentialy 2 more simulation */ public List current = new ArrayList<>(); /** contains all experience done */ public List history = new ArrayList<>(); /** * return experience requested, if this experience doesn't exist * create it. * * @param i simulation number * @return experience with simulation number fixed if new experience * is returned */ public Experience getExperience(int i) { Experience result; if (i < history.size()) { result = history.get(i); } else { result = new Experience(); result.simNumber = i; history.add(i, result); current.add(result); } return result; }//fin de la definition de getExperience /** * @return the history */ public List getHistory() { return this.history; } }// fin de la creation des listes experiences static public class Experience implements Comparable { public int simNumber; public double criteria; public double q1; public double q2; /** * Permit to order experience, first is experience with smallest criteria */ public int compareTo(Experience other) { int result = Double.compare(this.criteria, other.criteria); return result; } } }