[Isis-fish-devel] scipt simatrix

25 May 2007

      Salut Benjamin,
J'ai essayé de triturer deux ou 3 trucs sur le scipt simatrix. Ca gagne un peu de temps sur la base de Sigrid (4min contre 5min30), mais j'ai pas vérifié les résultats ( :-) ). Si ça t'intéresse je te joins le script. En gros l'astuce consiste à tester (dans totalFishingMortality et CatchRatePerStrategyMet) si effortPerZonePop  est non nul avant de lancer FishingMortality donc la boucle sur les groups sautent et on évite de refaire le calcul pleins de fois de l'effort (même si le résultat doit être dans le cache)
A+
Hilaire

-- 

Hilaire Drouineau
IFREMER, département EMH
Centre de Nantes - BP 21105 - 44311 Nantes Cedex 3
tel : +33 (0)2 40 37 42 12 - fax : +33 (0)2 40 37 40 75
email : hilaire.drouineau@ifremer.fr - web : http://www.ifremer.fr/emh
**********************************************************************

/* *##%
 * Copyright (C) 2006
 *     Code Lutin, Benjamin Poussin
 *
 * 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, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *##%*/

/* *
 * SiMatrix.java
 *
 * Created: 21 août 2006 15:53:01
 *
 * @author poussin
 * @version $Revision: 1.13 $
 *
 * Last update: $Date: 2007/03/01 17:16:12 $
 * by : $Author: bpoussin $
 */

package scripts;

import static org.codelutin.i18n.I18n._;
import static org.codelutin.i18n.I18n.n_;

import java.util.ArrayList;
import java.util.Collection;
import java.util.HashSet;
import java.util.Iterator;
import java.util.List;

import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;
import org.codelutin.math.matrix.MatrixFactory;
import org.codelutin.math.matrix.MatrixIterator;
import org.codelutin.math.matrix.MatrixND;
import org.codelutin.topia.TopiaContext;
import org.codelutin.topia.TopiaException;

import fr.ifremer.isisfish.IsisFishDAOHelper;
import fr.ifremer.isisfish.IsisFishException;
import fr.ifremer.isisfish.entities.Cell;
import fr.ifremer.isisfish.entities.EffortDescription;
import fr.ifremer.isisfish.entities.Gear;
import fr.ifremer.isisfish.entities.Metier;
import fr.ifremer.isisfish.entities.MetierSeasonInfo;
import fr.ifremer.isisfish.entities.Population;
import fr.ifremer.isisfish.entities.PopulationGroup;
import fr.ifremer.isisfish.entities.PopulationSeasonInfo;
import fr.ifremer.isisfish.entities.Selectivity;
import fr.ifremer.isisfish.entities.SetOfVessels;
import fr.ifremer.isisfish.entities.Strategy;
import fr.ifremer.isisfish.entities.StrategyMonthInfo;
import fr.ifremer.isisfish.entities.Zone;
import fr.ifremer.isisfish.entities.ZoneDAO;
import fr.ifremer.isisfish.simulator.SimulationContext;
import fr.ifremer.isisfish.types.Date;
import fr.ifremer.isisfish.types.Month;

/**
 * @author poussin
 *
 */

public class SiMatrix {

    /** to use log facility, just put in your code: log.info("..."); */
    static private Log log = LogFactory.getLog(SiMatrix.class);

    protected SimulationContext context = null;
    protected TopiaContext db = null;

    /**
     * Method used to get SiMatrix used for simulation
     * @param context context simulation
     * @return SiMatrix or null if no SiMatrix created for simulation
     */
    public static SiMatrix getSiMatrix(SimulationContext context) {
        SiMatrix result = (SiMatrix)context.getValue(SiMatrix.class.getName());
        return result;
    }

    private static void setSiMatrix(SimulationContext context, SiMatrix siMatrix) {
        context.setValue(SiMatrix.class.getName(), siMatrix);
    }

    /**
     *
     * @param context Simulation context
     * @param db TopiaContext with transaction opened. You must used this
     * TopiaContext and not used
     * SimulationContext.getSimulationStorage().getStorage()
     * @throws TopiaException
     */
    public SiMatrix(SimulationContext context) throws TopiaException {
        this.context = context;
        this.db = context.getDB();
        setSiMatrix(context, this);
    }

