MCS - Assigning random variables to indexed variable

daler6 · November 1, 2020, 3:06am

Hi,

I have been going back over Tutorial 5 as I have switched IDEs from Atom to VSCode. To simplify things, I am using the One-Region Model from Tutorial 4 rather than the Multi-Region Model. To further simplify things, I am just looking at two parameters - share and s. The latter is indexed by time, so I have been trying to just change the values for the time periods between 2015 and 2050. I did try to use a line like:

           s[2015:5:2020] = Uniform(0.2, 0.3)

This seems to work fine at first, but when I try to run the following:

si = run(sd, m, 10; trials_output_filename = “./Results/trialdata.csv”, results_output_dir="./Results")

I get this error “ArgumentError: indexed assignment with a single value to many locations is not supported; perhaps use broadcasting .= instead?”

So, I went back and modified the earlier line to read:

         s[2015:5:2020] .= Uniform(0.2, 0.3)

but, then I get this error: “LoadError: LoadError: Unrecognized expression ‘s[2015:5:2050] .= Uniform(0.2, 0.3)’ in @defsim”

I have tried lots of other things, but to no avail. I finally ended up with the following code to define the Monte Carlo Simulation, which worked with my si = run(. . . ) function.

sd = @defsim begin

# assign RVs to model Parameters
share = Uniform(0.2, 0.8)
s[2015] = Uniform(0.2, 0.3)
s[2020] = Uniform(0.2, 0.3)
s[2025] = Uniform(0.2, 0.3)
s[2030] = Uniform(0.2, 0.3)
s[2035] = Uniform(0.2, 0.3)
s[2040] = Uniform(0.2, 0.3)
s[2045] = Uniform(0.2, 0.3)
s[2050] = Uniform(0.2, 0.3)

# Indicate which variables to save for each model run.
# The syntax is: component_name.variable_name
save(grosseconomy.share, grosseconomy.s, grosseconomy.K, grosseconomy.YGROSS,
     emissions.E)

end

Surely, there must be a more elegant way to do this.

In addition, this now randomly assigns different values for s for each time step from 2015 to 2050 in each trial, rather than having the same values for s for all time steps from 2015 to 2050 within each trial, but different ones from trial to trial.
I tried many things, e.g.,

s[2015:5:2050] = s[2015:5:2050] = ones(8) .* Uniform(0.2, 0.3)

but to no avail.

p.s. I also played with the *= function like in the line from Tutorial 5. My specific line was:

s[2020:5:2050] *= Uniform(0.9, 1.1)

When I did this, I noted that instead of multiplying by the default value for s[2020:5:2050], which is 0.22, in each trial, it was multiplying by the value of s[2020:5:2050] from the previous trial. So, across the trials, instead of s[2020:5:2050] being a uniform distribution around 0.22, it took on the characteristic of a random walk starting at 0.22.

lrennels · November 2, 2020, 6:24pm

Hi Dale, thanks for the note. I will try to do some clarification below and then make sure I understand what functionality is missing so that we can fill the need. Some of the syntax you indicate trying:

s[2015:5:2020] = Uniform(0.2, 0.3)

Since this is all inside a Julia macro, this isn’t going to cause any issues until runtime, when Julia tries to run the expressions the macro creates and hits the broadcast error you encountered. If this was straight Julia instead of syntax within a macro, then the below would be the correct fix …

s[2015:5:2020] .= Uniform(0.2, 0.3)

As stated above this would in theory fix the problem, but we are actually taking these expressions and converting them into Julia, so you are hitting an error because we have not implemented understanding the .= operator within our macro. I’ll come back to this but for now this operator is not supported.

s[2020:5:2050] *= Uniform(0.9, 1.1)

As you point out, the *= syntax is not intended to work like .= but instead like it’s syntactically similar += that you might have seen in other languages. It operates as a random walk as you pointed out. That is explained in How-To Guide 3 as follows, but should be clarified in the text of the tutorial.

