src/mathb.h

/*****************************************************************************
* This file is provided under the Creative Commons Attribution 3.0 license.
*
* You are free to share, copy, distribute, transmit, or adapt this work
* PROVIDED THAT you attribute the work to the authors listed below.
* For more information, please see the following web page:
* http://creativecommons.org/licenses/by/3.0/
*
* This file is a component of the Sleipnir library for functional genomics,
* authored by:
* Curtis Huttenhower (chuttenh@princeton.edu)
* Mark Schroeder
* Maria D. Chikina
* Olga G. Troyanskaya (ogt@princeton.edu, primary contact)
*
* If you use this library, the included executable tools, or any related
* code in your work, please cite the following publication:
* Curtis Huttenhower, Mark Schroeder, Maria D. Chikina, and
* Olga G. Troyanskaya.
* "The Sleipnir library for computational functional genomics"
*****************************************************************************/
#ifndef MATH_H
#define MATH_H

#include <math.h>

#include "mathbi.h"

namespace Sleipnir {

class CMath : CMathImpl {
public:
    static double LogFact( size_t iN );

    static double Sigmoid( double dHeight, double dShift, double dSlope, double dVertical, double dX ) {

        return ( dVertical + ( dHeight / ( 1 + exp( -dSlope * ( dX - dShift ) ) ) ) ); }

    static size_t GCD( size_t iA, size_t iB ) {
        size_t  i;

        while( iB ) {
            i = iA;
            iA = iB;
            iB = i % iB; }

        return iA; }

    static size_t Round( double d ) {

        return (size_t)floor( d + 0.5 ); }

    template<class tType, class tIter>
    static bool LeastSquares( const tIter BeginY, const tIter EndY, const tIter BeginX, const tIter EndX,
        tType& Alpha, tType& Beta, bool fAlpha = true ) {
        tIter   CurX, CurY;
        tType   AveX, AveY, AveXY, AveX2, Tmp;
        size_t  iN;

        AveX = AveY = AveXY = AveX2 = 0;
        for( iN = 0,CurX = BeginX,CurY = BeginY; ( CurX != EndX ) && ( CurY != EndY ); ++CurX,++CurY,++iN ) {
            AveX += *CurX;
            AveY += *CurY;
            AveXY += *CurX * *CurY;
            AveX2 += *CurX * *CurX; }
        if( !iN )
            return false;
        if( !fAlpha ) {
            Beta = AveXY / AveX2;
            return ( AveX2 != 0 ); }
        AveX /= iN;
        AveX2 /= iN;
        Tmp = AveX * AveX;
        if( Tmp == AveX2 )
            return false;
        AveY /= iN;
        AveXY /= iN;

        Beta = ( AveXY - ( AveX * AveY ) ) / ( AveX2 - Tmp );
        Alpha = AveY - ( Beta * AveX );
        return true; }
};

}

#endif // MATH_H