blast/gumbel_params/njn_localmaxstatutil.cpp

/* $Id: njn_localmaxstatutil.cpp 462397 2015-03-18 16:23:53Z vasilche $
* ===========================================================================
*
*                            PUBLIC DOMAIN NOTICE
*               National Center for Biotechnology Information
*
*  This software/database is a "United States Government Work" under the
*  terms of the United States Copyright Act.  It was written as part of
*  the author's offical duties as a United States Government employee and
*  thus cannot be copyrighted.  This software/database is freely available
*  to the public for use. The National Library of Medicine and the U.S.
*  Government have not placed any restriction on its use or reproduction.
*
*  Although all reasonable efforts have been taken to ensure the accuracy
*  and reliability of the software and data, the NLM and the U.S.
*  Government do not and cannot warrant the performance or results that
*  may be obtained by using this software or data. The NLM and the U.S.
*  Government disclaim all warranties, express or implied, including
*  warranties of performance, merchantability or fitness for any particular
*  purpose.
*
*  Please cite the author in any work or product based on this material.
*
* ===========================================================================*/

/*****************************************************************************

File name: njn_localmaxstatutil.cpp

Author: John Spouge

Contents:

******************************************************************************/

#include <ncbi_pch.hpp>

#include "njn_localmaxstatutil.hpp"
#include "njn_approx.hpp"
#include "njn_dynprogproblim.hpp"
#include "njn_integer.hpp"
#include "njn_memutil.hpp"
#include "njn_root.hpp"

#include "sls_alp_data.hpp"

// Kludge: limit optimization by ICC 10.x to avoid undesired references to
// __svml_exp2 (problematic to use from DLLs on x86_64 or at all on ia32).
#if defined(NCBI_COMPILER_ICC)  &&  defined(__OPTIMIZE__) \
    &&  NCBI_COMPILER_VERSION >= 1000  &&  NCBI_COMPILER_VERSION < 1100
#  define NEED_ICC_OPTIMIZATION_LIMITS 1
#endif

USING_NCBI_SCOPE;
USING_SCOPE(blast);
USING_SCOPE(Njn);

void LocalMaxStatUtil::flatten ( // allocates memory for linear probabilities and scores
size_t dimension_, // dimension of equilProb_
const Int4 *const *scoreMatrix_, // packed scoring matrix [0...dimension_)[0...dimension2_)
const double *const *prob_, // prob_ [0...dimension_)[0...dimension2_) : distribution of scores sum to 1.0
size_t *dim_, // dimension of p_
Int4 **score_, // score [0...dim_) in increasing order
double **p_, // linear p_ [0...dim_) : distribution of scores
size_t dimension2_) // dimension2 of equilProb_
{
    if (dimension2_ == 0) dimension2_ = dimension_;

    size_t i = 0;
    size_t j = 0;

    double sum = 0.0;

    for (i = 0; i < dimension_; i++)
    {
      for (j = 0; j < dimension2_; j++)
      {
         sum += prob_ [i][j];
      }
    }

#ifdef _DEBUG /* to eliminate compiler warning in release mode */
    const double FUDGE = 20.0;
#endif
    _ASSERT(Approx::relApprox (sum, 1.0, FUDGE * REL_TOL));

    Int4 s = 0;
    Int4 min = scoreMatrix_[0][0];
    Int4 max = scoreMatrix_[0][0];

    for (i = 0; i < dimension_; i++)
    {
      for (j = 0; j < dimension2_; j++)
      {
         if (scoreMatrix_ [i][j] < min)
             min = scoreMatrix_ [i][j];
         else if (max < scoreMatrix_ [i][j])
             max = scoreMatrix_ [i][j];
      }
    }

    assert (min <= max);

    size_t dim = static_cast <size_t> (max - min + 1);
    double *p = new double [dim];
    for (i = 0; i < dim; i++) p [i] = 0.0;

    for (i = 0; i < dimension_; i++)
    {
      for (j = 0; j < dimension2_; j++)
      {
         p [scoreMatrix_ [i][j] - min] += prob_ [i][j];
      }
    }

    *dim_ = 0;

    for (s = min; s <= max; s++)
    {
      if (0.0 < p [s - min]) ++*dim_;
    }

    *p_ = new double [*dim_];
    *score_ = new Int4 [*dim_];
    *dim_ = 0;

    for (s = min; s <= max; s++)
    {
      if (0.0 < p [s - min]) {
         (*score_) [*dim_] = s;
         (*p_) [*dim_] = p [s - min];
         ++*dim_;
      }
    }

    delete [] p; p = 0;
}

