/*
 * #%L
 * IsisFish data
 * %%
 * Copyright (C) 2006 - 2011 Ifremer, 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, see
 * <http://www.gnu.org/licenses/gpl-2.0.html>.
 * #L%
 */

package scripts;

import static org.nuiton.i18n.I18n.n_;

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

import org.apache.commons.logging.Log;
import org.apache.commons.logging.LogFactory;
import org.nuiton.math.matrix.MatrixFactory;
import org.nuiton.math.matrix.MatrixIterator;
import org.nuiton.math.matrix.MatrixND;
import org.nuiton.topia.TopiaContext;
import org.nuiton.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.TimeStep;
import fr.ifremer.isisfish.types.Month;
import fr.ifremer.isisfish.datastore.ResultStorage;

/**
 * SiMatrix.java
 *
 * Created: 21 aout 2006 15:53:01
 *
 * @author poussin
 * @version $Revision: 1.18 $
 *
 * Last update: $Date: 2007-11-02 17:53:20 $
 * by : $Author: bpoussin $
 */
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
     * @throws TopiaException
     */
    public static SiMatrix getSiMatrix(SimulationContext context)
            throws TopiaException {
        SiMatrix result = (SiMatrix) context.getValue(SiMatrix.class.getName());
        if (result == null) {
            result = new SiMatrix(context);
        }
        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(TimeStep step) throws TopiaException {
        ZoneDAO dao = IsisFishDAOHelper.getZoneDAO(db);
        List<Zone> result = dao.findAll();
        return result;
    }

    /**
     * @return
     * @throws TopiaException
     */
    public List<Population> getPopulations(TimeStep step) 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(TimeStep step) 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(TimeStep step) throws TopiaException {
        //        if (metiers == null) {
        List<Metier> metiers = new ArrayList<Metier>();
        HashSet<Metier> tmp = new HashSet<Metier>();
        for (Strategy str : getStrategies(step)) {
            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 pratiques par une Strategie a une step donnee Un
     * metier est pratique si le PropStrMet est different de 0
     * 
     * @param str
     * @param step
     * @return
     */
    public List<Metier> getMetiers(Strategy str, TimeStep step) {
        StrategyMonthInfo info = str.getStrategyMonthInfo(step.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 utilise par un
     * metier pour une step donnee. Si la valeur de la matrice est 1 alors la
     * zone est utilise par le metier, si elle vaut 0 alors elle n'est pas
     * utilisee.
     * 
     * @param step
     * @return
     * @throws TopiaException
     */
    public MatrixND getMetierZone(TimeStep step) throws TopiaException {
        List<Metier> metiers = getMetiers(step);
        List<Zone> zones = getZones(step);

        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(
                    step.getMonth()).getZone();
            for (Zone zone : zoneMetier) {
                result.setValue(metier, zone, 1);
            }
        }
        return result;
    }

    public MatrixND matrixPrice(TimeStep step, 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
    //
    ///////////////////////////////////////////////////////////////////////////

    public MatrixND matrixCatchWeightPerStrategyMetPerZoneMet(TimeStep step,
            Population pop, MatrixND matrixCatchPerStrategyMetPerZoneMet)
            throws TopiaException, IsisFishException {
        List<PopulationGroup> groups = pop.getPopulationGroup();

        MatrixND result = matrixCatchPerStrategyMetPerZoneMet.copy();
        result.setName(ResultName.MATRIX_CATCH_WEIGHT_PER_STRATEGY_MET_PER_ZONE_MET);

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

        return result;
    }

    public MatrixND matrixCatchPerStrategyMetPerZoneMet(MatrixND N,
            Population pop, TimeStep step) throws TopiaException, IsisFishException {

        MatrixND matrixFishingMortalityPerCell = matrixFishingMortalityPerCell(
                step, pop);
        MatrixND matrixCatchRatePerStrategyMetPerCell = matrixCatchRatePerStrategyMetPerCell(
                pop, step, matrixFishingMortalityPerCell);
        MatrixND matrixCatchPerStrategyMetPerCell = matrixCatchPerStrategyMetPerCell(
                N, pop, step, matrixCatchRatePerStrategyMetPerCell);

        List<Strategy> strategies = getStrategies(step);
        List<Metier> metiers = getMetiers(step);
        List<PopulationGroup> groups = (List<PopulationGroup>)matrixCatchPerStrategyMetPerCell
                .getSemantic(2);
        List<Zone> zones = getZones(step);

        Set<Cell> cellPops = new HashSet(matrixCatchPerStrategyMetPerCell
                .getSemantic(4));

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

        // matrice temporaire ou les zones pops sont sommees
        MatrixND tmp = matrixCatchPerStrategyMetPerCell.sumOverDim(3);
        tmp = tmp.reduceDims(3);

        for (int s = 0; s < strategies.size(); s++) {
            Strategy str = strategies.get(s);
            for (int g = 0; g < groups.size(); g++) {
                PopulationGroup group = groups.get(g);
                for (int m = 0; m < metiers.size(); m++) {
                    Metier metier = metiers.get(m);
                    MetierSeasonInfo infoMet = metier.getMetierSeasonInfo(step
                            .getMonth());
                    Collection<Zone> zoneMet = infoMet.getZone();
                    for (Zone z : zoneMet) {
                        double value = 0;
                        List<Cell> cells = z.getCell();
                        for (int c = 0; c < cells.size(); c++) {
                            Cell cell = cells.get(c);
                            if (cellPops.contains(cell)) {
                                // les cells de la matrice sont les cells des 
                                // zones pops, donc seul les intersections avec
                                // les cells des metiers sont des cells valides
                                value += tmp.getValue(str, metier, group, cell);
                            }
                        }
                        result.setValue(str, metier, group, z, value);
                    }
                }
            }
        }

        return result;
    }

    public MatrixND matrixCatchPerStrategyMetPerZonePop(MatrixND N,
            Population pop, TimeStep step) throws TopiaException, IsisFishException {
        MatrixND matrixFishingMortalityPerCell = matrixFishingMortalityPerCell(
                step, pop);
        MatrixND matrixCatchRatePerStrategyMetPerCell = matrixCatchRatePerStrategyMetPerCell(
                pop, step, matrixFishingMortalityPerCell);
        MatrixND matrixCatchPerStrategyMetPerCell = matrixCatchPerStrategyMetPerCell(
                N, pop, step, matrixCatchRatePerStrategyMetPerCell);

        // on somme sur les cellules
        MatrixND result = matrixCatchPerStrategyMetPerCell.sumOverDim(4);
        result = result.reduceDims(4);
        result.setName(ResultName.MATRIX_CATCH_PER_STRATEGY_MET_PER_ZONE_POP);

        return result;
    }

    public MatrixND matrixCatchWeightPerStrategyMetPerZonePop(TimeStep step,
            Population pop, MatrixND matrixCatchPerStrategyMetPerZonePop)
            throws TopiaException, IsisFishException {
        List<PopulationGroup> groups = pop.getPopulationGroup();

        MatrixND result = matrixCatchPerStrategyMetPerZonePop.copy();
        result.setName(ResultName.MATRIX_CATCH_WEIGHT_PER_STRATEGY_MET_PER_ZONE_POP);

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

        return result;
    }

    /**
     * Utilise pour le calcul en Zone
     * 
     * 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 step
     * @param matrixCatchRatePerStrategyMet
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    public MatrixND matrixCatchPerStrategyMetPerZone(MatrixND N,
            Population pop, TimeStep step, 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, step);
        MatrixND result = matrixCatchRatePerStrategyMet.copy();
        result.setName(ResultName.MATRIX_CATCH_PER_STRATEGY_MET_PER_ZONE_POP);

        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;
    }

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

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

        // Optimisation Hilaire
        for (int s = 0; s < strategies.size(); s++) {
            Strategy str = strategies.get(s);
            metiers = getMetiers(str, step);
            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, step, zone);
                    if (effort > 0) {
                        for (int g = 0; g < groups.size(); g++) {
                            PopulationGroup group = groups.get(g);
                            double value = catchRatePerStrategyMet(str, metier,
                                    step, group, zone, matrixFishingMortality);
                            result.setValue(str, metier, group, zone, value);
                        }
                    }
                }
            }
        }
        return result;
    }

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

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

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

        double result = effortPerStrategyPerCell * inter;
        return result;
    }

    // Optimisation Hilaire
    private double catchRatePerStrategyMet(Strategy str, Metier metier,
            TimeStep step, PopulationGroup group, Zone zone,
            MatrixND matrixFishingMortality) throws TopiaException,
            IsisFishException {
        double totalFishingMortality = totalFishingMortality(step,
                matrixFishingMortality).getValue(group, zone);

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

        double fishingMortality = matrixFishingMortality.getValue(str, metier,
                group, zone);
        double totalCatchRate = totalCatchRate(step, group, zone,
                totalFishingMortality);

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

        return result;
    }

    /**
     * @param step
     * @param group
     * @param zone
     * @param totalFishingMortality
     * @return
     * @throws TopiaException
     */
    private double totalCatchRate(TimeStep step, 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 mortalite 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;
    }

    /**
     * Returne une matrice de mortalite group x zone, donc somme sur str et
     * metier
     * 
     * @param step
     * @param matrixFishingMortality
     * @param group
     * @param zone
     * @return
     */
    private MatrixND totalFishingMortality(TimeStep step,
            MatrixND matrixFishingMortality) {
        MatrixND result = matrixFishingMortality.sumOverDim(0);
        result = result.sumOverDim(1);
        result = result.reduceDims(0, 1);
        return result;
    }

    /**
     * Matrice fishing mortality dim [ Strategy x Metier x Classe x zonePop ]
     * 
     * @param pop
     * @param step
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    public MatrixND matrixFishingMortality(TimeStep step, Population pop)
            throws TopiaException, IsisFishException {
        List<Strategy> strategies = getStrategies(step);
        List<Metier> metiers = getMetiers(step);
        List<PopulationGroup> groups = pop.getPopulationGroup();
        List<Zone> zones = pop.getPopulationZone();

        // default value in matrix is 0
        MatrixND result = MatrixFactory.getInstance().create(
                ResultName.MATRIX_FISHING_MORTALITY,
                new List[] { strategies, metiers, groups, zones },
                new String[] { n_("Strategies"), n_("Metiers"), n_("Groups"),
                        n_("Zones") });

        Month month = step.getMonth();
        PopulationSeasonInfo infoPop = pop.getPopulationSeasonInfo(month);

        for (int g = 0; g < groups.size(); g++) {
            PopulationGroup group = groups.get(g);

            // getCapturability is check matrix validity and create matrix if needed and get one value in
            double capturability = infoPop.getCapturability(group);
            if (capturability != 0) { // check 0, this prevent next call, for default value

                for (int m = 0; m < metiers.size(); m++) {
                    Metier metier = metiers.get(m);

                    MetierSeasonInfo infoMet = metier
                            .getMetierSeasonInfo(month);
                    // getTargetFactor seem to be simple
                    double ciblage = infoMet.getTargetFactor(group);

                    if (ciblage != 0) { // check 0, this prevent next call, for default value

                        Gear gear = metier.getGear();
                        Selectivity selectivity = gear
                                .getPopulationSelectivity(pop);
                        if (selectivity != null) {

                            // getCoefficient is equation evaluation
                            double coeff = selectivity.getCoefficient(pop,
                                    group, metier);
                            if (coeff != 0) { // check 0, this prevent next call, for default value

                                for (int s = 0; s < strategies.size(); s++) {
                                    Strategy str = strategies.get(s);
                                    for (int z = 0; z < zones.size(); z++) {
                                        Zone zone = zones.get(z);
                                        double effort = effortPerZonePop(str,
                                                metier, step, zone);
                                        if (effort > 0) { // put value only if <> 0
                                            double value = coeff
                                                    * capturability * ciblage
                                                    * effort;
                                            result.setValue(str, metier, group,
                                                    zone, value);
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        return result;
    }

    /**
     * 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 step
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    private MatrixND matrixCatchPerStrategyMetPerCell(MatrixND N,
            Population pop, TimeStep step,
            MatrixND matrixCatchRatePerStrategyMetPerCell)
            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, step);
        MatrixND result = matrixCatchRatePerStrategyMetPerCell.copy();
        result.setName("matrixCatchPerStrategyMetPerCell");

        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) / (double) zone.sizeCell();
                subsub.mults(val);
            }
        }
        return result;
    }

    /**
     * Matrice des captures en poids dim [ Strategy x Metier x Classe x zonePop ]
     * 
     * @param pop
     * @param step
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    private MatrixND matrixCatchRatePerStrategyMetPerCell(Population pop,
            TimeStep step, MatrixND matrixFishingMortalityPerCell)
            throws TopiaException, IsisFishException {
        List<Strategy> strategies = getStrategies(step);
        List<Metier> metiers = getMetiers(step);
        List<PopulationGroup> groups = pop.getPopulationGroup();
        List<Zone> zones = pop.getPopulationZone();
        List<Cell> cells = getCells(zones);

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

        // Optimisation Hilaire
        MatrixND matrixFishingMortalityPerCellSumOverGroup = matrixFishingMortalityPerCell
                .sumOverDim(2);
        matrixFishingMortalityPerCellSumOverGroup = matrixFishingMortalityPerCellSumOverGroup
                .reduceDims(2);

        for (int s = 0; s < strategies.size(); s++) {
            Strategy str = strategies.get(s);
            metiers = getMetiers(str, step);
            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);
                    List<Cell> cellZones = zone.getCell();
                    for (int c = 0; c < cellZones.size(); c++) {
                        Cell cell = cellZones.get(c);
                        double effort = matrixFishingMortalityPerCellSumOverGroup
                                .getValue(str, metier, zone, cell);
                        if (effort > 0) {
                            for (int g = 0; g < groups.size(); g++) {
                                PopulationGroup group = groups.get(g);
                                double value = catchRatePerStrategyMetPerCell(
                                        str, metier, step, group, zone, cell,
                                        matrixFishingMortalityPerCell);
                                result.setValue(new Object[] { str, metier,
                                        group, zone, cell }, value);
                            }
                        }
                    }
                }
            }
        }

        return result;
    }

    /**
     * catchRatePerStrategyMetPerCell.
     * 
     * Optimisation Hilaire
     * 
     * @param str
     * @param metier
     * @param step
     * @param group
     * @param zone
     * @param cell
     * @param matrixFishingMortalityPerCell
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    private double catchRatePerStrategyMetPerCell(Strategy str, Metier metier,
            TimeStep step, PopulationGroup group, Zone zone, Cell cell,
            MatrixND matrixFishingMortalityPerCell) throws TopiaException,
            IsisFishException {

        MatrixND totalFishingMortalityPerCell = totalFishingMortalityPerCell(
                step, matrixFishingMortalityPerCell);
        double totalFishingMortality = totalFishingMortalityPerCell.getValue(
                group, zone, cell);

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

        double fishingMortalityPerCell = matrixFishingMortalityPerCell
                .getValue(new Object[] { str, metier, group, zone, cell });
        double totalCatchRatePerCell = totalCatchRatePerCell(step, group, zone,
                totalFishingMortality);

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

        return result;
    }

    /**
     * @param step
     * @param group
     * @param zone
     * @param totalFishingMortality
     * @return
     * @throws TopiaException
     */
    private double totalCatchRatePerCell(TimeStep step, PopulationGroup group,
            Zone zone, double totalFishingMortalityPerCell)
            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 mortalite naturelle pour: " + group);
            }
        }
        double F = totalFishingMortalityPerCell;

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

        return result;
    }

    /**
     * Returne une matrice de mortalite group x zone, donc somme sur str et
     * metier
     * 
     * @param step
     * @param matrixFishingMortality
     * @param group
     * @param zone
     * @return
     */
    private MatrixND totalFishingMortalityPerCell(TimeStep step,
            MatrixND matrixFishingMortalityPerCell) {
        MatrixND result = matrixFishingMortalityPerCell.sumOverDim(0);
        result = result.sumOverDim(1);
        result = result.reduceDims(0, 1);
        return result;
    }

    /**
     * Matrice fishing mortality dim [ Strategy x Metier x Classe x zonePop x
     * cellPop]
     * 
     * Nouvelle implantation suite a demande Steph et Sigrid: 20080208 (NouvellesEquationsfev2008ParCelluleAvecAlgo.odt)
     * 
     * @param pop
     * @param step
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    public MatrixND matrixFishingMortalityPerCell(TimeStep step, Population pop)
            throws TopiaException, IsisFishException {
        List<Strategy> strategies = getStrategies(step);
        List<Metier> metiers = getMetiers(step);
        List<PopulationGroup> groups = pop.getPopulationGroup();
        List<Zone> zones = pop.getPopulationZone();
        List<Cell> cells = getCells(zones);

        // default value in matrix is 0
        MatrixND result = MatrixFactory.getInstance().create(
                "matrixFishingMortalityPerCell",
                new List[] { strategies, metiers, groups, zones, cells },
                new String[] { n_("Strategies"), n_("Metiers"), n_("Groups"),
                        n_("Zones"), n_("Cells") });

        Month month = step.getMonth();
        PopulationSeasonInfo infoPop = pop.getPopulationSeasonInfo(month);

        for (int g = 0; g < groups.size(); g++) {
            PopulationGroup group = groups.get(g);

            // getCapturability is check matrix validity and create matrix
            // if needed and get one value in
            double capturability = infoPop.getCapturability(group);
            if (capturability != 0) { // check 0, this prevent next call, for default value
                for (int m = 0; m < metiers.size(); m++) {
                    Metier metier = metiers.get(m);

                    MetierSeasonInfo infoMet = metier
                            .getMetierSeasonInfo(month);
                    // getTargetFactor seem to be simple
                    double ciblage = infoMet.getTargetFactor(group);

                    if (ciblage != 0) { // check 0, this prevent next call, for default value

                        Gear gear = metier.getGear();
                        Selectivity selectivity = gear.getPopulationSelectivity(pop);
                        if (selectivity != null) {

                            // getCoefficient is equation evaluation
                            double coeff = selectivity.getCoefficient(pop,group, metier);
                            if (coeff != 0) { // check 0, this prevent next call, for default value

                                for (int s = 0; s < strategies.size(); s++) {
                                    Strategy str = strategies.get(s);
                                    // l'effort d'une cellule du metier
                                    double effort = effortPerStrategyPerCell(
                                            str, metier, step);
                                    if (effort > 0) { // put value only if <> 0
                                        for (int z = 0; z < zones.size(); z++) {
                                            Zone zone = zones.get(z);
                                            Set<Cell> cellPops = new HashSet<Cell>(zone.getCell());
                                            for (Cell cellMet : infoMet.getCells()) {
                                                if (cellPops.contains(cellMet)) {
                                                    double value = coeff * capturability * ciblage * effort;
                                                    result.setValue(
                                                            new Object[] { str,
                                                                    metier,
                                                                    group,
                                                                    zone,
                                                                    cellMet },
                                                            value);
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        return result;
    }

    /**
     * @param str
     * @param metier
     * @param step
     * @return
     */
    private double effortPerStrategyPerCell(Strategy str, Metier metier,
            TimeStep step) {
        Month month = step.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 pour le metier " + metier
                        + " pour le mois " + month
                        + " avec une zone sans maille: " + zones);
            return 0;
        }

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

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

        double result = effortPerStrategy / nbCell;

        return result;
    }

    /**
     * @param str
     * @param metier
     * @param step
     * @return
     */
    private double effortNominalPerStrategyMet(Strategy str, Metier metier, TimeStep step) {
        Month month = step.getMonth();
        StrategyMonthInfo smi = str.getStrategyMonthInfo(month);
    
        double propSetOfVessels = str.getProportionSetOfVessels();
        int nbOfVessels = str.getSetOfVessels().getNumberOfVessels();
        double propStrMet = smi.getProportionMetier(metier);
        double effortNominalPerVessel = effortNominalPerStrategyPerVessel(str, metier, step);

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

        double result = propSetOfVessels * nbOfVessels * propStrMet * effortNominalPerVessel;

        return result;
    }
    
    /**
     * @param str
     * @param metier
     * @param step
     * @return
     */
    private double effortNominalPerStrategyPerVessel(Strategy str, Metier metier, TimeStep step) {
        Month month = step.getMonth();
        StrategyMonthInfo smi = str.getStrategyMonthInfo(month);
        double nbTrips = smi.getNumberOfTrips();
        double fishingTime = fishingTimePerTrip(str, metier, step);
   

        if(log.isDebugEnabled()) {
            log.debug(
                    " strategy=" + str +
                    " metier=" + metier +
                    " nbTrips=" + nbTrips +
                    " fishingTime=" + fishingTime 
                 
            );
        }
        // nominal timeAtSea = nbTrips * fishingTime;
        //
        double result = nbTrips * fishingTime;

        return result;
    }
    
    /**
     * @param str
     * @param metier
     * @param step
     * @return
     */
    private double effortPerStrategyMet(Strategy str, Metier metier, TimeStep step) {
        Month month = step.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, step);

        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 step
     * @return
     */
    private double effortPerStrategyPerVessel(Strategy str, Metier metier,
            TimeStep step) {
        Month month = step.getMonth();
        StrategyMonthInfo smi = str.getStrategyMonthInfo(month);
        double nbTrips = smi.getNumberOfTrips();
        double fishingTime = fishingTimePerTrip(str, metier, step);
        double stdEffortPerHour = stdEffortPerHour(step, 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 step
     * @return
     */
    protected double fishingTimePerTrip(Strategy str, Metier metier, TimeStep step) {
        Month month = step.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(TimeStep step, SetOfVessels sov, Metier metier) {
        double result = 0;
        EffortDescription ed = sov.getPossibleMetiers(metier);
        if (ed != null) {
            double fstd = metier.getGear().getStandardisationFactor();
            double etp = sov.getTechnicalEfficiency(metier);
            double val = fstd * etp * 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 step
     * @return
     * @throws IsisFishException
     * @throws TopiaException
     */
    private MatrixND matrixAbundancePerCell(MatrixND N, Population pop,
            TimeStep step, MatrixND matrixFishingMortality) throws TopiaException,
            IsisFishException {

        List<PopulationGroup> groups = pop.getPopulationGroup();
        List<Zone> zones = pop.getPopulationZone();
        List<Cell> allCells = getCells(zones);

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

        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);
                List<Cell> cells = zone.getCell();
                for (int c = 0; c < cells.size(); c++) {
                    Cell cell = cells.get(c);
                    double value = survivalRatePerCell(step, group, zone, cell,
                            matrixFishingMortality);
                    double n = N.getValue(g, z) / zone.sizeCell();
                    value *= n;
                    result.setValue(g, z, c, value);
                }
            }
        }

        return result;
    }

    /**
     * Utilise pour le calcule en Cell.
     * 
     * @param N
     * @param pop
     * @param step
     * @return
     * @throws IsisFishException
     * @throws TopiaException
     */
    public MatrixND matrixAbundance(MatrixND N, Population pop, TimeStep step)
            throws TopiaException, IsisFishException {

        MatrixND matrixFishingMortalityPerCell = matrixFishingMortalityPerCell(
                step, pop);

        MatrixND result = matrixAbundancePerCell(N, pop, step,
                matrixFishingMortalityPerCell);
        result = result.sumOverDim(2);
        result = result.reduceDims(2);
        result.setName(ResultName.MATRIX_ABUNDANCE);

        return result;
    }

    /**
     * Utilise pour le calcule en Zone.
     * 
     * @param N
     * @param pop
     * @param step
     * @param matrixFishingMortality
     * @return
     * @throws IsisFishException
     * @throws TopiaException
     */
    public MatrixND matrixAbundance(MatrixND N, Population pop, TimeStep step,
            MatrixND matrixFishingMortality) 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 = survivalRatePerZone(step, group, zone,
                        matrixFishingMortality);
                double n = N.getValue(g, z);
                value *= n;
                result.setValue(g, z, value);
            }
        }

        return result;
    }

    /**
     * Utilise pour la mortalite des poissons, lorsqu'aucune strategie n'est
     * selectionnee.
     * 
     * @param N
     * @param pop
     * @param step
     * @return
     * @throws IsisFishException
     * @throws TopiaException
     */
    public MatrixND matrixAbundanceSsF(MatrixND N, Population pop, TimeStep step)
            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 M = group.getNaturalDeathRate(zone)
                        / (double) Month.NUMBER_OF_MONTH;
                double value = (double) Math.exp(-M);
                
                double n = N.getValue(g, z);
                value *= n;
                result.setValue(g, z, value);
            }
        }

        return result;
    }

    /**
     * @param step
     * @param group
     * @param zone
     * @return
     * @throws IsisFishException
     * @throws TopiaException
     */
    private double survivalRatePerZone(TimeStep step, PopulationGroup group,
            Zone zone, MatrixND matrixFishingMortality) throws TopiaException,
            IsisFishException {
        double F = totalFishingMortality(step, matrixFishingMortality)
                .getValue(group, zone); //totalFishingMortality(step, group, zone);  // rem perf: totalFishingMortality a deja ete calcule
				double M = group.getNaturalDeathRate(zone)
                / (double) Month.NUMBER_OF_MONTH;
        double result = Math.exp(-(F + M));

        return result;
    }

    /**
     * @param step
     * @param group
     * @param zone
     * @return
     * @throws IsisFishException
     * @throws TopiaException
     */
    private double survivalRatePerCell(TimeStep step, PopulationGroup group,
            Zone zone, Cell cell, MatrixND matrixFishingMortalityPerCell)
            throws TopiaException, IsisFishException {

        double F = totalFishingMortalityPerCell(step,
                matrixFishingMortalityPerCell).getValue(group, zone, cell); //totalFishingMortality(step, group, zone);  // rem perf: totalFishingMortality a deja ete calcule
        double M = group.getNaturalDeathRate(zone)
        / (double) Month.NUMBER_OF_MONTH;
        double result = (double) Math.exp(-(F + M));

        return result;
    }

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

    /**
     * @param N
     * @param pop
     * @param step
     * @return
     */
    public MatrixND matrixBiomass(MatrixND N, Population pop, TimeStep step) {
        List<PopulationGroup> groups = (List<PopulationGroup>)N.getSemantic(0);
        List<Zone> zones = (List<Zone>)N.getSemantic(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);
            }
        }

        return result;
    }

    public MatrixND matrixBiomassBeginMonth(MatrixND N, Population pop,
            TimeStep step) {
        List<PopulationGroup> groups = (List<PopulationGroup>)N.getSemantic(0);
        List<Zone> zones = (List<Zone>)N.getSemantic(1);

        MatrixND result = MatrixFactory.getInstance().create(
                ResultName.MATRIX_BIOMASS_BEGIN_MONTH, 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);
            }
        }
        return result;
    }

    public MatrixND matrixAbondanceBeginMonth(MatrixND N, Population pop,
            TimeStep step) {
        List<PopulationGroup> groups = (List<PopulationGroup>)N.getSemantic(0);
        List<Zone> zones = (List<Zone>)N.getSemantic(1);

        MatrixND result = MatrixFactory.getInstance().create(
                ResultName.MATRIX_ABUNDANCE_BEGIN_MONTH, 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 = N.getValue(group, zone);
                result.setValue(group, zone, value);
            }
        }
        return result;
    }

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

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

        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, step);
            for (int m = 0; m < metiers.size(); m++) {
                Metier metier = metiers.get(m);
                double value = effortPerStrategyMet(str, metier, step);
                result.setValue(str, metier, value);
            }
        }

        return result;
    }

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

    /**
     * 
     * 
     * @param pop
     * @param step
     * @param matrixDiscardPerStrategyMetPerZonePop
     * @return
     */
    public MatrixND matrixDiscardWeightPerStrategyMetPerZonePop(Population pop,
            TimeStep step, MatrixND matrixDiscardPerStrategyMetPerZonePop) {
        List<PopulationGroup> groups = pop.getPopulationGroup();

        MatrixND result = matrixDiscardPerStrategyMetPerZonePop.copy();
        result.setName(ResultName.MATRIX_DISCARDS_WEIGHT_PER_STR_MET_PER_ZONE_POP);

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

        return result;
    }

    ///////////////////////////////////////////////////////////////////////////
    //
    //
    //
    /////////////////////////////////////////////////////////////////////////// 
        public MatrixND matrixEffortNominalPerStrategyMet(TimeStep step) throws TopiaException {
      
       List<Strategy> strategies = getStrategies(step);
        List<Metier> metiers = getMetiers(step);

        MatrixND result = MatrixFactory.getInstance().create(
                ResultName.MATRIX_EFFORT_NOMINAL_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, step);
            for (int m=0; m < metiers.size(); m++) {
                Metier metier = metiers.get(m);
                double value = effortNominalPerStrategyMet(str, metier, step);
                result.setValue(str, metier, value);
            }
        }

        return result;
    }


    /**
     * Ce morceau de script sert a calculer la mortalite par peche par espece et par groupe d'age sur l'ensemble de la zone d'etude,
     * en minimisant la difference entre les captures d'ISIS et celles calculees a���������partir de l'equation de Baranov.
     */
	public MatrixND TotalFishingMortality (TimeStep step, ResultStorage resManager)throws TopiaException { // Zone zone, ou,  Population pop, 
	double Fgroup=0;
	double Ftemp=0;
	double Fpop=0;
	//double Ftot=0;
	double Cgroup=0;
	double Ctot=0;
	double SumMtot=0;
	double Mpop=0;
	double Mgroup=0;
	double Mtot=0;
	double NgroupJan=0;
	double NgroupDec=0;
	double DeltaN=0;
	int  dim;
		
	List<Population> populations = getPopulations(step);	
	//List<PopulationGroup> groups = pop.getPopulationGroup();
			
	MatrixND tfmMatrix = MatrixFactory.getInstance().create(
                ResultName.MATRIX_TOTAL_FISHING_MORTALITY, 
				new List[]{populations}, //, step
                new String[]{n_("Populations")}); //, n_("step")
	
		for (Population pop : populations) {
	
			if (step.getMonth() == Month.DECEMBER){
	
				List<PopulationGroup> groups = pop.getPopulationGroup();
				//List<Zone> zone = pop.getPopulationZone();
			
				for (PopulationGroup group : groups) { 
					
                MatrixND catchPerStrategy = resManager.getMatrix(pop, ResultName.MATRIX_CATCH_PER_STRATEGY_MET_PER_ZONE_POP);
				catchPerStrategy = catchPerStrategy.sumOverDim(0,12); // timeStep, avec un pas de 12
				catchPerStrategy = catchPerStrategy.sumOverDim(1); // Strategy
				catchPerStrategy = catchPerStrategy.sumOverDim(2); // Metier
				catchPerStrategy = catchPerStrategy.sumOverDim(4); // Zone : une pop peut avoir plusieurs zonespop dans ISIS
				catchPerStrategy = catchPerStrategy.reduce(); // Enleve les dimensions de taille 1
				//MatrixND catchPerYear = catchPerStrategy.getSubMatrix(1,step.getYear());
				dim = catchPerStrategy.getDimCount();
				log.info("Dim = " + dim);
				Cgroup = catchPerStrategy.getValue(pop.getName(), step.getYear(), group.getAge()); // pop.getName() devrait disparaitre avec le .reduce() d'avant ?
				log.info("Cgroup = " + Cgroup + "Year=" + step.getYear());
				//log.info("total catch per year, per group  = " + catchPerStrategy.getValue(step.getYear(),pop.getPopulationGroup()));
				//log.info("total catch per year, per group  = " + catchPerStrategy.getValue(step.getYear()));
				//Ctot = Ctot+catchPerStrategy.getValue(step.getYear(),pop.getPopulationGroup()); // On va chercher les valeurs du vecteur obtenu (10 valeurs par an) qui correspondent 
				
				MatrixND naturalDeathRatePopGroup = pop.getNaturalDeathRateMatrix();
				Mgroup = naturalDeathRatePopGroup.getValue(group.getAge());
				log.info("Mgroup= " + Mgroup + "Year=" + step.getYear());
				//naturalDeathRatePop = naturalDeathRatePop.meanOverDim(3); // Zone
				// naturalDeathRatePop = naturalDeathRatePop.reduce(); // Enleve les dimensions de taille 1
				//log.info("naturalDeathRatePopGroup= " + naturalDeathRatePopGroup.getValue(step.getYear()) + " step=" + step);
		
				//MatrixND abundancePop = resManager.getMatrix(pop, ResultName.MATRIX_ABUNDANCE);
				//MatrixND abundancePopBeginMonth = resManager.getMatrix(pop, ResultName.MATRIX_ABUNDANCE_BEGIN_MONTH);
				MatrixND abundancePopJan = resManager.getMatrix(new TimeStep(12*step.getYear()), pop, ResultName.MATRIX_ABUNDANCE);
				NgroupJan = abundancePopJan.getValue(group.getAge());
				//abundancePopGroupJan = abundancePopGroupJan.sumOverDim(2); // Zone : (idem)
				log.info("NgroupJan = " + NgroupJan + "Year=" + step.getYear());
				MatrixND abundancePopDec = resManager.getMatrix(step, pop, ResultName.MATRIX_ABUNDANCE);
				NgroupDec = abundancePopDec.getValue(group.getAge());
				//abundancePopGroupDec = abundancePopGroupDec.sumOverDim(2); // Zone : (idem)
				log.info("NgroupDec= " + NgroupDec + "Year=" + step.getYear());
				DeltaN = NgroupJan - NgroupDec ; 
				//DeltaN = abundancePopGroupJan.getValue(step.getYear(),pop.getPopulationGroup())-abundancePopGroupDec.getValue(step.getYear(),pop.getPopulationGroup());
				log.info("DeltaN= " + DeltaN + "Year=" + step.getYear());
				//log.info("DeltaN= " + DeltaN.getValue(step.getYear()));
				//DeltaN = DeltaN+(abundancePopGroupJan.sumAll()-abundancePopDec.sumAll());
				
				//Mtot = SumMtot / getPopulations(step).size();
				Fgroup=fmin(step,0.0,2.0,1.0e-10,Cgroup,Mgroup,DeltaN);
				log.info("Fgroup = " + Fgroup + "Year=" + step.getYear());
				log.info("Cgroup = " + Cgroup + "Year=" + step.getYear());
				log.info("Mgroup= " + Mgroup + "Year=" + step.getYear());
				log.info("DeltaN= " + DeltaN + "Year=" + step.getYear());
				
				Ftemp = Ftemp + Fgroup;
				}
			
			Fpop = Ftemp / groups.size(); // Moyenne sur tous les groupes : corriger pour ne prendre que les groupes representatifs
			log.info("Fpop= " + Fpop + "pop=" + pop +"Year=" + step.getYear());
			}	
			
			else {
			Fpop = 0;			
			}
			
		tfmMatrix.setValue(pop, Fpop); // Bien faire attention �� l'endroit o�� on met cette ��tape (quelle boucle) ?
		return tfmMatrix;	
		}
			
	}
	
	public double deriveF(TimeStep step, double xx, double C, double M, double N) { // Derives Fsq;
		double ff;
		ff = Math.pow((C - (xx/(xx+M))*(1-Math.exp(-(xx+M)))*N),2);
		return ff;  
		}
	
	public double fmin (TimeStep step, double a, double b, double tol, double C, double M, double N) {

		double c,d,e,eps,xm,p,q,r,tol1,t2,u,v,w,fu,fv,fw,fx,x,tol3;

		c = .5*(3.0 - Math.sqrt(5.0));
		d = 0.0;
		eps = 1.2e-16; // 1.1102e-16 is machine precision
		tol1 = eps + 1.0;
		eps = Math.sqrt(eps);
		v = a + c*(b-a);
		//v = 0.2;
		log.info("Finit=" + v);
		w = v;
		x = v;
		e = 0.0;
		fx = 0.0;
		fu = 0.0;
		
		fx = deriveF(step, x, C, M, N);
		
		fv = fx;
		fw = fx;
		tol3 = tol/3.0;
		xm = .5*(a + b);
		tol1 = eps*Math.abs(x) + tol3;
		t2 = 2.0*tol1;

		while (Math.abs(x-xm) > (t2 - .5*(b-a))) { // main loop;
			p = q = r = 0.0;
			if (Math.abs(e) > tol1) { // fit the parabola;
				r = (x-w)*(fx-fv);
				q = (x-v)*(fx-fw);
				p = (x-v)*q - (x-w)*r;
				q = 2.0*(q-r);
				if (q > 0.0) {
					p = -p;
				}
				else {
					q = -q;
				}
				r = e;
				e = d;         
			}
			if ((Math.abs(p) < Math.abs(.5*q*r)) && (p > q*(a-x)) && (p < q*(b-x))) { // parabolic interpolation step;
				d = p/q;
				u = x+d;
				if (((u-a) < t2) || ((b-u) < t2)) { // f must not be evaluated too close to a or b;
					d = tol1;
					if (x >= xm) d = -d;
				}
			}
			else { // a golden-section step;
				if (x < xm) {
					e = b-x;
				}
				else {
					e = a-x;
				}
				d = c*e;
			}
			if (Math.abs(d) >= tol1) { // f must not be evaluated too close to x;
				u = x+d;
			}
			else {
				if (d > 0.0) {
					u = x + tol1;
				}
				else {
					u = x - tol1;
				}
			}
			
			fu = deriveF(step, u, C, M, N);
			

			// Update a, b, v, w, and x
			if (fx <= fu) {
				if (u < x) {
					a = u;
				}
				else {
					b = u;
				}
			}
			if (fu <= fx) {
				if (u < x) {
					b = x;
				}
				else {
					a = x;
				}
				v = w;
				fv = fw;
				w = x;
				fw = fx;
				x = u;
				fx = fu;
				xm = .5*(a + b);
				tol1 = eps*Math.abs(x) + tol3;
				t2 = 2.0*tol1;
			}
			else {
				if ((fu <= fw) || (w == x)) {
					v = w;
					fv = fw;
					w = u;
					fw = fu;
					xm = .5*(a + b);
					tol1 = eps*Math.abs(x) + tol3;
					t2 = 2.0*tol1;
				}
				else if ((fu > fv) && (v != x) && (v != w)) {
					xm = .5*(a + b);
					tol1 = eps*Math.abs(x) + tol3;t2 = 2.0*tol1;
				}
				else {
					v = u;
					fv = fu;
					xm = .5*(a + b);
					tol1 = eps*Math.abs(x) + tol3;
					t2 = 2.0*tol1;
				}
			}
		}
		return x;
	}   
}