数値計算界に伝えたいオブジェクト指向構文とその用法

動機

• 数値計算する人のコードは、与式と実装式にギャップがある
• そのため、理論における与式が実装コードのどこで計算されているか理解に時間がかかる
• そこで、クラス変数やインスタンスメソッドなどオブジェクト指向構文を用いて、与式と計算式のギャップを小さくすることを示す

---

事例: 緯度経度と地図座標

• GeoLocation APIで得られた位置は地図のどこ?
  • 地表は楕円体曲面であり、楕円体で定義される座標である緯度経度(lat, lng)が、地図上の平面直角座標において具体的に縦x横yどういう値になるのか?
  • 注: 地図においては北がx軸で東がy軸
• 理論的な与式: 国土地理院/測量計算サイト/平面直角座標への換算...(A)

• 実装コード例: 同サイト引用論文/Gauss-Krüger投影における経緯度座標及び平面直角座標相互間の座標換算についてのより簡明な計算方法 ...(B)

罠

• 理論的な与式(A)と、引用論文の実装コード(B)が元にしている理論的な与式とが、一致してないので要注意
• (B)の理論式を５次近似した式が(A)
• (B)の実装コードは、5次近似する前の理論式を実装コードに落としてそのまま計算している

---

そもそも数値計算におけるプログラミング技法とは

• 基本
  • 各式共通に使い回す値を一時変数に格納する
  • n^5など、累乗の計算は、nnnnnと素直に計算するより、nsq=nn, nsqnsq*nとした方が、乗算回数が4回ではなく3回と回数が減る
    • n^2, n^4, n^8 などを用意しておく
  • 級数は小さい値から大きな値へ足していく
    • 大きな値から小さな値を足すと、計算誤差が大きくなるため
• 計算誤差が少なさと、解までの収束の速さが肝

---

オブジェクト指向でリファクタリング

1. 与式は一時変数ではなく、メソッドで定義
2. 複数式から参照される場合は、再計算を回避
3. for文の代わりに、コレクションメソッド
4. 定数パラメータはクラス定数

class GeoLocation {
  // 4. 定数パラメータはクラス定数
  // javascriptではクラス定数構文がないのでstaticで代用
  static get a() {return 6378137 }
  static get RF(){return 298.257222101}
  static get m0(){return 0.9999 }
  static get n() {return this._n || (this._n = 1/(2*GeoLocation.RF -1));} //参照ごとに何度も計算する無駄を省く
  static get nsq(){return GeoLocation.n * GeoLocation.n}
  
  get x(){ return GeoLocation.ra * this.xi - this.rs_phi0}
  get y(){ return GeoLocation.ra * this.eta};
  get gmm(){
    const tau = this.tau, cchi = this.cchi, cdlmbd = this.cdlmbd;
    const sgm = this.sgm, tchi = this.tchi, sdlmbd = this.sdlmbd;
    return Math.atan2(tau*cchi*cdlmbd+sgm*tchi*sdlmbd, sgm*cchi*cdlmbd-tau*tchi*sdlmbd)
  }
  get m(){
    const sgm = this.sgm, tau = this.tau, tchi = this.tchi, cdlmbd = this.cdlmbd, nphi = this.nphi;
    return GeoLocation.ra/GeoLocation.a * Math.sqrt((sgm*sgm+tau*tau)/(tchi*tchi+cdlmbd*cdlmbd)*(1+nphi*nphi));
  }  

  get xi(){
    return this.alpha.inject(this.xip, (prev, alp_j, j)=>{
      return j < 1 ? prev: (prev + alp_j.sin * cosh(2*j*this.etap))
    })
  }
  get eta(){
    return this.alpha.inject(this.etap, (prev, alp_j, j)=>{
      return j < 1 ? prev: (prev + alp_j.cos * sinh(2*j*this.etap))
    })
  }
  get xip(){
    if (this._xip == null)
      this._xip = Math.atan2(this.tchi, this.cdlmbd);
    return this._xip;
  }
  get etap(){
    if (this._etap == null)
      this._etap = arctanh(this.sdlmbd/this.cchi);
    return this._etap;
  }

