matrix1d.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:     Carlos Oscar S. Sorzano (coss@cnb.csic.es)
 *
 * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 *
 * 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 2 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, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 * 02111-1307  USA
 *
 *  All comments concerning this program package may be sent to the
 *  e-mail address 'xmipp@cnb.csic.es'
 ***************************************************************************/

#ifndef CORE_MATRIX1D_H_
#define CORE_MATRIX1D_H_

#include <iostream>
#include <string>
#include <vector>
#include "xmipp_error.h"
#include "xmipp_macros.h"

class FileName;
template<typename T>
class Matrix2D;

#define MATRIX1D_ARRAY(v) ((v).vdata)

#define FOR_ALL_ELEMENTS_IN_MATRIX1D(v) \
    for (size_t i=0; i<v.vdim; i++)

#define VEC_XSIZE(m) ((m).vdim)

#define XX(v) (v).vdata[0]

#define YY(v) (v).vdata[1]

#define ZZ(v) (v).vdata[2]

#define VECTOR_R2(v, x, y) { \
        XX(v) = x; YY(v) = y; }

#define VECTOR_R3(v, x, y, z) { \
        XX(v) = x; YY(v) = y; ZZ(v) = z;}

#define V2_PLUS_V2(a, b, c) { \
        XX(a) = XX(b) + XX(c); \
        YY(a) = YY(b) + YY(c); }

#define V2_MINUS_V2(a, b, c) { \
        XX(a) = XX(b) - XX(c); \
        YY(a) = YY(b) - YY(c); }

#define V2_PLUS_CT(a, b, k) { \
        XX(a) = XX(b) + (k); \
        YY(a) = YY(b) + (k); }

#define V2_BY_CT(a, b, k) { \
        XX(a) = XX(b) * (k); \
        YY(a) = YY(b) * (k); }

#define V3_PLUS_V3(a, b, c) { \
        XX(a) = XX(b) + XX(c); \
        YY(a) = YY(b) + YY(c); \
        ZZ(a) = ZZ(b) + ZZ(c); }

#define V3_MINUS_V3(a, b, c) { \
        XX(a) = XX(b) - XX(c); \
        YY(a) = YY(b) - YY(c); \
        ZZ(a) = ZZ(b) - ZZ(c); }

#define V3_PLUS_CT(a, b, c) { \
        XX(a) = XX(b) + (c); \
        YY(a) = YY(b) + (c); \
        ZZ(a) = ZZ(b) + (c); }

#define V3_BY_CT(a, b, c) { \
        XX(a) = XX(b) * (c); \
        YY(a) = YY(b) * (c); \
        ZZ(a) = ZZ(b) * (c); }

#define VEC_ELEM(v,i) ((v).vdata[(i)])
#define dMi(v, i) ((v).vdata[(i)])

#define VEC_SWAP(v, i, j, aux) {aux = dMi(v, i); dMi(v, i) = dMi(v, j); dMi(v, j) = aux; }


template<typename T>
class Matrix1D
{
public:
    T* vdata;

    bool destroyData;

    size_t vdim;

    bool row;



    Matrix1D(bool column = true)
    {
        coreInit();
        row = !column;
    }

    template<class T1>
    Matrix1D(T1 dim, bool column = true)
    {
        coreInit();
        row = !column;
        initZeros(dim);
    }

    Matrix1D(const Matrix1D<T>& v)
    {
        coreInit();
        *this = v;
        row = v.row;
    }

    ~Matrix1D()
    {
        coreDeallocate();
    }

    Matrix1D<T>& operator=(const Matrix1D<T>& op1);

    bool operator==(const Matrix1D<T>& op1) const;

    Matrix1D<T>& operator=(const std::vector<T>& op1);

    void clear();

    void coreInit();

    inline void coreAllocate(int _vdim)
    {
        if (_vdim <= 0)
        {
            clear();
            return;
        }

        vdim = _vdim;
        vdata = new T[vdim];
        memset(vdata, 0, vdim * sizeof(T));
        if (vdata == NULL
           )
            REPORT_ERROR(ERR_MEM_NOTENOUGH, "Allocate: No space left");
    }