double LocalMaxStatUtil::lambda (
size_t dimension_, // dimension of equilProb_
const Int4 *const *scoreMatrix_, // packed scoring matrix [0...dimension_)[0...dimension_)
const double *q_) // q_ [0...dimension_) : distribution of independent letters
{
    size_t i = 0;
    size_t j = 0;

   double **prob = MemUtil::newMatrix <double> (dimension_, dimension_);

    for (i = 0; i < dimension_; i++)
    {
        for (j = 0; j < dimension_; j++)
        {
            prob [i][j] = q_ [i] * q_ [j];
        }
    }

    size_t dim = 0;
    Int4 *score = 0;
    double *p = 0;

    flatten (dimension_, scoreMatrix_, prob, &dim, &score, &p);

    MemUtil::deleteMatrix <double> (prob, dimension_, dimension_); prob = 0;

    double lambdaHat = LocalMaxStatUtil::lambda (dim, score, p);

    delete [] p; p = 0;
    delete [] score; score = 0;

    return lambdaHat;
}


   size_t n_dimension = 0; // dimension of matrices
   const Int4 *n_score = 0; // score_ [0...dimension_ - 1]
   const double *n_prob = 0; // prob_ [0...dimension_ - 1]
   Int4 n_morgue = 0; // score_ [0] - 1
   Int4 n_entry = 0; // n_entry = 0 : weak descending ladder epoch ; n_entry = -1 : strict descending ladder epoch

   void n_setParameters (
   size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
   const Int4 *score_, // scores
   const double *prob_, // probabilities
   Int4 entry_ = 0) // entry_ = 0 : weak descending ladder epoch ; entry_ = -1 : strict descending ladder epoch
   {
      n_dimension = dimension_;
      n_score = score_;
      n_prob = prob_;
      n_morgue = score_ [0] - 1;
      n_entry = entry_;
   }

#ifdef NEED_ICC_OPTIMIZATION_LIMITS
#  pragma optimization_level 1
#endif
   double n_totalProbAssoc (double x_)
   {
      double sum = 0.0;
      for (size_t i = 0; i < n_dimension; i++) {
         sum += n_prob [i] * exp (x_ * static_cast <double> (n_score [i]));
      }
      return sum;
   }

   double n_meanPowerAssoc (double x_, Int4 power_ = 1L)
   {
      double sum = 0.0;
      for (size_t i = 0; i < n_dimension; i++) {
          sum += Integer::integerPower (static_cast <double> (n_score [i]), power_) *
            n_prob [i] * exp (x_ * static_cast <double> (n_score [i]));
      }
      return sum;
   }

   double n_meanAssoc (double x_)
   {
      return n_meanPowerAssoc (x_);
   }

   void n_bracket (double *p_, double *q_)
   {
      const double FACTOR = 0.5;
      *p_ = -log (n_prob [n_dimension - 1]) / static_cast <double> (n_score [n_dimension - 1]);
      while (1.0 <= n_totalProbAssoc (*p_)) {
         *p_ *= FACTOR;
      }
      *q_ = *p_ / FACTOR;
   }


double LocalMaxStatUtil::mu (
size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
const Int4 *score_, // scores in increasing order
const double *prob_) // corresponding probabilities
{
   double mu = 0.0;
   for (size_t i = 0; i < dimension_; i++) {
      mu += static_cast <double> (score_ [i]) * prob_ [i];
   }
   return mu;
}

double LocalMaxStatUtil::lambda (
size_t dimension_, // #(distinct values)
const Int4 *score_, // values
const double *prob_) // probability of corresponding value
{
   n_setParameters (dimension_, score_, prob_);

   double p = 0.0;
   double q = 0.0;

   n_bracket (&p, &q);

   return Root::bisection (1.0, n_totalProbAssoc, p, q, REL_TOL * fabs (p - q));
}

