nmrPolynomialTermPowerIndex.h

Go to the documentation of this file.

/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-    */
/* ex: set filetype=cpp softtabstop=4 shiftwidth=4 tabstop=4 cindent expandtab: */

/*

  Author(s):    Ofri Sadowsky
  Created on:   2001-10-21

  (C) Copyright 2001-2007 Johns Hopkins University (JHU), All Rights
  Reserved.

--- begin cisst license - do not edit ---

This software is provided "as is" under an open source license, with
no warranty.  The complete license can be found in license.txt and
http://www.cisst.org/cisst/license.txt.

--- end cisst license ---
*/

#ifndef _nmrPolynomialTermPowerIndex_h
#define _nmrPolynomialTermPowerIndex_h

#include <vector>
#include <stdexcept>
#include <cisstCommon/cmnPortability.h>
#include <cisstNumerical/nmrExport.h>

class nmrPolynomialBase;

#ifndef DEFINE_FORMATTED_OUTPUT
#define DEFINE_FORMATTED_OUTPUT 0
#endif

class CISST_EXPORT nmrPolynomialTermPowerIndex
{
public:
    typedef int VariableIndexType;
    typedef int PowerType;
    typedef unsigned long MultinomialCoefficientType;


    nmrPolynomialTermPowerIndex(VariableIndexType numVariables = 1, PowerType minDegree = 0, PowerType maxDegree = 4);

    nmrPolynomialTermPowerIndex(const nmrPolynomialBase & p);

    int GetDegree() const
    {
        return Degree;
    }

    PowerType GetMaxDegree() const
    {
        return MaxDegree;
    }

    PowerType GetMinDegree() const
    {
        return MinDegree;
    }

    void SetMaxDegree(PowerType newMax) throw(std::runtime_error)
    {
        if (newMax < this->Degree) {
            throw std::runtime_error("nmrPolynomialTermPowerIndex: Attempt to set max degree below current degree");
        }
        if (newMax < this->MinDegree) {
            throw std::runtime_error("nmrPolynomialTermPowerIndex: Attempt to set max degree below set minimum");
        }
        MaxDegree = newMax;
    }

    void SetMinDegree(PowerType newMin) throw(std::runtime_error)
    {
        if (newMin > this->Degree) {
            throw std::runtime_error("nmrPolynomialTermPowerIndex: Attempt to set min degree above current degree");
        }
        if (newMin > this->MaxDegree) {
            throw std::runtime_error("nmrPolynomialTermPowerIndex: Attempt to set max degree above set maximum");
        }
        MinDegree = newMin;
    }

    bool IsValid() const
    {
        return (MinDegree <= Degree) && (Degree <= MaxDegree);
    }
        
    VariableIndexType GetNumVariables() const
    {
        return Powers.size();
    }
    
    PowerType GetPower(VariableIndexType variable) const
    {
        return Powers[variable];
    }

    const PowerType * GetPowers() const
    {
        return &(Powers[0]); 
    }

    const int * GetPowersAsSigned() const
    {
        return reinterpret_cast<const int *>(GetPowers());
    }

    void SetPower(VariableIndexType variable, PowerType power)
    { 
        Degree += power - Powers[variable];
        Powers[variable] = power; 
    }

    void SetPowers(const PowerType powers[])
    {
        VariableIndexType i;
        Degree = 0;
        for (i = 0; i < GetNumVariables(); i++) {
            Powers[i] = powers[i];
            Degree += powers[i];
        }
    }

    void CopyPowers(const nmrPolynomialTermPowerIndex & other);

    void SetDegree(PowerType degree);

    void GoBegin();
    void GoEnd();

    void Increment();

    void Decrement();  

    bool IsBegin() const
    {
        return (GetDegree() == GetMinDegree()) && (GetPower(0) == GetMinDegree());
    }
    
    bool IsEnd() const
    {
        return (GetDegree() == GetMaxDegree()) && (GetPower(GetNumVariables() - 1) == GetMaxDegree());
    }

    int Compare(const nmrPolynomialTermPowerIndex & other) const;

    bool operator< (const nmrPolynomialTermPowerIndex& other) const
    {
        return (this->Compare(other) < 0);
    }

    bool operator > (const nmrPolynomialTermPowerIndex& other) const
    {
        return (other < *this);
    }

    bool operator<= (const nmrPolynomialTermPowerIndex & other) const
    {
        return (this->Compare(other) <= 0);
    }

    bool operator>= (const nmrPolynomialTermPowerIndex & other) const
    {
        return (other <= *this);
    }

    bool operator== (const nmrPolynomialTermPowerIndex& other) const;
    // This is an incorrect definition of operator==
    //{ return (this->Compare(other) == 0); }


    // For a set of powers 'p1',...,'pn' whose sum is 'd', return
    // d! / (p1! p2! ... pn!).
    // In this function, d is given as an argument, so there is assumed to
    // be an implicit (non-independent) variable in the polynomial.
    static MultinomialCoefficientType 
        EvaluateMultinomialCoefficient(VariableIndexType numIndices, 
        PowerType chooseFrom, const PowerType indices[]);

