Merge dev for preparing release v0.4.0

Project.toml

@@ -14,8 +14,8 @@ PlotlyJS = "f0f68f2c-4968-5e81-91da-67840de0976a"
 [compat]
 ColorSchemes = "3.0 - 3.19"
 ColorTypes = "0.7 - 0.11"
-DifferentialEquations = "7.1 - 7.5"
+DifferentialEquations = "7.1 - 7.6"
 IfElse = "0.1"
-ModelingToolkit = "8.26 - 8.27"
+ModelingToolkit = "8.26 - 8.29"
 PlotlyJS = "0.13 - 0.18"
 julia = "1.6 - 1.8"

docs/src/source.md

@@ -1,5 +1,12 @@
 # Source code documentation
 
+## Contents
+
+```@contents
+Pages = ["source.md"]
+Depth = 3
+```
+
 ## WorldDynamics constants
 ```@autodocs
 Modules = [WorldDynamics]
@@ -34,6 +41,15 @@ Order   = [:function]
 Pages   = ["solvesystems.jl"]
 ```
 
+## Reproducing World2 figures
+
+### World2 system
+
+```@autodocs
+Modules = [WorldDynamics.World2]
+Pages   = ["plots.jl"]
+```
+
 ## Reproducing World3 figures
 
 ### Agriculture system
@@ -57,6 +73,13 @@ Modules = [WorldDynamics.World3.NonRenewable]
 Pages   = ["plots.jl"]
 ```
 
+### Pollution system
+
+```@autodocs
+Modules = [WorldDynamics.World3.Pollution]
+Pages   = ["plots.jl"]
+```
+
 ### One level population system
 
 ```@autodocs

docs/src/tutorial.md

@@ -1,6 +1,6 @@
 # A WorldDynamics tutorial
 
-`WorldDynamics` allows the user to *play* with the World3 model introduced in the book *Dynamics of Growth in a Finite World* (1974). Informally speaking, this model is formed by five systems, each containg one or more subsystems. The following picture shows the structure of the model and the connections between the subsystems which share a common variable.
+`WorldDynamics` allows the user to *play* with the World3 model introduced in the book *Dynamics of Growth in a Finite World* (1974). Informally speaking, this model is formed by five systems, each containing one or more subsystems. The following picture shows the structure of the model and the connections between the subsystems which share a common variable.
 
 ![The World3 model](img/world3.png)
 