  get phi(){ // 1. 与式は一時変数ではなく、メソッドで定義
    if (this._phi == null){
      this._phi = deg2rad(this.coords.latitude);
    }
    return this._phi;
  }
  get lmbd(){
    return deg2rad(this.coords.longitude);
  }
  phi0(num=2){
    // 平面直角座標の座標系原点の緯度を度単位で、経度を分単位で格納                                                              
    const phi0 =[0,33,33,36,33,36,36,36,36,36,40,44,44,44,26,26,26,26,20,26] // 1-19                                             
    return deg2rad(phi0[num]).rad;
  }
  lmbd0(num=2){
    // 平面直角座標の座標系原点の緯度を度単位で、経度を分単位で格納                                                              
    const lmbd0=[0,7770,7860,7930,8010,8060,8160,8230,8310,8390,8450,8415,   // 1-19                                             
                 8535,8655,8520,7650,7440,7860,8160,9240]
    return deg2rad(lmbd0[num]/60).rad
  }

  get tchi(){
    if (this._tchi == null){
      const n    = GeoLocation.n;
      const e2n  = 2*Math.sqrt(n)/(1+n);
      const sphi = Math.sin(this.phi.rad);
      this._tchi = sinh(arctanh(sphi) - e2n*arctanh(e2n*sphi));
    }
    return this._tchi;
  }
  get cchi(){
    if (this._cchi == null)
      this._cchi = Math.sqrt(1+this.tchi*this.tchi);
    return this._cchi;
  }
  get nphi(){
    const n = GeoLocation.n;
    return (1-n)/(1+n)*Math.tan(this.phi.rad);
  }

  get cdlmbd(){
    if (this._cdlmbd == null)
      this._cdlmbd = Math.cos(this.dlmbd);
    return this._cdlmbd;
  }
  get sdlmbd(){
    if (this._sdlmbd == null)
      this._sdlmbd = Math.sin(this.dlmbd);
    return this._sdlmbd;
  }
  get dlmbd(){
    if (this._dlmbd == null) // 2. 複数式から参照される場合は、再計算を回避
      this._dlmbd = this.lmbd.rad - this.lmbd0(this.num);
    return this._dlmbd;
  }

  get sgm(){
    if (this._sgm == null)
      this._sgm = this.alpha.inject(1, (prev, alp_j, j)=>{
        return j < 1 ? prev: (prev + 2*j*alp_j.cos*cosh(2*j*this.etap))
      })
    return this._sgm;
  }
  get tau(){
    if (this._tau == null)
      this._tau = this.alpha.inject(0, (prev, alp_j, j)=>{
        return j < 1 ? prev: (prev + 2*j*alp_j.sin*sinh(2*j*this.etap))
      })
    return this._tau;
  }

  get alpha(){ 
    if (this._alp == null){
      const n = GeoLocation.n, nsq = GeoLocation.nsq;
      this._alp = [null,
            (1/2+(-2/3+(5/16+(41/180-127/288*n)*n)*n)*n)*n,
            (13/48+(-3/5+(557/1440+281/630*n)*n)*n)*nsq,
            (61/240+(-103/140+15061/26880*n)*n)*n*nsq,
            (49561/161280-179/168*n)*nsq*nsq,
            34729/80640*n*nsq*nsq
           ].map((e_j, j)=>{ // 3. for文の代わりにコレクションメソッド
             return j < 1 ? e_j: {sin: e_j*Math.sin(2*j*this.xip), cos: e_j*Math.cos(2*j*this.xip)};
           })
    }
    return this._alp;
  }
  get rs_phi0(){   
    const a   = GeoLocation.a; 
    const m0  = GeoLocation.m0;
    const n   = GeoLocation.n;
    const nsq = GeoLocation.nsq;
    const phi0= this.phi0(this.num);
    return [null,
            -3/2*(1-nsq/8-nsq*nsq/64)*n,
            15/16*(1-nsq/4)*nsq,
            -35/48*(1-5*nsq/16)*n*nsq,
            315/512*nsq*nsq,
            -693/1280*n*nsq*nsq
           ].inject((1+nsq/4+nsq*nsq/64)*phi0, (prev, ra_i, i)=>{
             return i < 1 ? prev : (prev + ra_i * Math.sin(2*i*phi0));
           })*(m0*a/(1+n))
  }
  static get ra(){
    const nsq = GeoLocation.nsq;
    return GeoLocation.m0*(1+nsq/4+nsq*nsq/64)*GeoLocation.a/(1+GeoLocation.n);
  }

