[dennou-ruby:002287] Re: Ruby-ISPACK 作成中

[Takeshi Horinouchi <horinout@xxxxxxxxxxxxxxxxxx>]

堀之内です．

＃ 亀レスご容赦を

> 北村です．
> 
> > 堀之内さんに言ったら作業変数を隠蔽するような上位のラッパーもほ 
> > しい
> > っていうことですのでそれもおいおいできればいいかなと思います。
> 
> 自分が利用しているN2PACKの一部だけですが，ラップしてみました．
> 作業変数の隠蔽と，グリッドデータ，スペクトルデータが1次元 
> 配列でなくても
> 動くことを目標にしています．以下に実装例とテストプログラムを添付します．
... 
> どういう設計にするのがいいのかがよく分からないので，もっと適切な 
> 方法が
> ありましたら教えていただければ幸いです．

私もまさにこの設計がいいと思います．ポイントは，

  * 変換はクラス化し，new (initialize) でパラメターを与えて
    初期化する．
  * 作業領域は内部変数化し，利用者からは隠蔽されるが，不必要な
    メモリー再割り当ては行われない．
  * 変換の実行の際は変換対象だけを引数にすればすむ．

ただし，作業配列の長さは関数によって違い得ます．よって手持ちの作
業配列の長さが足りない場合は，より長くするといった工夫が要るでしょ
う．private method 化するといいのでは？

   def _alloc_wkfft( length )
     if !@wkfft || @wkfft.length < length
	@wkfft = NArray.float( length )
     end
   end
   private _alloc_wkfft

これは，initialize 内でも使えます．

＃ 水田さま：こうするためには，Cラッパーにおいて，作業領域の長さ
   は必要な分「以上」あれば良いという判定基準である必要がありま
   すが，そうなってますか？

以下，細かな点に関しちょっとコメントします．

> require "narray"
> require "ispack_raw"
> 
> class N2Pack
>    def initialize(im, jm, nx, ny)
>      @im, @jm, @nx, @ny = im, jm, nx, ny
>      @itj = NArray.int(5)
>      @iti = NArray.int(5)
>      @tj = NArray.float(2*jm)
>      @ti = NArray.float(2*im)

この上の４行は不要です．（@itj, @tj, @iti, @ti はこの2行下で
上書きされます．中身のサイズは上の通りでしょうが，オブジェクトと
してはまったく別物で，上で作ったのはゴミ集めの対象になるだけ．）

あとで，その変換のサイズを調べられるよう，
@im, @jm, @nx, @ny については値を調べられるようにしたほうが
便利と思います．

   attr_accessor :im, :jm, :nx, :ny

とか．

>      @wkfft = NArray.float(im*jm)
>      @itj, @tj, @iti, @ti = ISPACK::n2init(@jm, @im)
>    end
> 
>    def forward(grid)
>      shape = grid.shape
>      grid.reshape!(grid.total)

この１次元化はしなくても動くようになってないでしょうか．
（もしもなってないなら，Cラッパーのほうを改訂したほうが
いいと思います）

>      obj = ISPACK::n2g2sa(@ny ,@nx, @jm, @im, grid,
>                   @wkfft, @itj, @tj, @iti, @ti)
>      grid.reshape!(*shape)
>      return obj
>    end
> 
>    def backward(spectrum)
>      shape = spectrum.shape
>      spectrum.reshape!(spectrum.total)
>      obj = ISPACK::n2s2ga(@ny ,@nx, @jm, @im, spectrum,
>               @wkfft, @itj, @tj, @iti, @ti)
>      spectrum.reshape!(*shape)
>      return obj
>    end
> 
>    def dx(spectrum)
>      tmp = spectrum.reshape(2*@ny+1, 2*@nx+1)

ここは reshape! でいいでしょう．下でも中身は上書きしてないので．
(tmp = tmp[-1..0, -1..0] では新オブジェクトができます）
なお，もとに戻すのが面倒なら，spectrum.refer.reshape! という手も
あります．

>      factor = -(NArray.float(1,tmp.shape[1]).indgen! - (tmp.shape[1]  
> - 1)/2)
>      tmp = tmp[-1..0, -1..0]
>      return (tmp*factor).reshape!(tmp.total)
>    end
> 
>    def dy(spectrum)
>      tmp = spectrum.reshape(2*@ny+1, 2*@nx+1)
>      factor = -(NArray.float(tmp.shape[0],1).indgen! - (tmp.shape[0]  
> - 1)/2)
>      tmp = tmp[-1..0, -1..0]
>      return (tmp*factor).reshape!(tmp.total)
>    end
> 
> end
> 
> =begin
> ###
> ### test program
> ###
> 
> require "numru/dcl"
> 
> include NumRu
> include NMath
> 
> ## グラフの描画
> def draw(grid)
>    DCL.grfrm
>    DCL.grswnd(0.0, 1.0, 0.0, 1.0)
>    DCL.uspfit
>    DCL.grstrf
> 
>    DCL.ussttl('X', '', 'Y', '')
>    DCL.uelset('ltone', true)
>    DCL.uetone(grid)
>    DCL.usdaxs
>    DCL.udcntr(grid)
> end
> 
> ## 初期設定
> i = 64
> j = 64
> nx = 31
> ny = 31
> x = NArray.float(1,i).indgen!/i
> y = NArray.float(j,1).indgen!/j
> grid = exp(-100.0*((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5)))
> 
> DCL.gropn(1)
> 
> draw(grid)
> ## Fourier  変換 + 演算
> n2pack = N2Pack.new(i, j, nx, ny)
> grid = n2pack.backward(n2pack.dx(n2pack.forward(grid))).
>         reshape!(j, i).transpose(1,0)
> 
> draw(grid)
> 
> DCL.grcls
> =end
> 
> ----
> Yuji Kitamura, Ph.D. <kitamura@xxxxxxxxxxxxxxxxxx>
> Division of Earth and Planetary Sciences,
> Graduate School of Science, Kyoto University
> Tel: +81-75-753-3933/Fax: +81-75-753-3715
>