/*
 * #%L
 * IsisFish data
 * %%
 * Copyright (C) 2006 - 2014 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 3 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-3.0.html>.
 * #L%
 */

package scripts;

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

import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
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.ResultManager;
import fr.ifremer.isisfish.simulator.SimulationContext;
import fr.ifremer.isisfish.types.Month;
import fr.ifremer.isisfish.types.TimeStep;
import fr.ifremer.isisfish.util.Nocache;

/**
 * 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 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());
        List<Metier> result = info.getMetierWithProportion();

//        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 {

        return matrixToWeightMatrix(step, 2,
                ResultName.MATRIX_CATCH_WEIGHT_PER_STRATEGY_MET_PER_ZONE_MET,
                matrixCatchPerStrategyMetPerZoneMet);

//        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);
                // on iter que sur les metiers qui ont une proportion non 0
                // car la matrice pour les autres metiers, contient 0.
                List<Metier> metiersNot0 = getMetiers(str, step);
                for (int m = 0; m < metiersNot0.size(); m++) {
                    Metier metier = metiersNot0.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 {

        return matrixToWeightMatrix(step, 2,
                ResultName.MATRIX_CATCH_WEIGHT_PER_STRATEGY_MET_PER_ZONE_POP,
                matrixCatchPerStrategyMetPerZonePop);

//        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 {

        int dimGroup = 2;
        int dimZone = 3;

        MatrixND result = matrixCatchRatePerStrategyMet.copy();
        result.setName(ResultName.MATRIX_CATCH_PER_STRATEGY_MET_PER_ZONE_POP);

        for (MatrixIterator i=result.iteratorNotZero(); i.next();) {
            Object[] posSem = i.getSemanticsCoordinates();
            double val = i.getValue();
            double n = N.getValue(posSem[dimGroup], posSem[dimZone]);
            i.setValue(val * n);
        }

        return result;
//        List<PopulationGroup> groups = pop.getPopulationGroup();
//        List<Zone> zones = pop.getPopulationZone();
//
//        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") });

        Month month = step.getMonth();
        // 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);
                MetierSeasonInfo msi = metier.getMetierSeasonInfo(month);
                for (int z = 0; z < zones.size(); z++) {
                    Zone zone = zones.get(z);
                    double effort = effortPerZonePop(str, metier, msi, 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, MetierSeasonInfo infoMet, TimeStep step,
            Zone zonePop) {
        Collection<Zone> zoneMet = infoMet.getZone();
        double inter = nbCellInter(zoneMet, zonePop);

        double effortPerStrategyPerCell = effortPerStrategyPerCell(str, metier,
                infoMet, 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, infoMet, 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 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 (Strategy str : strategies) {
            StrategyMonthInfo smi = str.getStrategyMonthInfo(month);
            double nbTrip = smi.getNumberOfTrips();
            double propSetOfVessels = str.getProportionSetOfVessels();
            int nbOfVessels = str.getSetOfVessels().getNumberOfVessels();

            // dans le calcul de l'effort si une de ces donnees est 0
            // alors le result dans la matrice est 0
            // donc autant ne pas faire tous les autres calculs couteux
            if (nbTrip * propSetOfVessels * nbOfVessels != 0) {
                for (PopulationGroup group : groups) {

                    // 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

                        // et de la meme facon si la proportion est 0 alors
                        // le resultat dans la matrice sera 0, donc on boucle que
                        // sur ceux qui ont une proportion
                        for (Metier metier : getMetiers(str, step)) {

//                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 (Zone zone : zones) {
                                            double effort = effortPerZonePop(str,
                                                    metier, infoMet, 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 {
        int groupDim = 2;
        int zoneDim = 3;

        int NGroupDim = 0;
        int NZoneDim = 1;

        MatrixND result = matrixCatchRatePerStrategyMetPerCell.copy();
        result.setName("matrixCatchPerStrategyMetPerCell");

        List<PopulationGroup> groups = (List<PopulationGroup>)result.getSemantic(groupDim);
        List<Zone> zones = (List<Zone>)result.getSemantic(zoneDim);

        double[] sizeCell = new double[zones.size()];
        int cpt = 0;
        for (Zone z : zones) {
            sizeCell[cpt++] = z.sizeCell();
        }

        // on suppose que les semantics de N et matrixCatchRatePerStrategyMetPerCell
        // pour les PopulationGroup et Zone sont les memes
        // mais on verifie tout de meme
        if (groups.equals(N.getSemantic(NGroupDim)) &&
                zones.equals(N.getSemantic(NZoneDim))) {

            // on place les effectifs dans un tableau pour un acces rapide
            // on divise directement par le nombre de cell de la zone pour eviter
            // de le faire beaucoup plus souvent dans l'autre boucle
            double[][] NArray = new double[groups.size()][zones.size()];
            for (MatrixIterator i=N.iteratorNotZero(); i.next();) {
                int[] pos = i.getCoordinates();
                double val = i.getValue();

                NArray[pos[NGroupDim]][pos[NZoneDim]] = val / sizeCell[pos[NZoneDim]];
            }

            // calcul de la matrice resultat
            for (MatrixIterator i=result.iteratorNotZero(); i.next();) {
                int[] pos = i.getCoordinates();
                double val = i.getValue();
                val *= NArray[pos[groupDim]][pos[zoneDim]];
                i.setValue(val);
            }

        } else {
            // N et matrixCatchRatePerStrategyMetPerCell non pas les meme semantics
            // on fait moins rapide
            for (MatrixIterator i=result.iteratorNotZero(); i.next();) {
                int[] pos = i.getCoordinates();
                Object[] sems = i.getSemanticsCoordinates();
                PopulationGroup group = (PopulationGroup)sems[groupDim];
                Zone zone = (Zone)sems[zoneDim];
                double val = i.getValue();
                val *= N.getValue(group, zone) / sizeCell[pos[zoneDim]];
                i.setValue(val);
            }
        }

        return result;

//        List<PopulationGroup> groups = pop.getPopulationGroup();
//        List<Zone> zones = pop.getPopulationZone();
//
//        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 {

        MatrixND totalFishingMortalityPerCell = totalFishingMortalityPerCell(
                step, matrixFishingMortalityPerCell);

        MatrixND result = matrixFishingMortalityPerCell.copy();
        result.setName("matrixCatchRatePerStrategyMetPerCell");

        for (MatrixIterator i=result.iteratorNotZero(); i.next();) {
            Object[] sems = i.getSemanticsCoordinates();

            PopulationGroup group = (PopulationGroup)sems[2];
            Zone zone = (Zone)sems[3];
            Cell cell = (Cell)sems[4];

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

            if (totalFishingMortality != 0) {

                double fishingMortalityPerCell = i.getValue();
                double totalCatchRatePerCell = totalCatchRatePerCell(step, group, zone,
                        totalFishingMortality);

                double value = fishingMortalityPerCell / totalFishingMortality
                        * totalCatchRatePerCell;

                i.setValue(value);
            }
        }

        return result;


//        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;
    }
// directement dans matrixCatchRatePerStrategyMetPerCell pour pouvoir mieux optimiser
//    /**
//     * 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, infoMet, 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;
//    }
//

    /**
     * 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 (Strategy str : strategies) {
            StrategyMonthInfo smi = str.getStrategyMonthInfo(month);
            double nbTrip = smi.getNumberOfTrips();
            double propSetOfVessels = str.getProportionSetOfVessels();
            int nbOfVessels = str.getSetOfVessels().getNumberOfVessels();

            // dans le calcul de l'effort si une de ces donnees est 0
            // alors le result dans la matrice est 0
            // donc autant ne pas faire tous les autres calculs couteux
            if (nbTrip * propSetOfVessels * nbOfVessels != 0) {
                for (PopulationGroup group : groups) {

                    // 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

                        // et de la meme facon si la proportion est 0 alors
                        // le resultat dans la matrice sera 0, donc on boucle que
                        // sur ceux qui ont une proportion
                        List<Metier> metiersNot0 = getMetiers(str, step);
                        for (Metier metier : metiersNot0) {

                            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

                                        // l'effort d'une cellule du metier
                                        double effort = effortPerStrategyPerCell(
                                                str, metier, infoMet, step);
                                        if (effort > 0) { // put value only if <> 0
                                            for (Zone zone : zones) {
                                                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,
            MetierSeasonInfo infoMet, TimeStep step) {
        //StrategyMonthInfo smi = str.getStrategyMonthInfo(month);
        Collection<Zone> zones = infoMet.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()) {
                Month month = step.getMonth();
                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 sov
     * @param zoneMetier
     * @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) {
        return matrixToWeightMatrix(step, 0, ResultName.MATRIX_BIOMASS, N);
//
//        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) {
        return matrixToWeightMatrix(step, 0, ResultName.MATRIX_BIOMASS_BEGIN_MONTH, N);

//        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) {
        return matrixToWeightMatrix(step, 2,
                ResultName.MATRIX_DISCARDS_WEIGHT_PER_STR_MET_PER_ZONE_POP,
                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.
     * 
     * @param step pas de temps
     * @param pop population
     * @return Fishing mortality per group per year (computed in December)
     */
    public MatrixND fishingMortalityPerGroup(TimeStep step, Population pop, ResultManager resManager) throws TopiaException {   
        double Fgroup = 0;
        double Cgroup = 0;
        double Mgroup = 0;
        double NgroupJan = 0;

        List<Population> populations = Collections.singletonList(pop);
        List<PopulationGroup> groups = pop.getPopulationGroup();

        MatrixND tfgMatrix = MatrixFactory.getInstance().create(
                ResultName.MATRIX_FISHING_MORTALITY_PER_GROUP, 
                new List[]{populations, groups},  // On travaille sur les pops ET les groupes
                new String[]{n("Population"), n("Group")});

        for (PopulationGroup group : groups) {

            if (step.getMonth() == Month.DECEMBER) {

                MatrixND catchPerStrategy = null;

                for (TimeStep loopstep = new TimeStep(step.getYear() * 12); loopstep.beforeOrEquals(step); loopstep=loopstep.next()) { 
                    // On fait cette boucle pour contourner les aspects de cache qui font que les resultats ne sont pas recuperes 
                    //s'ils ont deja ete calcules une fois (meme s'ils ont change depuis)
                    // beforeOrEquals sert a bien prendre Decembre aussi

                    MatrixND catchPerStrategyTemp = resManager.getMatrix(loopstep, pop, ResultName.MATRIX_CATCH_PER_STRATEGY_MET_PER_ZONE_POP);
                    if (catchPerStrategy == null) {
                         catchPerStrategy = catchPerStrategyTemp.clone();
                         // On clone la matrice car si on fait les operations sur celle contenue dans le cache on la modifie et donc on recupere des resultats faux.
                    } else {
                         catchPerStrategy = catchPerStrategy.add(catchPerStrategyTemp); // Pour avoir la somme des captures sur les 12 mois
                    }
                }

                //log.info("catchPerStrategy = " + catchPerStrategy);
                catchPerStrategy = catchPerStrategy.sumOverDim(0); // Strategy
                catchPerStrategy = catchPerStrategy.sumOverDim(1); // Metier
                catchPerStrategy = catchPerStrategy.sumOverDim(3); // Zone : une pop peut avoir plusieurs zonespop dans ISIS
                List semgroup = catchPerStrategy.getSemantic(2);
                catchPerStrategy = catchPerStrategy.reduce(); // Enleve les dimensions de taille 1
                catchPerStrategy.setSemantic(0, semgroup); // Ne plait pas a Eric
                Cgroup = catchPerStrategy.getValue(group);
                //log.info("Cgroup = " + Cgroup + "Year=" + step.getYear());
                //log.info("catchPerStrategy = " + catchPerStrategy + "Year=" + step.getYear());

                MatrixND naturalDeathRatePop = pop.getNaturalDeathRateMatrix();
                naturalDeathRatePop = naturalDeathRatePop.meanOverDim(1); // moyenne sur Zone
                naturalDeathRatePop = naturalDeathRatePop.reduce(); // Enleve les dimensions de taille 1
                Mgroup = naturalDeathRatePop.getValue(group);
                //log.info("Mgroup= " + Mgroup + "Year=" + step.getYear());

                MatrixND abundancePopJan = resManager.getMatrix(new TimeStep(12*step.getYear()), pop, ResultName.MATRIX_ABUNDANCE); // Le timestep 0 correspond a janvier de la premiere annee et les annees sont numerotees a partir de zero
                abundancePopJan = abundancePopJan.sumOverDim(1); // somme sur Zone
                abundancePopJan = abundancePopJan.reduce();
                NgroupJan = abundancePopJan.getValue(group);
                //log.info("NgroupJan = " + NgroupJan + "Year=" + step.getYear());

                ObjectiveFunction f = new ObjectiveFunctionBaranov(Cgroup, Mgroup, NgroupJan);
                Fgroup = MinimisationUtil.fmin(0.0,2.0,1.0e-10, f); // step ??

                //log.info("Fgroup = " + Fgroup);

                tfgMatrix.setValue(pop, group, Fgroup); // Bien faire attention a l'endroit ou on met cette etape (quelle boucle) ?

            } else {
               //Fgroup = 0; // Ce n'est plus une valeur unique mais une matrice, est-ce que cette notation peut fonctionner ?
               tfgMatrix.setValue(pop, group, 0); // Bien faire attention a l'endroit ou on met cette etape (quelle boucle) ?
            }    
        }     

        //log.info("tfgMatrix = " + tfgMatrix);
        //log.info("tfg.DimensionNames = " + Arrays.toString(tfgMatrix.getDimensionNames())); 
        //log.info("tfg.Semantics = " + Arrays.toString(tfgMatrix.getSemantics())); 

        return tfgMatrix;          
    }

    /**
     * Ce morceau de script sert a calculer la mortalite par peche par espece et
     * par groupe representatif a partir de la mortalite par  peche par groupe
     * calculee precedemment.
     */ 
    public MatrixND totalFishingMortality(TimeStep step, Population pop, MatrixND fishingMortalityPerGroup) throws TopiaException {
        MatrixND tfmMatrix = fishingMortalityPerGroup.copy();
        tfmMatrix.setName(ResultName.MATRIX_TOTAL_FISHING_MORTALITY);
        //log.info("tfmMatrix = " + tfmMatrix);

        List<PopulationGroup> groups = pop.getPopulationGroup();

        int groupMin = pop.getGroupMin();
        int groupMax = pop.getGroupMax();
        int Nbre = (int)groupMax - (int)groupMin + 1;
        //log.info("Nbre = " + Nbre);

        for (PopulationGroup group : groups) {
            if (group.getId() == groupMin) {
                // MatrixND getSubMatrix(int dim, Object, int nb)
                tfmMatrix = tfmMatrix.getSubMatrix(1, group, Nbre); 
                tfmMatrix = tfmMatrix.meanOverDim(1);
                tfmMatrix = tfmMatrix.reduce();
                //log.info("tfm.DimensionNames = " + Arrays.toString(tfmMatrix.getDimensionNames())); 
                //log.info("tfm.Semantics = " + Arrays.toString(tfmMatrix.getSemantics())); 
            }
        }

        return tfmMatrix;
    }

    /**
     * Converti une matrix contenant des effectifs en une matrice en poids
     *
     * @param step
     * @param dimGroup
     * @param resultName
     * @return 
     */
    @Nocache // no need to cache it, because caller method is cached
    public MatrixND matrixToWeightMatrix(TimeStep step, int dimGroup, String resultName,
            MatrixND matrix) {

        MatrixND result = matrix.copy();
        result.setName(resultName);

        List<PopulationGroup> groups = (List<PopulationGroup>)result.getSemantic(dimGroup);
        double[] meanWeight = new double[groups.size()];
        int cpt = 0;
        for (PopulationGroup group : groups) {
            meanWeight[cpt++] = group.getMeanWeight();
        }

        for (MatrixIterator i=result.iteratorNotZero(); i.next();) {
            double val = i.getValue();
            int[] pos = i.getCoordinates();
            i.setValue(val * meanWeight[pos[dimGroup]]);
        }

        return result;
    }
}
