Merge branch 'master' into myb/dd

YingboMa · YingboMa · commit 8c3e05378919 · 2024-10-04T17:18:17.000-04:00
diff --git a/docs/src/basics/MTKLanguage.md b/docs/src/basics/MTKLanguage.md
@@ -63,14 +63,13 @@ end
     @structural_parameters begin
         f = sin
         N = 2
-        M = 3
     end
     begin
         v_var = 1.0
     end
     @variables begin
         v(t) = v_var
-        v_array(t)[1:N, 1:M]
+        v_array(t)[1:2, 1:3]
         v_for_defaults(t)
     end
     @extend ModelB(; p1)
@@ -311,10 +310,10 @@ end
   - `:defaults`: Dictionary of variables and default values specified in the `@defaults`.
   - `:extend`: The list of extended unknowns, name given to the base system, and name of the base system.
   - `:structural_parameters`: Dictionary of structural parameters mapped to their metadata.
-  - `:parameters`: Dictionary of symbolic parameters mapped to their metadata. Metadata of
-    the parameter arrays is, for now, omitted.
-  - `:variables`: Dictionary of symbolic variables mapped to their metadata. Metadata of
-    the variable arrays is, for now, omitted.
+  - `:parameters`: Dictionary of symbolic parameters mapped to their metadata. For
+    parameter arrays, length is added to the metadata as `:size`.
+  - `:variables`: Dictionary of symbolic variables mapped to their metadata. For
+    variable arrays, length is added to the metadata as `:size`.
   - `:kwargs`: Dictionary of keyword arguments mapped to their metadata.
   - `:independent_variable`: Independent variable, which is added while generating the Model.
   - `:equations`: List of equations (represented as strings).
@@ -325,10 +324,10 @@ For example, the structure of `ModelC` is:
 julia> ModelC.structure
 Dict{Symbol, Any} with 10 entries:
   :components            => Any[Union{Expr, Symbol}[:model_a, :ModelA], Union{Expr, Symbol}[:model_array_a, :ModelA, :(1:N)], Union{Expr, Symbol}[:model_array_b, :ModelA, :(1:N)]]
-  :variables             => Dict{Symbol, Dict{Symbol, Any}}(:v=>Dict(:default=>:v_var, :type=>Real), :v_for_defaults=>Dict(:type=>Real))
+  :variables             => Dict{Symbol, Dict{Symbol, Any}}(:v=>Dict(:default=>:v_var, :type=>Real), :v_array=>Dict(:type=>Real, :size=>(2, 3)), :v_for_defaults=>Dict(:type=>Real))
   :icon                  => URI("https://github.com/SciML/SciMLDocs/blob/main/docs/src/assets/logo.png")
-  :kwargs                => Dict{Symbol, Dict}(:f => Dict(:value => :sin), :N => Dict(:value => 2), :M => Dict(:value => 3), :v => Dict{Symbol, Any}(:value => :v_var, :type => Real), :v_for_defaults => Dict{Symbol, Union{Nothing, DataType}}(:value => nothing, :type => Real), :p1 => Dict(:value => nothing)),
-  :structural_parameters => Dict{Symbol, Dict}(:f => Dict(:value => :sin), :N => Dict(:value => 2), :M => Dict(:value => 3))
+  :kwargs                => Dict{Symbol, Dict}(:f=>Dict(:value=>:sin), :N=>Dict(:value=>2), :v=>Dict{Symbol, Any}(:value=>:v_var, :type=>Real), :v_array=>Dict{Symbol, Union{Nothing, UnionAll}}(:value=>nothing, :type=>AbstractArray{Real}), :v_for_defaults=>Dict{Symbol, Union{Nothing, DataType}}(:value=>nothing, :type=>Real), :p1=>Dict(:value=>nothing))
+  :structural_parameters => Dict{Symbol, Dict}(:f=>Dict(:value=>:sin), :N=>Dict(:value=>2))
   :independent_variable  => t
   :constants             => Dict{Symbol, Dict}(:c=>Dict{Symbol, Any}(:value=>1, :type=>Int64, :description=>"Example constant."))
   :extend                => Any[[:p2, :p1], Symbol("#mtkmodel__anonymous__ModelB"), :ModelB]
diff --git a/docs/src/tutorials/programmatically_generating.md b/docs/src/tutorials/programmatically_generating.md
@@ -1,42 +1,44 @@
 # [Programmatically Generating and Scripting ODESystems](@id programmatically)
 
