Rework epid model

Closes #117

doc/application/plot_metamodel.py

@@ -51,7 +51,6 @@
 
 path_fmu = otfmi.example.utility.get_path_fmu("epid")
 mesh = ot.RegularGrid(0.0, 0.05, 20)
-meshSample = mesh.getVertices()
 
 function = otfmi.FMUPointToFieldFunction(
     mesh,
@@ -67,7 +66,7 @@
 # simulate the FMU.
 # The simulation inputs and outputs will be used to train the metamodel.
 
-inputLaw = ot.Uniform(0.001, 0.01)
+inputLaw = ot.Uniform(1.5, 2.5)
 inputSample = inputLaw.getSample(30)
 outputFMUSample = function(inputSample)
 
@@ -201,7 +200,7 @@ def globalMetamodel(sample):
 ot.Show(graph)
 
 # %%
-# As the epidemiological model considers a population size of 700, the residual
+# As the epidemiological model considers a population size of 763, the residual
 # mean error on the field is acceptable.
 
 # %%

doc/example/dynamic/plot_dyn_basics.py

@@ -29,7 +29,7 @@
 # Define the time grid for the FMU's output. The last value of the time grid,
 # here 10., will define the FMU stop time for simulation.
 
-mesh = ot.RegularGrid(0.0, 0.1, 100)
+mesh = ot.RegularGrid(0.0, 0.1, 2000)
 meshSample = mesh.getVertices()
 print(meshSample)
 
@@ -47,7 +47,7 @@
     inputs_fmu=["infection_rate"],
     outputs_fmu=["infected"],
     start_time=0.0,
-    final_time=10.0,
+    final_time=200.0,
 )
 print(type(function))
 
@@ -61,7 +61,7 @@
 # Simulate the function on an input :py:class:`openturns.Point` yields an output
 # :py:class:`openturns.Sample`, corresponding to the output evolution over time:
 
-inputPoint = ot.Point([0.007])
+inputPoint = ot.Point([2.0])
 outputSample = function(inputPoint)
 
 plt.xlabel("FMU simulation time (s)")
@@ -73,7 +73,7 @@
 # Simulate the function on a input :py:class:`openturns.Sample` yields a set of
 # fields called :py:class:`openturns.ProcessSample`:
 
-inputSample = ot.Sample([[0.007], [0.005], [0.003]])
+inputSample = ot.Sample([[2.0], [2.25], [2.5]])
 outputProcessSample = function(inputSample)
 print(outputProcessSample)
 

doc/example/dynamic/plot_dyn_init.py

@@ -42,7 +42,7 @@
 temporary_file = "initialization.mos"
 with open(temporary_file, "w") as f:
     f.write("total_pop = 500;\n")
-    f.write("healing_rate = 0.02;\n")
+    f.write("healing_rate = 0.5;\n")
 
 # %%
 # If no initial value is provided for an input / parameter, it is set to its
@@ -79,8 +79,8 @@
 # function input variable ``infection_rate`` to propagate its uncertainty
 # through the model:
 
-lawInfected = ot.Normal(0.01, 0.003)
-inputSample = lawInfected.getSample(10)
+law_infection_rate = ot.Normal(2.0, 0.25)
+inputSample = law_infection_rate.getSample(10)
 outputProcessSample = function(inputSample)
 
 # %%

doc/fmus/epid.rst

@@ -24,8 +24,8 @@ time writes:
 .. math::
 
    \begin{aligned}
-   \frac{\partial S}{\partial t}(t) &= - \beta S(t) I(t) \\
-   \frac{\partial I}{\partial t}(t) &= \beta S(t) I(t) - \gamma I(t) \\
+   \frac{\partial S}{\partial t}(t) &= - \frac{\beta}{N} S(t) I(t) \\
+   \frac{\partial I}{\partial t}(t) &= \frac{\beta}{N} S(t) I(t) - \gamma I(t) \\
    \frac{\partial R}{\partial t}(t) &= \gamma I(t)
    \end{aligned}
 
@@ -35,25 +35,27 @@ This model is implemented in Modelica language. The default simulation time is 5
 
    model epid
 
-   parameter Real total_pop = 700;
+   parameter Real total_pop = 763;
+   parameter Real infection_rate = 2.0;
+   parameter Real healing_rate = 0.5;
+
    Real infected;
    Real susceptible;
    Real removed;