double LocalMaxStatUtil::muPowerAssoc (
size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
const Int4 *score_, // scores in increasing order
const double *prob_, // corresponding probabilities
double lambda_, // lambda
Int4 power_) // power
{
   n_setParameters (dimension_, score_, prob_);

   if (lambda_ == 0.0) lambda_ = lambda (dimension_, score_, prob_);

   return n_meanPowerAssoc (lambda_, power_);
}

double LocalMaxStatUtil::muAssoc (
size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
const Int4 *score_, // scores in increasing order
const double *prob_, // corresponding probabilities
double lambda_) // lambda
{
   return muPowerAssoc (dimension_, score_, prob_, lambda_);
}

double LocalMaxStatUtil::thetaMin (
size_t dimension_, // #(distinct values)
const Int4 *score_, // values
const double *prob_, // probability of corresponding value
double lambda_) // lambda
// assumes logarithmic regime
{
   n_setParameters (dimension_, score_, prob_);

   if (lambda_ == 0.0) lambda_ = lambda (dimension_, score_, prob_);

   double p = 0.0;
   double q = 0.0;
   n_bracket (&p, &q);
   return Root::bisection (0.0, n_meanAssoc, 0.0, lambda_, REL_TOL * fabs (p - q));
}

double LocalMaxStatUtil::rMin (
size_t dimension_, // #(distinct values)
const Int4 *score_, // values
const double *prob_, // probability of corresponding value
double lambda_, // lambda
double thetaMin_) // argument of rate
// assumes logarithmic regime
{
   n_setParameters (dimension_, score_, prob_);

   if (thetaMin_ == 0.0) thetaMin_ = thetaMin (dimension_, score_, prob_, lambda_);

   return n_totalProbAssoc (thetaMin_);
}

#ifdef NEED_ICC_OPTIMIZATION_LIMITS
#  pragma optimization_level 1
#endif
double LocalMaxStatUtil::r ( // r (theta)
size_t dimension_, // #(distinct values)
const Int4 *score_, // scores in increasing order
const double *prob_, // probability of corresponding value
double theta_) // argument of rate
// assumes logarithmic regime
{
   double sum = 0.0;
   for (size_t i = 0; i < dimension_; i++) {
      sum += prob_ [i] * exp (theta_ * static_cast <double> (score_ [i]));
   }
   return sum;
}

Int4 LocalMaxStatUtil::delta ( // theta [minus delta] for ungapped sequence comparison
size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
const Int4 *score_) // scores
{
   size_t i = 0;

   Int4 delta = 0;
   for (i = 0; i < dimension_; i++) {
      delta = Integer::euclidAlgorithm <Int4> (delta, score_ [i]);
   }
   return delta;
}

double LocalMaxStatUtil::thetaMinusDelta ( // theta [minus delta] for ungapped sequence comparison
double lambda_, // lambda, the exponential rate for the local maximum
size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
const Int4 *score_) // scores
{
   double del = static_cast <double> (delta (dimension_, score_));
   return (1.0 - exp (-lambda_ * del)) / del;
}


   Int4 n_step (Int4 oldValue_, size_t state_)
   {
      assert (state_ < n_dimension);
      return n_morgue < oldValue_ ? oldValue_ + n_score [state_] : oldValue_;
   }

   Int4 n_bury (Int4 oldValue_, size_t state_)
   {
      assert (state_ < n_dimension);
      return n_entry < oldValue_ ? oldValue_ : n_morgue;
   }


