# Help in making multi-dimensional equations a bit sexier

I need to write some equations that loop over several dimensions, and I created my own function `getData(parName,dim1,dim2)` and `setData(value,dim1,dim2)` that, using anonymous functions, I dynamically wrapped to `parName(dim1,dim2)` and `parName!(value,dim1,dim2)`.

So I can now write my model as:

``````[par1!( par2(d1,d2)+par3(d1,d2)   ,d1,d2,dfix3) for d1 in DIM1, d2 in DIM2]
``````

where the first argument of `par1!()` function is the value that I want to store.

Now, it’s still not very inspiring, so I am trying to make this syntax a bit more easy to write (and to read/maintain).

Specifically, I am trying to write a macro that let’s write the equation as:

``````@meq par1!(d1 in DIM1, d2 in DIM2, dfix3) =  par2(d1,d2)+par3(d1,d2)
``````

That’s sexy!

While I have used macros in julia, it’s the first time I am trying to write one.

So far, I managed to write:

``````macro meq(eq)
dump(eq)
lhs_par               = eq.args.args
rhs                   = eq.args
lhs_dims              = eq.args.args[2:end]
loop_counters         = [d.args for d in lhs_dims if typeof(d) == Expr]
loop_sets             = [d.args for d in lhs_dims if typeof(d) == Expr]
loop_wholeElements    = [d for d in lhs_dims if typeof(d) == Expr]
lhs_dims_placeholders = []
for d in lhs_dims
if typeof(d) == Expr
push!(lhs_dims_placeholders,d.args)
else
push!(lhs_dims_placeholders,d)
end
end
outExp =  quote
[\$(lhs_par)(\$(rhs),\$(lhs_dims_placeholders ...)) for \$(loop_wholeElements ...)  ]
end
return outExp
end
``````

While I think “I am almost there”, unfortunately the above code return a syntax error (“invalid iteration specification”) due to the `for \$(loop_wholeElements)` part… indeed I don’t know how to expand the expressions in lhs_dims_placeholders and loop_wholeElements in order to “assemble” the expanded expression…

I can’t help with writing macros, but you may want to check out cartesian indexing as a way to save you some time.

My instinct is that there may be a better way to do this. Are your your parameter structures actually multi-dimensional arrays? If so, have you considered making them subtypes of `AbstractArray`? You would just need to implement `size()`, `getindex()`, `setindex()` and `IndexStyle()` https://docs.julialang.org/en/stable/manual/interfaces/#man-interface-array-1 and you’d get multi-dimensional indexing, iteration, and lots of other features for free. In particular, you’d be able to use Base.Cartesian to generate multi-dimensional loops over your data.

Not sure about your particular application here, but have you considered that using a macro might not be appropriate? Have you tried using a regular function that takes an `Expr` as input and then returns an `Expr` as output, and then using `eval` on it? Perhaps that’s worth a try first, before going the macro route.

I think your problem might be the `...` in

``````\$(loop_wholeElements)
``````

Thank you. My actual data is in the form of a DataFrame parName|region|d1|d2|y|value.

The reason I am using my own get/setData() functions is that I implement my own application-specific query logic on the data.
For example some variables are known at regional level, but other only an national one, so if data is not available at regional level the query pick it up it at national one, or if a variable is not available for that year I pick it up for the last available year…

Still it doesn’t compile it without it

`show(loop_wholeElements)` returns `Expr[:(d1 in DIM1), :(d2 in DIM2)]` and, when I leave the outExp as

``````outExp =  quote
[\$(lhs_par)(\$(rhs),\$(lhs_dims_placeholders ...)) for d1 in DIM1, d2 in DIM2  ]
end
``````

(to let the macro compile) `show(outExp)` returns:

``````quote  # /home/lobianco/git/ffsm/test.jl, line 47:
[par1!(begin  # /home/lobianco/git/ffsm/test.jl, line 54:
par2(d1, d2) + par3(d1, d2)
end, d1, d2, dfix3) for d1 = DIM1, d2 = DIM2]
end
``````

I don’t understand why the output expression has (a) the begin/end parts wrapping the replaced lhs expression; (b) the code `for d1 in DIM1, d2 in DIM2` replaced with `for d1 = DIM1, d2 = DIM2`

Here, try this, I think this should get you closer to what you want, it doesn’t give you the error about iterators and makes the expression you are looking for, i believe

``````return :([\$(Expr(:generator,:(\$(Expr(:call,lhs_par,rhs,lhs_dims_placeholders...))),loop_wholeElements...))])
``````

Although, you might want to replace `rhs` with something like `rhs.args`, if you want to avoid having the extra quote block in your final expression.

1 Like

I’m wondering why the following two expressions behave differently on the scope of `x`:

``````julia> macro bar(ex)
Expr(:comprehension, 1, ex)
end
@bar (macro with 1 method)

julia> @macroexpand @bar(x in X)
:(\$(Expr(:comprehension, 1, :(Main.x in Main.X))))

