Add new dispersion function, milsol and climate models

src/tests/multithreading/main.cpp

@@ -107,7 +107,7 @@ double grid_multithreading(int cluster_number)
             paradevs::tests::boost_graph::PartitioningParameters >
         > rc(0, 10, "root", paradevs::common::NoParameters(),
              paradevs::tests::boost_graph::PartitioningParameters(
-                 cluster_number, "rand", 20, false, generator));
+                 cluster_number, "gggp", 20, false, generator));
 
     steady_clock::time_point t1 = steady_clock::now();
 
@@ -209,8 +209,8 @@ void tree(int n)
 int main(int argc, char** argv)
 {
     if (argc > 1) {
-        // grid(atoi(argv[1]));
-        tree(atoi(argv[1]));
+        grid(atoi(argv[1]));
+        // tree(atoi(argv[1]));
     }
     return 0;
 }

src/tests/plot/climate.hpp

@@ -0,0 +1,85 @@
+/**
+ * @file tests/plot/climate.hpp
+ * @author The PARADEVS Development Team
+ * See the AUTHORS or Authors.txt file
+ */
+
+/*
+ * PARADEVS - the multimodeling and simulation environment
+ * This file is a part of the PARADEVS environment
+ *
+ * Copyright (C) 2013-2015 ULCO http://www.univ-litoral.fr
+ *
+ * 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/>.
+ */
+
+#ifndef TESTS_PLOT_CLIMATE_HPP
+#define TESTS_PLOT_CLIMATE_HPP 1
+
+namespace paradevs { namespace tests { namespace plot {
+
+class Climate
+{
+public:
+    Climate()
+    {
+    }
+
+    virtual void get_humidities(typename common::DoubleTime::type t,
+                                std::vector < double >& humidities) const
+    {
+        humidities.clear();
+        for (unsigned int i = 0; i < 24; ++i) {
+            humidities.push_back(95.);
+        }
+    }
+
+    virtual void get_temperatures(typename common::DoubleTime::type t,
+                                  std::vector < double >& temperatures) const
+    {
+        temperatures.clear();
+        for (unsigned int i = 0; i < 9; ++i) {
+            temperatures.push_back(10.);
+        }
+        for (unsigned int i = 9; i < 19; ++i) {
+            temperatures.push_back(20.);
+        }
+        for (unsigned int i = 19; i < 24; ++i) {
+            temperatures.push_back(10.);
+        }
+    }
+
+    virtual void get_wind_direction(
+        typename common::DoubleTime::type t,
+        std::vector < double >& wind_direction) const
+    {
+        wind_direction.clear();
+        for (unsigned int i = 0; i < 24; ++i) {
+            wind_direction.push_back(0.);
+        }
+    }
+
+    virtual void get_wind_speed(typename common::DoubleTime::type t,
+                                std::vector < double >& wind_speed) const
+    {
+        wind_speed.clear();
+        for (unsigned int i = 0; i < 24; ++i) {
+            wind_speed.push_back(20.);
+        }
+    }
+};
+
+} } } // namespace paradevs tests plot
+
+#endif

src/tests/plot/cohorte.hpp