void LocalMaxStatUtil::descendingLadderEpochRepeat (
size_t dimension_, // #(distinct values)
const Int4 *score_, // values
const double *prob_, // probability of corresponding value
double *eSumAlpha_, // expectation (sum [alpha])
double *eOneMinusExpSumAlpha_, // expectation [1.0 - exp (sum [alpha])]
bool isStrict_, // ? is this a strict descending ladder epoch
double lambda_, // lambda for repeats : default is lambda0_ below
size_t endW_, // maximum w plus 1
double *pAlphaW_, // probability {alpha = w} : pAlphaW_ [0, wEnd)
double *eOneMinusExpSumAlphaW_, // expectation [1.0 - exp (sum [alpha]); alpha = w] : eOneMinusExpSumAlphaW_ [0, wEnd)
double lambda0_, // lambda for flattened distribution (avoid recomputation)
double mu0_, // mean of flattened distribution (avoid recomputation)
double muAssoc0_, // mean of associated flattened distribution (avoid recomputation)
double thetaMin0_, // thetaMin of flattened distribution (avoid recomputation)
double rMin0_, // rMin of flattened distribution (avoid recomputation)
double time_, // get time for the dynamic programming computation
bool *terminated_) // ? Was the dynamic programming computation terminated prematurely ?
// assumes logarithmic regime
{
    // Start dynamic programming probability calculation using notation in
    //
    // Mott R. and Tribe R. (1999)
    // J. Computational Biology 6(1):91-112
    //
    // Karlin S. and Taylor H.M.(1981)
    // A Second Course in Stochastic Processes, p. 480
    //
    // Note there is an error in Eq (6.19) there, which is corrected in Eq (6.20)
    //
    // This program uses departure into (-Inf, 0] not (-Inf, 0)

    // avoid recomputation
#ifdef _DEBUG /* to eliminate compiler warning in release mode */
    double mu0 = 0.0 == mu0_ ? mu (dimension_, score_, prob_) : mu0_;
#endif
    _ASSERT(mu0 < 0.0);
    double lambda0 = 0.0 == lambda0_ ? lambda (dimension_, score_, prob_) : lambda0_;
    _ASSERT (0.0 < lambda0);
    if (lambda_ == 0.0) lambda_ = lambda0;
    _ASSERT (0.0 < lambda_);
#ifdef _DEBUG /* to eliminate compiler warning in release mode */
    double muAssoc0 = 0.0 == muAssoc0_ ? muAssoc (dimension_, score_, prob_, lambda0) : muAssoc0_;
#endif
    _ASSERT (0.0 < muAssoc0);
    double thetaMin0 = 0.0 == thetaMin0_ ? thetaMin (dimension_, score_, prob_, lambda0) : thetaMin0_;
    _ASSERT (0.0 < thetaMin0);
    double rMin0 = 0.0 == rMin0_ ? rMin (dimension_, score_, prob_, lambda0, thetaMin0) : rMin0_;
    _ASSERT (0.0 < rMin0 && rMin0 < 1.0);

    const Int4 ITER_MIN = static_cast <Int4> ((log (REL_TOL * (1.0 - rMin0)) / log (rMin0)));
    _ASSERT (0 < ITER_MIN);
    const Int4 ITER = static_cast <Int4> (endW_) < ITER_MIN ? ITER_MIN : static_cast <Int4> (endW_);
    _ASSERT (0 < ITER);
    const Int4 Y_MAX = static_cast <Int4> (-log (REL_TOL) / lambda0);

    Int4 entry = isStrict_ ? -1 : 0;
    n_setParameters (dimension_, score_, prob_, entry);


    double time0 = 0.0;
    double time1 = 0.0;
        if (time_ > 0.0) Sls::alp_data::get_current_time (time0);

    DynProgProbLim dynProgProb (n_step, dimension_, prob_, score_ [0] - 1, Y_MAX);

    if (pAlphaW_) pAlphaW_ [0] = 0.0;
    if (eOneMinusExpSumAlphaW_) eOneMinusExpSumAlphaW_ [0] = 0.0;

    dynProgProb.update (); // iterate random walk

    Int4 value = 0;

    if (eSumAlpha_) *eSumAlpha_ = 0.0;
    if (eOneMinusExpSumAlpha_) *eOneMinusExpSumAlpha_ = 0.0;

    for (size_t w = 1; w < static_cast <size_t> (ITER); w++) {

        if (w < endW_) { // sum pAlphaW_ [w] and eOneMinusExpSumAlphaW_ [w]

             if (pAlphaW_) pAlphaW_ [w] = 0.0;
             if (eOneMinusExpSumAlphaW_) eOneMinusExpSumAlphaW_ [w] = 0.0;

             for (value = score_ [0]; value <= entry; value++) {
                if (pAlphaW_) pAlphaW_ [w] += dynProgProb.getProb (value);
                if (eOneMinusExpSumAlphaW_) eOneMinusExpSumAlphaW_ [w] +=
                                               dynProgProb.getProb (value) *
                                               (1.0 - exp (lambda_ * static_cast <double> (value)));
             }
        }

        for (value = score_ [0]; value <= entry; value++) {
         if (eSumAlpha_) *eSumAlpha_ += dynProgProb.getProb (value) * static_cast <double> (value);
         if (eOneMinusExpSumAlpha_) *eOneMinusExpSumAlpha_ += dynProgProb.getProb (value) *
                                        (1.0 - exp (lambda_ * static_cast <double> (value)));
        }

        dynProgProb.setValueFct (n_bury);
        dynProgProb.update (); // put probability into the morgue

        dynProgProb.setValueFct (n_step);
        dynProgProb.update (); // iterate random walk

        if (time_ > 0.0)
        {
                        Sls::alp_data::get_current_time (time1);
            if (time1 - time0 > time_)
            {
                *terminated_ = true;
                return;
            }
        }

    }

    for (value = score_ [0]; value <= entry; value++) {
      if (eSumAlpha_) *eSumAlpha_ += dynProgProb.getProb (value) * static_cast <double> (value);
      if (eOneMinusExpSumAlpha_) *eOneMinusExpSumAlpha_ += dynProgProb.getProb (value) *
                                     (1.0 - exp (lambda_ * static_cast <double> (value)));
    }

    // check that not too much probability has been omitted
    double prob = 0.0;
    for (value = entry + 1; value < dynProgProb.getValueUpper (); value++) {
      prob += dynProgProb.getProb (value);
    }
    prob += dynProgProb.getProbLost ();

#ifdef _DEBUG /* to eliminate compiler warning in release mode */
    const double FUDGE = 2.0;
#endif
    _ASSERT (prob <= FUDGE * static_cast <double> (dimension_) * REL_TOL);
}