-   parameter Real infection_rate = 0.007;
-   parameter Real healing_rate = 0.02;
 
    initial equation
    infected = 1;
    removed = 0;
    total_pop = infected + susceptible + removed;
 
    equation
-   der(susceptible) = - infection_rate*infected*susceptible;
-   der(infected) = infection_rate*infected*susceptible - healing_rate*infected;
-   der(removed) = healing_rate*infected;
+   der(susceptible) = - infection_rate * infected * susceptible / total_pop;
+   der(infected) = infection_rate * infected * susceptible / total_pop -
+                   healing_rate * infected;
+   der(removed) = healing_rate * infected;
 
    annotation(
-       experiment(StartTime = 0, StopTime = 50, Tolerance = 1e-6, Interval = 0.1));
+       experiment(StartTime = 0.0, StopTime = 200.0, Tolerance = 1e-6, Interval = 0.1));
    end epid;
 
 We focus on the effect of the ``infection_rate`` and ``healing_rate`` on the evolution of the ``infected`` category.

otfmi/example/file/epid.mo

@@ -1,24 +1,24 @@
 model epid
 
-parameter Real total_pop = 700;
+parameter Real total_pop = 763;
+parameter Real infection_rate = 2.0;
+parameter Real healing_rate = 0.5;
 
 Real infected;
 Real susceptible;
 Real removed;
 
-parameter Real infection_rate = 0.007;
-parameter Real healing_rate = 0.02;
-
 initial equation
 infected = 1;
 removed = 0;
 total_pop = infected + susceptible + removed;
 
 equation
-der(susceptible) = - infection_rate*infected*susceptible;
-der(infected) = infection_rate*infected*susceptible - healing_rate*infected;
-der(removed) = healing_rate*infected;
+der(susceptible) = - infection_rate * infected * susceptible / total_pop;
+der(infected) = infection_rate * infected * susceptible / total_pop -
+                healing_rate * infected;
+der(removed) = healing_rate * infected;
 
 annotation(
-    experiment(StartTime = 0, StopTime = 50, Tolerance = 1e-6, Interval = 0.1));
+    experiment(StartTime = 0.0, StopTime = 200.0, Tolerance = 1e-6, Interval = 0.1));
 end epid;

test/test_epid.py

@@ -7,7 +7,9 @@
 import pytest
 
 
-input_value = [0.007, 0.02]
+# infection_rate, healing_rate
+input_value = [2.0, 0.5]
+final_time = 5
 
 
 def simulate_with_pyfmi(path_fmu, final_time=None):
@@ -24,7 +26,7 @@ def simulate_with_pyfmi(path_fmu, final_time=None):
     list_variable = list(check_model.get_model_variables().keys())
     infected_column = list_variable.index("infected") + 1
     # +1 stands for time column, inserted as first in the data matrix
-    list_last_infected_value = check_result.data_matrix[infected_column][-5:-1]
+    list_last_infected_value = check_result.data_matrix[infected_column][-5:]
     return list_last_infected_value
 
 
@@ -37,23 +39,13 @@ def path_fmu():
 def model_fmu(path_fmu):
     return otfmi.FMUFunction(path_fmu,
                              inputs_fmu=["infection_rate", "healing_rate"],
-                             outputs_fmu=["infected"])
+                             outputs_fmu=["infected"], final_time=final_time)
 
 
 def test_final_value(model_fmu):
     """Check reproducibility of the final value."""
     y = model_fmu(input_value)
-    assert abs(y[0] - 265.68) < 1e-2, "wrong value"
-
-
-def test_fmufunction_interpolation(path_fmu, model_fmu):
-    """Check that FMUFunction returns a point corresponding to Pyfmi's
-    interpolated output.
-    """
-    list_last_infected_value = simulate_with_pyfmi(path_fmu)
-    y = model_fmu(input_value)
-    assert y[0] < max(list_last_infected_value)
-    assert y[0] > min(list_last_infected_value)
+    assert abs(y[0] - 303.785) < 1e-2, "wrong value"
 
 
 def test_final_time(path_fmu):
@@ -69,8 +61,7 @@ def test_final_time(path_fmu):
         path_fmu, final_time=final_time
     )
     y = model_fmu_1(input_value)
-    assert y[0] < max(list_last_infected_value)
-    assert y[0] > min(list_last_infected_value)
+    assert abs(y[0] - list_last_infected_value[-1]) < 1e-5
 
 
 def test_keyword(path_fmu):