@@ -103,7 +103,7 @@ plotvariables(sol, (t, 1900, 2100), reference_variables, title="Fig. 7-10", show
 
 ### Modifying an interpolation table
 
-In order to reproduce Figure 7-13, in which the slope of the fraction of industrial output allocated to agriculture is increased, we can modify the two tables `FIOAA1` and `FIOAA2` by getting the table set of the agriculture sector, and by changing the value of these two tables. We then have to solve again the ODE system, by specifying which set of tables has to be used for the agriculture sector. Finally, we can plot the same seven variables of Figure 7-10. This is exactly what we do in the following code.
+In order to reproduce Figure 7-13, in which the slope of the fraction of industrial output allocated to agriculture is increased, we can modify the two tables `FIOAA1` and `FIOAA2` by getting the table set of the agriculture sector, and by changing the value of these two tables. We then have to solve the ODE system again, by specifying which set of tables has to be used for the agriculture sector. Finally, we can plot the same seven variables of Figure 7-10. This is exactly what we do in the following code.
 
 ```
 using WorldDynamics
@@ -115,3 +115,27 @@ system = World3.historicalrun(agriculture_tables=agriculture_tables_7_13);
 sol = WorldDynamics.solve(system, (1900, 2100));
 plotvariables(sol, (t, 1900, 2100), reference_variables, title="Fig. 7-13", showlegend=true, colored=true)
 ```
+
+### Updating the model with modern data
+The flexible structure of `WorldDynamics` allows the user to feed the model with modern data. For example, in the book, the variable `POP` of the pollution system is assigned the following interpolation table which corresponds to the population number (expressed in $10^8$) for a set of years between 1900 and 2100.
+
+```
+tables[:pop] = (16.0, 19.0, 22.0, 31.0, 42.0, 53.0, 67.0, 86.0, 109.0, 139.0, 176.0);
+ranges[:pop] = (1900, 2100)
+```
+
+Instead, we can extend the above set of years to a much larger one as well as replace any outdated estimations with more recent data available at open-source data catalogs. In the following, we consider past and future projections of the world population, taken from the recognized public database
+[Our World In Data](https://ourworldindata.org/grapher/population-past-future?time=1867..2100). We first have to modify the table `POP` by getting the table set of the pollution sector, and by changing its value. We then have to solve the ODE system again, by specifying which set of tables has to be used for the pollution sector. This is exactly what we do in the following code.
+
+
+```
+using WorldDynamics
+
+tables = Pollution.gettables();
+tables[:pop] = (16.47,16.59,16.73,16.87,17.02,17.16,17.31,17.47,17.62,17.77,17.93,18.05,18.18,18.30,18.43,18.56,18.69,18.82,18.95,19.09,19.26,19.40,19.56,19.73,19.90,20.08,20.26,20.44,20.63,20.82,21.04,21.22,21.44,21.66,21.88,22.10,22.33,22.57,22.80,23.03,23.27,23.45,23.64,23.82,24.00,24.17,24.35,24.54,24.75,25.01,24.99,25.43,25.90,26.40,26.92,27.46,28.01,28.58,29.16,29.70,30.19,30.68,31.27,31.96,32.67,33.37,34.06,34.75,35.47,36.21,36.95,37.70,38.45,39.20,39.96,40.69,41.43,42.16,42.90,43.66,44.44,45.25,46.08,46.92,47.76,48.62,49.50,50.41,51.32,52.24,53.16,54.06,54.93,55.77,56.61,57.43,58.25,59.06,59.87,60.68,61.49,62.31,63.12,63.94,64.76,65.58,66.41,67.26,68.12,68.98,69.86,70.73,71.62,72.51,73.39,74.27,75.13,76.00,76.84,77.65,78.41,79.09,79.75,80.45,81.19,81.92,82.64,83.36,84.07,84.77,85.46,86.15,86.82,87.49,88.15,88.79,89.43,90.06,90.68,91.29,91.88,92.47,93.04,93.60,94.14,94.68,95.19,95.69,96.18,96.65,97.09,97.53,97.94,98.34,98.72,99.08,99.43,99.76,100.08,100.39,100.68,100.96,101.22,101.48,101.73,101.96,102.18,102.40,102.60,102.79,102.97,103.14,103.30,103.45,103.59,103.71,103.82,103.92,104.01,104.08,104.15,104.20,104.24,104.27,104.29,104.31,104.31,104.30,104.29,104.27,104.24,104.20,104.15,104.09,104.03,103.96,103.89,103.80,103.70,103.60,103.49);
+
+system = World3.Pollution.historicalrun(tables=tables);
+sol = WorldDynamics.solve(system, (1900, 2100))
+```
+
+Finally, we can compare the model updated with new data against the one with outdated data by reproducing the figures from the book (as described within the [Replicating book runs](@ref Replicating-book-runs) section).

docs/src/world2.md

@@ -1,4 +1,4 @@
-# [World 2 equations, variables, and parameters](@id eqs_vars_pars_2) 
+# [World 2 equations, variables, and parameters](@id eqs_vars_pars_2)
 
 In this page we list the equations, the variables and the parameters of the World2 model as described in Chapter 3 and Appendix B of the book *World dynamics* (1973). Even if it is not said explicitily in the book, the World2 model consists of six systems containing several subsystems.
 
@@ -78,7 +78,7 @@ In this page we list the equations, the variables and the parameters of the Worl
 | Effective capital investment ratio | `ecir` | `capital_investment` |  |
 | Food from crowding multiplier | `fcm` | `agriculture_investment` |  |
 | Food potential from capital investment | `fpci` | `agriculture_investment` |  |
-| Food from pollution multiplier | `fpm` | `agriculture_investment` |  | 
+| Food from pollution multiplier | `fpm` | `agriculture_investment` |  |
 | Food ratio | `fr` | `agriculture_investment` | `birth_rate`, `death_rate`, `quality_life` |
 | Material standard of living | `msl` | `capital_investment` | `birth_rate`, `capital_investment_generation`, `death_rate`, `natural_resources_usage_rate`, `quality_life` |
 | Natural resources | `nr` | `natural_resources` |  |
@@ -139,7 +139,7 @@ In this page we list the equations, the variables and the parameters of the Worl
 | Switch time no. 3 for `nrun` | $\mathtt{swt3}$ | $1970$ | `Population` |
 | Switch time no. 4 for `cign` | $\mathtt{swt4}$ | $1970$ | `Capital investment` |
 | Switch time no. 5 for `cidn` | $\mathtt{swt5}$ | $1970$ | `Capital investment` |
-| Switch time no. 5 for `poln` | $\mathtt{swt6}$ | $1970$ | `Pollution` |
+| Switch time no. 6 for `poln` | $\mathtt{swt6}$ | $1970$ | `Pollution` |
 | Switch time no. 7 for `fc` | $\mathtt{swt7}$ | $1970$ | `Population` |
 
 ## Tables and ranges

src/World2/agricultureinvestment/subsystems.jl

@@ -2,25 +2,25 @@
 D = Differential(t)
 
 function agriculture_investment(; name, params=_params, inits=_inits, tables=_tables, ranges=_ranges)
-    @parameters fc = params[:fc]
-    @parameters fc1 = params[:fc1]
-    @parameters swt7 = params[:swt7]
-    @parameters fn = params[:fn]
-    @parameters ciaft = params[:ciaft]
+    @parameters fc = params[:fc] [description = "Food coefficient"]
+    @parameters fc1 = params[:fc1] [description = "Food coefficient no. 1"]
+    @parameters swt7 = params[:swt7] [description = "Switch time no. 7"]
+    @parameters fn = params[:fn] [description = "Food normal"]
+    @parameters ciaft = params[:ciaft] [description = "Capital investment in agriculture fraction adjustment time"]
 
-    @variables fr(t)
-    @variables fpci(t)
-    @variables fcm(t)
-    @variables fpm(t)
-    @variables ciaf(t) = inits[:ciaf]
-    @variables cfifr(t)
-    @variables ciqr(t)
+    @variables fr(t) [description = "Food ratio"]
+    @variables fcm(t) [description = "Food from crowding multiplier"]
+    @variables fpci(t) [description = "Food potential from capital investment"]
+    @variables fpm(t) [description = "Food from pollution multiplier"]
+    @variables ciaf(t) = inits[:ciaf] [description = "Capital investment in agriculture fraction"]
+    @variables cfifr(t) [description = "Capital fraction indicated by food ratio"]
+    @variables ciqr(t) [description = "Capital investment from quality ratio"]
 
-    @variables cr(t)
-    @variables cira(t)
-    @variables polr(t)
-    @variables qlm(t)
-    @variables qlf(t)
+    @variables cr(t) [description = "Crowding ratio"]
+    @variables cira(t) [description = "Capital investment ratio in agriculture"]
+    @variables polr(t) [description = "Pollution ratio"]
+    @variables qlm(t) [description = "Quality of life from material"]
+    @variables qlf(t) [description = "Quality of life from food"]
 
     eqs = [
         fr ~ fpci * fcm * fpm * clip(fc, fc1, swt7, t) / fn

src/World2/capitalinvestment/subsystems.jl

@@ -2,20 +2,20 @@
 D = Differential(t)
 
 function capital_investment(; name, params=_params, inits=_inits, tables=_tables, ranges=_ranges)
-    @parameters ecirn = params[:ecirn]
-    @parameters ciafn = params[:ciafn]
+    @parameters ecirn = params[:ecirn] [description = "Effective capital investment ratio normal"]
+    @parameters ciafn = params[:ciafn] [description = "Capital investment in agriculture fraction normal"]
 
-    @variables msl(t)
-    @variables ecir(t)
-    @variables cir(t)
-    @variables cira(t)
-    @variables ci(t) = inits[:ci]
+    @variables msl(t) [description = "Material standard of living"]
+    @variables ecir(t) [description = "Effective capital investment ratio"]
+    @variables cira(t) [description = "Capital investment ratio in agriculture"]
+    @variables cir(t) [description = "Capital investment ratio"]
+    @variables ci(t) = inits[:ci] [description = "Capital investment"]
 
-    @variables ciaf(t)
-    @variables nrem(t)
-    @variables p(t)
-    @variables cig(t)
-    @variables cid(t)
+    @variables ciaf(t) [description = "Capital investment in agriculture fraction"]
+    @variables nrem(t) [description = "Natural resource extraction multiplier"]
+    @variables p(t) [description = "Population"]
+    @variables cig(t) [description = "Capital investment generation"]
+    @variables cid(t) [description = "Capital investment discard"]
 
     eqs = [
         msl ~ ecir / ecirn
@@ -29,15 +29,15 @@ function capital_investment(; name, params=_params, inits=_inits, tables=_tables
 end
 
 function capital_investment_generation(; name, params=_params, inits=_inits, tables=_tables, ranges=_ranges)
-    @parameters cign = params[:cign]
-    @parameters cign1 = params[:cign1]
-    @parameters swt4 = params[:swt4]
+    @parameters cign = params[:cign] [description = "Capital investment generation normal"]
+    @parameters cign1 = params[:cign1] [description = "Capital investment generation normal no. 1"]
+    @parameters swt4 = params[:swt4] [description = "Switch time no. 4"]
 
-    @variables cig(t)
-    @variables cim(t)
+    @variables cig(t) [description = "Capital investment generation"]
+    @variables cim(t) [description = "Capital investment multiplier"]
 
-    @variables p(t)
-    @variables msl(t)
+    @variables p(t) [description = "Population"]
+    @variables msl(t) [description = "Material standard of living"]
 
     eqs = [
         cig ~ p * cim * clip(cign, cign1, swt4, t)
@@ -48,13 +48,13 @@ function capital_investment_generation(; name, params=_params, inits=_inits, tab
 end
 
 function capital_investment_discard(; name, params=_params, inits=_inits, tables=_tables, ranges=_ranges)
-    @parameters cidn = params[:cidn]
-    @parameters cidn1 = params[:cidn1]
-    @parameters swt5 = params[:swt5]
+    @parameters cidn = params[:cidn] [description = "Capital investment discard normal"]
+    @parameters cidn1 = params[:cidn1] [description = "Capital investment discard normal no. 1"]
+    @parameters swt5 = params[:swt5] [description = "Switch time no. 5"]
 
-    @variables cid(t)
+    @variables cid(t) [description = "Capital investment discard"]
 
-    @variables ci(t)
+    @variables ci(t) [description = "Capital investment"]
 
     eqs = [
         cid ~ ci * clip(cidn, cidn1, swt5, t)

src/World2/naturalresources/subsystems.jl

@@ -2,11 +2,11 @@
 D = Differential(t)
 
 function natural_resources(; name, params=_params, inits=_inits, tables=_tables, ranges=_ranges)
-    @variables nrem(t)
-    @variables nrfr(t)
-    @variables nr(t) = inits[:nr]
+    @variables nrem(t) [description = "Natural resource extraction multiplier"]
+    @variables nrfr(t) [description = "Natural resource fraction remaining"]
+    @variables nr(t) = inits[:nr] [description = "Natural resources"]
 
-    @variables nrur(t)
+    @variables nrur(t) [description = "Natural resource usage rate"]
 
     eqs = [
         nrem ~ interpolate(nrfr, tables[:nrem], ranges[:nrem])
@@ -18,15 +18,15 @@ function natural_resources(; name, params=_params, inits=_inits, tables=_tables,
 end
 
 function natural_resources_usage_rate(; name, params=_params, inits=_inits, tables=_tables, ranges=_ranges)
-    @parameters nrun = params[:nrun]
-    @parameters nrun1 = params[:nrun1]
-    @parameters swt2 = params[:swt2]
+    @parameters nrun = params[:nrun] [description = "Natural resource usage normal"]
+    @parameters nrun1 = params[:nrun1] [description = "Natural resource usage normal no. 1"]
+    @parameters swt2 = params[:swt2] [description = "Switch time no. 2"]
 
-    @variables nrur(t)
-    @variables nrmm(t)
+    @variables nrur(t) [description = "Natural resource usage rate"]
+    @variables nrmm(t) [description = "Natural resources from material multiplier"]
 
-    @variables p(t)
-    @variables msl(t)
+    @variables p(t) [description = "Population"]
+    @variables msl(t) [description = "Material standard of living"]
 
     eqs = [
         nrur ~ p * clip(nrun, nrun1, swt2, t) * nrmm

src/World2/pollution/subsystems.jl

@@ -2,13 +2,13 @@
 D = Differential(t)
 
 function pollution(; name, params=_params, inits=_inits, tables=_tables, ranges=_ranges)
-    @parameters pols = params[:pols]
+    @parameters pols = params[:pols] [description = "Pollution standard"]
 
-    @variables polr(t)
-    @variables pol(t) = inits[:pol]
+    @variables polr(t) [description = "Pollution ratio"]
+    @variables pol(t) = inits[:pol] [description = "Pollution"]
 
-    @variables polg(t)
-    @variables pola(t)
+    @variables polg(t) [description = "Pollution generation	"]
+    @variables pola(t) [description = "Pollution absorption"]
 
     eqs = [
         polr ~ pol / pols
@@ -19,11 +19,11 @@ function pollution(; name, params=_params, inits=_inits, tables=_tables, ranges=
 end
 
 function pollution_absorption(; name, params=_params, inits=_inits, tables=_tables, ranges=_ranges)
-    @variables pola(t)
-    @variables polat(t)
+    @variables pola(t) [description = "Pollution absorption"]
+    @variables polat(t) [description = "Pollution absorption time"]
 
-    @variables pol(t)
-    @variables polr(t)
+    @variables pol(t) [description = "Pollution"]
+    @variables polr(t) [description = "Pollution ratio"]
 
     eqs = [
         pola ~ pol / polat
@@ -34,15 +34,15 @@ function pollution_absorption(; name, params=_params, inits=_inits, tables=_tabl
 end
 
 function pollution_generation(; name, params=_params, inits=_inits, tables=_tables, ranges=_ranges)
-    @parameters poln = params[:poln]
-    @parameters poln1 = params[:poln1]
-    @parameters swt6 = params[:swt6]
+    @parameters poln = params[:poln] [description = "Pollution normal"]
+    @parameters poln1 = params[:poln1] [description = "Pollution normal no. 1"]
+    @parameters swt6 = params[:swt6] [description = "Switch time no. 6"]
 
-    @variables polg(t)
-    @variables polcm(t)
+    @variables polg(t) [description = "Pollution generation"]
+    @variables polcm(t) [description = "Pollution from capital multiplier"]
 
-    @variables p(t)
-    @variables cir(t)
+    @variables p(t) [description = "Population"]
+    @variables cir(t) [description = "Capital investment ratio"]
 
     eqs = [
         polg ~ p * clip(poln, poln1, swt6, t) * polcm