Hello
Voici un script de calibration de la capturabilité quand il n y a qu un
seul paramètre qui semble assez efficace (précision = 1/32 du pas
initial en 11 simulations).
A adapter aux fonctions objectif de chacun (ici captures mensuelles en
poids) mais la partie algo (le before simulation) est générique.
(déposé aussi sur le wiki).
A+
Sigrid
--
Sigrid LEHUTA
~ ><> ~
Doctorante
Département Ecologie et Modèles pour l'Halieutique
IFREMER, rue de l'ile d'Yeu BP 21105
44311 Nantes Cedex 03
Tél : +33 (0)2 40 37 41 23 (interne : 8123)
package analyseplans;
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.Date;
import fr.ifremer.isisfish.entities.*;
import fr.ifremer.isisfish.simulator.AnalysePlan;
import fr.ifremer.isisfish.simulator.AnalysePlanContext;
import fr.ifremer.isisfish.simulator.SimulationParameter;
import fr.ifremer.isisfish.datastore.SimulationStorage;
import fr.ifremer.isisfish.datastore.ResultStorage;
/**
* Calibration_1parameter.java
* this program computes successive values for the catchability parameter in order to minimise the sum of squarres difference between observed and simulated catches.
* not generic ! comments indicate where the code must be adaptated to each fishery
*
* Created: 07/06/10
*
* @author <>
* @version $Revision: 3.2 $
*
* Last update: $Date: 2007/05/24 09:29:18 $
* by : $Author: bpoussin $
*/
public class Calibration_1param implements AnalysePlan {
/** to use log facility, just put in your code: log.info("..."); */
static private Log log = LogFactory.getLog(Calibration_1param.class);
@Doc("path and name of the file where parameter values and performance criteria are saved")
public String param_nomFileHisto = "Inputs_Anchois/CalibrationSardine/histo.csv";
File exportHistoric = new File (param_nomFileHisto);
protected String exportHisto = "";
@Doc("population concerned")
public Population param_Population = null;
@Doc("values for the 5 initial experiments in increasing order with a constant increment")
public Double param_M1 = 0.98e-4;// devient un parametre du plan d analyse
public Double param_M2 = 1.0e-4;// devient un parametre du plan d analyse
public Double param_M3 = 1.02e-4;// devient un parametre du plan d analyse
public Double param_M4 = 1.04e-4;// devient un parametre du plan d analyse
public Double param_M5 = 1.06e-4;// devient un parametre du plan d analyse
@Doc("path for the file that contains observed values used for calibration")
public String param_nomfichier_debarquements = "Inputs_Anchois/CalibrationSardine/DebarquementsSardine.csv";
protected File debarquementsObserves;
protected MatrixND matrixDebarquement;
protected String[] strNomsNous;
protected List<Strategy> strNous;
protected List<Double> q;
protected List<Double> Fobj = new ArrayList<Double>();
// TO DO : matrix containing simulated values to compare to observed values
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 _("Single Factor Optimization from Walters et al., 1991");
}
/**
* 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(AnalysePlanContext context) throws Exception {
// chargement de la matrice de valeurs observées
if (param_nomfichier_debarquements==null || "".equals(param_nomfichier_debarquements)){
debarquementsObserves = FileUtil.getFile(".*.csv", "fichier csv");
} else {
debarquementsObserves = new File(param_nomfichier_debarquements);
}
// creation of the vector of catchability values
q = new ArrayList<Double>();
q.add(param_M1);
q.add(param_M2);
q.add(param_M3);
q.add(param_M4);
q.add(param_M5);
// Creation of the matrix of observed values : TO DO
// Here the observed values are landings (weight) per strategy per month
int nbYear = context.getParam().getNumberOfYear();
List<Strategy> allStrategies = context.getParam().getStrategies();
// TO DO : names of strategies in apparition order in the file of observed landings
// also used to keep only certain strategies if necessary
strNomsNous = new String[]{"PelProfil1","PelProfil2","BolBasques","BolBretons"};
strNous = new ArrayList<Strategy>();
for(String strNomCurrent : strNomsNous){
for (Strategy s : allStrategies){
if(s.getName().equals(strNomCurrent)){
strNous.add(s);
}
}
}
// create the vector of dates at which catches are observed
int dateFin = nbYear*12;
List<Date> dates = new ArrayList<Date>();
for(int i=0; i<dateFin;i++){
dates.add(new Date(i));
}
TopiaContext db = context.getParam().getRegion().getStorage().beginTransaction();
Population pop = (Population)db.findByTopiaId(param_Population.getTopiaId());
// create the matrix of observed values which correspondant dimension
matrixDebarquement = MatrixFactory.getInstance().create(
"matDebarq",
new List[]{dates,strNous},
new String[]{"Dates","strategies"});
matrixDebarquement.importCSV(new FileReader(debarquementsObserves),new int []{0,0});
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
*/
public boolean beforeSimulation(AnalysePlanContext context, SimulationStorage nextSimulation) throws Exception {
boolean doNext = true;
int number = nextSimulation.getParameter().getAnalysePlanNumber();
log.info("before simulation" + number);
// TO DO : limit the number of experiments
if(number <31){
/** the rest of the code is generic, after the 5 first experiments were run, two new experiments are generated with values of q
at intervals equal to one-half the original interval on either side of the best value found in the initial set. The procedure is repeted, the interval beeing halved each time.
**/
double pasInitial = q.get(1)-q.get(0);
if (number <5) {
log.info("number<5");
changeDB(q.get(number), nextSimulation);
}else if(number == 5){
int min = compareTo(Fobj);
double newq = q.get(min);
double pas = pasInitial/2;
newq = newq - pas;
q.add(newq);
}else {
if(!isPair(number)){
int min = compareTo(Fobj);
double newq = q.get(min);
int puiss = (number-3)/2;
double pas = pasInitial/Math.pow(2.0,puiss);
newq = newq - pas;
q.add(newq);
}else{
int min = compareTo(Fobj.subList(0,number-1));
double newq = q.get(min);
int puiss = (number - 4)/2;
double pas = pasInitial/Math.pow(2.0,puiss);
newq = newq + pas;
q.add(newq);
}
}
changeDB(q.get(number), nextSimulation);
}else doNext = false;
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(AnalysePlanContext context, SimulationStorage lastSimulation) throws Exception {
int number = lastSimulation.getParameter().getAnalysePlanNumber();
//int number = context.getNumber();
System.out.println("after simulation" + number);
ResultStorage result = lastSimulation.getResultStorage();
// Pour sommer sur certaines classes les resultats :
MatrixND L2 = result.getMatrix(param_Population, ResultName.MATRIX_CATCH_WEIGHT_PER_STRATEGY_MET_PER_ZONE_POP);
List<Strategy> str = (List<Strategy>)L2.getSemantic(1);
System.out.println("import de la matrice des debarquements");
//pour reccuperer les index des strategies francaises
int comptSt = 0;
int comptTab = 0;
int[] elem = new int[strNous.size()];
for(Strategy st : str){
if (strNous.contains(st)){
elem[comptTab] = comptSt;
comptTab +=1;
}
comptSt += 1;
}
MatrixND L = L2.getSubMatrix(1,elem).copy();
L = L.sumOverDim(3);// sum sur les groupes
L = L.sumOverDim(2);// sum sur les metiers
L = L.sumOverDim(4);// sum sur les zones
L = L.reduce(); // supprime les dimensions de taille 0
///////////////////Calcul de la fonction objectif//////////////////
log.info("calcul de la fobj");
System.out.println("dim de L" + " " + Arrays.toString(L.getDim()));
System.out.println("matrice de debarquements : "+ L);
double crit = 0;
for ( MatrixIterator g = L.iterator(); g.hasNext();){
g.next();
//boucle sur les dates et les flottilles
Strategy s = (Strategy)g.getSemanticsCoordinates()[1];
Date dat = (Date)g.getSemanticsCoordinates()[0];
// int [] dim = g.getCoordinates();
double obs = matrixDebarquement.getValue(dat,s);
double simules = g.getValue();
double sum_pond = obs-simules;
crit += Math.pow(sum_pond, 2); // crit = crit + (obs-simules)^2
}// fin du for
System.out.println("critere " + number + " = " + crit );
Fobj.add(crit);
//ecriture de la table historic
exportHisto += q.get(number) +";"+ Fobj.get(number) + "\n";
org.nuiton.util.FileUtil.writeString(exportHistoric, exportHisto);
return true;
}// fin du after simulation
/**
* Modify nextSimulation database with q1 and q2 in exp.
* @param exp
* @param nextSimulation
* @throws Exception
*/
protected void changeDB(Double q, 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Ã
for (MatrixIterator i = c.iterator(); i.hasNext();){
i.next();
i.setValue(q);
}//fin du for
db.commitTransaction(); // effectue la modification
db.closeContext(); // ferme le context
}
public int compareTo(List<Double> listCrit) {
int N = listCrit.size();
System.out.println("taille de la liste :"+N);
int minimum=0;
for(int i=0; i<N; i++) {
System.out.println("i = "+i);
if (listCrit.get(i)<listCrit.get(minimum)){minimum = i;}
}
return minimum;
}
public int compare2(Double f1,Double f2){
int minimum=0;
if (f2<f1) minimum = 1;
return minimum;
}
public boolean isPair(int nombre){
return (nombre%2 == 0) ;
}
}
