# Any concise way to express simulation equations in Julia (like algebraic modelling language but not necessarily linked to optimisation problems)?

**URL:** <https://discourse.julialang.org/t/any-concise-way-to-express-simulation-equations-in-julia-like-algebraic-modelling-language-but-not-necessarily-linked-to-optimisation-problems/52449>\
**Category:** Data\
**Created:** [December 27, 2020, 11:02am UTC](https://discourse.julialang.org/t/any-concise-way-to-express-simulation-equations-in-julia-like-algebraic-modelling-language-but-not-necessarily-linked-to-optimisation-problems/52449 "2020-12-27T11:02:20Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [December 27, 2020, 11:02am UTC](https://discourse.julialang.org/t/any-concise-way-to-express-simulation-equations-in-julia-like-algebraic-modelling-language-but-not-necessarily-linked-to-optimisation-problems/52449/1 "2020-12-27T11:02:20Z")

</div>

I know several [Algebraic Modelling Languages](https://en.wikipedia.org/wiki/Algebraic_modeling_language) in the context of optimisation, e.g. [GAMS](https://www.gams.com/), [Pyomo](https://www.pyomo.org/) (Python) or [JuMP](https://jump.dev/) (Julia). These allow you to express your objective function and your constraints with a syntax very concise and similar to your pen and paper equations.  
But… outside the domain of optimisation, when you have a model with hundred of equations like the following pseudocode, which is the simplest way to express and compute these ?

```julia
@meq nTrees!(r in reg, sp in species, dc in diameterClass[2-end], y in years) = nTrees_(r, sp, dc, y-1)*(1-mortRate_(r, sp, dc, y-1) - promotionRate_(r, sp, dc, y-1))) + nTrees_(r, sp, dc-1, y-1) * promotionRate_(r, sp, dc-1, y-1))
....
production!(r in reg, sp in secPr, y in years) = sum(trValues_(r,pp) * transfCoef_(r,pp) for pp in primPr)

```

I assume a matrix formulation of the equations or maybe using for loops/list comprehension/maps are the most efficient solutions, but they become cumbersome when you have many dimensions and you need to operate on a subset of them, or when you have to operates transforms like the summation above…

What is your approach? Is there a library that allows to write the kind of equations above in a concise way ? Is such library needed ? (I did wrote one some time ago based on IndexedTable.jl and have to decide if it is worth to update it to Julia 1.x or if there are better approaches)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [December 27, 2020, 1:19pm UTC](https://discourse.julialang.org/t/any-concise-way-to-express-simulation-equations-in-julia-like-algebraic-modelling-language-but-not-necessarily-linked-to-optimisation-problems/52449/2 "2020-12-27T13:19:37Z")

</div>

> **[GitHub - SciML/ModelingToolkit.jl: An acausal modeling framework for...](https://github.com/SciML/ModelingToolkit.jl)**
>
> An acausal modeling framework for automatically parallelized scientific machine learning (SciML) in Julia. A computer algebra system for integrated symbolics for physics-informed machine learning a...

So far it covers:

- ODEs
- SDEs
- Jump dynamics
- Nonlinear solving
- Nonlinear optimization
- Control
- Reaction systems
- PDEs

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [December 27, 2020, 3:04pm UTC](https://discourse.julialang.org/t/any-concise-way-to-express-simulation-equations-in-julia-like-algebraic-modelling-language-but-not-necessarily-linked-to-optimisation-problems/52449/3 "2020-12-27T15:04:16Z")

</div>

Thank you, but my needs are more on how to manipulate data (likely one of the two commonly used data structures in Julia, IndexedTables or DataFrames) in a concise way…

---

<div class="post-metadata">

**Author:** ![mikkoku](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikkoku/32/16274_2.png) [@mikkoku](https://discourse.julialang.org/u/mikkoku)\
**Post date:** [December 28, 2020, 10:26am UTC](https://discourse.julialang.org/t/any-concise-way-to-express-simulation-equations-in-julia-like-algebraic-modelling-language-but-not-necessarily-linked-to-optimisation-problems/52449/4 "2020-12-28T10:26:07Z")

</div>

I would use an [inner join](https://dataframes.juliadata.org/stable/man/joins/) for the first example and [split-apply-combine](https://dataframes.juliadata.org/stable/man/split_apply_combine/) for the second.  
There are also several data manipulation packages like [Query.jl](https://www.queryverse.org/Query.jl/stable/) that might allow a more concise syntax for the same operations.

---

<div class="post-metadata">

**Author:** ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)\
**Post date:** [December 28, 2020, 10:41am UTC](https://discourse.julialang.org/t/any-concise-way-to-express-simulation-equations-in-julia-like-algebraic-modelling-language-but-not-necessarily-linked-to-optimisation-problems/52449/5 "2020-12-28T10:41:10Z")

</div>

Thank you, but even using Query.jl, it would take lot of efforts to write every single equation/transformation.

Here is an extract of the model I am playing with, using the (old, not working) [MultiDimEquation.jl](https://github.com/sylvaticus/MultiDimEquations.jl) package I am considering to revive or looking if something similar (and better) exists already. Note that the functions are auto-generated with `Foo_(dim1,dim2)` a get\_like function and `Foo!(value,dim1,dim2)` a set-like function:

```julia
# ***Init market module***

# Setting da for prim products and sa for sec products for each region from the data at national level
@meq sa!(r in fr2, tp in secProducts,y2) = sa_("FRA",tp,y2)
@meq da!(r in fr2, pp in priProducts,y2) = da_("FRA",pp,y2)

# set local prices for transformed products
@meq pl!(r in fr2, tp in secProducts,y2) = sum(pl_(r,pp,y2)*a_("FRA",pp,y1,tp) for pp in priProducts) + m_(r,tp,y1)

# set local demand for primary products
# [dl!( sum(sl_(r,tp) * a_(r,pp,tp) for tp in secProducts) ,r, pp) for r in fr2, pp in priProducts]
@meq dl!(r in fr2, pp in priProducts, y2) = sum(sl_(r,tp,y2) * a_("FRA",pp,y1,tp) for tp in secProducts)

# set total supply (local + abroad[exports]) and total demand (local + abroad[imports])
@meq st!(r in fr2, p in products, y2) = sl_(r, p, y2) + sa_(r, p, y2)
@meq dt!(r in fr2, p in products, y2) = da_(r, p, y2) + dl_(r, p, y2)

# set initial values for the ces function (transf products)
@meq q1!(r in fr2, tp in secProducts,y1) = 1.0 / (1.0+((dl_(r,tp,y2)/da_(r,tp,y2))^(1/psi_(r,tp,y1)) ) * (pl_(r,tp,y2)/pl_(wld,tp,y2)) )
@meq p1!(r in fr2, tp in secProducts,y1) = 1 - q1_(r, tp,y1)

# composite demand/price (sec products)
@meq dc!(r in fr2, tp in secProducts,y2) = (q1_(r,tp,y1)*da_(r,tp,y2)^( (psi_(r,tp,y1)-1)/psi_(r,tp,y1) )+p1_(r,tp,y1)*dl_(r,tp,y1)^( (psi_(r,tp,y1)-1)/psi_(r,tp,y1) ) )^(psi_(r,tp,y1)/(psi_(r,tp,y1)-1))
@meq pc!(r in fr2, tp in secProducts,y2) = (da_(r,tp,y2)/dc_(r,tp,y2)) * pl_(wld,tp,y2) + (dl_(r,tp,y2)/dc_(r,tp,y2)) * pl_(r,tp,y2)

# Capacity parameter
@meq k!(r in fr2, tp in secProducts,y2) = k1_("FRA",tp,y1) * sl_(r,tp,y2)

# set initial values for the ces function (primary products)
@meq t1!(r in fr2, pp in priProducts,y1) = 1.0 / (1.0+((sl_(r,pp,y2)/sa_(r,pp,y2))^(1/eta_(r,pp,y1)) ) * (pl_(r,pp,y2)/pl_(wld,pp,y2)) )
@meq r1!(r in fr2, pp in priProducts,y1) = 1 - t1_(r, pp,y1)

# composite demand/price (prim products)
@meq sc!(r in fr2, pp in priProducts,y2) = (t1_(r,pp,y1)*sa_(r,pp,y2)^( (eta_(r,pp,y1)-1)/eta_(r,pp,y1) )+r1_(r,pp,y1)*sl_(r,pp,y2)^( (eta_(r,pp,y1)-1)/eta_(r,pp,y1) ) )^(eta_(r,pp,y1)/(eta_(r,pp,y1)-1))
@meq pc!(r in fr2, pp in priProducts, y2) = (sa_(r,pp,y2)/sc_(r,pp,y2)) * pl_(wld,pp,y2) + (sl_(r,pp,y2)/sc_(r,pp,y2)) * pl_(r,pp,y2)
# psi_("LO","hardWRoundW",y1)

# weighted prices
@meq pw!(r in fr2, tp in secProducts,y2) = (dl_(r,tp,y2)*pl_(r,tp,y2) + da_(r,tp,y2)*pl_(wld,tp,y2))/ (dl_(r,tp,y2)+da_(r,tp,y2))
# at this point there is an error, as some outputs are NaN!!
@meq pw!(r in fr2, pp in priProducts,y2) = (sl_(r,pp,y2)*pl_(r,pp,y2) + sa_(r,pp,y2)*pl_(wld,pp,y2))/(sl_(r,pp,y2)+sa_(r,pp,y2))

# setting regional trade to zero and positive transport costs between same region
@meq rt!(r in fr2, p in products, y2, r2 in fr2) = 0.0
[ct!(0.1, rFrom, p, y1, rTo ) for rFrom in fr2, p in products, rTo in fr2 if rFrom == rTo]

# available inventories
@meq in!(r in fr2, pp in priProducts, y2) = sum(app_(pp,ft,dc)* vol_(r,ft,dc,y1) for ft in fTypes, dc in invDClasses)
# [println("$r \t $p \t $(in_(r,p,y2))") for r in fr2, p in priProducts if r == "AQ"]

# ***Forest dynamic module***

# Harvesting rate
@meq hr!(r in fr2, ft in fTypes, dc in invDClasses, y2) = sum(app_(pp,ft,dc) * st_(r,pp,y2) / in_(r, pp, y2) for pp in priProducts)
@meq hV!(r in fr2, ft in fTypes, dc in invDClasses, y2) = hr_(r, ft, dc, y2) * vol_(r, ft, dc, y1)

# volumes.. first diameter class and then subsequent diameter classes..
@meq (vol!(r in fr2, ft in fTypes, invDClasses[1], y2) =
           vol_(r, ft, invDClasses[1], y1) * (1 - mortCoef_(r, ft, invDClasses[1], y1) - (1/tp_(r, ft, invDClasses[1], y1)) -hr_(r, ft, invDClasses[1], y2))
        + vReg_(r, ft, invDClasses[1], y1))

@meq (vol!(r in fr2, ft in fTypes, dc in invDClasses[2:end], y2) =
           vol_(r, ft, dc, y1) * (1-mortCoef_(r, ft, dc, y1) - 1/tp_(r, ft, dc, y1)-hr_(r, ft, dc, y2))
        + vol_(r, ft, dc-1, y1) * (1/tp_(r, ft, dc-1, y1)) * betaCoef_(r, ft, dc-1, y1))

```