/*
 Genome-wide Efficient Mixed Model Association (GEMMA)
 Copyright (C) 2011-2017, Xiang Zhou
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#include <iostream>
#include <fstream>
#include <sstream>

#include <iomanip>
#include <cmath>
#include <iostream>
#include <stdio.h>
#include <stdlib.h> 
#include <ctime>
#include <cstring>
#include <algorithm>

#include "gsl/gsl_vector.h"
#include "gsl/gsl_matrix.h"
#include "gsl/gsl_linalg.h"
#include "gsl/gsl_blas.h"
#include "gsl/gsl_eigen.h"
#include "gsl/gsl_randist.h"
#include "gsl/gsl_cdf.h"
#include "gsl/gsl_roots.h"
#include "Eigen/Dense"

#include "lapack.h"
#include "param.h"
#include "ldr.h"
#include "lm.h"
#include "mathfunc.h"

using namespace std;
using namespace Eigen;

void LDR::CopyFromParam (PARAM &cPar) {
	a_mode=cPar.a_mode;
	d_pace=cPar.d_pace;
	
	file_bfile=cPar.file_bfile;
	file_geno=cPar.file_geno;
	file_out=cPar.file_out;
	path_out=cPar.path_out;

	ni_total=cPar.ni_total;
	ns_total=cPar.ns_total;
	ni_test=cPar.ni_test;
	ns_test=cPar.ns_test;
	n_cvt=cPar.n_cvt;
	
	indicator_idv=cPar.indicator_idv;
	indicator_snp=cPar.indicator_snp;
	snpInfo=cPar.snpInfo;
	
	return;
}

void LDR::CopyToParam (PARAM &cPar) {
	return;
}

//X is a p by n matrix.
void LDR::VB (const vector<vector<unsigned char> > &Xt,
	      const gsl_matrix *W_gsl, const gsl_vector *y_gsl) {
  
  // Save gsl_vector and gsl_matrix into Eigen library formats.
  MatrixXd W(W_gsl->size1, W_gsl->size2);
  VectorXd y(y_gsl->size);
  VectorXd x_col(y_gsl->size);

  double d;
  for (size_t i=0; i<W_gsl->size1; i++) {
    d=gsl_vector_get(y_gsl, i);
    y(i)=d;
    for (size_t j=0; j<W_gsl->size2; j++) {
      W(i,j)=gsl_matrix_get(W_gsl, i, j);
    }
  }

  // Initial VB values by lm.
  cout<<indicator_snp[0]<<" "<<indicator_snp[1]<<" "<<indicator_snp[2]<<endl;
  uchar_matrix_get_row (Xt, 0, x_col);

  for (size_t j=0; j<10; j++) {
    cout<<x_col(j)<<endl;
  }

  // Run VB iterations.
  // TO DO.

  // Save results.
  // TO DO.
  
  return;
}