    /**
     * @return
     * @throws TopiaException
     */
    public List<Zone> getZones(Date date) throws TopiaException {
        ZoneDAO dao = IsisFishDAOHelper.getZoneDAO(db);
        List<Zone> result = dao.findAll();
        return result;
    }

    /**
     * @return
     * @throws TopiaException
     */
    public List<Population> getPopulations(Date date) throws TopiaException {
        List<Population> populations = new ArrayList<Population>();
        for (Population pop : context.getSimulationStorage().getParameter().getPopulations()) {
            Population tmp = (Population)db.findByTopiaId(pop.getTopiaId());
            populations.add(tmp);
        }
        return populations;
    }

    /**
     * @return
     * @throws TopiaException
     */
    public List<Strategy> getStrategies(Date date) throws TopiaException {
//        if (strategies == null) {
        List<Strategy> strategies = new ArrayList<Strategy>();
            for (Strategy str : context.getSimulationStorage().getParameter().getStrategies()) {
                Strategy tmp = (Strategy)db.findByTopiaId(str.getTopiaId());
                strategies.add(tmp);
            }
//        }
        return strategies;
    }

    public List<Metier> getMetiers(Date date) throws TopiaException {
//        if (metiers == null) {
            List<Metier> metiers = new ArrayList<Metier>();
            HashSet<Metier> tmp = new HashSet<Metier>();
            for (Strategy str : getStrategies(date)) {
                SetOfVessels sov = str.getSetOfVessels();
                for (EffortDescription effort : sov.getPossibleMetiers()) {
                    Metier metier = effort.getPossibleMetiers();
                    if (tmp.add(metier)) {
                        metiers.add(metier);
                    }
                }
            }
//        }
        return metiers;
    }

    /**
     * Retourne les metiers pratiqués par une Strategie à une date donnée
     * Un metier est pratiqué si le PropStrMet est différent de 0
     *
     * @param str
     * @param date
     * @return
     */
    public List<Metier> getMetiers(Strategy str, Date date) {
        StrategyMonthInfo info = str.getStrategyMonthInfo(date.getMonth());
        MatrixND props = info.getProportionMetier();

        List<Metier> result = new ArrayList<Metier>();

        for (MatrixIterator i=props.iterator(); i.hasNext();) {
            i.next();
            if (i.getValue() != 0) {
                Metier metier = (Metier)i.getSemanticsCoordinates()[0];
                result.add(metier);
            }
        }
        return result;
    }

    /**
     * Retourne la matrix Metier x Zone qui correspond au zone utilisé par
     * un métier pour une date donnée. Si la valeur de la matrice est 1 alors
     * la zone est utilisé par le métier, si elle vaut 0 alors elle n'est pas
     * utilisée.
     *
     * @param date
     * @return
     * @throws TopiaException
     */
    public MatrixND getMetierZone(Date date) throws TopiaException {
        List<Metier> metiers = getMetiers(date);
        List<Zone> zones = getZones(date);

        MatrixND result = MatrixFactory.getInstance().create(
                ResultName.MATRIX_METIER_ZONE,
                new List[]{metiers, zones},
                new String[]{n_("Metiers"), n_("Zones")});

        for (Metier metier : metiers) {
            Collection<Zone> zoneMetier = metier.getMetierSeasonInfo(date.getMonth()).getZone();
            for (Zone zone : zoneMetier) {
                result.setValue(metier, zone, 1);
            }
        }
        return result;
    }

    public MatrixND matrixPrice(Date date, Population pop) {
        List<PopulationGroup> groups = pop.getPopulationGroup();
        MatrixND result = MatrixFactory.getInstance().create(
                ResultName.MATRIX_PRICE,
                new List[]{groups},
                new String[]{n_("PopulationGroup")});

        for (PopulationGroup group : groups) {
            result.setValue(group, group.getPrice());
        }

        return result;
    }

   ///////////////////////////////////////////////////////////////////////////
    //
    // Toutes les methodes suivantes ne sont utiles que pour
    // matrixCatchPerStrategyMet
    //
    ///////////////////////////////////////////////////////////////////////////

