Asking for advice on development direction

[hyunjimoon]

Not knowing the current tool ecosystem of System dynamics community, I think I am wasting some time trying to figure out how to best deliver the tools. May I ask your advice on that aspect? @tomfid @tseyanglim

Being closer vs Being further to Vensim is the first decision I need to make. To be honest, blackboxness of Vensim and diversity of tools used in SD (anylogic, Stella) made me inclined to choose the latter at first. But after meeting Tom (with his support and overlapping goals in introducing Bayes Stats.) and getting to know vensim's advanced functions (sythesim, macro) I am leaning towards the first (still confused, though). There is no clear answer, but I am sure about one thing.

I feel this project is a long fight with a Goliath. Everytime I identify the stark gap between SD and Stats community. e.g. multiplicative vs additive noise, the number of parameters (300 vs 3), the existence of driving data, I am becoming certain this cannot be done alone but need some experienced collaborators. So may I ask for some help in plans for future development-side? I am not asking for a direct answer, but as you are in the frontier of this domain (code, Bayes + SD) your opinion matters greatly to me.

For the code I've been developing (which I walked through with Tom yesterday), currently it is on pysd repo, but as pysd and Stan use separate builder (Stan has its own block-based coding system which helps Bayesian modeling), the only part resorting to pysd is this submodule which builds Abstract Syntax Tree (AST). (many thanks to pysd's AST and its excellent documentation). However, considering different goals of pysd and BayesSD, James and I agreed pysd might not be the best home for hosting inferece codes. Functions like multiview-support and subscript will be needed for future Bayes+SD, but these functions, I can use pysd as a package. So second decision might be Being closer vs Being futher to PySD. The first (I discussed with James) is combining the pysd and BayesSD interface (data generation with Stan). The unified format makes comparing the two (from parameter to observable and from observable to parameter) easier. The second is to use PySD as a library. I think the first decision above (being close/further to vensim) would affect this decision. So to be honest, how vensim is seeing BayesSD (which @tomfid could share best) would be the most helpful information.

Regarding my Bayesian check code above, please don't run it yet as I haven't tested the full automation. However, this python module is the main whose concept is described here.

[tseyanglim]

At the risk of jumping into a discussion without knowing full context - and with the caveat that I've not worked with Stan myself -

@tomfid and Hazhir and I had an email thread recently discussing how to do basic Bayesian inference tasks in Vensim. Short version is that a lot of basic stuff can be done in Vensim, it's just quite fiddly. But what Vensim (and really Vengine) has going for it is the MCMC engine, which is fast, and allows working with larger, more complex models (e.g. the 300 vs. 3 parameters) without absurd computational requirements. My impression is that its speed and efficiency is hard to match without access to a HPC cluster. Computational constraints are a huge concern for Bayesian inference with any substantial model so this is a major advantage.

That's why for me personally my approach has been to retain Vensim/Vengine as the core analytical engine while trying to make it easier to do all upstream and downstream tasks outside Vensim, drawing on the wealth of other existing tools out there. Obviously it would be ideal if Vensim was better or more easily integrated with those ecosystems. But for now, the advantage of computing power is still (for me) worth the lack of straightforward integration.

[hyunjimoon]

This is tremendously helpful, thanks @tseyanglim. While waiting for Tom's comment, if I may share further context:

• Tom and I had discussed the possibility of introducing Stan library on autodiff(math), HMC, and convergence diagnostics codes in C++.

• I think adding SOPS/BATS/SDA on top of venim or its python extension (which I communicated with Tom) would have great value. I am making each libraries' introduction here as one diagram which may give you the idea of Stan as well.

[tomfid]

Here's a possible approach:

• Do this end to end in Python & Stan
• Document all the mappings to conventional SD
• Repeat in Vensim, or Vensim wrapped in venpy with Stan, and observe/document an differences, opportunities, etc.

[hyunjimoon]