julia> @macroexpand @bar(x = X)
:(\$(Expr(:comprehension, 1, :(#50#x = Main.X))))
``````

The result shows that the `x` of `x in X` is in the global scope, but the `x` of `x = X` is in the local scope, and both of the `X` are in the global scope.

Thank you… my knowledge about metaprogrammming is too limited to fully understand the issue… just looks odd that, for example, if I use this expression with your macro:

``````@meq  pl!(tp in secProducts, r in fr2,"",y2) = sum(pl_(r,pp,"",y2)*a_(r,pp,tp,y2) for pp in priProducts) + m_(r,tp,"",y2)
``````

I obtain:

`````` :(\$(Expr(:comprehension, :(pl!(begin  # /home/lobianco/git/ffsm/ffsm.jl, line 127:
sum((pl_(r, pp, "", y2) * a_(r, pp, tp, y2) for pp = priProducts)) + m_(r, tp, "", y2)
end, tp, r, "", y2)), :(tp = secProducts), :(r = fr2))))
``````

Or, if I use `rhs.args` instead of `args`, I obtain:

``````:(\$(Expr(:comprehension, :(pl!(sum((pl_(r, pp, "", y2) * a_(r, pp, tp, y2) for pp = priProducts)) + m_(r, tp, "", y2), tp, r, "", y2)), :(tp = secProducts), :(r = fr2))))
``````

(notice in both cases the transformation of `for pp in priProducts` to `for pp = priProducts`)

As an OT question, as SO and discourse.julialang.org are two separate but not disjoint communities, is it considered legit to post on SO and here ?

Many thanks to everyone that helped me with this question.

Double posting is discouraged, but if you choose to do it you should post the cross-reference both places.

2 Likes

Just wrap your data in a new type (possibly a subtype of `AbstractArray`), and define `getindex` and `setindex!` etcetera functions on that type. Probably don’t use macros for this. Then you can just do:

``````for d1 in DIM1, d2 in DIM2
par1[d1,d2,dfix3] = par2[d1,d2] + par3[d1,d2]
end
``````
1 Like

Yes, I see now my original solution has an extra parenthesis in it.

Here is the final corrected version, this should give you exactly what you need:

``````return Expr(:comprehension,Expr(:generator,Expr(:call,lhs_par,rhs.args,lhs_dims_placeholders...),loop_wholeElements...))
``````

Although, you need to be careful about using `rhs.args` if you are not entering your command from the REPL, in that case you need to replace `` with ``.

It returns

``````:([par1!(par2(d1, d2) + par3(d1, d2), d1, d2, dfix3) for d1 in DIM1, d2 in DIM2])
``````

Alright, so I went ahead and tested this. I was not able to make the `macro` work yet, but as a `function` it works (also, in order for it to work you need `=` not `in` and change the array element for `d.args`:

``````par1!(a,b,c,d) = 0; par2(a,b) = 0; par3(a,b) = 0; dfix3 = 0;
function meq(eq)
lhs_par               = eq.args.args
rhs                   = eq.args
lhs_dims              = eq.args.args[2:end]
loop_counters         = [d.args for d in lhs_dims if typeof(d) == Expr]
loop_sets             = [d.args for d in lhs_dims if typeof(d) == Expr]
loop_wholeElements    = [d for d in lhs_dims if typeof(d) == Expr]
lhs_dims_placeholders = []
for d in lhs_dims
if typeof(d) == Expr
push!(lhs_dims_placeholders,d.args)
else
push!(lhs_dims_placeholders,d)
end
end
return Expr(:comprehension,Expr(:generator,Expr(:call,lhs_par,rhs.args,lhs_dims_placeholders...),loop_wholeElements...))
end
``````

Then computing the result like this

``````julia> meq(:(par1!(d1 = 1:2, d2 = 1:2, 3) =  par2(d1,d2)+par3(d1,d2))) |> eval

2×2 Array{Int64,2}:
0  0
0  0
``````

Thank you… both versions now work I am using them to write a bio-economic model with a syntax similar to that of JuMP, so I can use JuMP for the optimisation part, and this syntax for everything else (but I am also testing the suggestion of stevengj):

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

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

# set total supply (local + abroad[exports])
@meq st!(r in fr2, p in products, "", y2) = sl_(r, p, "", y2) + sa_(r, p, "", y2)

# set total demand (local + abroad[imports])
@meq dt!(r in fr2, p in products, "", y2) = da_(r, p, "", y2) + dl_(r, p, "", y2)

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

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

# weighted price
@meq pw!(r in fr2, tp in secProducts) = (dl_(r,tp)*pl_(r,tp) + da_(r,tp)*pl_(wld,tp))/dl_(r,tp)+da_(r,tp)
``````

Just another general remark, for metaprogramming I have found that using the `Expr(::Symbol,args...)` constructor is overall much more robust and reliable than using the `\$` interpolation. The best way to get a first look at a how to construct an unusual expression correctly is to `dump` it and then inspect the `head` and `args` to see how the `Expr` constructor can be applied. With interpolation, sometimes weird things happen.

1 Like

for j in 1:nhubs, for k in 1:nplant
if supply_plant[j,k]==1
println(“supply_plant[\$j,\$k]=”,supply_plant[j,k])
else
end
end

i am getting this error…ERROR: syntax: invalid iteration specification

Please tell me how to resolve it

If you have a new question, please create a new thread rather than posting on an old thread.

1 Like

But since it’s here

`for j in 1:nhubs, k in 1:nplant`
should be fine