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 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:

```
# *** 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))
```