worlddynamics / worlddynamics.jl Goto Github PK

Inspired by the really nice overlap figure from @cvanwynsberghe in his Python implementation of World3, we can create similar ones manually. It would be nice to automate the process for every figure.

In Fig. 3.46, there is a slight difference for the sopc variable. Its corresponding plot should exhibit a higher peak around 1980.

The issue persists when changing the solver.

    The snippet below is an example on how to reproduce a figure currently implemented.

using WorldDynamics
fig = WorldDynamics.World3.fig_7()

Originally posted by @paulobruno in #63 (comment)

Once a system is solved, we should be able to save the solution in a file (JLD2, for example). Then, when we need the same solution, instead of solving again, we just load the file.

At the moment, all plots of a given subsystem have to be saved manually. Moreover, within a series of called-plots, only the last one is visualised.
We should be able to have them all stored (locally) at once, when first generated.

As mentioned in #18, the current way to change an equation is by printing the list of equations, getting the index of the one we would like to edit, and finally change it. This is obviously not user-friendly.

Here, the documentation mentions a picture which is not shown.
In fact, there is no link to any picture in the source of the tutorial page.

We should provide a bibtex entry to cite this work.

The recommended way of using an if-else check in ModellingToolkit is to use the IfElse.ifelse function.
However, this is not possible in our interpolate function, because we can potentially evaluate a position outside of the range. If this is the case, we get a BoundsError.

Currently, our workaround is to use Tuples, which allows the interpolate function to be called only after the variables were already evaluated.

World2 is already available in WorldDynamics.jl, but currently we only have fig 4.1 implemented.

In few figures of Chapter 2 (Population), the plot corresponding to the cbr variable fluctuates. This instability might be related to Issue #59 (e.g., iphst parameter).

In the following, we provide a minimal reproducible example.

Source code :

using WorldDynamics

include("../scenarios/pop1_historicalrun.jl")

@named pop = WorldDynamics.World3.Pop1.population()
@named br = WorldDynamics.World3.Pop1.birth_rate()
@named dr = WorldDynamics.World3.Pop1.death_rate()
@named so = WorldDynamics.World3.Pop1.service_output()
@named io = WorldDynamics.World3.Pop1.industrial_output()
@named f = WorldDynamics.World3.Pop1.food()

@variables t

fig_2_88_variables = [
    (so.sopc, 0,   1000,   "sopc"),
    (io.iopc, 0,   1000,   "iopc"),
    (f.fpc,   0,   1000,   "fpc"),
    (pop.pop, 0,   1.6e10, "pop"),
    (br.cbr,  0,   50,     "cbr"),
    (dr.cdr,  0,   50,     "cdr"),
    (dr.le,   0,   80,     "le"),
    (dr.fpu,  0,   1,      "fpu"),
    (br.fce,  0.5, 1,      "fce"),
]

parameters_2_88 = WorldDynamics.World3.Pop1.getparameters()
parameters_2_88[:iphst] = 4000
parameters_2_88[:lt2] = 1900
parameters_2_88[:cio] = 1000
parameters_2_88[:cso] = 1500
parameters_2_88[:cfood] = 2500

system = pop1_historicalrun(params=parameters_2_88)
sol_2_88 = WorldDynamics.solve(system, (1900, 2100))

plotvariables(sol_2_88, (t, 1900, 2100), fig_2_88_variables, name="Fig. 2.88a", showlegend=true, showaxis=true, colored=true)

Resulting plot :

Abstract types make the containers hold all possible subtypes. Not only this is very harmful to performance, since we are dealing with tuples, compilation time is also affected by them. This is also recommended by Julia Performance Tips.

My suggestion is to use Float64, which seems to be the common type used in SciML code.

It would be good to provide, as artifacts, pre-computed solutions that can e.g. be used to speed up the plotting of figures.

Related to #24.

We should also acknowledge support by:

• I3S, AAPI 02-2022,
• UCA Academy of Excellence RISE, Programme de chercheurs invités 2022,
• IDEX UCAjedi, Appel à candidatures pour Professeurs Invités UCA.

The logos should appear in the README.

A template will make it easier for external collaborators to create issues and pull requests containing all needed information.

