forked from rgcgithub/regenie
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathSKAT.hpp
118 lines (100 loc) · 11.7 KB
/
SKAT.hpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
/*
This file is part of the regenie software package.
Copyright (c) 2020-2024 Joelle Mbatchou, Andrey Ziyatdinov & Jonathan Marchini
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
*/
#ifndef SKAT_H
#define SKAT_H
// SKAT
void update_vc_gmat(SpMat&,Eigen::ArrayXd&,Eigen::ArrayXd&,ArrayXb&,const int&,const int&,struct param const&,const Eigen::Ref<const ArrayXb>&,Eigen::Ref<Eigen::MatrixXd>,std::vector<variant_block>&,const Eigen::Ref<MatrixXb>);
void update_vc_gmat(SpMat&,Eigen::ArrayXd&,Eigen::ArrayXd&,SpMat const&,const Eigen::Ref<const ArrayXb>&,const Eigen::Ref<const ArrayXb>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const ArrayXb>&,struct param const&);
bool get_custom_weights(std::string const&,Eigen::Ref<Eigen::ArrayXd>,std::vector<snp>&,std::vector<uint64> const&);
bool get_custom_weights(std::string const&,Eigen::Ref<Eigen::ArrayXd>,std::vector<snp>&,const Eigen::Ref<const Eigen::ArrayXi>&,std::vector<uint64> const&);
void compute_vc_masks(SpMat&,Eigen::Ref<Eigen::ArrayXd>,Eigen::Ref<Eigen::ArrayXd>,SpMat&,Eigen::Ref<MatrixXb>,const Eigen::Ref<const Eigen::MatrixXd>&, struct ests const&,struct f_ests const&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const MatrixXb>&,MatrixXb&,std::vector<variant_block>&,const Eigen::Ref<const ArrayXb>&,struct param const&,struct remeta_sumstat_writer&);
void prep_ultra_rare_mask(SpMat&,Eigen::Ref<Eigen::ArrayXd>,Eigen::Ref<Eigen::ArrayXd>,SpMat&,Eigen::Ref<MatrixXb>,MatrixXb&,const Eigen::Ref<const ArrayXb>&,struct param const&);
void compute_vc_masks_qt(SpMat&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const MatrixXb>&,std::vector<variant_block>&,struct param const&,struct remeta_sumstat_writer&);
void compute_vc_masks_qt_fixed_rho(SpMat&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const MatrixXb>&,std::vector<variant_block>&,double const&,double const&,double const&,uint const&,bool const&,struct param const&,struct remeta_sumstat_writer&);
void compute_vc_masks_qt(SpMat&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const MatrixXb>&,std::vector<variant_block>&,const Eigen::Ref<const Eigen::ArrayXd>&,double const&,double const&,uint const&,bool const&,struct param const&,struct remeta_sumstat_writer&);
void compute_vc_mats_qt(Eigen::Ref<Eigen::MatrixXd>,Eigen::Ref<Eigen::MatrixXd>,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<SpMat>);
void compute_skat_q(Eigen::MatrixXd&,Eigen::MatrixXd&,Eigen::Ref<Eigen::MatrixXd>,const Eigen::Ref<const Eigen::MatrixXd>&,Eigen::MatrixXd&,const Eigen::Ref<const ArrayXb>&,const Eigen::Ref<const MatrixXb>&,bool const&,bool const&);
void get_acatv_pv(int const&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,double&,double&,double const&,bool const&);
void get_single_pvs(Eigen::Ref<Eigen::MatrixXd>,const Eigen::Ref<const Eigen::MatrixXd>&);
void compute_vc_masks_bt(SpMat&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,struct ests const&,struct f_ests const&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const MatrixXb>&,const Eigen::Ref<const MatrixXb>&,std::vector<variant_block>&,struct param const&,struct remeta_sumstat_writer&);
void compute_vc_masks_bt_fixed_rho(SpMat&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,struct ests const&,struct f_ests const&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const MatrixXb>&,const Eigen::Ref<const MatrixXb>&,std::vector<variant_block>&,double const&,double const&,double const&,bool const&,uint const&,bool const&,struct param const&,struct remeta_sumstat_writer&);
void compute_vc_masks_bt(SpMat&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,struct ests const&,struct f_ests const&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const MatrixXb>&,const Eigen::Ref<const MatrixXb>&,std::vector<variant_block>&,const Eigen::Ref<const Eigen::ArrayXd>&,double const&,double const&,bool const&,uint const&,bool const&,struct param const&,struct remeta_sumstat_writer&);
void get_single_pvs_bt(Eigen::Ref<Eigen::ArrayXd>,const Eigen::Ref<const Eigen::ArrayXd>&);
Eigen::MatrixXd get_RsKRs(const Eigen::Ref<const Eigen::MatrixXd>&,const double&,const double&);
Eigen::MatrixXd get_RsKRs(const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const double&,const double&,const double&);
void get_lambdas(Eigen::VectorXd&,const Eigen::Ref<const Eigen::MatrixXd>&,const double&);
void compute_fixed_skato_p(double&,double&,double const&,double const&,double const&,Eigen::VectorXd&,const double&,bool const&);
void compute_fixed_skato_p(double&,double&,double&,double const&,Eigen::VectorXd&,const double&);
void compute_skat_pv(double&,double&,double const&,Eigen::VectorXd&,const double&);
double get_chisq_mix_pv(double const&,const Eigen::Ref<const Eigen::VectorXd>&);
double get_chisq_mix_logp(double const&,const Eigen::Ref<const Eigen::VectorXd>&,double&);
double get_davies_pv(double const&,Eigen::Ref<Eigen::VectorXd>,bool const&);
double get_kuonen_pv(const double&,const Eigen::Ref<const Eigen::VectorXd>&);
double get_liu_pv(const double&,const Eigen::Ref<const Eigen::VectorXd>&,const bool& lax = false);
double get_liu_pv(const double&,const Eigen::Ref<const Eigen::VectorXd>&,double&);
double get_tmin_lambda(const double&,const Eigen::Ref<const Eigen::ArrayXd>&);
double get_tmax_lambda(const Eigen::Ref<const Eigen::ArrayXd>&);
void solve_kp(bool&,double&,const double&,const double&,const double&,const Eigen::Ref<const Eigen::ArrayXd>&);
bool valid_bounds(double&,double const&,double const&,const double&,const Eigen::Ref<const Eigen::ArrayXd>&);
double K_lambda(const double&,const Eigen::Ref<const Eigen::ArrayXd>&);
double Kp_lambda(const double&,const Eigen::Ref<const Eigen::ArrayXd>&);
double Kpp_lambda(const double&,const Eigen::Ref<const Eigen::ArrayXd>&);
double get_spa_pv(const double&,const double&,const Eigen::Ref<const Eigen::ArrayXd>&);
void compute_vc_mats_bt(Eigen::Ref<Eigen::ArrayXd>,Eigen::Ref<Eigen::MatrixXd>,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,Eigen::Ref<SpMat>,SpMat&,Eigen::MatrixXd&);
void compute_skat_q(Eigen::VectorXd&,Eigen::VectorXd&,const Eigen::Ref<const Eigen::ArrayXd>&,Eigen::Ref<Eigen::MatrixXd>,const Eigen::Ref<const MatrixXb>&,bool const&);
void correct_vcov(const int&,const Eigen::Ref<const Eigen::ArrayXi>&,const Eigen::Ref<const Eigen::ArrayXd>&,Eigen::Ref<ArrayXb>,Eigen::Ref<Eigen::ArrayXd>,const Eigen::Ref<const Eigen::ArrayXd>&,Eigen::Ref<Eigen::MatrixXd>,SpMat const&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,SpMat const&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const ArrayXb>&,struct f_ests const&,struct param const&);
void apply_correction_cc(const int&,const Eigen::Ref<const Eigen::ArrayXi>&,const Eigen::Ref<const Eigen::ArrayXd>&,Eigen::Ref<Eigen::ArrayXd>,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,SpMat const&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,SpMat const&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const ArrayXb>&,struct f_ests const&,struct param const&,bool const&);
void apply_firth_snp(bool&,double&,double const&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::ArrayXi>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const ArrayXb>&,struct param const&);
bool correct_vcov_burden(const int&,double&,double const&,double const&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,SpMat const&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const ArrayXb>&,const Eigen::Ref<const Eigen::MatrixXd>&,struct param const&);
bool get_ztz_evals(const Eigen::Ref<const Eigen::MatrixXd>&,Eigen::MatrixXd&,double&,double&,double&,double const&,bool const&);
void get_skato_mom(double&,double&,double&,double&,Eigen::ArrayXd&,const Eigen::Ref<const Eigen::VectorXd>&,double const&,double const&,double const&,const Eigen::Ref<const Eigen::ArrayXd>&,bool const&);
void get_cvals(int const&,Eigen::Ref<Eigen::MatrixXd>,const Eigen::Ref<const Eigen::VectorXd>&);
void get_cvals(Eigen::Ref<Eigen::ArrayXd>,const Eigen::Ref<const Eigen::VectorXd>&);
void get_Qmin(int const&,double&,Eigen::Ref<Eigen::ArrayXd>,const Eigen::Ref<const Eigen::MatrixXd>&);
void get_skato_pv(double &,double&,double const&,int const&,double const&,bool const&);
void print_vc_sumstats(int const&,std::string const&,std::string const&,variant_block*,std::vector<snp> const&,struct in_files const&,struct param const*);
// for numerical integration with skat-o
#ifdef __cplusplus
extern "C"
{
#endif
extern void dqags_(double f(double*),double*,double*,double*,double*,double*,double*,int*,int*,int*,int*,int*,int*,double*);
double SKATO_integral_fn(double*);
double SKATO_integral_fn_liu(double*);
#ifdef __cplusplus
}
#endif
// declare global variables
extern Eigen::ArrayXd flipped_skato_rho;
extern Eigen::ArrayXd skato_Qmin_rho;
extern Eigen::ArrayXd skato_tau;
extern Eigen::VectorXd skato_lambdas;
extern double skato_muQ;
extern double skato_fdavies;
extern double skato_sdQ;
extern double skato_dfQ;
extern double skato_upper;
extern int skato_state; // positive if integration failed
void integrate(double f(double*),double&,int const&,bool const&);
// for lovo with bts
extern Eigen::MatrixXd vc_Rvec_start;
void check_cc_correction(SpMat&,const Eigen::Ref<const Eigen::ArrayXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,struct ests const&,struct f_ests const&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const Eigen::MatrixXd>&,const Eigen::Ref<const MatrixXb>&,struct param const&);
void check_sizes(SpMat const&, SpMat const&, const Eigen::Ref<const MatrixXb>&);
#endif