Cleanup data structures

[hdrake]

Before tackling documentation #2 and unit tests #3, the first priority is to cleanup the data structures (see branch https://github.com/hdrake/ClimateMARGO.jl/tree/data-structure-cleanup).

There are a few features which are particularly undesirable about the current data structures:

1. We hardly use any functions in the JuMP optimization code src/optimization/jl, instead choosing to write out the objective function in its entirety, which is excessively long and susceptible to bugs. My guess is that it is possible to write the objective condition in one line that looks something like (in pseudo-JuMP syntax)

for i in N; @constraint{ T(M,R,G,A) < Tstar }; end

where M, R, G, A are the control variables.
2. The few functions that are currently used in the JuMP optimization code have explicitly re-defined within the optimization code, probably for legacy reasons. Ideally, the diagnostic functions used for plotting / analysis should be re-used for the optimization to reduce the likelihood of discrepancies between the two versions of the code (which indeed caused breaking bugs in early versions of the implementation).
3. There are a lot of redundant functions that should be collapsed into a single function with keyword arguments and/or different methods.
4. We should rethink the ClimateModel, Physics, Controls, and Economics structs. They are fairly unwieldy because I did not really understand Julia constructors when I created them half a year ago. It might be useful to look at how more mature julia models like https://github.com/CliMA/Oceananigans.jl or https://github.com/CliMA/ClimateMachine.jl structure their submodules for inspiration.

[fonsp]

Hi! About 4: I think the current structure is easy to understand and work with - I do have some minor notes:

dt

I spent some time debugging why changing model.dt led to strange results and no change in runtime - I found that I also needed to set t. This is fine (I know it now 🙃) but you could fix it by taking dt out of the ClimeateModel struct and defining a function instead:

dt(model::ClimateModel) = model.domain[2] - model.domain[1]

A more Julian way would be to have domain be of the type AbstractRange (like 0:2:10). This would reflect the property that the model only works with constant dt. You would then define

dt(model::ClimateModel) = step(model.domain)

This does require some rewrites in the rest of the code and notebooks - so do with this what you like!

model defaults

You could write some of the model defaults directly in the constructor - this helps to "keep the scope clean" - it is not clear from the names of the variables in defaults.jl that they are default values. For example, discounting accidentally uses t (default) instead of model.domain (actual).

[hdrake]

Hi @fonsp, thanks for the suggestion. This was definitely just me be sloppy early on and now thinking of a better fix, such as what you propose above. I am busy for the next few hours but will try to make the necessary changes today.

A similar problem is the specification of the feedback parameter and climate sensitivity. Currently, the climate sensitivity is treated as a "diagnostic" quantity that is evaluated upon construction of the Physics() struct. However, since this is a mutable struct its value can be changed after construction, which is something that we do when "stepping forward" and re-elvaluating the parameter choices. The problem is that right now the only way to do this is to manually specify both the indepedent parameters (e.g. the feedback parameter) as well as the dependent parameter diagnostics (e.g. the climate sensitivity).

mutable struct Physics
    CO₂_init::Float64
    δT_init::Float64
    a::Float64
    B::Float64
    Cd::Float64
    κ::Float64
    r::Float64
    
    ECS::Float64
    τd::Float64
    function Physics(CO₂_init, δT_init, a, B, Cd, κ, r)
        FCO₂_2x = a*log(2) # Forcing due to doubling CO2 (Geoffrey 2013)
        sec_per_year = 60. * 60. * 24. * 365.25
        
        ECS = (FCO₂_2x*sec_per_year)/B # [degC]
        τd = (Cd/B) * (B+κ)/κ # [yr]
        return new(CO₂_init, δT_init, a, B, Cd, κ, r, ECS, τd)
    end
end

As in the case of t and dt, only one independent parameter should be required and the others should be automatically re-diagnosed whenever it changes.

[fonsp]

Right! You could make those dependent parameters (ECS and τd) functions of a Physics object, the parameter that is being changed dynamically is the independent B right?

Don't feel rushed to fix it today - I can work around these minor points. 👍

[hdrake]

About the model default values, eventually I want to move away from hard-coding this anywhere, whether in default.jl or directly in the constructor. Instead, I think it would be better to have the defaults in a .yml or .csv file. This would also be nice because it would make for easy sharing of parameter configurations.

For example, I imagine you could email me your run-time parameter fonsp_favorite.yml YAML file that has a structure like:

--- # Physical parameters
  a: 4.97
  B: 1.17
...

and I could just run:

ClimateModel(parameter_file = "fonsp_favorite.yml")

or something similar.

I think this would be really complimentary to have the defaults directly in the constructor as you suggest, because then we can set it up such that if any parameters are left out of the YAML / CSV file, then it can just take the hard-coded default value.

[fonsp]

Oh that's a cool idea! I'd recommend JSON over yaml/csv - since JSON can contain nested data structures:

The package JSON2.jl can directly read JSON files/strings into Julia types:

[hdrake]

That is really convenient! It would be great to have human-readable and easily modifiable input files that uniquely determine a MARGO configuration!

[hdrake]

I've started Julia-fying the package structure and submodel data structures in this branch https://github.com/hdrake/ClimateMARGO.jl/tree/updating-structure, which is still very much a work in progress.

Interpretable and convenient functions for diagnostics
My design philosophy has been that each key diagnostic, like temperature, should be a function with two key methods:

1. an interpretable method that "feels" more like a math equation and can be used outside of the context of the ClimateMARGO submodel structs

T(T0, F, Cd, κ, B, t, dt; A=0.) = sqrt.(1. .- A) .* (
    T0 .+
    T_fast(F, κ, B) .+
    T_slow(F, Cd, κ, B, t, dt)
)

and
2) a convenience method whose only argument is a ClimateModel instance, where the default behavior is the "baseline" and then each control can optionally be turned from false to true to turn it on

T(model::ClimateModel; M=false, R=false, G=false, A=false) = T(
    model.physics.T0,
    F(model, M=M, R=R, G=G),
    model.physics.Cd,
    model.physics.κ,
    model.physics.B,
    t(model),
    model.domain.dt,
    A=model.controls.adapt .* (1. .- .~future_mask(model) * ~A)
)

Re-tooling dependent parameters as diagnostic functions
As suggested, the t array for the time dimension is now a function, rather than a hard-coded variable in a struct. Similarly, the climate sensitivity can be diagnosed from the independent physical parameters.

ECS(a, B) = F2x(a)/B
ECS(model) = ECS(model.physics.a, model.physics.B)

@fonsp, what do you think? If you like it, I'll continue re-implementing all of the other diagnostics. I think it makes the code much more readable, user-friendly, and bug-repellent.

[hdrake]

A lot of this has now been implemented in #19, but it is still a work-in-progress. Will close this Issue once that is merged.