2016-08-11 4 views
2

NFFTパッケージには、ベクトルと行列のための高速で特殊な関数と一般的な配列の遅い関数があります。私は一般的な機能を速くし、特殊な機能は時代遅れにしたいと思っていますが、私は困っています。Julia NFFTでは、ベクトルと行列のメソッドを単一の配列メソッドで置き換えます。

myfunc!(f::Vector, g::Vector) 
    offset = 5 
    n = length(g) 
    N = length(f) 

    for l = 1:N 
     g[ ((l+offset)%n)+1 ] = f[l] 
    end 
end 

myfunc!(f::Matrix, g::Matrix) 
    offsetx = 5 
    offsety = 5 
    n1, n2 = size(g) 
    N1, N2 = size(f) 

    for ly = 1:N2 
     for lx = 1:N1 
      g[ ((lx+offsetx)%n1)+1, ((ly+offsety)%n2)+1 ] = f[lx, ly] 
     end 
    end 
end 

を:特化した機能の

一つの家族は(offsetは、実際のコードで計算されますが、そのは、この質問のために重要ではありません)このように、本質的に見えます一般的な関数は次のように書くことができます。

myfunc!{D}(f::Array{D}, g::Array{D}) 
    offset = ntuple(d -> 5, D) 
    n = size(g) 

    for l in CartesianRange(size(f)) 
     idx = CartesianIndex{D}(ntuple(d -> ((l[d]+offset[d])%n[d])+1, D)) 
     g[ idx ] = f[l] 
    end 
end 

残念ながら、これは非常に遅いです。大部分の時間はループ内のntupleで費やされます。 idxため

他の可能性はidx = Array{Int}(D)をさせると、内部ループが

for d = 1:D 
    idx[d] = ((l[d]+offset[d])%n[d])+1 
end 
g[idx...] = f[l] 

これも低速であるように見えるように含まれます。私が考えている

その次元D@generated機能がidxを計算するために行うことができるタイプの引数、ですが、私はそれを行う方法を見つけ出すことはできませんので、(またはより良い方法がある場合)。

私はJulia v0.4.5を使用しています。

答えて

1

Base.Cartesianヘルパーマクロを使用するには、生成された関数でこれを行う方法の答えは です。

using Base.Cartesian 

@generated function myfunc!{T,D}(f::Array{T,D}, g::Array{T,D}) 
    quote 
     @nexprs $D d->offset_d=5 #Declase offset_1=5, offset_2=5 etc 

     @nloops $D l f begin 
      (@nref $D g d->(l_d+offset_d)%size(g,d)+1) = @nref $D f l 
     end 
    end 
end 

私はこれが正しいことを2Dで確認しました。 私はそれをプロファイルするために読者に消費税として残します。 多かれ少なかれ割り当てられません。

+0

非常に良い解決策!このアプローチで私はOPを単純化しました。私はまた3番目の入力 'h :: Vector {Vector {Float64}}'を持っていて、内部ループは 'g [((lx + offsetx)%n1)+1、((ly + (lx + offsetx)%n1)+1、((ly + offsety)%n2)+1] = f [lx(0 + 2)+1)= f [1x、ly] 、ly] * h [1] [lx] * h [2] [ly] 'となります。 これらの用語は 'Base.Cartesian'を使って含めることもできます。それはネストされた '@ nref'のようなものはないと思われます。 – Robert

+1

確かに:最も内側のセクション( '@ nref')を: ' v = @nref $ D f l'に変更してください。 '@nexprs $ D d-> v * = h [d] [l_d]'; (g_d + offset_d)%size(g、d)+1)= v'; –

関連する問題