    inline void coreDeallocate()
    {
        if (vdata != NULL && destroyData)
            delete[] vdata;
        vdata = NULL;
    }



    inline void resize(size_t Xdim, bool copy = true)
    {
        if (Xdim == vdim)
            return;

        if (Xdim <= 0)
        {
            clear();
            return;
        }

        T * new_vdata;
        try
        {
            new_vdata = new T[(size_t) Xdim]; //TODO: Kino: Valgrind points here there is something wrong, but I could'nt resolve
            memset(new_vdata, 0, Xdim * sizeof(T));
        }
        catch (std::bad_alloc &)
        {
            REPORT_ERROR(ERR_MEM_NOTENOUGH, "Allocate: No space left");
        }

        // Copy needed elements, fill with 0 if necessary
        if (copy)
        {
            T zero=0; // Useful for complexes
            for (size_t j = 0; j < Xdim; j++)
            {
                T *val=NULL;
                if (j >= vdim)
                    val = &zero;
                else
                    val = &vdata[j];
                new_vdata[j] = *val;
            }
        }

        // deallocate old vector
        coreDeallocate();

        // assign *this vector to the newly created
        vdata = new_vdata;
        vdim = Xdim;

    }

    inline void resizeNoCopy(int Xdim)
    {
        resize(Xdim, false);
    }

    template<typename T1>
    void resize(const Matrix1D<T1> &v)
    {
        if (vdim != v.vdim)
            resize(v.vdim);
    }

    template<typename T1>
    void resizeNoCopy(const Matrix1D<T1> &v)
    {
        if (vdim != v.vdim)
            resize(v.vdim, false);
    }

    template<typename T1>
    bool sameShape(const Matrix1D<T1>& op) const
    {
        return (vdim == op.vdim);
    }

    inline size_t size() const
    {
        return vdim;
    }

    inline int isRow() const
    {
        return row;
    }

    inline int isCol() const
    {
        return !row;
    }

    inline void setRow()
    {
        row = true;
    }

    inline void setCol()
    {
        row = false;
    }



    void initConstant(T val);

    void initConstant(size_t Xdim, T val);


    void enumerate();

    inline void initZeros()
    {
        memset(vdata, 0, vdim * sizeof(T));
    }

    inline void initZeros(size_t Xdim)
    {
        if (vdim != Xdim)
            resizeNoCopy(Xdim);
        memset(vdata, 0, vdim * sizeof(T));
    }

    template<typename T1>
    inline void initZeros(const Matrix1D<T1>& op)
    {
        if (vdim != op.vdim)
            resize(op);
        memset(vdata, 0, vdim * sizeof(T));
    }

    bool isAnyNaN();



    Matrix1D<T> operator*(T op1) const;

    Matrix1D<T> operator/(T op1) const;

    Matrix1D<T> operator+(T op1) const;

    Matrix1D<T> operator-(T op1) const;

    friend Matrix1D<T> operator*(T op1, const Matrix1D<T>& op2)
    {
        Matrix1D<T> tmp(op2);
        T *ptr1 = &VEC_ELEM(tmp,0);
        const T *ptr2 = &VEC_ELEM(op2,0);
        size_t iBlockMax = op2.vdim / 4;
        for (size_t i = 0; i < iBlockMax; i++)
        {
            (*ptr1++) = (*ptr2++) * op1;
            (*ptr1++) = (*ptr2++) * op1;
            (*ptr1++) = (*ptr2++) * op1;
            (*ptr1++) = (*ptr2++) * op1;
        }
        for (size_t i = iBlockMax * 4; i < op2.vdim; ++i)
            (*ptr1++) = (*ptr2++) * op1;
        return tmp;
    }

    friend Matrix1D<T> operator/(T op1, const Matrix1D<T>& op2)
    {
        Matrix1D<T> tmp(op2);
        T *ptr1 = &VEC_ELEM(tmp,0);
        const T *ptr2 = &VEC_ELEM(op2,0);
        size_t iBlockMax = op2.vdim / 4;
        for (size_t i = 0; i < iBlockMax; i++)
        {
            (*ptr1++) = op1 / (*ptr2++);
            (*ptr1++) = op1 / (*ptr2++);
            (*ptr1++) = op1 / (*ptr2++);
            (*ptr1++) = op1 / (*ptr2++);
        }
        for (size_t i = iBlockMax * 4; i < op2.vdim; ++i)
            (*ptr1++) = op1 / (*ptr2++);
        return tmp;
    }