@@ -0,0 +1,308 @@
+/**
+ * @file tests/plot/cohorte.hpp
+ * @author The PARADEVS Development Team
+ * See the AUTHORS or Authors.txt file
+ */
+
+/*
+ * PARADEVS - the multimodeling and simulation environment
+ * This file is a part of the PARADEVS environment
+ *
+ * Copyright (C) 2013-2015 ULCO http://www.univ-litoral.fr
+ * Copyright (C) 2009 INRA
+ *
+ * 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/>.
+ */
+
+#ifndef TESTS_PLOT_COHORTE_HPP
+#define TESTS_PLOT_COHORTE_HPP 1
+
+namespace paradevs { namespace tests { namespace plot {
+
+class Cohorte
+{
+public:
+    Cohorte(int id, double InocPrim, double P1, double P2, double P3,
+            double P4, double P5, double P6, double P7, double P8,
+            double P9, double P10, double P11, double P12, double P13,
+            double P14, double d0, double d1, double dc, double cum0,
+            double cum1, double topt, double tmin, double fact, double ssa)
+    {
+        ID = id;
+        count_AFFAIB = 0;
+
+        TUDESPO = 0.0;
+        CUMDDS = 0.0;
+        SURVIE = 1.0;
+        SPORES = InocPrim;
+        GRAVI = 0.0;
+        POIDS = 0.0;
+        RET = 0.0;
+        TINCUB = 0.0;
+        AGE = 0.0;
+        KASPOTHEO = 0.0;
+        TRED = 0.0;
+        AFFAIB = 0.0;
+        KASPOREELLE = 0.0;
+        POSPO = 0.0;
+        SPOSPO = 0.0;
+        TUSPORU = 0.0;
+        CUMSPO = 0.0;
+        ACTISPO = 0.0;
+        SPORUL = 0.0;
+        SURFMIL = 0.0;
+
+        p1 = P1;
+        p2 = P2;
+        p3 = P3;
+        p4 = P4;
+        p5 = P5;
+        p6 = P6;
+        p7 = P7;
+        p8 = P8;
+        p9 = P9;
+        p10 = P10;
+        p11 = P11;
+        p12 = P12;
+        p13 = P13;
+        p14 = P14;
+        D0 = d0;
+        D1 = d1;
+        Dc = dc;
+        CUM0 = cum0;
+        CUM1 = cum1;
+        Topt = topt;
+        Tmin = tmin;
+        FACT = fact;
+        SSA = ssa;
+    }
+
+    virtual ~Cohorte()
+    { };
+
+    void compute(double Tmoy, double HR)
+    {
+        if (CUMDDS < D1) { //Compartiment Survie et Contamination
+            //TUDESPO
+            if (Tmoy < 18.0) {
+                TUDESPO = p1 * Tmoy + p2;
+            } else {
+                TUDESPO = p3;
+            }
+
+            //CUMDDS
+            CUMDDS += TUDESPO;
+
+            //SURVIE
+            if (HR > 90.0) {
+                SURVIE = 1.0;
+            } else if (CUMDDS <= Dc) {
+                SURVIE = 1 - CUMDDS / Dc;
+            } else {
+                SURVIE = 0.0;
+            }
+
+            //SPORES
+            SPORES *= SURVIE;
+
+            //GRAVI
+            if (HR < 90.0) {
+                GRAVI = 0.0;
+            } else {
+                GRAVI = (CUMDDS - D0) / (D1 - D0);
+            }
+
+            //POIDS
+            POIDS = SPORES * GRAVI;
+        } else { //Compartiment Incubation et Sporulation
+            //potentielle & réelle
+            //RET
+            if (Tmoy <= Topt) {
+                RET = 0.0;
+            } else {
+                RET = p5 * std::pow((Tmoy - Topt), p6);
+            }
+
+            //TINCUB
+            if (Tmoy <= 18.0) {
+                TINCUB = p4 * Tmoy;
+            } else {
+                TINCUB = p4 * Tmoy - RET;
+            }
+
+            //AGE
+            AGE += TINCUB;
+
+            if (AGE > p7)
+            {
+                //KASPOTHEO
+                if (AGE < p8) {
+                    KASPOTHEO = FACT * (AGE - p7) / p7;
+                } else if (AGE < p9) {
+                    KASPOTHEO = FACT * (p9 - AGE) / p9;
+                } else {
+                    KASPOTHEO = 0.0;
+                }
+
+                //TRED
+                if (Tmoy <= 18.0) {
+                    TRED = 0.0;
+                } else {
+                    TRED = FACT * p10 * std::pow(Tmoy - Topt, p6);
+                }
+
+                //AFFAIB
+                count_AFFAIB++;
+                if (count_AFFAIB % 12 != 1) {
+                    AFFAIB += TRED;
+                } else {
+                    AFFAIB = TRED; // Remise a 0 de AFFAIB toutes
+                    // les 12 heures
+                }
+                //KASPOREELLE
+                KASPOREELLE = KASPOTHEO - AFFAIB;
+
+                //POSPO
+                POSPO = POIDS * KASPOREELLE;
+
+                //SPOSPO
+                SPOSPO += POSPO;
+
+                //TUSPORU
+                if (HR <= 90.0) {
+                    TUSPORU = 0.0;
+                } else {
+                    TUSPORU = p11 * std::pow(Tmoy - Tmin, 2) *
+                        (p13 - p14 * (Tmoy - Tmin));
+                }
+
+                //CUMSPO
+                CUMSPO += TUSPORU;
+
+                //ACTISPO
+                if (CUMSPO < CUM0) {
+                    ACTISPO = 0.0;
+                } else if (CUMSPO < CUM1) {
+                    ACTISPO = (CUMSPO - CUM0) / (CUM1 - CUM0);
+                } else {
+                    ACTISPO = 1.0;
+                }
+
+                //SPORUL
+                SPORUL = ACTISPO * SPOSPO;
+
+                //SURFMIL
+                SURFMIL = SPORUL * SSA;
+            }
+        }
+    }
+
+    double age() const
+    { return AGE; }
+
+    double sporul() const
+    { return SPORUL; }
+
+private:
+    int ID;
+    unsigned int count_AFFAIB;
+
+    //paramètres du modèle de cohorte
+    double p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14;
+    double D0, D1, Dc;
+    double CUM0, CUM1;
+    double Topt, Tmin;
+    double FACT, SSA;
+
+    //variables d'état du modèle de cohorte
+    double TUDESPO; // Unité de développement des spores en 1
+                    // heure en fonction de la température horaire
+    double CUMDDS;  // Evolution du développement des spores depuis
+                    // l'instant initiale (arrivée de l'inoculum
+                    // primaire). Il faut 100 unités pour que
+                    // commence la germinantion des spores. A 150
+                    // unités, toutes les spores ont germé.
+    double SURVIE;  // Calcul de la survie des spores en fonction
+                    // de l'HR et de CUMDDS. Si HR > 90, toutes
+                    // les spores survivent ; si HR<90, la
+                    // mortalité des spores augmente au fur et à
+                    // mesure que  CUMDDS apporche de 100. Les
+                    // spores qui ont commencé à émettre un tube
+                    // germinantif meurent en 1 heure si HR < 90.
+    double SPORES;  // C'est le nombre de spores au cours du temps
+                    // en fonction de la survie. Cette variable
+                    // donne à la fin (CUMDDS = 150) le nombre de
+                    // spores réellement prêtes pour germer après
+                    // la mortatilté due aus conditions climatiques.
+    double GRAVI;   // Proportion de spores (parmi celles qui ont
+                    // survecu = SPORES) ayant germé en fonction de HR
+    double POIDS;   // Cette variable donne le nombre de spores
+                    // contaminatrices, càd celles qui entrent en
+                    // période d'incubation et qui sont
+                    // susceptibles de sporuler en fonction des
+                    // conditions de développement et climatiques
+    double RET;     // Les températures > 18 °C diminuent la
+                    // valeur de TINCUB. Cette variable permet de
+                    // calculer cette diminution.
+    double TINCUB;  // Unité de développement horaire des
+                    // mycéliums pour provoquer les lésions (en
+                    // fonction de  la témpérature horaire)
+    double AGE;     // C'est la valeur au cours du temps du
+                    // développement des mycéliums. 75 = nombre
+                    // d'unités nécessaire pour la période
+                    // d'incubation, ce qui correspond au début
+                    // sporulation si condition favorable 150 =
+                    // âge maximale de la lésion correspondant au
+                    // maximum de sporulation 225 = toute la
+                    // lésion est nécrosée, ce qui correspond à la
+                    // fin de la sporulation si condition favorable
+    double KASPOTHEO; // C'est la capacité de sporulation
+                      // théorique d'une spore qui commence
+                      // uniquement après qu'on ait atteint AGE =
+                      // 75. Il augmente et atteint son maximum à
+                      // AGE =150 et diminue jusqu'à 0 vers AGE = 225.
+    double TRED;    // Quand T > 18 °C, il y a réduction de KASPO,
+                    // cette variable quantifie cette diminution
+                    // en fonction de la température horaire.
+    double AFFAIB;  // C'est la somme de TRED qui est calculée par
+                    // période de 12 heures.
+    double KASPOREELLE; // C'est la capacité de sporulation réelle
+                        // après enlèvement des effets dus à AFFAIB
+    double POSPO;   // C'est la capacité de sporulation de
+                    // l'ensemble des spores contaminatrices (POIDS)
+    double SPOSPO;  // Potentiel de sporulation des lésions
+                    // (causées par toutes les spores
+                    // contaminatrices, tant que les conditions
+                    // sont favorables) jusqu'à ce qu'elles ne
+                    // sont plus en mesure de sporuler (necrose
+                    // des lésions).
+    double TUSPORU; // Unité de développement de la sporulation en
+                    // fonction de la température horaire
+    double CUMSPO;  // Degré de développement de la sporulation au
+                    // cours du temps.
+    double ACTISPO; // Proportion de spores effectivement
+                    // produites par rapport à sporulation
+                    // potentielle. Il faut CUMSPO = 6 pour
+                    // débuter la sporulation réelle CUMSPO = 10,
+                    // le potentiel est atteint
+    double SPORUL;  // C'est le nombre de spores réellement
+                    // produites et prêtes à être dispersées.
+    double SURFMIL; // Conversion du nombre de spores en une
+                    // taille de lésion. Elle correspond à la
+                    // surface mildiousée.
+};
+
+} } } // namespace paradevs tests plot
+
+#endif

