SNmultiGaussian.h

Go to the documentation of this file.

/*
Copyright 2017 Laurent Claessens
contact : laurent@claessens-donadello.eu

This program 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 3 of the License, or
(at your option) any later version.

This program 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, see <http://www.gnu.org/licenses/>.
*/

#ifndef __SNMULTIGAUSSIAN_H__105525
#define __SNMULTIGAUSSIAN_H__105525


#include "SNlowerTriangular.h"
#include "SNidentity.h"
#include "m_num.h"

template <class T,unsigned int tp_size>
class SNmultiGaussian : public SNgeneric<T,tp_size>
{
    
    private :
        SNlowerTriangular<T,tp_size> data_L;
        m_num data_last_column;       // the last non trivial column

        SpecialValue<T> checkForSpecialElements(const m_num&,const m_num&) const;
        T _get(m_num,m_num) const override;
        T& _at(m_num,m_num) override;
    public:
        SNmultiGaussian();
        SNmultiGaussian(const SNgeneric<T,tp_size>& A);
        SNmultiGaussian(const SNgaussian<T,tp_size>& A);
        SNmultiGaussian(const SNmultiGaussian<T,tp_size>& A);

        void swap(SNmultiGaussian<T,tp_size>&);
        SNmultiGaussian<T,tp_size>& operator=(SNmultiGaussian<T,tp_size>);

        void operator *=(const SNgaussian<T,tp_size>& other);

        SNmultiGaussian<T,tp_size> inverse() const;

        /*\brief Swap the lines \f$ i \f$ and \f$ j \f$.
         *
         * This only works when \f$ i \f$ and \f$ j \f$ are larger than
         * the last non trivial column.
         *
         * If \f$ i \f$ or \f$ j \f$ is lower or equal to the last non trivial line,
         * throw a `ProbablyNotWhatYouWantException`.
         *
         * In place replacement.
         */
        void swapLines(const m_num& i,const m_num& j);

        m_num getLastColumn() const;
        void setLastColumn(const m_num& lc);

        template <class U,unsigned int s>
        void setFirstLines(const SNgeneric<U,s>& other,const m_num max_l);
};

// CONSTRUCTORS -------------------------------------------------

// from nothing
template <class T,unsigned int tp_size>
SNmultiGaussian<T,tp_size>::SNmultiGaussian():
    data_L(SNidentity<T,tp_size>()),
    data_last_column(tp_size+1)     //force the user to initialize (see `_at`)
{ }

// from generic
template <class T,unsigned int tp_size>
SNmultiGaussian<T,tp_size>::SNmultiGaussian(const SNgeneric<T,tp_size>& A):
    data_L(A.getGaussian(0)),
    data_last_column(0)
{  }

// from multigaussian
template <class T,unsigned int tp_size>
SNmultiGaussian<T,tp_size>::SNmultiGaussian(const SNmultiGaussian<T,tp_size>& A):
    data_L(A.data_L),
    data_last_column(A.getLastColumn())
{ }

// from gaussian
template <class T,unsigned int tp_size>
SNmultiGaussian<T,tp_size>::SNmultiGaussian(const SNgaussian<T,tp_size>& A):
    data_last_column(A.getColumn())
{ 
    for (m_num c=0;c<A.getColumn();++c)
    {
        for (m_num l=c+1;l<tp_size;++l)
        {
            this->at(l,c)=A.get(l,c);
        }
    }
    for (m_num l=A.getColumn()+1;l<tp_size;++l)
    {
        this->at(l,A.getColumn())=A.get(l,A.getColumn());
    }
}

template <class T,unsigned int tp_size>
void SNmultiGaussian<T,tp_size>::swap(SNmultiGaussian<T,tp_size>& other)
{
    data_L.swap(other.data_L);
    data_last_column.swap(other.data_last_column);
}

// assignation
template <class T,unsigned int tp_size>
SNmultiGaussian<T,tp_size>& SNmultiGaussian<T,tp_size>::operator=(SNmultiGaussian<T,tp_size> other)
{
    swap(other);
    return *this;
}

// GETTER/SETTER METHODS  ---------------------------------------


template <class T,unsigned int tp_size>
m_num SNmultiGaussian<T,tp_size>::getLastColumn() const
{
    return data_last_column;
}