" param *= RV replaces the values in the parameter with the product of the original value and the value of the RV for the current trial "

Am I correct in saying that what you are hoping for is that you can (elegantly if possible ) set s at each time period 2015:5:2050 to the same value, pulled from a distribution, in this case Uniform(0.2, 0.3)? If so, I actually do think that your original attempt (#1) should work, as indicated by the ** Apply RVs to model parameters: Assigning to array slices** section in our tutorials, so will try to figure out now if there’s a bug making that fail.

lrennels · November 2, 2020, 7:44pm

Hi @daler6, so as it turns out you have found a bug, thank you! We have been struggling to deal with how our custom arrays handle broadcasting, it’s more complicated under the hood than expected, and this is just another case of that same issue. I will work on this (Issue is here: https://github.com/mimiframework/Mimi.jl/issues/779 with connected issue here: https://github.com/mimiframework/Mimi.jl/issues/765) and keep you in the loop, but for now the (not elegant) way to do what you want is this:

sd = @defsim begin

# assign RVs to model Parameters
share = Uniform(0.2, 0.8)
rv(s_rv) = Uniform(0.2, 0.3)

s[2015] = s_rv
s[2020] = s_rv
s[2025] = s_rv
s[2030] = s_rv
s[2035] = s_rv
s[2040] = s_rv
s[2045] = s_rv
s[2050] = s_rv

# Indicate which variables to save for each model run.
# The syntax is: component_name.variable_name
save(grosseconomy.share, grosseconomy.s, grosseconomy.K, grosseconomy.YGROSS,
     emissions.E)

end

daler6 · November 5, 2020, 5:53am

Lisa, thanks for looking into this. Two additional comments clarification in the Mimi Documentation:

You point to the text "param *= RV replaces the values in the parameter with the product of the original value and the value of the RV for the current trial " from the How-To-Guide 3. Since, as you agree, it operates as a random walk, it seems that the link should not say “product of the original value”, but rather “product of the value in the previous trial”? This, of course, leaves the question of whether there is a way to have it actually use the original value, which I think is more natural, i.e. to select a value in each trial that fall between 0.9 and 1.1 times the default value.
By the way, I patterned the line “s[2020:5:2050] *= Uniform(0.9, 1.1)” after the line “sigma[2020:5:2050, (Region2, Region3)] *= Uniform(0.8, 1.2)” in tutorial 5. For each trial, I believe this assigns the same value for sigma for all the years 2020, 2025, . . . , 2050 in Regions 2 and 3. It might be good to make this clear. Some people, e.g., me might accidentally interpret this line as implying that it chooses a separate value for each year and each region.

Anyway, I look forward to speaking tomorrow.

lrennels · November 10, 2020, 11:10pm

@daler6 I just wanted to follow up here and say I was incorrect, and the documentation was correct, on the *- and += operators. These will alter your original value, so you are getting variation around the original value. I did a small test for this that confirmed, by running the following and then looking at the trials values and the output model values. Let me know if you have any questions, since it seemed your results were indicating a random walk.

include("/Users/lisarennels/.julia/dev/Mimi/test/mcs/test-model-2/multi-region-model.jl")
using .MyModel
m = construct_MyModel()

N = 10

sd = @defsim begin

    share += Uniform(-0.3, 0.3)
    save(grosseconomy.share)

end

output_dir = joinpath(mypath, "sim")

# Run trials 
si = run(sd, m, N; trials_output_filename = joinpath(output_dir, "trialdata.csv"), results_output_dir=output_dir)

lrennels · December 16, 2020, 11:46pm

Hi @daler6 we have reinstated the functionality so that we broadcast under the hood when necessary (https://github.com/mimiframework/Mimi.jl/issues/779), as soon as we get this branch merged (https://github.com/mimiframework/Mimi.jl/pull/780) and a new version of Mimi tagged, you should now be a able to use your original syntax. Let me know if you have any questions, this will be done by the end of the week.