    friend Matrix1D<T> operator+(T op1, const Matrix1D<T>& op2)
    {
        Matrix1D<T> tmp(op2);
        T *ptr1 = &VEC_ELEM(tmp,0);
        const T *ptr2 = &VEC_ELEM(op2,0);
        size_t iBlockMax = op2.vdim / 4;
        for (size_t i = 0; i < iBlockMax; i++)
        {
            (*ptr1++) = (*ptr2++) + op1;
            (*ptr1++) = (*ptr2++) + op1;
            (*ptr1++) = (*ptr2++) + op1;
            (*ptr1++) = (*ptr2++) + op1;
        }
        for (size_t i = iBlockMax * 4; i < op2.vdim; ++i)
            (*ptr1++) = (*ptr2++) + op1;
        return tmp;
    }

    friend Matrix1D<T> operator-(T op1, const Matrix1D<T>& op2)
    {
        Matrix1D<T> tmp(op2);
        T *ptr1 = &VEC_ELEM(tmp,0);
        const T *ptr2 = &VEC_ELEM(op2,0);
        size_t iBlockMax = op2.vdim / 4;
        for (size_t i = 0; i < iBlockMax; i++)
        {
            (*ptr1++) = op1 - (*ptr2++);
            (*ptr1++) = op1 - (*ptr2++);
            (*ptr1++) = op1 - (*ptr2++);
            (*ptr1++) = op1 - (*ptr2++);
        }
        for (size_t i = iBlockMax * 4; i < op2.vdim; ++i)
            (*ptr1++) = op1 - (*ptr2++);
        return tmp;
    }

    void operator+=(const Matrix1D<T>& op1) const;

    void operator-=(const Matrix1D<T>& op1) const;

    void operator*=(T op1);

    void operator/=(T op1);

    void operator+=(T op1);

    void operator-=(T op1);

    Matrix1D<T> operator*(const Matrix1D<T>& op1) const;

    Matrix1D<T> operator/(const Matrix1D<T>& op1) const;
    Matrix1D<T> operator+(const Matrix1D<T>& op1) const;

    Matrix1D<T> operator-(const Matrix1D<T>& op1) const;

    void operator*=(const Matrix1D<T>& op1);

    void operator/=(const Matrix1D<T>& op1);

    void operator+=(const Matrix1D<T>& op1);

    void operator-=(const Matrix1D<T>& op1);

    Matrix1D<T> operator-() const;

    Matrix1D<T> operator*(const Matrix2D<T>& M);

    inline T& operator()(int i) const
    {
        return vdata[i];
    }

    inline T& operator [](int idx) {
        return vdata[idx];
    }

    inline T operator [](int idx) const {
        return vdata[idx];
    }



    inline T* adaptForNumericalRecipes() const
    {
        return MATRIX1D_ARRAY(*this) - 1;
    }

    void killAdaptationForNumericalRecipes(T* m) const
        {}

    void selfCEIL();

    void selfFLOOR();

    void selfROUND();

    Matrix1D<T> sort() const;

    void indexSort(Matrix1D<int> &indx) const;

    double computeMean() const;

    T computeMax() const;

    void computeMinMax(T &minval, T &maxval) const;

    void computeMeanAndStddev(double &mean, double &stddev) const;

    T computeMedian() const;

    int maxIndex() const;

    int minIndex() const;

    Matrix1D<T> transpose() const;

    inline void selfTranspose()
    {
        row = !row;
    }

    double sum(bool average = false) const;

    double sum2() const;

    double dotProduct(const Matrix1D<T> &op1) const;

    inline double module() const
    {
        return sqrt(sum2());
    }

    inline double angle() const
    {
        return atan2((double) YY(*this), (double) XX(*this));
    }

    void selfNormalize();

    void selfReverse();

    void showWithGnuPlot(const std::string& xlabel, const std::string& title);

    void numericalDerivative(Matrix1D<double> &result) const;

