/*
 * Copyright (c) 2009 The openGion Project.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND,
 * either express or implied. See the License for the specific language
 * governing permissions and limitations under the License.
 */
package org.opengion.penguin.math.statistics;

import java.util.Arrays;

/**
 * 独自実装の二次回帰計算クラスです。
 * f(x) = c1x^2 + c2x + c3
 * の曲線を求めます。
 */
public class HybsSquadraticRegression implements HybsRegression {
	private final double[] cnst = new double[3] ;		// 係数(０次、１次、２次)
	private double rsquare;		// 決定係数 今のところ求めていない

	/**
	 * コンストラクタ。
	 * 与えた二次元データを元に二次回帰を計算します。
	 *
	 * @param data xとyの組み合わせの配列
	 */
	public HybsSquadraticRegression( final double[][] data ) {
		//二次回帰曲線を求めるが、これはapacheにはなさそうなので自前で計算する。
		train( data );
	}

	/**
	 * 係数計算
	 *
	 *	c3Σ＋c2Σx＋c1Σx^2＝Σy
	 *	c3Σx＋c2Σ(x^2)＋c1Σx^3＝Σ(xy)
	 *	c3Σ(x^2)＋c2Σ(x^3)＋c1Σ(x^4)＝Σ(ｘ^2*y)
	 *	この三元連立方程式を解くことになる。
	 *
	 * @param data x,yの配列
	 */
	private void train( final double[][] data ) {
		// xの二乗等の総和用
		final int data_n=data.length;
		double sumx2	= 0;
		double sumx		= 0;
		double sumxy	= 0;
		double sumy		= 0;
		double sumx3	= 0;
		double sumx2y	= 0;
		double sumx4	= 0;

		// まずは計算に使うための和を計算
		for( int i=0; i < data_n; i++ ) {
			final double data_x	= data[i][0];
			final double data_y	= data[i][1];
			final double x2		= data_x*data_x;

			sumx	+= data_x;
			sumx2	+= x2;
			sumxy	+= data_x * data_y;
			sumy	+= data_y;
			sumx3	+= x2 * data_x;
			sumx2y	+= x2 * data_y;
			sumx4	+= x2 * x2;
		}

		// ガウス・ジョルダン法で係数計算
		final double diffx2 = sumx2 - sumx * sumx / data_n;
		final double diffxy = sumxy - sumx * sumy / data_n;
		final double diffx3 = sumx3 - sumx2 * sumx /data_n;
		final double diffx2y = sumx2y - sumx2 * sumy /data_n;
		final double diffx4 = sumx4 - sumx2 * sumx2 /data_n;
		final double diffd = diffx2 * diffx4 - diffx3 * diffx3;

		cnst[2] = ( diffx2y * diffx2 - diffxy * diffx3 ) / diffd;
		cnst[1] = ( diffxy * diffx4 - diffx2y * diffx3 ) / diffd;
		cnst[0] = sumy/data_n - cnst[1]*sumx/ data_n - cnst[2]*sumx2/data_n;

		rsquare = 0;		// 決定係数 今のところ求めていない
	}

	/**
	 * 係数(０次、１次、２次)の順にセットした配列を返します。
	 *
	 * @return 係数の配列
	 */
	@Override	// HybsRegression
	public double[] getCoefficient() {
		return Arrays.copyOf( cnst,cnst.length );
	}

	/**
	 * 決定係数の取得。
	 * @return 決定係数
	 */
	@Override	// HybsRegression
	public double getRSquare() {
		return rsquare;
	}

	/**
	 * c2*x^2 + c1*x + c0を計算。
	 *
	 * @param in_x 必要な大きさの変数配列
	 * @return 計算結果
	 */
	@Override	// HybsRegression
	public double predict( final double... in_x ) {
		return cnst[2] * in_x[0] * in_x[0] + cnst[1] * in_x[0] + cnst[0];
	}

	// ================ ここまでが本体 ================
	/**
	 * ここからテスト用mainメソッド 。
	 *
	 * @param args 引数
	 */
	public static void main( final String[] args ) {
		final double[][] data = {{1, 2.3}, {2, 5.1}, {3, 9.1}, {4, 16.2}};

		final HybsSquadraticRegression sr = new HybsSquadraticRegression(data);

		final double[] cnst = sr.getCoefficient();

		System.out.println(cnst[2]);
		System.out.println(cnst[1]);
		System.out.println(cnst[0]);

		System.out.println(sr.predict( 5 ));
	}
}

