#include <glmm.h>
Inheritance diagram for matvec::GLMM:
Public Types | |
| enum | ModelType { bad_model, fixed_model, mixed_model, reg_model } |
Public Methods | |
| virtual double | contrast (const Vector< double > &Kp, const double m=0.0) |
| virtual double | contrast (const doubleMatrix &Kp, const Vector< double > &M) |
| virtual double | contrast (const doubleMatrix &Kp) |
| virtual double | contrast (const std::string &termname, const doubleMatrix &Kp, const Vector< double > &M) |
| int | equation (const std::string &modelspecs) |
| virtual double | contrast (const double **kpme, const unsigned nr, const unsigned nc, const double *M=0, const int prt_flag=1, double **result=0) |
| GLMM () | |
| ~GLMM () | |
| void | SSQCP (void) |
| void ** | Param (void **par=0) |
| int | Large_Sample (int Sample=-1) |
| int | Pev_scale (int Sample=-1) |
| int | LR_Test (int LR_set=-1) |
| int | LM_Test (int LM_set=-1) |
| int | Data_static (int DS_set=-1) |
| int | Like_LR (int L_set=-1) |
| int | Westell (int W_set=-1) |
| int | AI_Sample (int S_set=-1) |
| void | Initial_eta (Vector< double > initial) |
| void | residual (int get_pev=0) |
| void | new_SinMat (void) |
| double | log_likelihood (int solve=1) |
| double | restricted_log_likelihood (int solve=1) |
| void | release_SinMat (void) |
| void | Build_SinMat (void) |
| int | print_level (int level=-99) |
| int | return_stat (int level=-99) |
| int | print_summary (int level=-99) |
| Observe * | get_obs (void) |
| doubleMatrix * | get_SinMat (void) |
| Vector< double > & | resid (int i) |
| Vector< double > & | Linear_Predictor (int i) |
| doubleMatrix & | pev (int i) |
| double | log_like (void) |
| doubleMatrix | INFO (void) |
| Vector< double > | SCORE (void) |
| Observe & | new_obs () |
| void | release_obs () |
| SparseMatrix & | ainv (void) |
| unsigned | Numrec (void) |
| doubleMatrix | AI_REML (int numiter=10, double tol=1.e-4, double info_scale=0) |
| virtual SparseMatrix * | setup_mme (Vector< double > *rhs=0) |
| void | link (void(*link_fn)(Observe &), int rvc=0) |
| unsigned | TotalNvc (void) const |
| int | Var2Vec (double *x) |
| void | Vec2Var (const double *x) |
| void | covariance_old (const std::string &termname, const std::string &termname2, const doubleMatrix &v) |
| Vector< double > * | glim (int iterations=10, double **ke=0, int nr=0, int nc=0, double *mean=0) |
| void | add_G_Sand (void) |
| void | add_IgSand (int t) |
| void | add_AgSand (int t) |
| void | add_Ag (int t, int ct) |
| void | add_Ag (int t, SparseMatrix &h) |
| void | add_Ag_old (int t) |
| void | add_Ig (int t, SparseMatrix &h) |
| void | add_G_1 (SparseMatrix &h) |
| void | add_G_1_old () |
| void | build_CorrVar (void) |
| void | build_hInv (void) |
| void | info (std::ostream &stream) |
| void | save (const std::string &fname, const int io_mode=std::ios::out) |
| void | normalLog (void) |
| doubleMatrix | getVarEstimates (const std::string &termname) |
| int | display (void) const |
| virtual int | prepare_data (const std::string &solver="ysmp") |
| unsigned | nonzero (void) const |
| void | nonzero (const unsigned nz) |
| unsigned | mmesize (void) const |
| unsigned | ntrait (void) const |
| unsigned | nterm (void) const |
| unsigned | totalnvc (void) const |
| void | trait_effect_level (const unsigned mmeaddr, std::string &retval, int info_flag=1) |
| void | initialize (void) |
| void | copyfrom (const Model &M) |
| void | release (void) |
| void | prior_dist (const std::string &termname, Pedigree &P, GeneticDist *D) |
| void | prior_dist (const std::string &termname, GeneticDist *D) |
| void | DGSamplerSetup (unsigned numLoci, Data *D) |
| void | DGSampler (unsigned numOfSamples, unsigned numOfSL, unsigned numOfHaplo, unsigned numOfSLCas) |
| void | RSamplerPrior_dist (const std::string &termname, Pedigree &P, Data *D, GeneticDist *G) |
| method to specify the prior distributions for each terms in the model for which the R-sampler is invoked; it also sets up the structures needed by the peeling process. | |
| void | RSampler (RSamplerParms rsParms) |
| void | RSampler (std::string fileName) |
| void | RSampler (unsigned pedBlockSize, unsigned numLoci, unsigned numOfSamples, const std::string &samplerType, const std::string &howToSample, const std::string &samplerUsed, const std::string &whatToCompute) |
| void | RSamplerInitialDG (const std::string &samplerType="genotypic", const std::string &whatToCompute="probDescHaplo") |
| method to generate the initial descent graph to be used by pop_graph | |
| void | RSamplerGibbs (unsigned pedBlockSize, unsigned numLoci, unsigned numOfSamples, const std::string &samplerType="genotypic", const std::string &whatToCompute="genotypeProb") |
| method to sample ordered genotypes using the Gibbs sampler to sample across loci ( the first sample is obtained by sequential imputation) | |
| void | RSamplerMH (unsigned pedBlockSize, unsigned numLoci, unsigned numOfSamples, const std::string &samplerType="genotypic", const std::string &whatToCompute="genotypeProb") |
| method to sample ordered genotypes using the Metropolis Hastings algorithm to accept samples proposed by sequential imputation. | |
| void | RSamplerMH (unsigned numLoci, unsigned numOfSamples, const std::string &samplerType="allelic", const std::string &whatToCompute="genotypeProb") |
| void | RSamplerGibbsMH (unsigned pedBlockSize, unsigned numLoci, unsigned numOfSamples, unsigned stepGibbsMH, const std::string &samplerType="genotypic", const std::string &whatToCompute="genotypeProb") |
| Matrix< double > | getIBDMatrix (void) |
| void | variance (const std::string &termname, const doubleMatrix &v) |
| void | variance (const std::string &termname, Pedigree &P, const doubleMatrix &v) |
| void | variance (const std::string &termname, const double v00,...) |
| void | variance (const std::string &termname, Pedigree &P, const double v00,...) |
| void | VarLink (const std::string &termname, const std::string &termname2) |
| void | covariance (const std::string &termname, const std::string &termname2, const doubleMatrix &v) |
| void | covariance (const std::string &termname, const std::string &termname2, const double v00,...) |
| void | weight (const std::string &wtname, const int r=0) |
| void | covariate (const std::string &factorname) |
| void | fitdata (Data &D) |
| void | info (std::ostream &stream) const |
| void | var2vec (Vector< double > &x) |
| void | vec2var (const Vector< double > &x) |
| void | release_mme (void) |
| void | vce_gibbs_sampler (const unsigned nvc, Vector< double > &vc, const Vector< double > &ratio, const Vector< unsigned > &nlevel, const double yy, double &ee, Vector< double > &uAu, Vector< double > &bu, Vector< double > &wspace, const Vector< double > &zz) |
| void | update_mme (const Vector< double > &ratio, const Vector< double > &oldratio, const Vector< double > &zz) |
| void | compute_xbzu (Vector< double > &bu) |
| Vector< double > | get_solutions (const std::string &termname) |
| void | sampleW (Vector< double > &sol, Vector< double > &newrhs) |
| void | vce_emreml_single_trait (const int miter=1000, const int intive=1, const double stoptol=1.0e-5) |
| void | vce_emreml_multi_trait (const int miter=1000, const int intive=1, const double stoptol=1.0e-5) |
| void | uAu_trCA (Vector< double > &tsol, Vector< double > &uAu, Vector< double > &trca, Vector< double > &ivect, Vector< double > &invect, const Vector< double > &ratio, const Vector< double > &zz) |
| Vector< double > * | blup (const std::string &method="ysmp", double stopval=0.001, double relax=1.0, unsigned mxiter=100000) |
| void | blup_pccg (double stopval=1.0e-8, unsigned int mxiter=100000) |
| Vector< double > * | blue (const std::string &method="ysmp", double stopval=0.001, double relax=1.0, unsigned mxiter=100000) |
| Vector< double > * | get_blupsol (const std::string &method="ysmp", double stopval=0.001, double relax=1.0, unsigned mxiter=100000) |
| Vector< double > * | vce_emreml (const int miter=1000, const int intive=1, const double stoptol=1.0e-4) |
| Vector< double > * | vce_gibbs (const unsigned warmup=2000, const unsigned gibbslen=5000) |
| Vector< double > * | vce_dfreml (const std::string &method="powell", int maxiter=500000, double funtol=1.0e-5) |
| Vector< double > * | genotypic_value_peeling (void) |
| Vector< double > * | genotypic_value_gibbs (const unsigned warmup=30, const unsigned gibbslen=1000) |
| Vector< double > * | genotypic_value_cat (const unsigned warmup, const unsigned gibbslen) |
| double | genotype_dist_gibbs (const unsigned warmup=30, const unsigned gibbslen=1000) |
| void | genotype_dist_peeling (const int prtlevel=1, const int estifreq=1) |
| double | restricted_log_likelihood (void) |
| double | log_likelihood_peeling (const unsigned maxit=3) |
| double | log_likelihood_gibbs (const unsigned wup=30, const unsigned glen=1000) |
| double | minfun (const Vector< double > &x, const int n) |
| double | minfun_peeling (const Vector< double > &x, const int n) |
| double | minfun_vce (const Vector< double > &x, const int n) |
| double | estimate_peeling (void) |
| double | ml_estimate_peeling (void) |
| double | estimate (const Vector< double > &Kp) |
| Vector< double > | estimate (const doubleMatrix &Kp) |
| Vector< double > | estimate (const std::string &termname, const doubleMatrix &Kp) |
| void | estimate (const double **kpme, const unsigned nr, const unsigned nc, double *kpb) |
| doubleMatrix | CovMat (doubleMatrix &Kp) |
| doubleMatrix | CovMat (const std::string &termname, doubleMatrix &Kp) |
| virtual double | contrast (const double **kpme, const int nr, const int nc, const double *M=0, const int prt_flag=1, double **result=0) |
| void | anova (double &ssm, double &sse, unsigned &rank_x, int &dfe) |
| void | lsmeans (const std::string &termnames, const std::string &filename="", const int io_mode=std::ios::out, const int savekp_flag=0) |
| SparseMatrix * | mmec (void) |
| Vector< double > * | rhs (void) |
| double | multipoint (int maxit) |
| void | multipoint_setup (void) |
| void | multipoint_setup (Fpmm &F, char map_fun='H', int map_par=2) |
| void | multipoint_tables (void) |
| void | multipoint_QTL_table (void) |
| void | multipoint_Rec_table (void) |
| void | multipoint_Rec_recalc (double distance) |
| double | MapF (double dist, int qr=0) |
| double | minfun_multipoint (const Vector< double > &x, const int n) |
| double | ProbGamete (Vector< int > value, int nloci) |
| void | multipoint_profile_distance (int method, int maxit, double Min, double Max, double step_size, int Print, const std::string &fname) |
| double | multipoint_ml_estimate (const int method=0, int niter=0, double tol=1.0e-16, double epsilon=1.0e-8, int printlevel=0) |
| void | multipoint_estimate_Ve (double Initial, double Max, double Min) |
| void | multipoint_estimate_Distance (double Initial, double Max, double Min) |
| void | multipoint_estimate_Qmeans (Vector< double > Initial, Vector< double > Max, Vector< double > Min, char mtype='G') |
| void | multipoint_estimate_Qfreq (double Initial, double Max, double Min) |
| void | multipoint_genot_probs (void) |
| void | multipoint_compute_Rec_table (void) |
| double | multipoint_mixture_approx (int maxit=0, double P_var=1.0) |
| void | multipoint_display_parms (void) |
| void | multipoint_estimate_PolyV (double Initial, double Max, double Min) |
| unsigned | get_popsize (void) const |
| void | descent_graph_peeling_initialisation (void) |
| Vector< int > | get_gamete_vector (unsigned index) |
| void | calc_gamete_prob (void) |
| void | map (void) |
| double | graph_log_likelihood_peeling (const unsigned maxiter) |
| std::string | label (const std::string &termname, const unsigned i) |
| void | min_prtlevel (const int pl) |
| double | praxis (Vector< double > &x, const int n, int &maxfun, const double tol=1.0e-16, const double epsilon=1.0e-8, const int pl=2) |
Public Attributes | |
| int | nrandom |
| int | num_glmm |
| int * | nt_vec2 |
| unsigned | num_ped |
| unsigned | act_numtrait |
| unsigned | data_static |
| RSamplerParms | myRSamplerParms |
| ModelType | type |
| Population * | pop |
| Data * | data |
| Matrix< double > * | rve |
| doubleMatrix | residual_var |
| Vector< double > | lng0vec |
| TermList | term |
| int | Q_pos |
| int | maxit |
| int | map_parm |
| char | map_function |
| double | new_distance |
| doubleMatrix | RecoTable |
| int | N_Parm_List |
| Vector< Parm_Elem > | Parm_List |
| Vector< double > | map_distance |
| doubleMatrix | gamete_prob_table |
| doubleMatrix | Cr_table |
Protected Methods | |
| virtual void | setup_ww (Vector< double > *rhs) |
| void | add_Ig (int t, Vector< double > *x_p=0) |
| void | add_Ig_diag (int t, double **M) |
| void | add_Ag (int t, Vector< double > *x_p=0) |
| void | add_Ag_diag (int t, double **M) |
| virtual void | add_G_1 (void) |
| virtual void | add_G_1_diag (double **M) |
| void | add_G_1_single_trait (const Vector< double > &ratio) |
| double | setup_ww_single_trait (const Vector< double > &ratio, const unsigned pt, const unsigned ibeg, Vector< double > *zz) |
| void | compute_rhs_diag_mme (double **M) |
| void | mme_times (Vector< double > &x) |
| int | build_model_term (const std::string &modelspecs) |
| void | inverse_residual_var (void) |
| void | assign_id_xact (const Vector< bool > &) |
| void | hashxact (const Vector< bool > &) |
| void | re_hash_data (Vector< bool > &) |
| virtual void | save_pos_val (Vector< bool > &, const std::string &solver="ysmp") |
| void | input_pos_val (std::istream &modelfile) |
| void | input_pos_val_iod (std::vector< pos_val_node >::iterator it) |
| void | get_lms_kp (const unsigned termth, const unsigned lvlth, const unsigned yth, Vector< double > &kp, const int *lsmtablevec) |
| void | out_lsmeans_to_stream (std::ostream &ofs, const unsigned nlsm, const Matrix< int > &lsmtable, const int savekp=0) |
Protected Attributes | |
| double | Marq_lambda |
| void ** | param |
| double | numy |
| double | ldet |
| int | aireml_called |
| int | Need_Ainv |
| int | Need_Residual |
| int | res_nvc |
| int | Return_Stat |
| int | do_west |
| int | num_west |
| int | AI_sample |
| int | Print_Level |
| int | PEV_Scale |
| int | Use_Like |
| int | Print_Summary |
| int | Update_estimate |
| int | Sand_Calc |
| int | LM |
| int | Need_PEV |
| Observe * | observation |
| doubleMatrix * | SinMat |
| doubleMatrix * | varinv |
| doubleMatrix * | Var |
| doubleMatrix ** | corrvar |
| doubleMatrix | CorrVar |
| int | ncorr |
| int * | corrmap |
| doubleMatrix | Info |
| Vector< double > | Score |
| Vector< double > | varcomp |
| Vector< double > | varest |
| Vector< int > | kvec |
| double | like_val |
| double | like_adj |
| doubleMatrix | Fmat |
| doubleMatrix | FRHS |
| doubleMatrix | FSol |
| SparseMatrix | hainv |
| SparseMatrix | hSand |
| SparseMatrix | hInv |
| SparseMatrix | itilde |
| SparseMatrix | ihat |
| SparseMatrix | jtilde |
| SparseMatrix | jhat |
| SparseMatrix | S |
| Vector< double > | sol_null |
| Vector< double > | qhat |
| void(* | link_function )(Observe &) |
| Vector< double > | Initial |
| int | Asym |
| int | LR |
| std::vector< pos_val_node > | pos_val_vector |
| std::string | modelfname |
| std::string | weightname |
| std::string | modelstring |
| std::string * | traitname |
| std::string | modelstringstr |
| char * | modelstr |
| int | modelpcount |
| unsigned | numtrait |
| unsigned | numterm |
| unsigned | numfact |
| unsigned | maxorder |
| unsigned | numrec |
| unsigned | numcol |
| unsigned | popsize |
| unsigned | numobs |
| unsigned | hmmesize |
| unsigned | non_zero |
| unsigned | max_nz |
| unsigned | npattern |
| unsigned | ntermGdist |
| unsigned * | rawpos |
| unsigned | numgroup |
| int | weightinverse |
| int | data_prepared |
| int | pop_created |
| int | categorical |
| std::vector< std::string > | pattern |
| std::set< std::string > | pattern_tree |
| UnknownDist * | unknown_prior |
| FieldStruct * | factor_struct |
| FieldStruct ** | trait_struct |
| HashTable * | xact_htable |
| HashTable ** | idlist |
| Vector< double > | lnr0vec |
| Vector< double > | blupsol |
| Vector< double > | rellrhs |
| Vector< double > | diag_mme |
| Vector< double > | mme_times_res |
| double | yry |
| unsigned | nfunk_in_dfreml |
| unsigned | dfreml_called |
| double | reml_value |
| std::string | dfreml_method |
| SparseMatrix | hmmec |
| Vector< int > | base_effect |
| Vector< int > | corr_link |
| Vector< int > | pos_vec |
| Vector< Vector< int > > | trait_map |
| Vector< int > | nt_vec |
| Vector< int > | var_link |
| unsigned * | pos_term |
| double * | xval_term |
| double * | trait_vec |
| unsigned * | rec_indid |
| int | minfun_indx |
|
|
Definition at line 173 of file model.h.
00173 {bad_model, fixed_model, mixed_model, reg_model};
|
|
|
Definition at line 117 of file glmm.h. References AI_sample, aireml_called, Asym, corrmap, do_west, like_val, link_function, LM, LR, Marq_lambda, ncorr, Need_Ainv, Need_PEV, Need_Residual, nrandom, nt_vec2, num_glmm, num_west, numy, param, PEV_Scale, Print_Level, Print_Summary, res_nvc, Return_Stat, Sand_Calc, Update_estimate, and Use_Like.
00117 {observation = 0;SinMat=0;Need_Ainv=1;varinv=0;like_val=0;Need_PEV=0;
00118 Var=0;link_function=0;Need_Residual=0;res_nvc=0;Print_Level=0;Print_Summary=0;Asym=0;LR=0;PEV_Scale=1;Use_Like=0;Sand_Calc=0;LM=0;Return_Stat=0;do_west=0;num_west=0;numy=0;
00119 param=0;Marq_lambda=1.;Update_estimate=1;ncorr=0;corrmap=0,corrvar=0;nrandom=0;nt_vec2=0;AI_sample=0;aireml_called=0;num_glmm=3;};
|
|
|
Definition at line 120 of file glmm.h. References corrmap, nt_vec2, matvec::Model::numterm, release_obs(), release_SinMat(), and matvec::Matrix< double >::resize().
00121 {release_obs();
00122 release_SinMat();
00123 int i;if(varinv){
00124 for(i=0;i<numterm;i++)varinv[i].resize(0,0);delete [] varinv; varinv=0;
00125 }
00126 if (Var){
00127 for(i=0;i<numterm;i++)Var[i].resize(0,0); delete [] Var;Var=0;
00128 }
00129 if(corrvar) {
00130 for(int i=0;i<numterm;i++) {
00131 delete [] corrvar[i];
00132 }
00133 delete corrvar;corrvar=0;;
00134 }
00135 if(corrmap) delete [] corrmap;corrmap=0;
00136 if(nt_vec2) delete [] nt_vec2;nt_vec2=0;
00137 };
|
|
||||||||||||
|
Add the inverse of additive relationship matrix to mixed model quation. Definition at line 474 of file model_blup.cpp. References matvec::Matrix< T >::begin(), matvec::Matrix< double >::begin(), matvec::Vector< T >::begin(), matvec::Vector< double >::begin(), matvec::Session::epsilon, matvec::Individual::father_id(), matvec::Individual::father_inbcoef(), matvec::getlambda(), matvec::ginverse1(), matvec::Model::hmmec, matvec::Individual::id(), matvec::SparseMatrix::insert(), matvec::Model::lng0vec, matvec::Population::member(), matvec::Model::mme_times_res, matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), matvec::Model::nt_vec, matvec::Model::numgroup, matvec::Model::pop, matvec::Model::popsize, matvec::SESSION, matvec::Model::term, and matvec::Model::var_link. Referenced by matvec::Model::add_G_1(), and matvec::Model::mme_times().
00475 {
00476 unsigned t1, t2, k, i, j, ii, jj, iii, jjj;
00477 double dii,val,value;
00478 double *x,*r;
00479 r = mme_times_res.begin() - 1;
00480 if (x_p) x = x_p->begin() - 1;
00481
00482 double **var = 0;
00483 doubleMatrix *tmp = term[var_link[t]].prior->var_matrix();
00484 if (!tmp) throw exception("Model::add_Ag(): you probably forgot to use M.variance(...)\n"
00485 " to set variance for each random effect\n");
00486 var = tmp->begin();
00487 int nt=nt_vec[var_link[t]];
00488 Matrix<double> rvarg(nt,nt);
00489 for (i=0; i<nt;i++) for (j=0; j<nt; j++) rvarg[i][j]= var[i][j];
00490 ginverse1(rvarg.begin(),nt,lng0vec[t],1,SESSION.epsilon);
00491
00492 Matrix<double> lambda(3,3);
00493 getlambda(lambda.begin(),3);
00494 unsigned startaddr = term[t].start + 1;
00495 unsigned asd[3];
00496 Individual *I;
00497 unsigned nanim=popsize-numgroup;
00498 double f1,f2;
00499 for (k=0; k<nanim; k++) {
00500 I = pop->member(k);
00501 // cout << popsize << " " << " "<<nanim << " " << k << I << "\n";
00502 //cout << I->id() << " " << I->father_id() << " " << I->mother_id() << "\n";
00503 asd[0] = I->id(); asd[1] = I->father_id(); asd[2] = I->mother_id();
00504 f1= I->father_inbcoef();
00505 f2= I->mother_inbcoef();
00506
00507 if(asd[1]> nanim) f1=-1.;
00508 if(asd[2]> nanim) f2=-1.;
00509 dii = 4.0/(2.0 -f1 -f2);
00510
00511 for (i=0; i<3; i++) {
00512 if (asd[i] == 0) continue;
00513 ii = startaddr + (asd[i]-1)*nt;
00514 for (j=0; j<3; j++) {
00515 if (asd[j] == 0) continue;
00516 jj = startaddr + (asd[j]-1)*nt;
00517 val = lambda[i][j]* dii;
00518 for (t1=0; t1<nt; t1++) {
00519 iii = ii + t1;
00520 for (t2=0; t2<nt; t2++) {
00521 jjj = jj + t2;
00522 if (jjj>=iii) {
00523 value = val*rvarg[t1][t2];
00524 if (x_p) { // iteration on data
00525 r[iii] += value*x[jjj];
00526 if (jjj > iii) r[jjj] += value*x[iii];
00527 } else {
00528 hmmec.insert(iii,jjj,value);
00529 }
00530 }
00531 }
00532 }
00533 }
00534 }
00535 }
00536 }
|
|
||||||||||||
|
|
|
||||||||||||
|
Definition at line 985 of file glmm2.cpp. References matvec::Matrix< T >::begin(), matvec::Matrix< double >::begin(), corrvar, EPSILON, matvec::Individual::father_id(), matvec::Individual::father_inbcoef(), matvec::getlambda(), matvec::ginverse1(), matvec::Model::hmmec, matvec::Individual::id(), matvec::SparseMatrix::insert(), matvec::Model::lng0vec, matvec::Population::member(), matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), matvec::Model::numtrait, matvec::Model::pop, matvec::Model::popsize, and matvec::Model::term.
00986 {
00987 unsigned t1, t2, k, i, j, ii, jj, iii, jjj;
00988 double dii,val,value;
00989 double **var = (double **)NULL;
00990 var = corrvar[t][ct].begin();
00991 Matrix<double> rvarg(numtrait,numtrait);
00992 for (i=0; i<numtrait;i++) for (j=0; j<numtrait; j++) rvarg[i][j]= var[i][j];
00993 ginverse1(rvarg.begin(),numtrait,lng0vec[t],1,EPSILON);
00994
00995 Matrix<double> lambda(3,3);
00996 getlambda(lambda.begin(),3);
00997 unsigned startaddr = term[t].start + 1;
00998 unsigned startaddrct = term[ct].start + 1;
00999 unsigned asd[3];
01000 Individual *I;
01001 if(startaddr > startaddrct){
01002 jjj=startaddr;
01003 startaddr =startaddrct;
01004 startaddrct=jjj;
01005 }
01006 for (k=0; k<popsize; k++) {
01007 I = pop->member(k);
01008 asd[0] = I->id(); asd[1] = I->father_id(); asd[2] = I->mother_id();
01009 dii = 4.0/(2.0 - I->father_inbcoef() - I->mother_inbcoef());
01010
01011 for (i=0; i<3; i++) {
01012 if (asd[i] != 0) {
01013 ii = startaddr + (asd[i]-1)*numtrait;
01014 for (j=0; j<3; j++) {
01015 if (asd[j] != 0) {
01016 jj = startaddrct + (asd[j]-1)*numtrait;
01017 val = lambda[i][j]* dii;
01018 for (t1=0; t1<numtrait; t1++) {
01019 iii = ii + t1;
01020 for (t2=0; t2<numtrait; t2++) {
01021 jjj = jj + t2;
01022 if (jjj>=iii) {
01023 value = val*rvarg[t1][t2];
01024 hmmec.insert(iii,jjj,value);
01025 }
01026 }
01027 }
01028 }
01029 }
01030 }
01031 }
01032 }
01033 }
|
|
||||||||||||
|
add the diagonals of inverse of additive relationship for pccg Definition at line 542 of file model_blup.cpp. References matvec::Matrix< T >::begin(), matvec::Matrix< double >::begin(), matvec::Session::epsilon, matvec::Individual::father_id(), matvec::Individual::father_inbcoef(), matvec::getlambda(), matvec::ginverse1(), matvec::Individual::id(), matvec::lidx(), matvec::Model::lng0vec, matvec::Population::member(), matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), matvec::Model::nt_vec, matvec::Model::numgroup, matvec::Model::pop, matvec::Model::popsize, matvec::SESSION, matvec::Model::term, and matvec::Model::var_link. Referenced by matvec::Model::add_G_1_diag().
00543 {
00544 unsigned t1, t2,i,j,k,kk,rec;
00545 double dii,val,value;
00546 double **var = 0;
00547 doubleMatrix *tmp = term[var_link[t]].prior->var_matrix();
00548 if (!tmp) throw exception("Model::add_Ag(): you probably forgot to use M.variance(...)\n"
00549 " to set variance for each random effect\n");
00550 var = tmp->begin();
00551 int nt=nt_vec[var_link[t]];
00552 Matrix<double> rvarg(nt,nt);
00553 for (i=0; i<nt;i++) for (j=0; j<nt; j++) rvarg[i][j]= var[i][j];
00554 ginverse1(rvarg.begin(),nt,lng0vec[t],1,SESSION.epsilon);
00555
00556 Matrix<double> lambda(3,3);
00557 getlambda(lambda.begin(),3);
00558 unsigned startaddr = term[t].start + 1;
00559 unsigned asd[3];
00560 Individual *I;
00561 unsigned nanim=popsize-numgroup;
00562 double f1,f2;
00563
00564 for (k=0,i=0; i<t; ++i) k += term[i].nlevel();
00565 for (rec=0; rec<nanim; rec++) {
00566 I = pop->member(rec);
00567 asd[0] = I->id(); asd[1] = I->father_id(); asd[2] = I->mother_id();
00568 f1= I->father_inbcoef();
00569 f2= I->mother_inbcoef();
00570
00571 if(asd[1]> nanim) f1=-1.;
00572 if(asd[2]> nanim) f2=-1.;
00573 dii = 4.0/(2.0 -f1 -f2);
00574
00575 for (i=0; i<3; i++) {
00576 if (asd[i] == 0) continue;
00577 kk = k + asd[i] - 1;
00578 val = lambda[i][i]* dii;
00579 for (t1=0; t1<nt; ++t1) {
00580 for (t2=t1; t2<nt; ++t2) {
00581 M[kk][lidx(t2,t1,nt)+1] += val*rvarg[t1][t2];
00582 }
00583 }
00584 }
00585 }
00586 }
|
|
|
Definition at line 1205 of file glmm2.cpp. References matvec::Model::act_numtrait, matvec::Matrix< T >::begin(), matvec::Matrix< double >::begin(), build_CorrVar(), CorrVar, EPSILON, matvec::Individual::father_id(), matvec::Individual::father_inbcoef(), matvec::getlambda(), matvec::ginverse1(), matvec::Model::hmmec, matvec::Individual::id(), matvec::SparseMatrix::insert(), matvec::Model::lng0vec, matvec::Population::member(), matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), matvec::Model::numgroup, matvec::Model::numterm, matvec::Model::numtrait, matvec::Model::pop, matvec::Model::popsize, matvec::Vector< T >::resize(), and matvec::Model::term.
01206 {
01207 unsigned t1, t2, k, i, j, ii, jj, iii, jjj,icnt;
01208 double dii,val,value;
01209 double **var = (double **)NULL;
01210 Vector<int> tmap;
01211 build_CorrVar();
01212 tmap.resize(act_numtrait);
01213 icnt=0;
01214 for(i=t;i<numterm;i++){
01215 if(term[i].classi()=='P') {
01216 for(j=0;j<numtrait;j++,icnt++){
01217 tmap[icnt]=i;
01218 }
01219 }
01220 }
01221 // std::cout <<tmap;
01222 doubleMatrix *tmp = &CorrVar;
01223 if (!tmp) throw exception("Model::add_Ag(): you probably forgot to use M.variance(...)\n"
01224 " to set variance for each random effect");
01225 var = tmp->begin();
01226 Matrix<double> rvarg(act_numtrait,act_numtrait);
01227 for (i=0; i<act_numtrait;i++) for (j=0; j<act_numtrait; j++) rvarg[i][j]= var[i][j];
01228
01229 ginverse1(rvarg.begin(),act_numtrait,lng0vec[t],1,EPSILON);
01230
01231
01232 double f1,f2;
01233 Matrix<double> lambda(3,3);
01234 getlambda(lambda.begin(),3);
01235 unsigned startaddr = term[t].start + 1;
01236 unsigned asd[3];
01237 Individual *I;
01238 unsigned nanim=popsize-numgroup;
01239 for (k=0; k<nanim; k++) {
01240 I = pop->member(k);
01241 asd[0] = I->id(); asd[1] = I->father_id(); asd[2] = I->mother_id();
01242 f1= I->father_inbcoef() ;
01243 f2= I->mother_inbcoef();
01244 if(asd[1]> nanim) f1=-1.;
01245 if(asd[2]> nanim) f2=-1.;
01246 dii = 4.0/(2.0 -f1-f2);
01247
01248 for (i=0; i<3; i++) {
01249 if (asd[i] != 0) {
01250 ii = startaddr + (asd[i]-1)*act_numtrait;
01251 for (j=0; j<3; j++) {
01252 if (asd[j] != 0) {
01253 jj = startaddr + (asd[j]-1)*act_numtrait;
01254 val = lambda[i][j]* dii;
01255 for (t1=0; t1<act_numtrait; t1++) {
01256 iii = ii + t1;
01257 for (t2=0; t2<act_numtrait; t2++) {
01258 jjj = jj + t2;
01259 if (jjj>=iii) {
01260 value = val*rvarg[t1][t2];
01261 //if(term[tmap[t1]].trait[t1%numtrait] && term[tmap[t2]].trait[t2%numtrait])
01262 hmmec.insert(iii,jjj,value);
01263 }
01264 }
01265 }
01266 }
01267 }
01268 }
01269 }
01270 }
01271 }
|
|
|
Definition at line 1036 of file glmm2.cpp. References matvec::Matrix< T >::begin(), matvec::Matrix< double >::begin(), EPSILON, matvec::Individual::father_id(), matvec::Individual::father_inbcoef(), matvec::getlambda(), matvec::ginverse1(), hSand, matvec::Individual::id(), matvec::SparseMatrix::insert(), matvec::Model::lng0vec, matvec::Population::member(), matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), matvec::Model::numtrait, matvec::Model::pop, matvec::Model::popsize, and matvec::Model::term. Referenced by add_G_Sand().
01037 {
01038 unsigned t1, t2, k, i, j, ii, jj, iii, jjj;
01039 double dii,val,value;
01040 double **var = (double **)NULL;
01041 doubleMatrix *tmp = term[t].prior->var_matrix();
01042 if (!tmp) throw exception("Model::add_Ag(): you probably forgot to use M.variance(...)\n"
01043 " to set variance for each random effect");
01044 var = tmp->begin();
01045 Matrix<double> rvarg(numtrait,numtrait);
01046 for (i=0; i<numtrait;i++) for (j=0; j<numtrait; j++) rvarg[i][j]= var[i][j];
01047 ginverse1(rvarg.begin(),numtrait,lng0vec[t],1,EPSILON);
01048
01049 Matrix<double> lambda(3,3);
01050 getlambda(lambda.begin(),3);
01051 unsigned startaddr = term[t].start + 1;
01052 unsigned asd[3];
01053 Individual *I;
01054
01055 for (k=0; k<popsize; k++) {
01056 I = pop->member(k);
01057 asd[0] = I->id(); asd[1] = I->father_id(); asd[2] = I->mother_id();
01058 dii = 4.0/(2.0 - I->father_inbcoef() - I->mother_inbcoef());
01059
01060 for (i=0; i<3; i++) {
01061 if (asd[i] != 0) {
01062 ii = startaddr + (asd[i]-1)*numtrait;
01063 for (j=0; j<3; j++) {
01064 if (asd[j] != 0) {
01065 jj = startaddr + (asd[j]-1)*numtrait;
01066 val = lambda[i][j]* dii;
01067 for (t1=0; t1<numtrait; t1++) {
01068 iii = ii + t1;
01069 for (t2=0; t2<numtrait; t2++) {
01070 jjj = jj + t2;
01071 if (jjj>=iii) {
01072 value = val*rvarg[t1][t2];
01073 hSand.insert(iii,jjj,value);
01074 }
01075 }
01076 }
01077 }
01078 }
01079 }
01080 }
01081 }
01082 }
|
|
|
Add G-inverse Definition at line 279 of file model_blup.cpp. References matvec::Model::add_Ag(), matvec::Model::add_Ig(), matvec::Model::numterm, and matvec::Model::term. Referenced by matvec::Model::setup_mme().
|
|
|
|
|
|
add diagonals of G-inverse for pccg Definition at line 297 of file model_blup.cpp. References matvec::Model::add_Ag_diag(), matvec::Model::add_Ig_diag(), matvec::Model::numterm, and matvec::Model::term. Referenced by matvec::Model::compute_rhs_diag_mme().
00298 {
00299 int done_ped=0;
00300 for (int t=0; t<numterm; t++) {
00301 if (term[t].classi() == 'P') {
00302 add_Ag_diag(t,M);
00303 //hmmec.display(12,24,12,24);
00304 }
00305 else if (term[t].classi() == 'R') {
00306 add_Ig_diag(t,M);
00307 }
00308 }
00309 }
|
|
|
Definition at line 1190 of file glmm2.cpp. References matvec::Model::numterm, and matvec::Model::term.
|
|
|
Definition at line 332 of file model_vce.cpp. References matvec::Matrix< T >::begin(), matvec::Individual::father_id(), matvec::Individual::father_inbcoef(), matvec::getlambda(), matvec::Model::hmmec, matvec::Individual::id(), matvec::SparseMatrix::insert(), matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), matvec::Model::numterm, matvec::Model::pop, matvec::Population::popmember, and matvec::Model::term. Referenced by matvec::Model::setup_ww_single_trait(), and matvec::Model::vce_emreml_multi_trait().
00333 {
00334 unsigned i,j,ii,jj,t,k,rec,nl,startaddr;
00335 double val,rr,xval;
00336 for (k=0,t=0; t<numterm; t++) {
00337 startaddr = term[t].startaddr();
00338 nl = term[t].nlevel();
00339 if (term[t].classi() == 'P') {
00340 Individual *I;
00341 unsigned asd[3];
00342 val = ratio[k++];
00343 Matrix<double> lambda(3,3);
00344 getlambda(lambda.begin(),3);
00345 for (rec=0; rec<nl; rec++) {
00346 I = pop->popmember[rec];
00347 asd[0]=I->id(); asd[1]=I->father_id(); asd[2] = I->mother_id();
00348 rr = val*4.0/(2.0-I->father_inbcoef()- I->mother_inbcoef());
00349 for (i=0; i<3; i++) {
00350 if (asd[i] != 0) {
00351 ii = startaddr + asd[i];
00352 for (j=0; j<3; j++) {
00353 if (asd[j] != 0) {
00354 jj = startaddr + asd[j];
00355 if (jj >= ii) {
00356 xval = lambda[i][j]*rr;
00357 hmmec.insert(ii,jj,xval);
00358 }
00359 }
00360 }
00361 }
00362 }
00363 }
00364 }
00365 else if (term[t].classi() == 'R') {
00366 val = ratio[k++];
00367 for (ii=startaddr+1,i=0; i<nl; i++,ii++) hmmec.insert(ii,ii,val);
00368 }
00369 }
00370 }
|
|
|
Definition at line 942 of file glmm2.cpp. References add_AgSand(), add_IgSand(), matvec::Model::numterm, and matvec::Model::term. Referenced by setup_mme().
00943 {
00944 for (int t=0; t<numterm; t++) {
00945 if (term[t].classi() == 'P') {
00946 add_AgSand(t);
00947 }
00948 else if (term[t].classi() == 'R') {
00949 add_IgSand(t);
00950 }
00951 }
00952 }
|
|
||||||||||||
|
Definition at line 399 of file model_blup.cpp. References matvec::Matrix< T >::begin(), matvec::Matrix< double >::begin(), matvec::Vector< T >::begin(), matvec::Vector< double >::begin(), matvec::Session::epsilon, matvec::ginverse1(), matvec::Model::hmmec, matvec::SparseMatrix::insert(), matvec::Model::lng0vec, matvec::Model::mme_times_res, matvec::Model::nt_vec, matvec::SESSION, matvec::Model::term, and matvec::Model::var_link. Referenced by matvec::Model::add_G_1(), and matvec::Model::mme_times().
00400 {
00401 unsigned k,i,ii,jj,nt;
00402 int t1,t2;
00403 double *x,*r,value;
00404 r = mme_times_res.begin() - 1;
00405 if (x_p) x = x_p->begin() - 1;
00406
00407 double **v = 0;
00408 doubleMatrix *tmp = term[var_link[t]].prior->var_matrix();
00409 if (!tmp) throw exception("Model::add_Ig(): you probably forgot to use M.variance(...)\n"
00410 " to set variance for each random effect\n");
00411 nt=nt_vec[var_link[t]];
00412 v = tmp->begin();
00413 Matrix<double> rv(nt,nt);
00414 for (t1=0; t1<nt; t1++) {
00415 for(t2=0; t2<nt; t2++) rv[t1][t2]=v[t1][t2];
00416 }
00417 ginverse1(rv.begin(),nt,lng0vec[t],1,SESSION.epsilon);
00418
00419 unsigned na = term[t].nlevel();
00420 unsigned startaddr = term[t].start + 1;
00421 for (k=0; k<na; k++) {
00422 i = startaddr + k*nt;
00423 for (t1=0; t1<nt; t1++) {
00424 ii = i + t1;
00425 for (t2=0; t2<nt; t2++) {
00426 jj = i + t2;
00427 if (jj>=ii) {
00428 value = rv[t1][t2];
00429 if(x_p) {// iteration on data
00430 r[ii] += value*x[jj];
00431 if (jj > ii) r[jj] += value*x[ii];
00432 } else {
00433 hmmec.insert(ii,jj,value);
00434 }
00435 }
00436 }
00437 }
00438 }
00439 }
|
|
||||||||||||
|
|
|
||||||||||||
|
Definition at line 441 of file model_blup.cpp. References matvec::Matrix< T >::begin(), matvec::Matrix< double >::begin(), matvec::Session::epsilon, matvec::ginverse1(), matvec::lidx(), matvec::Model::lng0vec, matvec::Model::nt_vec, matvec::SESSION, matvec::Model::term, and matvec::Model::var_link. Referenced by matvec::Model::add_G_1_diag().
00442 {
00443 unsigned rec,k,i,ii,jj,kk,nt;
00444 int t1,t2;
00445 double **v = 0;
00446 doubleMatrix *tmp = term[var_link[t]].prior->var_matrix();
00447 if (!tmp) throw exception("Model::add_Ig(): you probably forgot to use M.variance(...)\n"
00448 " to set variance for each random effect\n");
00449 nt=nt_vec[var_link[t]];
00450 v = tmp->begin();
00451 Matrix<double> rv(nt,nt);
00452 for (t1=0; t1<nt; t1++) {
00453 for(t2=0; t2<nt; t2++) rv[t1][t2]=v[t1][t2];
00454 }
00455 ginverse1(rv.begin(),nt,lng0vec[t],1,SESSION.epsilon);
00456
00457 unsigned na = term[t].nlevel();
00458 for (k=0,i=0; i<t; ++i) k += term[i].nlevel();
00459 for (rec=0; rec<na; rec++) {
00460 kk = k + rec;
00461 for (t1=0; t1<nt; ++t1) {
00462 for (t2=t1; t2<nt; ++t2) {
00463 M[kk][lidx(t2,t1,nt)+1] += rv[t1][t2];
00464 }
00465 }
00466 }
00467 }
|
|
|
Definition at line 954 of file glmm2.cpp. References matvec::Matrix< T >::begin(), matvec::Matrix< double >::begin(), EPSILON, matvec::ginverse1(), hSand, matvec::SparseMatrix::insert(), matvec::Model::lng0vec, matvec::Model::numtrait, and matvec::Model::term. Referenced by add_G_Sand().
00955 {
00956 unsigned k,i,ii,jj;
00957 int t1,t2;
00958 double **v = (double **)NULL;
00959 doubleMatrix *tmp = term[t].prior->var_matrix();
00960 if (!tmp) throw exception("Model::add_Ig(): you probably forgot to use M.variance(...)\n"
00961 " to set variance for each random effect");
00962 v = tmp->begin();
00963 Matrix<double> rv(numtrait,numtrait);
00964 for (t1=0; t1<numtrait; t1++) {
00965 for(t2=0; t2<numtrait; t2++) rv[t1][t2]=v[t1][t2];
00966 }
00967 ginverse1(rv.begin(),numtrait,lng0vec[t],1,EPSILON);
00968
00969 unsigned na = term[t].nlevel();
00970 unsigned startaddr = term[t].start + 1;
00971 for (k=0; k<na; k++) {
00972 i = startaddr + k*numtrait;
00973 for (t1=0; t1<numtrait; t1++) {
00974 ii = i + t1;
00975 for (t2=0; t2<numtrait; t2++) {
00976 jj = i + t2;
00977 if (jj>=ii) hSand.insert(ii,jj,rv[t1][t2]);
00978 }
00979 }
00980 }
00981 }
|
|
||||||||||||||||
|
Definition at line 1017 of file glmm1.cpp. References aireml_called, matvec::Matrix< double >::assign(), matvec::Vector< T >::begin(), DO_LOG, matvec::doubleMatrix::eigen(), matvec::Session::epsilon, EPSILON, matvec::doubleMatrix::ginv0(), glim(), matvec::hadjoin(), matvec::doubleMatrix::identity(), Info, matvec::Vector< T >::max(), nt_vec2, num_glmm, matvec::Model::numrec, matvec::Model::numtrait, observation, matvec::Session::output_precision, Print_Level, matvec::ranf(), res_nvc, matvec::Vector< T >::resize(), matvec::Matrix< double >::resize(), matvec::Model::restricted_log_likelihood(), Score, matvec::SESSION, SSQCP(), TotalNvc(), matvec::Matrix< double >::transpose(), Update_estimate, Var2Vec(), varest, and Vec2Var().
01018 {
01019 int Converged,Need_Like;
01020 aireml_called=-1;
01021 Converged=0;
01022 Need_Like=1;
01023 Vector<double> varold,varnew;
01024 int Update_estimate_AI;
01025 int nvc,iteration,numran,k,t1,t2,raneff,kidx,kend;
01026 Update_estimate=1;
01027 // info_scale=1.1;
01028 glim(num_glmm);
01029 // glim(5)->display();
01030 Update_estimate=1;
01031 nvc=TotalNvc();
01032 doubleMatrix varbound;
01033 varbound.resize(numtrait,numtrait);
01034 Vector<double> eig,flg,fix;
01035 eig.resize(numtrait);
01036 flg.resize(numtrait);
01037 fix.resize(nvc);
01038 varold.resize(nvc);
01039 varnew.resize(nvc);
01040 numran=Var2Vec(varold.begin()); //Have Numran return Number af random
01041 // effects plus one for P if it exists
01042
01043 doubleMatrix Adj_Cov;
01044 if(numran) P_Mat = new doubleMatrix [numran];
01045 doubleMatrix Ip,Info_orig;
01046 double log_like_old,log_like_new;
01047 if(Need_Like){
01048 log_like_old=restricted_log_likelihood(0);
01049 std::cout << "\nOriginal Residual Log Likelihood:" << std::setprecision(SESSION.output_precision) << log_like_old <<"\n\n";
01050 Need_Like=0;
01051 }
01052 for(k=0;k<numran;k++) P_Mat[k].identity(nt_vec2[k]);
01053 //if(Print_Level == 0) std::cout << "Initial Estimates " << varold;
01054 long subit;
01055 for(iteration=0;(iteration<numiter) && !Converged;iteration++){
01056
01057 Update_estimate=1;
01058 //log_like_old=restricted_log_likelihood(1);
01059 //std::cout << "\nOld Residual Log Likelihood:" << log_like_old <<"\n";
01060 SSQCP();
01061 if(Need_Like) {
01062 //glim(2);
01063 log_like_old=restricted_log_likelihood(0);
01064 Need_Like=0;
01065 }
01066 if(Print_Level > 0) std::cout << "\n Old Residual Log Likelihood:" << std::setprecision(SESSION.output_precision)<< log_like_old <<"\n";
01067 Var2Vec(varold.begin());
01068 subit=0;
01069
01070 info_scale=0.;
01071 Info_orig=Info;
01072 do{
01073 Info=Info_orig;
01074 // std::cout << "\nOrig Info" << Info_orig;
01075 if(info_scale ){
01076 double info_mult=std::exp(info_scale);
01077 //std::cout << "\n Info Mult " << info_mult;
01078 for(k=0;k<nvc;k++) Info[k][k]*=info_mult;
01079 }
01080 //std::cout <<"\n Fix " << fix;
01081 for(k=0;k<nvc;k++) {
01082 if(Info[k][k] < EPSILON || fix[k]) {
01083 for(int k2=0;k2<nvc;k2++) {
01084 Info[k][k2]=0;
01085 Info[k2][k]=0;
01086 }
01087 Info[k][k]=10.*EPSILON;
01088 //fix[k]=1;
01089 }
01090 }
01091 #ifdef DO_CHOL_NOT
01092 //varbound=Info;
01093 // std::cout << "Eigen" << varbound.eigen();
01094 /*
01095 for(k=0;k<nvc;k++) {
01096 if(fix[k] ) {
01097 Score[k]=0;
01098 for(int k2=0;k2<nvc;k2++) {
01099 Info[k2][k]=0;
01100 Info[k][k2]=0;
01101 }
01102 Info[k][k]=1.e-10
01103 }
01104 }
01105 */
01106 // Info_orig=Info;
01107 //std::cout << "Score " << Score;
01108 //std::cout << "Info " << Info;
01109 // std::cout << "ascov 1" << Info_orig.ginv1();
01110 // std::cout << "ascov 0" << Info_orig.ginv0();
01111 /*Info_orig=Info;*/
01112 #endif
01113 // std::cout << "Score " << Score;
01114 // std::cout << "Info " << Info;
01115 Vector<double> Delt;
01116 Delt=Info.ginv0()*Score;
01117 varnew=varold+Delt;
01118 //std::cout << "\nVar " << varnew << varold << Delt;
01119 // if(Print_Level >= 0) {
01120 // std::cout << " New Estimates before Adjustment" << varnew;
01121 //}
01122 #ifdef DO_LOG
01123 doubleMatrix P,L,D;
01124
01125 Adj_Cov.identity(nvc);
01126 int offset=nvc-res_nvc;
01127 for(k=0;k<res_nvc;k++)
01128 for(int k2=0;k2<res_nvc;k2++)
01129 Adj_Cov[offset+k][offset+k2]=observation[numrec].Adj_Cov[k][k2];
01130
01131
01132 Vector<double> d;
01133 int nt;
01134 for(k=0,raneff=0;raneff<numran;raneff++) {
01135 nt=nt_vec2[raneff];
01136 varbound.resize(nt,nt);
01137 for(kidx=k,t1=0;t1<nt;t1++) for(t2=t1;t2<nt;t2++,kidx++) {
01138 varbound[t1][t2]=varold[kidx];
01139 varbound[t2][t1]=varold[kidx];
01140 }
01141 // std::cout << "varbound "<< varbound;
01142 P=varbound;
01143 //std::cout << "\nP " << P;
01144 d=P.eigen();
01145 D.resize(nt,nt);
01146 for(t1=0;t1<nt;t1++){
01147 d[t1]=std::log(d[t1]);
01148 D[t1][t1]=d[t1];
01149 }
01150
01151
01152 doubleMatrix part;
01153
01154 for(kidx=k,t1=0;t1<nt;t1++) for(t2=t1;t2<nt;t2++,kidx++) {
01155 part.resize(nt,nt);
01156 double D_avg,D_delt;
01157 D_avg=.5*(d[t1]+d[t2]);
01158 D_delt=.5*(d[t1]-d[t2]);
01159 if(fabs(D_delt) > 1.e-8) {
01160
01161 //std::cout << "how did I get here!!!\n";
01162 part[t1][t2]=std::exp(D_avg)*(sinh(D_delt))/(D_delt);
01163 part[t2][t1]=std::exp(D_avg)*(sinh(D_delt))/(D_delt);
01164 }
01165 else{
01166 part[t1][t2]=std::exp(D_avg);
01167 part[t2][t1]=std::exp(D_avg);
01168 }
01169 //std::cout << "part\n" << part <<"d\n"<< d;
01170 part=P*part*P.transpose();
01171 if(d[t1] <= -9. && d[t2] <= -9.) part.assign(0.0);
01172 int kkidx,tt1,tt2;
01173 for( kkidx=k,tt1=0;tt1<nt;tt1++) for(tt2=tt1;tt2<nt;tt2++,kkidx++) {
01174 // std::cout <<"\n" << kkidx << " " << kidx << " " << tt1 << " " << tt2 << " " << nt << " " << raneff << " " << k <<"\n",
01175 Adj_Cov[kkidx][kidx]=part[tt1][tt2];
01176 }
01177 }
01178
01179 /*
01180 double dmax=0.;
01181 for(kidx=k,t1=0;t1<nt;t1++) {
01182 for(t2=t1;t2<nt;t2++,kidx++) {
01183 if(fabs(Delt[kidx])> dmax) dmax=fabs(Delt[kidx]);
01184 }
01185 }
01186
01187
01188 #define MAXCHG 1.
01189 if(dmax > MAXCHG){
01190 std::cout << "\ndmax " << dmax <<"\n";
01191 for(kidx=k,t1=0;t1<nt;t1++) {
01192 for(t2=t1;t2<nt;t2++,kidx++) {
01193 Delt[kidx]/=dmax/MAXCHG;
01194 }
01195 }
01196 }
01197 */
01198 //std::cout << "D" << D << "P" << P;
01199 for(kidx=k,t1=0;t1<nt;t1++) {
01200 for(t2=t1;t2<nt;t2++,kidx++) {
01201 D[t1][t2]+=Delt[kidx];
01202 D[t2][t1]=D[t1][t2];
01203 }
01204 }
01205
01206 //std::cout << "D new" << D;
01207
01208 D=P*D*P.transpose();
01209 doubleMatrix Q;
01210 doubleMatrix DE;
01211 Q=D;
01212 d=Q.eigen();
01213 //std::cout << " d " << d << " Q " << Q << " D" << D;
01214 double dmax=d.max();
01215 //if(dmax < -8.) std::cerr << "\nD\n" << D << "\nP\n" << P<< "\n";
01216 if(dmax > 15.) dmax=15.;//
01217 if(dmax < 0.) dmax=0;
01218 double dmin=std::log(10.*SESSION.epsilon)+dmax;
01219 for(t1=0;t1<nt;t1++) {
01220 if(d[t1] < dmin) d[t1]=dmin+std::log(.5);
01221 if(d[t1] > dmax) d[t1]=15;
01222 }
01223 DE.resize(nt,nt);
01224 for(t1=0;t1<nt;t1++)DE[t1][t1]=std::exp(d[t1]);
01225 D=Q*DE*Q.transpose();
01226 Q=D;
01227 d=Q.eigen();
01228 //std::cout << " d exp " << d << " Q " << Q;
01229 varbound=D;//P*D*P.t();
01230 //std::cout << " d " << d;
01231 // std::cout << " New Estimates before Adjustment" << varbound;
01232 for(t1=0;t1<nt;t1++) {
01233 // varbound[t1][t1]+=1.0e-4;
01234 if(varbound[t1][t1] < std::exp(-11.)) {
01235 for(t2=0;t2<nt;t2++){
01236 varbound[t1][t2]=0.;
01237 varbound[t2][t1]=0.;
01238 }
01239 varbound[t1][t1]=std::exp(-11.5);
01240 }
01241 }
01242 //std::cout << "varbound new" << varbound;
01243 for(t1=0;t1<nt;t1++) for(t2=t1;t2<nt;t2++,k++) {
01244 varnew[k]=varbound[t1][t2];
01245 }
01246 }
01247
01248 #endif
01249
01250
01251 // ****
01252 // If the variance covariance matrix is not pd make it slighty pd
01253 // ****
01254 //Print_Level=2;
01255 #ifdef DO_CHOL_2
01256 int nt;
01257 for(k=0,raneff=0;raneff<numran;raneff++) {
01258 nt=nt_vec2[raneff];
01259 flg.zeros(nt);
01260
01261 for(t1=0;t1<nt;t1++) for(t2=t1;t2<nt;t2++,k++) {
01262 if(t1==t2 && varnew[k]<1.e-8 && !fix[k]) {
01263 Score[k]=0;;
01264 varnew[k]=1.e-4;
01265 flg[t1]=1;
01266 fix[k]=1;
01267 }
01268 else if(flg[t1] && t1!=t2) {
01269 Score[k]=0.;
01270 varnew[k]=0.;
01271 fix[k]=1;
01272 }
01273 }
01274 }
01275
01276 // for(k=0;k<nvc;k++) {
01277 // if(fix[k]){
01278 // for(int k2=0;k2<nvc;k2++) {
01279 // Info[k][k2]=0;
01280 // Info[k2][k]=0;
01281 // }
01282 // Info[k][k]=1.e-4;
01283 // }
01284 // }
01285 // varnew=varold+Info.ginv0()*Score;
01286 #endif
01287
01288 #ifdef DO_PMAT
01289 int nt;
01290 for(k=0,raneff=0;raneff<numran;raneff++) {
01291 nt=nt_vec2[raneff];
01292 varbound.resize(nt,nt);
01293 eig.resize(nt);
01294 flg.resize(nt);
01295 for(kidx=k,t1=0;t1<nt;t1++) for(t2=t1;t2<nt;t2++,kidx++) {
01296 varbound[t1][t2]=varnew[kidx];
01297 varbound[t2][t1]=varnew[kidx];
01298 }
01299
01300 for(t1=0;t1<nt;t1++){
01301 if(varbound[t1][t1] < 1.e-8){
01302 for(t2=0;t2<nt;t2++) {
01303 varbound[t1][t2]=0;
01304 varbound[t2][t1]=0;
01305 }
01306 }
01307 }
01308
01309 // P_Mat[raneff]=P_Mat[raneff]*P_Mat[raneff].transpose()*varbound*P_Mat[raneff]*P_Mat[raneff].transpose();
01310
01311
01312 P_Mat[raneff]=varbound;
01313
01314 eig=P_Mat[raneff].eigen();
01315 if(Print_Level > 1) {
01316 // std::cout << "varbound" << varbound << "\n";
01317 // std::cout << "Eigen Values " << eig << "\n";
01318 // std::cout << "P_Mat" << P_Mat[raneff] << "\n";
01319 }
01320
01321 for(t1=0;t1<nt;t1++){
01322 flg[t1]=0;
01323 if(eig[t1]<1.e-8){
01324 eig[t1]=.99e-3;
01325 flg[t1]=1;
01326 }
01327 if(eig[t1]>1.e10){
01328 eig[t1]=1.01e10;
01329 flg[t1]=0;
01330 }
01331 }
01332 if(Print_Level > 1) {
01333 // std::cout << "Adjusted Eigen Values " << eig << "\n";
01334 }
01335 varbound=P_Mat[raneff]*eig.diag()*P_Mat[raneff].transpose();
01336 for(t1=0;t1<nt;t1++){
01337 if(flg[t1]){
01338 for(t2=0;t2<nt;t2++) {
01339 P_Mat[raneff][t1][t2]=0;
01340 }
01341 }
01342 }
01343
01344 for(t1=0;t1<nt;t1++) for(t2=t1;t2<nt;t2++,k++) {
01345 varnew[k]=varbound[t1][t2];
01346
01347 }
01348 eig=varbound.eigen();
01349 if(Print_Level > 1) {
01350 // std::cout << " Eigen Values bounded Matrix " << eig << "\n";
01351 }
01352
01353 }
01354 #endif /* DO_PMAT */
01355 Vec2Var(varnew.begin());
01356 //std::cout << "Varnew " <<varnew;
01357 glim(num_glmm);
01358 log_like_new=restricted_log_likelihood(0);
01359 std::cout << " Iteration " << (iteration+1) <<"."<< subit << " Res Log Like "<< std::setprecision(SESSION.output_precision)<< log_like_new << " Change " << std::setprecision(SESSION.output_precision)<<(log_like_new-log_like_old) << "\n";
01360 info_scale+=1.5*ranf();
01361 }while(tol && log_like_new-log_like_old < -tol && subit++<= 10);
01362 if(subit == 0 && std::abs(log_like_new-log_like_old) < std::max(tol,EPSILON) ) Converged=1;
01363 if(Converged) aireml_called=1;
01364 log_like_old=log_like_new;
01365 info_scale-=.75;
01366 varold=varnew;
01367
01368 // if(iteration%3 == 2) {
01369 // Update_estimate=1;
01370 // glim(2);
01371 // Update_estimate=1;
01372 // }
01373
01374 // Print_Level=2;
01375
01376 if(Print_Level > 0) {
01377 std::cout << " Iteration " << iteration << " Subit " << subit << "\n";
01378 std::cout << " Residual log likelihood "
01379 << std::setprecision(SESSION.output_precision)
01380 << log_like_new << "\n";
01381 std::cout << " New Estimates \n" << varnew;
01382 if(Print_Level > 1) {
01383 std::cout << " Score \n" << Score;
01384 std::cout << " Info \n" << Info;
01385 std::cout << " Asy Cov \n" << Info.ginv0();
01386 }
01387 }
01388
01389
01390 Vec2Var(varnew.begin());
01391 varold=varnew;
01392 }
01393 if(Print_Level >= 0) {
01394 if(Converged){
01395 std::cout << "\n Iteration " << iteration << " Converged\n";
01396 }
01397 else{
01398 std::cout << "\n Iteration " << iteration << " Subit " << subit << " Failed to Converge\n";
01399 }
01400 std::cout << " Final Estimates \n" << varnew;
01401 std::cout << " Last Change \n" << Info.ginv0()*Score;
01402 Info=Info_orig;
01403 if(info_scale) std::cout << " Unscaled Last Change \n" << Info.ginv0()*Score;
01404 std::cout << " Residual log likelihood "
01405 << std::setprecision(SESSION.output_precision)
01406 << restricted_log_likelihood(0) << "\n" ;
01407 Info=Info_orig;
01408 // std::cout << " Adj Matrix " << Adj_Cov;
01409 std::cout << " Asy Covariance Matrix\n" <<
01410 #ifdef DO_LOG
01411 Adj_Cov*
01412 #endif
01413 Info.ginv0()
01414 #ifdef DO_LOG
01415 *Adj_Cov.transpose()
01416 #endif
01417 ;
01418 }
01419 Update_estimate=1;
01420 //glim();
01421 // for(k=0;k<numran;k++) P_Mat[k].resize(0,0);
01422 if(P_Mat)delete [] P_Mat;
01423 P_Mat=NULL;
01424 varest=varnew;
01425 doubleMatrix ans;
01426 ans = hadjoin(varest,
01427 #ifdef DO_LOG
01428 Adj_Cov*
01429 #endif
01430 Info.ginv0()
01431 #ifdef DO_LOG
01432 *Adj_Cov.transpose()
01433 #endif
01434 );
01435 return ans;
01436 }
|
|
|
Definition at line 147 of file glmm.h. References AI_sample.
00147 {if(S_set != -1) AI_sample=S_set;return(AI_sample);}
|
|
|
Definition at line 966 of file glmm1.cpp. References matvec::Matrix< T >::begin(), matvec::SparseMatrix::close(), matvec::Individual::father_id(), matvec::Individual::father_inbcoef(), matvec::getlambda(), hainv, matvec::Individual::id(), matvec::SparseMatrix::insert(), ldet, matvec::SparseMatrix::logdet(), matvec::Population::member(), matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), Need_Ainv, matvec::Model::numgroup, matvec::Model::pop, matvec::Model::popsize, and matvec::SparseMatrix::resize(). Referenced by contrast(), residual(), and SSQCP().
00967 {
00968 if(!Need_Ainv) return(hainv);
00969 Need_Ainv=0;
00970 //TW int do_west;
00971 //TW Vector<double> West;
00972 //TW do_west=pop->do_west;
00973 //TW West=pop->West;
00974 unsigned nanim=popsize-numgroup;
00975 hainv.resize(popsize,4*popsize);
00976 double dii,val;
00977
00978 Matrix<double> lambda(3,3);
00979 getlambda(lambda.begin(),3);
00980 unsigned asd[3];
00981 int ii,jj,j;
00982 Individual *I;
00983 double f1,f2;
00984 ldet=0.;
00985 for (int k=0; k<nanim; k++) {
00986 I = pop->member(k);
00987 asd[0] = I->id(); asd[1] = I->father_id(); asd[2] = I->mother_id();
00988 dii = 4.0/(2.0 - I->father_inbcoef() - I->mother_inbcoef());;
00989 f1= I->father_inbcoef() ;
00990 f2= I->mother_inbcoef();
00991 if(asd[1]> nanim) f1=-1.;
00992 if(asd[2]> nanim) f2=-1.;
00993 dii = 4.0/(2.0 -f1-f2);
00994 ldet+=std::log(dii);
00995 for (int i=0; i<3; i++) {
00996 if (asd[i] != 0) {
00997 ii = asd[i];
00998 for (j=0; j<3; j++) {
00999 if (asd[j] != 0) {
01000 jj = asd[j];
01001 val = lambda[i][j]* dii;
01002 // std::cout << "\n ii " << ii << " jj " << jj << " val " << val << " dii " << dii;
01003 if(ii <= jj ) hainv.insert(ii,jj,val);
01004 }
01005 }
01006 }
01007 }
01008 }
01009 ldet=-2.*ldet;
01010 //std::cout << "ldet " << ldet << "\n";
01011 //std::cout << hainv.dense();
01012 hainv.logdet();
01013 hainv.close();
01014
01015 return(hainv);
01016 }
|
|
||||||||||||||||||||
|
Definition at line 426 of file model_hypo.cpp. References matvec::Model::blupsol, matvec::Model::hmmec, matvec::Vector< double >::inner_product(), matvec::SparseMatrix::logdet(), matvec::Model::numgroup, matvec::Model::numterm, matvec::Model::numtrait, matvec::Model::rellrhs, matvec::Model::residual_var, matvec::Model::term, and matvec::Model::yry.
00427 {
00428 ////////////////////////////////////////////////////////
00429 // calculates sums of squares for single trait model
00430 // the rank of the random part:
00431 ////////////////////////////////////////////////////////
00432 unsigned rank_mme=0,rank_z=0;
00433 char cl;
00434 for (unsigned i=0; i<numterm; i++) {
00435 cl = term[i].classi();
00436 if (cl == 'P' || cl == 'R') rank_z += term[i].nlevel();
00437 }
00438 rank_z = (rank_z - numgroup)*numtrait;
00439 hmmec.logdet(&rank_mme);
00440 double ve = residual_var[0][0];
00441 double sst = yry*ve;
00442 ssm = blupsol.inner_product(rellrhs)*ve;
00443 sse = sst - ssm;
00444 rank_x = rank_mme - rank_z;
00445 dfe = static_cast<int>(numobs) - static_cast<int>(rank_x);
00446 }
|
|
|
Definition at line 967 of file model.cpp. References matvec::Model::data, matvec::Model::hashxact(), matvec::Data::num_rows(), matvec::Model::numterm, matvec::HashTable::resize(), matvec::HashTable::size(), matvec::Model::term, and matvec::Model::xact_htable. Referenced by matvec::Model::prepare_data().
00968 {
00969 // **********************************************************
00970 // get sequantial integer ids for interation term in model
00971 // ***********************************************************
00972 if (!data) return;
00973 unsigned nord, nrow_in_data = data->num_rows();
00974
00975 if (xact_htable) {
00976 delete [] xact_htable;
00977 xact_htable = 0;
00978 }
00979 if(numterm>0){
00980 xact_htable = new HashTable [numterm];
00981 }
00982 else {
00983 xact_htable = 0;
00984 }
00985 int t;
00986 for (t=0; t<numterm; t++) {
00987 if (term[t].order() > 1) {
00988 xact_htable[t].resize(nrow_in_data,term[t].order()*sizeof(unsigned));
00989 }
00990 }
00991 hashxact(record_legality);
00992 for (t=0;t<numterm;t++) {
00993 nord = term[t].order();
00994 if (nord>1) {
00995 term[t].numlevel = xact_htable[t].size();
00996 xact_htable[t].resize(term[t].nlevel(),nord*sizeof(unsigned));
00997 }
00998 }
00999 hashxact(record_legality);
01000 }
|
|
||||||||||||||||||||
|
Return BLUE solution for a fixed model.
Definition at line 238 of file model_blup.cpp. References matvec::Model::blup(), matvec::Model::fixed_model, and matvec::Model::type.
00240 {
00241 Vector<double> *retval = 0;
00242 if (type != fixed_model) throw exception("Model::blue(): it must be a fixed model");
00243 retval = blup(method,stopval,relax,mxiter);
00244 return retval;
00245 }
|
|
||||||||||||||||||||
|
Return the BLUP solution.
Definition at line 206 of file model_blup.cpp. References matvec::Model::bad_model, matvec::Model::blup_pccg(), matvec::Model::blupsol, matvec::Model::hmmec, matvec::SparseMatrix::nz(), matvec::Model::rellrhs, matvec::Model::setup_mme(), matvec::SparseMatrix::solve(), and matvec::Model::type. Referenced by matvec::Model::blue(), glim(), and matvec::Model::lsmeans().
00208 {
00209 Vector<double> *retval = 0;
00210 if (type == bad_model) throw exception("Model::blup(): bad model");
00211 if (method == "iod") {
00212 blup_pccg(stopval,mxiter);
00213 return &blupsol;
00214 }
00215 if (hmmec.nz() == 0) {
00216 try {
00217 setup_mme(&rellrhs);
00218 } catch (exception &er) {
00219 throw;
00220 }
00221 }
00222 if (hmmec.solve(blupsol,rellrhs,method,relax,stopval,mxiter) == 1) {
00223 retval = &blupsol;
00224 //hmmec.close();
00225 }
00226 // note that after BLUP, HMME remains intack
00227 // hmmec.resize(0,0);
00228 return retval;
00229 }
|
|
||||||||||||
|
Compute blup by iteration on data using preconditioned conjugate gradient see J.Anim.Sci (2001) 79:1166-1172 Definition at line 618 of file model_blup.cpp. References matvec::Model::bad_model, matvec::Vector< double >::begin(), matvec::Vector< T >::begin(), matvec::Model::blupsol, matvec::Vector< T >::clear(), matvec::Vector< double >::clear(), matvec::Model::compute_rhs_diag_mme(), matvec::Model::data_prepared, matvec::Model::fixed_model, matvec::ginverse2(), matvec::Model::hmmesize, matvec::doubleMatrix::identity(), matvec::Vector< T >::inner_product(), matvec::Model::inverse_residual_var(), matvec::Model::mixed_model, matvec::Model::mme_times(), matvec::Model::mme_times_res, matvec::Model::nt_vec, matvec::Matrix< double >::num_rows(), matvec::Model::numterm, matvec::Model::numtrait, matvec::preconditioning(), matvec::Model::prepare_data(), matvec::Model::rellrhs, matvec::Vector< T >::reserve(), matvec::Vector< double >::reserve(), matvec::Model::residual_var, matvec::Vector< double >::resize(), matvec::Model::term, and matvec::Model::type. Referenced by matvec::Model::blup().
00618 {
00619
00620 if (type == bad_model) throw exception("Model::blup_pccg: bad model");
00621 if (!data_prepared) {
00622 if (!prepare_data("iod")) {
00623 type = bad_model;
00624 throw exception("Model::blup_pccg: bad model");
00625 }
00626 }
00627 if (type == mixed_model) {
00628 if (residual_var.num_rows() != numtrait) {
00629 type = bad_model;
00630 throw exception("Model::blup_pccg: residual variance has not yet assigned");
00631 }
00632 }
00633 else if (type == fixed_model) {
00634 if (residual_var.num_rows() != numtrait) {
00635 residual_var.identity(numtrait,numtrait);
00636 }
00637 }
00638 else {
00639 throw exception("Model::setup_mme(): inappropriate model");
00640 }
00641
00642 inverse_residual_var();
00643 blupsol.resize(hmmesize,0.0);
00644 rellrhs.reserve(hmmesize);
00645 unsigned i,j,k;
00646 unsigned n;
00647
00648 unsigned nt = (numtrait+1)*numtrait/2;
00649 for (n=0,i=0; i<numterm; ++i) n += term[i].nlevel();
00650
00651 double **M;
00652 M = new(nothrow) double *[n];
00653 for (i=0; i<n; ++i) {
00654 M[i] = new(nothrow) double[nt+1]; //M[i][0] is reserved to store the size
00655 memset(&(M[i][1]),'\0',sizeof(double)*nt);
00656 }
00657
00658 k = 0;
00659 for (i=0; i<numterm; ++i) {
00660 nt = nt_vec[i];
00661 for (j=0; j<term[i].nlevel(); ++j) M[k++][0] = nt;
00662 }
00663 compute_rhs_diag_mme(M);
00664
00665 for (i=0; i<n; ++i) {
00666 // nt = (unsigned)M[i][0];
00667 // nt = nt*(nt+1)/2;
00668 // for (j=1; j<=nt; ++j) cout << " " << M[i][j];
00669 // cout << "\n";
00670 ginverse2(&(M[i][1]),(unsigned)M[i][0],1.0e-15);
00671 }
00672 mme_times_res.reserve(hmmesize);
00673 blupsol.resize(hmmesize);
00674 Vector<double> r;
00675 r = rellrhs;
00676
00677 Vector<double> d;
00678 d.reserve(hmmesize);
00679 ///////// d = M*r ///////////////
00680 preconditioning(d.begin(),(const double **)M,r.begin(),n);
00681
00682 double delta_new, delta_old, convcr = 1.0;
00683 unsigned iter = 0;
00684 double alpha,beta;
00685 delta_new = r.inner_product(d);
00686 double delta_0 = delta_new*tol*tol;
00687 bool done = false;
00688 while (!done){
00689 iter++;
00690 std::cout <<"iteration: " << iter <<std::endl;
00691 /////// q = Ad /////////
00692 mme_times(d);// the result is stored in mme_times_res
00693 alpha = delta_new/(d.inner_product(mme_times_res));
00694 blupsol += alpha*d;
00695 if(n%50){
00696 //////// r = r - alpha*q //////////
00697 r -= alpha*mme_times_res;
00698 }
00699 else{
00700 ///// r = b - Ax /////////
00701 mme_times(blupsol);
00702 r = rellrhs - mme_times_res;
00703 }
00704 ///////// q = M*r ///////////////
00705 preconditioning(mme_times_res.begin(),(const double **)M,r.begin(),n);
00706
00707 delta_old = delta_new;
00708 delta_new = r.inner_product(mme_times_res);
00709 beta = delta_new/delta_old;
00710 d = mme_times_res + beta*d;
00711
00712 if (delta_new > delta_0 && iter < maxiter){
00713 done = false;
00714 }
00715 else if (delta_new <= delta_0){
00716 ///// r = b - Ax /////////
00717 mme_times(blupsol);
00718 r = rellrhs - mme_times_res;
00719
00720 ///////// q = M*r ///////////////
00721 preconditioning(mme_times_res.begin(),(const double **)M,r.begin(),n);
00722
00723 delta_new = r.inner_product(mme_times_res);
00724 if(delta_new > delta_0 && iter < maxiter){
00725 done = false;
00726 }
00727 else {
00728 done = true;
00729 }
00730 }
00731 else {
00732 done = true;
00733 }
00734 }
00735 // std::cout <<"r = "<<r;
00736 mme_times_res.clear();
00737 r.clear();
00738 }
|
|
|
Definition at line 1151 of file glmm2.cpp. References matvec::Model::act_numtrait, corrmap, corrvar, CorrVar, ncorr, matvec::Model::numterm, matvec::Model::numtrait, matvec::Matrix< double >::resize(), and matvec::Model::term. Referenced by add_Ag_old().
01151 {
01152 CorrVar.resize(act_numtrait,act_numtrait);
01153 int ipos,jpos,i,j;
01154 doubleMatrix *var;
01155 Matrix<int> posmat(numterm,act_numtrait);
01156 for(i=0,ipos=0;i<numterm;i++) {
01157 if(term[i].classi() == 'P' && corrmap[i]==0) {
01158 corrmap[i]=1;
01159 ncorr++;
01160 }
01161 if(corrmap[i]) {
01162 for(int t=0;t<numtrait;t++) {
01163 posmat[i][t]=ipos++;
01164 }
01165 }
01166 }
01167 for(i=0;i<numterm;i++) {
01168 if(corrmap[i]) {
01169 var = term[i].prior->var_matrix();
01170 for(int t1=0;t1<numtrait;t1++) {
01171 for(int t2=0;t2<numtrait;t2++) {
01172 CorrVar[posmat[i][t1]][posmat[i][t2]]=(*var)[t1][t2];
01173 }
01174 }
01175 for(j=(i+1);j<numterm;j++) {
01176 if(corrmap[j]) {
01177 var = &corrvar[i][j];
01178 for(int t1=0;t1<numtrait;t1++) {
01179 for(int t2=0;t2<numtrait;t2++) {
01180 CorrVar[posmat[i][t1]][posmat[j][t2]]=(*var)[t1][t2];
01181 CorrVar[posmat[j][t2]][posmat[i][t1]]=(*var)[t1][t2];
01182 }
01183 }
01184 }
01185 }
01186 }
01187 }
01188 }
|
|
|
Definition at line 1750 of file glmm2.cpp. References matvec::SparseMatrix::a(), matvec::SparseMatrix::aav, matvec::Vector< T >::begin(), matvec::Vector< unsigned >::begin(), matvec::SparseMatrix::dense(), matvec::SparseMatrix::dim, matvec::SparseMatrix::factor(), matvec::SparseMatrix::factor_done, hInv, matvec::Model::hmmec, matvec::Model::hmmesize, matvec::SparseMatrix::hsize, matvec::hsort_ija(), matvec::SparseMatrix::ia(), matvec::SparseMatrix::iiv, matvec::SparseMatrix::insert(), matvec::SparseMatrix::inv_order, matvec::SparseMatrix::ja(), matvec::SparseMatrix::jjv, matvec::Model::max_nz, matvec::SparseMatrix::order, pt_mutex_destroy, pt_mutex_init, matvec::SparseMatrix::resize(), matvec::Vector< T >::resize(), matvec::SparseMatrix::solvrow(), and matvec::squeeze(). Referenced by residual().
01751 {
01752
01753 unsigned row=0;
01754 Vector<double> InvCol,wrkVec;
01755 int *lhs_ia,*lhs_ja,*lhs_sh,nz;
01756 double *lhs_a;
01757 unsigned hsize=hmmesize;
01758 InvCol.resize(hmmesize);
01759 wrkVec.resize(hmmesize);
01760 hInv.resize(hmmesize,max_nz);
01761 #ifdef DEBUG
01762 doubleMatrix LHS;
01763 LHS=hmmec.dense();
01764 std::cout << "\nbuild_hInv LHS\n"<<LHS <<"\n";
01765 #endif
01766 if (!hmmec.factor_done) hmmec.factor();
01767 unsigned int *order=hmmec.order.begin();
01768
01769 #ifdef DO_THREADS
01770 pthread_mutex_t mdone_mutex;
01771 pt_mutex_init(&mdone_mutex,"init mutex"); /* initialize the mutex */
01772
01773 thread_arg mult_arg;
01774 mult_arg.mdone_mutex=&mdone_mutex;
01775 mult_arg.A=&hmmec;
01776 mult_arg.Ainv=&hInv;
01777 mult_arg.hmmesize=hmmesize;
01778 mult_arg.rowdone=&row;
01779 mult_arg.order=hmmec.order.begin();
01780 #endif
01781 int *ia = hmmec.iiv-1;
01782 int *ja = hmmec.jjv-1;
01783 double *a = hmmec.aav-1;
01784 unsigned nrow = 0;
01785 unsigned oldrow = 0;
01786 nz=squeeze(hmmec.hsize,ia,ja,a);
01787 for(int i=1;i<=nz;i++) {
01788 ia[i]=hmmec.inv_order[ia[i]-1];
01789 ja[i]=hmmec.inv_order[ja[i]-1];
01790 if(ja[i] < ia[i]) {
01791 int itmp=ia[i];
01792 ia[i]=ja[i];
01793 ja[i]=itmp;
01794 }
01795 }
01796 hsort_ija(nz,ia,ja,a);
01797 for (int k=1; k<=nz ; k++) {
01798 if (ia[k] != oldrow) {
01799 nrow++;
01800 ia[nrow] = k;
01801 oldrow = ia[k];
01802 }
01803 }
01804 ia[hmmec.dim+1] = nz+1;
01805
01806 lhs_ia=hmmec.ia();
01807 lhs_ja=hmmec.ja();
01808 lhs_a=hmmec.a();
01809 #ifdef DO_THREADS
01810 mult_arg.ia=ia;
01811 mult_arg.ja=ja;
01812 mult_arg.lhs_ja=lhs_ja;
01813 mult_arg.lhs_ia=lhs_ia;
01814 int NTHREADS;
01815 char *MP_NUM;
01816 NTHREADS=1;
01817 MP_NUM=getenv("MP_NUM_THREADS");
01818 if(MP_NUM) NTHREADS=atol(MP_NUM);
01819 if(NTHREADS >1) {
01820 NTHREADS--;
01821 std::cout << "Nthreads "<< NTHREADS << "\n";
01822 //thread_code(&mult_arg[0]);
01823 // std::cout << &mult_arg << "\n";
01824 //TW pt_fork(NTHREADS,thread_code,(pt_addr_t) &mult_arg,NULL); /* run code in parallel */
01825
01826
01827 //TW pt_mutex_destroy(&mdone_mutex,"del mutex");
01828 // hInv.close();
01829 }
01830 else{
01831 #endif
01832 Vector<double> solwrk,temp_rhswrk;
01833 Vector<unsigned int> temp_col,temp_row;
01834 int m,jj;
01835 #pragma omp parallel private(jj,m,wrkVec,solwrk,temp_rhswrk,InvCol,temp_col,temp_row)
01836 {
01837 #pragma omp critical
01838 {
01839 wrkVec.resize(hmmesize);
01840 InvCol.resize(hmmesize);
01841 solwrk.resize(hmmesize);
01842 temp_rhswrk.resize(hmmesize);
01843 temp_row.resize(hmmesize);
01844 temp_col.resize(hmmesize);
01845 }
01846 #pragma omp for schedule(dynamic,50) nowait
01847 for(m=1;m<=hmmesize;m++) {
01848 wrkVec(order[m-1])=1;
01849 hmmec.solvrow(InvCol.begin(),wrkVec.begin(),m,solwrk,temp_rhswrk,temp_row.begin(),temp_col.begin());
01850 wrkVec(order[m-1])=0;
01851 {
01852 for(jj=lhs_ia[m];jj<lhs_ia[m+1];jj++) {
01853 // hInv.insert(order[m-1],order[lhs_ja[jj]-1],InvCol[order[lhs_ja[jj]-1]-1]);
01854 // hInv.insert(order[m-1],order[lhs_ja[jj]-1],solwrk[order[lhs_ja[jj]-1]-1]);
01855 hInv.insert(order[m-1],order[lhs_ja[jj]-1],InvCol[lhs_ja[jj]-1]);
01856 }
01857 }
01858 }
01859 }
01860 #ifdef DO_THREADS
01861 pt_mutex_destroy(&mdone_mutex,"del mutex");
01862 }
01863 #endif
01864 }
|
|
|
Definition at line 2158 of file model.cpp. References matvec::Model::base_effect, matvec::Model::categorical, matvec::check_ptr(), matvec::FieldStruct::classi(), matvec::Model::corr_link, matvec::Model::factor_struct, matvec::fieldindex(), matvec::FieldStruct::name(), matvec::nest_xaction(), matvec::Model::nt_vec, matvec::Model::numfact, matvec::Model::numterm, matvec::Model::numtrait, matvec::Model::pos_vec, matvec::Vector< Vector< int > >::reserve(), matvec::Vector< int >::reserve(), matvec::Vector< Vector< int > >::resize(), matvec::Vector< int >::resize(), matvec::TermList::resize(), matvec::sort(), matvec::split(), matvec::Model::term, matvec::Model::trait_map, matvec::Model::trait_struct, matvec::Model::traitname, matvec::Model::unknown_prior, and matvec::Model::var_link. Referenced by matvec::Model::equation().
02159 {
02160 int retval = 1;
02161 unsigned i,j,k;
02162
02163 std::string mspec = modelspec; // modelspec: y1 = A B + A (B), y2 = B + X*B
02164 nest_xaction(mspec); // now mspec: y1 = A B + A**B , y2 = B + X*B
02165
02166 std::string sep("=+,* ");
02167 std::string TMP = mspec;
02168 //////// TMP.unique("=+,* "); // TMP ->y1 y2 A B X
02169 std::vector<std::string> tmpvec;
02170 std::vector<std::string>::iterator it;
02171 std::string::size_type begidx;
02172 split(TMP,sep,&tmpvec);
02173
02174 sort(tmpvec.begin(),tmpvec.end());
02175 it = unique(tmpvec.begin(),tmpvec.end());
02176 tmpvec.erase(it,tmpvec.end());
02177 // for(it = tmpvec.begin(); it != tmpvec.end(); it++){
02178 // std::cout << *it << endl;
02179 // }
02180
02181 categorical = 0;
02182 numfact = tmpvec.size();
02183 if(factor_struct){
02184 delete [] factor_struct;
02185 factor_struct = 0;
02186 }
02187 if(numfact>0){
02188 factor_struct = new FieldStruct [numfact];
02189 }
02190 else {
02191 factor_struct = 0;
02192 }
02193 check_ptr(factor_struct);
02194 for (i=0; i<numfact; i++) factor_struct[i].name(tmpvec[i]);
02195
02196 TMP = mspec;
02197 numtrait = 0;
02198 sep = "=";
02199 begidx = TMP.find(sep,0);
02200 while(begidx != std::string::npos) {
02201 numtrait++;
02202 begidx = TMP.find(sep,begidx+1);
02203 }
02204 if (numtrait == 0) {
02205 std::cerr << "Model::build_model_term(): there is no trait in the model\n";
02206 return 0;
02207 }
02208 if(trait_struct){
02209 delete [] trait_struct;
02210 trait_struct = 0;
02211 }
02212 if(numtrait>0){
02213 trait_struct = new FieldStruct* [numtrait];
02214 }
02215 else {
02216 trait_struct = 0;
02217 }
02218 check_ptr(trait_struct);
02219 if(traitname){
02220 delete [] traitname;
02221 traitname = 0;
02222 }
02223 if(numtrait>0){
02224 traitname = new std::string [numtrait];
02225 }
02226 else {
02227 traitname = 0;
02228 }
02229 check_ptr(traitname);
02230 std::vector<std::string> model_trait;
02231 std::vector<std::string> trtnm;
02232 sep = ",";
02233 split(TMP,sep,&model_trait);
02234 if (numtrait != model_trait.size()) throw exception("Model::build_model_term(): model for each trait is not separated properly");
02235
02236 sep = " +";
02237 unsigned n = 0;
02238 for (i=0; i<numtrait; i++) {
02239 // trait name should be before "="
02240 begidx = model_trait[i].find("=");
02241 // bbbtrtnmbb is the trait name with possible surrounding blanks
02242 std::string bbbtrtnmbbb = model_trait[i].substr(0,begidx);
02243 // getting rid of the surrounding blanks
02244 split(bbbtrtnmbbb,sep,&trtnm);
02245 traitname[i] = trtnm[i];
02246 model_trait[i] = model_trait[i].substr(begidx+1);
02247 n += split(model_trait[i],sep); /// n += model_trait[i].ntoken("+ ");
02248 k = fieldindex(traitname[i],factor_struct,numfact);
02249 if (k >= 0 && k < numfact) {
02250 factor_struct[k].classi('T');
02251 trait_struct[i] = &(factor_struct[k]);
02252 }
02253 else {
02254 throw exception("Model::build_model_term(): you have probably found a bug!");
02255 }
02256 }
02257 TermList T(n,numtrait);
02258 tmpvec.clear();
02259 sep = " +";
02260 for (numterm=0,i=0; i<numtrait; i++) {
02261 tmpvec.clear();
02262 n = split(model_trait[i],sep,&tmpvec); // tmpvec = model_trait[i].split(n,"+ ");
02263 for (k=0; k<n; k++) {
02264 for (j=0; j<numterm; j++) {
02265 if (T[j].match(tmpvec[k],"* ",factor_struct)) break;
02266 }
02267 if (j==numterm) {
02268 T[numterm++].input(tmpvec[k].c_str(),factor_struct,numfact);
02269 }
02270 T[j].trait[i]= 1;
02271 }
02272 }
02273
02274 term.resize(numterm,numtrait);
02275 if (unknown_prior) {
02276 delete [] unknown_prior; // note that unknown_prior[i]
02277 unknown_prior = 0;
02278 }
02279 if(numterm>0){
02280 unknown_prior = new UnknownDist[numterm]; // needs to resize(numtrait)
02281 }
02282 else {
02283 unknown_prior = 0;
02284 }
02285 check_ptr(unknown_prior);
02286 nt_vec.resize(numterm,numtrait);
02287 base_effect.reserve(numterm);
02288 corr_link.resize(numterm,-1);
02289 pos_vec.resize(numterm);
02290 trait_map.reserve(numterm);
02291 var_link.reserve(numterm);
02292
02293
02294 for (i=0; i<numterm; i++) {
02295 base_effect[i]=i;
02296 var_link[i]=i;
02297 term[i] = T[i];
02298 term[i].prior = &(unknown_prior[i]);
02299 trait_map[i].resize(numtrait,i);
02300 }
02301 return retval;
02302 }
|
|
|
Definition at line 406 of file glmm1.cpp. References matvec::Matrix< double >::assign(), matvec::doubleMatrix::eigen(), matvec::Session::epsilon, matvec::doubleMatrix::ginv0(), kvec, link_function, matvec::log(), matvec::Vector< T >::max(), new_SinMat(), matvec::Model::nt_vec, matvec::Model::numterm, matvec::Model::numtrait, matvec::Model::residual_var, matvec::Vector< int >::resize(), matvec::Matrix< double >::resize(), matvec::SESSION, SinMat, matvec::Model::term, matvec::Matrix< double >::transpose(), Var, matvec::Model::var_link, and varinv. Referenced by SSQCP().
00407 {
00408 int i,k,t1,t2,raneff;
00409 #ifdef DO_CHOL
00410 doubleMatrix L;
00411 #endif
00412 #ifdef DO_LOG
00413 doubleMatrix P,L;
00414 Vector<double> D;
00415 #endif
00416 if(!SinMat)new_SinMat();
00417 if(!varinv)varinv = new doubleMatrix [numterm];
00418 if(!Var)Var = new doubleMatrix [numterm];
00419 int done_ped=0;
00420 doubleMatrix var,part;
00421 part.resize(numtrait,numtrait);
00422 kvec.resize(numterm);
00423 int kk=0;
00424 for(k=0,i=0,raneff=0;i<numterm;i++) {
00425 if((term[i].classi() == 'P'&& !done_ped) || term[i].classi() == 'R'){
00426 int nt=nt_vec[i];
00427 kvec[i]=k;
00428 if(var_link[i]!=i) {
00429 kvec[i]=kvec[var_link[i]];
00430 }
00431 else{
00432 raneff++;
00433 kk+=(nt*(nt+1))/2;
00434 }
00435 k=kvec[i];
00436 Var[i]=*term[var_link[i]].prior->var_matrix();
00437
00438
00439 #ifdef DO_CHOL
00440 L=Var[i];
00441 L.chol();
00442 L=L.transpose();
00443 #endif
00444 #ifdef DO_LOG
00445 P=Var[i];
00446 D=P.eigen();
00447 double dmax=D.max();
00448 if(dmax < 1) dmax=1;
00449 double dmin=10.*dmax*SESSION.epsilon;
00450 for(int ii=0;ii<nt;ii++) {
00451 if(D[ii] < dmin) D[ii]=.5*dmin;
00452 }
00453 D=log(D);
00454 #endif
00455 var=Var[i].ginv0();
00456 varinv[i]=var;
00457 #ifdef DO_PMAT
00458 var*=P_Mat[raneff]*P_Mat[raneff].transpose();
00459 #endif
00460 for(t1=0;t1<nt;t1++) for(t2=t1;t2<nt;t2++,k++) {
00461 part.resize(nt,nt);
00462 // std::cout << "PART\n";
00463 part[t1][t2]=1.;
00464 part[t2][t1]=1.;
00465 #ifdef DO_CHOL
00466 part[t1][t2]=0.;
00467 part[t2][t1]=0.;
00468 for(int j=0;j<nt;j++) {
00469 part[t2][j]+=L[t1][j];
00470 if(j != t2) part[j][t2]+=L[t1][j];
00471 }
00472 #endif
00473 #ifdef DO_LOG
00474 part[t1][t2]=0.;
00475 part[t2][t1]=0.;
00476 //std::cout << "DO LOG\n";
00477 double D_avg,D_delt;
00478 D_avg=.5*(D[t1]+D[t2]);
00479 D_delt=.5*(D[t1]-D[t2]);
00480 if(fabs(D_delt) > 1.e-8) {
00481
00482 // Sinh(x)=(e^x-e^{-x})/2
00483 part[t1][t2]=std::exp(D_avg)*sinh(D_delt)/(D_delt);
00484 part[t2][t1]=std::exp(D_avg)*sinh(D_delt)/(D_delt);
00485 }
00486 else{
00487 part[t1][t2]=std::exp(D_avg);
00488 part[t2][t1]=std::exp(D_avg);
00489 }
00490
00491 /**/
00492 if(D[t1]<-11 && D[t2]<-11) {
00493 part[t2][t1]=0;
00494 part[t1][t2]=0;
00495 }
00496 /**/
00497
00498
00499
00500 //std::cout << "part" << part;
00501 part=P*part*P.transpose();
00502 // if(std::exp(D[t1]) <= dmin || std::exp(D[t2]) <= dmin) part.assign(0.0);
00503 #endif
00504 SinMat[k]=var*part*var.transpose();
00505 //std::cout << "SinMat" << SinMat[k];
00506 }
00507
00508 }
00509
00510 }
00511 k=kk;
00512 if(!link_function) {
00513 var=residual_var.ginv0();
00514 var*=P_Mat[raneff]*P_Mat[raneff].transpose();
00515 for(t1=0;t1<numtrait;t1++) for(t2=t1;t2<numtrait;t2++,k++) {
00516 part.assign(0.0);
00517 part[t1][t2]=1.;
00518 part[t2][t1]=1.;
00519 SinMat[k]=var*part*var.transpose();
00520 }
00521 }
00522 }
|
|
|
|
|
|
Compute rhs and diagonals of mme Definition at line 742 of file model_blup.cpp. References matvec::Model::add_G_1_diag(), matvec::Vector< double >::begin(), matvec::Model::hmmesize, matvec::Model::input_pos_val_iod(), matvec::lidx(), matvec::Model::nt_vec, matvec::Model::numterm, matvec::Model::numtrait, matvec::Model::pos_val_vector, matvec::Model::rellrhs, matvec::Model::rve, and matvec::Model::term. Referenced by matvec::Model::blup_pccg().
00742 {
00743 unsigned i,j,k,kk,ii,jj,t1,t2,tt,rec;
00744 std::vector<pos_val_node>::iterator it;
00745 double **vep,val,vy,xval;
00746 memset(rellrhs.begin(),'\0',sizeof(double)*hmmesize);
00747 double *rhstmp = rellrhs.begin()-1;
00748 Matrix<double> ve(numtrait,numtrait);
00749
00750 for (it=pos_val_vector.begin(); it!=pos_val_vector.end(); it++) {
00751 if(!it->pos_term.empty()) {
00752 input_pos_val_iod(it);
00753 }
00754 else {
00755 continue;
00756 }
00757
00758 vep = rve[it->pos_term[numterm]].begin();
00759 val = it->xval_term[numterm]; // val = weight variable
00760 for (t1=0; t1<numtrait; t1++) {
00761 for (t2=0; t2<numtrait; t2++) ve[t1][t2] = val*vep[t1][t2];
00762 }
00763
00764 for (t1=0; t1<numtrait; t1++) {
00765 vy = 0.0;
00766 for(t2=0; t2<numtrait; t2++) vy += ve[t1][t2]*it->trait_vec[t2];
00767 for (i=0; i<numterm; i++) {
00768 ii = term[i].addr[t1];
00769 if (ii == 0) continue;
00770 k = (ii-1)/numtrait;
00771 xval = it->xval_term[i];
00772 for (t2=t1; t2<numtrait; t2++) {
00773 jj = term[i].addr[t2];
00774 if (jj) {
00775 val = xval* ve[t1][t2];
00776 M[k][lidx(t2,t1,nt_vec[i])+1] += val*xval;
00777 }
00778 }
00779 rhstmp[ii] += vy*xval;
00780 }
00781 }
00782 }
00783
00784 add_G_1_diag(M);
00785 }
|
|
|
Definition at line 204 of file model_gibbs.cpp. Referenced by matvec::Model::genotype_dist_gibbs(), matvec::Model::genotype_dist_peeling(), matvec::Model::genotypic_value_gibbs(), matvec::Model::genotypic_value_peeling(), matvec::Model::graph_log_likelihood_peeling(), and matvec::Model::log_likelihood_peeling().
00205 {
00206 unsigned i,ii,rec;
00207 double val;
00208 double *sol = bu.begin()-1;
00209 // for gcc-3.2 replaced
00210 // std::istrstream modelfile(modelstr, modelpcount);
00211 // with
00212 std::istringstream modelfile(modelstringstr);
00213 if (!modelfile) throw exception("Model::compute_xbzu(): cannot open file");
00214 for (rec=0; rec<numrec; rec++) {
00215 input_pos_val(modelfile);
00216 val = 0.0;
00217 for (i=0; i<numterm; i++) {
00218 ii = term[i].addr[0]; // only one trait allowed now
00219 if (ii) val += sol[pos_term[i]]*xval_term[i];
00220 }
00221 ii = rec_indid[rec];
00222 if (ii > 0) pop->popmember[ii-1]->xbzu(val);
00223 }
00224 //BRS modelfile.close();
00225 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 276 of file model_hypo.cpp. References matvec::warning().
00278 {
00279 ////////////////////////////////////////////////////
00280 // H0: K'*b = M
00281 // trailing zeros in K' can be omitted
00282 // REMINDER: blupsol & rellrhs remain intact
00283 //
00284 // if M=NULL, then it acts as a vector of zeros.
00285 // meanp must be declared as double meanp[nr][4]
00286 /////////////////////////////////////////////////////
00287 if (!kpme) {
00288 warning("Model::contrast(kp): kp is null, thus it is ignored");
00289 return 0.0;
00290 }
00291 if (!get_blupsol()) throw exception("Model::contrast(): bad model");
00292 if (nc > hmmesize) throw exception("Model::contrast(): size not conformable");
00293 int est=1;
00294 doubleMatrix kvk(nr,nr);
00295 doubleMatrix kvkori(nr,nr);
00296 Vector<double> kpdiff(nr);
00297 Vector<double> kpb_m(nr);
00298 Vector<double> xy(hmmesize);
00299 Vector<double> tmpsol(hmmesize);
00300 double *sol, kpd;
00301 const double *kp;
00302 unsigned i,j,k;
00303 for (i=0; i<nr; i++) {
00304 if (nc < hmmesize) {
00305 memset(xy.begin(),'\0',sizeof(double)*hmmesize);
00306 memcpy(xy.begin(),kpme[i],sizeof(double)*nc);
00307 if (!hmmec.solve(tmpsol,xy,"ysmp")) return 0.0;
00308 }
00309 else {
00310 if (!hmmec.solve(tmpsol.begin(),kpme[i],"ysmp")) return 0.0;
00311 }
00312 hmmec.mv(tmpsol,xy);
00313 kp = kpme[i]; sol = xy.begin();
00314 for (kpd=0.0,j=0; j<nc; j++) {
00315 kpd = std::max(kpd,fabs(*kp - *sol));
00316 kp++;
00317 sol++;
00318 }
00319 for (j=nc; j<hmmesize; j++) {
00320 kpd = std::max(kpd,fabs(*sol));
00321 sol++;
00322 }
00323 kpdiff[i] = kpd;
00324 for (k=0; k<nr; k++) {
00325 kp = kpme[k]; sol = tmpsol.begin();
00326 for (kpd=0.0,j=0; j<nc; j++) kpd += *kp++ * *sol++;
00327 kvk[k][i] = kpd;
00328 }
00329 }
00330 kvkori = kvk;
00331 kvk.ginv1(&k);
00332 if (k < nr) throw exception("Model::contrast(K'): K'V{-1}K is singular. Possible reason:\n"
00333 " (1) K is non-estimable, and (2) K isn't of full column rank");
00334 // K'b-M
00335 for (i=0; i<nr; i++) kpb_m[i] = blupsol.inner_product(kpme[i]);
00336 if (M) for (i=0; i<nr; i++) kpb_m[i] -= M[i];
00337 double quq = kvk.quadratic(kpb_m,kpb_m); // (K'b-m)'(K'(X'V-1X)-K)-1(K'b-m)
00338 double ssm,sse;
00339 unsigned rank_x;
00340 int dfe = 0;
00341 anova(ssm,sse,rank_x,dfe);
00342 double res_var_est = sse/dfe;
00343 int degen = 0;
00344 int errcode = 0;
00345 double f_stat=0.0, prob=0.0;
00346 double sigma_e = residual_var[0][0];
00347 int chisq = 0;
00348 if (chisq) {
00349 f_stat = quq;
00350 prob = ChiSquare_cdf(f_stat,static_cast<double>(nr),0.0,errcode);
00351 }
00352 else {
00353 if (dfe <= 0 || sse <= SESSION.epsilon) {
00354 warning("Model::contrast(): either dfe or sse is zero");
00355 degen = 1;
00356 }
00357 else {
00358 f_stat = (quq*sigma_e/res_var_est)/nr;
00359 prob = F_cdf(f_stat,static_cast<double>(nr),static_cast<double>(dfe),0.0,errcode);
00360 }
00361 // in univariate cases only: the estimated variance of estimator K'B
00362 if (!degen) kvkori *= res_var_est/sigma_e;
00363 }
00364 std::string raw_data_code;
00365 double p_value,var,tcal=0.0;
00366 int W = SESSION.output_precision;
00367 std::cout.precision(W);
00368 if (prt_flag) {
00369 std::cout << "\n RESULTS FROM CONTRAST(S)\n";
00370 std::cout << " ----------------------------------------------------------\n";
00371 std::cout << " Contrast MME_addr K_coef Raw_data_code\n";
00372 std::cout << " ---------------------------------------------\n";
00373 }
00374 for (i=0; i<nr; i++) {
00375 if (prt_flag) {
00376 kp = kpme[i];
00377 for (j=0; j<nc; j++) {
00378 if (kp[j] == 0.0) continue;
00379 trait_effect_level(j,raw_data_code);
00380 std::cout << std::setw(4) << i+1 <<std::setw(11) << j+1 << std::setw(13) << kp[j]
00381 << " " << raw_data_code << "\n";
00382 }
00383 }
00384 kpd = kpdiff[i];
00385 if (result) {
00386 result[i][0] = kpb_m[i];
00387 result[i][3] = kpd;
00388 }
00389 var = kvkori[i][i];
00390 if (var <= 0.0) throw exception(" Model::contrast(): you have probably found a bug!");
00391 var = std::sqrt(var);
00392 tcal = fabs(kpb_m[i]/var);
00393 p_value = 2.0*(1.0-t_cdf(tcal,static_cast<double>(dfe),0.0));
00394 if (result) {
00395 result[i][1] = var;
00396 result[i][2] = p_value;
00397 }
00398 if (!prt_flag) continue;
00399 if (kpd > 1.0e-5) {
00400 est = 0;
00401 std::cout << " **** NON-ESTIMABLE ****\n\n";
00402 }
00403 else {
00404 if (kpd > 1.0e-10) {
00405 std::cout << " ***WARNING***: it may or mayn't be estimable\n";
00406 std::cout << " max(k'*ginv(mme)*mme - k') = " << kpd << "\n";
00407 }
00408 std::cout << " estimated value (K'b-M) = " << kpb_m[i]
00409 << " +- " << var << "\n";
00410 std::cout << " Prob(|t| > " << tcal << ") = "
00411 << p_value << " (p_value)\n";
00412 if (i != nr-1) std::cout << "\n";
00413 }
00414 }
00415 if (prt_flag) {
00416 std::cout << " ----------------------------------------------------------\n";
00417 if (nr > 1 && est) {
00418 std::cout << " joint hypothesis test H: K'b = M\n";
00419 std::cout << " Prob(F > " << f_stat << ") = "
00420 << 1.0-prob << " (p_value)\n";
00421 }
00422 }
00423 return 1.0-prob; // p_value
00424 }
|
|
||||||||||||||||||||||||||||
|
Definition at line 250 of file glmm2.cpp. References matvec::SparseMatrix::a(), ainv(), matvec::Matrix< double >::assign(), matvec::Vector< T >::begin(), matvec::Model::blupsol, matvec::ChiSquare_cdf(), EPSILON, matvec::F_cdf(), matvec::Model::get_blupsol(), matvec::doubleMatrix::ginv0(), matvec::doubleMatrix::ginv1(), glim(), matvec::Observe::H, hainv, matvec::Model::hmmec, matvec::Model::hmmesize, matvec::SparseMatrix::ia(), Info, matvec::Vector< double >::inner_product(), matvec::doubleMatrix::inv(), matvec::SparseMatrix::ja(), like_val, link_function, LR, matvec::SparseMatrix::mv(), matvec::Model::nt_vec, matvec::SparseMatrix::num_rows(), matvec::Matrix< double >::num_rows(), matvec::Model::numrec, matvec::Model::numterm, matvec::Model::numtrait, observation, matvec::Session::output_precision, Print_Level, Print_Summary, matvec::doubleMatrix::quadratic(), res_nvc, matvec::Observe::Resid_Sin_Mat, residual(), matvec::Model::residual_var, matvec::Observe::rve, matvec::SESSION, SinMat, matvec::SparseMatrix::solve(), matvec::Matrix< double >::submat(), matvec::t_cdf(), matvec::Model::term, TotalNvc(), matvec::Model::trait_effect_level(), matvec::Matrix< double >::transpose(), matvec::warning(), and matvec::Observe::weight.
00252 {
00253 ////////////////////////////////////////////////////
00254 // H0: K'*b = M
00255 // trailing zeros in K' can be omitted
00256 // REMINDER: blupsol & rellrhs remain intact
00257 //
00258 // meanp must be declared as double meanp[nr][4]
00259 /////////////////////////////////////////////////////
00260
00261
00262 if(!link_function){
00263 return(Model::contrast(kpme,nr,nc,M, prt_flag,result));
00264 }
00265 int do_varcomp=Info.num_rows();
00266 int Print_Flag=1;
00267 if(Print_Level < 0 || !prt_flag) Print_Flag=0;
00268
00269 if (!kpme) {
00270 warning("GLMM::contrast(kp): kp is null, thus it is ignored");
00271 return 0.0;
00272 }
00273 if (!get_blupsol()) throw exception("GLMM::contrast(): bad model");
00274
00275 if (nc > hmmesize) throw exception("GLMM::contrast(): size not conformable");
00276 int est=1;
00277 doubleMatrix kvk(nr,nr);
00278 doubleMatrix kvkori(nr,nr);
00279 doubleMatrix kv(nr,hmmesize);
00280 doubleMatrix kGk(nr,nr);
00281 doubleMatrix kRk(nr,nr);
00282 Vector<double> kpdiff(nr);
00283 Vector<double> kpb_m(nr);
00284 Vector<double> xy(hmmesize);
00285 Vector<double> tmpsol(hmmesize);
00286 double *sol, kpd;
00287 double log_LR;
00288 log_LR=2.*like_val;
00289 const double *kp;
00290 unsigned i,j,k;
00291 int t1,t2,ii,*ainv_ia,*ainv_ja,ainv_nrow;
00292 int ipos,jpos,je,jj,jaddr;
00293 double *ainv_a;
00294 if(LR==2) glim(10,(double **)kpme,nr,nc,(double *)M);//Lagrange Multiplier uses KVK under Ho
00295
00296 for (i=0; i<nr; i++) {
00297 if (nc < hmmesize) {
00298 memset(xy.begin(),'\0',sizeof(double)*hmmesize);
00299 memcpy(xy.begin(),kpme[i],sizeof(double)*nc);
00300 hmmec.solve(tmpsol,xy,"ysmp");
00301 }
00302 else {
00303 hmmec.solve(tmpsol.begin(),kpme[i],"ysmp");
00304 }
00305 hmmec.mv(tmpsol,xy);
00306 memcpy(kv[i],tmpsol.begin(),sizeof(double)*hmmesize);
00307 kp = kpme[i];
00308 sol = xy.begin();
00309 for (kpd=0.0,j=0; j<nc; j++) {
00310 kpd = std::max(kpd,fabs(kpme[i][j] - xy[j]));
00311 kp++;
00312 sol++;
00313 }
00314 for (j=nc; j<hmmesize; j++) {
00315 kpd = std::max(kpd,fabs(xy[j]));
00316 sol++;
00317 }
00318 kpdiff[i] = kpd;
00319 for (k=0; k<nr; k++) {
00320 kp = kpme[k]; sol = tmpsol.begin();
00321 for (kpd=0.0,j=0; j<nc; j++) kpd += kpme[k][j] * tmpsol[j];
00322 kvk[k][i] = kpd;
00323 }
00324 }
00325 for(k=1;k<nr;k++) for(i=0;i<k;i++) kvk[k][i]=kvk[i][k];
00326 kvkori=kvk;
00327 kvk = kvkori;
00328 kvk.ginv1(&k);
00329 if (k < nr) {
00330 warning("GLMM::contrast(K'): K'V{-1}K is singular. Possible reason:\n"
00331 " (1) K is non-estimable, and (2) K isn't of full column rank");
00332 //return 0.0;
00333 }
00334 // Get Lambda's for vc calculations
00335 int nvc=TotalNvc();
00336 unsigned startaddr,nlevels,iaddr,t;
00337 Vector<double> Lambda(nvc);
00338 double addj=0.;
00339 doubleMatrix *smat;
00340 int k_vc,i_vc;
00341 doubleMatrix wmat(numtrait,numtrait);
00342 doubleMatrix rmat(numtrait,numtrait);
00343 for(k_vc=0,i_vc=0;do_varcomp && i_vc<numterm;i_vc++) {
00344 if(term[i_vc].classi() == 'P' || term[i_vc].classi() == 'R'){
00345 startaddr=term[i_vc].startaddr()+1;
00346 int nt=nt_vec[i_vc];
00347 for(t1=0;t1<nt;t1++) for(t2=t1;t2<nt;t2++,k_vc++) {
00348
00349 rmat=*(term[i_vc].prior->var_matrix());
00350 rmat.ginv1();
00351 wmat=SinMat[k_vc];
00352 kGk.assign(0.0);
00353 kRk.assign(0.0);
00354 if(term[i_vc].classi() == 'R') {
00355 nlevels=term[i_vc].nlevel();
00356 for(ii=1;ii<=nlevels;ii++) {
00357 iaddr=startaddr+(ii-1)*numtrait;
00358 kGk+=kv.submat(0,iaddr-1,nr,nt)*wmat
00359 *(kv.submat(0,iaddr-1,nr,nt)).transpose();
00360 if(t1==0 && t2==0) kRk+=kv.submat(0,iaddr-1,nr,nt)*rmat
00361 *(kv.submat(0,iaddr-1,nr,nt)).transpose();
00362 }
00363 }
00364 if(term[i_vc].classi() == 'P') {
00365 ainv();
00366 ainv_ia=hainv.ia();
00367 ainv_ja=hainv.ja();
00368 ainv_a=hainv.a();
00369 ainv_nrow=hainv.num_rows();
00370 for(ii=1;ii<=ainv_nrow;ii++) {
00371 ipos=ii;
00372 iaddr=startaddr+(ipos-1)*nt;
00373 je=ainv_ia[ii+1];
00374 for(jj=ainv_ia[ii];jj<je;jj++) {
00375 jpos=ainv_ja[jj];
00376 double offdiag;
00377 offdiag=2.;
00378 if(ipos==jpos) offdiag=1.;
00379 jaddr=startaddr+(jpos-1)*nt;
00380 kGk+=ainv_a[jj]*offdiag*
00381 kv.submat(0,iaddr-1,nr,nt)*wmat
00382 *(kv.submat(0,iaddr-1,nr,nt)).transpose();
00383 if(t1==0 && t2==0) kRk+=kv.submat(0,iaddr-1,nr,nt)*rmat
00384 *(kv.submat(0,iaddr-1,nr,nt)).transpose();
00385 }
00386 }
00387 }
00388 Lambda[k_vc]=(kvk*kGk).trace();
00389 if(t1 == 0 && t2 == 0) addj+=(kvk*kRk).trace();
00390 }
00391 }
00392 }
00393
00394 // ****
00395 // Build Residual SSs and CP
00396 // ****
00397 unsigned kvc;
00398
00399 if(link_function){
00400
00401 for(kvc=0;do_varcomp && kvc<res_nvc;kvc++,k_vc++){
00402 wmat.assign(0.0);
00403 rmat.assign(0.0);
00404 for(ii=0;ii<numrec;ii++) {
00405 wmat+=observation[ii].H.transpose()*observation[ii].Resid_Sin_Mat[kvc]*
00406 observation[ii].H*observation[ii].weight;
00407 rmat+=observation[ii].rve*observation[ii].weight;
00408 //
00409 // At this point I am not sure if it is better to adjust
00410 // for prediction errors or leave it alone.
00411 // For "most" of the cases I typically run accross it is
00412 // probably not much of an issue.
00413 //
00414
00415 }
00416 Lambda[k_vc]=(static_cast<double>(nr) - addj)*(wmat*rmat.inv()).trace();
00417
00418 }
00419
00420
00421 }
00422 else{
00423 rmat.assign(0.0);
00424 for(ii=0;ii<numrec;ii++) {
00425 rmat+=observation[ii].rve*observation[ii].weight;
00426 }
00427 for(t1=0;t1<numtrait;t1++) for(t2=t1;t2<numtrait;t2++,k_vc++) {
00428 wmat.assign(0.0);
00429 for(ii=0;ii<numrec;ii++) {
00430 wmat+=SinMat[k_vc]*observation[ii].weight;
00431
00432 }
00433 Lambda[k_vc]=(static_cast<double>(nr) - addj)*(wmat*rmat.inv()).trace();
00434 }
00435 }
00436
00437 if(LR==2) glim();
00438 double dfe;
00439
00440
00441 // K'b-M
00442 for (i=0; i<nr; i++) kpb_m[i] = blupsol.inner_product(kpme[i]);
00443 if (M) for (i=0; i<nr; i++) kpb_m[i] -= M[i];
00444 double quq;
00445 if(LR!=1) {
00446 quq = kvk.quadratic(kpb_m,kpb_m); // (K'b-m)'(K'(X'V-1X)-K)-1(K'b-m)
00447 }
00448 else{
00449 // tmpsol=blupsol-((kvk*kpb_m)*kv); //tmpsol now contains the restricted solutions
00450 {
00451 //double *blupve;
00452 //blupve=blupsol.ve;
00453 //blupsol.ve=tmpsol.ve;
00454 glim(10,(double **)kpme,nr,nc,(double *)M);
00455 residual();
00456 quq= -like_val;
00457 //blupsol.ve=blupve;
00458 glim();
00459 residual();
00460 quq+=like_val;
00461 quq*=2.;
00462 }
00463 }
00464 double ssm,sse;
00465 unsigned rank_x;
00466 int degen = 0;
00467 int errcode = 0;
00468 double f_stat=0.0, prob=0.0,vart=0.0;
00469 double sigma_e = residual_var[0][0];
00470 int not_chisq = Info.num_rows();
00471 if(Asym) not_chisq=0;
00472
00473 if(not_chisq) vart=Info.ginv0().quadratic(Lambda,Lambda);
00474 if (!not_chisq || vart <= EPSILON ) {
00475 f_stat = quq;
00476
00477 prob = ChiSquare_cdf(f_stat,static_cast<double>(nr),0.0,errcode);
00478 f_stat /= static_cast<double>(nr);
00479 dfe=2000.;
00480 }
00481 else {
00482
00483 dfe=(2.0*nr*nr)/vart;
00484 if(dfe < 1.) dfe=1.;
00485 if(dfe > 2000.) {
00486 dfe=2000.;
00487 }
00488 f_stat = (quq)/nr;
00489 // if(f_stat > 1000.) f_stat=1000.;
00490 prob = F_cdf(f_stat,static_cast<double>(nr),static_cast<double>(dfe),0.0,errcode);
00491
00492 // in univariate cases only: the estimated variance of estimator K'B
00493 }
00494 std::string raw_data_code;
00495 double p_value,var,tcal=0.0;
00496 int W = SESSION.output_precision;
00497 std::cout.precision(W);
00498 if(Print_Level < 0) Print_Flag=0;
00499 if (Print_Flag) {
00500 std::cout << "\n RESULTS FROM CONTRAST(S)\n";
00501 std::cout << " ----------------------------------------------------------\n";
00502 std::cout << " Contrast MME_addr K_coef Raw_data_code\n";
00503 std::cout << " ---------------------------------------------\n";
00504 }
00505 for (i=0; i<nr; i++) {
00506 if (Print_Flag) {
00507 kp = kpme[i];
00508 for (j=0; j<nc; j++) {
00509 if (kp[j] == 0.0) continue;
00510 trait_effect_level(j,raw_data_code);
00511 std::cout << std::setw(4) << i+1 <<std::setw(11) << j+1 << std::setw(13) << kp[j]
00512 << " " << raw_data_code << "\n";
00513 }
00514 }
00515 kpd = kpdiff[i];
00516 if (result) {
00517 result[i][0] = kpb_m[i];
00518 result[i][3] = kpd;
00519 }
00520 var = kvkori[i][i];
00521 if (var <= 0.0) throw exception("GLMM::contrast(): you have probably found a bug!");
00522 var = std::sqrt(var);
00523 // std::cout << var << " \n";
00524 // std::cout << kpb_m[i] << " xxxxxx \n";
00525 tcal = fabs(kpb_m[i]/var);
00526 // std::cout << tcal << " "<< dfe << " \n";
00527 p_value = 2.0*(1.0-t_cdf(tcal,static_cast<double>(dfe),0.0));
00528
00529 if (result) {
00530 result[i][1] = var;
00531 result[i][2] = p_value;
00532 }
00533 if (!Print_Flag) continue;
00534 if (kpd > 1.0e-5) {
00535 est = 0;
00536 std::cout << " **** NON-ESTIMABLE ****\n\n";
00537 }
00538 else {
00539 if (kpd > 1.0e-10) {
00540 std::cout << " ***WARNING***: it may or mayn't be estimable\n";
00541 std::cout << " max(k'*ginv(mme)*mme - k') = " << kpd << "\n";
00542 }
00543 std::cout << " estimated value (K'b-M) = " << kpb_m[i]
00544 << " +- " << var << "\n";
00545 std::cout << " Prob(|t| > " << tcal << ") = "
00546 << p_value << " (p_value)\n";
00547 if(nr == 1 && not_chisq)
00548 std::cout << " " << dfe << " (Error degrees of freedom)\n";
00549 if (i != nr-1) std::cout << "\n";
00550 }
00551 }
00552 if (Print_Flag || Print_Summary) {
00553 std::cout << " ----------------------------------------------------------\n";
00554 if ((nr > 1 || Print_Summary) && est) {
00555 std::cout << " joint hypothesis test H: K'b = M\n";
00556 std::cout << " Prob(F > " << f_stat << ") = "
00557 << 1.0-prob << " (p_value)\n";
00558 if(not_chisq) std::cout << " " << dfe << " (Error degrees of freedom)\n";
00559 }
00560
00561 }
00562 if(Return_Stat) return f_stat;
00563 return 1.0-prob; // p_value
00564 }
|
|
||||||||||||||||
|
Reimplemented from matvec::Model. Definition at line 623 of file glmm2.cpp. References matvec::Matrix< T >::begin(), matvec::Vector< T >::begin(), matvec::Matrix< double >::begin(), contrast(), matvec::Model::factor_struct, matvec::Model::get_blupsol(), matvec::Model::hmmesize, matvec::TermList::index(), matvec::Model::nt_vec, matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), matvec::Vector< T >::size(), matvec::Model::term, and matvec::warning().
00624 {
00625 if (!(Kp.begin())) {
00626 warning("GLMM::contrast(kp): kp is null, thus it is ignored");
00627 return 0.0;
00628 }
00629
00630 double retval = 0.0;
00631 if (!get_blupsol()) throw exception("GLMM::contrastranspose(): bad model");
00632
00633 int k = term.index(termname,factor_struct);
00634 if (k < 0) {
00635 warning("GLMM::contrast(): no such term in the model");
00636 return retval;
00637 }
00638
00639 unsigned nr = Kp.num_rows();
00640 unsigned nc = Kp.num_cols();
00641 if (M.size() != nr) {
00642 warning("GLMM::contrast(): unconformable size");
00643 return retval;
00644 }
00645 unsigned i,startaddr;
00646 int nt=nt_vec[k];
00647 if (nc <= term[k].nlevel()*nt) {
00648 startaddr = term[k].startaddr();
00649 Vector<double> fullsize_M(hmmesize);
00650 Matrix<double>fullsize_Kp(nr,hmmesize);
00651 for (i=0; i<nr; i++) {
00652 memcpy(&(fullsize_Kp[i][startaddr]),Kp[i], sizeof(double)*nc);
00653 }
00654 memcpy(&(fullsize_M[startaddr]),M.begin(), sizeof(double)*nc);
00655 retval=contrast((const double **)fullsize_Kp.begin(),nr,hmmesize,fullsize_M.begin());
00656 }
00657 else {
00658 warning("GLMM::contrast(%s,Kp): too many columns in Kp",termname.c_str());
00659 }
00660 return retval;
00661 }
|
|
|
Reimplemented from matvec::Model. Definition at line 593 of file glmm2.cpp. References matvec::Matrix< double >::begin(), contrast(), matvec::Matrix< double >::num_cols(), and matvec::Matrix< double >::num_rows().
00594 {
00595 ///////////////////////////////////////////
00596 // trailing zeros in K' can be omitted
00597 ///////////////////////////////////////////
00598
00599 return contrast((const double **)Kp.begin(), Kp.num_rows(), Kp.num_cols());
00600 }
|
|
||||||||||||
|
Reimplemented from matvec::Model. Definition at line 602 of file glmm2.cpp. References matvec::Matrix< T >::begin(), matvec::Vector< T >::begin(), matvec::Model::hmmesize, matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), and matvec::Vector< T >::size().
00603 {
00604 ///////////////////////////////////////////
00605 // trailing zeros in K' can be omitted
00606 ///////////////////////////////////////////
00607
00608
00609 double retval = 0.0;
00610 int nr = Kp.num_rows();
00611 unsigned nc = Kp.num_cols();
00612 if (M.size() != nr ) throw exception(" GLMM::contrast: size is incomformable");
00613 Vector<double> fullsize_M(hmmesize);
00614 Matrix<double>fullsize_Kp(nr,hmmesize);
00615 for (unsigned i=0; i<nr; i++) {
00616 memcpy(fullsize_Kp[i],Kp[i], sizeof(double)*nc);
00617 }
00618 memcpy(fullsize_M.begin(),M.begin(), sizeof(double)*nc);
00619 retval=GLMM::contrast((const double **)fullsize_Kp.begin(),nr,hmmesize,fullsize_M.begin());
00620 return retval;
00621 }
|
|
||||||||||||
|
Reimplemented from matvec::Model. Definition at line 573 of file glmm2.cpp. References matvec::Vector< T >::begin(), and matvec::Vector< T >::size(). Referenced by contrast().
00574 {
00575 ///////////////////////////////////////////
00576 // trailing zeros in K' can be omitted
00577 ///////////////////////////////////////////
00578
00579 double retval, **kpme = new double *[1];
00580 kpme[0] = Kp.begin();
00581 unsigned nc = Kp.size();
00582 double* M = (double *)NULL;
00583 if (m != 0.0) {
00584 M = new double [1];
00585 M[0] = m;
00586 }
00587 retval = contrast((const double **)kpme,1,nc,M);
00588 delete [] kpme;
00589 if (M) delete [] M;
00590 return retval;
00591 }
|
|
|
Definition at line 83 of file model.cpp. References matvec::Model::blupsol, matvec::check_ptr(), matvec::UnknownDist::copyfrom(), matvec::Model::data, matvec::Model::data_prepared, matvec::Model::dfreml_called, matvec::Model::factor_struct, matvec::fieldindex(), matvec::Data::hashtable, matvec::Model::hmmec, matvec::Model::hmmesize, matvec::Model::idlist, matvec::Model::initialize(), matvec::Model::kvec, matvec::Model::lng0vec, matvec::Model::lnr0vec, matvec::Model::max_nz, matvec::Model::maxorder, matvec::Model::modelfname, matvec::Model::modelstring, matvec::GeneticDist::name(), matvec::Model::nfunk_in_dfreml, matvec::Model::non_zero, matvec::Model::npattern, matvec::Model::ntermGdist, matvec::Model::num_ped, matvec::Model::numcol, matvec::Model::numfact, matvec::Model::numgroup, matvec::Model::numobs, matvec::Model::numrec, matvec::Model::numterm, matvec::Model::numtrait, matvec::Model::pattern, matvec::Model::pattern_tree, matvec::Model::pop, matvec::Model::pop_created, matvec::Model::popsize, matvec::Model::pos_term, matvec::ModelTerm::prior, matvec::Model::rawpos, matvec::Model::rec_indid, matvec::Model::rellrhs, matvec::Model::reml_value, matvec::Matrix< double >::reserve(), matvec::Model::residual_var, matvec::SparseMatrix::resize(), matvec::Model::rve, matvec::Model::term, matvec::TermList::termlist, matvec::Model::trait_struct, matvec::Model::trait_vec, matvec::Model::traitname, matvec::Model::type, matvec::Model::unknown_prior, matvec::Model::weightinverse, matvec::Model::weightname, matvec::Model::xact_htable, matvec::Model::xval_term, and matvec::Model::yry. Referenced by matvec::Model::Model().
00084 {
00085 initialize();
00086 num_ped = M.num_ped;
00087 yry = M.yry;
00088 data = M.data;
00089 term = M.term;
00090 max_nz = M.max_nz;
00091 numrec = M.numrec;
00092 numobs = M.numobs;
00093 numcol = M.numcol;
00094 max_nz = M.max_nz;
00095 popsize = M.popsize;
00096 numterm = M.numterm;
00097 numfact = M.numfact;
00098 lng0vec = M.lng0vec;
00099 lnr0vec = M.lnr0vec;
00100 blupsol = M.blupsol;
00101 rellrhs = M.rellrhs;
00102 popsize = M.popsize;
00103 hmmesize = M.hmmesize;
00104 numtrait = M.numtrait;
00105 npattern = M.npattern;
00106 hmmesize = M.hmmesize;
00107 numgroup = M.numgroup;
00108 maxorder = M.maxorder;
00109 non_zero = M.non_zero;
00110 rec_indid = M.rec_indid;
00111 type = M.type;
00112 modelfname = M.modelfname;
00113 weightname = M.weightname;
00114 ntermGdist = M.ntermGdist;
00115 modelstring = M.modelstring;
00116 pattern_tree = M.pattern_tree;
00117 residual_var = M.residual_var;
00118 data_prepared = M.data_prepared;
00119 weightinverse = M.weightinverse;
00120 nfunk_in_dfreml= M.nfunk_in_dfreml;
00121 reml_value = M.reml_value;
00122 dfreml_called = M.dfreml_called;
00123 hmmec.resize(hmmesize,max_nz);
00124
00125 int i,k;
00126 // pos_term, xval_term, trait_vec are working spaces
00127 if (pos_term) {delete [] pos_term; pos_term = 0;}
00128 if (M.pos_term) pos_term = new unsigned [numterm+1];
00129 if (xval_term) {delete [] xval_term; xval_term = 0;}
00130 if (M.xval_term) xval_term = new double [numterm+1];
00131 if (trait_vec) {delete [] trait_vec; trait_vec= 0;}
00132 if (M.trait_vec) trait_vec = new double [numterm+1];
00133
00134 if (unknown_prior) {
00135 delete [] unknown_prior;unknown_prior = 0;
00136 }
00137 if (M.unknown_prior) {
00138 if(numterm>0){
00139 unknown_prior = new UnknownDist[numterm];
00140 }
00141 else {
00142 unknown_prior = 0;
00143 }
00144 for (i=0; i<numterm; i++) {
00145 unknown_prior[i].copyfrom(M.unknown_prior[i]);
00146 if (strcmp(M.term.termlist[i].prior->name(),"UnknownDist") == 0) {
00147 term.termlist[i].prior = &(unknown_prior[i]);
00148 }
00149 else {
00150 term.termlist[i].prior = M.term.termlist[i].prior;
00151 }
00152 }
00153 }
00154 if (traitname) {delete [] traitname; traitname = 0;}
00155 if (M.traitname) {
00156 if(numtrait>0){
00157 traitname = new std::string [numtrait];
00158 }
00159 else {
00160 traitname = 0;
00161 }
00162 for (i=0; i<numtrait; i++) traitname[i] = M.traitname[i];
00163 }
00164 if (kvec) {delete [] kvec; kvec = 0;}
00165 if (M.kvec) {
00166 if(numtrait>0){
00167 kvec = new unsigned [npattern];
00168 }
00169 else {
00170 kvec = 0;
00171 }
00172 for (i=0; i<npattern; i++) kvec[i] = M.kvec[i];
00173 }
00174 if (rawpos) {delete [] rawpos; rawpos = 0;}
00175 if (M.rawpos) {
00176 if(maxorder>0){
00177 rawpos = new unsigned [maxorder]; // this is just a working space
00178 }
00179 else {
00180 rawpos = 0;
00181 }
00182 }
00183 pattern.resize(npattern);
00184 std::set<std::string>::iterator pos;
00185 for (i=0, pos=pattern_tree.begin(); pos != pattern_tree.end(); ++pos,++i) {
00186 pattern[i] = *pos;
00187 }
00188 pop_created = M.pop_created;
00189 if (M.pop_created) {
00190 pop = new Population(*(M.pop));
00191 }
00192 else {
00193 pop = M.pop;
00194 }
00195 if (rve) {
00196 delete [] rve;
00197 rve = 0;
00198 }
00199 if (M.rve) {
00200 if(npattern>0){
00201 rve = new Matrix<double> [npattern];
00202 }
00203 else {
00204 rve = 0;
00205 }
00206 check_ptr(rve);
00207 for (i=0; i<npattern; i++) rve[i].reserve(numtrait,numtrait);
00208 }
00209 if (factor_struct) {
00210 delete [] factor_struct; factor_struct = 0;
00211 }
00212 if (M.factor_struct) {
00213 if(numfact>0){
00214 factor_struct = new FieldStruct [numfact];
00215 }
00216 else {
00217 factor_struct = 0;
00218 }
00219 for (i=0; i<numfact; i++) factor_struct[i] = M.factor_struct[i];
00220 }
00221 if (trait_struct) {
00222 delete [] trait_struct; trait_struct = 0;
00223 }
00224 if (M.trait_struct) {
00225 if(numtrait>0){
00226 trait_struct = new FieldStruct* [numtrait];
00227 }
00228 else {
00229 trait_struct = 0;
00230 }
00231 for (i=0; i<numtrait; i++) {
00232 k = fieldindex(traitname[i],factor_struct,numfact);
00233 if (k >= 0 && k < numfact) {
00234 trait_struct[i] = &(factor_struct[k]);
00235 }
00236 else {
00237 throw exception("Model::copyfrom(): you have probably found a bug!");
00238 }
00239 }
00240 }
00241 if (xact_htable) {delete [] xact_htable; xact_htable = 0;}
00242 if (M.xact_htable) {
00243 if(numterm>0){
00244 xact_htable = new HashTable [numterm];
00245 }
00246 else {
00247 xact_htable = 0;
00248 }
00249 for (i=0; i<numterm; i++) xact_htable[i] = M.xact_htable[i];
00250 }
00251 if (idlist) {delete [] idlist; idlist = 0;}
00252 if (M.idlist) {
00253 if(numcol>0){
00254 idlist = new HashTable* [numcol];
00255 }
00256 else {
00257 idlist = 0;
00258 }
00259 check_ptr(idlist);
00260 if (data) for (i=0; i<numcol; i++) idlist[i] = data->hashtable[i];
00261 }
00262 }
|
|
||||||||||||||||||||
|
Definition at line 2059 of file model.cpp.
02063 :
02064 // *
02065 // * M.variance("animal",1.0,...)
02066 // *********************************************************************
02067 {
02068 if (numtrait<1) {
02069 std::cerr << "Model::variance(args): no trait(dependent-variable) in the model\n";
02070 return;
02071 }
02072 doubleMatrix var(numtrait,numtrait);
02073 int t1,t2;
02074 va_list param_pt; // an object param_pt
02075 va_start(param_pt,v00); // call the setup macro
02076 var[0][0] = v00;
02077 for (t2=1; t2<numtrait; t2++) var[0][t2] = va_arg(param_pt,double);
02078
02079 for (t1=1; t1<numtrait; t1++) for (t2=0; t2<numtrait; t2++) {
02080 var[t1][t2] = va_arg(param_pt,double);
02081 }
02082 va_end(param_pt);
02083 covariance(termname,termname2,var);
02084 }
|
|
||||||||||||||||
|
Definition at line 1947 of file model.cpp. References matvec::Matrix< double >::clear(), matvec::Matrix< double >::copy(), matvec::Matrix< double >::me, matvec::Vector< T >::reserve(), and matvec::Vector< T >::resize().
01952 {
01953 if (!factor_struct) {
01954 std::cerr << "Model::covariance(args): no equation(s) in the model\n";
01955 return;
01956 }
01957 if (v.num_rows() != v.num_cols() || v.num_rows() != numtrait) {
01958 std::cerr << "Model::covariance(): invalid variance matrix size\n";
01959 return;
01960 }
01961 int i,k,k2;
01962 k = term.index(termname,factor_struct);
01963
01964 k2 = term.index(termname2,factor_struct);
01965
01966 if( k < 0 ){
01967 std::cerr << "Model::covariance(..): not in model: " << termname << "\n";
01968 return;
01969 }
01970 if( k2 < 0 ){
01971 std::cerr << "Model::covariance(..): not in model: " << termname2 << "\n";
01972 return;
01973 }
01974 int bek=base_effect[k];
01975 if(term[bek].classi() != 'R' && term[bek].classi() != 'P'){
01976 std::cerr << "Model::covariance(..): not a random effect: " << termname << "\n";
01977 return;
01978 }
01979 int bek2=base_effect[k2];
01980 if(term[bek2].classi() != 'R' && term[bek2].classi() != 'P'){
01981 std::cerr << "Model::covariance(..): not a random effect: " << termname2 << "\n";
01982 return;
01983 }
01984 if(term[bek].classi() != term[bek2].classi() ){
01985 std::cerr << "Model::covariance(..): different types of random effects: \n";
01986 return;
01987 }
01988
01989
01990 int be;
01991 if(base_effect[k]==base_effect[k2]){
01992 be=base_effect[k];
01993 }
01994 else {
01995 be=base_effect[k];
01996 int be2=base_effect[k2];
01997 if(be2 <be) {
01998 int tmp;
01999 tmp=k;
02000 k=k2;
02001 k2=tmp;
02002 tmp=be;
02003 be=be2;
02004 be2=tmp;
02005 }
02006 base_effect[k2]=be;
02007 doubleMatrix orig1,orig2;
02008 orig1.copy(*term[be].prior->var_matrix());
02009 orig2.copy(*term[be2].prior->var_matrix());
02010 int newsize=nt_vec[be]+nt_vec[be2];
02011 Vector<int> orig=trait_map[be];
02012 trait_map[be].reserve(newsize);
02013 int ii=0;
02014 for(int i=0;i<nt_vec[be];i++,ii++) trait_map[be][ii]=orig[i];
02015 for(int i=0;i<nt_vec[be2];i++,ii++) trait_map[be][ii]=trait_map[be2][i];
02016 term[be].prior->var_matrix()->resize(newsize,newsize);
02017
02018 for(int i=0;i<nt_vec[be];i++){
02019 for(int j=0;j<nt_vec[be];j++){
02020 term[be].prior->var_matrix()->me[i][j]=orig1[i][j];
02021 }
02022 }
02023
02024 for(int i=0;i<nt_vec[be2];i++){
02025 for(int j=0;j<nt_vec[be2];j++){
02026 term[be].prior->var_matrix()->me[nt_vec[be]+i][nt_vec[be]+j]=orig2[i][j];
02027 }
02028 }
02029 nt_vec[be]=newsize;
02030 term[be2].prior->var_matrix()->clear();
02031 nt_vec[be2]=0;
02032 term[be2].classi('F');
02033 int last_eff=be;
02034 for(;corr_link[last_eff] != -1;)last_eff=corr_link[last_eff];
02035 corr_link[last_eff]=be2;
02036 int last_pos=pos_vec[last_eff];
02037 for(;corr_link[last_eff] != -1;){
02038 last_eff=corr_link[last_eff];
02039 pos_vec[last_eff]=pos_vec[last_eff]+last_pos+1;
02040 }
02041 }
02042
02043 int srow,scol;
02044 srow=pos_vec[k]*numtrait;
02045 scol=pos_vec[k2]*numtrait;
02046 for(int irow=0;irow<numtrait;irow++){
02047 for(int jcol=0;jcol<numtrait;jcol++){
02048 term[be].prior->var_matrix()->me[srow+irow][scol+jcol]=v[irow][jcol];
02049 term[be].prior->var_matrix()->me[scol+jcol][srow+irow]=v[irow][jcol];
02050 }
02051 }
02052
02053
02054 // std::cout << "Cov " << *term[be].prior->var_matrix();
02055
02056
02057 }
|
|
||||||||||||||||
|
Definition at line 1098 of file glmm2.cpp. References matvec::warning().
01103 {
01104 if (!factor_struct) throw exception("Model::variance(args): no equation(s) in the model");
01105
01106 int k = term.index(termname,factor_struct);
01107 int k2 = term.index(termname2,factor_struct);
01108 if(k2<k) {
01109 int kt=k;
01110 k=k2;
01111 k2=kt;
01112 }
01113 if ( k < 0) {
01114 warning("GLMM::covariance(): %s: not in the model",termname.c_str());
01115 type = bad_model;
01116 return;
01117 }
01118 if ( k2< 0) {
01119 warning("GLMM::covariance(): %s: not in the model",termname2.c_str());
01120 type = bad_model;
01121 return;
01122 }
01123 if (term[k].classi() != 'P') {
01124 warning("GLMM::covariance(): %s: no associated Pedigree",termname.c_str());
01125 type = bad_model;
01126 return;
01127 } // factor classi will override
01128 if (term[k2].classi() != 'P') {
01129 warning("GLMM::covariance(): %s: no associated Pedigree",termname2.c_str());
01130 type = bad_model;
01131 return;
01132 } // factor classi will override
01133 if(!corrmap[k]) {
01134 corrmap[k]=1;
01135 ncorr++;
01136 }
01137 if(!corrmap[k2]) {
01138 corrmap[k2]=1;
01139 ncorr++;
01140 }
01141 if (v.num_rows() == v.num_cols() && v.num_rows() == numtrait) {
01142 corrvar[k][k2]=v;
01143 }
01144 else {
01145 warning("Model::variance(): invalid variance matrix size");
01146 type = bad_model;
01147 }
01148 }
|
|
|
Definition at line 2097 of file model.cpp. References matvec::FieldStruct::classi(), matvec::Model::data, matvec::Model::factor_struct, matvec::fieldindex(), matvec::Model::numfact, matvec::split(), and matvec::Model::term.
02098 {
02099 // *********************************************************************
02100 // it specifies a list of effect (factor) as covariates.
02101 // note that an effect is not a model-term, which consists of effects,
02102 // instead an effect is associated with a data column. It is required
02103 // that the name of an effect must be the same as the associated
02104 // data-column.
02105 // **********************************************************************
02106 if (!factor_struct) {
02107 std::cerr << "Model::covariate(): no equation(s) in the model\n";
02108 return;
02109 }
02110
02111 int j,k,t;
02112 unsigned i,nc;
02113 std::vector<std::string> efflist;
02114 nc = split(covariate_names," ,",&efflist); // efflist = covariate_names.split(nc,", ");
02115
02116 for (i=0; i<nc; i++) {
02117 k = fieldindex(efflist[i],factor_struct,numfact);
02118 if (k < 0) throw exception("Model::covariate(): not in the model, thus it's ignored");
02119 // need to check if classi is changed
02120 if (factor_struct[k].classi() != 'F') data= 0;
02121 factor_struct[k].classi('C');
02122
02123 ////////////////////////////////////////////////////////////////////
02124 // if a term contains at least one covariate, then its classi = 'C'
02125 // else leave it alone
02126 //////////////////////////////////////////////////////////////////
02127 for (t=0; t<numterm; t++) {
02128 j = term[t].partial_match(efflist[i].c_str(),factor_struct);
02129 if (j >=0 ) {
02130 term[t].classi('C');
02131 break;
02132 }
02133 }
02134 }
02135 }
|
|
||||||||||||
|
Definition at line 128 of file model_hypo.cpp. References matvec::Model::CovMat(), matvec::Model::factor_struct, matvec::Model::get_blupsol(), matvec::Model::hmmesize, matvec::TermList::index(), matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), matvec::Model::numtrait, and matvec::Model::term.
00128 {
00129 doubleMatrix retval;
00130 if (!get_blupsol()) throw exception("Model::CovMat(): bad model");
00131
00132 int k = term.index(termname,factor_struct);
00133 if (k < 0) throw exception("Model::CovMat(): no such term in the model");
00134
00135 unsigned nr = Kp.num_rows();
00136 unsigned nc = Kp.num_cols();
00137 unsigned i,startaddr;
00138 if (nc <= term[k].nlevel()*numtrait) {
00139 startaddr = term[k].startaddr();
00140 doubleMatrix fullsize_Kp(nr,hmmesize);
00141 for (i=0; i<nr; i++) {
00142 memcpy(&(fullsize_Kp[i][startaddr]),Kp[i], sizeof(double)*nc);
00143 }
00144 retval=CovMat(fullsize_Kp);
00145 }
00146 else {
00147 throw exception("Model::CovMat(): too many columns in Kp");
00148 }
00149 return retval;
00150 }
|
|
|
Definition at line 69 of file model_hypo.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< double >::begin(), matvec::Model::get_blupsol(), matvec::Model::hmmec, matvec::Model::hmmesize, matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), and matvec::SparseMatrix::solve(). Referenced by matvec::Model::CovMat().
00069 {
00070 unsigned nr=Kp.num_rows();
00071 unsigned nc = Kp.num_cols();
00072 double *kp,*sol;
00073 doubleMatrix KvK(nr,nr);
00074 if (!get_blupsol("ysmp")) throw exception("Model::CovMat(): bad model");
00075 if (nc > hmmesize) throw exception("Model::CovMat(): size not conformable");
00076
00077 Vector<double> xy(hmmesize);
00078 Vector<double> tmpsol(hmmesize);
00079 double kpd,**kpme=Kp.begin();
00080 unsigned j;
00081 for (unsigned i=0; i<nr; i++) {
00082 if (nc < hmmesize) {
00083 memset(xy.begin(),'\0',sizeof(double)*hmmesize);
00084 memcpy(xy.begin(),kpme[i],sizeof(double)*nc);
00085 if (!hmmec.solve(tmpsol,xy,"ysmp")) throw exception("Model::CovMat(): hmmec.solve failed");
00086 }
00087 else {
00088 if (!hmmec.solve(tmpsol.begin(),kpme[i],"ysmp")) throw exception("Model::CovMat(): hmmec.solve failed");
00089 }
00090 for (unsigned k=0; k<nr; k++) {
00091 kp = Kp[k]; sol = tmpsol.begin();
00092 for (kpd=0.0,j=0; j<nc; j++) kpd += *kp++ * *sol++;
00093 KvK[k][i] = kpd;
00094 }
00095 }
00096 return(KvK);
00097 }
|
|
|
Definition at line 144 of file glmm.h. References matvec::Model::data_static.
00144 {if(DS_set != -1) data_static=DS_set;return(data_static);}
|
|
|
|
|
||||||||||||||||||||
|
Definition at line 1176 of file model.cpp. References matvec::Individual::display_freq_haplotype(), matvec::Population::member(), matvec::Population::MH_ibd_sample_map(), matvec::Model::pop, matvec::Population::size(), and matvec::Population::sum_descentState_map().
01176 {
01177 // Authors: L. Radu Totir
01178 // (August, 2004)
01179 // Contributors:
01180 double l_hood = 0.0;
01181 unsigned k = 0, j = 0, option = 0;
01182 unsigned k1 = numOfSL;
01183 unsigned k2 = numOfHaplo;
01184 unsigned k3 = numOfSLCas;
01185 cout<< endl << "Descent graph iterations = " << numOfSamples << endl;
01186 for(int sample=1; sample <= numOfSamples; sample++) {
01187 if(sample%5000 == 0) cout<<sample<<".";
01188 // cout<<endl<< "SAMPLE.................................. " << sample <<endl;
01189 if (k >= k1 + k2 + k3) {k = 0;}
01190 if (k < k1) {option = 1;}
01191 else if ((k >= k1) && (k < k1 + k2)) {option = 2;}
01192 else if (k <= k1 + k2 + k3) {option = 3;}
01193 k++;
01194 j++;
01195
01196 l_hood = pop->MH_ibd_sample_map(l_hood,option);
01197 pop->sum_descentState_map(); // update pdq's for all marker loci
01198 if (j >= 100000) {
01199 j = 0;
01200 cout << "SI_PDQ sample = " << sample << endl;
01201 }
01202 }
01203 cout << endl << "Haplotypes origins mm mp pm pp / mother - father " << endl;
01204 for (int i=0; i<pop->size(); i++) {
01205 pop->member(i)->display_freq_haplotype(numOfSamples);
01206 }
01207 cout << endl;
01208 }
|
|
||||||||||||
|
Definition at line 1153 of file model.cpp. References matvec::Population::descent_graph_setup(), matvec::Vector< unsigned >::empty(), matvec::Population::founder_allele_counter, matvec::Individual::m_counter_map, matvec::Population::member(), matvec::Individual::p_counter_map, matvec::Model::pop, matvec::Model::popsize, matvec::Vector< unsigned >::resize(), matvec::Individual::search_heteroz(), matvec::Individual::set_founder_alleles(), and matvec::Population::size().
01153 {
01154 // Authors: L. Radu Totir
01155 // (August, 2004)
01156 // Contributors:
01157 pop->descent_graph_setup(D);
01158 pop->founder_allele_counter = 0;
01159 int popsize = pop->size();
01160 int zero = 0;
01161 for (int i=0;i<popsize;i++) {
01162 pop->member(i)->set_founder_alleles(0);
01163 }
01164 unsigned lastLocus = numLoci;
01165 for (int i=0; i<popsize; i++) {
01166 if (pop->member(i)->m_counter_map.empty()){
01167 int zero = 0;
01168 pop->member(i)->m_counter_map.resize(numLoci,zero);
01169 pop->member(i)->p_counter_map.resize(numLoci,zero);
01170 }
01171 pop->member(i)->search_heteroz(lastLocus);
01172 //cout << pop->member(i)-> id() << " " << pop->member(i)-> ord_heter << endl;
01173 }
01174 }
|
|
|
Definition at line 324 of file model.cpp. References matvec::Model::modelstring.
00325 {
00326 std::cout << " model: " << modelstring << "\n";
00327 return 1;
00328 }
|
|
|
Reimplemented from matvec::Model. Definition at line 1507 of file glmm1.cpp.
01507 {
01508
01509
01510 int retval;
01511 retval=Model::equation(modelspecs);
01512
01513
01514 return(retval);
01515 }
|
|
||||||||||||||||||||
|
Definition at line 154 of file model_hypo.cpp. References matvec::warning().
00155 {
00156 ///////////////////////////////////////////
00157 // trailing zeros in K' can be omitted
00158 ///////////////////////////////////////////
00159 if (!kpme) {
00160 warning("Model::estimate(kp): kp is null, thus it is ignored");
00161 return;
00162 }
00163 if (!get_blupsol("ysmp")) throw exception("Model::estimate(): bad model");
00164
00165 Vector<double> xy;
00166 if (nc < hmmesize) xy.resize(hmmesize);
00167 for (unsigned i=0; i<nr; i++) {
00168 if (nc < hmmesize) {
00169 memcpy(xy.begin(),kpme[i],sizeof(double)*nc);
00170 kpb[i] = blupsol.inner_product(xy);
00171 }
00172 else {
00173 kpb[i] = blupsol.inner_product(kpme[i]);
00174 }
00175 }
00176 }
|
|
||||||||||||
|
Definition at line 99 of file model_hypo.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< T >::begin(), matvec::Model::estimate(), matvec::Model::factor_struct, matvec::Model::get_blupsol(), matvec::Model::hmmesize, matvec::TermList::index(), matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), matvec::Model::numtrait, matvec::Vector< T >::resize(), and matvec::Model::term.
00100 {
00101 Vector<double> retval;
00102 if (!get_blupsol()) throw exception("Model::estimate(): bad model");
00103
00104 int k = term.index(termname,factor_struct);
00105 if (k < 0) throw exception("Model::estimate(): no such term in the model");
00106
00107 unsigned nr = Kp.num_rows();
00108 unsigned nc = Kp.num_cols();
00109 unsigned i,startaddr;
00110 if (nc <= term[k].nlevel()*numtrait) {
00111 startaddr = term[k].startaddr();
00112 Matrix<double> fullsize_Kp(nr,hmmesize);
00113 for (i=0; i<nr; i++) {
00114 memcpy(&(fullsize_Kp[i][startaddr]),Kp[i], sizeof(double)*nc);
00115 }
00116 retval.resize(nr);
00117 estimate((const double **)fullsize_Kp.begin(),nr,hmmesize,retval.begin());
00118 }
00119 else {
00120 throw exception("Model::estimate(): too many columns in Kp");
00121 }
00122 return retval;
00123 }
|
|
|
Definition at line 55 of file model_hypo.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< double >::begin(), matvec::Model::estimate(), matvec::Matrix< double >::num_cols(), and matvec::Matrix< double >::num_rows().
00056 {
00057 ///////////////////////////////////////////
00058 // trailing zeros in K' can be omitted
00059 ///////////////////////////////////////////
00060
00061 unsigned nr = Kp.num_rows();
00062 unsigned nc = Kp.num_cols();
00063 Vector<double> kpb(nr);
00064 estimate((const double **)Kp.begin(),nr,nc,kpb.begin());
00065 return kpb;
00066 }
|
|
|
Definition at line 37 of file model_hypo.cpp. References matvec::Vector< T >::begin(), and matvec::Vector< T >::size(). Referenced by matvec::Model::estimate().
00038 {
00039 ///////////////////////////////////////////
00040 // trailing zeros in K' can be omitted
00041 ///////////////////////////////////////////
00042
00043 double **kpme = new double *[1];
00044 kpme[0] = Kp.begin();
00045 unsigned nc = Kp.size();
00046 double kpb[1];
00047 estimate((const double **)kpme,1,nc,kpb);
00048 if(kpme) {
00049 delete [] kpme;
00050 kpme = 0;
00051 }
00052 return kpb[0];
00053 }
|
|
|
Definition at line 113 of file model_peeling.cpp. References matvec::ChromStruct::locus.
00115 {
00116 if (type == bad_model) throw exception("Model::estimate_peeling(): bad model");
00117 double max_llhood = 0.0;
00118 if (!data_prepared) {
00119 if (!prepare_data()) {
00120 type = bad_model;
00121 return max_llhood;
00122 }
00123 }
00124
00125 minfun_indx = 2;
00126 pop->residual_var = residual_var;
00127 pop->input_data(data);
00128 pop->build_pop_gamete();
00129 unsigned i,n,m,k,j;
00130 m = pop->tn_gamete - 1; // because p(A1)+ p(A2) +...+ p(At) = 1
00131 n = 1 + m; // n = residual_var + allele_freq
00132 // n += pop->tn_genotype; // + genotypic_value // uncommented in old version
00133 Vector<double> x(n);
00134
00135 x[0] = pop->residual_var[0][0]; // residual_variance
00136
00137 ChromStruct *chrom = pop->prior->chrom();
00138 double* ve = chrom[0].locus[0].allele_freq.begin();
00139 for (i=0; i<m; i++) x[i+1] = ve[i]; // 1 is for residual_variance
00140
00141 // uncommented in old version
00142 // k = 1 + m;
00143 // const double** me = pop->prior->genotypic_val(0,0);
00144 // for (i=0; i<tn_genotype; i++) for (j=0; j<=i; j++) x[k++] = me[i][j];
00145 int iter = 0;
00146 max_llhood = praxis(x,n,iter);
00147 minfun_indx = 0;
00148 return max_llhood;
00149 }
|
|
|
Definition at line 304 of file model.cpp. References matvec::Model::bad_model, matvec::Model::data, matvec::Model::data_prepared, matvec::Model::hmmec, matvec::Model::hmmesize, matvec::Model::max_nz, matvec::Model::npattern, matvec::Data::num_rows(), matvec::Model::numcol, matvec::Model::numobs, matvec::Model::numrec, matvec::SparseMatrix::resize(), matvec::Model::type, and matvec::warning(). Referenced by matvec::Population::setupRSampler().
00305 {
00306 if (D.num_rows() > 0) {
00307 data = &D;
00308 }
00309 else {
00310 data = 0;
00311 type = bad_model;
00312 warning("Model::fitdata(D): Data is empty");
00313 }
00314 max_nz = 0;
00315 numrec = 0;
00316 numobs = 0;
00317 numcol = 0;
00318 npattern = 0;
00319 hmmesize = 0;
00320 data_prepared = 0;
00321 hmmec.resize(0,0);
00322 }
|
|
||||||||||||
|
Definition at line 109 of file model_gibbs.cpp. References matvec::Vector< T >::begin(), matvec::Vector< double >::begin(), matvec::Model::blupsol, matvec::Model::compute_xbzu(), matvec::Population::cond_genotype_config(), matvec::Model::data, matvec::Model::data_prepared, matvec::Population::gibbs_iterate(), matvec::SparseMatrix::gibbs_iterate(), matvec::Model::hmmec, matvec::Model::hmmesize, matvec::Population::input_data(), matvec::Model::ntermGdist, matvec::Model::pop, matvec::Model::prepare_data(), matvec::Model::rellrhs, matvec::Vector< T >::reserve(), matvec::Model::residual_var, matvec::Population::residual_var, matvec::Vector< double >::resize(), matvec::Population::resize_genotype_counter(), matvec::Model::setup_mme(), and matvec::warning().
00110 {
00111 if (ntermGdist==0) {
00112 warning(" Model::genotype_dist_gibbs(): no GeneticDist term");
00113 return 0.0;
00114 }
00115 if (!data_prepared) prepare_data();
00116 unsigned i,it;
00117 blupsol.resize(hmmesize,0.0);
00118 Vector<double> wspace; wspace.reserve(hmmesize);
00119 pop->resize_genotype_counter();
00120 pop->residual_var = residual_var;
00121 pop->input_data(data);
00122 pop->cond_genotype_config(); // an initial configuration
00123
00124 if (numterm) setup_mme(&rellrhs);
00125
00126 /////////////////////////////////////////////////////
00127 // Gibbs Sampler warm up
00128 /////////////////////////////////////////////
00129 for (it=0; it<warmup; it++) {
00130 if (numterm) hmmec.gibbs_iterate(blupsol.begin(),rellrhs.begin(),wspace.begin());
00131 for (i=0; i<ntermGdist; i++) {
00132 compute_xbzu(blupsol); // to adjust trait
00133 pop->gibbs_iterate();
00134 }
00135 }
00136
00137 /////////////////////////////////////////////////////
00138 // samples from Gibbs Sampler now will be used
00139 ////////////////////////////////////////////////////
00140
00141 for (it=0; it<gibbslen; it++) {
00142 if (numterm) hmmec.gibbs_iterate(blupsol.begin(),rellrhs.begin(),wspace.begin());
00143 for (i=0; i<ntermGdist; i++) {
00144 compute_xbzu(blupsol); // to adjust trait
00145 pop->gibbs_iterate(1);
00146 }
00147 }
00148 hmmec.resize(0,0);
00149 ntermGdist = 0;
00150 return 0.0; // no useful value could be returned so far
00151 }
|
|
||||||||||||
|
Definition at line 29 of file model_peeling.cpp. References matvec::Model::bad_model, matvec::Model::blupsol, matvec::Model::compute_xbzu(), matvec::Model::data, matvec::Model::data_prepared, matvec::Population::genotype_dist_peeling(), matvec::Model::hmmec, matvec::Model::hmmesize, matvec::Population::input_data(), matvec::Model::ntermGdist, matvec::Model::pop, matvec::Model::prepare_data(), matvec::Model::rellrhs, matvec::Model::residual_var, matvec::Population::residual_var, matvec::SparseMatrix::resize(), matvec::Vector< double >::resize(), matvec::Model::setup_mme(), matvec::SparseMatrix::solve(), and matvec::Model::type.
00030 {
00031 //////////////////////////////////////////////////////////////////////////
00032 // this works for arbitrary pedigrees
00033 // based on iterative peeling algorithm: Arendonk(1989), Fernando(1993)
00034 // prtl = 0, print nothing
00035 // 1, default, print genotype distribution and genotypic values
00036 // estifreq = 0, not estimate allele frequencies
00037 // 1, default, estimate allele frequencies
00038 //////////////////////////////////////////////////////////////////////////
00039 if (type == bad_model) return;
00040 if (!data_prepared) {
00041 if (!prepare_data()) {
00042 type = bad_model;
00043 return;
00044 }
00045 }
00046
00047 pop->residual_var = residual_var;
00048 pop->input_data(data);
00049
00050 blupsol.resize(hmmesize);
00051
00052 if (numterm) {
00053 setup_mme(&rellrhs);
00054 hmmec.solve(blupsol,rellrhs);
00055 }
00056 if (ntermGdist) {
00057 compute_xbzu(blupsol); // to adjust trait
00058 pop->genotype_dist_peeling(prtl,estifreq);
00059 }
00060
00061 hmmec.resize(0,0);
00062 ntermGdist = 0;
00063 }
|
|
||||||||||||
|
Definition at line 77 of file model_gibbs.cpp. References matvec::Vector< double >::begin(), matvec::Vector< T >::begin(), matvec::Model::blupsol, matvec::Model::data_prepared, matvec::SparseMatrix::gibbs_iterate(), matvec::Model::hmmec, matvec::Model::hmmesize, matvec::Model::prepare_data(), matvec::Model::rellrhs, matvec::Vector< T >::reserve(), matvec::Vector< double >::resize(), matvec::Model::sampleW(), and matvec::Model::setup_mme(). Referenced by matvec::Model::genotypic_value_gibbs().
00079 {
00080 if (!data_prepared) prepare_data();
00081
00082 blupsol.resize(hmmesize,0.0);
00083 Vector<double> tmpsol(hmmesize);
00084 Vector<double> wspace; wspace.reserve(hmmesize);
00085 setup_mme(&rellrhs); // relrhs is used here
00086
00087 /////////////////////////////////////////////////////
00088 // Gibbs Sampler warm up
00089 /////////////////////////////////////////////
00090 unsigned it;
00091 for (it=0; it<warmup; it++) {
00092 sampleW(tmpsol,rellrhs);
00093 hmmec.gibbs_iterate(tmpsol.begin(),rellrhs.begin(),wspace.begin());
00094 }
00095
00096 /////////////////////////////////////////////////////
00097 // samples from Gibbs Sampler now will be used
00098 ////////////////////////////////////////////////////
00099 for (it=0; it<gibbslen; it++) {
00100 sampleW(tmpsol,rellrhs);
00101 hmmec.gibbs_iterate(tmpsol.begin(),rellrhs.begin(),wspace.begin());
00102 blupsol += tmpsol;
00103 }
00104 blupsol /= gibbslen;
00105 hmmec.resize(0,0);
00106 return &blupsol;
00107 }
|
|
||||||||||||
|
Definition at line 153 of file model_gibbs.cpp. References matvec::Vector< T >::assign(), matvec::Vector< double >::begin(), matvec::Vector< T >::begin(), matvec::Model::blupsol, matvec::Model::compute_xbzu(), matvec::Population::cond_genotype_config(), matvec::Model::data, matvec::Model::data_prepared, matvec::Model::genotypic_value_cat(), matvec::Population::gibbs_iterate(), matvec::SparseMatrix::gibbs_iterate(), matvec::Model::hmmec, matvec::Model::hmmesize, matvec::Population::input_data(), matvec::Model::ntermGdist, matvec::Model::pop, matvec::Model::popsize, matvec::Model::prepare_data(), matvec::Model::rellrhs, matvec::Vector< T >::reserve(), matvec::Model::residual_var, matvec::Population::residual_var, matvec::Vector< double >::resize(), and matvec::Model::setup_mme().
00155 {
00156 if (categorical) return genotypic_value_cat(warmup,gibbslen);
00157 if (!data_prepared) prepare_data();
00158 unsigned i,it;
00159 if (numterm) setup_mme(&rellrhs);
00160 blupsol.resize(hmmesize + ntermGdist*popsize, 0.0);
00161 Vector<double> tmpsol(hmmesize);
00162 tmpsol.assign(0.0);
00163 Vector<double> wspace; wspace.reserve(hmmesize);
00164 double * ve = &(blupsol[hmmesize]);
00165 if (ntermGdist) {
00166 pop->residual_var = residual_var;
00167 pop->input_data(data);
00168 pop->cond_genotype_config();
00169 }
00170
00171 /////////////////////////////////////////////////////
00172 // Gibbs Sampler warm up
00173 /////////////////////////////////////////////
00174 for (it=0; it<warmup; it++) {
00175 if (numterm) hmmec.gibbs_iterate(tmpsol.begin(),rellrhs.begin(),wspace.begin());
00176 for (i=0; i<ntermGdist; i++) {
00177 compute_xbzu(tmpsol); // to adjust trait
00178 pop->gibbs_iterate();
00179 }
00180 }
00181
00182 /////////////////////////////////////////////////////
00183 // samples from Gibbs Sampler now will be used
00184 ////////////////////////////////////////////////////
00185
00186 for (it=0; it<gibbslen; it++) {
00187 if (numterm) {
00188 hmmec.gibbs_iterate(tmpsol.begin(),rellrhs.begin(),wspace.begin());
00189 for (i=0; i<hmmesize; i++) blupsol[i] += tmpsol[i];
00190 }
00191 for (i=0; i<ntermGdist; i++) {
00192 compute_xbzu(tmpsol); // to adjust trait
00193 pop->gibbs_iterate();
00194 for (i=0; i<popsize; i++) {
00195 ve[i] += pop->popmember[i]->genotypic_val();
00196 }
00197 }
00198 }
00199 blupsol /= gibbslen;
00200 hmmec.resize(0,0);
00201 return &(blupsol);
00202 }
|
|
|
Definition at line 65 of file model_peeling.cpp. References matvec::Model::bad_model, matvec::Model::blupsol, matvec::Model::compute_xbzu(), matvec::Model::data, matvec::Model::data_prepared, matvec::Individual::est_GV, matvec::Population::genotype_dist_peeling(), matvec::Model::get_blupsol(), matvec::Model::hmmec, matvec::Model::hmmesize, matvec::Population::input_data(), matvec::Model::ntermGdist, matvec::Model::pop, matvec::Population::popmember, matvec::Model::popsize, matvec::Model::prepare_data(), matvec::Vector< double >::print(), matvec::Model::rellrhs, matvec::Model::residual_var, matvec::Population::residual_var, matvec::Vector< double >::resize(), matvec::Model::setup_mme(), matvec::SparseMatrix::solve(), and matvec::Model::type.
00066 {
00067 if (type == bad_model) return 0;
00068 if (!data_prepared) {
00069 if (!prepare_data()) {
00070 type = bad_model;
00071 return 0;
00072 }
00073 }
00074
00075 blupsol.resize(hmmesize + ntermGdist*popsize);
00076 unsigned i;
00077
00078 if (ntermGdist ) {
00079 pop->residual_var = residual_var;
00080 pop->input_data(data);
00081 if (numterm) {
00082 get_blupsol();
00083 compute_xbzu(blupsol); // to adjust trait
00084 }
00085 else {
00086 pop->genotype_dist_peeling(0);
00087 for (i=0; i<popsize; i++) blupsol[i] = pop->popmember[i]->est_GV;
00088 return &(blupsol);
00089 }
00090 }
00091 else {
00092 setup_mme(&rellrhs);
00093 hmmec.solve(blupsol,rellrhs);
00094 return &(blupsol);
00095 }
00096
00097 ///////////////////////////////////////////
00098 // iterate begins now
00099 ////////////////////////////////////////////
00100
00101 double * ve = &(blupsol[hmmesize]);
00102 unsigned niterate = 0;
00103 while (niterate++ < 4) {
00104 pop->genotype_dist_peeling(0);
00105 for (i=0; i<popsize; i++) ve[i] = pop->popmember[i]->est_GV;
00106 blupsol.print();
00107 setup_mme(&rellrhs);
00108 hmmec.solve(blupsol,rellrhs);
00109 }
00110 return &(blupsol);
00111 }
|
|
||||||||||||||||||||
|
Return the BLUP solution The difference between blup() and get_blupsol() is that blup() will always recompute blup again, while get_blupsol() gets whatsoever in blupsol, Of course, if it's empty then it will compute blup solution
Definition at line 257 of file model_blup.cpp. References matvec::Model::bad_model, matvec::Model::blupsol, matvec::Model::hmmec, matvec::SparseMatrix::nz(), matvec::Model::rellrhs, matvec::Model::setup_mme(), matvec::SparseMatrix::solve(), and matvec::Model::type. Referenced by matvec::Model::contrast(), contrast(), matvec::Model::CovMat(), matvec::Model::estimate(), matvec::Model::genotypic_value_peeling(), and matvec::Model::get_solutions().
00259 {
00260 Vector<double> *retval = 0;
00261 if (type == bad_model) throw exception("Model::blup(): bad model");
00262 if (hmmec.nz() == 0) {
00263 setup_mme(&rellrhs);
00264 if (type == bad_model) throw exception("Model::blup(): bad model");
00265 if (hmmec.solve(blupsol,rellrhs,method,relax,stopval,mxiter)==1) {
00266 retval = &blupsol;
00267 //hmmec.close();
00268 }
00269 }
00270 else {
00271 retval = &blupsol;
00272 }
00273 return retval;
00274 }
|
|
|
|
|
||||||||||||||||||||||||
|
Definition at line 448 of file model_hypo.cpp. References matvec::SparseMatrix::a(), matvec::Vector< T >::assign(), matvec::Vector< T >::begin(), matvec::Session::epsilon, matvec::Model::hmmec, matvec::Model::hmmesize, matvec::SparseMatrix::ia(), matvec::SparseMatrix::ja(), matvec::Model::numterm, matvec::Model::numtrait, matvec::SESSION, and matvec::Model::term. Referenced by matvec::Model::out_lsmeans_to_stream().
00450 {
00451 int i,j,k,h;
00452 unsigned l,t,jb,n,nh,nl,term_start,term_nlevel,endaddr;
00453 int mcell;
00454
00455 k = term[termid].startaddr() + lvlth*numtrait + yth;
00456 int *ia = hmmec.ia();
00457 int *ja = hmmec.ja();
00458 double *a = hmmec.a();
00459 kp.assign(0.0);
00460 double *kpve = kp.begin();
00461 kpve--;
00462 mcell = 0;
00463 k++;
00464 endaddr = ia[k+1];
00465 for (j=ia[k]; j<endaddr; j++) {
00466 if (ja[j] != k) continue;
00467 if (a[j] <= SESSION.epsilon) {
00468 mcell = 1;
00469 return;
00470 }
00471 else {
00472 break;
00473 }
00474 }
00475
00476 ////////////////////////////////////////////////////////////////
00477 // set element kp[i] = 1.0 if MME[k][i] is non-zero element
00478 //
00479 ////////////////////////////////////////////////////////////////
00480 for (j=ia[k]; j<endaddr; j++) kpve[ja[j]] = 1.0;
00481 for (i=1; i<=hmmesize; i++) {
00482 endaddr = ia[i+1];
00483 for (j=ia[i]; j<endaddr; j++) if (ja[j] == k) kpve[i] = 1.0;
00484 }
00485
00486 for (t=0; t<numterm; t++) {
00487 if (!(term[t].classi() == 'F' || term[t].classi() == 'C')) continue;
00488 if (lsmtablevec[t] >= 1) continue;
00489 term_start = term[t].startaddr() + 1;
00490 term_nlevel = term[t].nlevel();
00491 l = term_start + yth;
00492 for (i=0; i<term_nlevel; i++) {
00493 for (n=ia[l+1], j=ia[l]; j<n; j++) {
00494 if (ja[j] == l) {
00495 if (fabs(a[j]) > SESSION.epsilon) kpve[l] = 1.0;
00496 }
00497 }
00498 l += numtrait;
00499 }
00500 }
00501
00502 ////////////////////////////////////////////////////////////////////
00503 // if term[i] is random, then set all elements kp[j] for term[i]
00504 // to be zeros,
00505 // else set elements kp[j] for non-relevent trait to be zero
00506 ////////////////////////////////////////////////////////////////
00507 for (i=0; i<numterm; i++) {
00508 term_start = term[i].startaddr() + 1;
00509 term_nlevel = term[i].nlevel();
00510 if (term[i].classi() == 'F' || term[i].classi() == 'C') {
00511 for (t=0; t<numtrait; t++) {
00512 if (t != yth) {
00513 l = term_start;
00514 for ( j=0; j<term_nlevel; j++) {
00515 kpve[l+t] = 0.0;
00516 l += numtrait;
00517 }
00518 }
00519 }
00520 }
00521 else {
00522 memset(&(kpve[term_start]),'\0',numtrait*sizeof(double)*term_nlevel);
00523 }
00524 }
00525
00526 double dx,xn;
00527 DataNode dnode;
00528 for (j=0; j<numterm; j++) {
00529 if (lsmtablevec[j] == 2) continue; // j = lsm_term
00530 term_nlevel = term[j].nlevel();
00531 term_start = term[j].startaddr() + 1;
00532 dx = 1.0;
00533 if (term[j].classi() == 'C') {
00534 k = factor_struct[term[j].factorindx[0]].index();
00535 dnode = data->datasheet[k].mean(0);
00536 if (dnode.missing) throw exception("get_lms_kp(): no valid values in the covariate column");
00537 dx = dnode.double_val();
00538 }
00539 if (!(term[j].classi() == 'F' || term[j].classi() == 'C')) continue;
00540 if (lsmtablevec[j] == 1) { // j is related to lsm_term
00541 n = term_start + term_nlevel*numtrait;
00542 for (xn=0.0,k=term_start; k<n; k++) xn += kpve[k];
00543 if (xn >= SESSION.epsilon*10) {
00544 for (k=term_start; k<n; k++) kpve[k] *= dx/xn;
00545 }
00546 }
00547 else {
00548 if (term[j].order() == 1) {
00549 n = term_start + term_nlevel*numtrait;
00550 for (xn=0.0, k=term_start; k<n; k++) xn += kpve[k];
00551 if (xn >= SESSION.epsilon*10) {
00552 for (k=term_start; k<n; k++) kpve[k] *= dx/xn;
00553 }
00554 }
00555 else {
00556 nl = factor_struct[term[j].factorindx[0]].nlevel();
00557 nh = term_nlevel/nl;
00558 for (k=0; k<nh; k++) {
00559 jb = term_start + k;
00560 for (xn=0.0,h=0; h<nl; h++) {
00561 xn += kpve[jb];
00562 jb += nh*numtrait;
00563 }
00564 if (xn < SESSION.epsilon*10) continue;
00565 jb = term_start + k;
00566 for (h=0; h<nl; h++) {
00567 kpve[jb] *= dx/(xn*nh);
00568 jb += nh*numtrait;
00569 }
00570 }
00571 }
00572 }
00573 }
00574 }
|
|
|
Definition at line 158 of file glmm.h.
00158 {return observation;}
|
|
|
Definition at line 339 of file model.h. References matvec::Model::popsize.
00339 {return popsize;};
|
|
|
Definition at line 159 of file glmm.h.
00159 {return SinMat;}
|
|
|
Definition at line 2390 of file model.cpp. References matvec::Model::blupsol, matvec::Model::factor_struct, matvec::Model::get_blupsol(), matvec::TermList::index(), matvec::Model::numtrait, matvec::Vector< T >::resize(), and matvec::Model::term.
02391 {
02392 // Authors: Liviu R. Totir and Rohan L. Fernando (November, 2002)
02393 // Contributors:
02394 Vector<double> retval;
02395 if (!get_blupsol()) throw exception("Model::get_solutions(): get_blupsol() failed");
02396
02397 int k = term.index(termname,factor_struct);
02398 if (k < 0) throw exception("Model::get_solutions(): no such term in the model");
02399 int nlevels = term[k].nlevel()*numtrait;
02400 unsigned i,startaddr;
02401 startaddr = term[k].startaddr();
02402 retval.resize(nlevels);
02403 for (i=0; i<nlevels; i++) {
02404 retval[i] = blupsol[startaddr+i];
02405 }
02406 return retval;
02407 }
|
|
|
Definition at line 1713 of file model.cpp. References matvec::Model::pop, matvec::Population::popmember, matvec::Model::popsize, matvec::Population::popsize, matvec::Matrix< T >::resize(), and matvec::Matrix< T >::transpose(). Referenced by matvec::Model::RSamplerGibbs(), matvec::Model::RSamplerGibbsMH(), and matvec::Model::RSamplerMH().
01713 {
01714 // Author: L. Radu Totir
01715 // (July, 2004)
01716 // Contributors:
01717 matvec::Matrix<double> geneticValVec;
01718 matvec::Matrix<double> CovMatrix;
01719 CovMatrix.resize(pop->popsize,pop->popsize,0.0);
01720 matvec::Vector<double> qtlMu;
01721 qtlMu.resize(4,0.0);
01722 qtlMu[0]=-1.41421,qtlMu[3]=1.41421;
01723 geneticValVec.resize(pop->popsize,1,0.0);
01724 for (unsigned t=0;t<pop->popsize;t++){
01725 // cout << (pop->popmember[t])->id() << " "
01726 // << (pop->popmember[t])->genotNodeVector[1].acceptedGenotypeState
01727 // << endl;
01728 int geno = (pop->popmember[t])->genotNodeVector[1].acceptedGenotypeState;
01729 geneticValVec[t][0] = qtlMu[geno];
01730 }
01731 CovMatrix =geneticValVec*geneticValVec.transpose();
01732 return CovMatrix;
01733 }
|
|
|
Definition at line 1866 of file glmm2.cpp. References matvec::Model::bad_model, matvec::ModelTerm::classi(), matvec::Model::factor_struct, matvec::TermList::index(), matvec::ModelTerm::prior, matvec::Model::term, matvec::Model::type, matvec::GeneticDist::var_matrix(), and matvec::warning().
01867 {
01868 // Authors: Liviu R. Totir and Rohan L. Fernando (March, 2004)
01869 // Contributors:
01870 doubleMatrix* retval = 0;
01871 if (type == bad_model) {
01872 warning("Model::info(stream): model is too bad");
01873 return *retval;
01874 }
01875 int k = term.index(termname,factor_struct);
01876 if (k < 0) throw exception("GLMM::getVarEstimates: no such term in the model");
01877 const ModelTerm *T;
01878 T = &(term(k));
01879 if (T->classi() == 'R' || (T->classi() == 'P')) {
01880 retval=T->prior->var_matrix();
01881 }
01882 return *retval;
01883 }
|
|
||||||||||||||||||||||||
|
Definition at line 142 of file glmm2.cpp. References matvec::Vector< T >::begin(), matvec::Model::blup(), matvec::Model::blupsol, matvec::SparseMatrix::dense(), matvec::doubleMatrix::ginv1(), matvec::Model::hmmec, matvec::Model::hmmesize, link_function, matvec::Model::mmec(), matvec::Model::rellrhs, matvec::Matrix< double >::resize(), matvec::Vector< T >::resize(), setup_mme(), matvec::SparseMatrix::solve(), and matvec::Vector< double >::subvec(). Referenced by AI_REML(), contrast(), and SSQCP().
00143 {
00144 int iter,i,j,k;
00145 doubleMatrix kvk;
00146 Vector<double> kpm,wrkvec;
00147 double old_like,new_like,*krow,*m,*tmpsol,kpd;
00148 double *kp,*sol;
00149 Vector<double> old_sol;
00150 Vector<double> *retval;
00151 if(!link_function) {
00152 blup("ysmp1");
00153 if(ke)
00154 {
00155 tmpsol=new double [hmmesize];
00156 krow=new double [hmmesize];
00157 }
00158 if(ke)
00159 {
00160 kpm.resize(nr);
00161 kvk.resize(nr,nr);
00162 if(mean) memcpy(kpm.begin(),mean,sizeof(double)*nr);
00163 for(i=0;i<nr;i++){
00164 for(j=0;j<nc;j++) kpm[i]-=ke[i][j]*blupsol[j];
00165 memset(krow,'\0',sizeof(double)*hmmesize);
00166 memcpy(krow,ke[i],sizeof(double)*nc);
00167 hmmec.solve(tmpsol,krow,"ysmp");
00168 for (k=0; k<nr; k++) {
00169 kp = ke[k]; sol = tmpsol;
00170 for (kpd=0.0,j=0; j<nc; j++) kpd += *kp++ * *sol++;
00171 kvk[k][i] = kpd;
00172 }
00173 }
00174 kvk.ginv1();
00175 kpm=kvk*kpm;
00176 memset(krow,'\0',sizeof(double)*hmmesize);
00177 for(i=0;i<nr;i++) {
00178 for(j=0;j<nc;j++) {
00179 krow[j]+=ke[i][j]*kpm[i];
00180 }
00181 }
00182 hmmec.solve(tmpsol,krow,"ysmp");
00183 for(i=0;i<hmmesize;i++) blupsol[i]+=tmpsol[i];
00184
00185 delete [] krow;
00186 delete [] tmpsol;
00187 }
00188
00189
00190
00191 return(&blupsol);
00192 }
00193 for(iter=0;iter< iterations;iter++) {
00194 setup_mme(&rellrhs);
00195 #ifdef DEBUG
00196 doubleMatrix LHS;
00197 LHS=mmec()->dense();
00198 std::cout << "\nRHS\n"<<rellrhs.subvec()<<"\n";
00199 std::cout << "\nLHS\n"<<LHS <<"\n";
00200 #endif
00201 blup("ysmp1");
00202
00203 if(ke && iter==0)
00204 {
00205 tmpsol=new double [hmmesize];
00206 krow=new double [hmmesize];
00207 }
00208 if(ke)
00209 {
00210 kpm.resize(nr);
00211 kvk.resize(nr,nr);
00212 if(mean) memcpy(kpm.begin(),mean,sizeof(double)*nr);
00213 for(i=0;i<nr;i++){
00214 for(j=0;j<nc;j++) kpm[i]-=ke[i][j]*blupsol[j];
00215 memset(krow,'\0',sizeof(double)*hmmesize);
00216 memcpy(krow,ke[i],sizeof(double)*nc);
00217 hmmec.solve(tmpsol,krow,"ysmp");
00218 for (k=0; k<nr; k++) {
00219 kp = ke[k]; sol = tmpsol;
00220 kpd=0.0;
00221 for (j=0; j<nc; j++) kpd += ke[k][j] * tmpsol[j];
00222 kvk[k][i] = kpd;
00223 }
00224 }
00225 kvk.ginv1();
00226 kpm=kvk*kpm;
00227 memset(krow,'\0',sizeof(double)*hmmesize);
00228 for(i=0;i<nr;i++) {
00229 for(j=0;j<nc;j++) {
00230 krow[j]+=ke[i][j]*kpm[i];
00231 }
00232 }
00233 hmmec.solve(tmpsol,krow,"ysmp");
00234 for(i=0;i<hmmesize;i++) blupsol[i]+=tmpsol[i];
00235 }
00236
00237 }
00238
00239
00240 retval=&blupsol;
00241 if(ke)
00242 {
00243 delete [] krow;
00244 delete [] tmpsol;
00245 }
00246 return retval;
00247 }
|
|
|
Definition at line 1155 of file model_peeling.cpp. References matvec::Population::analysis_type, matvec::Model::bad_model, |