template <class T,unsigned int tp_size>
void SNmultiGaussian<T,tp_size>::setLastColumn(const m_num& lc)
{
    if (lc>tp_size-1)   // makes no sense to have a gaussian 
                        // behaviour on the last line.
    {
        throw OutOfRangeColumnNumber("The specified column number is larger than the size of the matrix.");
    }
    data_last_column=lc;
}

template <class T,unsigned int tp_size>
template <class U,unsigned int s>
void SNmultiGaussian<T,tp_size>::setFirstLines
        (const SNgeneric<U,s>& other,const m_num max_l)
{
    for (m_num line=0;line<max_l+1;line++)
    {
        for (m_num col=0; col < line ;++col)
        {
            this->at(line,col)=other.get(line,col);
        }
    }
}

// OPERATORS  ---------------------------------------

template <class T,unsigned int tp_size>
void SNmultiGaussian<T,tp_size>::operator *=(const SNgaussian<T,tp_size>& other)
{
    checkSizeCompatibility(*this,other);
    if (other.getColumn()!=data_last_column+1)
    {
        throw ProbablyNotWhatYouWantException("You are trying to multiply a multi-Gaussian matrix by a gaussian matrix whose column is not the next one. This is mathematically possible, but probably not what you want. However; this situation is not yet implemented.");
    }
    data_last_column++;
    for (m_num l=other.getColumn()+1;l<tp_size;++l)
    {
        this->at(l,other.getColumn())+=other.get(l,other.getColumn());
    }
}

// MATHEMATICS -------------------------------------------


template <class T,unsigned int tp_size>
SNmultiGaussian<T,tp_size> SNmultiGaussian<T,tp_size>::inverse() const
{
    auto ans=*this;
    for (m_num col=0;col<getLastColumn()+1;++col)
    {
        // one can parallelize with respect to the columns,
        // not with respect to the lines.
        for (m_num line=col+1;line<tp_size;++line)
        {
            T acc=0;
            for (m_num k=col;k<line;++k)
            {
                acc+=this->get(line,k)*ans.get(k,col);
            }
            ans.at(line,col)= - acc;
        }
    }
    return ans;
}


template <class T,unsigned int tp_size>
void SNmultiGaussian<T,tp_size>::swapLines(const m_num& i,const m_num& j)
{
    if (i<=getLastColumn() or j<=getLastColumn())
    {
        throw ProbablyNotWhatYouWantException("You are trying to swap the lines "+std::to_string(i)+" and "+std::to_string(j)+" while the last non trivial column is "+std::to_string(getLastColumn()) );
    }
    for (m_num col=0;col<=getLastColumn();++col)
    {
        std::swap( this->at(i,col),this->at(j,col)  );
    }
}

// UTILITIES  ---------------------------------------


template <class T,unsigned int tp_size>
SpecialValue<T> SNmultiGaussian<T,tp_size>::checkForSpecialElements(const m_num& i,const m_num& j) const
{
    if (i==j)
    {
        return SpecialValue<T>(1,true);
    }
    if (j>data_last_column)
    {
        return SpecialValue<T>(0,true);
    }
    if (i<j)
    {
        return SpecialValue<T>(0,true);
    }
    return SpecialValue<T>(0,false);      // In this case, the value "0" is dummy.
}


// _GET AND _AT METHODS ---------------------------------------

template <class T,unsigned int tp_size>
T SNmultiGaussian<T,tp_size>::_get(m_num i,m_num j) const
{
    SpecialValue<T> sv=checkForSpecialElements(i,j);
    if (sv.special)
    {
        return sv.value;
    }
    return data_L.get(i,j);  //if you change here, you have to change _at
}

template <class T,unsigned int tp_size>
T& SNmultiGaussian<T,tp_size>::_at(m_num i,m_num j) 
{
    if (data_last_column==tp_size+1)
    {
        throw NotInitializedMemberException("You are trying to populate a 'SNmultiGaussian' before to initialize the member 'data_last_column'. Use setLastColumn().");
    }
    SpecialValue<T> sv=checkForSpecialElements(i,j);
    if (sv.special)
    {
        throw SNchangeNotAllowedException(i,j);
    }
    return data_L.at(i,j);  //if you change here, you have to change _get
}



#endif