    /**
     * Matrice des captures en nombre
     * dim [ Strategy x Metier x Classe x zonePop ]
     *
     * @param N l'abondance sous forme de matrice 2D [class x zone]
     * @param pop
     * @param date
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    public MatrixND matrixCatchPerStrategyMet(MatrixND N, Population pop,
            Date date, MatrixND matrixCatchRatePerStrategyMet) throws TopiaException, IsisFishException {
        List<PopulationGroup> groups = pop.getPopulationGroup();
        List<Zone> zones = pop.getPopulationZone();

        // on le passe en argument ce qui evite de le calculer 2 fois
//        MatrixND matrixCatchRatePerStrategyMet = matrixCatchRatePerStrategyMet(pop, date);
        MatrixND result = matrixCatchRatePerStrategyMet.copy();
        result.setName(ResultName.MATRIX_CATCH_PER_STRATEGY_MET);

        for(PopulationGroup group : groups) {
            MatrixND sub = result.getSubMatrix(2, group, 1);
            for(Zone zone : zones){
                MatrixND subsub = sub.getSubMatrix(3, zone, 1);
                double val = N.getValue(group, zone);
                subsub.mults(val);
            }
        }
        return result;
    }

    /**
     * Matrice des captures en poids
     * dim [ Strategy x Metier x Classe x zonePop ]
     *
     * @param pop
     * @param date
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    public MatrixND matrixCatchRatePerStrategyMet(Population pop, Date date) throws TopiaException, IsisFishException {
        List<Strategy> strategies = getStrategies(date);
        List<Metier> metiers = getMetiers(date);
        List<PopulationGroup> groups = pop.getPopulationGroup();
        List<Zone> zones = pop.getPopulationZone();

        MatrixND result = MatrixFactory.getInstance().create(
                ResultName.MATRIX_CATCH_RATE_PER_STRATEGY_MET,
                new List[]{strategies, metiers, groups, zones},
                new String[]{n_("Strategies"), n_("Metiers"), n_("Groups"), n_("Zones")});

        for (int s=0; s < strategies.size(); s++) {
            Strategy str = strategies.get(s);
            metiers = getMetiers(str, date);
            for (int m=0; m < metiers.size(); m++) {
                Metier metier = metiers.get(m);
                for (int z=0; z < zones.size(); z++) {
                    Zone zone = zones.get(z);
                    double effort=effortPerZonePop(str,metier,date,zone);
                    if (effort>0){
                      for (int g=0; g < groups.size(); g++) {
                          PopulationGroup group = groups.get(g);
                          double value = catchRatePerStrategyMet(str, metier, date, group, zone,effort);
                          result.setValue(str, metier, group, zone, value);
                      }
                    }
                    else {
                      for (int g=0; g < groups.size(); g++) {
                          PopulationGroup group = groups.get(g);
                          result.setValue(str, metier, group, zone, 0);
                      }
                    }
                }
            }
        }

//        for (Strategy str : strategies) {
//            List<Metier> metierStr = getMetiers(str, date);
//            for (Metier metier : metierStr) {
//                for (PopulationGroup group : groups) {
//                    for (Zone zone : zones) {
//                        double val = catchRatePerStrategyMet(str, metier, date, group, zone);
//                        result.setValue(str, metier, group, zone, val);
//                    }
//                }
//            }
//        }
        return result;
    }

    /**
     * @param str
     * @param metier
     * @param date
     * @param group
     * @param zone
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    private double catchRatePerStrategyMet(Strategy str, Metier metier, Date date, PopulationGroup group, Zone zone,double effort) throws TopiaException, IsisFishException {
        double totalFishingMortality = totalFishingMortality(date, group, zone);

        if(totalFishingMortality == 0){
            if(log.isDebugEnabled()) {log.debug("pas de totalFishingMortality pour (" + group + ", " + zone +")");}
            return 0;
        }

        double fishingMortality = fishingMortality(str, metier, date, group, zone,effort);
        double totalCatchRate = totalCatchRate(date, group, zone, totalFishingMortality);

        if(log.isDebugEnabled()) {
            log.debug(
            " totalFishingMortality=" + totalFishingMortality +
            " fishingMortality=" + fishingMortality +
            " totalCatchRate=" + totalCatchRate);
        }
        double result = fishingMortality / totalFishingMortality * totalCatchRate;

        return result;
    }

    /**
     * @param date
     * @param group
     * @param zone
     * @param totalFishingMortality
     * @return
     * @throws TopiaException
     */
    private double totalCatchRate(Date date, PopulationGroup group,
            Zone zone, double totalFishingMortality) throws TopiaException {
        double M = group.getNaturalDeathRate(zone) / Month.NUMBER_OF_MONTH;
        if(M == 0){
            // normalement il devrait y avoir de la mortalite naturelle
            if (log.isWarnEnabled()) {
                log.warn("Pas de mortalité naturelle pour: " + group);
            }
        }
        double F = totalFishingMortality;

        double result = 0;
        if( M != 0 || F != 0){
            result =  F/(F+M) * (1 - Math.exp(-(F+M)));
        }

        return result;
    }

