数値計算界に伝えたいオブジェクト指向構文とその用法
数値計算界に伝えたいオブジェクト指向構文とその用法:
これが
このように分解されている。 xi=xipとループ中に直接xi+=...という記述順序から 以下の式がパッと想像できるだろうか?
Javascriptの場合は、代入文ではなく代入式なので、次のように短縮してもよい。
動機
- 数値計算する人のコードは、与式と計算式にギャップがある
- そのため、理論における与式が実装コードのどこで計算されているか理解に時間がかかる
- そこで、クラス変数やインスタンスメソッドなどオブジェクト指向構文を用いて、与式と計算式のギャップを小さくすることを示す
事例: 緯度経度と地図座標
- GeoLocation APIで得られた位置は地図のどこ?
- 地表は楕円体曲面であり、楕円体で定義される座標である緯度経度(lat, lng)が、地図上の平面直角座標において具体的に縦x横yどういう値になるのか?
-
注: 地図においては北がx軸で東がy軸
- 理論的な与式: 国土地理院/測量計算サイト/平面直角座標への換算
- 実装コード例:[同サイト引用論文/Gauss-Krüger投影における経緯度座標及び平面直角座標相互間の座標換算についてのより簡明な計算方法](http://www.gsi.go.jp/common/000061216.pdf）
function sinh(x) { return 0.5*(Math.exp(x)-Math.exp(-x)) }
function cosh(x) { return 0.5*(Math.exp(x)+Math.exp(-x)) }
function arctanh(x) { return 0.5*Math.log((1+x)/(1-x)) }
a=6378137 ; rf=298.257222101 ; m0=0.9999 ; s2r=Math.PI/648000 ; n=0.5/(rf-0.5)
n15=1.5*n ; anh=0.5*a/(1+n) ; nsq=n*n
e2n=2*Math.sqrt(n)/(1+n) ; ra=2*anh*m0*(1+nsq/4+nsq*nsq/64)
jt=5 ; jt2=2*jt ; ep=1.0 ; e=[] ; s=[0.0] ; t=[] ; alp=[]
for(k=1; k<=jt; k++) { ep*=e[k]=n15/k-n ; e[k+jt]=n15/(k+jt)-n }
// 展開パラメータの事前入力
alp[1]=(1/2+(-2/3+(5/16+(41/180-127/288*n)*n)*n)*n)*n
alp[2]=(13/48+(-3/5+(557/1440+281/630*n)*n)*n)*nsq
alp[3]=(61/240+(-103/140+15061/26880*n)*n)*n*nsq
alp[4]=(49561/161280-179/168*n)*nsq*nsq
alp[5]=34729/80640*n*nsq*nsq
// 平面直角座標の座標系原点の緯度を度単位で、経度を分単位で格納
phi0=[0,33,33,36,33,36,36,36,36,36,40,44,44,44,26,26,26,26,20,26]
lmbd0=[0,7770,7860,7930,8010,8060,8160,8230,8310,8390,8450,8415,8535,8655,8520,7
650,7440,7860,8160,9240]
// 該当緯度の 2 倍角の入力により赤道からの子午線弧長を求める関数
function Merid(phi2) {
dc=2.0*Math.cos(phi2) ; s[1]=Math.sin(phi2)
for(i=1; i<=jt2; i++) { s[i+1]=dc*s[i]-s[i-1] ; t[i]=(1.0/i-4.0*i)*s[i] }
sum=0.0 ; c1=ep ; j=jt
while(j) {
c2=phi2 ; c3=2.0 ; l=j ; m=0
while(l) { c2+=(c3/=e[l--])*t[++m]+(c3*=e[2*j-l])*t[++m] }
sum+=c1*c1*c2 ; c1/=e[j--]
}
return anh*(sum+phi2)
}
// 与件入力
num=eval(prompt("座標系番号を入力してください。"))
phi=eval(prompt("緯度を ddmmss.ssss 形式で入力してください。"))
lmbd=eval(prompt("経度を dddmmss.ssss 形式で入力してください。"))
phideg=Math.floor(phi/10000) ; phimin=Math.floor((phi-phideg*10000)/100)
phirad=(phideg*3600+phimin*60+phi-phideg*10000-phimin*100)*s2r
lmbddeg=Math.floor(lmbd/10000) ; lmbdmin=Math.floor((lmbd-lmbddeg*10000)/100)
lmbdsec=lmbddeg*3600+lmbdmin*60+lmbd-lmbddeg*10000-lmbdmin*100
// 実際の計算実行部分
sphi=Math.sin(phirad) ; nphi=(1-n)/(1+n)*Math.tan(phirad)
dlmbd=(lmbdsec-lmbd0[num]*60)*s2r
sdlmbd=Math.sin(dlmbd) ; cdlmbd=Math.cos(dlmbd)
tchi=sinh(arctanh(sphi)-e2n*arctanh(e2n*sphi)) ; cchi=Math.sqrt(1+tchi*tchi)
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)
}
x=ra*xi-m0*Merid(2*phi0[num]*3600*s2r) ; y=ra*eta
gmm=Math.atan2(tau*cchi*cdlmbd+sgm*tchi*sdlmbd, sgm*cchi*cdlmbd-tau*tchi*sdlmbd)
m=ra/a*Math.sqrt((sgm*sgm+tau*tau)/(tchi*tchi+cdlmbd*cdlmbd)*(1+nphi*nphi))
// ラジアン → 度分秒変換
sgn=(gmm<0)
gdo=Math.floor(gmm/s2r/3600)+sgn
gfun=Math.floor((gmm/s2r-gdo*3600)/60)+sgn
gbyou=gmm/s2r-gdo*3600-gfun*60
// 結果表示
document.write("<h2>座標系番号： " + num + " 緯度： " + phi + " 経度： " + lmbd +
"<br/><br/>")
document.write("Ｘ＝" + x + "，Ｙ＝" + y + "<br/>")
document.write("γ＝" + (sgn?"－":"＋") + Math.abs(gdo) + "°" + Math.abs(gfun) + "′" + Math.abs(gbyou) + "″，m＝" + m + "<br/></h2>")
オブジェクト指向でリファクタリング
- 一時変数ではなく、与式ごとにメソッド定義
- 複数式から参照される場合は、再計算を回避
- for文の代わりに、コレクションメソッド
- 定数パラメータはクラス定数
class GeoLocation {
// 1. 定数パラメータはクラス定数
// 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}
static get ra(){
const nsq = GeoLocation.nsq;
return GeoLocation.m0*(1+nsq/4+nsq*nsq/64)*GeoLocation.a/(1+GeoLocation.n);
}
get phi(){ // 2. 一時変数ではなく、与式ごとにメソッド定義
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 sphi(){
return Math.sin(this.phi.rad);
}
get nphi(){
const n = GeoLocation.n;
return (1-n)/(1+n)*Math.tan(this.phi.rad);
}
get dlmbd(){
if (this._dlmbd == null) // 3. 複数式から参照される場合は、再計算を回避
this._dlmbd = this.lmbd.rad - this.lmbd0(this.num);
return this._dlmbd;
}
get sdlmbd(){
if (this._sdlmbd == null)
this._sdlmbd = Math.sin(this.dlmbd);
return this._sdlmbd;
}
get cdlmbd(){
if (this._cdlmbd == null)
this._cdlmbd = Math.cos(this.dlmbd);
return this._cdlmbd;
}
get tchi(){
if (this._tchi == null){
const n = GeoLocation.n;
const e2n = 2*Math.sqrt(n)/(1+n);
const sphi = this.sphi;
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 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 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 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 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 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))
}
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));
}
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;
}
}
}
元コードを参照しつつ解説
- 一時変数ではなく、与式ごとにメソッド定義
- 複数式から参照される場合は、再計算を回避
- for文の代わりに、コレクションメソッド
- 定数パラメータはクラス定数
1. 一時変数ではなく、与式ごとにメソッドを定義すべし
- 理由: 理論上の与式と実装コードの対応を見つけやすくするため
- 元コードをアプリ界からみると、変数(xi, eta, sgm, tau)の使い方に違和感
- 考えないとわからない
// 実際の計算実行部分
// ...中略
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'
\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;
}
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がついている
4. 定数パラメータはクラス定数
- 理由: インスタンスオブジェクト毎に定数値を共通に持つ必要がないから
- JavaScriptの場合は、クラス定数の構文がないので、static get で代用
まとめ
- 数値計算する人のコードは、与式と計算式にギャップがある
- これは、数式はトップダウン記述であるのに対して
- 実装コードである計算式は、処理順序記述になっており
- 理解の順序と記述順序が真逆になっているためである
- そのため、理論における与式が実装コードのどこで計算されているか理解に時間がかかっていた
- そこで、クラス変数やインスタンスメソッドなどオブジェクト指向構文を用いて、与式と計算式のギャップを小さくすることを示した
- 具体的には次の4点を施した
- 一時変数ではなく、与式ごとにメソッド定義
- 複数式から参照される場合は、再計算を回避
- for文の代わりに、コレクションメソッド
- 定数パラメータはクラス定数
発展的な課題
- 添字が1始まりの配列を表現する、文法要素、があるともう少し記述が自然になる
- 値がnullなら値を代入してしかも、再設定できることを示す、演算子
- オブジェクト指向的には、シングルトンパターンと類似したイメージ
- 変数には、未定義状態と値が設定された状態のどちらかしかない、erlangのような変数の考え方なら、必要ない
// SwiftのNil Coalescing Operatorに似たようなもの
// 使い方はイメージこんな感じで記述できるといいかな
get xip(){
return this._xip ??= Math.atan2(this.tchi, this.cdlmbd);
}
コメント
コメントを投稿