/*
 * 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の二乗等の総和用
                // 8.5.5.1 (2024/02/29) PMD 7.0.0 LocalVariableNamingConventions
//              final int data_n=data.length;
                final int count=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++ ) {
                for( int i=0; i < count; i++ ) {
                        // 8.5.5.1 (2024/02/29) PMD 7.0.0 LocalVariableNamingConventions
//                      final double data_x     = data[i][0];
//                      final double data_y     = data[i][1];
                        final double dataX      = data[i][0];
                        final double dataY      = data[i][1];
//                      final double x2         = data_x*data_x;
                        final double x2         = dataX*dataX;

//                      sumx    += data_x;
                        sumx    += dataX;
                        sumx2   += x2;
//                      sumxy   += data_x * data_y;
                        sumxy   += dataX * dataY;
                        sumy    += dataY;
//                      sumx3   += x2 * data_x;
//                      sumx2y  += x2 * data_y;
                        sumx3   += x2 * dataX;
                        sumx2y  += x2 * dataY;
                        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 diffx2 = sumx2 - sumx * sumx / count;
                final double diffxy = sumxy - sumx * sumy / count;
                final double diffx3 = sumx3 - sumx2 * sumx /count;
                final double diffx2y = sumx2y - sumx2 * sumy /count;
                final double diffx4 = sumx4 - sumx2 * sumx2 /count;
                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;
                cnst[0] = sumy/count - cnst[1]*sumx/count - cnst[2]*sumx2/count;

                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 ));
        }
}