    /**
     * @param date
     * @param group
     * @param zone
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    private double totalFishingMortality(Date date, PopulationGroup group, Zone zone) throws TopiaException, IsisFishException {
        List<Strategy> strategies = getStrategies(date);

        float result = 0;

        for(Strategy str : strategies){
            List<Metier> metierStr = getMetiers(str, date);
            for (Metier metier : metierStr) {
              double effort=effortPerZonePop(str,metier,date,zone);
                // TODO peut etre ne pas le faire si classe.pop n'est pas
                /// peche par le metier
                if (effort>0)
                  result += fishingMortality(str, metier, date, group, zone,effort);
            }
        }

        return result;
    }

    /**
     * @param str
     * @param metier
     * @param date
     * @param group
     * @param zone
     * @return
     * @throws IsisFishException
     */
    private double fishingMortality(Strategy str, Metier metier, Date date, PopulationGroup group, Zone zone,double effort) throws IsisFishException {
        Month month = date.getMonth();
        Population pop = group.getPopulation();
        Gear gear = metier.getGear();

        Selectivity selectivity = gear.getPopulationSelectivity(pop);

        double result = 0;

        if (selectivity != null) {
            MetierSeasonInfo infoMet = metier.getMetierSeasonInfo(month);
            double ciblage = infoMet.getTargetFactor(group);

            PopulationSeasonInfo infoPop = pop.getPopulationSeasonInfo(month);
            double capturability = infoPop.getCapturability(group);

            double coeff = selectivity.getCoefficient(pop, group, metier);

            if(log.isDebugEnabled()) {
                log.debug(
                        " strategy=" + str +
                        " metier=" + metier +
                        " ciblage=" + ciblage +
                        " capturabilite=" + capturability +
                        " selectivity=" + coeff +
                        " effort=" + effort);
            }

            result = coeff * capturability * ciblage * effort;
        }

        return result;
    }

    /**
     * @param str
     * @param metier
     * @param date
     * @param zone
     * @return
     */
    private double effortPerZonePop(Strategy str, Metier metier, Date date, Zone zonePop) {
        Month month = date.getMonth();
        Collection<Zone> zoneMet = metier.getMetierSeasonInfo(month).getZone();
        double inter = nbCellInter(zoneMet, zonePop);

        double effortPerStrategyPerCell =  effortPerStrategyPerCell(str, metier, date);

        if(log.isDebugEnabled()) {
            log.debug(
                    " strategy=" + str +
                    " metier=" + metier +
                    " inter=" + inter +
                    " effortPerStrategyPerCell=" + effortPerStrategyPerCell
            );
        }

        double result = effortPerStrategyPerCell * inter;
        return result;
    }

    /**
     * @param str
     * @param metier
     * @param date
     * @return
     */
    private double effortPerStrategyPerCell(Strategy str, Metier metier, Date date) {
        Month month = date.getMonth();
        StrategyMonthInfo smi = str.getStrategyMonthInfo(month);
        Collection<Zone> zones = metier.getMetierSeasonInfo(month).getZone();
        double nbCell = getCells(zones).size();

        if(nbCell == 0){
            // normalement il devrait y avoir des mailles, mais pour les
            // ancienne zone AuPort, il n'y en avait pas
            if(log.isWarnEnabled()) log.warn("Calcul d'une distance avec une zone sans maille");
            return 0;
        }

        double effortPerStrategy =  effortPerStrategyMet(str, metier, date);

        if(log.isDebugEnabled()) {
            log.debug(
                    " strategy=" + str +
                    " metier=" + metier +
                    " nbCell=" + nbCell +
                    " effortPerStrategy=" + effortPerStrategy
            );
        }

        double result = effortPerStrategy/nbCell;

        return result;
    }