  async getCurrentPosition(opts){
      try {
        const api = navigator.geolocation || google.gears.factory.create('beta.geolocation');
        api.getCurrentPosition(this.success, this.failure, opts);
      }catch(err){
        this.failure({code:-1, message:"ご使用のブラウザはGeoLocation APIに対応していないため、位置情報が得られません"})
      }
  }

  success(pos){
    this.coords = pos.coords;
    ["phi", "dlmbd", "sdlmbd", "cdlmbd", "tchi", "cchi", "xip", "etap", "alpha", "sgm", "tau"].forEach((k)=>{
      this["_"+k] = null;
    })
    console.log(this.x, this.y);
  }
  failure(err){
    const map = document.getElementById('map');
    switch(err.code){
    case 1:
      map.innerHTML= "位置情報の利用が許可されていません";
      break;
    case 2:
      map.innerHTML= "デバイスの位置が判定できません";
      break;
    case 3:
      map.innerHTML= "タイムアウトしました";
    default:
      map.innerHTML= err.message;
    }
  }


}

---

元コードを参照しつつ解説

1. 与式は一時変数ではなく、メソッドで定義
2. 複数式から参照される場合は、再計算を回避
3. for文の代わりに、コレクションメソッド
4. 定数パラメータはクラス定数

---

1. 与式は一時変数ではなく、メソッドで定義

• 理由: 理論上の与式と実装コードの対応を見つけやすくするため
• 元コードをアプリ界からみると、変数(xi, eta, sgm, tau)の使い方に違和感
  • 共通の係数alsinとalcosがあるからって、まとめすぎ
  • 与式と実装コードにギャップが生じていて、考えないと元式がわからない

元コード

// 実際の計算実行部分
// ...中略
xi=xip=Math.atan2(tchi, cdlmbd); 
eta=etap=arctanh(sdlmbd/cchi); 
sgm=1;
tau=0;
for(j=alp.length; --j; ) {
  alsin=alp[j]*Math.sin(2*j*xip) ; alcos=alp[j]*Math.cos(2*j*xip);
  xi+=alsin*cosh(2*j*etap);
  eta+=alcos*sinh(2*j*etap);
  sgm+=2*j*alcos*cosh(2*j*etap);
  tau+=2*j*alsin*sinh(2*j*etap);
}

例えば変数xiに注目する

• xi は理論式に出てこない、何を意味するかよくわからない変数(に見える)
  • (1) 唐突に出てくる 変数xi
  • (2) xiはループで足されている
  • (3) ここまできてようやく意味がわかる

// xi に関係ある文だけを抽出
xi=xip=Math.atan2(tchi, cdlmbd);    // ...(1)
for(j=alp.length; --j;){
 alsin=alp[j]*Math.sin(2*j*xip);
 xi+=alsin*cosh(2*j*etap);          // ...(2)
}
x= ra * xi - rs_phi0;               //...(3)

---

理論式を理解する流れが、計算の順序と真逆

• 理論式は宣言的な記述
• 実装コードは計算の順序であって手続き的な記述になっている
• 記述順序が真逆なのが、理解をしにくさの最大の要因