void LocalMaxStatUtil::descendingLadderEpoch (
size_t dimension_, // #(distinct values)
const Int4 *score_, // values
const double *prob_, // probability of corresponding value
double *eSumAlpha_, // expectation (sum [alpha])
double *eOneMinusExpSumAlpha_, // expectation [1.0 - exp (sum [alpha])]
bool isStrict_, // ? is this a strict descending ladder epoch
double lambda0_, // lambda for flattened distribution (avoid recomputation)
double mu0_, // mean of flattened distribution (avoid recomputation)
double muAssoc0_, // mean of associated flattened distribution (avoid recomputation)
double thetaMin0_, // thetaMin of flattened distribution (avoid recomputation)
double rMin0_, // rMin of flattened distribution (avoid recomputation)
double time_, // get time for the dynamic programming computation
bool *terminated_) // ? Was the dynamic programming computation terminated prematurely ?
{
   descendingLadderEpochRepeat (dimension_, score_, prob_,
      eSumAlpha_, eOneMinusExpSumAlpha_, isStrict_, 0.0, 0, 0, 0,
      lambda0_, mu0_, muAssoc0_, thetaMin0_, rMin0_, time_, terminated_);
}

bool LocalMaxStatUtil::isProbDist (
size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
const double *prob_) // corresponding probabilities
{
   double sum = 0.0;
   for (size_t i = 0; i < dimension_; i++) {
      if (prob_ [i] < 0.0 || 1.0 < prob_ [i]) return false;
      sum += prob_ [i];
   }
   return Approx::relApprox (sum, 1.0, REL_TOL);
}

bool LocalMaxStatUtil::isScoreIncreasing (
size_t dimension_, // #(distinct values)
const Int4 *score_) // scores in increasing order
{
   for (size_t i = 1; i < dimension_; i++) {
      if (score_ [i] <= score_ [i - 1]) return false;
   }
   return true;
}

bool LocalMaxStatUtil::isLogarithmic (
size_t dimension_, // #(distinct values) of scores & probabilities (which are paired)
const Int4 *score_, // scores in increasing order
const double *prob_) // corresponding probabilities
{
   _ASSERT (score_);
   _ASSERT (prob_);
   if (! isScoreIncreasing (dimension_, score_)) return false;
   if (! isProbDist (dimension_, prob_)) return false;
   if (0.0 <= mu (dimension_, score_, prob_)) return false;
   if (score_ [dimension_ - 1] <= 0.0) return false;
   return true;
}