ttquanti.h

Go to the documentation of this file.

/*
*
*  @file ttquanti.h
*  Nemo2
*
*   Copyright (C) 2006-2023 Frederic Guillaume
*   frederic.guillaume@helsinki.fi
*
*   This file is part of Nemo
*
*   Nemo is free software; you can redistribute it and/or modify
*   it under the terms of the GNU General Public License as published by
*   the Free Software Foundation; either version 2 of the License, or
*   (at your option) any later version.
*
*   Nemo is distributed in the hope that it will be useful,
*   but WITHOUT ANY WARRANTY; without even the implied warranty of
*   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*   GNU General Public License for more details.
*
*   You should have received a copy of the GNU General Public License
*   along with this program; if not, write to the Free Software
*   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*
*  created on @date 14.11.2005
* 
*  @author fred
*/
 
#ifndef TTQUANTI_H
#define TTQUANTI_H
 
#include <cmath>
#include <vector>
#include <set>
#include "ttrait_with_map.h"
#include "bitstring.h"
#include "filehandler.h"
#include "stathandler.h"
#include "metapop.h"
#include "datatable.h"
#include "Uniform.h"
 
#ifdef HAS_GSL
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#endif
 
class TTQuantiSH;
class TTQuantiFH;
class TTQFreqExtractor;
class TTQOhtaStats;
class TProtoQuanti;
 
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuanti : public TTrait {
  
public:
  TTQuanti () : _phenotypes(0), _genotypic_values(0), _myProto(0)
  {
//    message("CONSTRUCTOR::TTQuanti::CONSTRUCTOR\n");
  }
 
  TTQuanti (const TTQuanti& T) : _phenotypes(0), _genotypic_values(0), _myProto(T._myProto){}
//  {  message("COPY CONSTRUCTOR::TTQuanti::COPY CONSTRUCTOR\n"); }
 
  virtual ~TTQuanti () {}
 
  virtual   trait_t         get_type () const                       {return QUANT;}
  virtual   void            mutate ();
  virtual   void            inherit (const TTrait* mother, const TTrait* father);
  virtual   void*           set_trait (void* value)                 {return value;}
  virtual   void            set_value ();
  virtual   void*           getValue () const                       {return _phenotypes;}
 
  // TTQuanti interface:
  virtual double            get_additive_genotype (const unsigned int trait) const = 0;
  virtual double            get_dominant_genotype (const unsigned int trait) const = 0;
  virtual double            get_full_genotype     (unsigned int trait) = 0;
  virtual void              copy_sequence_block   (sex_t SEX, unsigned int strand, unsigned int from_pos,
                                                   unsigned int to_pos, const TTQuanti *parent) = 0;
  virtual void              copy_sequence_1locus  (sex_t SEX, unsigned int strand, unsigned int at,
                                                   const TTQuanti *parent) = 0;
  virtual void              mutate_add            (unsigned int position, unsigned int allele, double value) = 0;
  virtual void              mutate_inplace        (unsigned int position, unsigned int allele, double value) = 0;
 
  // interface to di-allelic architecture
  virtual bool              get_allele_bit (unsigned int position, unsigned int allele) const = 0;
  virtual void              set_allele_bit (unsigned int position, unsigned int allele, bool value) = 0;
 
  //accessors
  void                      set_proto  (TProtoQuanti* proto)        {_myProto = proto;}
  TProtoQuanti*             get_proto  ()                           {return _myProto;}
  double                    get_phenotype         (unsigned int trai);
  void                      set_phenotype         (unsigned int trait, double value);
  void                      reset_phenotype_to_genotypic_value ();
 
protected:
 
  double *_phenotypes;
  double *_genotypic_values;
  TProtoQuanti* _myProto;
};
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuanti_continuous : public TTQuanti  {
 
public:
 
  TTQuanti_continuous () : TTQuanti(), _sequence(0)
  {
//    message("CONSTRUCTOR::TTQuanti_continuous::CONSTRUCTOR\n");
  }
 
  TTQuanti_continuous (const TTQuanti& T) : TTQuanti(T), _sequence(0)
  {
//    message("COPY CONSTRUCTOR::TTQuanti_continuous::COPY CONSTRUCTOR\n");
  }
 
  virtual ~TTQuanti_continuous () {
//    message("DESTRUCTOR::TTQuanti_continuous::DESTRUCTOR\n");
    reset();
  }
  
  //implements TTrait:
  virtual   void            reset ();
  virtual   void            init  ();
  virtual   void            set_sequence (void** seq);
  virtual   void**          get_sequence () const                   {return (void**)_sequence;}
  virtual   unsigned int    get_allele (int loc, int all) const ;
  virtual   double          get_allele_value (int loc, int all) const;
  virtual   void            set_allele_value (unsigned int locus, unsigned int allele, double value);
 
  virtual   TTQuanti_continuous&       operator= (const TTrait& T);
  virtual   bool            operator== (const TTrait& T);
  virtual   bool            operator!= (const TTrait& T);
  
  //implements StorableComponent:
  virtual   void            store_data    ( BinaryStorageBuffer* saver  );
  virtual   bool            retrieve_data ( BinaryStorageBuffer* reader );
  
  //implements TTQuanti:
  virtual double            get_full_genotype     (unsigned int trait);
  virtual void              set_allele            (int locus, int allele, double value) {_sequence[allele][locus] = value;}
  virtual void              mutate_add            (unsigned int position, unsigned int allele, double value)
                               {_sequence[allele][position] += value;}
  virtual void              mutate_inplace        (unsigned int position, unsigned int allele, double value)
                               {_sequence[allele][position] = value;}
  virtual bool              get_allele_bit (unsigned int position, unsigned int allele) const {return false;}
  virtual void              set_allele_bit (unsigned int position, unsigned int allele, bool value) {}
 
  
protected:
    
  double **_sequence;
};
// Deriving specialized classes
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuanti_continuous_full_pleio : public TTQuanti_continuous{
 
public:
  TTQuanti_continuous_full_pleio() : TTQuanti_continuous() {}
 
  TTQuanti_continuous_full_pleio(const TTQuanti_continuous_full_pleio& TT) : TTQuanti_continuous(TT) {}
 
  virtual ~TTQuanti_continuous_full_pleio() {}
 
  //TTtrait implementation:
  virtual void            init_sequence ();
  virtual TTQuanti_continuous_full_pleio*      clone ()     {return new TTQuanti_continuous_full_pleio(*this);}
  virtual void            show_up();
 
  //TTQuanti implementation:
  virtual double          get_additive_genotype (const unsigned int trait) const;
  virtual double          get_dominant_genotype (const unsigned int trait) const;
  virtual void            copy_sequence_block   (sex_t SEX, unsigned int strand, unsigned int from_locus,
                                                 unsigned int to_locus, const TTQuanti *parent);
  virtual void            copy_sequence_1locus  (sex_t SEX, unsigned int strand, unsigned int at,
                                                   const TTQuanti *parent);
};
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuanti_continuous_var_pleio : public TTQuanti_continuous{
 
public:
  TTQuanti_continuous_var_pleio() : TTQuanti_continuous() {}
 
  TTQuanti_continuous_var_pleio(const TTQuanti_continuous_var_pleio& TT) : TTQuanti_continuous(TT) {}
 
  virtual ~TTQuanti_continuous_var_pleio() {}
 
  unsigned int            get_sequence_block_size (unsigned int from, unsigned int to);
 
  //TTtrait implementation:
  virtual void            init_sequence ();
  virtual TTQuanti_continuous_var_pleio*      clone ()     {return new TTQuanti_continuous_var_pleio(*this);}
  virtual void            show_up();
 
  //TTQuanti implementation:
  virtual double          get_additive_genotype (const unsigned int trait) const;
  virtual double          get_dominant_genotype (const unsigned int trait) const;
  virtual void            copy_sequence_block   (sex_t SEX, unsigned int strand, unsigned int from_locus,
      unsigned int to_locus, const TTQuanti *parent);
  virtual void            copy_sequence_1locus  (sex_t SEX, unsigned int strand, unsigned int at,
                                                   const TTQuanti *parent);
};
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuanti_continuous_no_pleio : public TTQuanti_continuous {
 
public:
  TTQuanti_continuous_no_pleio() : TTQuanti_continuous() {}
 
  TTQuanti_continuous_no_pleio(const TTQuanti_continuous_no_pleio& TT) : TTQuanti_continuous(TT) {}
 
  virtual ~TTQuanti_continuous_no_pleio() {}
 
  //TTtrait implementation:
  virtual void init_sequence ();
  virtual TTQuanti_continuous_no_pleio*      clone ()     {return new TTQuanti_continuous_no_pleio(*this);}
  virtual void show_up();
 
  //TTQuanti implementation:
  virtual double          get_additive_genotype (const unsigned int trait) const;
  virtual double          get_dominant_genotype (const unsigned int trait) const;
  virtual void            copy_sequence_block   (sex_t SEX, unsigned int strand, unsigned int from_locus,
      unsigned int to_locus, const TTQuanti *parent);
  virtual void            copy_sequence_1locus  (sex_t SEX, unsigned int strand, unsigned int at,
                                                   const TTQuanti *parent);
};
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuanti_continuous_single : public TTQuanti_continuous {
 
public:
  TTQuanti_continuous_single () : TTQuanti_continuous() {}
 
  TTQuanti_continuous_single (const TTQuanti_continuous_single& TT) : TTQuanti_continuous(TT) {}
 
  virtual ~TTQuanti_continuous_single () {}
 
  //TTtrait implementation:
  virtual void init_sequence ();
  virtual TTQuanti_continuous_single*      clone ()     {return new TTQuanti_continuous_single(*this);}
  virtual void show_up();
 
  //TTQuanti implementation:
  virtual double get_additive_genotype (const unsigned int trait) const;
  virtual double get_dominant_genotype (const unsigned int trait) const;
  virtual void   copy_sequence_block   (sex_t SEX, unsigned int strand, unsigned int from_locus,
      unsigned int to_locus, const TTQuanti *parent);
  virtual void            copy_sequence_1locus  (sex_t SEX, unsigned int strand, unsigned int at,
                                                   const TTQuanti *parent);
};
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuanti_diallelic :  public TTQuanti {
 
public:
 
  TTQuanti_diallelic () : TTQuanti(), _sequence(0) {}
 
  TTQuanti_diallelic (const TTQuanti& T) : TTQuanti(T), _sequence(0) {}
 
  virtual ~TTQuanti_diallelic () {reset();}
 
  virtual bool              get_allele_bit (unsigned int position, unsigned int allele) const;
  virtual void              set_allele_bit (unsigned int position, unsigned int allele, bool value);
  void                      inherit_free (const TTrait* mother, const TTrait* father);
 
  //implements TTrait:
  virtual   void            reset ();
  virtual   void            init  ();
  virtual   void            set_sequence (void** seq);
  virtual   void**          get_sequence () const                             {return (void**)_sequence;}
  virtual   unsigned int    get_allele (int loc, int all) const ;
  virtual   double          get_allele_value (int locus, int allele) const;
  virtual   void            set_allele_value (unsigned int locus, unsigned int allele, double value);
 
  virtual   TTQuanti_diallelic&       operator= (const TTrait& T);
  virtual   bool            operator== (const TTrait& T);
  virtual   bool            operator!= (const TTrait& T);
 
  //implements StorableComponent:
  virtual   void            store_data    ( BinaryStorageBuffer* saver  );
  virtual   bool            retrieve_data ( BinaryStorageBuffer* reader );
 
  //implements TTQuanti:
  virtual double            get_full_genotype     (unsigned int trait);
  virtual void              set_allele   (int locus, int allele, double value) {_sequence[allele][locus] = bool(value);}
  virtual void              mutate_add            (unsigned int position, unsigned int allele, double value)
                               {_sequence[allele][position] = bool(value);}
  virtual void              mutate_inplace        (unsigned int position, unsigned int allele, double value)
                               {_sequence[allele][position] ^= 1;}
 
protected:
 
  unsigned char** _sequence;
};
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuanti_diallelic_no_pleio : public TTQuanti_diallelic {
 
public:
  TTQuanti_diallelic_no_pleio() : TTQuanti_diallelic() {}
 
  TTQuanti_diallelic_no_pleio(const TTQuanti_diallelic_no_pleio& TT) : TTQuanti_diallelic(TT) {}
 
  virtual ~TTQuanti_diallelic_no_pleio() {}
 
  //TTtrait implementation:
  virtual void init_sequence ();
  virtual TTQuanti_diallelic_no_pleio*      clone ()     {return new TTQuanti_diallelic_no_pleio(*this);}
  virtual void show_up();
 
  //TTQuanti implementation:
  virtual double get_additive_genotype (const unsigned int trait) const;
  virtual double get_dominant_genotype (const unsigned int trait) const;
  virtual void   copy_sequence_block   (sex_t SEX, unsigned int strand, unsigned int from_locus,
      unsigned int to_locus, const TTQuanti *parent);
  virtual void   copy_sequence_1locus  (sex_t SEX, unsigned int strand, unsigned int at,
                                                   const TTQuanti *parent);
};
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuanti_diallelic_full_pleio : public TTQuanti_diallelic {
 
public:
  TTQuanti_diallelic_full_pleio () : TTQuanti_diallelic() {}
 
  TTQuanti_diallelic_full_pleio (const TTQuanti_diallelic_full_pleio& TT) : TTQuanti_diallelic(TT) {}
 
  virtual ~TTQuanti_diallelic_full_pleio () {}
 
  //TTtrait implementation:
  virtual void init_sequence ();
  virtual TTQuanti_diallelic_full_pleio*      clone ()     {return new TTQuanti_diallelic_full_pleio(*this);}
  virtual void show_up();
 
  //TTQuanti implementation:
  virtual double get_additive_genotype (const unsigned int trait) const;
  virtual double get_dominant_genotype (const unsigned int trait) const;
  virtual void   copy_sequence_block   (sex_t SEX, unsigned int strand, unsigned int from_locus,
      unsigned int to_locus, const TTQuanti *parent);
  virtual void   copy_sequence_1locus  (sex_t SEX, unsigned int strand, unsigned int at,
                                                   const TTQuanti *parent);
};
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuanti_diallelic_var_pleio : public TTQuanti_diallelic {
 
public:
  TTQuanti_diallelic_var_pleio () : TTQuanti_diallelic() {}
 
  TTQuanti_diallelic_var_pleio (const TTQuanti_diallelic_var_pleio& TT) : TTQuanti_diallelic(TT) {}
 
  virtual ~TTQuanti_diallelic_var_pleio () {}
 
  unsigned int get_sequence_block_size (unsigned int from, unsigned int to);
 
  //TTtrait implementation:
  virtual void init_sequence ();
  virtual TTQuanti_diallelic_var_pleio*      clone ()     {return new TTQuanti_diallelic_var_pleio(*this);}
  virtual void show_up();
 
  //TTQuanti implementation:
  virtual double get_additive_genotype (const unsigned int trait) const;
  virtual double get_dominant_genotype (const unsigned int trait) const;
  virtual void   copy_sequence_block   (sex_t SEX, unsigned int strand, unsigned int from_locus,
      unsigned int to_locus, const TTQuanti *parent);
 
  virtual void   copy_sequence_1locus  (sex_t SEX, unsigned int strand, unsigned int at,
                                                   const TTQuanti *parent);
};
 
 
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TProtoQuanti : public TTProtoWithMap {
 
public:
  
  TProtoQuanti ();
  TProtoQuanti (const TProtoQuanti& T);
  virtual ~TProtoQuanti ();
  
  unsigned int              get_num_traits()                    {return _num_traits;}
  unsigned int              get_num_locus()                     {return _num_locus;}
  unsigned int              get_num_locus (unsigned int trait)  {return _trait_table[trait].size();}
  unsigned int              get_pleiotropy_type()               {return _pleio_type;}
  unsigned int              get_seq_length ()                   {return _seq_length;}
  size_t                    get_size_locus_type()               {return _sizeofLocusType;}
  size_t                    get_locus_byte_size()               {return _locusByteSize;}
  vector<double>            get_env_var ()                      {return _eVariance;}
  vector<double>            get_heritability ()                 {return _h2;}
  unsigned int              get_h2_setTime   ()                 {return _h2_setTime;}
  bool                      get_h2_isBroad   ()                 {return _h2_isBroad;}
  unsigned int              get_allele_model ()                 {return _allele_model;}
  const TMatrix&            get_diallele_values ()              {return _allele_value;}
  double                    get_diallele_value  (unsigned int locus, unsigned int allele)
                                    {return _allele_value.get(locus, allele);}
  double                    get_seq_diallele_value (unsigned int position, unsigned int allele)
                                    {return _sequence_diallele_values[allele][position];}
  bool                      has_equal_domCoeff()                {return _equal_dom_coeff;}
  bool                      has_equal_effects ()                {return _equal_effects;}
  double                    get_equal_val_0 ()                  {return _equal_val_0;}
  double                    get_equal_val_1 ()                  {return _equal_val_1;}
  double                    get_equal_val_diff ()               {return _equal_val_diff;}
  const bitstring&          get_trait_mask (unsigned int trait) const
                                                                {return _trait_masks[trait];}
  double                    get_init_value (unsigned int i)     {return _init_value.get(0,i); }
  double                    get_init_variance (unsigned int i)  {return _init_variance.get(0,i);}
  unsigned int              get_doInitMutation ()               {return _doInitMutation;}
  TTQuantiSH*               get_stater  ()                      {return _stats;}
 
  const TMatrix&            get_pleio_matrix ()                 {return _pleio_matx;}
  unsigned int              get_locus_seq_pos (unsigned int loc, unsigned int trait)
                                                                {return _trait_table[trait][loc];}
  unsigned int              get_locus_ID  (unsigned int locus, unsigned int trait)
                                                                {return _trait_locus_table[trait][locus];}
  unsigned int              get_locus_PD        (unsigned int locus)  {return _locus_table[locus][1];}
  unsigned int              get_locus_start_pos (unsigned int locus)  {return _locus_table[locus][0];}
  unsigned int              get_sequence_block_size (unsigned int from, unsigned int to) {
                              return (to < _num_locus ? _locus_table[to][0] : _seq_length)
                                     - _locus_table[from][0];
                            }
 
  void                      set_eVarianceSD                     (unsigned int trait, double SD);
  void                      set_init_values                     (const double *values, unsigned int nval);
  void                      set_trait_value_func_ptr            (bool withVe);
  double                    set_trait_value_VE                  (const TTQuanti* ind, const unsigned int trait);
  double                    set_trait_value_noVE                (const TTQuanti* ind, const unsigned int trait);
 
  double                    set_genotype_value_additive         (const TTQuanti* ind, const unsigned int trait);
  double                    set_genotype_value_dominance        (const TTQuanti* ind, const unsigned int trait);
 
  int                       get_allele_position                 (const unsigned int locus, const unsigned int trait);
  double                    get_allele_value                    (const TTQuanti* ind, const unsigned int allele, const unsigned int locus, const unsigned int trait);
  double                    get_genotype_with_dominance         (const double a1, const double a2, const unsigned int locus, const unsigned int trait);
  double                    get_genotype_dominance_h            (double a1, double a2, double h);
  double                    get_genotype_dominance_k            (double a1, double a2, double k);
  unsigned int              get_dominance_model                 () {return _dominance_model;}
  double                    get_dominance                       (unsigned int locus, unsigned int trait)
                                                                   {return _dominance_effects.get(trait, locus); }
  const TMatrix&            get_dominance_effects               () {return _dominance_effects;}
 
  double                    get_trait_mutation_variance         (unsigned int trait);
 
  double                    get_genotypic_value                 (const TTQuanti* ind, const unsigned int trait) {
     return (this->*_set_genotype_func_ptr)(ind, trait);
   }
 
  double                    get_phenotypic_value                 (const TTQuanti* ind, const unsigned int trait) {
     return (this->*_set_trait_value_func_ptr)(ind, trait);
  }
 
 
  // Parameter setters
  bool                      setHeritabilityParams             ();
  bool                      setInitialValuesParams            ();
  bool                      setGeneticMapParams               ();
  unsigned int              setAlleleModel                    ();
  bool                      setMutationModel_no_pleio         ();
  bool                      setMutationModel_full_pleio       ();
  bool                      setMutationModel_var_pleio        ();
  void                      setTraitAndLocusTables_full_pleio ();
  void                      setTraitAndLocusTables_no_pleio   (TMatrix& mat);
  bool                      setMutationCorrelation            ();
  bool                      setDominanceParameters            ();
 
  bool                      setDiallelicMutationModel         ();
 
  bool                      setMutationSigmaFromQuantiMutationVariance   ();
  bool                      setMutationSigmaFromQuantiMutationVariance_no_pleio ();
  bool                      readMatrixFromQuantiMutationMatrix           (vector<vector<double>>& varmat);
  bool                      setContinuousMutationModel_no_pleio          ();
  bool                      setContinuousMutationModel_full_pleio        ();
  bool                      setContinuousMutationModel_var_pleio         ();
  void                      allocate_gsl_mutation_matrix_space           (unsigned int num_locus);
  void                      deallocate_gsl_mutation_matrix_space         ();
  gsl_matrix*               set_gsl_mutation_matrix_from_sigma           (unsigned int loc, unsigned int pleio_deg);
  gsl_matrix*               set_gsl_mutation_matrix                      (unsigned int pleio_deg, const vector<double>& varcov);
  void                      set_mutation_matrix_decomposition            (unsigned int loc, unsigned int pleio_deg);
  
  // mutation generation functions
  // first set for continuum-of-allele model with a single mutation matrix across all loci:
  double*                   getMutationEffectMultivariateGaussian (unsigned int loc);
  double*                   getMutationEffectBivariateGaussian    (unsigned int loc);
  double*                   getMutationEffectUnivariateGaussian   (unsigned int loc);
  // different implementation for locus-specific mutation matrices:
  double*                   getMutationEffectMultivariateGaussianLocSpec (unsigned int loc);
  double*                   getMutationEffectBivariateGaussianLocSpec    (unsigned int loc);
  double*                   getMutationEffectUnivariateGaussianLocSpec   (unsigned int loc);
 
//  double*                   getMutationEffectUnivariateDiallelic  (unsigned int loc);
  double*                   getMutationEffectBivariateDiallelic   (unsigned int loc);
  double*                   getMutationEffects                    (unsigned int loc) {
    return (this->* _getMutationValues)(loc);
  }
  double*                   getMutationEffectsVarPleio            (unsigned int loc) {
    return (this->* _getMutationValuesVarPleio[loc]) (loc);
  }
  void                      inherit                               (sex_t SEX, TTQuanti* ind, const TTQuanti* parent);
  void                      inherit_free                          (sex_t SEX, TTQuanti* ind, const TTQuanti* parent);
  void                      inherit_low                           (sex_t SEX, TTQuanti* ind, const TTQuanti* parent);
 
  unsigned int              get_num_mutations                     () {
    return (unsigned int)RAND::Binomial(_mutation_rate, _2L);
  }
  void                      mutate                                (TTQuanti* ind);
  void                      mutate_nill                           (TTQuanti* ind);
  void                      mutate_full_pleio                     (TTQuanti* ind);
  void                      mutate_var_pleio                      (TTQuanti* ind);
  void                      mutate_no_pleio                       (TTQuanti* ind);
  void                      mutate_inplace_full_pleio             (TTQuanti* ind);
  void                      mutate_inplace_var_pleio              (TTQuanti* ind);
  void                      mutate_inplace_no_pleio               (TTQuanti* ind);
 
  void                      mutate_diallelic_pleio                (TTQuanti* ind);
  void                      mutate_diallelic_var_pleio            (TTQuanti* ind);
  void                      mutate_diallelic_no_pleio             (TTQuanti* ind);
    
  // EPISTASIS
  bool                      setEpistasisParameters               ();
//  double                    get_epistatic_mean                   () {return _epistatic_mean;}
//  double                    get_epistatic_variance               () {return _epistatic_variance;}
//  double                    get_epistatic_prop_large             () {return _epistatic_prop_large;}
//  double                    get_epistatic_large_mean             () {return _epistatic_large_mean;}
//  double                    get_epistatic_large_variance         () {return _epistatic_large_variance;}
  unsigned int              get_num_epi_coefs                    () {return _num_epi_coefs;}
  bool                      do_epistasis                         () {return _epistasis;}
  const TMatrix&            get_epi_coefs                        ()  const  {return _epistatic_coefs_matrix;}
  const TMatrix&            get_epi_coef_index                   ()  const  {return _epistatic_coefs_indices;}
 
  virtual   void            reset ()                            {TTProtoWithMap::reset();}
  virtual   TTQuanti*       hatch ();
  virtual   TProtoQuanti*   clone ()                            {return new TProtoQuanti(*this);}
  virtual   trait_t         get_type () const                   {return QUANT;}
  virtual   int                     get_phenotype_dimension () {return _num_traits;}
  virtual   int                     get_allele_number       () {return (_allele_model > 2 ? 3: 2);} //3 just means more than di-allelic
  virtual   int                     get_locus_number        () {return _num_locus;}
  virtual   void            store_data    ( BinaryStorageBuffer* saver  )
    {saver->store(&_seq_length,sizeof(int));}
  
  virtual   bool            retrieve_data ( BinaryStorageBuffer* reader )
    {reader->read(&_seq_length,sizeof(int));return true;}
  virtual   bool            setParameters();
  virtual   void            loadFileServices ( FileServices* loader );
  virtual   void            loadStatServices ( StatServices* loader );
  virtual   bool            resetParameterFromSource (std::string param, SimComponent* cmpt) {return false;}
 
protected:
  
  unsigned int    _num_locus;
  unsigned int    _2L;
  unsigned int    _num_traits;
  unsigned int    _seq_length;
  unsigned int    _allele_model;
  string          _diallele_datatype;
 
  //mutations:
  double          _mutation_rate;
  double          _genomic_mutation_rate;
 
  //Pleiotropy type: 0 = "no pleio", 1 = "full", 2 = "variable"
  unsigned int    _pleio_type;
 
private:
 
  TMatrix         _allele_value;  // num locus x 2
  double*         _sequence_diallele_values[2];
  bool            _equal_effects;    // true when all loci share the same (val_0, val_1)
  double          _equal_val_0;      // shared ancestral allele value
  double          _equal_val_1;      // shared derived allele value
  double          _equal_val_diff;   // val_1 - val_0, precomputed
  vector<bitstring> _trait_masks;    // one per trait, bit i set iff position i belongs to that trait
  TMatrix         _init_value;    // initial value of each trait
  TMatrix         _init_variance; // initial variance of each trait
//  double*         _phenotypes;   // num traits
  
  //mutations:
  TMatrix         _mutation_correlation;  // 1 x n_locus, one correlation per locus, applied to all trait pairs
  TMatrix         _mutation_sigma;        // n_locus x n_trait, variance per locus, used in the bivariate case, and to initialize the m-matrix
  gsl_matrix      **_gsl_mutation_matrix; // one mutation matrix per locus
  gsl_matrix      **_evect;
  gsl_vector      **_eval;
  gsl_vector      **_effects_multivar;
  gsl_vector      **_ws;
  double          _effects_bivar[2];
  unsigned int    _doInitMutation;
  bool            _mutationVarianceIsLocusSpecific; // a single mutation matrix for all loci
  bool            _mutationEffectIsFixedDiAllele;   // bitstring sequence not storing effects
 
  // dominance
  unsigned int    _dominance_model;
  TMatrix         _dominance_effects;
  bool            _equal_dom_coeff;
  
  //recombination:
  bool*           _all_chooser;
  size_t          _locusByteSize;
  size_t          _sizeofLocusType;
  
  //heritability
  vector<double>  _eVariance;
  vector<double>  _h2;
  unsigned int    _h2_setTime;
  bool            _h2_isBroad;
 
  //variable pleiotropy
 
  TMatrix                        _pleio_matx;
 
  vector< vector<unsigned int> > _trait_table;
 
  vector< vector<unsigned int> > _trait_locus_table;
 
  vector< vector<unsigned int> > _locus_table;
 
 
  // EPISTASIS
  bool                      _epistasis;
  unsigned int              _num_epi_coefs;
  TMatrix                   _epistatic_coefs_matrix;
  TMatrix                   _epistatic_coefs_indices;
 
 
  // function pointers
 
  void    (TProtoQuanti::* _inherit_fun_ptr)   (sex_t, TTQuanti*, const TTQuanti*);
 
  void    (TProtoQuanti::* _mutation_func_ptr) (TTQuanti*);
 
  double* (TProtoQuanti::* _getMutationValues) (unsigned int);
 
  vector< double* (TProtoQuanti::* ) (unsigned int) > _getMutationValuesVarPleio;
 
  double  (TProtoQuanti::* _getGenotypeWithDominance) (double, double, double);
 
 
  double  (TProtoQuanti::* _set_trait_value_func_ptr) (const TTQuanti*, const unsigned int);
 
 
  double  (TProtoQuanti::* _set_genotype_func_ptr)    (const TTQuanti*, const unsigned int);
  
  friend class TTQuanti;
  friend class TTQuanti_continuous;
  friend class TTQuanti_continuous_full_pleio;
  friend class TTQuanti_continuous_var_pleio;
  friend class TTQuanti_continuous_no_pleio;
  friend class TTQuanti_continuous_single;
  
  TTQuantiSH* _stats;
  TTQuantiFH* _writer;
  TTQuantiFH* _reader;
  TTQFreqExtractor* _freqExtractor;
  TTQOhtaStats*     _ohtaStats;
};
 
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuantiSH : public TraitStatHandler<TProtoQuanti, TTQuantiSH> {
 
  double *_meanP, *_meanG, *_Va, *_Vg, *_Vb, *_Vp, *_covar;
  double **_pmeanP, **_pmeanG, **_pVa, **_pVp, **_pcovar, **_peigval, **_peigvect;
  bool _eVar;
 
  unsigned int _num_locus, _num_trait, _patchNbr;
 
  gsl_matrix *_G, *_evec;
  gsl_vector *_eval;
  gsl_eigen_symmv_workspace *_ws;
  
  DataTable< double > _phenoTable, _genoTable;
  unsigned int _table_set_gen, _table_set_age, _table_set_repl;
 
public:
 
    TTQuantiSH(TProtoQuanti* TP) 
    : TraitStatHandler<TProtoQuanti, TTQuantiSH> (TP), 
    _meanP(0), _meanG(0), _Va(0), _Vg(0), _Vb(0), _Vp(0), _covar(0),
    _pmeanP(0), _pmeanG(0), _pVa(0), _pVp(0), _pcovar(0), _peigval(0), _peigvect(0),
    _eVar(0), _num_locus(0), _num_trait(0),_patchNbr(0), _G(0), _evec(0),_eval(0),_ws(0),
    _table_set_gen(999999), _table_set_age(999999), _table_set_repl(999999)
    {}
  
  virtual ~TTQuantiSH() {resetPtrs();}
  
  void  resetPtrs();
  
  virtual void init ( );
  
  virtual bool setStatRecorders    (std::string& token);
 
  void addQuanti                   (age_t AGE);
  void addQuantiPerPatch           (age_t AGE);
  void addAvgPerPatch              (age_t AGE);
  void addGenotPerPatch            (age_t AGE);
  void addVarPerPatch              (age_t AGE);
  void addCovarPerPatch            (age_t AGE);
  void addEigenPerPatch            (age_t AGE);
  void addEigenValuesPerPatch      (age_t AGE);
  void addEigenVect1PerPatch       (age_t AGE);
  void addSkewPerPatch             (age_t AGE);
 
  void   setDataTables             (age_t AGE);
  void   setAdultStats             ( ) {setStats(ADULTS);}
  void   setOffsprgStats           ( ) {setStats(OFFSPRG);}
  void   setStats                  (age_t AGE);
  double getMeanGenot              (unsigned int i) {return _meanG[i];}
  double getMeanPhenot             (unsigned int i) {return _meanP[i];}
  double getVa                     (unsigned int i) {return _Va[i];}
  double getVg                     (unsigned int i) {return _Vg[i];}
  double getVb                     (unsigned int i) {return _Vb[i];}
  double getVp                     (unsigned int i) {return _Vp[i];}
  double getQst                    (unsigned int i) {return _Vb[i]/(_Vb[i]+2*_Va[i]);}
  double getCovar                  (unsigned int i) {return _covar[i];}
 
  double getMeanGenotPerPatch      (unsigned int i, unsigned int p) {return _pmeanG[i][p];}
  double getMeanPhenotPerPatch     (unsigned int i, unsigned int p) {return _pmeanP[i][p];}
  double getVaPerPatch             (unsigned int i, unsigned int p) {return _pVa[i][p];}
  double getVpPerPatch             (unsigned int i, unsigned int p) {return _pVp[i][p];}
  double getEigenValuePerPatch     (unsigned int i, unsigned int p) {return _peigval[i][p];}
  double getCovarPerPatch          (unsigned int p, unsigned int i) {return _pcovar[p][i];}
  double getEigenVectorEltPerPatch (unsigned int p, unsigned int v) {return _peigvect[p][v];}
  double getSkewPerPatch           (unsigned int i, unsigned int p);
 
  vector<double> getSNPalleleFreqInPatch  (Patch* patch, const age_idx AGE);
  vector<double> getVaWithDominance (Patch* curPop, const age_idx AGE);
  vector<double> getVaNoDominance   (Patch* curPop, const age_idx AGE);
 
};
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQuantiFH : public TraitFileHandler<TProtoQuanti> {
 
  string _output_option;
  bool   _has_genetic_map;
  
public:
  TTQuantiFH(TProtoQuanti* T) : TraitFileHandler<TProtoQuanti>(T,".quanti"),
    _has_genetic_map(0) {}
  virtual ~TTQuantiFH(){}
  
  void setOutputOption (string opt);
  
  virtual void  FHwrite ();
  void write_TABLE ();
  void print(ofstream& FH, Patch* patch, sex_t SEX, age_idx Ax, unsigned int print_gene, bool print_genotype, bool print_additive_genotype);
  void write_PLINK ();
  void print_PLINK_PED(ofstream& FH, age_idx Ax, Patch *patch);
  void print_PLINK_FAM(ofstream& FH, age_idx Ax, Patch *patch);
 
  virtual void FHread (string& filename);
};
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQFreqExtractor : public TraitFileHandler<TProtoQuanti> {
 
  vector< string > _records;
 
public:
  TTQFreqExtractor(TProtoQuanti* T) : TraitFileHandler<TProtoQuanti>(T,".qfreq") {}
 
  void resetTable();
 
  virtual ~TTQFreqExtractor () {}
  virtual void FHwrite ();
  virtual void FHread (string& filename) {}
};
 
// ------------------------------------------------------------------------------
// ------------------------------------------------------------------------------
class TTQOhtaStats : public TraitFileHandler<TProtoQuanti> {
 
  TMatrix _pairwiseCombs;
 
public:
  TTQOhtaStats(TProtoQuanti* T) : TraitFileHandler<TProtoQuanti>(T,".qohta") {}
 
  virtual ~TTQOhtaStats () {}
  virtual void FHwrite ();
  virtual void FHread (string& filename) {}
};
#endif