    /**
     * @param str
     * @param metier
     * @param date
     * @return
     */
    private double effortPerStrategyMet(Strategy str, Metier metier, Date date) {
        Month month = date.getMonth();
        StrategyMonthInfo smi = str.getStrategyMonthInfo(month);

        double propSetOfVessels = str.getProportionSetOfVessels();
        int nbOfVessels = str.getSetOfVessels().getNumberOfVessels();
        double propStrMet = smi.getProportionMetier(metier);
        double effortPerVessel = effortPerStrategyPerVessel(str, metier, date);

        if(log.isDebugEnabled()) {
            log.debug(
                    " strategy=" + str +
                    " metier=" + metier +
                    " propSetOfVessels=" + propSetOfVessels +
                    " nbOfVessels=" + nbOfVessels +
                    " propStrMet=" + propStrMet +
                    " effortPerVessel=" + effortPerVessel
            );
        }

        double result = propSetOfVessels * nbOfVessels * propStrMet * effortPerVessel;

        return result;
    }

    /**
     * @param str
     * @param metier
     * @param date
     * @return
     */
    private double effortPerStrategyPerVessel(Strategy str, Metier metier, Date date) {
        Month month = date.getMonth();
        StrategyMonthInfo smi = str.getStrategyMonthInfo(month);
        int nbTrips = smi.getNumberOfTrips();
        double fishingTime = fishingTimePerTrip(str, metier, date);
        double stdEffortPerHour = stdEffortPerHour(str.getSetOfVessels(), metier);

        if(log.isDebugEnabled()) {
            log.debug(
                    " strategy=" + str +
                    " metier=" + metier +
                    " nbTrips=" + nbTrips +
                    " fishingTime=" + fishingTime +
                    " stdEffortPerHour=" + stdEffortPerHour
            );
        }

        double result = nbTrips * fishingTime * stdEffortPerHour;

        return result;
    }

    /**
     * Used in GravityModel too
     *
     * @param str
     * @param metier
     * @param date
     * @return
     */
    protected double fishingTimePerTrip(Strategy str, Metier metier, Date date) {
        Month month = date.getMonth();
        StrategyMonthInfo smi = str.getStrategyMonthInfo(month);
        Collection<Zone> zone = metier.getMetierSeasonInfo(month).getZone();

        if (zone == null) {
            if(log.isWarnEnabled()) log.warn(
                    "missing zone for metier =" + metier +
                    " for month" + month
            );
        }

        double tripDuration =  smi.getTripType().getTripDuration().getHour();
        double travelTime = travelTimePerTrip(str.getSetOfVessels(), zone);

        double result = tripDuration - travelTime;

        if (result < 0 ) {
            if(log.isWarnEnabled()) log.warn(
            " strategy=" + str +
            " metier=" + metier +
            " tripDuration=" + tripDuration +
            " travelTime=" + travelTime
            );
        }
        return result;
    }

    /**
     *
     * @param setOfVessels
     * @param zone
     * @return
     */
    protected double travelTimePerTrip(SetOfVessels sov, Collection<Zone> zoneMetier) {
        Cell maille = sov.getPort().getCell();
        double result =
            2 * distance(zoneMetier, maille) / sov.getVesselType().getSpeed();

        return result;
    }

    /**
     * @param zoneMetier
     * @param maille
     * @return
     */
    private double distance(Collection<Zone> zones, Cell cell) {
        double result = 0;
        List<Cell> cells = getCells(zones);

        if(cells.size() == 0){
            // normalement il devrait y avoir des mailles, mais pour les
            // ancienne zone AuPort, il n'y en avait pas
            if(log.isWarnEnabled()) {
                log.warn("Calcul d'une distance avec une zone sans maille");
            }
            return 0;
        }

        for(Cell c : cells){
            result += distance(c, cell);
        }

        if(log.isDebugEnabled()) {
            log.debug(" result=" + result + " nbMaille="+cells.size());
        }

        result = result / (double)cells.size();
        return result;
    }

