UMath.h

Go to the documentation of this file.

/*
*  utilite is a cross-platform library with
*  useful utilities for fast and small developing.
*  Copyright (C) 2010  Mathieu Labbe
*
*  utilite is free library: you can redistribute it and/or modify
*  it under the terms of the GNU Lesser General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  utilite 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 Lesser General Public License for more details.
*
*  You should have received a copy of the GNU Lesser General Public License
*  along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#ifndef UMATH_H
#define UMATH_H
 
#include <cmath>
#include <list>
#include <vector>
 
#if _MSC_VER
        #undef min
        #undef max
#endif
 
template<class T>
inline bool uIsNan(const T & value)
{
#if _MSC_VER
        return _isnan(value) != 0;
#else
        return std::isnan(value);
#endif
}
 
template<class T>
inline bool uIsFinite(const T & value)
{
#if _MSC_VER
        return _finite(value) != 0;
#else
        return std::isfinite(value);
#endif
}
 
template<class T>
inline T uMin3( const T& a, const T& b, const T& c)
{
        float m=a<b?a:b;
        return m<c?m:c;
}
 
template<class T>
inline T uMax3( const T& a, const T& b, const T& c)
{
        float m=a>b?a:b;
        return m>c?m:c;
}
 
template<class T>
inline T uMax(const T * v, unsigned int size, unsigned int & index)
{
        T max = 0;
        index = 0;
        if(!v || size == 0)
        {
                return max;
        }
        max = v[0];
        for(unsigned int i=1; i<size; ++i)
        {
                if(uIsNan(max) || (max < v[i] && !uIsNan(v[i])))
                {
                        max = v[i];
                        index = i;
                }
        }
 
        return max;
}
 
template<class T>
inline T uMax(const std::vector<T> & v, unsigned int & index)
{
        return uMax(v.data(), v->size(), index);
}
 
template<class T>
inline T uMax(const T * v, unsigned int size)
{
        unsigned int index;
        return uMax(v, size, index);
}
 
template<class T>
inline T uMax(const std::vector<T> & v)
{
        return uMax(v.data(), v.size());
}
 
template<class T>
inline T uMin(const T * v, unsigned int size, unsigned int & index)
{
        T min = 0;
        index = 0;
        if(!v || size == 0)
        {
                return min;
        }
        min = v[0];
        for(unsigned int i=1; i<size; ++i)
        {
                if(uIsNan(min) || (min > v[i] && !uIsNan(v[i])))
                {
                        min = v[i];
                        index = i;
                }
        }
 
        return min;
}
 
template<class T>
inline T uMin(const std::vector<T> & v, unsigned int & index)
{
        return uMin(v.data(), v.size(), index);
}
 
template<class T>
inline T uMin(const T * v, unsigned int size)
{
        unsigned int index;
        return uMin(v, size, index);
}
 
template<class T>
inline T uMin(const std::vector<T> & v)
{
        return uMin(v.data(), v.size());
}
 
template<class T>
inline void uMinMax(const T * v, unsigned int size, T & min, T & max, unsigned int & indexMin, unsigned int & indexMax)
{
        min = 0;
        max = 0;
        indexMin = 0;
        indexMax = 0;
        if(!v || size == 0)
        {
                return;
        }
        min = v[0];
        max = v[0];
 
        for(unsigned int i=1; i<size; ++i)
        {
                if(uIsNan(min) || (min > v[i] && !uIsNan(v[i])))
                {
                        min = v[i];
                        indexMin = i;
                }
 
                if(uIsNan(max) || (max < v[i] && !uIsNan(v[i])))
                {
                        max = v[i];
                        indexMax = i;
                }
        }
}
 
template<class T>
inline void uMinMax(const std::vector<T> & v, T & min, T & max, unsigned int & indexMin, unsigned int & indexMax)
{
        uMinMax(v.data(), v.size(), min, max, indexMin, indexMax);
}
 
template<class T>
inline void uMinMax(const T * v, unsigned int size, T & min, T & max)
{
        unsigned int indexMin;
        unsigned int indexMax;
        uMinMax(v, size, min, max, indexMin, indexMax);
}
 
template<class T>
inline void uMinMax(const std::vector<T> & v, T & min, T & max)
{
        uMinMax(v.data(), v.size(), min, max);
}
 
template<class T>
inline int uSign(const T & v)
{
        if(v < 0)
        {
                return -1;
        }
        else
        {
                return 1;
        }
}
 
template<class T>
inline T uSum(const std::list<T> & list)
{
        T sum = 0;
        for(typename std::list<T>::const_iterator i=list.begin(); i!=list.end(); ++i)
        {
                sum += *i;
        }
        return sum;
}
 
template<class T>
inline T uSum(const T * v, unsigned int size)
{
        T sum = 0;
        if(v && size)
        {
                for(unsigned int i=0; i<size; ++i)
                {
                        sum += v[i];
                }
        }
        return sum;
}
 
template<class T>
inline T uSum(const std::vector<T> & v)
{
        return uSum(v.data(), (int)v.size());
}
 
template<class T>
inline T uSumSquared(const T * v, unsigned int size, T subtract = T())
{
        T sum = 0;
        if(v && size)
        {
                for(unsigned int i=0; i<size; ++i)
                {
                        sum += (v[i]-subtract)*(v[i]-subtract);
                }
        }
        return sum;
}
 
template<class T>
inline T uSumSquared(const std::vector<T> & v, T subtract = T())
{
        return uSumSquared(v.data(), v.size(), subtract);
}
 
template<class T>
inline T uMean(const T * v, unsigned int size)
{
        T buf = 0;
        if(v && size)
        {
                for(unsigned int i=0; i<size; ++i)
                {
                        buf += v[i];
                }
                buf /= size;
        }
        return buf;
}
 
template<class T>
inline T uMean(const std::list<T> & list)
{
        T m = 0;
        if(list.size())
        {
                for(typename std::list<T>::const_iterator i=list.begin(); i!=list.end(); ++i)
                {
                        m += *i;
                }
                m /= list.size();
        }
        return m;
}
 
template<class T>
inline T uMean(const std::vector<T> & v)
{
        return uMean(v.data(), v.size());
}
 
template<class T>
inline T uMeanSquaredError(const T * x, unsigned int sizeX, const T * y, unsigned int sizeY)
{
        T sum = 0;
        if(x && y && sizeX == sizeY)
        {
                for(unsigned int i=0; i<sizeX; ++i)
                {
                        T diff = x[i]-y[i];
                        sum += diff*diff;
                }
                return sum/(T)sizeX;
        }
        return (T)-1;
}
 
template<class T>
inline T uMeanSquaredError(const std::vector<T> & x, const std::vector<T> & y)
{
        return uMeanSquaredError(x.data(), x.size(), y.data(), y.size());
}
 
template<class T>
inline T uVariance(const T * v, unsigned int size, T meanV)
{
        T buf = 0;
        if(v && size>1)
        {
                float sum = 0;
                for(unsigned int i=0; i<size; ++i)
                {
                        sum += (v[i]-meanV)*(v[i]-meanV);
                }
                buf = sum/(size-1);
        }
        return buf;
}
 
template<class T>
inline T uVariance(const std::list<T> & list, const T & m)
{
        T buf = 0;
        if(list.size()>1)
        {
                float sum = 0;
                for(typename std::list<T>::const_iterator i=list.begin(); i!=list.end(); ++i)
                {
                        sum += (*i-m)*(*i-m);
                }
                buf = sum/(list.size()-1);
        }
        return buf;
}
 
template<class T>
inline T uVariance(const T * v, unsigned int size)
{
        T m = uMean(v, size);
        return uVariance(v, size, m);
}
 
template<class T>
inline T uVariance(const std::vector<T> & v, const T & m)
{
        return uVariance(v.data(), v.size(), m);
}
 
template<class T>
inline T uNormSquared(const std::vector<T> & v)
{
        float sum = 0.0f;
        for(unsigned int i=0; i<v.size(); ++i)
        {
                sum += v[i]*v[i];
        }
        return sum;
}
 
template<class T>
inline T uNorm(const std::vector<T> & v)
{
        return std::sqrt(uNormSquared(v));
}
 
template<class T>
inline T uNormSquared(const T & x1, const T & x2)
{
        return x1*x1 + x2*x2;
}
 
template<class T>
inline T uNorm(const T & x1, const T & x2)
{
        return std::sqrt(uNormSquared(x1, x2));
}
 
template<class T>
inline T uNormSquared(const T & x1, const T & x2, const T & x3)
{
        return x1*x1 + x2*x2 + x3*x3;
}
 
template<class T>
inline T uNorm(const T & x1, const T & x2, const T & x3)
{
        return std::sqrt(uNormSquared(x1, x2, x3));
}
 
template<class T>
inline std::vector<T> uNormalize(const std::vector<T> & v)
{
        float norm = uNorm(v);
        if(norm == 0)
        {
                return v;
        }
        else
        {
                std::vector<T> r(v.size());
                for(unsigned int i=0; i<v.size(); ++i)
                {
                        r[i] = v[i]/norm;
                }
                return r;
        }
}
 
template<class T>
inline std::list<unsigned int> uLocalMaxima(const T * v, unsigned int size)
{
        std::list<unsigned int> maxima;
        if(size)
        {
                for(unsigned int i=0; i<size; ++i)
                {
                        if(i == 0)
                        {
                                // first item
                                if((i+1 < size && v[i] > v[i+1]) ||
                                        i+1 >= size)
                                {
                                        maxima.push_back(i);
                                }
                        }
                        else if(i == size - 1)
                        {
                                //last item
                                if((i >= 1 && v[i] > v[i-1]) ||
                                        i == 0)
                                {
                                        maxima.push_back(i);
                                }
                        }
                        else
                        {
                                //all others, check previous and next
                                if(v[i] > v[i-1] && v[i] > v[i+1])
                                {
                                        maxima.push_back(i);
                                }
                        }
                }
        }
        return maxima;
}
 
template<class T>
inline std::list<unsigned int> uLocalMaxima(const std::vector<T> & v)
{
        return uLocalMaxima(v.data(), v.size());
}
 
enum UXMatchMethod{UXCorrRaw, UXCorrBiased, UXCorrUnbiased, UXCorrCoeff, UXCovRaw, UXCovBiased, UXCovUnbiased, UXCovCoeff};
 
template<class T>
inline std::vector<T> uXMatch(const T * vA, const T * vB, unsigned int sizeA, unsigned int sizeB, UXMatchMethod method)
{
        if(!vA || !vB || sizeA == 0 || sizeB == 0)
        {
                return std::vector<T>();
        }
 
        std::vector<T> result(sizeA + sizeB - 1);
 
        T meanA = 0;
        T meanB = 0;
        if(method > UXCorrCoeff)
        {
                meanA = uMean(vA, sizeA);
                meanB = uMean(vB, sizeB);
        }
 
        T den = 1;
        if(method == UXCorrCoeff || method == UXCovCoeff)
        {
                den = std::sqrt(uSumSquared(vA, sizeA, meanA) * uSumSquared(vB, sizeB, meanB));
        }
        else if(method == UXCorrBiased || method == UXCovBiased)
        {
                den = (T)std::max(sizeA, sizeB);
        }
 
        if(sizeA == sizeB)
        {
                T resultA;
                T resultB;
 
                int posA;
                int posB;
                unsigned int j;
 
                // Optimization, filling two results at once
                for(unsigned int i=0; i<sizeA; ++i)
                {
                        if(method == UXCorrUnbiased || method == UXCovUnbiased)
                        {
                                den = 0;
                        }
 
                        posA = sizeA - i - 1;
                        posB = sizeB - i - 1;
                        resultA = 0;
                        resultB = 0;
                        for(j=0; (j + posB) < sizeB && (j + posA) < sizeA; ++j)
                        {
                                resultA += (vA[j] - meanA) * (vB[j + posB] - meanB);
                                resultB += (vA[j + posA] - meanA) * (vB[j] - meanB);
                                if(method == UXCorrUnbiased || method == UXCovUnbiased)
                                {
                                        ++den;
                                }
                        }
 
                        result[i] = resultA / den;
                        result[result.size()-1 -i] = resultB / den;
                }
        }
        else
        {
                for(unsigned int i=0; i<result.size(); ++i)
                {
                        if(method == UXCorrUnbiased || method == UXCovUnbiased)
                        {
                                den = 0;
                        }
 
                        int posB = sizeB - i - 1;
                        T r = 0;
                        if(posB >= 0)
                        {
                                for(unsigned int j=0; (j + posB) < sizeB && j < sizeA; ++j)
                                {
                                        r += (vA[j] - meanA) * (vB[j + posB] - meanB);
                                        if(method == UXCorrUnbiased || method == UXCovUnbiased)
                                        {
                                                ++den;
                                        }
                                }
                        }
                        else
                        {
                                int posA = posB*-1;
                                for(unsigned int i=0; (i+posA) < sizeA && i < sizeB; ++i)
                                {
                                        r += (vA[i+posA] - meanA) * (vB[i] - meanB);
                                        if(method == UXCorrUnbiased || method == UXCovUnbiased)
                                        {
                                                ++den;
                                        }
                                }
                        }
 
                        result[i] = r / den;
                }
        }
 
        return result;
}
 
template<class T>
inline std::vector<T> uXMatch(const std::vector<T> & vA, const std::vector<T> & vB, UXMatchMethod method)
{
        return uXMatch(vA.data(), vB.data(), vA.size(), vB.size(), method);
}
 
template<class T>
inline T uXMatch(const T * vA, const T * vB, unsigned int sizeA, unsigned int sizeB, unsigned int index, UXMatchMethod method)
{
        T result = 0;
        if(!vA || !vB || sizeA == 0 || sizeB == 0)
        {
                return result;
        }
 
        T meanA = 0;
        T meanB = 0;
        if(method > UXCorrCoeff)
        {
                meanA = uMean(vA, sizeA);
                meanB = uMean(vB, sizeB);
        }
        unsigned int size = sizeA + sizeB - 1;
 
        T den = 1;
        if(method == UXCorrCoeff || method == UXCovCoeff)
        {
                den = std::sqrt(uSumSquared(vA, sizeA, meanA) * uSumSquared(vB, sizeB, meanB));
        }
        else if(method == UXCorrBiased || method == UXCovBiased)
        {
                den = (T)std::max(sizeA, sizeB);
        }
        else if(method == UXCorrUnbiased || method == UXCovUnbiased)
        {
                den = 0;
        }
 
        if(index < size)
        {
                int posB = sizeB - index - 1;
                unsigned int i;
                if(posB >= 0)
                {
                        for(i=0; (i + posB) < sizeB && i < sizeA; ++i)
                        {
                                result += (vA[i] - meanA) * (vB[i + posB] - meanB);
                                if(method == UXCorrUnbiased || method == UXCovUnbiased)
                                {
                                        ++den;
                                }
                        }
                }
                else
                {
                        int posA = posB*-1;
                        for(i=0; (i+posA) < sizeA && i < sizeB; ++i)
                        {
                                result += (vA[i+posA] - meanA) * (vB[i] - meanB);
                                if(method == UXCorrUnbiased || method == UXCovUnbiased)
                                {
                                        ++den;
                                }
                        }
                }
        }
        return result / den;
}
 
template<class T>
inline T uXMatch(const std::vector<T> & vA, const std::vector<T> & vB, unsigned int index, UXMatchMethod method)
{
        return uXMatch(vA.data(), vB.data(), vA.size(), vB.size(), index, method);
}
 
inline std::vector<float> uHamming(unsigned int L)
{
        std::vector<float> w(L);
        unsigned int N = L-1;
        float pi = 3.14159265f;
        for(unsigned int n=0; n<N; ++n)
        {
                w[n] = 0.54f-0.46f*std::cos(2.0f*pi*float(n)/float(N));
        }
        return w;
}
 
template <typename T>
bool uIsInBounds(const T& value, const T& low, const T& high)
{
        return uIsFinite(value) && !(value < low) && !(value >= high);
}
 
#endif // UMATH_H