src/tree/parameter.h

// $Id: parameter.h,v 1.45 2011/10/31 22:04:51 ewalkup Exp $

/*
  Copyright 2002  Peter Beerli, Mary Kuhner, Jon Yamato and Joseph Felsenstein

  This software is distributed free of charge for non-commercial use
  and is copyrighted.  Of course, we do not guarantee that the software
  works, and are not responsible for any damage you may cause or have.
*/

/*************************************************************************
   The Parameter class and its helper ResultStruct are used to store
   information about individual parameters of the force model.  For
   example, a single Parameter object might represent the migration
   rate from population 2 to population 1.

   Parameters contain information needed to describe a parameter in
   input/output, and also store information about results for that
   parameter, such as its MLE, profiles, and perhaps someday plots.

   Because everything about a parameter except its results is fixed
   after the user-input phase, only the results have set methods

   A Parameter has an IsValidParameter() method establishing its validity.
   No other field in an invalid Parameter should be used!  Invalid
   Parameters are used as placeholders, for example representing the
   non-existent "migration rate from 1 to 1".

   **

  The ParamVector class is a structure that holds a "checked out" set
  of all Parameters.  This is done so that parts of the program which
  need rapid, read-write access to all Parameters can get it.  One
  creates a ParamVector to "check out" the Parameters and destroys it
  (or lets it go out of scope) to automatically "check in" your changes.
  It is a an error, and will throw a runtime exception, to check out
  a second ParamVector while the first is still out.

  If you want only read access to a ParamVector, construct one with an
  argument of "true" (meaning read-only) (and make it const to keep you
  honest).  This ParamVector will not attempt a checkin.  It is fine to
  have as many read-only ParamVectors as you like.  However, bear in mind
  that if you check out a writable one and write to it, the read-only ones
  will present a stale view (they are never updated).

  NB:  There may be invalid Parameters in a ParamVector (such as
  the nonexistent migration rates along the diagonal).  It is a fatal
  mistake to try to extract results, names, etc. from such a Parameter.

Written by Mary Kuhner October 1 2001
****************************************************************************/

#ifndef PARAMETER_H
#define PARAMETER_H

#include <cassert>
#include <vector>
#include <string>
#include <map>

#include "plotstat.h"
#include "prior.h"
#include "constants.h"
#include "defaults.h"
#include "types.h"
#include "paramstat.h"

//------------------------------------------------------------------------------------

class ForceSummary;
class Random;

//------------------------------------------------------------------------------------

struct ResultStruct
{
  private:
    DoubleVec1d mles;                // dim: regions
    DoubleVec1d overallmle;          // vector of one element
    vector<ProfileStruct> profiles;        // dim: regions
    vector<ProfileStruct> overallprofile;  // vector of one element

  public:

    // Creation and Destruction
    // ResultStruct();                         // we accept the default
    // ResultStruct(const ResultStruct& src);  // we accept the default
    // ~ResultStruct();                        // we accept the default
    // ResultStruct& operator=(const ResultStruct& src); // we accept the default

    // Getters (mostly not inline due to error checking code)
    double GetMLE(long region) const;
    DoubleVec1d GetRegionMLEs() const { return mles; };
    double GetOverallMLE() const;
    DoubleVec1d GetAllMLEs() const;

    const ProfileStruct& GetProfile(long region) const;
    const vector<ProfileStruct>& GetRegionProfiles() const { return profiles; };
    const ProfileStruct& GetOverallProfile() const;
    vector<ProfileStruct> GetAllProfiles() const;
    const ProfileLineStruct& GetProfileFor(double centile, long reg) const;

    // Setters
    void AddMLE(double mle, long region);
    void AddOverallMLE(double mle);

    void AddProfile(const ProfileStruct& profile)
    { profiles.push_back(profile); };
    void AddOverallProfile(const ProfileStruct& profile)
    { overallprofile.push_back(profile); };

}; // ResultStruct definition

//------------------------------------------------------------------------------------