    /**
     * @param c
     * @param cell
     * @return
     */
    private double distance(Cell m1, Cell m2) {
        double earthRadius = 6378.388;
        double p = 180/Math.PI;

        if(log.isDebugEnabled()) log.debug("p: " + p);

        double m1lat = m1.getLatitude();
        double m2lat = m2.getLatitude();
        double m1lon = m1.getLongitude();
        double m2lon = m2.getLongitude();

        if(log.isDebugEnabled()) log.debug(
        " m1lat=" + m1lat +
        " m2lat=" + m2lat +
        " m1lon=" + m1lon +
        " m2lonx=" + m2lon
        );

        double m1lat_div_p = m1lat/p;
        double m2lat_div_p = m2lat/p;
        double m1lon_div_p = m1lon/p;
        double m2lon_div_p = m2lon/p;

        if(log.isDebugEnabled()) log.debug(
        " m1lat_div_p=" + m1lat_div_p +
        " m2lat_div_p=" + m2lat_div_p +
        " m1lon_div_p=" + m1lon_div_p +
        " m2lon_div_p=" + m2lon_div_p
        );

        double sin_m1lat_div_p = Math.sin(m1lat_div_p);
        double sin_m2lat_div_p = Math.sin(m2lat_div_p);
        double cos_m1lat_div_p = Math.cos(m1lat_div_p);
        double cos_m2lat_div_p = Math.cos(m2lat_div_p);

        if(log.isDebugEnabled()) log.debug(
        " sin_m1lat_div_p=" + sin_m1lat_div_p +
        " sin_m2lat_div_p=" + sin_m2lat_div_p +
        " cos_m1lat_div_p=" + cos_m1lat_div_p +
        " cos_m2lat_div_p=" + cos_m2lat_div_p
        );

        double cos_m1lon_div_p_minus_m2lon_div_p = Math.cos(m1lon_div_p - m2lon_div_p);

        if(log.isDebugEnabled()) log.debug(
        " cos_m1lon_div_p_minus_m2lon_div_p=" + cos_m1lon_div_p_minus_m2lon_div_p);

        double acos = Math.acos(
                    sin_m1lat_div_p
                    * sin_m2lat_div_p
                    +
                    cos_m1lat_div_p
                    * cos_m2lat_div_p
                    * cos_m1lon_div_p_minus_m2lon_div_p
                    );

        if(log.isDebugEnabled()) log.debug(" acos=" + acos);

        double result = earthRadius * acos;

        return result;
    }

    /**
     * @param setOfVessels
     * @param metier
     * @return
     */
    private double stdEffortPerHour(SetOfVessels sov, Metier metier) {
        double result = 0;
        EffortDescription ed = sov.getPossibleMetiers(metier);
        if(ed != null){
            double fstd = metier.getGear().getStandardisationFactor();
            double val = fstd * ed.getFishingOperation() * ed.getGearsNumberPerOperation();
            result = val;
        }
        result = result/24; // 24 heures

        return result;
    }

    /**
     * used here and in Rule (CantonnementPreSimu)
     *
     * @param zones
     * @return
     */
    public List<Cell> getCells(Collection<Zone> zones) {
        List<Cell> result = new ArrayList<Cell>();
        for (Zone zone : zones) {
            result.addAll(zone.getCell());
        }
        return result;
    }

    /**
     * used here and in Rule (CantonnementPreSimu)
     *
     * @param zoneMet
     * @param zonePop
     * @return
     */
    public int nbCellInter(Collection<Zone> zoneMet, Zone zonePop) {
        List<Cell> cells = getCells(zoneMet);
        List<Cell> tmp = new ArrayList<Cell>(cells);
        tmp.retainAll(zonePop.getCell());
        return tmp.size();
    }

    ///////////////////////////////////////////////////////////////////////////
    //
    //
    //
    ///////////////////////////////////////////////////////////////////////////

    /**
     * @param N
     * @param pop
     * @param date
     * @return
     * @throws IsisFishException
     * @throws TopiaException
     */
    public MatrixND matrixAbundance(MatrixND N, Population pop, Date date) throws TopiaException, IsisFishException {
        List<PopulationGroup> groups = pop.getPopulationGroup();
        List<Zone> zones = pop.getPopulationZone();

        MatrixND result = MatrixFactory.getInstance().create(
                ResultName.MATRIX_ABUNDANCE,
                new List[]{groups, zones},
                new String[]{n_("Groups"), n_("Zones")});

        for (int g=0; g < groups.size(); g++) {
            PopulationGroup group = groups.get(g);
            for (int z=0; z < zones.size(); z++) {
                Zone zone = zones.get(z);
                double value = survivalRate(date, group, zone);
                double n = N.getValue(g, z);
                value *= n;
                result.setValue(g, z, value);
            }
        }

//        for(PopulationGroup group : groups){
//            for(Zone zone : zones){
//                double val = survivalRate(date, group, zone);
//                val *= N.getValue(group, zone);
//                result.setValue(group, zone, val);
//            }
//        }

        return result;
    }

