ベータテスト終了　VfucとEqSysとMatとVec

ベータテスト終了　VfuncとEqSysとMatとVec

妄想中のVec, Mat、VfuncとEqSysの実装は済んで、残っていたベータテストも終わった。
　　（最新仕様は各ページ更新済み）

アルファテストの対象は
　１．MatのLU分解、２．一元一次方程式、３．一元三次方程式、４．三元一次連立方程式
複数の解法で得られた解の比較でテストした。

で、今回は
　５．三元非線形連立方程式（４．の延長）の求解
　　　→例のあれが解ける事が確認できた。
　６．線形単回帰関数(y=ax+b)の観測に基づく係数の推定　（省略）
　７．線形重回帰関数(y=ax+by+cz+d)の観測に基づく係数の推定　（省略）

やっぱ、面倒だった。特に８番（結局、完璧じゃあ無い・・・別記事参照）
　８．非線形関数の観測に基づくパラメータ推定（別記事）

とりあえずパッケージlll.Locに含めた形でベータ公開準備中
ただし、使い易さ重視なので、 メモリ使用効率・計算効率・収束性・頑健性・計算精度・デバッグ効率については一切保障はないことに注意
==== ベータテストドライバ（５）====

//
//  Test for Vfunc, EqSys, Mat, Vec, and their solver methods.
//  Created by Classiclll on 06/03/19.
//

import lll.Loc.*;
public void setup() {
// 5. multi-value nonlinear equation test
  //	<Formulation　for 3D triangulation with two wiiremotes>
  //		 [Pi]：the i-th reference point
  Loc[] pts
  	= new Loc[]{new Loc(.2f,.5f,.1f), new Loc(.12f,.44f,.2f), new Loc(.6f,.3f,.4f)};
  //		 |[Pj-Pi]|=Lij, (i!=j)
  double[] l2ij
  	= new double[]{pts[0].dist2(pts[1]), pts[0].dist2(pts[2]), pts[1].dist2(pts[2])};
  //		　[di]：the directional vector to [Pi]
  Loc[] dir = new Loc[]{pts[0].unit(), pts[1].unit(), pts[2].unit()};
  //		then each reference point is on the line
  //		     [Pi]=ti*[di]
  //		let cij=([dj]*[di]), then we get the followings
  double[] cij
  	= new double[]{dir[0].dot(dir[1]), dir[0].dot(dir[2]), dir[1].dot(dir[2])};
  //		　tj^2 - 2*cij*tj*ti + ti^2 - Lij^2 = 0
  Triangular trig = new Triangular(l2ij,cij);
  Vec x0 = new Vec(new double[]{1,1,1});
  Vec root;

  println("[ nonlinear equation system Test ]");
  println(">>by Newton\n"+ trig.jacobAt(x0));
  println(">>by Newton\n"+ (root=trig.solveByNewton(x0)));
  println(">>   F(x) = "+ trig.valueAt(root));
  println(">>by Simplex\n"+ (root=trig.solveBySimplex(x0,1000)));
  println(">>   F(x) = "+ trig.valueAt(root));

テスト結果

[ nonlinear equation system Test ]
>>by Newton
Mat{{0.06425809860229492, 0.06425809860229492, 0.0},
    {0.0, 0.5506737232208252, 0.5506737232208252},
    {0.5396480560302734, 0.0, 0.5396480560302734}}
>>by Newton
Vec(0.5579702401804243 ,0.5578155284121784 ,0.02707795585202451)
>>   F(x) = Vec(6.743494651573201E-13 ,3.0631813197068425E-9 ,2.9968867942820054E-9)
>>by Simplex
Vec(0.4979142256587791 ,0.5474806070350684 ,0.7809789539962118)
>>   F(x) = Vec(-2.6521510391508407E-5 ,-2.6521527641598652E-5 ,-2.652153279614211E-5)

Vfunc、EqSysのテスト用実装クラス

import lll.Loc.*;
//model of equation, tj^2 - 2*cij*tj*ti + ti^2 - Lij^2 = 0
public class Triangular extends EqSys {
  private double[] l2ij;
  private double[] cij;
  public Triangular(double[] l2, double[] c) {
	super(3,3);
	l2ij = l2;
	cij = c;
  }
  public Vec valueAt(Vec t) {	// fij(t) = tj^2 - 2*cij*tj*ti + ti^2 - Lij^2
  return new Vec(new double[]{
	t.elem(0) * (t.elem(0)-2*cij[0]*t.elem(1)) + t.elem(1)*t.elem(1) - l2ij[0],
	t.elem(1) * (t.elem(1)-2*cij[1]*t.elem(2)) + t.elem(2)*t.elem(2) - l2ij[1],
	t.elem(2) * (t.elem(2)-2*cij[2]*t.elem(0)) + t.elem(0)*t.elem(0) - l2ij[2]});
  }
  public Mat jacobAt(Vec t) {	// Jacobian of f(X)={f0(X), f1(X), f2(X)}
  double[][] jacobi =	{
// { {		dxf0(X),	dyf0(X),	dzf0(X)}
	{2*t.elem(0)-2*cij[0]*t.elem(1), -2*cij[0]*t.elem(0)+2*t.elem(1), 0 },
//   {		dxf1(X),	dyf1(X),	dzf1(X)}
	{ 0,	2*t.elem(1)-2*cij[1]*t.elem(2),    -2*cij[1]*t.elem(1)+2*t.elem(2)},
//   {		dxf2(X),	dyf2(X), 	dzf2(X)} }
	{ -2*cij[2]*t.elem(2)+2*t.elem(0), 0, 2*t.elem(2)-2*cij[2]*t.elem(0)    } };
	return new Mat(jacobi);
 }
 public Mat gradParamAt(Vec x) {
	return Mat.NaN.copy();
 }
}