src/tests/plot/dispersion.hpp

@@ -0,0 +1,76 @@
+/**
+ * @file tests/plot/dispersion.hpp
+ * @author The PARADEVS Development Team
+ * See the AUTHORS or Authors.txt file
+ */
+
+/*
+ * PARADEVS - the multimodeling and simulation environment
+ * This file is a part of the PARADEVS environment
+ *
+ * Copyright (C) 2013-2015 ULCO http://www.univ-litoral.fr
+ *
+ * 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/>.
+ */
+
+#ifndef TESTS_PLOT_DISPERSION_HPP
+#define TESTS_PLOT_DISPERSION_HPP 1
+
+namespace paradevs { namespace tests { namespace plot {
+
+#define SCALE 4
+
+class KleinDispersionFunction
+{
+public:
+    KleinDispersionFunction()
+    {
+    }
+
+    virtual double operator()(
+        typename common::DoubleTime::type t,
+        const paradevs::tests::boost_graph::Point& source,
+        const paradevs::tests::boost_graph::Point& destination,
+        double wind_speed,
+        double wind_direction) const
+    {
+        double distance = std::sqrt(
+            (destination._x - source._x) * (destination._x - source._x) +
+            (destination._y - source._y) * (destination._y - source._y));
+        double r, rt;
+
+        r = distance * SCALE / 100;
+        if (r <= 1.5) {
+            rt = 0.340 - 0.405 * r + 0.128 * std::pow (r, 2.0);
+        } else if (r <= 50) {
+            rt = 0.03985 / (1.0 + std::pow (r, 3.12) / 3.80);
+        } else { // r > 50m
+            double rr, gamma, K;
+
+            rr = 50;
+            gamma = -2.29;
+            //valeur la plus lourde -2.14, + legere - 2.56
+            K = (0.03985 / (1.0 + std::pow (rr, 3.12) / 3.80)) / std::pow (1.0 + rr, gamma);
+            rt = K * std::pow (1.0 + r, gamma);
+
+            // std::cout << _index << "\t" << t << "\t" << r << "\t" << K << "\t" << rr << "\t" << rt << "\t" << value << std::endl;
+
+        }
+        return rt;
+    }
+};
+
+} } } // namespace paradevs tests plot
+
+#endif