    /**
     * @param date
     * @param group
     * @param zone
     * @return
     * @throws IsisFishException
     * @throws TopiaException
     */
    private double survivalRate(Date date, PopulationGroup group, Zone zone) throws TopiaException, IsisFishException {
        double F = totalFishingMortality(date, group, zone);  // rem perf: totalFishingMortality a deja ete calculé
        double M = group.getNaturalDeathRate(zone)/(double)Month.NUMBER_OF_MONTH;
        double result = (double)Math.exp(-(F+M));

        return result;
    }

    ///////////////////////////////////////////////////////////////////////////
    //
    //
    //
    ///////////////////////////////////////////////////////////////////////////

    /**
     * @param n
     * @param pop
     * @param date
     * @return
     */
    public MatrixND matrixBiomass(MatrixND N, Population pop, Date date) {
        List<PopulationGroup> groups = N.getSemantics(0);
        List<Zone> zones = N.getSemantics(1);

        MatrixND result = MatrixFactory.getInstance().create(
                ResultName.MATRIX_BIOMASS,
                new List[]{groups, zones},
                new String[]{n_("Groups"), n_("Zones")});

        for (int g=0; g < groups.size(); g++) {
            PopulationGroup group = groups.get(g);
            double meanWeight = group.getMeanWeight();
            for (int z=0; z < zones.size(); z++) {
                Zone zone = zones.get(z);
                double n = N.getValue(group, zone);
                double value = n * meanWeight;
                result.setValue(group, zone, value);
            }
        }

//        for(PopulationGroup group : groups){
//            double meanWeight = group.getMeanWeight();
//            for(Zone zone : zones){
//                double val = N.getValue(group, zone) * meanWeight;
//                result.setValue(group, zone, val);
//            }
//        }

        return result;
    }

    ///////////////////////////////////////////////////////////////////////////
    //
    //
    //
    ///////////////////////////////////////////////////////////////////////////

    /**
     * @param date
     * @return
     * @throws TopiaException
     */
    public MatrixND matrixEffortPerStrategyMet(Date date) throws TopiaException {
        List<Strategy> strategies = getStrategies(date);
        List<Metier> metiers = getMetiers(date);

        MatrixND result = MatrixFactory.getInstance().create(
                ResultName.MATRIX_EFFORT_PER_STRATEGY_MET,
                new List[]{strategies, metiers},
                new String[]{n_("Strategies"), n_("Metiers")});

        for (int s=0; s < strategies.size(); s++) {
            Strategy str = strategies.get(s);
            metiers = getMetiers(str, date);
            for (int m=0; m < metiers.size(); m++) {
                Metier metier = metiers.get(m);
                double value = effortPerStrategyMet(str, metier, date);
                result.setValue(str, metier, value);
            }
        }

//        for(Strategy str : strategies){
//            List<Metier> metierStr = getMetiers(str, date);
//            for(Metier metier : metierStr) {
//                double val = effortPerStrategyMet(str, metier, date); // rem perf: effortPerStrategyMet a deja ete calculé
//                result.setValue(str, metier, val);
//            }
//        }

        return result;
    }

    ///////////////////////////////////////////////////////////////////////////
    //
    //
    //
    ///////////////////////////////////////////////////////////////////////////

    /**
     * @param pop
     * @param date
     * @return
     */
    public MatrixND matrixCatchWeightPerStrategyMet(Population pop,
            Date date, MatrixND matrixCatchPerStrategyMet) {
      List<PopulationGroup> groups = pop.getPopulationGroup();

      MatrixND result = matrixCatchPerStrategyMet.copy();
      result.setName(ResultName.MATRIX_CATCH_WEIGHT_PER_STRATEGY_MET);

      for(PopulationGroup group : groups){
          MatrixND sub = result.getSubMatrix(2, group, 1);
          double meanWeight = group.getMeanWeight();
          sub.mults(meanWeight);
      }

      return result;
    }

    /**
     * @param pop
     * @param date
     * @return
     */
    public MatrixND matrixDiscardWeightPerStrategyMet(Population pop,
            Date date, MatrixND matrixDiscardPerStrategyMet) {
      List<PopulationGroup> groups = pop.getPopulationGroup();

      MatrixND result = matrixDiscardPerStrategyMet.copy();
      result.setName(ResultName.MATRIX_DISCARDS_WEIGHT_PER_STR_MET);

      for(PopulationGroup group : groups){
          MatrixND sub = result.getSubMatrix(2, group, 1);
          double meanWeight = group.getMeanWeight();
          sub.mults(meanWeight);
      }

      return result;
    }

}

[Isis-fish-devel] scipt simatrix

Hilaire DROUINEAU