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_in;		// データ次元
	private final int n_out;	// ラベル種別数

	// 写像変数ベクトル f(x) = Wx + b
	private double[][] vW;
	private double[] vb;

	/**
	 * コンストラクタ。
	 * 
	 * 学習もしてしまう。
	 * 
	 * xはデータセット各行がn次元の説明変数となっている。
	 * trainはそれに対する{0,1,0},{1,0,0}のようなラベルを示すベクトルとなる。
	 * 学習率は通常、0.1程度を設定する。
	 * このロジックではループ毎に0.95をかけて徐々に学習率が下がるようにしている。
	 * 全データを利用すると時間がかかる場合があるので、確率的勾配降下法を利用しているが、
	 * 選択個数はデータに対する割合を与える。
	 * データ個数が少ない場合は1をセットすればよい。
	 * 
	 * @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_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できないため
		for( int i=0; i<n_N; i++) {
			random_index[i] = i; 
		}
		final List<Integer> indexList = Arrays.asList( random_index );

		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++) {
	//			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)を利用して
	 * 拡張したものとなる。（以下はソフトマックス関数利用）
	 * 
	 * @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 ) {
		final double[] p_y_given_x = 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;
			for(int j=0; j<n_in; j++) {
				p_y_given_x[i] += vW[i][j] * in_x[j];
			}
			p_y_given_x[i] += vb[i];
		}
		softmax( p_y_given_x );

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

			for(int j=0; j<n_in; j++) {
				vW[i][j] += lr * dy[i] * in_x[j] / n_N;
			}

			vb[i] += lr * dy[i] / n_N;
		}

		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) {
		final double[] out_y = new double[n_out];

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

		softmax(out_y);

		return out_y;
	}

	//************** ここまでが本体 **************
	/**
	 * ここからテスト用mainメソッド 。
	 *
	 * @param args 引数
	 */
	public static void main( final String[] args ) {
		// ３つの分類で分ける
		final double[][] train_X = {
				{-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 = {
				{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 = {
				{-2.5, 2.0}
				,{0.1, -0.1}
				,{1.5,-2.5}
		};

		double[][] test_Y = new double[test_X.length][train_Y[0].length];

		final HybsLogisticRegression hlr = new HybsLogisticRegression( train_X, train_Y, 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]) );
		}
	}
}