class Parameter
// EWCOMMENT -- keep an eye on how much "special case" code is added
// to this class. If we get a lot, then we need to think about
// using subclasses.  We don't do it now because we stuff these
// into vectors.
{
  private:
    ParamStatus m_status;               // if this is "invalid" nothing else matters
    unsigned long m_paramvecIndex;
    string m_shortname;
    string m_name;
    force_type m_force;                 // parameter knows to which force it belongs

    ResultStruct m_results;             // used to be "mutable", but no longer needed.
    method_type m_method;               // method used to calculate this parameter
    proftype m_profiletype;
    double m_truevalue;                 // true value of variable

    // bayesian variables
    Prior  m_prior;                     // a 'Prior' object (can be linear, log, ...?)

    Parameter();                        // not defined

  public:

    // Creation and destruction
    Parameter(const ParamStatus& status,
              unsigned long paramIndex,
              const string sname,
              const string lname,
              force_type thisforce,
              method_type meth,
              proftype prof,
              const UIVarsPrior & uiprior,
              double truevalue
              );

    //Use only to create invalid parameters.
    Parameter(const ParamStatus& status, unsigned long paramIndex);

    // We accept default copy ctor, op=, dtor

    // Getters
    ParamStatus GetStatus()   const { return m_status; };
    unsigned long GetParamVecIndex()   const { return m_paramvecIndex; };
    bool IsValidParameter()   const { return (m_status.Valid()); };
    bool IsVariable()         const { return (m_status.Inferred()); };
    bool IsForce(force_type thistag)  const { return (m_force == thistag); };
    bool IsEasilyBayesianRearrangeable() const;
    force_type WhichForce()       const { return m_force; };
    string GetShortName()     const { assert (IsValidParameter()); return m_shortname; };
    string GetName()          const { assert (IsValidParameter()); return m_name; };
    string GetUserName()      const;
    proftype GetProfileType() const { return m_profiletype; };
    bool IsProfiled()         const { return (m_profiletype != profile_NONE); };
    method_type GetMethod()   const { return m_method; };
    double GetTruth()         const { return m_truevalue; };

    // The Prior interface (for Bayesian analysis)
    Prior      GetPrior()         const { return m_prior; };

    void       SetShortName(string n) { m_shortname = n; };
    void       SetName(string n) { m_name = n; };
    std::pair<double, double> DrawFromPrior() const;
    bool IsZeroTrueMin();

    // MLEs
    double GetMLE(long region) const { assert (IsValidParameter()); return m_results.GetMLE(region); };
    DoubleVec1d GetRegionMLEs() const { assert (IsValidParameter()); return m_results.GetRegionMLEs(); };
    double GetOverallMLE() const { assert (IsValidParameter()); return m_results.GetOverallMLE(); };
    DoubleVec1d GetAllMLEs() const { assert (IsValidParameter()); return m_results.GetAllMLEs(); };

    // Profiles
    vector<vector<centilepair> > GetProfiles(long region) const;
    vector<vector<centilepair> > GetOverallProfile() const;
    vector<centilepair> GetPriorLikes(long region) const;
    vector<centilepair> GetOverallPriorLikes() const;
    vector<centilepair> GetCIs(long region) const;
    vector<centilepair> GetOverallCIs() const;

    bool CentileIsExtremeLow (double centile, long reg) const;
    bool CentileIsExtremeHigh(double centile, long reg) const;
    bool CentileHadWarning(double centile, long reg) const;

    // Setters
    //   no setters for most things outside of m_results, since everything
    //   is set in the constructor.  The m_results setters can be const
    //   since m_results is mutable.  (The prior can be set, but in that
    //   case the Parameter must be non-const).

    // MLEs
    void AddMLE(double mle, long region) { assert (IsValidParameter()); m_results.AddMLE(mle, region); };
    void AddOverallMLE(double mle) { assert (IsValidParameter()); m_results.AddOverallMLE(mle); };

    // Profiles
    void AddProfile(const ProfileStruct& prof, likelihoodtype like);
};

class ParamVector
{
  private:

    // This class is a lock.  NO COPYING ALLOWED.  Don't define these.
    ParamVector(const ParamVector&);              // not defined
    ParamVector& operator=(const ParamVector&);   // not defined

    static bool s_locked;
    bool m_readonly;
    ForceSummary& forcesum;      // to allow check-in in destructor

    // not very private, as we hand out references....
    vector<Parameter> parameters;

  public:

    // The constructor checks out the Parameters and the destructor checks them
    // back in.  If you check out a read-only version, you will do yourself a
    // favor if you also make it const (i.e. 'const ParamVector pvec(true)' )
    ParamVector(bool rd);
    ~ParamVector();

    // vector emulation routines
    typedef vector<Parameter>::iterator iterator;
    typedef vector<Parameter>::const_iterator const_iterator;
    Parameter&     operator[](long index);
    const Parameter&     operator[](long index) const;
    unsigned long  size()  const { return parameters.size(); };
    bool           empty() const { return parameters.empty(); };
    iterator       begin()       { return parameters.begin(); };
    const_iterator begin() const { return parameters.begin(); };
    iterator       end()         { return parameters.end(); };
    const_iterator end()   const { return parameters.end(); };

    paramlistcondition CheckCalcProfiles() const;
    paramlistcondition CheckCalcPProfiles() const;
    long           NumProfiledParameters() const;
    long           NumVariableParameters() const;
    long           ChooseSampleParameterIndex(Random * randomSource) const;
};

#endif // PARAMETER_H

//____________________________________________________________________________________