---

matvec::Model Class Reference

#include <model.h>

Inheritance diagram for matvec::Model:

List of all members.

---

Detailed Description

a generalized linear mixed model

term[t].classi =

'F'	a fixed term, including intercept
'C'	a covariate term
'R'	a regular random term
'P'	a random term associated with a pedigree

  factor_struct[i].classi = 'F', a fixed factor
                          = 'I', intercept
                          = 'C', covariate factor
                          = 'P', a factor associated with a pedigree

  weightinverse  = 0, no inverse is needed (default)
                 = 1, inverse is needed
 * 

See also:

Data

Definition at line 88 of file model.h.

Public Types
enum	ModelType { bad_model, fixed_model, mixed_model, reg_model }
Public Methods
	Model (void)
	Model (const std::string &modelspecs)
	Model (const Model &M)
virtual	~Model (void)
int	display (void) const
virtual int	prepare_data (const std::string &solver="ysmp")
unsigned	nonzero (void) const
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	nonzero (const unsigned nz)
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	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	RSampler (unsigned pedBlockSize, unsigned numLoci, unsigned numOfSamples, const std::string &samplerType, const std::string &howToSample, const std::string &samplerUsed, const std::string &whatToCompute)
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)
int	equation (const std::string &modelspecs)
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	save (const std::string &fname, const int io_mode=std::ios::out)
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 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)
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)
virtual SparseMatrix *	setup_mme (Vector< double > *rhs=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
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
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 *	kvec
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

---

Member Enumeration Documentation

enum matvec::Model::ModelType

Enumeration values: bad_model  fixed_model  mixed_model  reg_model  Definition at line 173 of file model.h. {bad_model, fixed_model, mixed_model, reg_model};

---

Constructor & Destructor Documentation

matvec::Model::Model (  void    )  [inline]

Definition at line 184 of file model.h. References initialize(). {initialize();}

matvec::Model::Model (  const std::string &    modelspecs )  [inline]

Definition at line 185 of file model.h. References equation(), and initialize(). {initialize(); equation(modelspecs);}

matvec::Model::Model (  const Model &    M )  [inline]

Definition at line 186 of file model.h. References copyfrom(), and initialize(). :Minimizer() {initialize(); copyfrom(M);}

virtual matvec::Model::~Model (  void    )  [inline, virtual]

Definition at line 187 of file model.h. References release(). {release();}

---

Member Function Documentation

void matvec::Model::add_Ag (  int    t, Vector< double > *    x_p = 0 )  [protected]

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(), hmmec, matvec::Individual::id(), matvec::SparseMatrix::insert(), lng0vec, matvec::Population::member(), mme_times_res, matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), nt_vec, numgroup, pop, popsize, matvec::SESSION, term, and var_link. Referenced by add_G_1(), and mme_times(). { unsigned t1, t2, k, i, j, ii, jj, iii, jjj; double dii,val,value; double *x,*r; r = mme_times_res.begin() - 1; if (x_p) x = x_p->begin() - 1; double **var = 0; doubleMatrix *tmp = term[var_link[t]].prior->var_matrix(); if (!tmp) throw exception("Model::add_Ag(): you probably forgot to use M.variance(...)\n" " to set variance for each random effect\n"); var = tmp->begin(); int nt=nt_vec[var_link[t]]; Matrix<double> rvarg(nt,nt); for (i=0; i<nt;i++) for (j=0; j<nt; j++) rvarg[i][j]= var[i][j]; ginverse1(rvarg.begin(),nt,lng0vec[t],1,SESSION.epsilon); Matrix<double> lambda(3,3); getlambda(lambda.begin(),3); unsigned startaddr = term[t].start + 1; unsigned asd[3]; Individual *I; unsigned nanim=popsize-numgroup; double f1,f2; for (k=0; k<nanim; k++) { I = pop->member(k); // cout << popsize << " " << " "<<nanim << " " << k << I << "\n"; //cout << I->id() << " " << I->father_id() << " " << I->mother_id() << "\n"; asd[0] = I->id(); asd[1] = I->father_id(); asd[2] = I->mother_id(); f1= I->father_inbcoef(); f2= I->mother_inbcoef(); if(asd[1]> nanim) f1=-1.; if(asd[2]> nanim) f2=-1.; dii = 4.0/(2.0 -f1 -f2); for (i=0; i<3; i++) { if (asd[i] == 0) continue; ii = startaddr + (asd[i]-1)*nt; for (j=0; j<3; j++) { if (asd[j] == 0) continue; jj = startaddr + (asd[j]-1)*nt; val = lambda[i][j]* dii; for (t1=0; t1<nt; t1++) { iii = ii + t1; for (t2=0; t2<nt; t2++) { jjj = jj + t2; if (jjj>=iii) { value = val*rvarg[t1][t2]; if (x_p) { // iteration on data r[iii] += value*x[jjj]; if (jjj > iii) r[jjj] += value*x[iii]; } else { hmmec.insert(iii,jjj,value); } } } } } } } }

void matvec::Model::add_Ag_diag (  int    t, double **    M )  [protected]

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(), lng0vec, matvec::Population::member(), matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), nt_vec, numgroup, pop, popsize, matvec::SESSION, term, and var_link. Referenced by add_G_1_diag(). { unsigned t1, t2,i,j,k,kk,rec; double dii,val,value; double **var = 0; doubleMatrix *tmp = term[var_link[t]].prior->var_matrix(); if (!tmp) throw exception("Model::add_Ag(): you probably forgot to use M.variance(...)\n" " to set variance for each random effect\n"); var = tmp->begin(); int nt=nt_vec[var_link[t]]; Matrix<double> rvarg(nt,nt); for (i=0; i<nt;i++) for (j=0; j<nt; j++) rvarg[i][j]= var[i][j]; ginverse1(rvarg.begin(),nt,lng0vec[t],1,SESSION.epsilon); Matrix<double> lambda(3,3); getlambda(lambda.begin(),3); unsigned startaddr = term[t].start + 1; unsigned asd[3]; Individual *I; unsigned nanim=popsize-numgroup; double f1,f2; for (k=0,i=0; i<t; ++i) k += term[i].nlevel(); for (rec=0; rec<nanim; rec++) { I = pop->member(rec); asd[0] = I->id(); asd[1] = I->father_id(); asd[2] = I->mother_id(); f1= I->father_inbcoef(); f2= I->mother_inbcoef(); if(asd[1]> nanim) f1=-1.; if(asd[2]> nanim) f2=-1.; dii = 4.0/(2.0 -f1 -f2); for (i=0; i<3; i++) { if (asd[i] == 0) continue; kk = k + asd[i] - 1; val = lambda[i][i]* dii; for (t1=0; t1<nt; ++t1) { for (t2=t1; t2<nt; ++t2) { M[kk][lidx(t2,t1,nt)+1] += val*rvarg[t1][t2]; } } } } }

void matvec::Model::add_G_1 (  void    )  [protected, virtual]

Add G-inverse Definition at line 279 of file model_blup.cpp. References add_Ag(), add_Ig(), numterm, and term. Referenced by setup_mme(). { int done_ped=0; for (int t=0; t<numterm; t++) { if (term[t].classi() == 'P') { add_Ag(t); //hmmec.display(12,24,12,24); } else if (term[t].classi() == 'R') { add_Ig(t); } } //hmmec.display(12,24,12,24); }

void matvec::Model::add_G_1_diag (  double **    M )  [protected, virtual]

add diagonals of G-inverse for pccg Definition at line 297 of file model_blup.cpp. References add_Ag_diag(), add_Ig_diag(), numterm, and term. Referenced by compute_rhs_diag_mme(). { int done_ped=0; for (int t=0; t<numterm; t++) { if (term[t].classi() == 'P') { add_Ag_diag(t,M); //hmmec.display(12,24,12,24); } else if (term[t].classi() == 'R') { add_Ig_diag(t,M); } } }

void matvec::Model::add_G_1_single_trait (  const Vector< double > &    ratio )  [protected]

Definition at line 332 of file model_vce.cpp. References matvec::Matrix< T >::begin(), matvec::Individual::father_id(), matvec::Individual::father_inbcoef(), matvec::getlambda(), hmmec, matvec::Individual::id(), matvec::SparseMatrix::insert(), matvec::Individual::mother_id(), matvec::Individual::mother_inbcoef(), numterm, pop, matvec::Population::popmember, and term. Referenced by setup_ww_single_trait(), and vce_emreml_multi_trait(). { unsigned i,j,ii,jj,t,k,rec,nl,startaddr; double val,rr,xval; for (k=0,t=0; t<numterm; t++) { startaddr = term[t].startaddr(); nl = term[t].nlevel(); if (term[t].classi() == 'P') { Individual *I; unsigned asd[3]; val = ratio[k++]; Matrix<double> lambda(3,3); getlambda(lambda.begin(),3); for (rec=0; rec<nl; rec++) { I = pop->popmember[rec]; asd[0]=I->id(); asd[1]=I->father_id(); asd[2] = I->mother_id(); rr = val*4.0/(2.0-I->father_inbcoef()- I->mother_inbcoef()); for (i=0; i<3; i++) { if (asd[i] != 0) { ii = startaddr + asd[i]; for (j=0; j<3; j++) { if (asd[j] != 0) { jj = startaddr + asd[j]; if (jj >= ii) { xval = lambda[i][j]*rr; hmmec.insert(ii,jj,xval); } } } } } } } else if (term[t].classi() == 'R') { val = ratio[k++]; for (ii=startaddr+1,i=0; i<nl; i++,ii++) hmmec.insert(ii,ii,val); } } }

void matvec::Model::add_Ig (  int    t, Vector< double > *    x_p = 0 )  [protected]

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(), hmmec, matvec::SparseMatrix::insert(), lng0vec, mme_times_res, nt_vec, matvec::SESSION, term, and var_link. Referenced by add_G_1(), and mme_times(). { unsigned k,i,ii,jj,nt; int t1,t2; double *x,*r,value; r = mme_times_res.begin() - 1; if (x_p) x = x_p->begin() - 1; double **v = 0; doubleMatrix *tmp = term[var_link[t]].prior->var_matrix(); if (!tmp) throw exception("Model::add_Ig(): you probably forgot to use M.variance(...)\n" " to set variance for each random effect\n"); nt=nt_vec[var_link[t]]; v = tmp->begin(); Matrix<double> rv(nt,nt); for (t1=0; t1<nt; t1++) { for(t2=0; t2<nt; t2++) rv[t1][t2]=v[t1][t2]; } ginverse1(rv.begin(),nt,lng0vec[t],1,SESSION.epsilon); unsigned na = term[t].nlevel(); unsigned startaddr = term[t].start + 1; for (k=0; k<na; k++) { i = startaddr + k*nt; for (t1=0; t1<nt; t1++) { ii = i + t1; for (t2=0; t2<nt; t2++) { jj = i + t2; if (jj>=ii) { value = rv[t1][t2]; if(x_p) {// iteration on data r[ii] += value*x[jj]; if (jj > ii) r[jj] += value*x[ii]; } else { hmmec.insert(ii,jj,value); } } } } } }

void matvec::Model::add_Ig_diag (  int    t, double **    M )  [protected]

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(), lng0vec, nt_vec, matvec::SESSION, term, and var_link. Referenced by add_G_1_diag(). { unsigned rec,k,i,ii,jj,kk,nt; int t1,t2; double **v = 0; doubleMatrix *tmp = term[var_link[t]].prior->var_matrix(); if (!tmp) throw exception("Model::add_Ig(): you probably forgot to use M.variance(...)\n" " to set variance for each random effect\n"); nt=nt_vec[var_link[t]]; v = tmp->begin(); Matrix<double> rv(nt,nt); for (t1=0; t1<nt; t1++) { for(t2=0; t2<nt; t2++) rv[t1][t2]=v[t1][t2]; } ginverse1(rv.begin(),nt,lng0vec[t],1,SESSION.epsilon); unsigned na = term[t].nlevel(); for (k=0,i=0; i<t; ++i) k += term[i].nlevel(); for (rec=0; rec<na; rec++) { kk = k + rec; for (t1=0; t1<nt; ++t1) { for (t2=t1; t2<nt; ++t2) { M[kk][lidx(t2,t1,nt)+1] += rv[t1][t2]; } } } }

void matvec::Model::anova (  double &    ssm, double &    sse, unsigned &    rank_x, int &    dfe )

Definition at line 426 of file model_hypo.cpp. References blupsol, hmmec, matvec::Vector< double >::inner_product(), matvec::SparseMatrix::logdet(), numgroup, numterm, numtrait, rellrhs, residual_var, term, and yry. { //////////////////////////////////////////////////////// // calculates sums of squares for single trait model // the rank of the random part: //////////////////////////////////////////////////////// unsigned rank_mme=0,rank_z=0; char cl; for (unsigned i=0; i<numterm; i++) { cl = term[i].classi(); if (cl == 'P' || cl == 'R') rank_z += term[i].nlevel(); } rank_z = (rank_z - numgroup)*numtrait; hmmec.logdet(&rank_mme); double ve = residual_var[0][0]; double sst = yry*ve; ssm = blupsol.inner_product(rellrhs)*ve; sse = sst - ssm; rank_x = rank_mme - rank_z; dfe = static_cast<int>(numobs) - static_cast<int>(rank_x); }

void matvec::Model::assign_id_xact (  const Vector< bool > &    )  [protected]

Definition at line 967 of file model.cpp. References data, hashxact(), matvec::Data::num_rows(), numterm, matvec::HashTable::resize(), matvec::HashTable::size(), term, and xact_htable. Referenced by prepare_data(). { // ********************************************************** // get sequantial integer ids for interation term in model // *********************************************************** if (!data) return; unsigned nord, nrow_in_data = data->num_rows(); if (xact_htable) { delete [] xact_htable; xact_htable = 0; } if(numterm>0){ xact_htable = new HashTable [numterm]; } else { xact_htable = 0; } int t; for (t=0; t<numterm; t++) { if (term[t].order() > 1) { xact_htable[t].resize(nrow_in_data,term[t].order()*sizeof(unsigned)); } } hashxact(record_legality); for (t=0;t<numterm;t++) { nord = term[t].order(); if (nord>1) { term[t].numlevel = xact_htable[t].size(); xact_htable[t].resize(term[t].nlevel(),nord*sizeof(unsigned)); } } hashxact(record_legality); }

Vector< double > * matvec::Model::blue (  const std::string &    method = "ysmp", double    stopval = 0.001, double    relax = 1.0, unsigned    mxiter = 100000 )

Return BLUE solution for a fixed model. Parameters: method  choices are ioc,iod,direct stopval  stop criteria for solver ioc and iod relax  relaxation factor for solver ioc mxiter  maximum number of iterations allowed for solver ioc,ios Definition at line 238 of file model_blup.cpp. References blup(), fixed_model, and type. { Vector<double> *retval = 0; if (type != fixed_model) throw exception("Model::blue(): it must be a fixed model"); retval = blup(method,stopval,relax,mxiter); return retval; }

Vector< double > * matvec::Model::blup (  const std::string &    method = "ysmp", double    stopval = 0.001, double    relax = 1.0, unsigned    mxiter = 100000 )

Return the BLUP solution. Parameters: method  choices are ioc,iod,direct stopval  stop criteria for solver ioc and iod relax  relaxation factor for solver ioc mxiter  maximum number of iterations allowed for solver ioc,ios Definition at line 206 of file model_blup.cpp. References bad_model, blup_pccg(), blupsol, hmmec, matvec::SparseMatrix::nz(), rellrhs, setup_mme(), matvec::SparseMatrix::solve(), and type. Referenced by blue(), matvec::GLMM::glim(), and lsmeans(). { Vector<double> *retval = 0; if (type == bad_model) throw exception("Model::blup(): bad model"); if (method == "iod") { blup_pccg(stopval,mxiter); return &blupsol; } if (hmmec.nz() == 0) { try { setup_mme(&rellrhs); } catch (exception &er) { throw; } } if (hmmec.solve(blupsol,rellrhs,method,relax,stopval,mxiter) == 1) { retval = &blupsol; //hmmec.close(); } // note that after BLUP, HMME remains intack // hmmec.resize(0,0); return retval; }

void matvec::Model::blup_pccg (  double    tol = 1.0e-8, unsigned int    mxiter = 100000 )

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 bad_model, matvec::Vector< double >::begin(), matvec::Vector< T >::begin(), blupsol, matvec::Vector< T >::clear(), matvec::Vector< double >::clear(), compute_rhs_diag_mme(), data_prepared, fixed_model, matvec::ginverse2(), hmmesize, matvec::doubleMatrix::identity(), matvec::Vector< T >::inner_product(), inverse_residual_var(), mixed_model, mme_times(), mme_times_res, nt_vec, matvec::Matrix< double >::num_rows(), numterm, numtrait, matvec::preconditioning(), prepare_data(), rellrhs, matvec::Vector< T >::reserve(), matvec::Vector< double >::reserve(), residual_var, matvec::Vector< double >::resize(), term, and type. Referenced by blup(). { if (type == bad_model) throw exception("Model::blup_pccg: bad model"); if (!data_prepared) { if (!prepare_data("iod")) { type = bad_model; throw exception("Model::blup_pccg: bad model"); } } if (type == mixed_model) { if (residual_var.num_rows() != numtrait) { type = bad_model; throw exception("Model::blup_pccg: residual variance has not yet assigned"); } } else if (type == fixed_model) { if (residual_var.num_rows() != numtrait) { residual_var.identity(numtrait,numtrait); } } else { throw exception("Model::setup_mme(): inappropriate model"); } inverse_residual_var(); blupsol.resize(hmmesize,0.0); rellrhs.reserve(hmmesize); unsigned i,j,k; unsigned n; unsigned nt = (numtrait+1)*numtrait/2; for (n=0,i=0; i<numterm; ++i) n += term[i].nlevel(); double **M; M = new(nothrow) double *[n]; for (i=0; i<n; ++i) { M[i] = new(nothrow) double[nt+1]; //M[i][0] is reserved to store the size memset(&(M[i][1]),'\0',sizeof(double)*nt); } k = 0; for (i=0; i<numterm; ++i) { nt = nt_vec[i]; for (j=0; j<term[i].nlevel(); ++j) M[k++][0] = nt; } compute_rhs_diag_mme(M); for (i=0; i<n; ++i) { // nt = (unsigned)M[i][0]; // nt = nt*(nt+1)/2; // for (j=1; j<=nt; ++j) cout << " " << M[i][j]; // cout << "\n"; ginverse2(&(M[i][1]),(unsigned)M[i][0],1.0e-15); } mme_times_res.reserve(hmmesize); blupsol.resize(hmmesize); Vector<double> r; r = rellrhs; Vector<double> d; d.reserve(hmmesize); ///////// d = M*r /////////////// preconditioning(d.begin(),(const double **)M,r.begin(),n); double delta_new, delta_old, convcr = 1.0; unsigned iter = 0; double alpha,beta; delta_new = r.inner_product(d); double delta_0 = delta_new*tol*tol; bool done = false; while (!done){ iter++; std::cout <<"iteration: " << iter <<std::endl; /////// q = Ad ///////// mme_times(d);// the result is stored in mme_times_res alpha = delta_new/(d.inner_product(mme_times_res)); blupsol += alpha*d; if(n%50){ //////// r = r - alpha*q ////////// r -= alpha*mme_times_res; } else{ ///// r = b - Ax ///////// mme_times(blupsol); r = rellrhs - mme_times_res; } ///////// q = M*r /////////////// preconditioning(mme_times_res.begin(),(const double **)M,r.begin(),n); delta_old = delta_new; delta_new = r.inner_product(mme_times_res); beta = delta_new/delta_old; d = mme_times_res + beta*d; if (delta_new > delta_0 && iter < maxiter){ done = false; } else if (delta_new <= delta_0){ ///// r = b - Ax ///////// mme_times(blupsol); r = rellrhs - mme_times_res; ///////// q = M*r /////////////// preconditioning(mme_times_res.begin(),(const double **)M,r.begin(),n); delta_new = r.inner_product(mme_times_res); if(delta_new > delta_0 && iter < maxiter){ done = false; } else { done = true; } } else { done = true; } } // std::cout <<"r = "<<r; mme_times_res.clear(); r.clear(); }

int matvec::Model::build_model_term (  const std::string &    modelspecs )  [protected]

Definition at line 2158 of file model.cpp. References base_effect, categorical, matvec::check_ptr(), matvec::FieldStruct::classi(), corr_link, factor_struct, matvec::fieldindex(), matvec::FieldStruct::name(), matvec::nest_xaction(), nt_vec, numfact, numterm, numtrait, 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(), term, trait_map, trait_struct, traitname, unknown_prior, and var_link. Referenced by equation(). { int retval = 1; unsigned i,j,k; std::string mspec = modelspec; // modelspec: y1 = A B + A (B), y2 = B + X*B nest_xaction(mspec); // now mspec: y1 = A B + A**B , y2 = B + X*B std::string sep("=+,* "); std::string TMP = mspec; //////// TMP.unique("=+,* "); // TMP ->y1 y2 A B X std::vector<std::string> tmpvec; std::vector<std::string>::iterator it; std::string::size_type begidx; split(TMP,sep,&tmpvec); sort(tmpvec.begin(),tmpvec.end()); it = unique(tmpvec.begin(),tmpvec.end()); tmpvec.erase(it,tmpvec.end()); // for(it = tmpvec.begin(); it != tmpvec.end(); it++){ // std::cout << *it << endl; // } categorical = 0; numfact = tmpvec.size(); if(factor_struct){ delete [] factor_struct; factor_struct = 0; } if(numfact>0){ factor_struct = new FieldStruct [numfact]; } else { factor_struct = 0; } check_ptr(factor_struct); for (i=0; i<numfact; i++) factor_struct[i].name(tmpvec[i]); TMP = mspec; numtrait = 0; sep = "="; begidx = TMP.find(sep,0); while(begidx != std::string::npos) { numtrait++; begidx = TMP.find(sep,begidx+1); } if (numtrait == 0) { std::cerr << "Model::build_model_term(): there is no trait in the model\n"; return 0; } if(trait_struct){ delete [] trait_struct; trait_struct = 0; } if(numtrait>0){ trait_struct = new FieldStruct* [numtrait]; } else { trait_struct = 0; } check_ptr(trait_struct); if(traitname){ delete [] traitname; traitname = 0; } if(numtrait>0){ traitname = new std::string [numtrait]; } else { traitname = 0; } check_ptr(traitname); std::vector<std::string> model_trait; std::vector<std::string> trtnm; sep = ","; split(TMP,sep,&model_trait); if (numtrait != model_trait.size()) throw exception("Model::build_model_term(): model for each trait is not separated properly"); sep = " +"; unsigned n = 0; for (i=0; i<numtrait; i++) { // trait name should be before "=" begidx = model_trait[i].find("="); // bbbtrtnmbb is the trait name with possible surrounding blanks std::string bbbtrtnmbbb = model_trait[i].substr(0,begidx); // getting rid of the surrounding blanks split(bbbtrtnmbbb,sep,&trtnm); traitname[i] = trtnm[i]; model_trait[i] = model_trait[i].substr(begidx+1); n += split(model_trait[i],sep); /// n += model_trait[i].ntoken("+ "); k = fieldindex(traitname[i],factor_struct,numfact); if (k >= 0 && k < numfact) { factor_struct[k].classi('T'); trait_struct[i] = &(factor_struct[k]); } else { throw exception("Model::build_model_term(): you have probably found a bug!"); } } TermList T(n,numtrait); tmpvec.clear(); sep = " +"; for (numterm=0,i=0; i<numtrait; i++) { tmpvec.clear(); n = split(model_trait[i],sep,&tmpvec); // tmpvec = model_trait[i].split(n,"+ "); for (k=0; k<n; k++) { for (j=0; j<numterm; j++) { if (T[j].match(tmpvec[k],"* ",factor_struct)) break; } if (j==numterm) { T[numterm++].input(tmpvec[k].c_str(),factor_struct,numfact); } T[j].trait[i]= 1; } } term.resize(numterm,numtrait); if (unknown_prior) { delete [] unknown_prior; // note that unknown_prior[i] unknown_prior = 0; } if(numterm>0){ unknown_prior = new UnknownDist[numterm]; // needs to resize(numtrait) } else { unknown_prior = 0; } check_ptr(unknown_prior); nt_vec.resize(numterm,numtrait); base_effect.reserve(numterm); corr_link.resize(numterm,-1); pos_vec.resize(numterm); trait_map.reserve(numterm); var_link.reserve(numterm); for (i=0; i<numterm; i++) { base_effect[i]=i; var_link[i]=i; term[i] = T[i]; term[i].prior = &(unknown_prior[i]); trait_map[i].resize(numtrait,i); } return retval; }

void matvec::Model::calc_gamete_prob (  void    )

void matvec::Model::compute_rhs_diag_mme (  double **    M )  [protected]

Compute rhs and diagonals of mme Definition at line 742 of file model_blup.cpp. References add_G_1_diag(), matvec::Vector< double >::begin(), hmmesize, input_pos_val_iod(), matvec::lidx(), nt_vec, numterm, numtrait, pos_val_vector, rellrhs, rve, and term. Referenced by blup_pccg(). { unsigned i,j,k,kk,ii,jj,t1,t2,tt,rec; std::vector<pos_val_node>::iterator it; double **vep,val,vy,xval; memset(rellrhs.begin(),'\0',sizeof(double)*hmmesize); double *rhstmp = rellrhs.begin()-1; Matrix<double> ve(numtrait,numtrait); for (it=pos_val_vector.begin(); it!=pos_val_vector.end(); it++) { if(!it->pos_term.empty()) { input_pos_val_iod(it); } else { continue; } vep = rve[it->pos_term[numterm]].begin(); val = it->xval_term[numterm]; // val = weight variable for (t1=0; t1<numtrait; t1++) { for (t2=0; t2<numtrait; t2++) ve[t1][t2] = val*vep[t1][t2]; } for (t1=0; t1<numtrait; t1++) { vy = 0.0; for(t2=0; t2<numtrait; t2++) vy += ve[t1][t2]*it->trait_vec[t2]; for (i=0; i<numterm; i++) { ii = term[i].addr[t1]; if (ii == 0) continue; k = (ii-1)/numtrait; xval = it->xval_term[i]; for (t2=t1; t2<numtrait; t2++) { jj = term[i].addr[t2]; if (jj) { val = xval* ve[t1][t2]; M[k][lidx(t2,t1,nt_vec[i])+1] += val*xval; } } rhstmp[ii] += vy*xval; } } } add_G_1_diag(M); }

void matvec::Model::compute_xbzu (  Vector< double > &    bu )

Definition at line 204 of file model_gibbs.cpp. Referenced by genotype_dist_gibbs(), genotype_dist_peeling(), genotypic_value_gibbs(), genotypic_value_peeling(), graph_log_likelihood_peeling(), and log_likelihood_peeling(). { unsigned i,ii,rec; double val; double *sol = bu.begin()-1; // for gcc-3.2 replaced // std::istrstream modelfile(modelstr, modelpcount); // with std::istringstream modelfile(modelstringstr); if (!modelfile) throw exception("Model::compute_xbzu(): cannot open file"); for (rec=0; rec<numrec; rec++) { input_pos_val(modelfile); val = 0.0; for (i=0; i<numterm; i++) { ii = term[i].addr[0]; // only one trait allowed now if (ii) val += sol[pos_term[i]]*xval_term[i]; } ii = rec_indid[rec]; if (ii > 0) pop->popmember[ii-1]->xbzu(val); } //BRS modelfile.close(); }

double matvec::Model::contrast (  const double **    kpme, const int    nr, const int    nc, const double *    M = 0, const int    prt_flag = 1, double **    result = 0 )  [virtual]

Definition at line 276 of file model_hypo.cpp. References matvec::warning(). { //////////////////////////////////////////////////// // H0: K'*b = M // trailing zeros in K' can be omitted // REMINDER: blupsol & rellrhs remain intact // // if M=NULL, then it acts as a vector of zeros. // meanp must be declared as double meanp[nr][4] ///////////////////////////////////////////////////// if (!kpme) { warning("Model::contrast(kp): kp is null, thus it is ignored"); return 0.0; } if (!get_blupsol()) throw exception("Model::contrast(): bad model"); if (nc > hmmesize) throw exception("Model::contrast(): size not conformable"); int est=1; doubleMatrix kvk(nr,nr); doubleMatrix kvkori(nr,nr); Vector<double> kpdiff(nr); Vector<double> kpb_m(nr); Vector<double> xy(hmmesize); Vector<double> tmpsol(hmmesize); double *sol, kpd; const double *kp; unsigned i,j,k; for (i=0; i<nr; i++) { if (nc < hmmesize) { memset(xy.begin(),'\0',sizeof(double)*hmmesize); memcpy(xy.begin(),kpme[i],sizeof(double)*nc); if (!hmmec.solve(tmpsol,xy,"ysmp")) return 0.0; } else { if (!hmmec.solve(tmpsol.begin(),kpme[i],"ysmp")) return 0.0; } hmmec.mv(tmpsol,xy); kp = kpme[i]; sol = xy.begin(); for (kpd=0.0,j=0; j<nc; j++) { kpd = std::max(kpd,fabs(*kp - *sol)); kp++; sol++; } for (j=nc; j<hmmesize; j++) { kpd = std::max(kpd,fabs(*sol)); sol++; } kpdiff[i] = kpd; for (k=0; k<nr; k++) { kp = kpme[k]; sol = tmpsol.begin(); for (kpd=0.0,j=0; j<nc; j++) kpd += *kp++ * *sol++; kvk[k][i] = kpd; } } kvkori = kvk; kvk.ginv1(&k); if (k < nr) throw exception("Model::contrast(K'): K'V{-1}K is singular. Possible reason:\n" " (1) K is non-estimable, and (2) K isn't of full column rank"); // K'b-M for (i=0; i<nr; i++) kpb_m[i] = blupsol.inner_product(kpme[i]); if (M) for (i=0; i<nr; i++) kpb_m[i] -= M[i]; double quq = kvk.quadratic(kpb_m,kpb_m); // (K'b-m)'(K'(X'V-1X)-K)-1(K'b-m) double ssm,sse; unsigned rank_x; int dfe = 0; anova(ssm,sse,rank_x,dfe); double res_var_est = sse/dfe; int degen = 0; int errcode = 0; double f_stat=0.0, prob=0.0; double sigma_e = residual_var[0][0]; int chisq = 0; if (chisq) { f_stat = quq; prob = ChiSquare_cdf(f_stat,static_cast<double>(nr),0.0,errcode); } else { if (dfe <= 0 || sse <= SESSION.epsilon) { warning("Model::contrast(): either dfe or sse is zero"); degen = 1; } else { f_stat = (quq*sigma_e/res_var_est)/nr; prob = F_cdf(f_stat,static_cast<double>(nr),static_cast<double>(dfe),0.0,errcode); } // in univariate cases only: the estimated variance of estimator K'B if (!degen) kvkori *= res_var_est/sigma_e; } std::string raw_data_code; double p_value,var,tcal=0.0; int W = SESSION.output_precision; std::cout.precision(W); if (prt_flag) { std::cout << "\n RESULTS FROM CONTRAST(S)\n"; std::cout << " ----------------------------------------------------------\n"; std::cout << " Contrast MME_addr K_coef Raw_data_code\n"; std::cout << " ---------------------------------------------\n"; } for (i=0; i<nr; i++) { if (prt_flag) { kp = kpme[i]; for (j=0; j<nc; j++) { if (kp[j] == 0.0) continue; trait_effect_level(j,raw_data_code); std::cout << std::setw(4) << i+1 <<std::setw(11) << j+1 << std::setw(13) << kp[j] << " " << raw_data_code << "\n"; } } kpd = kpdiff[i]; if (result) { result[i][0] = kpb_m[i]; result[i][3] = kpd; } var = kvkori[i][i]; if (var <= 0.0) throw exception(" Model::contrast(): you have probably found a bug!"); var = std::sqrt(var); tcal = fabs(kpb_m[i]/var); p_value = 2.0*(1.0-t_cdf(tcal,static_cast<double>(dfe),0.0)); if (result) { result[i][1] = var; result[i][2] = p_value; } if (!prt_flag) continue; if (kpd > 1.0e-5) { est = 0; std::cout << " **** NON-ESTIMABLE ****\n\n"; } else { if (kpd > 1.0e-10) { std::cout << " ***WARNING***: it may or mayn't be estimable\n"; std::cout << " max(k'*ginv(mme)*mme - k') = " << kpd << "\n"; } std::cout << " estimated value (K'b-M) = " << kpb_m[i] << " +- " << var << "\n"; std::cout << " Prob(|t| > " << tcal << ") = " << p_value << " (p_value)\n"; if (i != nr-1) std::cout << "\n"; } } if (prt_flag) { std::cout << " ----------------------------------------------------------\n"; if (nr > 1 && est) { std::cout << " joint hypothesis test H: K'b = M\n"; std::cout << " Prob(F > " << f_stat << ") = " << 1.0-prob << " (p_value)\n"; } } return 1.0-prob; // p_value }

double matvec::Model::contrast (  const std::string &    termname, const doubleMatrix &    Kp, const Vector< double > &    M )  [virtual]

Reimplemented in matvec::GLMM. Definition at line 227 of file model_hypo.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< T >::begin(), matvec::Matrix< double >::begin(), contrast(), factor_struct, get_blupsol(), hmmesize, matvec::TermList::index(), matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), numtrait, matvec::Vector< T >::size(), term, and matvec::warning(). { if (!(Kp.begin())) { warning("Model::contrast(kp): kp is null, thus it is ignored"); return 0.0; } double retval = 0.0; if (!get_blupsol()) throw exception("Model::contrast(): bad model"); int k = term.index(termname,factor_struct); if (k < 0) throw exception("Model::contrast(): no such term in the model"); unsigned nr = Kp.num_rows(); unsigned nc = Kp.num_cols(); if (M.size() != nr) throw exception("Model::contrast(): unconformable size"); unsigned i,startaddr; if (nc <= term[k].nlevel()*numtrait) { startaddr = term[k].startaddr(); Matrix<double> fullsize_Kp(nr,hmmesize); for (i=0; i<nr; i++) { memcpy(&(fullsize_Kp[i][startaddr]),Kp[i],sizeof(double)*nc); } contrast((const double **)fullsize_Kp.begin(),nr,hmmesize,(const double*)M.begin()); } else { throw exception("Model::contrast(): too many columns in Kp"); } return retval; }

double matvec::Model::contrast (  const doubleMatrix &    Kp )  [virtual]

Reimplemented in matvec::GLMM. Definition at line 207 of file model_hypo.cpp. References matvec::Matrix< double >::begin(), contrast(), matvec::Matrix< double >::num_cols(), and matvec::Matrix< double >::num_rows(). { /////////////////////////////////////////// // trailing zeros in K' can be omitted /////////////////////////////////////////// return contrast((const double **)Kp.begin(), Kp.num_rows(), Kp.num_cols()); }

double matvec::Model::contrast (  const doubleMatrix &    Kp, const Vector< double > &    M )  [virtual]

Reimplemented in matvec::GLMM. Definition at line 216 of file model_hypo.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< double >::begin(), contrast(), matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), and matvec::Vector< T >::size(). { /////////////////////////////////////////// // trailing zeros in K' can be omitted /////////////////////////////////////////// int nr = Kp.num_rows(); if (M.size() != nr ) throw exception(" Model::contrast(): size is incomformable"); return contrast((const double **)Kp.begin(),nr,Kp.num_cols(),M.begin()); }

double matvec::Model::contrast (  const Vector< double > &    Kp, const double    m = 0.0 )  [virtual]

Reimplemented in matvec::GLMM. Definition at line 181 of file model_hypo.cpp. References matvec::Vector< T >::begin(), and matvec::Vector< T >::size(). Referenced by contrast(), and out_lsmeans_to_stream(). { /////////////////////////////////////////// // trailing zeros in K' can be omitted /////////////////////////////////////////// double retval, **kpme = new double *[1]; kpme[0] = Kp.begin(); unsigned nc = Kp.size(); double* M = 0; if (m != 0.0) { M = new double [1]; M[0] = m; } retval = contrast((const double **)kpme,1,nc,M); if(kpme){ delete [] kpme; kpme = 0; } if (M) { delete [] M; M = 0; } return retval; }

void matvec::Model::copyfrom (  const Model &    M )

Definition at line 83 of file model.cpp. References blupsol, matvec::check_ptr(), matvec::UnknownDist::copyfrom(), data, data_prepared, dfreml_called, factor_struct, matvec::fieldindex(), matvec::Data::hashtable, hmmec, hmmesize, idlist, initialize(), kvec, lng0vec, lnr0vec, max_nz, maxorder, modelfname, modelstring, matvec::GeneticDist::name(), nfunk_in_dfreml, non_zero, npattern, ntermGdist, num_ped, numcol, numfact, numgroup, numobs, numrec, numterm, numtrait, pattern, pattern_tree, pop, pop_created, popsize, pos_term, matvec::ModelTerm::prior, rawpos, rec_indid, rellrhs, reml_value, matvec::Matrix< double >::reserve(), residual_var, matvec::SparseMatrix::resize(), rve, term, matvec::TermList::termlist, trait_struct, trait_vec, traitname, type, unknown_prior, weightinverse, weightname, xact_htable, xval_term, and yry. Referenced by Model(). { initialize(); num_ped = M.num_ped; yry = M.yry; data = M.data; term = M.term; max_nz = M.max_nz; numrec = M.numrec; numobs = M.numobs; numcol = M.numcol; max_nz = M.max_nz; popsize = M.popsize; numterm = M.numterm; numfact = M.numfact; lng0vec = M.lng0vec; lnr0vec = M.lnr0vec; blupsol = M.blupsol; rellrhs = M.rellrhs; popsize = M.popsize; hmmesize = M.hmmesize; numtrait = M.numtrait; npattern = M.npattern; hmmesize = M.hmmesize; numgroup = M.numgroup; maxorder = M.maxorder; non_zero = M.non_zero; rec_indid = M.rec_indid; type = M.type; modelfname = M.modelfname; weightname = M.weightname; ntermGdist = M.ntermGdist; modelstring = M.modelstring; pattern_tree = M.pattern_tree; residual_var = M.residual_var; data_prepared = M.data_prepared; weightinverse = M.weightinverse; nfunk_in_dfreml= M.nfunk_in_dfreml; reml_value = M.reml_value; dfreml_called = M.dfreml_called; hmmec.resize(hmmesize,max_nz); int i,k; // pos_term, xval_term, trait_vec are working spaces if (pos_term) {delete [] pos_term; pos_term = 0;} if (M.pos_term) pos_term = new unsigned [numterm+1]; if (xval_term) {delete [] xval_term; xval_term = 0;} if (M.xval_term) xval_term = new double [numterm+1]; if (trait_vec) {delete [] trait_vec; trait_vec= 0;} if (M.trait_vec) trait_vec = new double [numterm+1]; if (unknown_prior) { delete [] unknown_prior;unknown_prior = 0; } if (M.unknown_prior) { if(numterm>0){ unknown_prior = new UnknownDist[numterm]; } else { unknown_prior = 0; } for (i=0; i<numterm; i++) { unknown_prior[i].copyfrom(M.unknown_prior[i]); if (strcmp(M.term.termlist[i].prior->name(),"UnknownDist") == 0) { term.termlist[i].prior = &(unknown_prior[i]); } else { term.termlist[i].prior = M.term.termlist[i].prior; } } } if (traitname) {delete [] traitname; traitname = 0;} if (M.traitname) { if(numtrait>0){ traitname = new std::string [numtrait]; } else { traitname = 0; } for (i=0; i<numtrait; i++) traitname[i] = M.traitname[i]; } if (kvec) {delete [] kvec; kvec = 0;} if (M.kvec) { if(numtrait>0){ kvec = new unsigned [npattern]; } else { kvec = 0; } for (i=0; i<npattern; i++) kvec[i] = M.kvec[i]; } if (rawpos) {delete [] rawpos; rawpos = 0;} if (M.rawpos) { if(maxorder>0){ rawpos = new unsigned [maxorder]; // this is just a working space } else { rawpos = 0; } } pattern.resize(npattern); std::set<std::string>::iterator pos; for (i=0, pos=pattern_tree.begin(); pos != pattern_tree.end(); ++pos,++i) { pattern[i] = *pos; } pop_created = M.pop_created; if (M.pop_created) { pop = new Population(*(M.pop)); } else { pop = M.pop; } if (rve) { delete [] rve; rve = 0; } if (M.rve) { if(npattern>0){ rve = new Matrix<double> [npattern]; } else { rve = 0; } check_ptr(rve); for (i=0; i<npattern; i++) rve[i].reserve(numtrait,numtrait); } if (factor_struct) { delete [] factor_struct; factor_struct = 0; } if (M.factor_struct) { if(numfact>0){ factor_struct = new FieldStruct [numfact]; } else { factor_struct = 0; } for (i=0; i<numfact; i++) factor_struct[i] = M.factor_struct[i]; } if (trait_struct) { delete [] trait_struct; trait_struct = 0; } if (M.trait_struct) { if(numtrait>0){ trait_struct = new FieldStruct* [numtrait]; } else { trait_struct = 0; } for (i=0; i<numtrait; i++) { k = fieldindex(traitname[i],factor_struct,numfact); if (k >= 0 && k < numfact) { trait_struct[i] = &(factor_struct[k]); } else { throw exception("Model::copyfrom(): you have probably found a bug!"); } } } if (xact_htable) {delete [] xact_htable; xact_htable = 0;} if (M.xact_htable) { if(numterm>0){ xact_htable = new HashTable [numterm]; } else { xact_htable = 0; } for (i=0; i<numterm; i++) xact_htable[i] = M.xact_htable[i]; } if (idlist) {delete [] idlist; idlist = 0;} if (M.idlist) { if(numcol>0){ idlist = new HashTable* [numcol]; } else { idlist = 0; } check_ptr(idlist); if (data) for (i=0; i<numcol; i++) idlist[i] = data->hashtable[i]; } }

void matvec::Model::covariance (  const std::string &    termname, const std::string &    termname2, const double    v00, ...    )

Definition at line 2059 of file model.cpp. : // * // * M.variance("animal",1.0,...) // ********************************************************************* { if (numtrait<1) { std::cerr << "Model::variance(args): no trait(dependent-variable) in the model\n"; return; } doubleMatrix var(numtrait,numtrait); int t1,t2; va_list param_pt; // an object param_pt va_start(param_pt,v00); // call the setup macro var[0][0] = v00; for (t2=1; t2<numtrait; t2++) var[0][t2] = va_arg(param_pt,double); for (t1=1; t1<numtrait; t1++) for (t2=0; t2<numtrait; t2++) { var[t1][t2] = va_arg(param_pt,double); } va_end(param_pt); covariance(termname,termname2,var); }

void matvec::Model::covariance (  const std::string &    termname, const std::string &    termname2, const doubleMatrix &    v )

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(). { if (!factor_struct) { std::cerr << "Model::covariance(args): no equation(s) in the model\n"; return; } if (v.num_rows() != v.num_cols() || v.num_rows() != numtrait) { std::cerr << "Model::covariance(): invalid variance matrix size\n"; return; } int i,k,k2; k = term.index(termname,factor_struct); k2 = term.index(termname2,factor_struct); if( k < 0 ){ std::cerr << "Model::covariance(..): not in model: " << termname << "\n"; return; } if( k2 < 0 ){ std::cerr << "Model::covariance(..): not in model: " << termname2 << "\n"; return; } int bek=base_effect[k]; if(term[bek].classi() != 'R' && term[bek].classi() != 'P'){ std::cerr << "Model::covariance(..): not a random effect: " << termname << "\n"; return; } int bek2=base_effect[k2]; if(term[bek2].classi() != 'R' && term[bek2].classi() != 'P'){ std::cerr << "Model::covariance(..): not a random effect: " << termname2 << "\n"; return; } if(term[bek].classi() != term[bek2].classi() ){ std::cerr << "Model::covariance(..): different types of random effects: \n"; return; } int be; if(base_effect[k]==base_effect[k2]){ be=base_effect[k]; } else { be=base_effect[k]; int be2=base_effect[k2]; if(be2 <be) { int tmp; tmp=k; k=k2; k2=tmp; tmp=be; be=be2; be2=tmp; } base_effect[k2]=be; doubleMatrix orig1,orig2; orig1.copy(*term[be].prior->var_matrix()); orig2.copy(*term[be2].prior->var_matrix()); int newsize=nt_vec[be]+nt_vec[be2]; Vector<int> orig=trait_map[be]; trait_map[be].reserve(newsize); int ii=0; for(int i=0;i<nt_vec[be];i++,ii++) trait_map[be][ii]=orig[i]; for(int i=0;i<nt_vec[be2];i++,ii++) trait_map[be][ii]=trait_map[be2][i]; term[be].prior->var_matrix()->resize(newsize,newsize); for(int i=0;i<nt_vec[be];i++){ for(int j=0;j<nt_vec[be];j++){ term[be].prior->var_matrix()->me[i][j]=orig1[i][j]; } } for(int i=0;i<nt_vec[be2];i++){ for(int j=0;j<nt_vec[be2];j++){ term[be].prior->var_matrix()->me[nt_vec[be]+i][nt_vec[be]+j]=orig2[i][j]; } } nt_vec[be]=newsize; term[be2].prior->var_matrix()->clear(); nt_vec[be2]=0; term[be2].classi('F'); int last_eff=be; for(;corr_link[last_eff] != -1;)last_eff=corr_link[last_eff]; corr_link[last_eff]=be2; int last_pos=pos_vec[last_eff]; for(;corr_link[last_eff] != -1;){ last_eff=corr_link[last_eff]; pos_vec[last_eff]=pos_vec[last_eff]+last_pos+1; } } int srow,scol; srow=pos_vec[k]*numtrait; scol=pos_vec[k2]*numtrait; for(int irow=0;irow<numtrait;irow++){ for(int jcol=0;jcol<numtrait;jcol++){ term[be].prior->var_matrix()->me[srow+irow][scol+jcol]=v[irow][jcol]; term[be].prior->var_matrix()->me[scol+jcol][srow+irow]=v[irow][jcol]; } } // std::cout << "Cov " << *term[be].prior->var_matrix(); }

void matvec::Model::covariate (  const std::string &    factorname )

Definition at line 2097 of file model.cpp. References matvec::FieldStruct::classi(), data, factor_struct, matvec::fieldindex(), numfact, matvec::split(), and term. { // ********************************************************************* // it specifies a list of effect (factor) as covariates. // note that an effect is not a model-term, which consists of effects, // instead an effect is associated with a data column. It is required // that the name of an effect must be the same as the associated // data-column. // ********************************************************************** if (!factor_struct) { std::cerr << "Model::covariate(): no equation(s) in the model\n"; return; } int j,k,t; unsigned i,nc; std::vector<std::string> efflist; nc = split(covariate_names," ,",&efflist); // efflist = covariate_names.split(nc,", "); for (i=0; i<nc; i++) { k = fieldindex(efflist[i],factor_struct,numfact); if (k < 0) throw exception("Model::covariate(): not in the model, thus it's ignored"); // need to check if classi is changed if (factor_struct[k].classi() != 'F') data= 0; factor_struct[k].classi('C'); //////////////////////////////////////////////////////////////////// // if a term contains at least one covariate, then its classi = 'C' // else leave it alone ////////////////////////////////////////////////////////////////// for (t=0; t<numterm; t++) { j = term[t].partial_match(efflist[i].c_str(),factor_struct); if (j >=0 ) { term[t].classi('C'); break; } } } }

doubleMatrix matvec::Model::CovMat (  const std::string &    termname, doubleMatrix &    Kp )

Definition at line 128 of file model_hypo.cpp. References CovMat(), factor_struct, get_blupsol(), hmmesize, matvec::TermList::index(), matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), numtrait, and term. { doubleMatrix retval; if (!get_blupsol()) throw exception("Model::CovMat(): bad model"); int k = term.index(termname,factor_struct); if (k < 0) throw exception("Model::CovMat(): no such term in the model"); unsigned nr = Kp.num_rows(); unsigned nc = Kp.num_cols(); unsigned i,startaddr; if (nc <= term[k].nlevel()*numtrait) { startaddr = term[k].startaddr(); doubleMatrix fullsize_Kp(nr,hmmesize); for (i=0; i<nr; i++) { memcpy(&(fullsize_Kp[i][startaddr]),Kp[i], sizeof(double)*nc); } retval=CovMat(fullsize_Kp); } else { throw exception("Model::CovMat(): too many columns in Kp"); } return retval; }

doubleMatrix matvec::Model::CovMat (  doubleMatrix &    Kp )

Definition at line 69 of file model_hypo.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< double >::begin(), get_blupsol(), hmmec, hmmesize, matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), and matvec::SparseMatrix::solve(). Referenced by CovMat(). { unsigned nr=Kp.num_rows(); unsigned nc = Kp.num_cols(); double *kp,*sol; doubleMatrix KvK(nr,nr); if (!get_blupsol("ysmp")) throw exception("Model::CovMat(): bad model"); if (nc > hmmesize) throw exception("Model::CovMat(): size not conformable"); Vector<double> xy(hmmesize); Vector<double> tmpsol(hmmesize); double kpd,**kpme=Kp.begin(); unsigned j; for (unsigned i=0; i<nr; i++) { if (nc < hmmesize) { memset(xy.begin(),'\0',sizeof(double)*hmmesize); memcpy(xy.begin(),kpme[i],sizeof(double)*nc); if (!hmmec.solve(tmpsol,xy,"ysmp")) throw exception("Model::CovMat(): hmmec.solve failed"); } else { if (!hmmec.solve(tmpsol.begin(),kpme[i],"ysmp")) throw exception("Model::CovMat(): hmmec.solve failed"); } for (unsigned k=0; k<nr; k++) { kp = Kp[k]; sol = tmpsol.begin(); for (kpd=0.0,j=0; j<nc; j++) kpd += *kp++ * *sol++; KvK[k][i] = kpd; } } return(KvK); }

void matvec::Model::descent_graph_peeling_initialisation (  void    )

void matvec::Model::DGSampler (  unsigned    numOfSamples, unsigned    numOfSL, unsigned    numOfHaplo, unsigned    numOfSLCas )

Definition at line 1176 of file model.cpp. References matvec::Individual::display_freq_haplotype(), matvec::Population::member(), matvec::Population::MH_ibd_sample_map(), pop, matvec::Population::size(), and matvec::Population::sum_descentState_map(). { // Authors: L. Radu Totir // (August, 2004) // Contributors: double l_hood = 0.0; unsigned k = 0, j = 0, option = 0; unsigned k1 = numOfSL; unsigned k2 = numOfHaplo; unsigned k3 = numOfSLCas; cout<< endl << "Descent graph iterations = " << numOfSamples << endl; for(int sample=1; sample <= numOfSamples; sample++) { if(sample%5000 == 0) cout<<sample<<"."; // cout<<endl<< "SAMPLE.................................. " << sample <<endl; if (k >= k1 + k2 + k3) {k = 0;} if (k < k1) {option = 1;} else if ((k >= k1) && (k < k1 + k2)) {option = 2;} else if (k <= k1 + k2 + k3) {option = 3;} k++; j++; l_hood = pop->MH_ibd_sample_map(l_hood,option); pop->sum_descentState_map(); // update pdq's for all marker loci if (j >= 100000) { j = 0; cout << "SI_PDQ sample = " << sample << endl; } } cout << endl << "Haplotypes origins mm mp pm pp / mother - father " << endl; for (int i=0; i<pop->size(); i++) { pop->member(i)->display_freq_haplotype(numOfSamples); } cout << endl; }

void matvec::Model::DGSamplerSetup (  unsigned    numLoci, Data *    D )

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, pop, popsize, matvec::Vector< unsigned >::resize(), matvec::Individual::search_heteroz(), matvec::Individual::set_founder_alleles(), and matvec::Population::size(). { // Authors: L. Radu Totir // (August, 2004) // Contributors: pop->descent_graph_setup(D); pop->founder_allele_counter = 0; int popsize = pop->size(); int zero = 0; for (int i=0;i<popsize;i++) { pop->member(i)->set_founder_alleles(0); } unsigned lastLocus = numLoci; for (int i=0; i<popsize; i++) { if (pop->member(i)->m_counter_map.empty()){ int zero = 0; pop->member(i)->m_counter_map.resize(numLoci,zero); pop->member(i)->p_counter_map.resize(numLoci,zero); } pop->member(i)->search_heteroz(lastLocus); //cout << pop->member(i)-> id() << " " << pop->member(i)-> ord_heter << endl; } }

int matvec::Model::display (  void    )  const

Definition at line 324 of file model.cpp. References modelstring. { std::cout << " model: " << modelstring << "\n"; return 1; }

int matvec::Model::equation (  const std::string &    modelspecs )

Reimplemented in matvec::GLMM. Definition at line 264 of file model.cpp. References bad_model, build_model_term(), fixed_model, lng0vec, matvec::TermList::maxorder(), maxorder, modelstring, num_ped, numterm, rawpos, matvec::Vector< double >::resize(), term, and type. Referenced by Model(). { //****************************************************************** // take the folowing model as example to show how I process model // // model M("y1 = A + B + A * B," // "y2 = B X(B)" // // ***************************************************************** int retval = 1; num_ped=0; if (modelspecs == "") throw exception("Model::equation(): equation is empty"); modelstring = modelspecs; if (build_model_term(modelstring)) { maxorder = term.maxorder(); if (rawpos) { delete [] rawpos; rawpos=0; } if(maxorder>0){ rawpos = new unsigned [maxorder]; // this is just a working space } else { rawpos = 0; } lng0vec.resize(numterm); } else { modelstring = ""; retval = 0; } if (retval == 1) { type = fixed_model; } else { type = bad_model; } return retval; }

void matvec::Model::estimate (  const double **    kpme, const unsigned    nr, const unsigned    nc, double *    kpb )

Definition at line 154 of file model_hypo.cpp. References matvec::warning(). { /////////////////////////////////////////// // trailing zeros in K' can be omitted /////////////////////////////////////////// if (!kpme) { warning("Model::estimate(kp): kp is null, thus it is ignored"); return; } if (!get_blupsol("ysmp")) throw exception("Model::estimate(): bad model"); Vector<double> xy; if (nc < hmmesize) xy.resize(hmmesize); for (unsigned i=0; i<nr; i++) { if (nc < hmmesize) { memcpy(xy.begin(),kpme[i],sizeof(double)*nc); kpb[i] = blupsol.inner_product(xy); } else { kpb[i] = blupsol.inner_product(kpme[i]); } } }

Vector< double > matvec::Model::estimate (  const std::string &    termname, const doubleMatrix &    Kp )

Definition at line 99 of file model_hypo.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< T >::begin(), estimate(), factor_struct, get_blupsol(), hmmesize, matvec::TermList::index(), matvec::Matrix< double >::num_cols(), matvec::Matrix< double >::num_rows(), numtrait, matvec::Vector< T >::resize(), and term. { Vector<double> retval; if (!get_blupsol()) throw exception("Model::estimate(): bad model"); int k = term.index(termname,factor_struct); if (k < 0) throw exception("Model::estimate(): no such term in the model"); unsigned nr = Kp.num_rows(); unsigned nc = Kp.num_cols(); unsigned i,startaddr; if (nc <= term[k].nlevel()*numtrait) { startaddr = term[k].startaddr(); Matrix<double> fullsize_Kp(nr,hmmesize); for (i=0; i<nr; i++) { memcpy(&(fullsize_Kp[i][startaddr]),Kp[i], sizeof(double)*nc); } retval.resize(nr); estimate((const double **)fullsize_Kp.begin(),nr,hmmesize,retval.begin()); } else { throw exception("Model::estimate(): too many columns in Kp"); } return retval; }

Vector< double > matvec::Model::estimate (  const doubleMatrix &    Kp )

Definition at line 55 of file model_hypo.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< double >::begin(), estimate(), matvec::Matrix< double >::num_cols(), and matvec::Matrix< double >::num_rows(). { /////////////////////////////////////////// // trailing zeros in K' can be omitted /////////////////////////////////////////// unsigned nr = Kp.num_rows(); unsigned nc = Kp.num_cols(); Vector<double> kpb(nr); estimate((const double **)Kp.begin(),nr,nc,kpb.begin()); return kpb; }

double matvec::Model::estimate (  const Vector< double > &    Kp )

Definition at line 37 of file model_hypo.cpp. References matvec::Vector< T >::begin(), and matvec::Vector< T >::size(). Referenced by estimate(). { /////////////////////////////////////////// // trailing zeros in K' can be omitted /////////////////////////////////////////// double **kpme = new double *[1]; kpme[0] = Kp.begin(); unsigned nc = Kp.size(); double kpb[1]; estimate((const double **)kpme,1,nc,kpb); if(kpme) { delete [] kpme; kpme = 0; } return kpb[0]; }

double matvec::Model::estimate_peeling (  void    )

Definition at line 113 of file model_peeling.cpp. References matvec::ChromStruct::locus. { if (type == bad_model) throw exception("Model::estimate_peeling(): bad model"); double max_llhood = 0.0; if (!data_prepared) { if (!prepare_data()) { type = bad_model; return max_llhood; } } minfun_indx = 2; pop->residual_var = residual_var; pop->input_data(data); pop->build_pop_gamete(); unsigned i,n,m,k,j; m = pop->tn_gamete - 1; // because p(A1)+ p(A2) +...+ p(At) = 1 n = 1 + m; // n = residual_var + allele_freq // n += pop->tn_genotype; // + genotypic_value // uncommented in old version Vector<double> x(n); x[0] = pop->residual_var[0][0]; // residual_variance ChromStruct *chrom = pop->prior->chrom(); double* ve = chrom[0].locus[0].allele_freq.begin(); for (i=0; i<m; i++) x[i+1] = ve[i]; // 1 is for residual_variance // uncommented in old version // k = 1 + m; // const double** me = pop->prior->genotypic_val(0,0); // for (i=0; i<tn_genotype; i++) for (j=0; j<=i; j++) x[k++] = me[i][j]; int iter = 0; max_llhood = praxis(x,n,iter); minfun_indx = 0; return max_llhood; }

void matvec::Model::fitdata (  Data &    D )

Definition at line 304 of file model.cpp. References bad_model, data, data_prepared, hmmec, hmmesize, max_nz, npattern, matvec::Data::num_rows(), numcol, numobs, numrec, matvec::SparseMatrix::resize(), type, and matvec::warning(). Referenced by matvec::Population::setupRSampler(). { if (D.num_rows() > 0) { data = &D; } else { data = 0; type = bad_model; warning("Model::fitdata(D): Data is empty"); } max_nz = 0; numrec = 0; numobs = 0; numcol = 0; npattern = 0; hmmesize = 0; data_prepared = 0; hmmec.resize(0,0); }

double matvec::Model::genotype_dist_gibbs (  const unsigned    warmup = 30, const unsigned    gibbslen = 1000 )

Definition at line 109 of file model_gibbs.cpp. References matvec::Vector< T >::begin(), matvec::Vector< double >::begin(), blupsol, compute_xbzu(), matvec::Population::cond_genotype_config(), data, data_prepared, matvec::Population::gibbs_iterate(), matvec::SparseMatrix::gibbs_iterate(), hmmec, hmmesize, matvec::Population::input_data(), ntermGdist, pop, prepare_data(), rellrhs, matvec::Vector< T >::reserve(), residual_var, matvec::Population::residual_var, matvec::Vector< double >::resize(), matvec::Population::resize_genotype_counter(), setup_mme(), and matvec::warning(). { if (ntermGdist==0) { warning(" Model::genotype_dist_gibbs(): no GeneticDist term"); return 0.0; } if (!data_prepared) prepare_data(); unsigned i,it; blupsol.resize(hmmesize,0.0); Vector<double> wspace; wspace.reserve(hmmesize); pop->resize_genotype_counter(); pop->residual_var = residual_var; pop->input_data(data); pop->cond_genotype_config(); // an initial configuration if (numterm) setup_mme(&rellrhs); ///////////////////////////////////////////////////// // Gibbs Sampler warm up ///////////////////////////////////////////// for (it=0; it<warmup; it++) { if (numterm) hmmec.gibbs_iterate(blupsol.begin(),rellrhs.begin(),wspace.begin()); for (i=0; i<ntermGdist; i++) { compute_xbzu(blupsol); // to adjust trait pop->gibbs_iterate(); } } ///////////////////////////////////////////////////// // samples from Gibbs Sampler now will be used //////////////////////////////////////////////////// for (it=0; it<gibbslen; it++) { if (numterm) hmmec.gibbs_iterate(blupsol.begin(),rellrhs.begin(),wspace.begin()); for (i=0; i<ntermGdist; i++) { compute_xbzu(blupsol); // to adjust trait pop->gibbs_iterate(1); } } hmmec.resize(0,0); ntermGdist = 0; return 0.0; // no useful value could be returned so far }

void matvec::Model::genotype_dist_peeling (  const int    prtlevel = 1, const int    estifreq = 1 )

Definition at line 29 of file model_peeling.cpp. References bad_model, blupsol, compute_xbzu(), data, data_prepared, matvec::Population::genotype_dist_peeling(), hmmec, hmmesize, matvec::Population::input_data(), ntermGdist, pop, prepare_data(), rellrhs, residual_var, matvec::Population::residual_var, matvec::SparseMatrix::resize(), matvec::Vector< double >::resize(), setup_mme(), matvec::SparseMatrix::solve(), and type. { ////////////////////////////////////////////////////////////////////////// // this works for arbitrary pedigrees // based on iterative peeling algorithm: Arendonk(1989), Fernando(1993) // prtl = 0, print nothing // 1, default, print genotype distribution and genotypic values // estifreq = 0, not estimate allele frequencies // 1, default, estimate allele frequencies ////////////////////////////////////////////////////////////////////////// if (type == bad_model) return; if (!data_prepared) { if (!prepare_data()) { type = bad_model; return; } } pop->residual_var = residual_var; pop->input_data(data); blupsol.resize(hmmesize); if (numterm) { setup_mme(&rellrhs); hmmec.solve(blupsol,rellrhs); } if (ntermGdist) { compute_xbzu(blupsol); // to adjust trait pop->genotype_dist_peeling(prtl,estifreq); } hmmec.resize(0,0); ntermGdist = 0; }

Vector< double > * matvec::Model::genotypic_value_cat (  const unsigned    warmup, const unsigned    gibbslen )

Definition at line 77 of file model_gibbs.cpp. References matvec::Vector< double >::begin(), matvec::Vector< T >::begin(), blupsol, data_prepared, matvec::SparseMatrix::gibbs_iterate(), hmmec, hmmesize, prepare_data(), rellrhs, matvec::Vector< T >::reserve(), matvec::Vector< double >::resize(), sampleW(), and setup_mme(). Referenced by genotypic_value_gibbs(). { if (!data_prepared) prepare_data(); blupsol.resize(hmmesize,0.0); Vector<double> tmpsol(hmmesize); Vector<double> wspace; wspace.reserve(hmmesize); setup_mme(&rellrhs); // relrhs is used here ///////////////////////////////////////////////////// // Gibbs Sampler warm up ///////////////////////////////////////////// unsigned it; for (it=0; it<warmup; it++) { sampleW(tmpsol,rellrhs); hmmec.gibbs_iterate(tmpsol.begin(),rellrhs.begin(),wspace.begin()); } ///////////////////////////////////////////////////// // samples from Gibbs Sampler now will be used //////////////////////////////////////////////////// for (it=0; it<gibbslen; it++) { sampleW(tmpsol,rellrhs); hmmec.gibbs_iterate(tmpsol.begin(),rellrhs.begin(),wspace.begin()); blupsol += tmpsol; } blupsol /= gibbslen; hmmec.resize(0,0); return &blupsol; }

Vector< double > * matvec::Model::genotypic_value_gibbs (  const unsigned    warmup = 30, const unsigned    gibbslen = 1000 )

Definition at line 153 of file model_gibbs.cpp. References matvec::Vector< T >::assign(), matvec::Vector< double >::begin(), matvec::Vector< T >::begin(), blupsol, compute_xbzu(), matvec::Population::cond_genotype_config(), data, data_prepared, genotypic_value_cat(), matvec::Population::gibbs_iterate(), matvec::SparseMatrix::gibbs_iterate(), hmmec, hmmesize, matvec::Population::input_data(), ntermGdist, pop, popsize, prepare_data(), rellrhs, matvec::Vector< T >::reserve(), residual_var, matvec::Population::residual_var, matvec::Vector< double >::resize(), and setup_mme(). { if (categorical) return genotypic_value_cat(warmup,gibbslen); if (!data_prepared) prepare_data(); unsigned i,it; if (numterm) setup_mme(&rellrhs); blupsol.resize(hmmesize + ntermGdist*popsize, 0.0); Vector<double> tmpsol(hmmesize); tmpsol.assign(0.0); Vector<double> wspace; wspace.reserve(hmmesize); double * ve = &(blupsol[hmmesize]); if (ntermGdist) { pop->residual_var = residual_var; pop->input_data(data); pop->cond_genotype_config(); } ///////////////////////////////////////////////////// // Gibbs Sampler warm up ///////////////////////////////////////////// for (it=0; it<warmup; it++) { if (numterm) hmmec.gibbs_iterate(tmpsol.begin(),rellrhs.begin(),wspace.begin()); for (i=0; i<ntermGdist; i++) { compute_xbzu(tmpsol); // to adjust trait pop->gibbs_iterate(); } } ///////////////////////////////////////////////////// // samples from Gibbs Sampler now will be used //////////////////////////////////////////////////// for (it=0; it<gibbslen; it++) { if (numterm) { hmmec.gibbs_iterate(tmpsol.begin(),rellrhs.begin(),wspace.begin()); for (i=0; i<hmmesize; i++) blupsol[i] += tmpsol[i]; } for (i=0; i<ntermGdist; i++) { compute_xbzu(tmpsol); // to adjust trait pop->gibbs_iterate(); for (i=0; i<popsize; i++) { ve[i] += pop->popmember[i]->genotypic_val(); } } } blupsol /= gibbslen; hmmec.resize(0,0); return &(blupsol); }

Vector< double > * matvec::Model::genotypic_value_peeling (  void    )

Definition at line 65 of file model_peeling.cpp. References bad_model, blupsol, compute_xbzu(), data, data_prepared, matvec::Individual::est_GV, matvec::Population::genotype_dist_peeling(), get_blupsol(), hmmec, hmmesize, matvec::Population::input_data(), ntermGdist, pop, matvec::Population::popmember, popsize, prepare_data(), matvec::Vector< double >::print(), rellrhs, residual_var, matvec::Population::residual_var, matvec::Vector< double >::resize(), setup_mme(), matvec::SparseMatrix::solve(), and type. { if (type == bad_model) return 0; if (!data_prepared) { if (!prepare_data()) { type = bad_model; return 0; } } blupsol.resize(hmmesize + ntermGdist*popsize); unsigned i; if (ntermGdist ) { pop->residual_var = residual_var; pop->input_data(data); if (numterm) { get_blupsol(); compute_xbzu(blupsol); // to adjust trait } else { pop->genotype_dist_peeling(0); for (i=0; i<popsize; i++) blupsol[i] = pop->popmember[i]->est_GV; return &(blupsol); } } else { setup_mme(&rellrhs); hmmec.solve(blupsol,rellrhs); return &(blupsol); } /////////////////////////////////////////// // iterate begins now //////////////////////////////////////////// double * ve = &(blupsol[hmmesize]); unsigned niterate = 0; while (niterate++ < 4) { pop->genotype_dist_peeling(0); for (i=0; i<popsize; i++) ve[i] = pop->popmember[i]->est_GV; blupsol.print(); setup_mme(&rellrhs); hmmec.solve(blupsol,rellrhs); } return &(blupsol); }

Vector< double > * matvec::Model::get_blupsol (  const std::string &    method = "ysmp", double    stopval = 0.001, double    relax = 1.0, unsigned    mxiter = 100000 )

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 See also: blup,blue Definition at line 257 of file model_blup.cpp. References bad_model, blupsol, hmmec, matvec::SparseMatrix::nz(), rellrhs, setup_mme(), matvec::SparseMatrix::solve(), and type. Referenced by contrast(), matvec::GLMM::contrast(), CovMat(), estimate(), genotypic_value_peeling(), and get_solutions(). { Vector<double> *retval = 0; if (type == bad_model) throw exception("Model::blup(): bad model"); if (hmmec.nz() == 0) { setup_mme(&rellrhs); if (type == bad_model) throw exception("Model::blup(): bad model"); if (hmmec.solve(blupsol,rellrhs,method,relax,stopval,mxiter)==1) { retval = &blupsol; //hmmec.close(); } } else { retval = &blupsol; } return retval; }

Vector<int> matvec::Model::get_gamete_vector (  unsigned    index )

void matvec::Model::get_lms_kp (  const unsigned    termth, const unsigned    lvlth, const unsigned    yth, Vector< double > &    kp, const int *    lsmtablevec )  [protected]

Definition at line 448 of file model_hypo.cpp. References matvec::SparseMatrix::a(), matvec::Vector< T >::assign(), matvec::Vector< T >::begin(), matvec::Session::epsilon, hmmec, hmmesize, matvec::SparseMatrix::ia(), matvec::SparseMatrix::ja(), numterm, numtrait, matvec::SESSION, and term. Referenced by out_lsmeans_to_stream(). { int i,j,k,h; unsigned l,t,jb,n,nh,nl,term_start,term_nlevel,endaddr; int mcell; k = term[termid].startaddr() + lvlth*numtrait + yth; int *ia = hmmec.ia(); int *ja = hmmec.ja(); double *a = hmmec.a(); kp.assign(0.0); double *kpve = kp.begin(); kpve--; mcell = 0; k++; endaddr = ia[k+1]; for (j=ia[k]; j<endaddr; j++) { if (ja[j] != k) continue; if (a[j] <= SESSION.epsilon) { mcell = 1; return; } else { break; } } //////////////////////////////////////////////////////////////// // set element kp[i] = 1.0 if MME[k][i] is non-zero element // //////////////////////////////////////////////////////////////// for (j=ia[k]; j<endaddr; j++) kpve[ja[j]] = 1.0; for (i=1; i<=hmmesize; i++) { endaddr = ia[i+1]; for (j=ia[i]; j<endaddr; j++) if (ja[j] == k) kpve[i] = 1.0; } for (t=0; t<numterm; t++) { if (!(term[t].classi() == 'F' || term[t].classi() == 'C')) continue; if (lsmtablevec[t] >= 1) continue; term_start = term[t].startaddr() + 1; term_nlevel = term[t].nlevel(); l = term_start + yth; for (i=0; i<term_nlevel; i++) { for (n=ia[l+1], j=ia[l]; j<n; j++) { if (ja[j] == l) { if (fabs(a[j]) > SESSION.epsilon) kpve[l] = 1.0; } } l += numtrait; } } //////////////////////////////////////////////////////////////////// // if term[i] is random, then set all elements kp[j] for term[i] // to be zeros, // else set elements kp[j] for non-relevent trait to be zero //////////////////////////////////////////////////////////////// for (i=0; i<numterm; i++) { term_start = term[i].startaddr() + 1; term_nlevel = term[i].nlevel(); if (term[i].classi() == 'F' || term[i].classi() == 'C') { for (t=0; t<numtrait; t++) { if (t != yth) { l = term_start; for ( j=0; j<term_nlevel; j++) { kpve[l+t] = 0.0; l += numtrait; } } } } else { memset(&(kpve[term_start]),'\0',numtrait*sizeof(double)*term_nlevel); } } double dx,xn; DataNode dnode; for (j=0; j<numterm; j++) { if (lsmtablevec[j] == 2) continue; // j = lsm_term term_nlevel = term[j].nlevel(); term_start = term[j].startaddr() + 1; dx = 1.0; if (term[j].classi() == 'C') { k = factor_struct[term[j].factorindx[0]].index(); dnode = data->datasheet[k].mean(0); if (dnode.missing) throw exception("get_lms_kp(): no valid values in the covariate column"); dx = dnode.double_val(); } if (!(term[j].classi() == 'F' || term[j].classi() == 'C')) continue; if (lsmtablevec[j] == 1) { // j is related to lsm_term n = term_start + term_nlevel*numtrait; for (xn=0.0,k=term_start; k<n; k++) xn += kpve[k]; if (xn >= SESSION.epsilon*10) { for (k=term_start; k<n; k++) kpve[k] *= dx/xn; } } else { if (term[j].order() == 1) { n = term_start + term_nlevel*numtrait; for (xn=0.0, k=term_start; k<n; k++) xn += kpve[k]; if (xn >= SESSION.epsilon*10) { for (k=term_start; k<n; k++) kpve[k] *= dx/xn; } } else { nl = factor_struct[term[j].factorindx[0]].nlevel(); nh = term_nlevel/nl; for (k=0; k<nh; k++) { jb = term_start + k; for (xn=0.0,h=0; h<nl; h++) { xn += kpve[jb]; jb += nh*numtrait; } if (xn < SESSION.epsilon*10) continue; jb = term_start + k; for (h=0; h<nl; h++) { kpve[jb] *= dx/(xn*nh); jb += nh*numtrait; } } } } } }

unsigned matvec::Model::get_popsize (  void    )  const [inline]

Definition at line 339 of file model.h. References popsize. {return popsize;};

Vector< double > matvec::Model::get_solutions (  const std::string &    termname )

Definition at line 2390 of file model.cpp. References blupsol, factor_struct, get_blupsol(), matvec::TermList::index(), numtrait, matvec::Vector< T >::resize(), and term. { // Authors: Liviu R. Totir and Rohan L. Fernando (November, 2002) // Contributors: Vector<double> retval; if (!get_blupsol()) throw exception("Model::get_solutions(): get_blupsol() failed"); int k = term.index(termname,factor_struct); if (k < 0) throw exception("Model::get_solutions(): no such term in the model"); int nlevels = term[k].nlevel()*numtrait; unsigned i,startaddr; startaddr = term[k].startaddr(); retval.resize(nlevels); for (i=0; i<nlevels; i++) { retval[i] = blupsol[startaddr+i]; } return retval; }

matvec::Matrix< double > matvec::Model::getIBDMatrix (  void    )

Definition at line 1713 of file model.cpp. References pop, matvec::Population::popmember, popsize, matvec::Population::popsize, matvec::Matrix< T >::resize(), and matvec::Matrix< T >::transpose(). Referenced by RSamplerGibbs(), RSamplerGibbsMH(), and RSamplerMH(). { // Author: L. Radu Totir // (July, 2004) // Contributors: matvec::Matrix<double> geneticValVec; matvec::Matrix<double> CovMatrix; CovMatrix.resize(pop->popsize,pop->popsize,0.0); matvec::Vector<double> qtlMu; qtlMu.resize(4,0.0); qtlMu[0]=-1.41421,qtlMu[3]=1.41421; geneticValVec.resize(pop->popsize,1,0.0); for (unsigned t=0;t<pop->popsize;t++){ // cout << (pop->popmember[t])->id() << " " // << (pop->popmember[t])->genotNodeVector[1].acceptedGenotypeState // << endl; int geno = (pop->popmember[t])->genotNodeVector[1].acceptedGenotypeState; geneticValVec[t][0] = qtlMu[geno]; } CovMatrix =geneticValVec*geneticValVec.transpose(); return CovMatrix; }

double matvec::Model::graph_log_likelihood_peeling (  const unsigned    maxiter )

Definition at line 1155 of file model_peeling.cpp. References matvec::Population::analysis_type, bad_model, matvec::Vector< double >::begin(), blupsol, compute_xbzu(), data, data_prepared, matvec::Individual::est_GV, matvec::SparseMatrix::factorization(), matvec::Population::genotype_dist_peeling(), hmmec, hmmesize, matvec::SparseMatrix::info_vec, matvec::Population::input_data(), max_nz, ntermGdist, numobs, numterm, pop, matvec::Population::popmember, popsize, prepare_data(), matvec::SparseMatrix::q(), rellrhs, matvec::SparseMatrix::reorder(), matvec::Population::residual_var, residual_var, matvec::Vector< double >::resize(), matvec::SparseMatrix::resize(), setup_ww_single_trait(), matvec::SparseMatrix::solve(), term, matvec::Population::tn_qtl, totalnvc(), type, var2vec(), and matvec::Individual::xbzu_val. { if (ntermGdist == 0) { throw exception("Model::graph_log_likelihood_peeling(): ntermGdist == 0"); return 0.0; } if (type == bad_model) { throw exception("Model::graph_log_likelihood_peeling(): bad model"); return 0.0; } if (!data_prepared) { if (!prepare_data()) { type = bad_model; return 0.0; } } int pt,i; Individual *I; double sigma_mu,retval = 0.0; double *sol, rr; unsigned nl,rank_mme,startaddr,iter = 0; pop->analysis_type = "not-multi"; pop->tn_qtl=4; if (numterm) { retval += numobs*std::log(residual_var[0][0]); hmmec.resize(hmmesize,max_nz); rellrhs.resize(hmmesize); blupsol.resize(hmmesize); pop->genotype_dist_peeling(0); for (i=0; i<popsize; i++) { I = pop->popmember[i]; I->xbzu_val = I->est_GV; } unsigned nvc = totalnvc() - 1; // number of ratio's Vector<double> zz(popsize); Vector<double> ratio(nvc); Vector<double> vc(nvc+1); var2vec(vc); for (i=0; i<nvc; i++) ratio[i] = vc[nvc]/vc[i]; for (pt=-1,i=0; i<numterm; i++) { if (term[i].classi() == 'P') { pt = i; startaddr = term[pt].start + 1; nl = term[pt].nlevel(); break; } } setup_ww_single_trait(ratio,pt,startaddr,&zz); hmmec.reorder(); rank_mme = hmmec.factorization(5); hmmec.solve(blupsol,rellrhs,"ysmp1"); retval += hmmec.info_vec[0]- rank_mme*std::log(residual_var[0][0]); if (pt >= 0) { sol = &(blupsol[startaddr-1]); for (rr=0.0,i=0; i<nl; i++) { rr += (*sol * zz[i] * *sol); sol++; } sigma_mu = *(term[pt].prior->var_matrix())[0][0]; retval += nl*std::log(sigma_mu); sigma_mu = hmmec.q(blupsol.begin(),blupsol.begin(),startaddr,startaddr+nl-1) - rr; retval += std::log(sigma_mu); } compute_xbzu(blupsol); // to adjust trait } retval *= -0.5; pop->residual_var = residual_var; pop->input_data(data); // retval += pop->graph_log_likelihood_peeling(maxiter); return retval; }

void matvec::Model::hashxact (  const Vector< bool > &    )  [protected]

Definition at line 1002 of file model.cpp. References matvec::Field::classi(), data, matvec::Data::datasheet, matvec::DataNode::double_val(), factor_struct, matvec::HashTable::get_id(), idlist, matvec::TermList::index(), matvec::Vector< T >::initialize(), matvec::HashTable::insert(), matvec::Data::num_rows(), numcol, numterm, rawpos, matvec::Data::row(), term, matvec::DataNode::unsigned_val(), and xact_htable. Referenced by assign_id_xact(). { if (!data) return; char CL; int i,j,t,k,nord; int *tempval = new int; Vector<unsigned> pos_datacol(numcol); DataNode* recd; if(numcol>0){ recd = new DataNode [numcol]; } else { recd = 0; } unsigned nrow_in_data = data->num_rows(); for (i=0; i<nrow_in_data; i++) { if (record_legality[i]) continue; data->row(i,recd); pos_datacol.initialize(numcol,0); for (t=0; t<numterm; t++) { nord = term[t].order(); if (nord > 1) { for (j=0; j<nord; j++) { k = factor_struct[term[t].factorindx[j]].index(); if (pos_datacol[k] == 0) { CL = data->datasheet[k].classi(); if (CL == 'F') { *tempval = static_cast<int>(recd[k].double_val()); pos_datacol[k] = idlist[k]->get_id(tempval); } else if (CL == 'U' ||CL == 'P') { pos_datacol[k] = recd[k].unsigned_val(); } else if (CL =='C') { pos_datacol[k] = 1; } } rawpos[j] = pos_datacol[k]; } xact_htable[t].insert(rawpos); } } } if (recd) { delete [] recd; recd = 0; } if(tempval){ delete tempval; tempval = 0; } }

void matvec::Model::info (  std::ostream &    stream )  const

Definition at line 1017 of file model_vce.cpp. References bad_model, matvec::ModelTerm::classi(), dfreml_method, factor_struct, matvec::ModelTerm::factorindx, hmmec, hmmesize, matvec::FieldStruct::mean(), matvec::FieldStruct::name(), nfunk_in_dfreml, matvec::FieldStruct::nlevel(), numterm, numtrait, matvec::SparseMatrix::nz(), matvec::Session::output_precision, matvec::Matrix< double >::print(), matvec::ModelTerm::prior, reml_value, residual_var, matvec::SESSION, matvec::FieldStruct::std(), term, trait_struct, type, and matvec::GeneticDist::var_matrix(). Referenced by matvec::funk(), lsmeans(), save(), vce_dfreml(), vce_emreml_multi_trait(), vce_emreml_single_trait(), and vce_gibbs(). { if (type == bad_model) throw exception("Model::info(stream): model is too bad"); stream.precision(SESSION.output_precision); stream <<"\n some extra information in the model\n"; stream << " --------------------------------------------------------\n"; if (dfreml_called) { stream << "optimizor used in dfreml = " << dfreml_method <<"\n"; stream <<"# of function evaluations used in dfreml = " << nfunk_in_dfreml<<"\n"; stream << "maximum log restricted likelihood in dfreml = " << reml_value <<"\n"; } else { ; } const ModelTerm *T; int i; for (i=0; i<numterm; i++) { T = &(term(i)); if (T->classi() == 'R' || T->classi() == 'P') { stream << " variance for " << factor_struct[T->factorindx[0]].name() << " = \n"; T->prior->var_matrix()->print(stream); } } stream << " residual variance =\n" << residual_var; if (data_prepared) { stream << "\n"; stream << " MME dimension : " << hmmesize << "\n"; stream << " non-zeros in MME: " << hmmec.nz() << "\n\n"; stream << " basic statistics for dependent variables\n"; stream << " -----------------------------------------------------\n"; stream << " trait-name n mean std\n"; unsigned n; for (i=0; i<numtrait; i++) { n = trait_struct[i]->nlevel(); stream << " " << std::setw(16) << trait_struct[i]->name() << " " << std::setw(10) << n << " " << std::setw(12) << trait_struct[i]->mean(); if (n==1) { stream << " " << std::setw(12) << "." << "\n"; } else { stream<< " " << std::setw(12) << trait_struct[i]->std() << "\n"; } } } stream << " --------------------------------------------------------\n"; }

void matvec::Model::initialize (  void    )

Definition at line 35 of file model.cpp. References bad_model, data, data_prepared, dfreml_called, factor_struct, hmmesize, idlist, kvec, max_nz, maxorder, modelpcount, modelstr, nfunk_in_dfreml, non_zero, npattern, ntermGdist, num_ped, numcol, numfact, numgroup, numobs, numrec, numterm, numtrait, pop, pop_created, popsize, pos_term, rawpos, rec_indid, reml_value, rve, trait_struct, trait_vec, traitname, type, unknown_prior, weightinverse, xact_htable, and xval_term. Referenced by copyfrom(), and Model(). { num_ped = 0; max_nz = 0; numrec = 0; numobs = 0; numcol = 0; numterm = 0; popsize = 0; numfact = 0; maxorder = 0; numtrait = 0; npattern = 0; hmmesize = 0; numgroup = 0; non_zero = 0; type = bad_model; // it becomes good mode after having equations - was model_type ntermGdist = 0; pop_created = 0; data_prepared = 0; weightinverse = 0; modelpcount =0; modelstr = 0; pop = 0; rve = 0; data = 0; kvec = 0; rawpos = 0; idlist = 0; pos_term = 0; xval_term = 0; traitname = 0; trait_vec = 0; rec_indid = 0; xact_htable = 0; factor_struct = 0; trait_struct = 0; unknown_prior = 0; // ********************************************** // information variables for Model::info(stream) // ********************************************** nfunk_in_dfreml = 0; reml_value = 0.0; dfreml_called = 0; }

void matvec::Model::input_pos_val (  std::istream &    modelfile )  [protected]

Input position and value for each term in the model Note that term[k].addr[t] starts from 1 rather than from 0 Definition at line 367 of file model_blup.cpp. References base_effect, nt_vec, numterm, numtrait, pos_term, term, trait_vec, and xval_term. Referenced by matvec::GLMM::residual(), sampleW(), setup_ww(), matvec::GLMM::setup_ww(), setup_ww_single_trait(), and matvec::GLMM::SSQCP(). { int t,k,nt; unsigned startaddr; modelfile.read((char *)pos_term, (numterm+1)*sizeof(unsigned)); modelfile.read((char *)xval_term, (numterm+1)*sizeof(double)); modelfile.read((char *)trait_vec, numtrait*sizeof(double)); for (k=0; k<numterm; k++) { nt=nt_vec[base_effect[k]]; startaddr = term[k].start + 1 +(pos_term[k]-1)*nt ; for (t=0; t<numtrait; t++) { if (term[k].trait[t]) term[k].addr[t] = startaddr + t; else term[k].addr[t] = 0; } } }

void matvec::Model::input_pos_val_iod (  std::vector< pos_val_node >::iterator    it )  [protected]

Definition at line 385 of file model_blup.cpp. References base_effect, nt_vec, numterm, numtrait, and term. Referenced by compute_rhs_diag_mme(), and mme_times(). { int t,k,nt; unsigned startaddr; for (k=0; k<numterm; k++) { nt=nt_vec[base_effect[k]]; startaddr = term[k].start + 1 + (it->pos_term[k]-1)*nt ; for (t=0; t<numtrait; t++) { if (term[k].trait[t]) term[k].addr[t] = startaddr + t; else term[k].addr[t] = 0; } } }

void matvec::Model::inverse_residual_var (  void    )  [protected]

Compute the inverse of residual variance matrix Definition at line 594 of file model_blup.cpp. References matvec::Matrix< double >::begin(), matvec::Session::epsilon, matvec::ginverse1(), lnr0vec, npattern, numtrait, pattern, residual_var, rve, and matvec::SESSION. Referenced by blup_pccg(), matvec::GLMM::residual(), and setup_mme(). { unsigned i,t,t1,t2; double **ve; for (i=0; i<npattern; i++) { ve = rve[i].begin(); for (t1=0; t1<numtrait; t1++) for (t2=0; t2<numtrait; t2++) { ve[t1][t2] = residual_var[t1][t2]; } for (t=0; t<numtrait; t++) { if (pattern[i][t]=='0') { for (t1=0; t1<numtrait; t1++) ve[t1][t] = 0.0; // cross col t for (t2=0; t2<numtrait; t2++) ve[t][t2] = 0.0; // cross row t } } ginverse1(ve,numtrait,lnr0vec[i],1,SESSION.epsilon); } }

std::string matvec::Model::label (  const std::string &    termname, const unsigned    i )

Definition at line 260 of file model_hypo.cpp. References factor_struct, matvec::TermList::index(), term, and trait_effect_level(). { std::string rawcode; if(termname==""){ trait_effect_level(i,rawcode); } else{ int k = term.index(termname,factor_struct); if (k < 0) throw exception("Model::label(): no such term in the model"); unsigned startaddr = term[k].startaddr(); // cout << startaddr+i << std::endl; trait_effect_level(startaddr+i,rawcode); } return(rawcode); }

double matvec::Model::log_likelihood_gibbs (  const unsigned    wup = 30, const unsigned    glen = 1000 )

Definition at line 71 of file model_gibbs.cpp. References matvec::warning(). { warning("Model.log_likelihood_gibbs(): not inplemented yet"); return 0; }

double matvec::Model::log_likelihood_peeling (  const unsigned    maxit = 3 )

Definition at line 195 of file model_peeling.cpp. References matvec::Population::analysis_type, bad_model, matvec::Vector< double >::begin(), blupsol, compute_xbzu(), data, data_prepared, matvec::Individual::est_GV, matvec::SparseMatrix::factorization(), matvec::Population::genotype_dist_peeling(), hmmec, hmmesize, matvec::SparseMatrix::info_vec, matvec::Population::input_data(), matvec::Population::log_likelihood_peeling(), max_nz, ntermGdist, numobs, numterm, pop, matvec::Population::popmember, popsize, prepare_data(), matvec::SparseMatrix::q(), rellrhs, matvec::SparseMatrix::reorder(), matvec::Population::residual_var, residual_var, matvec::Vector< double >::resize(), matvec::SparseMatrix::resize(), restricted_log_likelihood(), setup_ww_single_trait(), matvec::SparseMatrix::solve(), term, matvec::Population::tn_qtl, totalnvc(), type, var2vec(), and matvec::Individual::xbzu_val. { if (type == bad_model) throw exception("Model::log_likelihood_peeling(): bad model"); if (ntermGdist == 0) { return restricted_log_likelihood(); } if (!data_prepared) { if (!prepare_data()) { type = bad_model; return 0.0; } } int pt,i; Individual *I; double sigma_mu,retval = 0.0; double *sol, rr; unsigned nl,rank_mme,startaddr,iter = 0; pop->analysis_type = "not-multi"; pop->tn_qtl=4; if (numterm) { retval += numobs * std::log(residual_var[0][0]); hmmec.resize(hmmesize,max_nz); rellrhs.resize(hmmesize); blupsol.resize(hmmesize); pop->genotype_dist_peeling(0); for (i=0; i<popsize; i++) { I = pop->popmember[i]; I->xbzu_val = I->est_GV; } unsigned nvc = totalnvc() - 1; // number of ratio's Vector<double> zz(popsize); Vector<double> ratio(nvc); Vector<double> vc(nvc+1); var2vec(vc); for (i=0; i<nvc; i++) ratio[i] = vc[nvc]/vc[i]; for (pt=-1,i=0; i<numterm; i++) { if (term[i].classi() == 'P') { pt = i; startaddr = term[pt].start + 1; nl = term[pt].nlevel(); break; } } setup_ww_single_trait(ratio,pt,startaddr,&zz); hmmec.reorder(); rank_mme = hmmec.factorization(5); hmmec.solve(blupsol,rellrhs,"ysmp1"); retval += hmmec.info_vec[0]- rank_mme * std::log(residual_var[0][0]); if (pt >= 0) { sol = &(blupsol[startaddr-1]); for (rr=0.0,i=0; i<nl; i++) { rr += (*sol * zz[i] * *sol); sol++; } sigma_mu = *(term[pt].prior->var_matrix())[0][0]; retval += nl * std::log(sigma_mu); sigma_mu = hmmec.q(blupsol.begin(),blupsol.begin(),startaddr,startaddr+nl-1) - rr; retval += std::log(sigma_mu); } compute_xbzu(blupsol); // to adjust trait } retval *= -0.5; pop->residual_var = residual_var; pop->input_data(data); retval += pop->log_likelihood_peeling(maxiter); return retval; }

void matvec::Model::lsmeans (  const std::string &    termnames, const std::string &    filename = "", const int    io_mode = std::ios::out, const int    savekp_flag = 0 )

Definition at line 576 of file model_hypo.cpp. References blup(), factor_struct, matvec::TermList::index(), info(), numterm, out_lsmeans_to_stream(), term, and matvec::warning(). { /////////////////////////////////////////////////////////////////////// // currently not working together with group-model // A*B works, but not A * B, we will fix it later // // REMINDER: blupsol & rellrhs remain intact /////////////////////////////////////////////////////////////////////// // lsm_terms.split(nlsm," "); std::string sep(" "); std::vector<std::string> tmpstr; std::string::size_type begidx,endidx; begidx = termnames.find_first_not_of(sep); while(begidx != std::string::npos) { endidx = termnames.find_first_of(sep,begidx); if (endidx == std::string::npos) endidx = termnames.length(); tmpstr.push_back(termnames.substr(begidx,endidx - begidx)); begidx = termnames.find_first_not_of(sep,endidx); } unsigned nlsm = tmpstr.size(); if (nlsm == 0) { warning("Model::lsmeans(termnames): termnames is empty"); return; } if (!blup("ysmp")) throw exception("Model::lsmeans(): bad model"); Matrix<int> lsmtable(nlsm,numterm+1); unsigned term_nlevel,i,j,t; unsigned nord = 0; int k,kk; sep = "*"; std::vector<std::string> effectnames; for (term_nlevel = 0, i=0; i<nlsm; i++) { k = term.index(tmpstr[i],factor_struct); if (k >= 0) { if (term[k].classi() != 'F') { warning("Model::lsmeans(%s): only fixed(discrete) term allowed", tmpstr[i].c_str()); lsmtable[i][numterm] = -1; // lsm_term is discarded continue; } } else { warning("Model::lsmeans(%s): no such term in the model, it's ignored", tmpstr[i].c_str()); lsmtable[i][numterm] = -1; // lsm_term can't be found in the model continue; } lsmtable[i][k] = 2; // term k is the lsmeans term lsmtable[i][numterm] = k; term_nlevel += term[k].nlevel(); effectnames.clear(); //////// effectnames = tmpstr[i].split(nord,"*"); begidx = tmpstr[i].find_first_not_of(sep); while(begidx != std::string::npos) { endidx = tmpstr[i].find_first_of(sep,begidx); if (endidx == std::string::npos) endidx = tmpstr[i].length(); effectnames.push_back(tmpstr[i].substr(begidx,endidx - begidx)); begidx = tmpstr[i].find_first_not_of(sep,endidx); } nord = effectnames.size(); for (t=0; t<numterm; t++) { if (term[t].order() > 2) { lsmtable[i][numterm] = -1; throw exception("Model::lsmeans(args): can't handle high order interaction"); } if (lsmtable[i][t] == 2) continue; for (j=0; j<nord; j++) { kk = term[t].partial_match(effectnames[j],factor_struct); if (kk>=0) break; } if (j<nord) lsmtable[i][t] = 1; } } if (filename == "") { std::ofstream ofs; ofs.open(filename.c_str(),(OpenModeType)io_mode); if (!ofs) throw exception("Model::lsmeans(): cann't open or already exist"); out_lsmeans_to_stream(ofs,nlsm,lsmtable,savekp_flag); info(ofs); ofs.close(); } else { out_lsmeans_to_stream(std::cout,nlsm,lsmtable,savekp_flag); } }

void matvec::Model::map (  void    )

double matvec::Model::MapF (  double    dist, int    qr = 0 )

Definition at line 527 of file model_peeling.cpp. References map_function, and map_parm. Referenced by matvec::Population::buildRecombinationMatrix(), matvec::Population::descent_graph_setup(), and multipoint_compute_Rec_table(). { // Computes different map functions that are multilocus feasible: // Model::map_method M is Morgans // Model::map_method H (default) is Haldanes // Model::map_method B Binomial with extra_parm // qr is a flag // 0 means theta is a recombination rate // 1 means theta is a distance in Morgans no centiMorgrans double temp; if ((map_function=='M') || (map_function=='m')) { //Morgan's function recombination rate equals distance return theta; } else if ((map_function=='B') || (map_function=='b')) { // Binomial Method if (qr) { temp=std::pow((1.0-2.0*theta),(1.0/map_parm)); return (0.5*map_parm*(1.0-temp)); } else { if (theta < map_parm/2.0) { temp=(1.0-2.0*theta/map_parm); return 0.5*(1.0-(std::pow(temp,map_parm))); } else return 0.5; } } else { // Haldane's function if method not set if (qr) { if (theta < 0.5) { temp= (-0.5 * std::log(1.0-2.0*theta)); return temp; } else return 5.0; } else { temp= (0.5*(1.0 - std::exp(-2.0*theta))); return temp; } } }

void matvec::Minimizer::min_prtlevel (  const int    pl )  [inline, inherited]

Definition at line 50 of file minimizer.h. References matvec::Minimizer::prtlevel. {prtlevel = pl;}

double matvec::Model::minfun (  const Vector< double > &    x, const int    n )  [virtual]

Implements matvec::Minimizer. Definition at line 2622 of file model.cpp. References matvec::Minimizer::minfun_indx, minfun_multipoint(), minfun_peeling(), minfun_vce(), and matvec::warning(). { double llhood = 0.0; switch (minfun_indx) { case 1: // dfreml vce, Model_vce.cc llhood = minfun_vce(x,n); break; case 2: // MLE peeling, Model_peeling llhood = minfun_peeling(x,n); break; // BRS case 3: // Multipoint llhood = minfun_multipoint(x,n); break; // BRS default: warning("Model::Model::minfun(): unknown minfun_indx: %d",minfun_indx); } return llhood; }

double matvec::Model::minfun_multipoint (  const Vector< double > &    x, const int    n )

Definition at line 770 of file model_peeling.cpp. References matvec::Population::analysis_type, matvec::GeneticDist::chrom(), matvec::ChromStruct::locus, maxit, matvec::Population::mean_for_genotype, matvec::Population::multi_m_log_likelihood_peeling(), matvec::Population::multipoint_likelihood(), multipoint_QTL_table(), multipoint_Rec_recalc(), N_Parm_List, new_distance, Parm_List, pop, matvec::Population::prior, and matvec::Population::Q_freq. Referenced by minfun(). { int i; // std::cout << "entered minfun_multipoint " << std::endl; // setup values from praxis into correct places for (i=0; i < N_Parm_List; i++) { *(Parm_List[i].Value) = x[i]; // std::cout << " Parameter " << i << " is " << x[i] << std::endl; if ((x[i] > Parm_List[i].Max) || (x[i] < Parm_List[i].Min)) { std::cout << " Parameter " << i << " is out of bounds " << x[i] << " Max is " << Parm_List[i].Max; std::cout << " Min is " << Parm_List[i].Min << std::endl; return 1.0e+25; } } // First need to find the new location and position of the QTL and recompute QTL table // std::cout << "New dist " << new_distance; // << "RecoTable Before change " << RecoTable; multipoint_Rec_recalc(new_distance); // std::cout << "After change " << RecoTable; // std::cout << " POP Q FREQS " << pop->Q_freq; // Recompute QTL genotype frequencies double QTL_freq= pop->prior->chrom()[0].locus[0].allele_freq[0]; pop->Q_freq(1) = QTL_freq*QTL_freq; pop->Q_freq(2) = QTL_freq*(1-QTL_freq); pop->Q_freq(3) = QTL_freq*(1-QTL_freq); pop->Q_freq(4) = (1-QTL_freq)*(1-QTL_freq); // End of old Stuff multipoint_QTL_table(); int N_means=5; // Find which means was used for (i=1; i <= N_Parm_List; i++) { // std::cout << "Start:" <<Parm_List[i].Name << ":End" << endl; if (Parm_List[i].Name == "Mean_Qq ") { N_means= N_means-2; } if (Parm_List[i].Name == "Dom_Qq ") { N_means--; } if (Parm_List[i].Name == "Mean_QQ ") { N_means--; } } //Dominance Means if (N_means == 3) { if (pop->mean_for_genotype[0] > pop->mean_for_genotype[3]){ if (pop->mean_for_genotype[1] > pop->mean_for_genotype[0]){ pop->mean_for_genotype[1] = pop->mean_for_genotype[0]; } else if (pop->mean_for_genotype[3] > pop->mean_for_genotype[1]){ pop->mean_for_genotype[1] = pop->mean_for_genotype[3]; } } else { if (pop->mean_for_genotype[1] > pop->mean_for_genotype[3]){ pop->mean_for_genotype[1] = pop->mean_for_genotype[3]; } else if (pop->mean_for_genotype[0] > pop->mean_for_genotype[1]){ pop->mean_for_genotype[1] = pop->mean_for_genotype[0]; } } } //Additive means else if (N_means == 4) { (pop->mean_for_genotype[1]) = 0.5*((pop->mean_for_genotype[0])+ (pop->mean_for_genotype[3])); } //Same means else if (N_means == 5) { (pop->mean_for_genotype[1]) = (pop->mean_for_genotype[0]); (pop->mean_for_genotype[3]) = (pop->mean_for_genotype[0]); } // General model don't do anything! //Force heterozygotes to have same mean (pop->mean_for_genotype[2]) = (pop->mean_for_genotype[1]); // cout << "mean 0 = " << pop->mean_for_genotype[0] << endl; // cout << "mean 1 = " << pop->mean_for_genotype[1] << endl; // cout << "mean 2 = " << pop->mean_for_genotype[2] << endl; // cout << "mean 3 = " << pop->mean_for_genotype[3] << endl; // std::cout << "parameters set " << std::endl; double log_llhood; if (pop->analysis_type == "multipoint") { // std::cout << "doing " << pop->analysis_type << " multipoint" << std::endl; log_llhood= -1.0*pop->multipoint_likelihood(maxit); } else { // std::cout << "doing " << pop->analysis_type << " multipoint" << std::endl; log_llhood= -1.0*pop->multi_m_log_likelihood_peeling(maxit); } // std::cout << "log likelihood value (its negative) = " << std::setprecision(14); // std::cout << log_llhood; // for (i=0; i < N_Parm_List; i++) // std::cout << "(" << i << ")" << " = " << x[i] << " "; // exit(1); return log_llhood; }

double matvec::Model::minfun_peeling (  const Vector< double > &    x, const int    n )

Definition at line 151 of file model_peeling.cpp. References matvec::Genome::chromosome, matvec::Chromosome::freq(), matvec::Population::log_likelihood_peeling(), matvec::Population::mean_for_genotype, pop, matvec::Population::pop_gamete, matvec::Population::residual_var, matvec::Population::tn_gamete, and matvec::Population::tn_qtl. Referenced by minfun(). { // requirement: build_pop_gamete() must have been called int i,m = pop->tn_gamete - 1; // because p(A1)+ p(A2) +...+ p(At) = 1 double p,pm; ////////////////////////////////////////////////////////////// // constraints for minimization // (0) residual_variance >= 0; // (1) 0.0 <= p(Ai) <= 1 ////////////////////////////////////////////////////////////// if (x[0] < 0.0) return 1.0e+25; for (pm=1.0, i=0; i<m; i++) { pm -= p = x[1+i]; if (p>1.0 || p<0.0) return 1.0e+25; } if (pm>1.0 || pm < 0.0 ) return 1.0e+25; ////////////////////////////////////////////////////////////// // now put back new values in x into appropriate places ////////////////////////////////////////////////////////////// //RLF /* pop->residual_var.me[0][0] = x[0]; for (i=0; i<m; i++) pop->pop_gamete[0].chromosome[i].freq(x[i+1]); pop->pop_gamete[0].chromosome[m].freq(pm); // pop->prior->genotypic_val(1,1,&x[1+m]); // get genotypic_values from x */ pop->residual_var[0][0] = x[0]; pop->pop_gamete[0].chromosome[0].freq(x[1]); pop->pop_gamete[0].chromosome[1].freq(x[2]); pop->pop_gamete[0].chromosome[2].freq(pm); pop->mean_for_genotype[0] = x[3]; pop->mean_for_genotype[1] = x[4]; pop->mean_for_genotype[2] = x[5]; pop->tn_qtl=4; //RLF double log_llhood = -1.0*pop->log_likelihood_peeling(); // std::cout << "log likelihood value (its negative) = " << log_llhood << std::endl; return log_llhood; }

double matvec::Model::minfun_vce (  const Vector< double > &    x, const int    n )

Definition at line 143 of file model_vce.cpp. References numterm, matvec::doubleMatrix::psd(), residual_var, restricted_log_likelihood(), term, and vec2var(). Referenced by minfun(). { //////////////////////////////////////////////////////////////////////// // it returns negative restricted log likelihood. // (it is for minimization program) // // x contains upper triangular part of variance matrices for residual // and other random terms in the model. // the hiden order in x is : // residual--randtermlist[0]--randtermlist[1] ..... // the order of randtermlist[i] is based on "random = A B A*B" given by user // user is totally responsible for size of vector x. // size of x == (1+nrandterm)*numtrait*(numtrait+1)/2 // // Important Notice: // because restricted_log_likelihood(D) will be called sequentially // under the same MME structure. // hmmec is kept after restricted_log_likelihood(D) call. // Call release_mme() to release hmmec from memory /////////////////////////////////////////////////////////////////////////// double mlhood; int i; vec2var(x); if (!residual_var.psd()) return (1.0e30); for (mlhood=0.0,i=0; i<numterm; i++) { if (term[i].classi() == 'R' || term[i].classi() == 'P') { if (!(term[i].prior->var_matrix()->psd())) { mlhood =1.0e30 + 1.0e10*(i+1); break; } } } if (mlhood == 0.0) mlhood = -1.0*restricted_log_likelihood(); return ( mlhood ); }

double matvec::Model::ml_estimate_peeling (  void    )

void matvec::Model::mme_times (  Vector< double > &    x )  [protected]

Compute mme*x without storing mme The result is stored in mme_times_res Definition at line 791 of file model_blup.cpp. References add_Ag(), add_Ig(), matvec::Vector< T >::begin(), matvec::Vector< double >::begin(), hmmesize, input_pos_val_iod(), mme_times_res, numterm, numtrait, pos_val_vector, rve, and term. Referenced by blup_pccg(). { std::vector<pos_val_node>::iterator it; unsigned i,j,ii,jj,t1,t2; double vy,xval,val,value; double *r = mme_times_res.begin(); memset(r,'\0',sizeof(double)*hmmesize); r--; double *xp = x.begin() - 1; double **vep; Matrix<double> ve(numtrait,numtrait); for (it=pos_val_vector.begin(); it!=pos_val_vector.end(); it++) { if(!it->pos_term.empty()){ input_pos_val_iod(it); } else{ continue; } if (ntermGdist) { throw exception("Model::mme_times: something does not seem to be compatible with iod"); } vep = rve[it->pos_term[numterm]].begin(); val = it->xval_term[numterm]; // val = weight variable for (t1=0; t1<numtrait; t1++) { for (t2=0; t2<numtrait; t2++) ve[t1][t2] = val*vep[t1][t2]; } for (t1=0; t1<numtrait; t1++) { for (i=0; i<numterm; i++) { ii = term[i].addr[t1]; if (ii == 0) continue; xval = it->xval_term[i]; for (t2=0; t2<numtrait; t2++) { val = xval * ve[t1][t2]; for (j=i; j<numterm; j++) { jj = term[j].addr[t2]; if (jj >= ii) { value = val*it->xval_term[j]; r[ii] += value*xp[jj]; if (jj > ii) r[jj] += value*xp[ii]; } } } } } } for (int t=0; t<numterm; t++) { if (term[t].classi() == 'P') { add_Ag(t,&x); } else if (term[t].classi() == 'R') { add_Ig(t,&x); } } }

SparseMatrix* matvec::Model::mmec (  void    )  [inline]

Definition at line 304 of file model.h. Referenced by matvec::GLMM::glim(). {return &hmmec;}

unsigned matvec::Model::mmesize (  void    )  const [inline]

Definition at line 192 of file model.h. References hmmesize. {return hmmesize;}

double matvec::Model::multipoint (  int    maxit )

Definition at line 322 of file model_peeling.cpp. References matvec::Population::analysis_type, matvec::Population::F, matvec::Population::multipoint_init_parm(), matvec::Population::multipoint_likelihood(), and pop. { double llh; pop->analysis_type = "multipoint"; // This is a bit excessive as means we have two versions of Farg when doing FPMM Fpmm Farg(0,0.5,1.0); // setup for 0 polygenic loci. if (pop->F == 0) { // Need to setup F as not using FPMM stuff pop->multipoint_init_parm(Farg); } llh=pop->multipoint_likelihood(maxiter); // std::cout << "QTL pos " << Q_pos << std::endl; std::cout << "Natural log llh= " << std::setprecision(10) << llh << " Log10 llh = " << std::setprecision(10) << (std::log10(std::exp(1.0))*llh) << std::endl; /* multipoint_Rec_recalc(6); multipoint_QTL_table(); llh=pop->multipoint_likelihood(maxit); // std::cout << "QTL pos " << Q_pos << std::endl; std::cout << "Natural log llh= " << std::setprecision(10) << llh << " Log10 llh = " << std::setprecision(10) << (log10(exp(1.0))*llh) << std::endl; */ return llh; }

void matvec::Model::multipoint_compute_Rec_table (  void    )

Definition at line 594 of file model_peeling.cpp. References MapF(), matvec::Population::n_markerLoci, pop, Q_pos, and RecoTable. Referenced by multipoint_Rec_recalc(), and multipoint_Rec_table(). { int nloci=pop->n_markerLoci+1, i,j; double temp; int Q_start = Q_pos; for(i=Q_start; i<nloci-1; i++) { if(RecoTable[i][i] > RecoTable[i+1][i+1]) { temp = RecoTable[i][i]; RecoTable[i][i] = RecoTable[i+1][i+1]; RecoTable[i+1][i+1] = temp; Q_pos++; } else { break; } } for(i=Q_start; i>0; i--) { if(RecoTable[i][i] < RecoTable[i-1][i-1]) { temp = RecoTable[i][i]; RecoTable[i][i] = RecoTable[i-1][i-1]; RecoTable[i-1][i-1] = temp; Q_pos--; } else { break; } } for (i=0; i < nloci-1; i++){ for (j=i+1; j < nloci; j++){ RecoTable[i][j]=MapF(RecoTable[j][j]-RecoTable[i][i]); RecoTable[j][i]=1.0-RecoTable[i][j]; // std::cout << "i " << i << " j " << j << " Reco " << RecoTable[i][j] << " RT[j] " << RecoTable[j][j] << " RT[i] " << RecoTable[i][i] << std::endl; } } // std::cout << RecoTable; }

void matvec::Model::multipoint_display_parms (  void    )

Definition at line 1079 of file model_peeling.cpp. References matvec::Population::analysis_type, N_Parm_List, matvec::Population::P_ndim, Parm_List, and pop. { int nl=0; if (pop->analysis_type == "multipoint") { if (pop->P_ndim > 1) { nl= static_cast<int>(0.5*(pop->P_ndim-1)); std::cout << "P_dim" << pop->P_ndim << " nl = " << nl << std::endl; std::cout << "Multipoint with FPMM with "; if (nl >1) { std::cout << nl << " loci." << std::endl; } else { std::cout << "1 locus." << std::endl; } } else std::cout << "Simple Multipoint" << std::endl; } else if (pop->analysis_type== "multipoint_m") { std::cout << "Multipoint with the PAP style mixture approximation" << std::endl; } else std::cout << "what did you run?" << std::endl; //Now pretty print the final values if (N_Parm_List > 0) { std::cout << " Parameter Final Value " << std::endl; int i; for (i=1; i <= N_Parm_List; i++) { if (Parm_List(i).Name != "Distance ") std::cout << Parm_List(i).Name << " " << *Parm_List(i).Value << std::endl; else std::cout << Parm_List(i).Name << " " << *Parm_List(i).Value-2.0 << std::endl; } } }

void matvec::Model::multipoint_estimate_Distance (  double    Initial, double    Max, double    Min )

Definition at line 895 of file model_peeling.cpp. References multipoint_QTL_table(), multipoint_Rec_recalc(), matvec::Population::n_markerLoci, N_Parm_List, new_distance, Parm_List, pop, RecoTable, and matvec::Vector< Parm_Elem >::resize(). { double T_dist,tmp; if (N_Parm_List == 0) { Parm_List.resize(7); } if (Min > Max) { tmp=Min; Min=Max; Max=tmp; } Initial +=2.0; Max +=2.0; Min +=2.0; T_dist = 2.0+RecoTable[pop->n_markerLoci][pop->n_markerLoci]; // CHECK parameter bounds; if (Max > T_dist) Max=T_dist; if (Min < 0.0) Min=0.0; if (Initial > Max) Initial=Max; if (Initial < Min) Initial=Min; new_distance=Initial; N_Parm_List++; Parm_List(N_Parm_List).Value = &new_distance; Parm_List(N_Parm_List).Max = Max; Parm_List(N_Parm_List).Min = Min; Parm_List(N_Parm_List).Name = "Distance "; // Need to recompute tables for the new starting location multipoint_Rec_recalc(new_distance); multipoint_QTL_table(); // std::cout << T_dist << " min " << Min << " Initial " << Initial << " Max " << Max << std::endl; // std::cout << RecoTable ; // exit(1); }

void matvec::Model::multipoint_estimate_PolyV (  double    Initial, double    Max, double    Min )

Definition at line 1115 of file model_peeling.cpp. References matvec::Population::F, N_Parm_List, Parm_List, pop, matvec::Population::residual_var, matvec::Vector< Parm_Elem >::resize(), and matvec::Fpmm::var. { // Maximize the polygenic variance under the mixture model approximation // not used by any other method! double tmp; if (N_Parm_List == 0) { Parm_List.resize(7); // Parm_list is empty therefore resize it } if (Min < 0.0) Min=0.0; if (Min > Max) { tmp=Min; Min=Max; Max=tmp; } if (Initial > Max) Initial=Max; if (Initial < Min) Initial=Min; pop->residual_var[0][0]=Initial; N_Parm_List++; Parm_List(N_Parm_List).Value = &(pop->F->var); Parm_List(N_Parm_List).Max = Max; Parm_List(N_Parm_List).Min = Min; Parm_List(N_Parm_List).Name = "Poly_Var "; }

void matvec::Model::multipoint_estimate_Qfreq (  double    Initial, double    Max, double    Min )

Definition at line 1030 of file model_peeling.cpp. References matvec::GeneticDist::chrom(), matvec::ChromStruct::locus, N_Parm_List, Parm_List, pop, matvec::Population::prior, and matvec::Vector< Parm_Elem >::resize(). { double QTL_allele_freq =Initial,tmp; if (N_Parm_List == 0) { Parm_List.resize(7); } if (Max > 1.0) Max=1.0; if (Min < 0.0) Min=0.0; if (Min > Max) { tmp=Min; Min=Max; Max=tmp; } if (Initial > Max) Initial=Max; if (Initial < Min) Initial=Min; N_Parm_List++; Parm_List(N_Parm_List).Value = &(pop->prior->chrom()[0].locus[0].allele_freq[0]); Parm_List(N_Parm_List).Max = Max; Parm_List(N_Parm_List).Min = Min; Parm_List(N_Parm_List).Name = "QTL_Freq "; }

void matvec::Model::multipoint_estimate_Qmeans (  Vector< double >    Initial, Vector< double >    Max, Vector< double >    Min, char    mtype = 'G' )

Definition at line 932 of file model_peeling.cpp. References matvec::Population::mean_for_genotype, min, N_Parm_List, Parm_List, pop, and matvec::Vector< Parm_Elem >::resize(). { // This maximizes the genotypic means for the genotypes qq, Qq and QQ // only two alleles at QTL supported // Currently force same heterozygotes mean_Qq is set the same as mean_qQ // Future flag maybe added to allow different heterozygote means // Current types // S is same mean : mean qq == mean Qq == mean QQ // A is additive model : mean Qq == average of (mean qq + mean QQ) // D is dominance model: mean qq <= mean Qq <= mean QQ // G or other value has unrestricted means so that overdominance is possible double tmp; int i, j=0; if (N_Parm_List == 0) { Parm_List.resize(7); // Parm_list is empty therefore resize it } // Need to check if any other type of QTL mean maximization have been defined for (i=1; i <= N_Parm_List; i++) { // cout << "Name " << Parm_List(i).Name << endl; if (Parm_List(i).Name == "Mean_qq ") { j=1; } } if (j) { std::cout << "Another estimate Qmeans has been previously defined" << std::endl; throw exception("Model::estimate_Qmeans: has been previously defined"); } else { N_Parm_List++; // For qq genotype pop->mean_for_genotype[0] = Initial[0]; Parm_List(N_Parm_List).Max = Max[0]; Parm_List(N_Parm_List).Min = Min[0]; Parm_List(N_Parm_List).Value = &(pop->mean_for_genotype[0]); Parm_List(N_Parm_List).Name = "Mean_qq "; if ((mtype == 'S') || (mtype=='s')){ // Same means (pop->mean_for_genotype[1]) = (pop->mean_for_genotype[0]); (pop->mean_for_genotype[2]) = (pop->mean_for_genotype[0]); (pop->mean_for_genotype[3]) = (pop->mean_for_genotype[0]); } else { // For QQ genotype N_Parm_List++; pop->mean_for_genotype[3] = Initial[2]; Parm_List(N_Parm_List).Max = Max[2]; Parm_List(N_Parm_List).Min = Min[2]; Parm_List(N_Parm_List).Value = &(pop->mean_for_genotype[3]); Parm_List(N_Parm_List).Name = "Mean_QQ "; if ((mtype=='A') || (mtype=='a') ) { // Additive (pop->mean_for_genotype[1]) = 0.5*((pop->mean_for_genotype[0])+(pop->mean_for_genotype[3])); (pop->mean_for_genotype[2]) = (pop->mean_for_genotype[1]); } else { // set heterozygotes to the same value depending on model N_Parm_List++; if ((mtype == 'D')|| (mtype == 'd')) {// Dominance model Parm_List(N_Parm_List).Name = "Dom_Qq "; // Now check that the Qq means fall within the bounds of QQ and qq means tmp= (double) std::max(Max[0],Max[2]); if (Max[1] > tmp) Max[1]=tmp; tmp=std::min(Min[0],Min[2]); if (Min[1] < tmp) Min[1]=tmp; tmp= (double) std::max(Initial[0],Initial[2]); if (Initial[1] > tmp) Initial[1]=tmp; tmp=std::min(Initial[0],Initial[2]); if (Initial[1] < tmp) Initial[1]=tmp; pop->mean_for_genotype[1] = Initial[1]; Parm_List(N_Parm_List).Max = Max[1]; Parm_List(N_Parm_List).Min = Min[1]; Parm_List(N_Parm_List).Value = &(pop->mean_for_genotype[1]); (pop->mean_for_genotype[2]) = (pop->mean_for_genotype[1]); } else { // General model Parm_List(N_Parm_List).Name = "Mean_Qq "; pop->mean_for_genotype[1] = Initial[1]; Parm_List(N_Parm_List).Max = Max[1]; Parm_List(N_Parm_List).Min = Min[1]; Parm_List(N_Parm_List).Value = &(pop->mean_for_genotype[1]); (pop->mean_for_genotype[2]) = (pop->mean_for_genotype[1]); } pop->mean_for_genotype[2] = Initial[1]; (pop->mean_for_genotype[2]) = (pop->mean_for_genotype[1]); } } } /* * for (i=1; i <= N_Parm_List; i++) { * cout << "Name " << Parm_List(i).Name << " Value= " << &Parm_List(i).Value * << std::endl; * } * */ }

void matvec::Model::multipoint_estimate_Ve (  double    Initial, double    Max, double    Min )

Definition at line 869 of file model_peeling.cpp. References matvec::Matrix< double >::assign(), N_Parm_List, Parm_List, pop, matvec::Population::residual_var, and matvec::Vector< Parm_Elem >::resize(). { double tmp; if (N_Parm_List == 0) { Parm_List.resize(7); } if (Min < 0.0) Min=0.0; if (Min > Max) { tmp=Min; Min=Max; Max=tmp; } if (Initial > Max) Initial=Max; if (Initial < Min) Initial=Min; pop->residual_var.assign(Initial); N_Parm_List++; Parm_List(N_Parm_List).Value = &(pop->residual_var(1,1)); Parm_List(N_Parm_List).Max = Max; Parm_List(N_Parm_List).Min = Min; Parm_List(N_Parm_List).Name = "Res_Var "; // std::cout << " value " << *Parm_List(N_Parm_List).Value << " Min " << Parm_List(N_Parm_List).Min << " Max " << Parm_List(N_Parm_List).Max << " name " << Parm_List(N_Parm_List).Name << std::endl; }

void matvec::Model::multipoint_genot_probs (  void    )

Definition at line 345 of file model_peeling.cpp. References matvec::Population::analysis_type, matvec::Population::build_nufamily(), matvec::Individual::initial_posterior(), matvec::Population::multi_geno_dist_peeling(), pop, matvec::Population::popmember, and matvec::Population::popsize. { int i; Individual *I; std::cout << "computing genot_probs" << std::endl; for (i=0;i<pop->popsize; i++){ I = pop->popmember[i]; //I->initial_anterior(genotype_freq.ve,tn_genotype); I->initial_posterior(3); } pop->analysis_type = "multipoint"; pop->build_nufamily(); pop->multi_geno_dist_peeling(0); }

double matvec::Model::multipoint_mixture_approx (  int    maxit = 0, double    P_var = 1.0 )

Definition at line 1056 of file model_peeling.cpp. References matvec::Population::analysis_type, matvec::Population::multi_m_log_likelihood_peeling(), matvec::Population::multipoint_init_parm(), pop, and Q_pos. { double llh; pop->analysis_type = "multipoint_m"; std::cout << "lets do multipoint with the mixture approximation!" << std::endl; Fpmm Farg(0,0.5,P_var); // setup for 0 polygenic loci. // Need to setup F as not using FPMM stuff pop->multipoint_init_parm(Farg); llh = pop->multi_m_log_likelihood_peeling(maxiter); std::cout << "QTL pos " << Q_pos << std::endl; std::cout << "Natural log llh= " << std::setprecision(10) << llh << " Log10 llh = " << std::setprecision(10) << (std::log10(std::exp(1.0))*llh) << std::endl; /* multipoint_Rec_recalc(14); multipoint_QTL_table(); llh = pop->multi_m_log_likelihood_peeling(maxiter); std::cout << "QTL pos " << Q_pos << std::endl; std::cout << "Natural log llh= " << std::setprecision(10) << llh << " Log10 llh = " << std::setprecision(10) << (log10(exp(1.0))*llh) << std::endl; */ return llh; }

double matvec::Model::multipoint_ml_estimate (  const int    method = 0, int    niter = 0, double    tol = 1.0e-16, double    epsilon = 1.0e-8, int    printlevel = 0 )

Definition at line 722 of file model_peeling.cpp. References matvec::Population::analysis_type, bad_model, data_prepared, maxit, matvec::Minimizer::minfun_indx, N_Parm_List, Parm_List, pop, matvec::Minimizer::praxis(), prepare_data(), and type. { /* * tol is the tolerance allowed for the precision of the * solution. praxis returns if the criterion * 2 * ||x[k]-x[k-1]|| <= sqrt(macheps) * ||x[k]|| + tol * is fulfilled more than ktm times. * the default value 1.0e-16 * tol is usually equal to EPISILON*EPISILON * * printlevel controls praxis printing output: * 0 -> no output * 1 -> print only starting and final values * 2 -> detailed map of the minimization process * 3 -> print also eigenvalues and vectors of the * search directions * the default value is 0 * */ if (type == bad_model) throw exception("Model::ml_estimate_peeling(): bad model"); // Setting up to run praxis to maximise the llh if (method == 0) pop->analysis_type = "multipoint"; else pop->analysis_type = "multipoint_m"; maxit=niter; unsigned n; double max_llhood = 0.0; if (!data_prepared) { if (!prepare_data()) { type = bad_model; return max_llhood; } } minfun_indx = 3; // std::cout << "get x " << N_Parm_List << std::endl; Vector<double> x(N_Parm_List); for (n=0; n < N_Parm_List; n++) x[n]= *(Parm_List[n].Value); int iter = 0; // std::cout << "Go to praxis " << std::endl; max_llhood = praxis(x,n,iter, tol, epsilon, printlevel); minfun_indx = 0; return max_llhood; }

void matvec::Model::multipoint_profile_distance (  int    method, int    maxit, double    Min, double    Max, double    step_size, int    Print, const std::string &    fname )

Definition at line 641 of file model_peeling.cpp. References matvec::Population::analysis_type, matvec::Population::F, matvec::Population::multi_m_log_likelihood_peeling(), matvec::Population::multipoint_init_parm(), matvec::Population::multipoint_likelihood(), multipoint_QTL_table(), multipoint_Rec_recalc(), matvec::Population::n_markerLoci, matvec::Population::P_ndim, pop, and matvec::Matrix< double >::resize(). { int i, nsteps; double distance=0.0, T_dist, temp=0.0; doubleMatrix Results; std::ofstream Outfile (fname.c_str(), std::ios::app); if (!Outfile) throw exception("Cannot open requested filename"); Fpmm Farg(0,0.5,1.0); // setup for 0 polygenic loci. if (pop->F == 0) { // Need to setup F as not using FPMM stuff pop->multipoint_init_parm(Farg); } std::cout << "Nloci " << pop->n_markerLoci << std::endl; T_dist = Max-Min; // T_dist = 2.0+RecoTable[pop->n_markerLoci][pop->n_markerLoci]; temp= (T_dist/step_size); nsteps= (int) (std::ceil(temp+1.0)); std::cout << "Total dist " << T_dist << " With steps " << nsteps << " Of size " << step_size << std::endl; Results.resize(nsteps,2); distance = Min -step_size + 2.0; // std::cout << distance << std::endl; if (method==0) { pop->analysis_type = "multipoint"; if (pop->P_ndim > 1) { std::cout << "Method used : multipoint FPMM" << std::endl; } else { std::cout << "Method used : Simple multipoint" << std::endl; } for (i=0; i < nsteps; i++) { distance += (step_size); Results[i][0] = (distance-2.0); // Find the new location and position of the QTL and recompute QTL table multipoint_Rec_recalc(distance); multipoint_QTL_table(); temp=(pop->multipoint_likelihood(maxiter)); Results[i][1]=temp; } if (Print==1) { for (i=0; i < nsteps; i++) { Outfile << std::setw(10) << std::setprecision(6) << Results[i][0]; Outfile << " " << std::setw(20) << std::setprecision(15); Outfile << Results[i][1] << std::endl; } } } else if (method==1) { pop->analysis_type = "multipoint_m"; std::cout << "Method used : multipoint Mixture approx" << std::endl; for (i=0; i < nsteps; i++) { distance += (step_size); Results[i][0] = (distance-2.0); // Find the new location and position of the QTL and recompute QTL table multipoint_Rec_recalc(distance); multipoint_QTL_table(); temp = pop->multi_m_log_likelihood_peeling(maxiter); Results[i][1]=temp; } if (Print) { for (i=0; i < nsteps; i++) { Outfile << std::setw(10) << std::setprecision(6) << Results[i][0]; Outfile << " " << std::setw(20) << std::setprecision(15); Outfile << Results[i][1] << std::endl; } } } else { std::cout << "Unknown method for multipoint_profile_distance" << std::endl; } if (Print==0) { std::cout << Results << std::endl; } Outfile.close(); }

void matvec::Model::multipoint_QTL_table (  void    )

Definition at line 438 of file model_peeling.cpp. References matvec::Population::n_markerLoci, pop, ProbGamete(), Q_pos, RecoTable, and matvec::Population::switch_table_prob. Referenced by minfun_multipoint(), multipoint_estimate_Distance(), multipoint_profile_distance(), and multipoint_setup(). { /* Try to create the marker and qtl gametic event */ int current, row, EL, TE,i,j,Nel,pow2,pow3; int eindex,sindex,etemp,stemp,eval,sval,locus; double epow,spow,count=0.0; double diff; pow2= (int) (1 << pop->n_markerLoci); // pow(2,pop->n_markerLoci)); pow3= (int) (pow((double) 3,(int)pop->n_markerLoci)); int qindex; doubleMatrix QT(pow3,3); Vector<int> gamete((pop->n_markerLoci)+1); for(eindex=0; eindex < pow3; eindex++) { etemp=eindex; for (locus=(pop->n_markerLoci-1); locus > -1; locus--){ epow=pow((double) 3,(int) locus); eval=(int) (etemp/epow); etemp -= ((int) (epow*eval)); if (((pop->n_markerLoci) - 1 - locus) < (Q_pos))