/****************************************************************************
*
*                                 M U E S L I   v 1.8
*
*
*     Copyright 2020 IMDEA Materials Institute, Getafe, Madrid, Spain
*     Contact: muesli.materials@imdea.org
*     Author: Ignacio Romero (ignacio.romero@imdea.org)
*
*     This file is part of MUESLI.
*
*     MUESLI 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.
*
*     MUESLI 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 MUESLI.  If not, see <http://www.gnu.org/licenses/>.
*
****************************************************************************/

#ifndef WITHEIGEN
#include "muesli/Math/mrealvector.h"

#include <iostream>
#include <iomanip>
#include <algorithm>
#include "cassert"

#include <cstdarg>
#include <valarray>
#include <vector>


using namespace muesli;


static double randomUniform(const double a, const double b);

double randomUniform(const double low, const double up)
{
    static bool initialized = false;

    if (!initialized)
    {
        srand((unsigned)time(0));
        initialized = true;
    }

    int     random_int = rand();
    double  random_dou = ( static_cast<double_t>(random_int))/ static_cast<double_t>(RAND_MAX);

    return random_dou * (up-low) + low;
}



realvector::realvector(const size_t dim) :
    std::valarray<double>(dim)
{
}




realvector::realvector(const size_t dim, const double v1, ...) :
std::valarray<double>(dim)
{
    std::vector<double> tmp;
    va_list components;
    va_start(components, v1);
    
    for (size_t i=0; i<dim; i++)
    {
        (*this)[i] = va_arg(components, double);
    }
    
    va_end(components);
}




double realvector::dot(const realvector& v) const
{
    assert( v.size() == this->size() );
    double sum=0.0;
    for (size_t i=0; i<size(); i++) sum += (*this)[i]*v[i];
    return sum;
}




double realvector::norm() const
{
    const double n2 = this->dot(*this);
    return sqrt(n2);
}




void realvector::normalize()
{
    double n = this->norm();
    if (n > 0.0) (*this) /= n;
}




void realvector::normalize(realvector& v)
{
    double n = v.norm();
    assert( n > 0.0 );
    
    v /= n;
}




realvector& realvector::operator=(const realvector &v)
{
    assert( v.size() == this->size() );
    for (size_t i=0; i<size(); i++) (*this)[i] = v[i];
    return *this;
}




realvector& realvector::operator+=(const realvector &v)
{
    assert( v.size() == this->size() );
    for (size_t i=0; i<size(); i++) (*this)[i] += v[i];
    return *this;
}




realvector& realvector::operator-=(const realvector &v)
{
    assert( v.size() == this->size() );
    for (size_t i=0; i<size(); i++) (*this)[i] -= v[i];
    return *this;
}




realvector& realvector::operator*=(const double a)
{
    for (size_t i=0; i<size(); i++) (*this)[i] *= a;
    return *this;
}




realvector& realvector::operator/=(const double a)
{
    for (size_t i=0; i<size(); i++) (*this)[i] /= a;
    return *this;
}



namespace muesli
{
    realvector operator+(const realvector& left, const realvector& right)
    {
        assert( left.size() == right.size() );
        realvector v( left.size() );

        for (size_t i=0; i< left.size(); i++)
            v[i] = left[i] + right[i];

        return v;
    }




    realvector operator-(const realvector &left, const realvector& right)
    {
        assert( left.size() == right.size() );
        realvector v( left.size() );

        for (size_t i=0; i< left.size(); i++)
            v[i] = left[i] - right[i];

        return v;
    }
}



void realvector::print(std::ostream& of) const
{
    of << std::scientific << std::setprecision(6);
    for (size_t i=0; i<size(); i++)
        of << (*this)[i] << "  ";
}



void realvector::setRandom()
{
    for (size_t i=0; i<size(); i++)
        (*this)[i] = randomUniform(0.0, 1.0);
}




void realvector::setZero()
{
    for (size_t i=0; i<size(); i++) (*this)[i] = 0.0;
}

#endif


