From 37a8308790c5af0a09a844e4560ed7223546736e Mon Sep 17 00:00:00 2001
From: abaucher <achille.baucher@inria.fr>
Date: Tue, 8 Feb 2022 17:37:40 +0100
Subject: [PATCH] Better compare, various plot, removed dynamo_all, cleaner
pydynamo_code
---
pydynamo/__init__.py | 4 +--
pydynamo/core/plot_system.py | 9 +++++-
pydynamo/core/system.py | 8 ++---
pydynamo/world3/code_pydynamo_w3.py | 47 +++++------------------------
requirements.txt | 2 ++
5 files changed, 24 insertions(+), 46 deletions(-)
diff --git a/pydynamo/__init__.py b/pydynamo/__init__.py
index 5206dcdc..ea3e9401 100644
--- a/pydynamo/__init__.py
+++ b/pydynamo/__init__.py
@@ -26,9 +26,9 @@ def plot_w2(w2, title=''):
def get_scenario(number):
w3s = get_w3()
- changes = scenarios[number]['changes']
+ changes = scenarios[number - 1]['changes']
for c, v in changes.items():
setattr(w3s, c, v)
- return w3s
+ return w3s, scenarios[number - 1]['title']
diff --git a/pydynamo/core/plot_system.py b/pydynamo/core/plot_system.py
index c51a11f0..8786ecd6 100644
--- a/pydynamo/core/plot_system.py
+++ b/pydynamo/core/plot_system.py
@@ -62,8 +62,15 @@ def plot_tabhl(s, name):
plt.title(title)
def compare_systems(s1, s2, v_names, scales=None, rescale=False, *args, **kwargs):
+ # remove tables
+ v_names = [v for v in v_names if not isinstance(getattr(s1,v), list)]
if not scales and rescale:
- scales = {v: max(max(getattr(s1, v)), max(getattr(s2, v))) for v in v_names}
+ scales= {}
+ for v in v_names:
+ try:
+ scales[v] = max(max(getattr(s1, v)), max(getattr(s2, v)))
+ except TypeError:
+ scales[v] = max(getattr(s1, v), getattr(s2, v))
plot_system(s1, v_names, *args, **kwargs, linestyle='-', scales=scales)
plt.gca().set_prop_cycle(None)
plot_system(s2, v_names, *args, **kwargs, linestyle='--', scales=scales)
diff --git a/pydynamo/core/system.py b/pydynamo/core/system.py
index 984fb20e..2565b580 100644
--- a/pydynamo/core/system.py
+++ b/pydynamo/core/system.py
@@ -365,8 +365,9 @@ class System:
raise Exception(f"Error in setting constants {cst}, used in:\n"
f"{use_eqs}\n"
f"Constant {k.args[0]} is never set.") from None
-
- self.set_cst(cst, args)
+
+ if not cst in dir(self):
+ self.set_cst(cst, args)
def generate_var(self, var, N):
setattr(self, var, np.full(N, np.NAN))
@@ -537,14 +538,13 @@ class System:
f"Returned value: {value}"))
- def prepare(self, initial_time=0):
+ def prepare(self):
self.assert_cst_defined()
self.assert_init_defined()
self.assert_update_defined()
assert nx.is_directed_acyclic_graph(self.get_cst_graph()),\
("Cycle detected for constant equations",
nx.find_cycle(self.get_update_graph()))
- self.initial_time = initial_time
self.set_all_funs()
self.set_all_csts()
diff --git a/pydynamo/world3/code_pydynamo_w3.py b/pydynamo/world3/code_pydynamo_w3.py
index f19f3d01..405fbaee 100644
--- a/pydynamo/world3/code_pydynamo_w3.py
+++ b/pydynamo/world3/code_pydynamo_w3.py
@@ -72,12 +72,12 @@ dcfs.k = clip(2.0, dcfsn*frsn.k*sfsn.k, time.k, zpgt)
zpgt = 4000
dcfsn = 3.8
sfsn.k = tabhl(sfsnt, diopc.k, 0, 800, 200)
-sfsnt = [1.25, 0.94, 0.715, 0.59, 0.5] # ! before wrld3 it was [1.25, 1, 0.9, 0.8, 0.75]
+sfsnt = [1.25, 0.94, 0.715, 0.59, 0.5] # Changed from Vensim. Before it was [1.25, 1, 0.9, 0.8, 0.75]
diopc.k = dlinf3(iopc.k, sad)
sad = 20
frsn.k = tabhl(frsnt, fie.k, -0.2, 0.2, 0.1)
frsnt = [0.5, 0.6, 0.7, 0.85, 1]
-# frsn.i = 0.82 # ! removed with wrld3
+# frsn.i = 0.82 # Removed from Vensim
fie.k = (iopc.k - aiopc.k)/aiopc.k
aiopc.k = smooth(iopc.j, ieat)
ieat = 3
@@ -93,7 +93,6 @@ fsafct = [0, 0.005, 0.015, 0.025, 0.03, 0.035]
# Capital sector
## Industrial subsector
iopc.k = io.k/pop.k
-# ATTENTION: from next line a small difference is induced
io.k = ic.k*(1-fcaor.k)*cuf.k/icor.k
icor.k = clip(icor2.k, icor1, time.k, pyear)
icor1 = 3
@@ -101,7 +100,7 @@ icor2.k = icormrrct.k*icormlyt.k*icormpt.k
icormrrct.k = tabhl(icormrrctt, nruf.j,0, 1, 0.1) # industrial capital output ratio multiplier from resource conservation technology
icormrrctt = [3.75, 3.6, 3.47, 3.36, 3.25, 3.16, 3.1, 3.06, 3.02, 3.01, 3] # industrial capital output ratio multiplier from resource table
icormlyt.k = tabhl(icormlytt, lymt.j, 1, 2, 0.2) # industrial capital output ratio multiplier from land yield technology
-# TRUC BIZARE !!(1,0.8)-(2,2)],
+# Note: Weird stuff in Vensim file (1,0.8)-(2,2)],
icormlytt = [1, 1.05, 1.12, 1.25, 1.35, 1.5] # industrial capital output ratio multiplier table
icormpt.k = tabhl(icormptt, ppgf.j, 0, 1, 0.1) # industrial capital output ratio multiplier from pollution technology
icormptt = [1.25, 1.2, 1.15, 1.11, 1.08, 1.05, 1.03, 1.02, 1.01, 1, 1] # industrial capital output ratio multiplier from pollution table
@@ -163,7 +162,7 @@ lf.k = (p2.k + p3.k)*lfpf
lfpf = 0.75
luf.k = j.k/lf.k
lufdi.k = clip(luf.j, 1, time.k, initial_time) # helper for luf
-lufd.k = smooth(lufdi.k, lufdt) #! added from wrld3
+lufd.k = smooth(lufdi.k, lufdt) # Added from Vensim
lufdt = 2
cuf.k = tabhl(cuft, lufd.k, 1, 11, 2)
cuf.i = 1
@@ -242,7 +241,7 @@ llmy.k = clip(0.95**((time.k - llmytm)/oy)*llmy1.k + (1 - 0.95**((time.k - llmyt
llmy1.k = tabhl(llmy1t, ly.k/ilf, 0, 9, 1)
llmy1t = [1.2, 1, 0.63, 0.36, 0.16, 0.055, 0.04, 0.025, 0.015, 0.01]
llmy2.k = tabhl(llmy2t, ly.k/ilf, 0, 9, 1)
-llmy2t = [1.2, 1, 0.63, 0.36, 0.29, 0.26, 0.24, 0.22, 0.21, 0.2] # WRLD3+03 !
+llmy2t = [1.2, 1, 0.63, 0.36, 0.29, 0.26, 0.24, 0.22, 0.21, 0.2] # Note: added from Vensim
ler.k = al.k/all.k
uilpc.k = tabhl(uilpct, iopc.k, 0, 1600, 200)
uilpct = [0.005, 0.008, 0.015, 0.025, 0.04, 0.055, 0.07, 0.08, 0.09]
@@ -267,7 +266,7 @@ ilf = 600
lfrt.k = tabhl(lfrtt, falm.k, 0, 0.10, 0.02)
lfrtt = [20.0, 13.0, 8.0, 4.0, 2.0, 2.0]
-## Loop 6: dDiscontinuing land maintinance
+## Loop 6: Discontinuing land maintinance
falm.k = tabhl(falmt, pfr.k, 0, 4, 1)
falmt = [0.0, 0.04, 0.07, 0.09, 0.1]
fr.k = fpc.k/sfpc
@@ -304,7 +303,6 @@ fcaor2t = [1.0, 0.2, 0.1, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05]
ppgr.k = (ppgio.k + ppgao.k)*ppgf.k
ppgf.k = clip(ppgf2.k, ppgf1, time.k, pyear)
ppgf1 = 1
-# Added for 03
ppgf2.k = smooth(ppt.k, tdd)
tdd = 20 # technology development delay
ppt.k = ppt.j + dt*pptcr.k # persistent pollution technology
@@ -314,8 +312,7 @@ pptcr.k = clip(ppt.j*pptcm.j, 0, time.k, pyear) # persistent pollution technolog
pptcm.k = tabhl(pptcmt, 1 - ppolx.k/dppolx, -1, 0, 1)# persistent pollution technology change multiplier
pptcmt = [0, 0] # persistent pollution technology change mult table
dppolx = 1.2 # desired persistent pollution index
-# industrial capital output ratio multiplier from persistent pollution technology
-# is not implemented in Vensim while on the schemas
+# Note: industrial capital output ratio multiplier from persistent pollution technology is not implemented in Vensim while on the schemas
ppgio.k = pcrum.k * pop.k * frpm*imef*imti
frpm = 0.02
imef = 0.1
@@ -336,10 +333,6 @@ ahlmt = [1.0, 11.0, 21.0, 31.0, 41.0]
ahl.k = ahl70*ahlm.k
ahl70 = 1.5
-# Supplementary equations sector
-# foa.k = 0.22*f.k/(0.22*f.k + so.k + io.k)
-# foi.k = io.k/(0.22*f.k + so.k + io.k)
-# fos.k = so.k/(0.22*f.k + so.k + io.k)
# World3 03 supplementary equations sector
cio.k = io.k*fioacc.k # consumed industrial output
ciopc.k = cio.k/pop.k # consumed industrial output per capita
@@ -377,31 +370,7 @@ ulgha.k = uil.k/hpgha # "Urban Land (GHA)"
algha.k = al.k/hpgha # "Arable Land in Gigahectares (GHA)"
# Control card for simulation
-# Some have been removed because not implemented this way in pydynamo
-# SPEC dt = 0.5
-# SPEC length = 2100
pyear = 1975
-# N time = 1900
initial_time = 1900
-# A pltper.k = step(plp, plit)
-# C plp = 5
-# C plit = 1900
-# C prp = 0
-# A prtper.k = step(prp, prit) + step(-prp, prtt)
-# C prit = 1900
-# C prtt = 2100
-# Take care to keep the last line !
-# Parameter and structural changes for limits to growth
-
-
-# Additional lines for limits to growth code
-# Capital sector
-## Industrial subsector
-# icir.k = clip(icir2.k, io.k*fioai.k, time.k, icet)
-# helpicir2.k = min(icdr.j, io.k* fioai.k) # Icir2 value if certain strategy is chosen
-# icir2.k = clip(helpicir2.k, io.k* fioai.k, diopc.k-diop.k, 0) # Icir value after icet
-# icet=4000 # Industrial capital economy time
-# diop.k = sample(iopc.k, dist.k, 0) # Desired industrial ouptput p
-# dist.k = step(4000, disi+1905)+disi # Desired industrial stop time
-# disi=4000 # Desired industrial stop interval
+# Take care to keep the last line :)
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
index 3b06f9ad..aa968044 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -27,3 +27,5 @@ pyvis==0.1.9
six==1.16.0
traitlets==5.1.1
wcwidth==0.2.5
+pandas==1.4.0
+jupyter-lab==3.2.9
--
GitLab