src/tests/plot/graph_builder.hpp

@@ -110,6 +110,15 @@ private:
 
         SHPGetInfo(hSHP, &nEntities, &nShapeType, adfMinBound, adfMaxBound);
 
+        // for (int i = 0; i < 4; ++i) {
+        //     std::cout << adfMinBound[i] << " ";
+        // }
+        // std::cout << std::endl;
+        // for (int i = 0; i < 4; ++i) {
+        //     std::cout << adfMaxBound[i] << " ";
+        // }
+        // std::cout << std::endl;
+
         PrecisionModel precisionModel(100., 0, 0);
         GeometryFactory geomFactory(&precisionModel);
         const CoordinateSequenceFactory* seqFactory =
@@ -166,21 +175,56 @@ private:
         }
 
         // compute neighbour
-        int buffer = 10;
+        // int buffer = 10;
+        // std::map < int, Polygon* >::const_iterator it = polygons.begin();
+
+        // while (it != polygons.end()) {
+        //     Geometry* bufferedPolygon = it->second->buffer(buffer);
+        //     std::map < int, Polygon* >::const_iterator it2 = polygons.begin();
+
+        //     while (it2 != polygons.end()) {
+        //         if (it->first != it2->first) {
+        //             if (bufferedPolygon->intersects(it2->second)) {
+        //                 boost::add_edge(vertices[it->first],
+        //                                 vertices[it2->first], 1., g);
+
+        //                 g[vertices[it2->first]]._neighbour_centroids.push_back(
+        //                     paradevs::tests::boost_graph::Point(
+        //                         g[vertices[it->first]]._centroid));
+        //             }
+        //         }
+        //         ++it2;
+        //     }
+        //     delete bufferedPolygon;
+        //     ++it;
+        // }
+
+        // compute neighbour
         std::map < int, Polygon* >::const_iterator it = polygons.begin();
 
         while (it != polygons.end()) {
-            Geometry* bufferedPolygon = it->second->buffer(buffer);
             std::map < int, Polygon* >::const_iterator it2 = polygons.begin();
 
             while (it2 != polygons.end()) {
-                if (bufferedPolygon->intersects(it2->second)) {
-                    boost::add_edge(vertices[it->first],
-                                    vertices[it2->first], 1., g);
+                if (it->first != it2->first) {
+                    // geos::algorithm::CentroidArea centroid;
+
+                    // centroid.add(dynamic_cast < Geometry* >(it->second));
+
+                    // std::cout << it->first << " " << it2->first << " " << it->second->distance(it2->second)
+                    //           << " " << centroid.getCentroid()->x << " " << centroid.getCentroid()->y << std::endl;
+
+                    if (it->second->distance(it2->second) < 10) {
+                        boost::add_edge(vertices[it->first],
+                                        vertices[it2->first], 1., g);
+
+                        g[vertices[it2->first]]._neighbour_centroids.push_back(
+                            paradevs::tests::boost_graph::Point(
+                                g[vertices[it->first]]._centroid));
+                    }
                 }
                 ++it2;
             }
-            delete bufferedPolygon;
             ++it;
         }
 

src/tests/plot/graph_defs.hpp

@@ -1,5 +1,5 @@
 /**
- * @file tests/boost_graph/graph_defs.hpp
+ * @file tests/plot/graph_defs.hpp
  * @author The PARADEVS Development Team
  * See the AUTHORS or Authors.txt file
  */
@@ -28,25 +28,44 @@
 #define __TESTS_PLOT_GRAPH_DEFS_HPP 1
 
 #include <boost/graph/adjacency_list.hpp>
+#include <vector>
 
 namespace paradevs { namespace tests { namespace boost_graph {
 
+struct Point
+{
+    double _x;
+    double _y;
+
+    Point() : _x(0), _y(0)
+    { }
+
+    Point(double x, double y) : _x(x), _y(y)
+    { }
+
+    Point(const Point& pt) : _x(pt._x), _y(pt._y)
+    { }
+};
+
+typedef std::vector < Point > Points;
+
 struct VertexProperties
 {
     int    _index;
     double _weight;
-    double _x;
-    double _y;
+    Point  _centroid;
+    Points _neighbour_centroids;
 
-    VertexProperties() : _index(0), _weight(0), _x(0.), _y(0)
+    VertexProperties() : _index(0), _weight(0)
     { }
 
     VertexProperties(int index, double weight, double x, double y) :
-        _index(index), _weight(weight), _x(x), _y(y)
+        _index(index), _weight(weight), _centroid(x, y)
     { }
 
     VertexProperties(const VertexProperties& vp) :
-        _index(vp._index), _weight(vp._weight), _x(vp._x), _y(vp._y)
+        _index(vp._index), _weight(vp._weight), _centroid(vp._centroid),
+        _neighbour_centroids(vp._neighbour_centroids)
     { }
 };
 

src/tests/plot/graph_manager.hpp

@@ -75,35 +75,17 @@ public:
     void build_flat_graph(const OrientedGraph& g, const InputEdges& inputs)
     {
         OrientedGraph::vertex_iterator vertexIt, vertexEnd;
-        std::map < int, int > neighbour_numbers;
-
-        boost::tie(vertexIt, vertexEnd) = boost::vertices(g);
-        for (; vertexIt != vertexEnd; ++vertexIt)
-        {
-            OrientedGraph::adjacency_iterator neighbourIt, neighbourEnd;
-
-            neighbour_numbers[g[*vertexIt]._index] = 0;
-            boost::tie(neighbourIt, neighbourEnd) =
-                boost::adjacent_vertices(*vertexIt, g);
-            for (; neighbourIt != neighbourEnd; ++neighbourIt) {
-                ++neighbour_numbers[g[*vertexIt]._index];
-            }
-        }
-
-        for (Edges::const_iterator it = inputs.begin(); it != inputs.end();
-             ++it) {
-            ++neighbour_numbers[it->second];
-        }
 
         boost::tie(vertexIt, vertexEnd) = boost::vertices(g);
         for (; vertexIt != vertexEnd; ++vertexIt) {
             std::ostringstream ss;
+            PlotParameters parameters(g[*vertexIt]._index,
+                                      g[*vertexIt]._centroid,
+                                      g[*vertexIt]._neighbour_centroids);
 
-            ss << "a" << g[*vertexIt]._index;
-            _simulators[g[*vertexIt]._index] =
-                new Simulator(ss.str(),
-                              PlotParameters(g[*vertexIt]._x, g[*vertexIt]._y,
-                                  neighbour_numbers[g[*vertexIt]._index]));
+            ss << "" << g[*vertexIt]._index;
+            _simulators[g[*vertexIt]._index] = new Simulator(ss.str(),
+                                                             parameters);
             _simulators[g[*vertexIt]._index]->add_out_port("out");
             _simulators[g[*vertexIt]._index]->add_in_port("in");
             FlatGraphManager < Parameters >::add_child(

src/tests/plot/main.cpp

@@ -44,7 +44,7 @@ double plot_monothreading()
             HierarchicalGraphManager,
             paradevs::common::NoParameters,
             GraphManagerParameters >
-        > rc(0, 10, "root", paradevs::common::NoParameters(),
+        > rc(0, 60, "root", paradevs::common::NoParameters(),
              GraphManagerParameters("/home/eric/tmp/parcelle/oye_plage.shp",
                                     4));
 

src/tests/plot/milsol.hpp

@@ -0,0 +1,150 @@
+/**
+ * @file tests/plot/milsol.hpp
+ * @author The PARADEVS Development Team
+ * See the AUTHORS or Authors.txt file
+ */
+
+/*
+ * PARADEVS - the multimodeling and simulation environment
+ * This file is a part of the PARADEVS environment
+ *
+ * Copyright (C) 2013-2015 ULCO http://www.univ-litoral.fr
+ * Copyright (C) 2009 INRA
+ *
+ * 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/>.
+ */
+
+#ifndef TESTS_PLOT_MILSOL_HPP
+#define TESTS_PLOT_MILSOL_HPP 1
+
+#include <tests/plot/climate.hpp>
+#include <tests/plot/cohorte.hpp>
+#include <vector>
+
+namespace paradevs { namespace tests { namespace plot {
+
+class Milsol
+{
+public:
+    Milsol() : count(0), TotalSporeReady(0)
+    {
+        p1 = 0.5;
+        p2 = 1.0;
+        p3 = 10.0;
+        p4 = 0.05;
+        p5 = 0.025;
+        p6 = 1.5;
+        p7 = 75.0;
+        p8 = 150.0;
+        p9 = 225.0;
+        p10 = 0.004;
+        p11 = 0.009;
+        p12 = 2.0;
+        p13 = 1.0;
+        p14 = 0.037;
+        D0 = 100.0;
+        D1 = 150.0;
+        Dc = 100.0;
+        CUM0 = 6.0;
+        CUM1 = 10.0;
+        Topt = 18.0;
+        Tmin = 3.0;
+        FACT = 2000.0;
+        SSA = 1. / 17700;
+
+        _total_ready_spore_number_t = 0.0;
+        _total_ready_spore_number_t2 = 0.0;
+        _total_ready_spore_number_t3 = 0.0;
+
+        _zoospore_number = 0.0;
+    }
+
+    double get_ready_spore_number() const
+    { return TotalSporeReady; }
+
+    virtual void operator()(typename common::DoubleTime::type t)
+    {
+        _total_ready_spore_number_t = 0.0;
+
+        _climate.get_temperatures(t, _temperatures);
+        _climate.get_humidities(t, _humidities);
+
+        createCohorte();
+        for (unsigned int i = 0; i < _cohortes.size(); ++i) {
+            for (unsigned int j = 0; j < 24 ; ++j) {
+                _cohortes[i].compute(_temperatures[j], _humidities[j]);
+                _total_ready_spore_number_t += _cohortes[i].sporul();
+                if (condEndCohorte(i)) {
+                    _cohortes.erase(_cohortes.begin() + i);
+                    i--;
+                    break;
+                }
+            }
+        }
+        _total_ready_spore_number_t2 = _total_ready_spore_number_t -
+            _total_ready_spore_number_t3;
+        _total_ready_spore_number_t3 = _total_ready_spore_number_t;
+        TotalSporeReady = _total_ready_spore_number_t2;
+    }
+
+    void add_zoospore_number(double n)
+    { _zoospore_number += n; }
+
+    double get_zoospore_number() const
+    { return _zoospore_number; }
+
+    void set_inoculum_primaire_number(double n)
+    { _zoospore_number = n; }
+
+private:
+    bool condEndCohorte(int k)
+    { return (_cohortes[k].age() > p9); }
+
+    bool condBeginCohorte()
+    { return _zoospore_number > 0.0; }
+
+    void createCohorte()
+    {
+        if (condBeginCohorte()) {
+            ++count;
+            _cohortes.push_back(
+                Cohorte(count, _zoospore_number, p1, p2, p3, p4, p5,
+                        p6, p7,
+                        p8, p9, p10, p11, p12, p13, p14, D0, D1, Dc, CUM0,
+                        CUM1, Topt, Tmin, FACT, SSA));
+        }
+    }
+
+    double p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12, p13, p14;
+    double D0, D1, Dc;
+    double CUM0, CUM1;
+    double Topt, Tmin;
+    double FACT, SSA;
+
+    std::vector < Cohorte > _cohortes;
+    std::vector < double > _temperatures;
+    std::vector < double > _humidities;
+    int count;
+    double TotalSporeReady;
+    double _total_ready_spore_number_t;
+    double _total_ready_spore_number_t2;
+    double _total_ready_spore_number_t3;
+    double _zoospore_number;
+
+    Climate _climate;
+};
+
+} } } // namespace paradevs tests plot
+
+#endif

src/tests/plot/models.hpp

@@ -28,33 +28,44 @@
 #define TESTS_PLOT_MODELS_HPP 1
 
 #include <paradevs/common/time/DoubleTime.hpp>
-
 #include <paradevs/kernel/pdevs/Dynamics.hpp>
 
+#include <tests/plot/graph_defs.hpp>
+#include <tests/plot/dispersion.hpp>
+#include <tests/plot/milsol.hpp>
+
+#include <shapefil.h>
+
 namespace paradevs { namespace tests { namespace plot {
 
 struct PlotParameters
 {
-    PlotParameters(double x, double y, int n) : x(x), y(y), neighbour_number(n)
+    PlotParameters(int index,
+                   const paradevs::tests::boost_graph::Point& centroid,
+                   const paradevs::tests::boost_graph::Points&
+                   neighbour_centroids)
+        : _index(index), _centroid(centroid),
+          _neighbour_centroids(neighbour_centroids)
     { }
 
-    double x;
-    double y;
-    int neighbour_number;
+    int _index;
+    paradevs::tests::boost_graph::Point _centroid;
+    paradevs::tests::boost_graph::Points _neighbour_centroids;
 };
 
 struct PlotData
 {
-    PlotData() : x(-1), y(-1), quantity(0)
+    PlotData() : index(-1), x(-1), y(-1), ready_spore_number(0)
     { }
 
-    PlotData(double x, double y, double quantity) : x(x), y(y),
-                                                    quantity(quantity)
+    PlotData(int index, double x, double y, double ready_spore_number) :
+        index(index), x(x), y(y), ready_spore_number(ready_spore_number)
     { }
 
+    int index;
     double x;
     double y;
-    double quantity;
+    double ready_spore_number;
 };
 
 class Plot : public paradevs::pdevs::Dynamics < common::DoubleTime,
@@ -64,50 +75,173 @@ public:
     Plot(const std::string& name, const PlotParameters& parameters) :
         paradevs::pdevs::Dynamics < common::DoubleTime, PlotParameters >(
             name, parameters),
-        _neighbour_number(parameters.neighbour_number)
+        _index(parameters._index),
+        _centroid(parameters._centroid),
+        _neighbour_centroids(parameters._neighbour_centroids)
     {
-        _data.x = parameters.x;
-        _data.y = parameters.y;
+        _data.index = _index;
+        _data.x = parameters._centroid._x;
+        _data.y = parameters._centroid._y;
+
+        // std::cout << "(" << _centroid._x << "," << _centroid._y << "): ";
+        // for (Points::const_iterator it = _neighbour_centroids.begin();
+        //      it != _neighbour_centroids.end(); ++it) {
+        //     double distance = std::sqrt(
+        //         (_centroid._x - it->_x) * (_centroid._x - it->_x) +
+        //         (_centroid._y - it->_y) * (_centroid._y - it->_y));
+
+        //     std::cout << "(" << it->_x << "," << it->_y << ") -> " << distance
+        //               << " ";
+        // }
+        // std::cout << std::endl;
+
+        _handle = DBFOpen("/home/eric/tmp/parcelle/test/parcellaire.dbf",
+                          "rb+");
+
+        if (!_handle) {
+            std::cout << "Error" << std::endl;
+        }
+
     }
 
     virtual ~Plot()
-    { }
+    {
+        DBFClose(_handle);
+    }
 
-    void dint(typename common::DoubleTime::type /* t */)
+    void dint(typename common::DoubleTime::type t)
     {
+        // std::cout << _index << " " << t << ": dint - BEFORE - "
+        //           << _phase << std::endl;
+
+        if (_phase == SEND) {
+            if (_neighbour_centroids.size() > 0) {
+                _phase = WAIT;
+                _sigma = common::DoubleTime::infinity;
+            } else {
+                _phase  = NEWSTATE;
+                _sigma = 0;
+            }
+        } else if (_phase == NEWSTATE) {
+            // cohorte update
+            _milsol(t);
+            _data.ready_spore_number = _milsol.get_ready_spore_number();
+
+            // compute
+            double sum = 0;
+
+            for (std::vector < PlotData >::const_iterator it =
+                     _neighbour_data.begin(); it != _neighbour_data.end();
+                 ++it) {
+                if (it->ready_spore_number > 0) {
+                    sum += _dispersion_function(
+                        t,
+                        paradevs::tests::boost_graph::Point(it->x, it->y),
+                        _centroid, 0., 0.) * it->ready_spore_number;
+
+                    // std::cout << (_index + 1) << " " << (it->index + 1) << " " << it->concentration << std::endl;
+
+                }
+            }
+            sum += _dispersion_function(t, _centroid, _centroid, 0., 0.) *
+                _data.ready_spore_number;
+            _milsol.add_zoospore_number(sum);
+
+            double p = 0;
+
+            if (_milsol.get_zoospore_number() > 0) {
+                p =  std::log10(_milsol.get_zoospore_number());
+            }
+
+            // std::cout << t << "\t" << get_name() << "\t"
+            //           << _data.ready_spore_number << "\t"
+            //           << _milsol.get_zoospore_number() << "\t"
+            //           << p << std::endl;
+
+            DBFWriteDoubleAttribute(_handle, _index, 0, p);
+
+            // std::cout << (_index + 1) << " " << _data.concentration << std::endl;
+
+            _phase  = SEND;
+            _sigma = 1;
+            _received = 0;
+            _neighbour_data.clear();
+        }
+
+        // std::cout << _index << " " << t << ": dint - AFTER - "
+        //           << _phase << std::endl;
+
     }
 
-    void dext(typename common::DoubleTime::type /* t */,
+    void dext(typename common::DoubleTime::type t,
               typename common::DoubleTime::type /* e */,
-              const common::Bag < common::DoubleTime >& /* bag */)
+              const common::Bag < common::DoubleTime >& bag)
     {
+        for (common::Bag < common::DoubleTime >::const_iterator it =
+                 bag.begin(); it != bag.end(); ++it) {
+            _neighbour_data.push_back(*(PlotData*)(it->get_content()));
+            ++_received;
+        }
+
+        if (_received == _neighbour_centroids.size()) {
+            _phase = NEWSTATE;
+            _sigma = 0;
+        } else {
+            _phase = WAIT;
+            _sigma = common::DoubleTime::infinity;
+        }
+
+        // std::cout << _index << " " << t << ": dext - "
+        //           << _received << "/"
+        //           << _neighbour_centroids.size()
+        //           << " " << _phase << std::endl;
     }
 
-    void dconf(typename common::DoubleTime::type /* t */,
-               typename common::DoubleTime::type /* e */,
-               const common::Bag < common::DoubleTime >& /* bag */)
+    void dconf(typename common::DoubleTime::type t,
+               typename common::DoubleTime::type e,
+               const common::Bag < common::DoubleTime >& bag)
     {
-        // dext(t, e, bag);
+
+        // std::cout << t << "[" << get_name() << "]: dconf" << std::endl;
+
+        dext(t, e, bag);
     }
 
     typename common::DoubleTime::type start(
         typename common::DoubleTime::type /* t */)
     {
-        _data.quantity = 0;
+        if (_index == 210) {
+            _milsol.set_inoculum_primaire_number(10);
+        } else {
+            _milsol.set_inoculum_primaire_number(0);
+        }
+        _phase = SEND;
+        if (_neighbour_centroids.size() > 0) {
+            _sigma = 0;
+        } else {
+            _sigma = 1;
+        }
+        _received = 0;
         return 0;
     }
 
     typename common::DoubleTime::type ta(
         typename common::DoubleTime::type /* t */) const
-    { return common::DoubleTime::infinity; }
+    { return _sigma; }
 
     common::Bag < common::DoubleTime > lambda(
-        typename common::DoubleTime::type /* t */) const
+        typename common::DoubleTime::type t) const
     {
         common::Bag < common::DoubleTime > bag;
 
-        bag.push_back(common::ExternalEvent < common::DoubleTime >(
-                          "out", (void*)(&_data)));
+        if (_phase == SEND) {
+            bag.push_back(common::ExternalEvent < common::DoubleTime >(
+                              "out", (void*)(&_data)));
+
+            // std::cout << t << "[" << get_name() << "]: lambda - "
+            //           << _phase << std::endl;
+
+        }
         return bag;
     }
 
@@ -115,11 +249,28 @@ public:
     { }
 
 private:
+    enum Phase { WAIT, SEND, NEWSTATE, NEXT };
+
     // parameters
-    unsigned int _neighbour_number;
+    int _index;
+    paradevs::tests::boost_graph::Point _centroid;
+    std::vector < paradevs::tests::boost_graph::Point > _neighbour_centroids;
+
+    DBFHandle _handle;
 
     // data
     PlotData _data;
+    std::vector < PlotData > _neighbour_data;
+
+    // submodels
+    KleinDispersionFunction _dispersion_function;
+    Milsol                  _milsol;
+    Climate                 _climate;
+
+    // state
+    Phase _phase;
+    common::DoubleTime::type _sigma;
+    int _received;
 };
 
 } } } // namespace paradevs tests plot