diff options
Diffstat (limited to 'vc.cpp')
| -rw-r--r-- | vc.cpp | 443 |
1 files changed, 443 insertions, 0 deletions
diff --git a/vc.cpp b/vc.cpp new file mode 100644 index 0000000..77cf746 --- /dev/null +++ b/vc.cpp @@ -0,0 +1,443 @@ +/* + Genome-wide Efficient Mixed Model Association (GEMMA) + Copyright (C) 2011 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 <bitset> +#include <cstring> + +#include "gsl/gsl_vector.h" +#include "gsl/gsl_matrix.h" +#include "gsl/gsl_linalg.h" +#include "gsl/gsl_blas.h" + +#include "gsl/gsl_cdf.h" +#include "gsl/gsl_multiroots.h" +#include "gsl/gsl_min.h" + +#include "io.h" +#include "lapack.h" +#include "gzstream.h" + +#ifdef FORCE_FLOAT +#include "lmm_float.h" +#include "vc_float.h" +#else +#include "lmm.h" +#include "vc.h" +#endif + + + +using namespace std; + + +//in this file, X, Y are already transformed (i.e. UtX and UtY) + + +void VC::CopyFromParam (PARAM &cPar) +{ + file_out=cPar.file_out; + + // v_sigma2=cPar.v_sigma2; + + time_UtX=0.0; + time_opt=0.0; + + v_traceG=cPar.v_traceG; + + return; +} + + +void VC::CopyToParam (PARAM &cPar) +{ + cPar.time_UtX=time_UtX; + cPar.time_opt=time_opt; + + cPar.v_sigma2=v_sigma2; + cPar.v_se_sigma2=v_se_sigma2; + cPar.v_pve=v_pve; + cPar.v_se_pve=v_se_pve; + cPar.v_traceG=v_traceG; + + cPar.v_beta=v_beta; + cPar.v_se_beta=v_se_beta; + + return; +} + + + +void UpdateParam (const gsl_vector *log_sigma2, VC_PARAM *p) +{ + size_t n1=(p->K)->size1, n_vc=log_sigma2->size-1, n_cvt=(p->W)->size2; + + gsl_matrix *K_temp=gsl_matrix_alloc(n1, n1); + gsl_matrix *HiW=gsl_matrix_alloc(n1, n_cvt); + gsl_matrix *WtHiW=gsl_matrix_alloc(n_cvt, n_cvt); + gsl_matrix *WtHiWi=gsl_matrix_alloc(n_cvt, n_cvt); + gsl_matrix *WtHiWiWtHi=gsl_matrix_alloc(n_cvt, n1); + + double sigma2; + //calculate H=\sum_i^{k+1} \sigma_i^2 K_i + gsl_matrix_set_zero (p->P); + for (size_t i=0; i<n_vc+1; i++) { + if (i==n_vc) { + gsl_matrix_set_identity (K_temp); + } else { + gsl_matrix_const_view K_sub=gsl_matrix_const_submatrix (p->K, 0, n1*i, n1, n1); + gsl_matrix_memcpy (K_temp, &K_sub.matrix); + } + + sigma2=exp(gsl_vector_get (log_sigma2, i) ); + gsl_matrix_scale(K_temp, sigma2); + gsl_matrix_add (p->P, K_temp); + } + + //calculate H^{-1} + int sig; + gsl_permutation * pmt1=gsl_permutation_alloc (n1); + LUDecomp (p->P, pmt1, &sig); + LUInvert (p->P, pmt1, K_temp); + gsl_permutation_free(pmt1); + + gsl_matrix_memcpy (p->P, K_temp); + + //calculate P=H^{-1}-H^{-1}W(W^TH^{-1}W)^{-1}W^TH^{-1} + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, p->P, p->W, 0.0, HiW); + gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, p->W, HiW, 0.0, WtHiW); + + gsl_permutation * pmt2=gsl_permutation_alloc (n_cvt); + LUDecomp (WtHiW, pmt2, &sig); + LUInvert (WtHiW, pmt2, WtHiWi); + gsl_permutation_free(pmt2); + + gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1.0, WtHiWi, HiW, 0.0, WtHiWiWtHi); + gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, -1.0, HiW, WtHiWiWtHi, 1.0, p->P); + + //calculate Py, KPy, PKPy + gsl_blas_dgemv(CblasNoTrans, 1.0, p->P, p->y, 0.0, p->Py); + + for (size_t i=0; i<n_vc+1; i++) { + gsl_vector_view KPy=gsl_matrix_column (p->KPy_mat, i); + gsl_vector_view PKPy=gsl_matrix_column (p->PKPy_mat, i); + + if (i==n_vc) { + gsl_vector_memcpy (&KPy.vector, p->Py); + } else { + gsl_matrix_const_view K_sub=gsl_matrix_const_submatrix (p->K, 0, n1*i, n1, n1); + gsl_blas_dgemv(CblasNoTrans, 1.0, &K_sub.matrix, p->Py, 0.0, &KPy.vector); + } + + gsl_blas_dgemv(CblasNoTrans, 1.0, p->P, &KPy.vector, 0.0, &PKPy.vector); + } + + gsl_matrix_free (K_temp); + gsl_matrix_free (HiW); + gsl_matrix_free (WtHiW); + gsl_matrix_free (WtHiWi); + gsl_matrix_free (WtHiWiWtHi); + + return; +} + + +//below are functions for AI algorithm +int LogRL_dev1 (const gsl_vector *log_sigma2, void *params, gsl_vector *dev1) +{ + VC_PARAM *p=(VC_PARAM *) params; + + size_t n1=(p->K)->size1, n_vc=log_sigma2->size-1; + + double tr, d; + + //update parameters + UpdateParam (log_sigma2, p); + + //calculate dev1=-0.5*trace(PK_i)+0.5*yPKPy + for (size_t i=0; i<n_vc+1; i++) { + if (i==n_vc) { + tr=0; + for (size_t l=0; l<n1; l++) { + tr+=gsl_matrix_get (p->P, l, l); + } + } else { + tr=0; + for (size_t l=0; l<n1; l++) { + gsl_vector_view P_row=gsl_matrix_row (p->P, l); + gsl_vector_const_view K_col=gsl_matrix_const_column (p->K, n1*i+l); + gsl_blas_ddot(&P_row.vector, &K_col.vector, &d); + tr+=d; + } + } + + gsl_vector_view KPy_i=gsl_matrix_column (p->KPy_mat, i); + gsl_blas_ddot(p->Py, &KPy_i.vector, &d); + + d=(-0.5*tr+0.5*d)*exp(gsl_vector_get(log_sigma2, i)); + + gsl_vector_set(dev1, i, d); + } + + return GSL_SUCCESS; +} + + + +int LogRL_dev2 (const gsl_vector *log_sigma2, void *params, gsl_matrix *dev2) +{ + VC_PARAM *p=(VC_PARAM *) params; + + size_t n_vc=log_sigma2->size-1; + + double d, sigma2_i, sigma2_j; + + //update parameters + UpdateParam (log_sigma2, p); + + //calculate dev2=0.5(yPKPKPy) + for (size_t i=0; i<n_vc+1; i++) { + gsl_vector_view KPy_i=gsl_matrix_column (p->KPy_mat, i); + sigma2_i=exp(gsl_vector_get(log_sigma2, i)); + + for (size_t j=i; j<n_vc+1; j++) { + gsl_vector_view PKPy_j=gsl_matrix_column (p->PKPy_mat, j); + + gsl_blas_ddot(&KPy_i.vector, &PKPy_j.vector, &d); + sigma2_j=exp(gsl_vector_get(log_sigma2, j)); + + d*=-0.5*sigma2_i*sigma2_j; + + gsl_matrix_set(dev2, i, j, d); + if (j!=i) {gsl_matrix_set(dev2, j, i, d);} + } + } + + gsl_matrix_memcpy (p->Hessian, dev2); + + return GSL_SUCCESS; +} + + + +int LogRL_dev12 (const gsl_vector *log_sigma2, void *params, gsl_vector *dev1, gsl_matrix *dev2) +{ + VC_PARAM *p=(VC_PARAM *) params; + + size_t n1=(p->K)->size1, n_vc=log_sigma2->size-1; + + double tr, d, sigma2_i, sigma2_j; + + //update parameters + UpdateParam (log_sigma2, p); + + //calculate dev1=-0.5*trace(PK_i)+0.5*yPKPy + //calculate dev2=0.5(yPKPKPy) + for (size_t i=0; i<n_vc+1; i++) { + if (i==n_vc) { + tr=0; + for (size_t l=0; l<n1; l++) { + tr+=gsl_matrix_get (p->P, l, l); + } + } else { + tr=0; + for (size_t l=0; l<n1; l++) { + gsl_vector_view P_row=gsl_matrix_row (p->P, l); + gsl_vector_const_view K_col=gsl_matrix_const_column (p->K, n1*i+l); + gsl_blas_ddot(&P_row.vector, &K_col.vector, &d); + tr+=d; + } + } + + gsl_vector_view KPy_i=gsl_matrix_column (p->KPy_mat, i); + gsl_blas_ddot(p->Py, &KPy_i.vector, &d); + + sigma2_i=exp(gsl_vector_get(log_sigma2, i)); + d=(-0.5*tr+0.5*d)*sigma2_i; + + gsl_vector_set(dev1, i, d); + + for (size_t j=i; j<n_vc+1; j++) { + gsl_vector_view PKPy_j=gsl_matrix_column (p->PKPy_mat, j); + gsl_blas_ddot(&KPy_i.vector, &PKPy_j.vector, &d); + + sigma2_j=exp(gsl_vector_get(log_sigma2, j)); + d*=-0.5*sigma2_i*sigma2_j; + + gsl_matrix_set(dev2, i, j, d); + if (j!=i) {gsl_matrix_set(dev2, j, i, d);} + } + + } + + gsl_matrix_memcpy (p->Hessian, dev2); + + return GSL_SUCCESS; +} + + + + +void VC::CalcVCreml (const gsl_matrix *K, const gsl_matrix *W, const gsl_vector *y) +{ + size_t n1=K->size1, n2=K->size2; + size_t n_vc=n2/n1; + gsl_vector *log_sigma2=gsl_vector_alloc (n_vc+1); + double d, s; + + //set up params + gsl_matrix *P=gsl_matrix_alloc (n1, n1); + gsl_vector *Py=gsl_vector_alloc (n1); + gsl_matrix *KPy_mat=gsl_matrix_alloc (n1, n_vc+1); + gsl_matrix *PKPy_mat=gsl_matrix_alloc (n1, n_vc+1); + gsl_vector *dev1=gsl_vector_alloc (n_vc+1); + gsl_matrix *dev2=gsl_matrix_alloc (n_vc+1, n_vc+1); + gsl_matrix *Hessian=gsl_matrix_alloc (n_vc+1, n_vc+1); + VC_PARAM params={K, W, y, P, Py, KPy_mat, PKPy_mat, Hessian}; + + //initialize sigma2/log_sigma2 + gsl_blas_ddot (y, y, &s); + s/=(double)n1; + for (size_t i=0; i<n_vc+1; i++) { + if (i==n_vc) { + d=s/((double)n_vc+1.0); + } else { + d=s/( ((double)n_vc+1.0)*v_traceG[i]); + } + + gsl_vector_set (log_sigma2, i, d); + } + // gsl_vector_set (log_sigma2, 0, 0.38); + // gsl_vector_set (log_sigma2, 1, -1.08); + + cout<<"iteration "<<0<<endl; + cout<<"sigma2 = "; + for (size_t i=0; i<n_vc+1; i++) { + cout<<exp(gsl_vector_get(log_sigma2, i))<<" "; + } + cout<<endl; + + //set up fdf + gsl_multiroot_function_fdf FDF; + FDF.n=n_vc+1; + FDF.params=¶ms; + FDF.f=&LogRL_dev1; + FDF.df=&LogRL_dev2; + FDF.fdf=&LogRL_dev12; + + //set up solver + int status; + int iter=0, max_iter=100; + + const gsl_multiroot_fdfsolver_type *T_fdf; + gsl_multiroot_fdfsolver *s_fdf; + T_fdf=gsl_multiroot_fdfsolver_hybridsj; + s_fdf=gsl_multiroot_fdfsolver_alloc (T_fdf, n_vc+1); + + gsl_multiroot_fdfsolver_set (s_fdf, &FDF, log_sigma2); + + do { + iter++; + status=gsl_multiroot_fdfsolver_iterate (s_fdf); + + if (status) break; + + cout<<"iteration "<<iter<<endl; + cout<<"sigma2 = "; + for (size_t i=0; i<n_vc+1; i++) { + cout<<exp(gsl_vector_get(s_fdf->x, i))<<" "; + } + cout<<endl; + cout<<"derivatives = "; + for (size_t i=0; i<n_vc+1; i++) { + cout<<gsl_vector_get(s_fdf->f, i)<<" "; + } + cout<<endl; + + status=gsl_multiroot_test_residual (s_fdf->f, 1e-3); + } + while (status==GSL_CONTINUE && iter<max_iter); + + //obtain Hessian inverse + int sig=LogRL_dev12 (s_fdf->f, ¶ms, dev1, dev2); + + gsl_permutation * pmt=gsl_permutation_alloc (n_vc+1); + LUDecomp (dev2, pmt, &sig); + LUInvert (dev2, pmt, Hessian); + gsl_permutation_free(pmt); + + //save data + v_sigma2.clear(); + for (size_t i=0; i<n_vc+1; i++) { + d=exp(gsl_vector_get(s_fdf->x, i)); + v_sigma2.push_back(d); + } + + v_se_sigma2.clear(); + for (size_t i=0; i<n_vc+1; i++) { + d=-1.0*v_sigma2[i]*v_sigma2[i]*gsl_matrix_get(Hessian, i, i); + v_se_sigma2.push_back(sqrt(d)); + } + + s=0; + for (size_t i=0; i<n_vc; i++) { + s+=v_traceG[i]*v_sigma2[i]; + } + s+=v_sigma2[n_vc]; + + v_pve.clear(); + for (size_t i=0; i<n_vc; i++) { + d=v_traceG[i]*v_sigma2[i]/s; + v_pve.push_back(d); + } + + v_se_pve.clear(); + for (size_t i=0; i<n_vc; i++) { + d=v_traceG[i]*(s-v_sigma2[i]*v_traceG[i])/(s*s)*v_se_sigma2[i]*v_se_sigma2[i]; + v_se_pve.push_back(sqrt(d) ); + } + + gsl_multiroot_fdfsolver_free(s_fdf); + + gsl_vector_free(log_sigma2); + gsl_matrix_free(P); + gsl_vector_free(Py); + gsl_matrix_free(KPy_mat); + gsl_matrix_free(PKPy_mat); + gsl_vector_free(dev1); + gsl_matrix_free(dev2); + gsl_matrix_free(Hessian); + + return; +} + + + + + + |