/*
 * #%L
 * IsisFish data
 * %%
 * Copyright (C) 2006 - 2015 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 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.*;
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;
import resultinfos.*;

/**
 * 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;
    protected ResultManager resManager;

    /**
     * Method used to get SiMatrix used for simulation, before call this method
     * new SiMatrix must called
     * 
     * @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());
        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();
        this.resManager = context.getResultManager();
        setSiMatrix(context, this);
    }

    /**
     * This method permit to add some specific computation for specific Simulator.
     * example use in SiMatrixEffortByCell
     *
     * @param step
     * @param pop
     * @param N
     * @throws IsisFishException 
     */
    public void computeMonthExtra(TimeStep step, Population pop, MatrixND N) throws TopiaException, IsisFishException {
    }

    /**
     * @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(
                MatrixMetierZone.NAME, 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(
                MatrixPrice.NAME, new List[] { groups },
                new String[] { n("PopulationGroup") });

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

        return result;
    }

    /**
     * Since 4.4, FishingMortalityOtherFleets must be added everewhere NaturalDeathRate was
     * used.
     * 
     * @param step current step
     * @param group group
     * @param zone zone
     * @return
     */
    public double getTotalDeathRate(TimeStep step, PopulationGroup group, Zone zone) {
        double n = group.getNaturalDeathRate(zone);
        double o = group.getFishingMortalityOtherFleets(zone);
        return n + o;
    }

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

    public MatrixND matrixCatchWeightPerStrategyMetPerZonePop(TimeStep step,
            Population pop, MatrixND matrixCatchPerStrategyMetPerZonePop)
            throws TopiaException, IsisFishException {

        return matrixToWeightMatrix(step, 2,
                MatrixCatchWeightPerStrategyMetPerZonePop.NAME,
                matrixCatchPerStrategyMetPerZonePop);
    }

    /**
     * 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
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    public MatrixND matrixCatchPerStrategyMetPerZone(MatrixND N,
            Population pop, TimeStep step)
            throws TopiaException, IsisFishException {

        MatrixND matrixCatchRatePerStrategyMet = matrixCatchRatePerStrategyMetPerZone(pop, step);
        
        int dimGroup = 2;
        int dimZone = 3;

        MatrixND result = matrixCatchRatePerStrategyMet.copy();
        result.setName(MatrixCatchPerStrategyMetPerZonePop.NAME);

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

    /**
     * Utilise pour le calcul en Zone Matrice des captures en poids dim [
     * Strategy x Metier x Classe x zonePop ]
     * 
     * @param pop
     * @param step
     * @return
     * @throws TopiaException
     * @throws IsisFishException
     */
    public MatrixND matrixCatchRatePerStrategyMetPerZone(Population pop, TimeStep step) throws TopiaException,
            IsisFishException {
        
        MatrixND matrixFishingMortality = matrixFishingMortality(step, pop);
        
        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(
                MatrixCatchRatePerStrategyMet.NAME,
                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);
                        }
                    }
                }
            }
        }
        resManager.addResult(step, pop, result);
        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;
    }

    /**
     * @param str
     * @param metier
     * @param infoMet 
     * @param step
     * @return
     */
    protected 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;
    }

    // Optimisation Hilaire
    protected 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 = getTotalDeathRate(step, group, 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(
                MatrixFishingMortality.NAME,
                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);
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        resManager.addResult(step, pop, result);
        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
     */
    protected 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();
    }

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

    /**
     * Utilise pour le calcule en Zone.
     * 
     * @param N
     * @param pop
     * @param step
     * @return
     * @throws IsisFishException
     * @throws TopiaException
     */
    public MatrixND matrixAbundance(MatrixND N, Population pop, TimeStep step) throws TopiaException,
            IsisFishException {
        
        MatrixND matrixFishingMortality = matrixFishingMortality(step, pop);
        
        List<PopulationGroup> groups = pop.getPopulationGroup();
        List<Zone> zones = pop.getPopulationZone();

        MatrixND result = MatrixFactory.getInstance().create(
                MatrixAbundance.NAME, 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(
                MatrixAbundance.NAME, 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 = getTotalDeathRate(step, group, 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
     * @param matrixFishingMortality 
     * @return
     * @throws IsisFishException
     * @throws TopiaException
     */
    protected 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 = getTotalDeathRate(step, group, zone)
                / (double) Month.NUMBER_OF_MONTH;
        double result = 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, MatrixBiomass.NAME, 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, MatrixBiomassBeginMonth.NAME, 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(
                MatrixAbundanceBeginMonth.NAME, 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(
                MatrixEffortPerStrategyMet.NAME,
                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,
                MatrixDiscardsWeightPerStrMetPerZonePop.NAME,
                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(
                MatrixEffortNominalPerStrategyMet.NAME,
                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 NgroupRef = 0;

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

        MatrixND tfgMatrix = MatrixFactory.getInstance().create(
                MatrixFishingMortalityPerGroup.NAME, 
                new List[]{groups},  // On travaille sur les groupes
                new String[]{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;
                    if (pop.isComputeFOnLandings()) {
                        catchPerStrategyTemp = resManager.getMatrix(loopstep, pop, MatrixLandingPerMet.NAME);
                    } else {
                        catchPerStrategyTemp = resManager.getMatrix(loopstep, pop, MatrixCatchPerStrategyMetPerZonePop.NAME);
                    }
                    if (catchPerStrategy == null) {
                        catchPerStrategy = catchPerStrategyTemp.copy();
                        // 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
                catchPerStrategy = catchPerStrategy.reduceDims(0, 1, 3); // Enleve les dimensions de taille 1
                Cgroup = catchPerStrategy.getValue(group);

                // compute total mortality by suming naturalDeathRate and fishingMortalityOtherFleetsMatrix
                MatrixND naturalDeathRatePop = pop.getNaturalDeathRateMatrix();
                naturalDeathRatePop = naturalDeathRatePop.meanOverDim(1); // moyenne sur Zone
                naturalDeathRatePop = naturalDeathRatePop.reduce(); // Enleve les dimensions de taille 1
                MatrixND fishingMortalityOther = pop.getFishingMortalityOtherFleetsMatrix();
                fishingMortalityOther = fishingMortalityOther.meanOverDim(1); // moyenne sur Zone
                fishingMortalityOther = fishingMortalityOther.reduce(); // Enleve les dimensions de taille 1
                Mgroup = naturalDeathRatePop.getValue(group) + fishingMortalityOther.getValue(group);
                //log.info("Mgroup= " + Mgroup + "Year=" + step.getYear());

                Month abundanceReferenceMonth = pop.getAbundanceReferenceMonth();
                if (abundanceReferenceMonth == null) {
                    abundanceReferenceMonth = Month.JANUARY;
                }
                TimeStep refStep = new TimeStep(12 * step.getYear() + abundanceReferenceMonth.getMonthNumber());
                // Le timestep 0 correspond a janvier de la premiere annee et les annees sont numerotees a partir de zero
                MatrixND abundancePopRef = resManager.getMatrix(refStep, pop, MatrixAbundanceBeginMonth.NAME);
                abundancePopRef = abundancePopRef.sumOverDim(1); // somme sur Zone
                abundancePopRef = abundancePopRef.reduce();
                NgroupRef = abundancePopRef.getValue(group);
                //log.info("NgroupRef = " + NgroupJan + "Year=" + step.getYear());

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

                tfgMatrix.setValue(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(group, 0); // Bien faire attention a l'endroit ou on met cette etape (quelle boucle) ?
            }    
        }     

        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(MatrixTotalFishingMortality.NAME);

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

        int groupMin = pop.getFbarGroupMin();
        int groupMax = pop.getFbarGroupMax();
        int groupCount = (int)groupMax - (int)groupMin + 1;

        for (PopulationGroup group : groups) {
            if (group.getId() == groupMin) {
                // MatrixND getSubMatrix(int dim, Object, int nb)
                tfmMatrix = tfmMatrix.getSubMatrix(0, group, groupCount); 
                tfmMatrix = tfmMatrix.meanOverDim(0);
                tfmMatrix = tfmMatrix.reduce();
            }
        }

        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;
    }
    
    /**
     * Permet de calculer la longueur moyenne (Lbar) d'une population.
     * 
     * @param step
     * @param pop
     * @param N
     * @return 
     */
    public MatrixND Lbar(TimeStep step, Population pop, MatrixND N) throws TopiaException {
        MatrixND LbarMatrix = N.copy();
        LbarMatrix.setName(MatrixLbar.NAME);
        //log.info("Lbar_1 = " + LbarMatrix);

        List<PopulationGroup> groups = pop.getPopulationGroup();
        List<Zone> zones = pop.getPopulationZone();
        
        // On multiplie l'abondance de chaque groupe par la taille de ce groupe
        
        // On somme sur les groupes
        
        // On utilise divise la valeur obtenue par l'abondance totale

        for (PopulationGroup group : groups) {
            for (Zone zn : zones) {
                LbarMatrix.setValue(group, zn, LbarMatrix.getValue(group, zn) * group.getLength());
            }
        }
        
        LbarMatrix = LbarMatrix.sumOverDim(0); // Somme sur les groupes
        LbarMatrix = LbarMatrix.reduce();
        //log.info("LbarMatrix_2 = " + LbarMatrix);
        
        MatrixND matAB = N.copy();
        matAB = matAB.sumOverDim(0);
        matAB = matAB.reduce();
        //log.info("matAB = " + matAB);
        
        for (Zone zn : zones) {
            LbarMatrix.setValue(zn, LbarMatrix.getValue(zn) / matAB.getValue(zn));
        }

        //log.info("LbarMatrix_3 = " + LbarMatrix);

        return LbarMatrix;
    }

    /**
     * Permet de calculer la longueur moyenne total (LbarTotal) d'une population.
     * 
     * @param step
     * @param pop
     * @param N
     * @return 
     */
    public MatrixND LbarTotal(TimeStep step, Population pop, MatrixND N) throws TopiaException {
        MatrixND LbarTotalMatrix = N.copy();
        LbarTotalMatrix = LbarTotalMatrix.sumOverDim(1).reduceDims(1); // reduction des zones
        double ab = LbarTotalMatrix.sumAll();

        LbarTotalMatrix.setName(MatrixLbarTotal.NAME);
        List<PopulationGroup> groups = pop.getPopulationGroup();

        for (PopulationGroup group : groups) {
            LbarTotalMatrix.setValue(group, LbarTotalMatrix.getValue(group) * group.getLength());
        }

        LbarTotalMatrix = LbarTotalMatrix.sumOverDim(0); // Somme sur les groupes
        LbarTotalMatrix = LbarTotalMatrix.reduce();

        LbarTotalMatrix = LbarTotalMatrix.divs(ab);

        return LbarTotalMatrix;
    }

    /**
     * Permet de sortir le recrutement par pas de temps.
     * @param step
     * @param pop
     * @return 
     */
    public MatrixND Recruitment(TimeStep step, Population pop) throws TopiaException {
        
        List<Population> populations = Collections.singletonList(pop);

        MatrixND RecruitmentMatrix = MatrixFactory.getInstance().create(
                MatrixRecruitment.NAME, 
                new List[]{populations},  // 
                new String[]{n("Population")});

        //log.info("Recruitment_1 = " + RecruitmentMatrix);
        
        // Pour chaque pas de temps on recupere le recrutement
        RecruitmentMatrix.setValue(pop, context.getPopulationMonitor().getRecruitment(step, pop).sumAll());

        //log.info("RecruitmentMatrix_2 = " + RecruitmentMatrix);

        return RecruitmentMatrix;
    }
}