At the moment, the tables are obtained from the book Dynamics of Growth in a Finite World (1974).
We should feed the model with more recent data.

In Figure 7.38, starting from around year 2050, the plot of the variable pop is under the plot of the variable iopc (cf. right subfigure attached below). Instead, following the book figure, pop should be over iopc (cf. left subfigure attached below).

This issue persists after considering different solvers.

In the following, we provide a minimum reproducing example.

using WorldDynamics
using DifferentialEquations

include("../scenarios/world3_historicalrun.jl")


system = world3_historicalrun()
sol = WorldDynamics.solve(system, (1900, 2100))


@named pop = World3.Pop15.population()
@named br = World3.Pop15.birth_rate()
@named dr = World3.Pop15.death_rate()
@named is = World3.Capital.industrial_subsector()
@named ss = World3.Capital.service_subsector()
@named ld = World3.Agriculture.land_development()
@named nr = World3.NonRenewable.non_renewable()
@named pp = World3.Pollution.persistent_pollution()
@named se = World3.SupplementaryEquations.supplementary_equations()
@named ai = World3.Agriculture.agricultural_inputs()
@named lfd = World3.Agriculture.land_fertility_degradation()

@variables t

fig_7_38_variables = [
    (nr.nrfr,  0, 1,    "nrfr"),
    (is.iopc,  0, 1000, "iopc"),
    (ld.fpc,   0, 1000, "fpc"),
    (pop.pop,  0, 16e9, "pop"),
    (pp.ppolx, 0, 32,   "ppolx"),
    (br.cbr,   0, 50,   "cbr"),
    (dr.cdr,   0, 50,   "cdr"),
]

cap_tables_7_38 = World3.Capital.gettables()
cap_tables_7_38[:isopc2] = (60, 450, 960, 1500, 1830, 2175, 2475, 2700, 3000)

agr_tables_7_38 = World3.Agriculture.gettables()
agr_tables_7_38[:ifpc2] = (345, 720, 1035, 1275, 1455, 1605, 1725, 1815, 1875)

pop_parameters_7_38 = World3.Pop15.getparameters()
pop_parameters_7_38[:zpgt] = 1975

cap_parameters_7_38 = World3.Capital.getparameters()
cap_parameters_7_38[:iet] = 1990
cap_parameters_7_38[:iopcd] = 320
cap_parameters_7_38[:alic2] = 21
cap_parameters_7_38[:alsc2] = 30

nr_parameters_7_38 = World3.NonRenewable.getparameters()
nr_parameters_7_38[:nruf2] = 0.125

pol_parameters_7_38 = World3.Pollution.getparameters()
pol_parameters_7_38[:ppgf21] = 0.25

system = world3_historicalrun(pop_params=pop_parameters_7_38, capital_params=cap_parameters_7_38, nonrenewable_params=nr_parameters_7_38, pollution_params=pol_parameters_7_38, capital_tables=cap_tables_7_38, agriculture_tables=agr_tables_7_38)
sol_7_38 = WorldDynamics.solve(system, (1900, 2100), solver = AutoVern9(Rodas5()))

plotvariables(sol_7_38, (t, 1900, 2100), fig_7_38_variables, name="Fig. 7.38", showlegend=true, showaxis=true, colored=true)

For sanity check, we also provide herafter the DYNAMO code snippet from the book chapter.

At the moment, the atol value in the switch function is arbitrary.
We should compute a proper value.

The initialization of falm and lfrt are incorrect in Figures 4.69d and 4.70d. This seems to be related to the usage of the smooth function when computing derivatives.

Following the example of other projects and packages, we can add the list of the project maintainers to the README.

This is a minor problem, I am creating this issue just to remember to improve the aesthetics of the plots.

So far, we don't test anything. The biggest challenge is to define how we can test if a system is working correctly, given that the expected result is a figure.

Currently, the logos are shown in a disorganized way.
Following some examples I saw in other repositories, I can put them in a Markdown table.

For solver that has keyword arguments associated with it, the user should be able to easily pass them throught WorldDynamics.solve function.
For example, in Issue #79, @natema created the system, simplified it, created the problem, and solved it using ModelingToolkit explicitly. This process should be simplified.

To generate links from one documentation page to an anchor in another documentation page .

I'm adding the usual badges to the README to link the documentation.

