/*
 * 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;
import java.util.Collections;
import java.util.List;

/**
 * 多項ロジスティック回帰の実装です。
 * 確率的勾配降下法(SGD)を利用します。
 *
 * ロジスティック回帰はn次元の情報からどのグループに所属するかの予測値を得るための手法の一つです。
 *
 * 実装は
 * http://nbviewer.jupyter.org/gist/mitmul/9283713
 * https://yusugomori.com/projects/deep-learning/
 * を参考にしています。
 */
public class HybsLogisticRegression {
//      private final int n_N;                  // データ個数
        /** データ個数 */
        private final int n_cnt;                // 8.5.4.2 (2024/01/12) PMD 7.0.0 FieldNamingConventions 対応 n_N ⇒ n_cnt

        /** データ次元 */
        private final int n_in;
        /** ラベル種別数 */
        private final int n_out;

        // 8.5.4.2 (2024/01/12) PMD 7.0.0 ImmutableField 対応
        /** 写像変数ベクトル f(x) = Wx + b */
        private final double[][] vW;
        private final double[] vb;

        /**
         * コンストラクタ。
         *
         * 学習もしてしまう。
         *
         * xはデータセット各行がn次元の説明変数となっている。
         * trainはそれに対する{0,1,0},{1,0,0}のようなラベルを示すベクトルとなる。
         * 学習率は通常、0.1程度を設定する。
         * このロジックではループ毎に0.95をかけて徐々に学習率が下がるようにしている。
         * 全データを利用すると時間がかかる場合があるので、確率的勾配降下法を利用しているが、
         * 選択個数はデータに対する割合を与える。
         * データ個数が少ない場合は1をセットすればよい。
         *
         * @og.rev 8.5.4.2 (2024/01/12) PMD 7.0.0 FieldNamingConventions 対応 n_N ⇒ n_cnt
         *
         * @param data データセット配列
         * @param label データに対応したラベルを示す配列
         * @param learning_rate 学習係数(0から1の間の数値)
         * @param loop 学習のループ回数（ミニバッチを作る回数）
         * @param minibatch_rate 全体に対するミニバッチの割合(0から1の間の数値)
         *
         */
        public HybsLogisticRegression(final double data[][], final int label[][], final double learning_rate ,final int loop, final double minibatch_rate ) {
        //      List<Integer> indexList; //シャッフル用

//              this.n_N = data.length;
                this.n_cnt = data.length;                                       // 8.5.4.2 (2024/01/12) PMD 7.0.0
                this.n_in = data[0].length;
                this.n_out = label[0].length;                           // ラベル種別

                vW = new double[n_out][n_in];
                vb = new double[n_out];

                // 確率勾配に利用するための配列インデックス配列
//              final Integer[] random_index = new Integer[n_N]; //プリミティブ型だとasListできないため
                // 8.5.5.1 (2024/02/29) PMD 7.0.0 LocalVariableNamingConventions
//              final Integer[] random_index = new Integer[n_cnt]; //プリミティブ型だとasListできないため
                final Integer[] randomIdx = new Integer[n_cnt]; //プリミティブ型だとasListできないため
//              for( int i=0; i<n_N; i++) {
                for( int i=0; i<n_cnt; i++) {                           // 8.5.4.2 (2024/01/12) PMD 7.0.0
//                      random_index[i] = i;
                        randomIdx[i] = i;
                }
//              final List<Integer> indexList = Arrays.asList( random_index );
                final List<Integer> indexList = Arrays.asList( randomIdx );

                double localRate = learning_rate;
                for(int epoch=0; epoch<loop; epoch++) {
                        Collections.shuffle( indexList );
        //              random_index = indexList.toArray(new Integer[indexList.size()]);

                        //random_indexの先頭からn_N*minibatch_rate個のものを対象に学習をかける（ミニバッチ）
//                      for(int i=0; i< n_N * minibatch_rate; i++) {
                        for(int i=0; i< n_cnt * minibatch_rate; i++) {          // 8.5.4.2 (2024/01/12) PMD 7.0.0
        //                      final int idx = random_index[i];
                                final int idx = indexList.get(i);
                                train(data[idx], label[idx], localRate);
                        }
                    localRate *= 0.95; //徐々に学習率を下げて振動を抑える。
                }
        }

