Help in making multi-dimensional equations a bit sexier

General Usage
,
1

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[1].args[1]
    rhs                   = eq.args[2]
    lhs_dims              = eq.args[1].args[2:end]
    loop_counters         = [d.args[2] for d in lhs_dims if typeof(d) == Expr]
    loop_sets             = [d.args[3] 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[2])
        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…

2

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

3

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.

4

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.

5

I think your problem might be the ... in

$(loop_wholeElements)
6

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…

7

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

8

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[2], if you want to avoid having the extra quote block in your final expression.

9

x-ref: Julia macro for transforming `f(dim1,dim2,..) = value` into `f(value,dim1,dim2,..)` - Stack Overflow

1 Like
10

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.

11

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[2] 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)

In any case your macro seems to still working, I accepted your answer on S.O.

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.

12

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

2 Likes
13

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
14

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[2],lhs_dims_placeholders...),loop_wholeElements...))

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

It returns

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

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[2]:

par1!(a,b,c,d) = 0; par2(a,b) = 0; par3(a,b) = 0; dfix3 = 0;
function meq(eq)
    lhs_par               = eq.args[1].args[1]
    rhs                   = eq.args[2]
    lhs_dims              = eq.args[1].args[2:end]
    loop_counters         = [d.args[2] for d in lhs_dims if typeof(d) == Expr]
    loop_sets             = [d.args[2] 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[2])
        else
            push!(lhs_dims_placeholders,d)
        end
    end
    return Expr(:comprehension,Expr(:generator,Expr(:call,lhs_par,rhs.args[2],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
16

Thank you… both versions now work :slight_smile:

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)
17

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
18

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

19

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

1 Like
20

But since it’s here

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