    // For a set of powers 'p1',...,'pn' whose sum is 'd', return
    // d! / (p1! p2! ... pn!).
    // In this function, d is the sum of all the powers stored in 'indices'
    static MultinomialCoefficientType 
        EvaluateMultinomialCoefficient(VariableIndexType numIndices, 
        const PowerType indices[]);


    // For a set of powers 'p1',...,'pn' whose sum is 'd', return
    // d! / (p1! p2! ... pn!).
    // In this function, d is the sum of all the powers stored in this term.
    MultinomialCoefficientType GetMultinomialCoefficient() const
    {
        return EvaluateMultinomialCoefficient(GetNumVariables(), GetPowers());
    }

    // For a set of powers 'p1',...,'pn' whose sum is 'd', return
    // d! / (p1! p2! ... pn!).
    // In this function, d is given as an argument, so there is assumed to
    // be an implicit (non-independent) variable in the polynomial.
    MultinomialCoefficientType GetMultinomialCoefficient(PowerType d) const
    {
        return EvaluateMultinomialCoefficient(GetNumVariables(), d, GetPowers());
    }

    // Return the combinatorial number of possible power sets for the term.
    //
    // For a fixed degree d, and n variables, the number of combinations
    // is the number of multisets of size d over n elements.
    // First we choose k = number of variables in the term. if d<n we can choose
    // k=1..d, otherwise we can choose k=1..n . But for simplicity and arithmetic
    // reasons we will assume k=1..d. We will see that it does not affect the 
    // final result.  Next, we have to choose the variables in the term. We have
    // \choose{n}{k} combinations.  Finally we choose the power for each variable.
    // We have \choose{d-1}{k-1} combinations. The entire expression is
    //
    //   \sum_{k=1}^{d} \choose{n}{k} \choose{d-1}{k-1}
    //
    // Note that for any k>n, \choose{n}{k} = 0, so we don't need to worry about
    // d>n.
    // The entire sum can be simplified into \choose{n+d-1}{d}. However, for a term
    // index of n variables, we have to sum over all possible degrees:
    //
    //   \sum_{d=MinDegree}^{MaxDegree} \choose{n+d-1}{d}.
    //
    // \choose{n+d}{d+1} = (n+d)! / ( (d+1)! (n-1)! ) = ((n+d)(n+d-1)) / ((d+1) d! (n-1)!) 
    //                   = ((n+d)/(d+1)) * ( (n+d-1)! d! (n-1)! ) = ((n+d)/(d+1)) * \choose{n+d-1}{d}
    //
    // Therefore, we only need to make a long computation for the first element of the
    // sum, and the rest are received incrementally.
    MultinomialCoefficientType CountPowerCombinations() const;

    MultinomialCoefficientType CountLowerOrderTerms() const;

#if DEFINE_FORMATTED_OUTPUT
    // Print the powers of the term index in sprintf format into a character
    // buffer. No assertion is made to ensure that the buffer has enough capacity
    // to store the formatted message. All powers are printed, included zeroes.
    // Parameters:
    //   buffer [i/o]   -- destination buffer for the output
    //   widthFormat [i] -- a format specification to enable equal column width.
    //
    // Return: a pointer to the character immediately following the formatted output
    char * FormatPowers(char * buffer, const char * widthFormat = "%d") const;
#endif


#if DEFINE_FORMATTED_OUTPUT
    // Print the entire term index in formatted output. The format is:
    // [numVariables minDegree maxDegree] powers...
    // See FormatPowers for parameter and return value specification.
    char * FormatTermIndex(char * buffer, const char * widthFormat = "%d") const;
#endif

#if DEFINE_FORMATTED_OUTPUT
    // Return the estimated number of characters required to store the formatted
    // powers (from FormatPowers())
    int CalculateFormatPowerLength(const char * widthFormat = "%d") const
    {
        char buff[10];
        sprintf( buff, widthFormat, GetMaxDegree() );
        int varLength = strlen(buff) + 1;
        return varLength * GetNumVariables();
    }
#endif

    void SerializeRaw(std::ostream & output) const;

    void DeserializeRaw(std::istream & input);

    void SerializeIndexRaw(std::ostream & output) const;

    void DeserializeIndexRaw(std::istream & input);

protected:
    std::vector<PowerType> Powers;      // the powers of all the variables in the term
    PowerType MinDegree;    // the minimal possible degree of the term
    PowerType MaxDegree;    // the maximal possible degree of the term

    // the degree of the term. Note that it may be negative, in which case the term is
    // invalid.  This member should always be equal to the sum of all powers.  It
    // is therefore redundant information is only serves as a cache for sorting
    // purposes.  It is not serialized.
    int Degree;

};

#endif // _nmrPolynomialTermPowerIndex_h