    void read(const FileName& fn);

    // This function must be explicitly implemented outside
    friend std::istream& operator>>(std::istream& in, Matrix1D<T>& v)
    {
        for (int i = 0; i < VEC_XSIZE(v); i++)
            in >> VEC_ELEM(v,i);
        return in;
    }

    friend std::ostream& operator<<(std::ostream& ostrm, const Matrix1D<T>& v)
    {
        if (v.vdim == 0)
            ostrm << "NULL Array\n";
        else
            ostrm << std::endl;

        double max_val = ABS(v.vdata[0]);

        for (size_t j = 0; j < v.vdim; j++)
        {
            max_val = XMIPP_MAX(max_val, v.vdata[j]);
        }

        int prec = bestPrecision(max_val, 10);

        if (v.row)
            for (size_t j = 0; j < v.vdim; j++)
            {
                ostrm << floatToString((double) v.vdata[j], 10, prec)
                << std::endl;
            }
        else
            for (size_t j = 0; j < v.vdim; j++)
            {
                ostrm << floatToString((double) v.vdata[j], 10, prec) << " ";
            }

        return ostrm;
    }

    void write(const FileName& fn) const;

    void edit();
};


 typedef Matrix1D<double> DVector;
 typedef Matrix1D<int> IVector;

Matrix1D<double> vectorR2(double x, double y);

Matrix1D<double> vectorR3(double x, double y, double z);

Matrix1D<int> vectorR3(int x, int y, int z);

template<typename T>
T dotProduct(const Matrix1D<T>& v1, const Matrix1D<T>& v2)
{
    if (!v1.sameShape(v2))
        REPORT_ERROR(ERR_MATRIX_SIZE,
                     "Dot product: vectors of different size or shape");

    T accumulate = 0;
    for (size_t j = 0; j < v1.vdim; j++)
        accumulate += v1.vdata[j] * v2.vdata[j];
    return accumulate;
}

template<typename T>
Matrix1D<T> vectorProduct(const Matrix1D<T>& v1, const Matrix1D<T>& v2)
{
    if (v1.vdim != 3 || v2.vdim != 3)
        REPORT_ERROR(ERR_MATRIX_SIZE, "Vector_product: vectors are not in R3");

    if (v1.isRow() != v2.isRow())
        REPORT_ERROR(ERR_MATRIX_SIZE,
                     "Vector_product: vectors are of different shape");

    Matrix1D<T> result(3);
    vectorProduct(v1, v2, result);
    return result;
}

template<typename T>
void vectorProduct(const Matrix1D<T>& v1, const Matrix1D<T>& v2,
                   Matrix1D<T> &result)
{
    XX(result) = YY(v1) * ZZ(v2) - ZZ(v1) * YY(v2);
    YY(result) = ZZ(v1) * XX(v2) - XX(v1) * ZZ(v2);
    ZZ(result) = XX(v1) * YY(v2) - YY(v1) * XX(v2);
}

template<typename T>
void sortTwoVectors(Matrix1D<T>& v1, Matrix1D<T>& v2)
{
    T temp;
    if (!v1.sameShape(v2))
        REPORT_ERROR(ERR_MATRIX_SIZE,
                     "sortTwoVectors: vectors are not of the same shape");

    for (size_t j = 0; j < v1.vdim; j++)
    {
        temp = XMIPP_MIN(v1.vdata[j], v2.vdata[j]);
        v2.vdata[j] = XMIPP_MAX(v1.vdata[j], v2.vdata[j]);
        v1.vdata[j] = temp;
    }
}

template<typename T1, typename T2>
void typeCast(const Matrix1D<T1>& v1, Matrix1D<T2>& v2)
{
    if (v1.vdim == 0)
    {
        v2.clear();
        return;
    }

    v2.resizeNoCopy(v1.vdim);
    for (size_t j = 0; j < v1.vdim; j++)
        v2.vdata[j] = static_cast<T2>(v1.vdata[j]);
}

template<typename T1>
void typeCast(const Matrix1D<T1>& v1, Matrix1D<T1>& v2)
{
    v2 = v1;
}


#endif /* MATRIX1D_H_ */