x = \overline{A}(\xi'+\sum_{j=1}^5 \alpha_j \sin 2j\xi' \cosh2j\eta') - \overline{S}_{{\varphi}_0}

これが

x = \overline{A}\xi- \overline{S}_{{\varphi}_0} \\

\xi = \xi'+\sum_{j=1}^5 \alpha_j \sin 2j\xi' \cosh2j\eta' 

このように分解されている。　

---

対策: メソッド化

• 計算順序を変えずに記述順序を変える
• 理解の流れ＝記述の順序にできる
• 見ればわかる
• ループで回さずにコレクションメソッドを使用することで、第1項がxipであることも明示的にわかる

  get x(){ 
    return GeoLocation.ra * this.xi - this.rs_phi0
  }

  get xi(){
    return this.alpha.inject(this.xip, (prev, alp_j, j)=>{
      return j < 1 ? prev: (prev + alp_j.sin * cosh(2*j*this.etap))
    })
  }  

  get alpha(){ 
    if (this._alp == null){
      const n = GeoLocation.n, nsq = GeoLocation.nsq;
      this._alp = [null,
            (1/2+(-2/3+(5/16+(41/180-127/288*n)*n)*n)*n)*n,
            (13/48+(-3/5+(557/1440+281/630*n)*n)*n)*nsq,
            (61/240+(-103/140+15061/26880*n)*n)*n*nsq,
            (49561/161280-179/168*n)*nsq*nsq,
            34729/80640*n*nsq*nsq
           ].map((e_j, j)=>{ // 4. for文の代わりにコレクションメソッド
             return j < 1 ? e_j: {sin: e_j*Math.sin(2*j*this.xip), cos: e_j*Math.cos(2*j*this.xip)};
           })
    }
    return this._alp;
  }

  get xip(){
    if (this._xip == null)
      this._xip = Math.atan2(this.tchi, this.cdlmbd);
    return this._xip;
  }
  get etap(){
    if (this._etap == null)
      this._etap = arctanh(this.sdlmbd/this.cchi);
    return this._etap;
  }

---

2. 複数式から参照される場合は、再計算を回避

• 理由: メソッド化することによって、複数式からコールされると、呼ばれた回数分同じ計算をしてしまう
• 対策(1): ガード条件を入れる

  get xip(){
    if (this._xip == null)
      this._xip = Math.atan2(this.tchi, this.cdlmbd);
    return this._xip;
  }

Javascriptの場合は、代入文ではなく代入式なので、次のように短縮してもよい。

  get xip(){
    return this._xip || (this._xip = Math.atan2(this.tchi, this.cdlmbd);
  }
  // this._xip がnullなら、代入式を評価する　　　　
  // 論理演算の結合規則にしたがって、代入式は括弧が必要

• 対策(2): 変数の初期化も忘れずに
  • 新しい位置情報が得られたら、一旦、null値で初期化

  success(pos){
    this.coords = pos.coords;
    ["phi", "dlmbd", "sdlmbd", "cdlmbd", "tchi", "cchi", "xip", "etap", "alpha", "sgm", "tau"].forEach((k)=>{
      this["_"+k] = null;
    })
    console.log(this.x, this.y);
  }

---

3. for文の代わりに、コレクションメソッド

• 理由: 処理を繰り返すのではなく、値を集めているから。意味論として
  • コレクションメソッドを使うのがループ文よりも適切と思われる

  get xi(){
    return this.alpha.inject(this.xip, (prev, alp_j, j)=>{
      return j < 1 ? prev: (prev + alp_j.sin * cosh(2*j*this.etap))
    })
  }
  Array.prototype.inject = function(init, cb){
    return this.reduce(cb, init);
  }  

• このケースの場合は、j=0という添字に対応する係数alp[0]が存在しないので、条件j<1がついている
• なお、FORTRANやCOBOLやBASIC（可変）の配列は添字が1始まり

---

4. 定数パラメータはクラス定数

• 理由: インスタンスオブジェクト毎に定数値を共通に持つ必要がないから
  • JavaScriptの場合は、クラス定数の構文がないので、static get で代用

---

#まとめ

• 数値計算する人のコードは、与式と計算式にギャップがある
  • これは、数式はトップダウンの宣言的な記述であるのに対して
  • 実装コードである計算式は、処理順序を指定する手続き的な記述になっており
  • 理解の順序と記述順序が真逆になっているためである
• そのため、理論における与式が実装コードのどこで計算されているか理解に時間がかかっていた
• そこで、クラス変数やインスタンスメソッドなどオブジェクト指向構文を用いて、与式と計算式のギャップを小さくすることを示した
  1. 与式ごとに一時変数ではなく、メソッドで定義
  2. 複数式から参照される場合は、再計算を回避
  3. for文の代わりに、コレクションメソッド
  4. 定数パラメータはクラス定数

発展的な課題

• 添字が1始まりの配列を表現する、文法要素、があるともう少し記述が自然になる
• 値がnullなら値を代入してしかも、再設定できることを示す、演算子
  • オブジェクト指向的には、シングルトンパターンと類似したイメージ
  • 変数には、未定義状態と値が設定された状態のどちらかしかない、erlangのような変数の考え方なら、必要ない

  // SwiftのNil Coalescing Operatorに似たようなもの
  // 使い方はイメージこんな感じで記述できるといいかな
  get xip(){
    return this._xip ??= Math.atan2(this.tchi, this.cdlmbd);
  }