-In the following tutorial we will discuss how to programmatically generate `ODESystem`s.
-This is for cases where one is writing functions that generating `ODESystem`s, for example
-if implementing a reader which parses some file format to generate an `ODESystem` (for example,
-SBML), or for writing functions that transform an `ODESystem` (for example, if you write a
-function that log-transforms a variable in an `ODESystem`).
+In the following tutorial, we will discuss how to programmatically generate `ODESystem`s.
+This is useful for functions that generate `ODESystem`s, for example
+when you implement a reader that parses some file format, such as SBML, to generate an `ODESystem`.
+It is also useful for functions that transform an `ODESystem`, for example
+when you write a function that log-transforms a variable in an `ODESystem`.
 
 ## The Representation of a ModelingToolkit System
 
 ModelingToolkit is built on [Symbolics.jl](https://symbolics.juliasymbolics.org/dev/),
 a symbolic Computer Algebra System (CAS) developed in Julia. As such, all CAS functionality
-is available on ModelingToolkit systems, such as symbolic differentiation, Groebner basis
+is also available to be used on ModelingToolkit systems, such as symbolic differentiation, Groebner basis
 calculations, and whatever else you can think of. Under the hood, all ModelingToolkit
 variables and expressions are Symbolics.jl variables and expressions. Thus when scripting
 a ModelingToolkit system, one simply needs to generate Symbolics.jl variables and equations
 as demonstrated in the Symbolics.jl documentation. This looks like:
 
 ```@example scripting
-using Symbolics
-using ModelingToolkit: t_nounits as t, D_nounits as D
-
-@variables x(t) y(t) # Define variables
+using ModelingToolkit # reexports Symbolics
+@variables t x(t) y(t) # Define variables
+D = Differential(t)
 eqs = [D(x) ~ y
        D(y) ~ x] # Define an array of equations
 ```
 
+However, ModelingToolkit has many higher-level features which will make scripting ModelingToolkit systems more convenient.
+For example, as shown in the next section, defining your own independent variables and differentials is rarely needed.
+
 ## The Non-DSL (non-`@mtkmodel`) Way of Defining an ODESystem
 
-Using `@mtkmodel` is the preferred way of defining ODEs with MTK. However, let us
-look at how we can define the same system without `@mtkmodel`. This is useful for
-defining PDESystem etc.
+Using `@mtkmodel`, like in the [getting started tutorial](@ref getting_started),
+is the preferred way of defining ODEs with MTK.
+However generating the contents of a `@mtkmodel` programmatically can be tedious.
+Let us look at how we can define the same system without `@mtkmodel`.
 
 ```@example scripting
 using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D
-
-@variables x(t)   # independent and dependent variables
-@parameters τ       # parameters
+@variables x(t) = 0.0  # independent and dependent variables
+@parameters τ = 3.0       # parameters
 @constants h = 1    # constants
 eqs = [D(x) ~ (h - x) / τ] # create an array of equations
 
@@ -45,10 +47,16 @@ eqs = [D(x) ~ (h - x) / τ] # create an array of equations
 
 # Perform the standard transformations and mark the model complete
 # Note: Complete models cannot be subsystems of other models!
-fol_model = structural_simplify(model)
+fol = structural_simplify(model)
+prob = ODEProblem(fol, [], (0.0, 10.0), [])
+using DifferentialEquations: solve
+sol = solve(prob)
+
+using Plots
+plot(sol)
 ```
 
-As you can see, generating an ODESystem is as simple as creating the array of equations
+As you can see, generating an ODESystem is as simple as creating an array of equations
 and passing it to the `ODESystem` constructor.
 
 ## Understanding the Difference Between the Julia Variable and the Symbolic Variable
diff --git a/src/systems/model_parsing.jl b/src/systems/model_parsing.jl
@@ -180,16 +180,6 @@ function update_kwargs_and_metadata!(dict, kwargs, a, def, indices, type, var,
     end
 end
 
-function unit_handled_variable_value(mod, y, varname)
-    meta = parse_metadata(mod, y)
-    varval = if meta isa Nothing || get(meta, VariableUnit, nothing) isa Nothing
-        varname
-    else
-        :($convert_units($(meta[VariableUnit]), $varname))
-    end
-    return varval
-end
-
 function parse_variable_def!(dict, mod, arg, varclass, kwargs, where_types;
         def = nothing, indices::Union{Vector{UnitRange{Int}}, Nothing} = nothing,
         type::Type = Real, meta = Dict{DataType, Expr}())
@@ -232,66 +222,6 @@ function parse_variable_def!(dict, mod, arg, varclass, kwargs, where_types;
                 varclass, where_types, meta)
             return var, def, Dict()
         end
-        Expr(:tuple, Expr(:(::), Expr(:ref, a, b...), type), y) || Expr(:tuple, Expr(:ref, a, b...), y) => begin
-            (@isdefined type) || (type = Real)
-            varname = Meta.isexpr(a, :call) ? a.args[1] : a
-            push!(kwargs, Expr(:kw, varname, nothing))
-            varval = unit_handled_variable_value(mod, y, varname)
-            if varclass == :parameters
-                var = :($varname = $first(@parameters $a[$(b...)]::$type = ($varval, $y)))
-            else
-                var = :($varname = $first(@variables $a[$(b...)]::$type = ($varval, $y)))
-            end
-            #TODO: update `dict` aka `Model.structure` with the metadata
-            (:($varname...), var), nothing, Dict()
-        end
-        Expr(:(=), Expr(:(::), Expr(:ref, a, b...), type), y) || Expr(:(=), Expr(:ref, a, b...), y) => begin
-            (@isdefined type) || (type = Real)
-            varname = Meta.isexpr(a, :call) ? a.args[1] : a
-            if Meta.isexpr(y, :tuple)
-                varval = unit_handled_variable_value(mod, y, varname)
-                val, y = (y.args[1], y.args[2:end])
-                push!(kwargs, Expr(:kw, varname, nothing))
-                if varclass == :parameters
-                    var = :($varname = $varname === nothing ? $val : $varname;
-                    $varname = $first(@parameters $a[$(b...)]::$type = (
-                        $varval, $(y...))))
-                else
-                    var = :($varname = $varname === nothing ? $val : $varname;
-                    $varname = $first(@variables $a[$(b...)]::$type = (
-                        $varval, $(y...))))
-                end
-            else
-                push!(kwargs, Expr(:kw, varname, nothing))
-                if varclass == :parameters
-                    var = :($varname = $varname === nothing ? $y : $varname;
-                    $varname = $first(@parameters $a[$(b...)]::$type = $varname))
-                else
-                    var = :($varname = $varname === nothing ? $y : $varname;
-                    $varname = $first(@variables $a[$(b...)]::$type = $varname))
-                end
-            end
-            #TODO: update `dict`` aka `Model.structure` with the metadata
-            (:($varname...), var), nothing, Dict()
-        end
-        Expr(:(::), Expr(:ref, a, b...), type) || Expr(:ref, a, b...) => begin
-            (@isdefined type) || (type = Real)
-            varname = a isa Expr && a.head == :call ? a.args[1] : a
-            push!(kwargs, Expr(:kw, varname, nothing))
-            if varclass == :parameters
-                var = :($varname = $first(@parameters $a[$(b...)]::$type = $varname))
-            elseif varclass == :variables
-                var = :($varname = $first(@variables $a[$(b...)]::$type = $varname))
-            else
-                throw("Symbolic array with arbitrary length is not handled for $varclass.
-                    Please open an issue with an example.")
-            end
-            dict[varclass] = get!(dict, varclass) do
-                Dict{Symbol, Dict{Symbol, Any}}()
-            end
-            # dict[:kwargs][varname] = dict[varclass][varname] = Dict(:size => b)
-            (:($varname...), var), nothing, Dict()
-        end
         Expr(:(=), a, b) => begin
             Base.remove_linenums!(b)
             def, meta = parse_default(mod, b)
@@ -338,6 +268,11 @@ function parse_variable_def!(dict, mod, arg, varclass, kwargs, where_types;
             end
             return var, def, Dict()
         end
+        Expr(:ref, a, b...) => begin
+            indices = map(i -> UnitRange(i.args[2], i.args[end]), b)
+            parse_variable_def!(dict, mod, a, varclass, kwargs, where_types;
+                def, indices, type, meta)
+        end
         _ => error("$arg cannot be parsed")
     end
 end
@@ -445,23 +380,14 @@ function parse_default(mod, a)
     end
 end
 
-function parse_metadata(mod, a::Expr)
+function parse_metadata(mod, a)
     MLStyle.@match a begin
-        Expr(:vect, b...) => Dict(parse_metadata(mod, m) for m in b)
-        Expr(:tuple, a, b...) => parse_metadata(mod, b)
+        Expr(:vect, eles...) => Dict(parse_metadata(mod, e) for e in eles)
         Expr(:(=), a, b) => Symbolics.option_to_metadata_type(Val(a)) => get_var(mod, b)
         _ => error("Cannot parse metadata $a")
     end
 end
 
-function parse_metadata(mod, metadata::AbstractArray)
-    ret = Dict()
-    for m in metadata
-        merge!(ret, parse_metadata(mod, m))
-    end
-    ret
-end
-
 function _set_var_metadata!(metadata_with_exprs, a, m, v::Expr)
     push!(metadata_with_exprs, m => v)
     a
@@ -719,7 +645,6 @@ function parse_variable_arg!(exprs, vs, dict, mod, arg, varclass, kwargs, where_
 end
 
 function convert_units(varunits::DynamicQuantities.Quantity, value)
-    value isa Nothing && return nothing
     DynamicQuantities.ustrip(DynamicQuantities.uconvert(
         DynamicQuantities.SymbolicUnits.as_quantity(varunits), value))
 end
@@ -731,7 +656,6 @@ function convert_units(
 end
 
 function convert_units(varunits::Unitful.FreeUnits, value)
-    value isa Nothing && return nothing
     Unitful.ustrip(varunits, value)
 end
 
@@ -750,50 +674,47 @@ end
 function parse_variable_arg(dict, mod, arg, varclass, kwargs, where_types)
     vv, def, metadata_with_exprs = parse_variable_def!(
         dict, mod, arg, varclass, kwargs, where_types)
-    if !(vv isa Tuple)
-        name = getname(vv)
-        varexpr = if haskey(metadata_with_exprs, VariableUnit)
-            unit = metadata_with_exprs[VariableUnit]
-            quote
-                $name = if $name === nothing
-                    $setdefault($vv, $def)
-                else
-                    try
-                        $setdefault($vv, $convert_units($unit, $name))
-                    catch e
-                        if isa(e, $(DynamicQuantities.DimensionError)) ||
-                           isa(e, $(Unitful.DimensionError))
-                            error("Unable to convert units for \'" * string(:($$vv)) * "\'")
-                        elseif isa(e, MethodError)
-                            error("No or invalid units provided for \'" * string(:($$vv)) *
-                                  "\'")
-                        else
-                            rethrow(e)
-                        end
+    name = getname(vv)
+
+    varexpr = if haskey(metadata_with_exprs, VariableUnit)
+        unit = metadata_with_exprs[VariableUnit]
+        quote
+            $name = if $name === nothing
+                $setdefault($vv, $def)
+            else
+                try
+                    $setdefault($vv, $convert_units($unit, $name))
+                catch e
+                    if isa(e, $(DynamicQuantities.DimensionError)) ||
+                       isa(e, $(Unitful.DimensionError))
+                        error("Unable to convert units for \'" * string(:($$vv)) * "\'")
+                    elseif isa(e, MethodError)
+                        error("No or invalid units provided for \'" * string(:($$vv)) *
+                              "\'")
+                    else
+                        rethrow(e)
                     end
                 end
             end
-        else
-            quote
-                $name = if $name === nothing
-                    $setdefault($vv, $def)
-                else
-                    $setdefault($vv, $name)
-                end
-            end
         end
-
-        metadata_expr = Expr(:block)
-        for (k, v) in metadata_with_exprs
-            push!(metadata_expr.args,
-                :($name = $wrap($set_scalar_metadata($unwrap($name), $k, $v))))
+    else
+        quote
+            $name = if $name === nothing
+                $setdefault($vv, $def)
+            else
+                $setdefault($vv, $name)
+            end
         end
+    end
 
-        push!(varexpr.args, metadata_expr)
-        return vv isa Num ? name : :($name...), varexpr
-    else
-        return vv
+    metadata_expr = Expr(:block)
+    for (k, v) in metadata_with_exprs
+        push!(metadata_expr.args,
+            :($name = $wrap($set_scalar_metadata($unwrap($name), $k, $v))))
     end
+
+    push!(varexpr.args, metadata_expr)
+    return vv isa Num ? name : :($name...), varexpr
 end
 
 function handle_conditional_vars!(
diff --git a/test/model_parsing.jl b/test/model_parsing.jl