I will be doing the packaging job early next week to deliver Stan version of vensim-python connection to seminar students. During this process, I wish to discuss the best location to host this service.

Agenda for week-unit

@tseyanglim shared this notebook with the following functions:

notebook showcases vst version 0.1.0 with updated functions and syntax to demonstrate a basic but complete model analysis workflow, using a simple toy model:
- Partial model calibrations using variant model files (generation of these is manual for now but could easily be automated)
- Full model calibration, building on the results of previous partial calibrations
- Calibration with MCMC to derive a sample of posterior parameter estimates
- Semi-automated generation of additional files to run projections
- Basic projections + parametric sensitivity (for uncertainty intervals) based on calibrated parameter estimates
- Projections for policy analysis + sensitivity analysis

To avoid the worst case of not introducing Stan connection, this week I am aiming to replicate the above with Stan-version first to
i) compare the value of state and estimated parameter sample
ii) get review from Hazhir, Tom, TY, Jason

Then I will move on to inventory management model. Meanwhile, @tseyanglim could help me by testing this inventory model in VST? Your VST result can make checking validity much easier in python; VST is only library that can do this. @hazhirr thinks pushing to this complexity is needed (which I agree)!

Agenda for month-unit

Restarting from Tom's reply above (proposal accepted!, thanks)

Do this end to end in Python & Stan
Document all the mappings to conventional SD
Repeat in Vensim, or Vensim wrapped in venpy with Stan, and observe/document an differences, opportunities, etc.

the following is my objective and constraints:

Binary-ified objective

• whether or not this service will be sustainable i.e. implementing reinforcement loop between stock wisdom and its inflow learning
• format: .vpd, model block of _data2draws.stan

Strong constraints

• must connect Stan with Vensim by including:
• Stan's convergence diagnostic, optimization engine (HMC, ADVI), penalty functions for diverse distribution
• Vensim's established data analysis framework e.g. dynamic hypothesis updating
• C++ level may be needed
• @tomfid could you help me by listing some service of Vensim that I can publicize to Stan users interested in dynamic model?
• format: parameter range in .voc, parameter block in _data2draws.stan

Weak constraints

Pooling information from four agents (MIT SD group, Stan team, SD community, Stan community) helped increase both stock level (literature recommendation) and flow rate (mail exchange). Delay was shorter than expected: witnessed co-evolve in less than one month in each.

• I think talking to pharmacokinetics/pharmacodynamics domain will be most helpful in the future e.g. one actively library: https://metrumresearchgroup.github.io/Torsten/
• format: regularization term to payoff in .vpd, as prior in model block of _data2draws.stan

[tseyanglim]

What exactly would you like me to test? estimating parameters for the inventory model w/ Vensim+VST MCMC? if so, will need data, payoff, and optimization control files

[hyunjimoon]

Thank you very much for quick reply and great question!
Could we compare the vector of production_rate_stocked and production_start_rate_stocked with the following setting?

This python file walks through model setting, generation, and estimation.

1. time_step:

n_t = 20

2. classification of parameter

setting_assumption = {
    "est_param_scalar" : ("inventory_adjustment_time", "minimum_order_processing_time"),
    "ass_param_scalar" : ("inventory_coverage", "manufacturing_cycle_time", "safety_stock_coverage", "time_to_average_order_rate", "wip_adjustment_time"),
    "target_simulated_vector" : ( "production_rate_stocked", "production_start_rate_stocked"),
    "driving_data_vector" : ("customer_order_rate", "process_noise_std_norm_data", "production_start_rate_m_noise_trun_norm_data", "production_rate_m_noise_trun_norm_data"),
    "model_name": "mngInven",
    "integration_times": list(range(1, n_t + 1)),
    "initial_time": 0.0
}

3. differential equation structure
   importing differential equation structures from vensim file

vf = VensimFile('vensim_models/mngChain/InventoryManagementWeek7woPN.mdl')

4. numeric value of driving data

