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.*; 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 calibration2paramtransfo3param implements SimulationPlan { /** to use log facility, just put in your code: log.info("..."); */ /////***class name to check here static private Log log = LogFactory.getLog(calibration2paramtransfo3param.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) public String param_nomfichier_debarquements = "C:/StephAncienPortableDisqueD/SimulationsISIS/Calibration/MerluLangoustineGdG/MerluCapturesFinalenkilos2001-2003SScol1.csv";//in row : time ; in columns : age or length group protected File debarquementsObserves; protected MatrixND matrixDebarquement; protected State state = State.STATE_INIT; public Experiences experiences = new Experiences(); public String [] necessaryResult = { ResultName.MATRIX_CATCH_WEIGHT_PER_STRATEGY_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 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_debarquements==null || "".equals(param_nomfichier_debarquements)){ debarquementsObserves = FileUtil.getFile(".*.csv", "fichier csv separateur ';'"); } else { debarquementsObserves = new File(param_nomfichier_debarquements); } 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 nbTrim = nbYear * 1; int[] dimMatrix = {nbTrim,nbGroup}; matrixDebarquement = MatrixFactory.getInstance().create(dimMatrix); //matrixDebarquement = MatrixFactory.getInstance().create(nbGroup); //matrixDebarquement = MatrixFactory.getInstance().create(new int[]{nbGroup}); // List groups = pop.getPopulationGroup(); // matrixDebarquement = MatrixFactory.getInstance().create(new List[]{groups}); //matrixDebarquement.importCSV(new FileReader(debarquementsObserves),new int []{0}); matrixDebarquement.importCSV(new FileReader(debarquementsObserves),new int []{0,0}); matrixDebarquement=matrixDebarquement.reduce(); log.info("MatrixDebarquement : " + matrixDebarquement); 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 MatrixND L2 = result.getMatrix(param_Population, ResultName.MATRIX_CATCH_WEIGHT_PER_STRATEGY_MET_PER_ZONE_POP); log.info("dim de L2" + " " + Arrays.toString(L2.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; //MatrixND L = L2.copy(); L = L.sumOverDim(0,12);//sum over dates 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 L2" + " " + Arrays.toString(L2.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(); ///////////////////Calcul du critere////////////////// log.info("calcul du critere"); log.info("dim de L" + " " + Arrays.toString(L.getDim())); log.info("dim de obs" + " " + Arrays.toString(matrixDebarquement.getDim())); double crit = 0; for ( MatrixIterator g = L.iterator(); g.hasNext();){ g.next(); //boucle sur les trimestres et les classes d age int [] dim = g.getCoordinates(); double obs = matrixDebarquement.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; } } }