Help with making multidimensional algorithm type generic

New to Julia
1

Hi Everyone!
first i’d like to say that i’m really enjoying using Julia, i’m using it in my PhD with increasing success!
lately i’ve been developing a MeanShift algorithm and i managed to make the code dimension generic but still not fully type generic

using StaticArrays,LinearAlgebra

@inline function Kern(x,h)
    max(zero(h),1-norm(x./h)^2)
end
function circledim(M,c,i)
    mod(c-1,size(M,i))+1
end
function makerange(p,h,i)
    round(Int,p[i]-h-1):round(Int,p[i]+h+1)
end
@generated function genms(M::Array{T,N},v,h::T) where {T,N}
    quote
        @inbounds begin
            p=SVector{$N}(v)
            c=SVector{$N}(v.*0)

            W=zero(T)
            xsq=zero(T)
            ls=1/(π*h^2)

            Base.Cartesian.@nexprs $N i -> R_i = makerange(p,h,i)
            Base.Cartesian.@nloops $N i i->R_i begin
                cp=SVector(Base.Cartesian.@ntuple $N k->i_k)
                sW=Kern(p.-cp,h)*((Base.Cartesian.@nref $N M k->circledim(M,i_k,k)) + ls)
                xsq+=norm(p.-cp)^2*sW
                c+=cp*sW
                W+=sW
            end
            c/=W
            xsq/=W
            c=2*c-p
            return SVector(Base.Cartesian.@ntuple $N k->circledim(M,c[k],k)),sqrt(3*xsq),W/xsq,norm(c.-p)
        end
    end
end
function MWE()
    T=[(x,y) for x in -5:0.1:5, y in -5:0.1:5]
    M=[exp(-(x-1)^2/2)* exp(-(y-2)^2/2) for x in -5:0.1:5, y in -5:0.1:5]
    v=@SVector [50.0,50.0]
    h=5.0
    del=1.0
    while del>1e-5
        v,h,W,del=genms(M,v,h)
    end
    M[round.(Int,v)...],T[round.(Int,v)...],h
end
MWE()

I struggled a lot to use the macro magic contained in Base.Cartesian, could someone help me with making the code fully generic?

the result of MWE() is (1.0, (1.0, 2.0), 3.3222927225309373)

2

As a start: given that you’re using an SVector for the base point, why not make M a vector of SVector points too? Then I think you can write this fairly easily in generic form without using a @generated function or Base.Cartesian.

2 Likes
3

Thank you Tim,
:: why not make M a vector of SVector points too?

M in my use case is a big multidimensional matrix obtained by loading images, it requires a lot of time to transform it into a SArray, so i’d prefer a way to avoid specifying the nature of the arrays, since in the future we plan to use GPU’s also (so i’m not sure what the final array type is going to be)

:: Then I think you can write this fairly easily in generic form without using a @generated function or Base.Cartesian.

what language constructs allow me to generate N nested loops at compile time other than Base.Cartesian?

the docs say:

“Starting in Julia 0.4-pre, the recommended approach is to use a @generated function .”

but also this

“The (non-exported) Cartesian module provides macros that facilitate writing multidimensional algorithms. It is hoped that Cartesian will not, in the long term, be necessary; however, at present it is one of the few ways to write compact and performant multidimensional code.”

So perhaps this module is no longer necessary, and there are more idiomatic/transparent ways to do the same thing, but i didn’t manage to find any

4

In your case it may be necessary, but much easier approaches support many constructs: Multidimensional algorithms and iteration

Which docs were you following? Sounds like we may need to update something, even if though there are some cases where it’s still necessary to use Base.Cartesian.

1 Like
5

i was looking here Base.Cartesian doc page
i attempted

using StaticArrays,LinearAlgebra

@inline function Kern(x,h)
    max(zero(h),1-norm(x./h)^2)
end
@inline function circledim(M,c,i)
    mod(c-1,size(M,i))+1
end
@inline function makerange(p,h,i)
    round(Int,p[i]-h-1):round(Int,p[i]+h+1)
end
function genms(M::Array{T,N},v::Vector{T},h::T) where {T,N}
        @inbounds begin
            c=zeros(T,N)
            W=xsq=zero(T)
            ls=1/(π*h^2)
            CI=CartesianIndices(Tuple(makerange(v,h,i) for i in 1:N))
            for cp in CI
                tcp=Tuple(cp)
                sW=Kern(v.-tcp,h)*(M[Tuple(circledim(M,tcp[i],i) for i in 1:N)...] + ls)
                xsq+=norm(v.-tcp)^2*sW
                c.+=tcp.*sW
                W+=sW
            end
            xsq/=W
            c.=2 .* c ./ W .-v
            return [circledim(M,c[i],i) for i in 1:N],sqrt(3*xsq),W/xsq,norm(c.-v)
        end
end
function MWE()
    T=[(x,y) for x in -5:0.1:5, y in -5:0.1:5]
    M=[exp(-(x-1)^2/2)* exp(-(y-2)^2/2) for x in -5:0.1:5, y in -5:0.1:5]
    v=[50.0,50.0]
    h=5.0
    del=1.0
    @time while del>1e-5
        v,h,W,del=genms(M,v,h)
    end
    M[round.(Int,v)...],T[round.(Int,v)...],h
end
MWE()

The only problem is that this new version runs ~500 times slower and allocates a lot…
while i had no allocation with the older approach

6

ProfileView.jl is your friend. Check out the original profile:

profileview

All that red on top is a sign of runtime dispatch due to non-inferrability (“type instability”). code_warntype/@code_warntype are then useful to diagnose the problem in more detail.

Try changing your two Tuple(f(i, x, ...) for i = 1:N)s into ntuple(i->f(i, x, ...), Val(N)). That will get rid of your type instability. Then I think you’ll find opportunities for further improvement (e.g., using StaticArrays, possibly more efficient alternatives to mod, etc). ProfileView can be your guide here too. Just make sure you run it long enough that you collect adequate numbers of samples.

With regards to making it generic, you should be able to replace the type declarations ::Array{T,N} and similar with ::AbstractArray{T,N}.

3 Likes