For some figures of Chapter 2 and Chapter 7, there is an incorrect bump during the initialization. This bump concerns most of the time the variable cbr.

While some major fluctuations of the cbr variable have been fixed by changing the solver (cf. #60), such initial bump resist and seems unrelated to the choice of the solver.

Examples :

Example affecting other variables in addition to the cbr variable :

I registered the repo on zenodo.org and we'll be able to add a badge to the README and reference how to cite it as soon as we have a first release.

The functions reproducing the plots internally solve the system and then only return the plot.
It would be handy to also have the solution available.
I'm not advancing any strong opinion on how this should happen (e.g. returning the solution together with the plot), but it's currently a waste to solve the system and then drop the solution away after plotting.
In trying to diagnose #67 for example, I would have liked to be easily able to inspect the solution on which the plot is based.

I have at least three good reasons to do so, listed below.

1. Performance improvement: functions are guaranteed to be type stable.
2. Clarity: users can clearly see which type of arguments they should pass and receive.
3. Avoid future bugs: since we are still making several changes in the code, is not uncommon to call functions passing wrong types.

We should list the authors in alphabetical order.

It is possible to do this without changing the package code, but it's not trivial. Currently, we have to look for the equation we want to change in the vector of system equations manually.

In MTK when we @register a function we tell Symbolics.jl not to recurse on it, but it still seem that we need to provide the derivative of the registered function.
So I'm opening this issue not to forget to have a closer look at this, at some point.

It would be good to add a Getting Started section to the README, showing how to add the package and then reproduce one figure, before linking to the tutorial.

At the moment, some standard figures from the historical runs of the book Dynamics of Growth in a Finite World (1974) are implemented.
We should implement all of the figures and offer a side-by-side comparison to converge towards asserting that the current code is indeed in pair with the original.

In the tutorial we should add using ModelingToolkit.

Currently, the tutorial shows how to reproduce a scenario by using the source files of the package.
We should allow users to reproduce a scenario without having to access the source files.

Although the standard figures are currently implemented and at least very close to the original ones, implementing all of them (including the parameter variations) would increase the confidence that the current code is indeed in pair with the original.

In the definition of pop1 and pop2 in odesystem of agricultural sector, the exponential function currently receives the exppop parameter which can be found in the parameters file.
Right now, the value of exppop is 0.011 and the code runs fine.
However, the code breaks if exppopv is defined as 0.012 (the actual value in World3 code).

I suggest to prepend to each fig_* function with a docstring so that we would get in the documentation a big list of all those function, with a brief line recalling what the figure is about in the book.

Related to discussion #83 .

We should be able to add equations. Currently, the only figures not implemented are those from Chapter 7 that require inserting new equations (namely, 23, 24, 26, 27, 30, 32, 39, and 41).

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

In the book "Limits to Growth: The 30 Year Update", the authors made a big revision of the World3 model, including new equations, variables, and parameters.
As such, it is a natural follow-up to add the revision into WorldDynamics.jl.

Although we are already able to do so without touching the package code, there is no example showing how to do this.

Instantiating the current version and running it gives us the following error:

ERROR: LoadError: UndefVarError: [VARIABLE] not defined

@JuliaRegistrator register

Related to issue #59, we have a problem with figures 2.97-100.
The changes of zpgt and fcest are inconsistent, which leads to wrong plots.

So far the tests are restricted to check if scenarios and figures codes run.

It might be useful to enhance the documentation of some parameters (with the book as a reference), that would act as a sanity check and help ensure that some inconsistencies (e.g., Issue #60) are not caused by a parameter discrepancy between the code and the book.

For example :

• In Chapter 2, instead of 1900, 1970, and 2100, the authors use 0, 70, and 200, respectively. Thus, the unit of lt is years.
• By looking at the source code in Chapter 7 (cf. screenshot below), we see that ppgf2 is a constant (= 1). Instead, if we look at Chapter 6, ppgf2 is a switch. This had led to an instability perceived in Figure 7.20, which was fixed (thanks to @paulobruno) by setting that the change on variable ppgf2 in the book corresponds to ppgf21 in the code.

pol_parameters_7_20 = World3.Pollution.getparameters()
pol_parameters_7_20[:ppgf21] = 0.1