        /**
         * データを与えて学習をさせます。
         * パラメータの1行を与えています。
         *
         * 0/1のロジスティック回帰の場合は
         * ラベルc(0or1)が各xに対して与えられている時
         * s(x)=σ(Wx+b)=1/(1+ exp(-Wx-b))として、
         * 確率の対数和L(W,b)の符号反転させたものの偏導関数
         * ∂L/∂w=-∑x(c-s(x))
         * ∂L/∂b=-∑=(c-s(x))
         * が最小になるようなW,bの値をパラメータを変えながら求める。
         * というのが実装になる。(=0を求められないため)
         * 多次元の場合はシグモイド関数σ(x)の代わりにソフトマックス関数π(x)を利用して
         * 拡張したものとなる。（以下はソフトマックス関数利用）
         *
         * @og.rev 8.5.4.2 (2024/01/12) PMD 7.0.0 FieldNamingConventions 対応 n_N ⇒ n_cnt
         *
         * @param in_x 1行分のデータ
         * @param in_y xに対するラベル
         * @param lr 学習率
         * @return 差分配列
         */
        private double[] train( final double[] in_x, final int[] in_y, final double lr ) {
                // 8.5.5.1 (2024/02/29) PMD 7.0.0 LocalVariableNamingConventions
//              final double[] p_y_given_x = new double[n_out];
                final double[] givenX   = new double[n_out];
                final double[] dy               = new double[n_out];

                for(int i=0; i<n_out; i++) {
//                      p_y_given_x[i] = 0;
                        givenX[i] = 0;
                        for(int j=0; j<n_in; j++) {
//                              p_y_given_x[i] += vW[i][j] * in_x[j];
                                givenX[i] += vW[i][j] * in_x[j];
                        }
//                      p_y_given_x[i] += vb[i];
                        givenX[i] += vb[i];
                }
//              softmax( p_y_given_x );
                softmax( givenX );

                // 勾配の平均で更新？
                for(int i=0; i<n_out; i++) {
//                      dy[i] = in_y[i] - p_y_given_x[i];
                        dy[i] = in_y[i] - givenX[i];

                        for(int j=0; j<n_in; j++) {
//                              vW[i][j] += lr * dy[i] * in_x[j] / n_N;
                                vW[i][j] += lr * dy[i] * in_x[j] / n_cnt;               // 8.5.4.2 (2024/01/12) PMD 7.0.0
                        }

//                      vb[i] += lr * dy[i] / n_N;
                        vb[i] += lr * dy[i] / n_cnt;                                            // 8.5.4.2 (2024/01/12) PMD 7.0.0
                }

                return dy;
        }

        /**
         * ソフトマックス関数。
         * π(xi) = exp(xi)/Σexp(x)
         * @param in_x 変数X
         */
        private void softmax( final double[] in_x ) {
                // double max = 0.0;
                double sum = 0.0;

                // for(int i=0; i<n_out; i++) {
                //      if(max < x[i]) {
                //              max = x[i];
                //      }
                // }

                for(int i=0; i<n_out; i++) {
                        //x[i] = Math.exp(x[i] - max); // maxとの差分を取ると利点があるのか分からなかった
                        in_x[i] = Math.exp(in_x[i]);
                        sum += in_x[i];
                }

                for(int i=0; i<n_out; i++) {
                        in_x[i] /= sum;
                }
        }

        /**
         * 写像式 Wx+b のW、係数ベクトル。
         * @return 係数ベクトル
         */
        public double[][] getW() {
                return vW;
        }

        /**
         * 写像式 Wx + bのb、バイアス。
         * @return バイアスベクトル
         */
        public double[] getB() {
                return vb;
        }

        /**
         * 出来た予測式に対して、データを入力してyを出力する。
         * (yは各ラベルに対する確率分布となる）
         * @param in_x 予測したいデータ
         * @return 予測結果
         */
        public double[] predict(final double[] in_x) {
                // 8.5.5.1 (2024/02/29) PMD 7.0.0 LocalVariableNamingConventions
//              final double[] out_y = new double[n_out];
                final double[] outY = new double[n_out];

                for(int i=0; i<n_out; i++) {
//                      out_y[i] = 0.;
                        outY[i] = 0.;
                        for(int j=0; j<n_in; j++) {
//                              out_y[i] += vW[i][j] * in_x[j];
                                outY[i] += vW[i][j] * in_x[j];
                        }
//                      out_y[i] += vb[i];
                        outY[i] += vb[i];
                }

//              softmax(out_y);
                softmax(outY);

//              return out_y;
                return outY;
        }

        // ================ ここまでが本体 ================

        /**
         * ここからテスト用mainメソッド 。
         *
         * @og.rev 8.5.4.2 (2024/01/12) PMD 7.0.0 LocalVariableNamingConventions 対応
         *
         * @param args 引数
         */
        public static void main( final String[] args ) {
                // ３つの分類で分ける
//              final double[][] train_X = {
                final double[][] trainX = {
                                {-2.0, 2.0}
                                ,{-2.1, 1.9}
                                ,{-1.8, 2.1}
                                ,{0.0, 0.0}
                                ,{0.2, -0.2}
                                ,{-0.1, 0.1}
                                ,{2.0, -2.0}
                                ,{2.2, -2.1}
                                ,{1.9, -2.0}
                };

//              final int[][] train_Y = {
                final int[][] trainY = {
                                {1, 0, 0}
                                ,{1, 0, 0}
                                ,{1, 0, 0}
                                ,{0, 1, 0}
                                ,{0, 1, 0}
                                ,{0, 1, 0}
                                ,{0, 0, 1}
                                ,{0, 0, 1}
                                ,{0, 0, 1}
                };

                 // test data
//              final double[][] test_X = {
                final double[][] testX = {
                                {-2.5, 2.0}
                                ,{0.1, -0.1}
                                ,{1.5,-2.5}
                };

//              final double[][] test_Y = new double[test_X.length][train_Y[0].length];
                final double[][] testY = new double[testX.length][trainY[0].length];

//              final HybsLogisticRegression hlr = new HybsLogisticRegression( train_X, train_Y, 0.1, 500, 1 );
                final HybsLogisticRegression hlr = new HybsLogisticRegression( trainX, trainY, 0.1, 500, 1 );

                // テスト
                // このデータでは2番目の条件には入りにくい？
//              for(int i=0; i<test_X.length; i++) {
//                       test_Y[i] = hlr.predict(test_X[i]);
//                       System.out.print( Arrays.toString(test_Y[i]) );
//              }
                for(int i=0; i<testX.length; i++) {
                         testY[i] = hlr.predict(testX[i]);
                         System.out.print( Arrays.toString(testY[i]) );
                }
        }
}