package simulationplans; import static org.nuiton.i18n.I18n._; import org.apache.commons.logging.Log; import org.apache.commons.logging.LogFactory; import scripts.ResultName; import java.io.*; import java.util.*; import org.nuiton.math.matrix.*; import org.nuiton.topia.*;// pour pouvoir utiliser la methode StringUtil.toDouble() import org.nuiton.util.*;// pour pouvoir utiliser la methode StringUtil.toDouble() import fr.ifremer.isisfish.*; import fr.ifremer.isisfish.types.*; import fr.ifremer.isisfish.simulator.SimulationContext; import fr.ifremer.isisfish.types.TimeStep; import fr.ifremer.isisfish.entities.Metier; import fr.ifremer.isisfish.entities.*; import fr.ifremer.isisfish.simulator.SimulationPlan; import fr.ifremer.isisfish.simulator.SimulationPlanContext; import fr.ifremer.isisfish.simulator.SimulationParameter; import fr.ifremer.isisfish.datastore.SimulationStorage; import fr.ifremer.isisfish.datastore.ResultStorage; /** * /////***File must be copied in isis-database-3/ analyseplans/ * File name : CalibrationEspeceq1q2Export.java * /////*** File name could be modified if needed * /////*** BUT class name must be identical to file name (without the extention ".java") see below * * Created: 17 septembre 2007 * * @author <> * @version $Revision: 1.27 $ * * Last update: $Date: 2007/05/24 09:29:18 $ * by : $Author: bpoussin $ */ /////***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 : "CalibrationEspeceq1q2Export" public class calibrationCaptDebMerlu2paramtransfo3param implements SimulationPlan { /** to use log facility, just put in your code: log.info("..."); */ /////***class name to check here static private Log log = LogFactory.getLog(calibrationCaptDebMerlu2paramtransfo3param.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") where q1, q2 and criteria computed at each simulation are written /////***Attention : before beginning a new calibration rename your eventual old file Historic.csv otherwise it will be lost File exportHistoric = new File ("C:/StephAncienPortableDisqueD/SimulationsISIS/Calibration/MerluLangoustineGdG/historicCalibMerluLangGdG.csv"); protected String exportHisto = ""; //initial points of the simplex public Population param_Population = null; public String param_M1 = "1.1e-5;9.5e-7";// devient un parametre du plan d analyse public String param_M2 = "1.3e-5;8.9e-7";// devient un parametre du plan d analyse public String param_M3 = "9.9e-6;9.9e-7";// devient un parametre du plan d analyse //public String param_pas = "1e-5";// devient un parametre du plan d analyse ///// ***put here the path and name of the file containing the data on which you calibrate your fichery ( here observed catches) // debarquements sur 2001-2002 pour tous les metiers public String param_nomfichier_debarquements200102 = "C:/StephAncienPortableDisqueD/SimulationsISIS/Calibration/MerluLangoustineGdG/MerluDebarquementsFinalenkilos2001-2002SScol1.csv";//in row : time ; in columns : age or length group protected File debarquementsObserves200102; protected MatrixND matrixDebarquement200102; // debarquements sur 2003 pour tous les metiers sauf langoustiniers et captures pour les langoustiniers public String param_nomfichier_capturesdebarquements2003 = "C:/StephAncienPortableDisqueD/SimulationsISIS/Calibration/MerluLangoustineGdG/MerluCapturesDebFinalenkilos2003SScol1.csv";//in row : time ; in columns : age or length group protected File capturesdebarquementsObserves2003; protected MatrixND matrixCapturesDebarquement2003; protected State state = State.STATE_INIT; public Experiences experiences = new Experiences(); public String [] necessaryResult = { ResultName.MATRIX_LANDING_PER_MET, ResultName.MATRIX_DISCARDS_WEIGHT_PER_STR_MET_PER_ZONE_POP }; 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 Merlu (q1,q2 en moyennenant q1 q2 pour obtenir un parametre q3 pour les classes intermediares, sur 2001-2002 avec les debarquements dans GdG et 2003 les debarquements auxquels on a ajoute les rejets des langoustiniers using variable step Simplex method (Walters): user gives a file of observations (here catches) by time step and group (.csv), 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 simulation La simulation pour lequel on utilise cette regle */ public void init(SimulationPlanContext context) throws Exception { if (param_nomfichier_debarquements200102==null || "".equals(param_nomfichier_debarquements200102)){ debarquementsObserves200102 = FileUtil.getFile(".*.csv", "fichier csv separateur ';'"); } else { debarquementsObserves200102 = new File(param_nomfichier_debarquements200102); } if (param_nomfichier_capturesdebarquements2003==null || "".equals(param_nomfichier_capturesdebarquements2003)){ capturesdebarquementsObserves2003 = FileUtil.getFile(".*.csv", "fichier csv separateur ';'"); } else { capturesdebarquementsObserves2003 = new File(param_nomfichier_capturesdebarquements2003); } int nbYear = context.getParam().getNumberOfYear(); TopiaContext db = context.getParam().getRegion().getStorage().beginTransaction(); Population pop = (Population)db.findByTopiaId(param_Population.getTopiaId()); /////*** specify dimention of the matrix containning observations (observed landings for instance) /////*** numbers of group/columns : could be equal to your number of classes in ISIS but may also be different if your had only aggregated data int nbGroup = 65 ; //int[] nbGroup = {55}; /////*** enter number of observation per year (if you have observation by quarter put 4) / lines of the observations file int nbyearDeb = 2; //ici 2001 et 2002 int[] dimMatrix = {nbyearDeb,nbGroup}; matrixDebarquement200102 = MatrixFactory.getInstance().create(dimMatrix); matrixDebarquement200102.importCSV(new FileReader(debarquementsObserves200102),new int []{0,0}); matrixDebarquement200102=matrixDebarquement200102.reduce(); log.info("matrixDebarquement200102 : " + matrixDebarquement200102); int nbyearDebCapt = 1; //ici 2003 int[] dimMatrix2 = {nbyearDebCapt,nbGroup}; matrixCapturesDebarquement2003 = MatrixFactory.getInstance().create(dimMatrix2); matrixCapturesDebarquement2003.importCSV(new FileReader(capturesdebarquementsObserves2003),new int []{0,0}); matrixCapturesDebarquement2003=matrixCapturesDebarquement2003.reduce(); log.info("matrixCapturesDebarquement2003 : " + matrixCapturesDebarquement2003); db.closeContext(); } /** * 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 */ double g1; double g2; double worst1; double worst2; public boolean beforeSimulation(SimulationPlanContext context, SimulationStorage nextSimulation) throws Exception { boolean doNext = true; boolean doBoucle = true; log.info("before simulation"); int number = nextSimulation.getParameter().getSimulationPlanNumber(); //if (number>4) { // doNext = false; //} else 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(1).q1 + experiences.current.get(2).q1) / 2.0 - experiences.current.get(3).q1; q2 = experiences.getExperience(number-1).q2 + (experiences.current.get(1).q2 + experiences.current.get(2).q2) / 2.0 - experiences.current.get(3).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(1).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 nextSimulation 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(); /////*** Simulated catches are cumulated over strategies, areas and metiers // pour les debarquements de 2001 et 2002 MatrixND L = result.getMatrix(param_Population, ResultName.MATRIX_LANDING_PER_MET); //log.info("dim de L" + " " + Arrays.toString(L.getDim())); /////*** If some strategies, metiers or areas must not be included in the cumulated catches see below else put "//" at the beginning of line 297 and replace "L2" by "L" at line 289 /////*** the strategies or so to exclude must be at the begining or end of the list you enter in the parameters interface /////*** fill the line MatrixND L = L2.getSubMatrix(a,b,c).copy(); as follow : /////*** a = 1 if you want to exclude strategies, 2 if metiers and 4 if areas /////*** b = indice of the first object considered /////*** c = number of object to keep after b, b included //MatrixND L = L2.getSubMatrix(1,2,4).copy(); // for instance : keep strategies 2,3,4,5 // si tu veux toutes les strategies: //MatrixND L = L2; L = L.getSubMatrix(0,0,24); //pour prendre les deux premieres annees de simu (2001 et 2002) L = L.sumOverDim(0,12);//sum over dates //pour agreger par annee (ici on calibre sur les captures totales de l'annee L = L.sumOverDim(1);// sum over strategies L = L.sumOverDim(2);// sum over metiers L = L.sumOverDim(4);// sum over zones // log.info("sommes sur les strategies, metiers et zones faites"); // log.info("sous matrice extraite"); // log.info("dim de L" + " " + Arrays.toString(L.getDim())); /////*** if observations are cumulated by groups else put "//" at the beginning of the following line //L = L.sumOverDim(3); /////*** else if number of groups in your database differs from number of groups in observations ... write me an e-mail ////*** if observations are cumulated over quarter a = 3, over year a = 12, else put "//" at the beginning of the following line // L = L.sumOverDim(0,a); L = L.reduce(); // pour 2003 : on prend tous les debarquements et on rajoute les rejets des langoustiniers en 2003 MatrixND L3 = result.getMatrix(param_Population, ResultName.MATRIX_LANDING_PER_MET); // log.info("dim de L3" + " " + Arrays.toString(L3.getDim())); L3 = L3.getSubMatrix(0,23,12); //pour prendre la troisieme annee 2003 //pour les rejets MatrixND L5 = result.getMatrix(param_Population, ResultName.MATRIX_DISCARDS_WEIGHT_PER_STR_MET_PER_ZONE_POP); // log.info("dim de L5" + " " + Arrays.toString(L5.getDim())); L5 = L5.getSubMatrix(0,23,12); ///pour prendre la troisieme annee 2003 // log.info("--------- Classe L5 : " +L3.getClass()); // log.info("--------- Classe L6 : " +L5.getClass()); // parcourir L4 et L6 pour ajouter dans L4 les valeurs de L6 si Metiers langoustiniers MatrixIterator j = L3.iterator(); for (MatrixIterator i = L5.iterator(); i.hasNext();) { i.next(); j.next(); Object[] coordonneesL5 = i.getSemanticsCoordinates(); Object[] coordonneesL3 = j.getSemanticsCoordinates(); // log.info("--------- coordonneesL5[2] :" +coordonneesL5[2]); String metL5 = (String) coordonneesL5[2]; String metL3 = (String) coordonneesL3[2]; if ((metL5.equals("metier lang simp Sud")&&metL3.equals("metier lang simp Sud"))||(metL5.equals("metier lang simp Nord")&&metL3.equals("metier lang simp Nord"))||(metL5.equals("metier lang simple")&&metL3.equals("metier lang simple"))||(metL5.equals("metier lang jum Nord")&&metL3.equals("metier lang jum Nord"))||(metL5.equals("metier lang jum Sud")&&metL3.equals("metier lang jum Sud"))||(metL5.equals("metier lang jum")&&metL3.equals("metier lang jum"))) { double value = i.getValue() + j.getValue(); // log.info("VALUE" + value); j.setValue(value); } } L3 = L3.sumOverDim(0,12);//sum over dates (agregation sur l'annee) L3 = L3.sumOverDim(1);// sum over strategies L3 = L3.sumOverDim(2);// sum over metiers L3 = L3.sumOverDim(4);// sum over zones // log.info("sommes sur les strategies, metiers et zones faites"); // log.info("sous matrice extraite"); // log.info("dim de L3" + " " + Arrays.toString(L3.getDim())); L3 = L3.reduce(); // ///////////////////Calcul du critere////////////////// log.info("calcul du critere"); double crit = 0; //sur 2001 et 2002 for ( MatrixIterator g = L.iterator(); g.hasNext();){ g.next(); //boucle sur les annees et les groupes int [] dim = g.getCoordinates(); double obs = matrixDebarquement200102.getValue(dim); double simules = g.getValue(); crit += Math.pow(obs-simules, 2); // crit = crit + (obs-simules)^2 }// fin du for // sur 2003 for ( MatrixIterator g = L3.iterator(); g.hasNext();){ g.next(); //boucle sur les annees (il n'y en a qu'une ... j'espere que ca ne posera pas de pb et les groupes int [] dim = g.getCoordinates(); double obs = matrixCapturesDebarquement2003.getValue(dim); double simules = g.getValue(); crit += Math.pow(obs-simules, 2); // crit = crit + (obs-simules)^2 }// fin du for log.info("critere " + number + " = " + crit ); //ajoute le critere dans la table experiences experiences.getExperience(number).criteria = crit; //ecriture de la table historic exportHisto += experiences.getExperience(number).q1 +";"+ experiences.getExperience(number).q2 +";"+ experiences.getExperience(number).criteria + "\n"; org.nuiton.util.FileUtil.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 capturabilite //log.info("Pour cette simulation : q1 = " + exp.q1 + ";" + "q2 = " + exp.q2 ); /////*** 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]; PopulationSeasonInfo season = (PopulationSeasonInfo)sem[1]; ////*** exemple when q2 corresponds to the 21 first groups (+ de 50% d'immatures) if (group.getId() < 16){ i.setValue(exp.q1); }else if (group.getId() < 30){ i.setValue((exp.q1+exp.q2)/2); }else { i.setValue(exp.q2); } /////*** exemple when it depends on seasons and groups /* 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 stoque 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 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(Object arg0) { Experience other = (Experience)arg0; int result = Double.compare(this.criteria, other.criteria); return result; } } }