numeric_assumption = {
    "n_t": n_t,
    "customer_order_rate": pd.read_csv('data/example_retail_sales.csv').iloc[:n_t, 1].values,
    "process_noise_std_norm_data": np.random.normal(0,1, size=n_t),
    "production_start_rate_m_noise_trun_norm_data": truncnorm.rvs(0, 2, size=n_t),
    "production_rate_m_noise_trun_norm_data": truncnorm.rvs(0, 2, size=n_t),
}

Important thing is the inclusion of random variable in the form of driving data. i.e. process_noise_std_norm_data production_start_rate_m_noise_trun_norm_data, production_rate_m_noise_trun_norm_data from above.

Two different implementations for random variable were identified with @jandraor: First is Vensim's implementation of generating random variable as along time. Second is my implementation of using single parameter value for each generator and treat aggregated synthetic data from M different generators as if they are from one one generator.

Vensim:

t = 1
- sample theta from N(2,1): theta_1 = 2
- sample y_1 from neg_binom(theta = theta_1, phi): y = 4
t = 2
- sample theta from N(2,1): theta_2 = 1
- sample y_2 from neg_binom(theta = theta_2, phi): y = 1
t = 3
- sample theta from N(2,1): theta_3 = 3
- sample y_3 from neg_binom(theta = theta_3, phi): y = 5

y_(11,12,13) = 4, 1, 5
Repeat the above three times: 
- sample $y_(21,22,23)$ = ...
- sample $y_(31,32,33)$ = ...

Current version of Stan connection

sample all three theta ~ N(2, 1): theta_1 = 2, theta_2 = 1, theta_3 = 3
- sample $y_(11,12,13)$ from neg_binom(theta_1 = 2, phi): y = 4,3,6
- sample $y_(21,22,23)$ from neg_binom(theta_1 = 1, phi): y = 1,1,2
- sample $y_(31,32,33)$ from neg_binom(theta_1 = 3, phi): y = 2,3,4

Let alone which implementation is correct, I wish to make sure the difference between the two is small enough as implementing the change of coefficients in differential equations is much complicated.

[hyunjimoon]

@tomfid,
@tseyanglim may be joining the 3pm meeting; I wish to verify vensim (from VST) and stan outputs (from stanify).
For this, I prepared https://github.com/Data4DM/stanify repository; you can run the demo.py within which returns plots for both stock values and estimated parameter samples.

[hyunjimoon]

@tseyanglim, thanks for you help yesterday.
#58 is what @tomfid and I thought is needed at this point which I wish to give demo to 879 students next Friday: verification framework in three different methods (Vensim, VST, stanify).

Could you offer some help by sharing VST code and result for the following by preferably Tuesday morning? Sorry, I couldn’t afford preparing three demos alone. Results, in short, are:

1. 30 different parameter sets from prior used in stanify demo
2. outcome array (xarray or inference datatype is the format in the end which can passed on to library arviz — diagnostic plot heaven)
3. 30 * 1000 number of estimated sample

Lmk if anything is unclear ;)

[tseyanglim]

@hyunjimoon Could you clarify exactly what the different outputs / runs you want are?

1. Predator-prey model outputs run with 30 different parameter draws from the priors specified:

model.set_prior("alpha", "normal", 0.8, 0.08, lower = 0)
model.set_prior("beta", "normal", 0.05, 0.005, lower = 0)
model.set_prior("delta", "normal", 0.05, 0.005, lower = 0)
model.set_prior("gamma", "normal", 0.8, 0.08, lower = 0)

or do you just want the prior draws themselves?

2. I'm not sure which outcome array you mean here

3. I assume the '* 1000' part is 1000 MCMC samples, but it's still not clear to me why you'd have 30 different sets of this (corresponding to the 30 parameter sets drawn from the prior?)

(Alternatively, instead of the three results you mentioned, you could outline in plain language what is the analytical procedure that you want replicated, step by step)