//
// MathVector.cpp
//
// This file contains function definitions for the Vector class.
// This code is part of the Bioinformatics Template Library (BTL).
//
// Copyright (C) 1997 Birkbeck College, Malet Street, London WC1E 7HX, U.K.
// (classlib@mail.cryst.bbk.ac.uk)
//
// This library is free software; you can
// redistribute it and/or modify it under the terms of
// the GNU Library General Public License as published by the Free
// Software Foundation; either version 2 of the License, or (at your
// option) any later version.  This library 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 Library General Public License for more details.
// You should have received a copy of the GNU Library General Public
// License along with this library; if not, write to the Free Software
// Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
////////////////////////////////////////////////////////////////////////////

#include "MathVector.h"
#include "NumericLimits.h"
#include "NumericUtilities.h"
#include <iomanip.h>

//.............................................................................
// Construct from an array of reals ( given value_type* )

Vector::Vector(const value_type* const array, const size_type& p)
{
    for (size_type i=0; i<p; i++)
    	vec.push_back(array[i]);
}
//.............................................................................
// Returns the element, i positions from the beginning of the Vector, 
// in constant time. (i >= 0).
    	    	    
Vector::value_type& Vector::operator[](const size_type& i)
{ 

#if defined(DEBUG_VERSION)
    // Bounds checks
    //
    if (i>=size()) FATAL_ERROR("Index i is out of range");
#endif
    return vec[i];
}

//.............................................................................
// Returns the element, i positions from the beginning of the Vector, 
// in constant time. (i >= 0).
    	    	    
Vector::value_type Vector::operator[](const size_type& i) const
{ 

#if defined(DEBUG_VERSION)
    // Bounds checks
    //
    if (i>=size()) FATAL_ERROR("Index i is out of range");
#endif

    return vec[i];
}

//.............................................................................
// Returns the i-th element in the Vector (i >= 1).
    	    	    
Vector::value_type& Vector::operator()(const size_type& i) 
{

#if defined(DEBUG_VERSION)
    // Bounds checks
    //
    if (i>size() || i<1) FATAL_ERROR("Index i is out of range");
#endif

    return vec[i-1];
}

//************************** VECTOR MODULUS ***********************************
Vector::value_type Vector::Modulus() const
{
    value_type sumsqrs = 0.0;
    for (const_iterator i=vec.begin(); i!=vec.end(); i++)
    	sumsqrs += (*i) * (*i);
    return sqrt(sumsqrs);
}

//************************** LARGEST ELEMENT DIFFERENCE ***********************
// Returns the largest difference between equivalent elements in this and
// Vector v.
    	    	    
Vector::value_type Vector::MaxDifference(const Vector& v) const
{
    if (v.size() != size())
    { 
    	WARNING("The input Vector must be the same size as the this Vector.");
    	return 99999.9;
    }

    double maxdiff = 0.0, diff;
    for (size_type i=0; i<size(); i++)
    {
    	diff = fabs(vec[i] - v[i]);
    	if (diff > maxdiff) maxdiff = diff;
    }
    return maxdiff;    
}

//************************** VECTOR DISTANCE **********************************
Vector::value_type Vector::Distance(const Vector& v) const
{
    if (v.size() != size())
    {
    	WARNING("The input Vector must be the same size as the this Vector.");
    	return 99999.9;
    }

    value_type dist=0.0, diff;
    for (size_type i=0; i<size(); i++)
    {
    	diff = vec[i]-v[i];
    	dist += diff*diff;
    }	
    return (sqrt(dist));
}    
//************************** VECTOR EQUALITY **********************************

bool Vector::operator==(const Vector& v) const
{
    
    if ( size() != v.size() )
    	return false;

    unsigned int error = IsNotEqualTo(v);
    if (error == 0) 
    	return true;
    else
    	return false;
}

//************************** VECTOR ADDITION **********************************
Vector Vector::operator+(const Vector &v) const
{
    if (v.size() != size())
    	FATAL_ERROR("The input Vector must be the same size as this Vector");
    	
    Vector result(*this);
    const_iterator i;
    iterator j;
    for (i=v.vec.begin(), j=result.vec.begin(); i!=v.vec.end(); i++, j++)
    	*j += *i;
    return result;
}

//************************** VECTOR SUBTRACTION *******************************
Vector Vector::operator-(const Vector &v) const
{ 
    if (v.size() != size())
    	FATAL_ERROR("The input Vector must be the same size as this Vector");

    Vector result(*this);
    const_iterator i;
    iterator j;
    for (i=v.vec.begin(), j=result.vec.begin(); i!=v.vec.end(); i++, j++)
    	*j -= *i;
    return result;
}
//************************** VECTOR NEGATION **********************************
// Vector negation (unary minus operator).
    	    	    
Vector Vector::operator-() const
{
    Vector negthis = *this;
    for (iterator i=negthis.begin(); i!=negthis.end(); i++)
    	*i = -(*i);
    return negthis;
}


//************************** VECTOR INCREMENT *********************************
Vector& Vector::operator+=(const Vector &v)
{ 
    if (v.size() != size())
    	FATAL_ERROR("The input Vector must be the same size as this Vector");

    iterator i;
    const_iterator j;
    for (i=vec.begin(), j=v.vec.begin(); i!=vec.end(); i++, j++)
    	*i += *j;
    return *this;
}

//************************** VECTOR DECREMENT *********************************
Vector& Vector::operator-=(const Vector &v)
{ 
    if (v.size() != size())
    	FATAL_ERROR("The input Vector must be the same size as this Vector");

    iterator i;
    const_iterator j;
    for (i=vec.begin(), j=v.vec.begin(); i!=vec.end(); i++, j++)
    	*i -= *j;

    return *this;
}

//*****************************************************************************
// Add a scalar to each element in the Vector.

Vector Vector::operator+(const value_type &s) const
{
    Vector v1 = *this;
    v1 += s;
    return v1;
}    	


//*****************************************************************************
// Subtract a scalar from each element in the Vector.
    	    	    
Vector Vector::operator-(const value_type &s) const
{
    Vector v1 = *this;
    v1 -= s;
    return v1;
}    	


//*****************************************************************************
// Add a scalar to each element in the Vector.
    	    	    
Vector& Vector::operator+=(const value_type &s)
{
    for (size_type i=0; i<size(); i++)
	vec[i] += s;
    return *this;
}    	

//*****************************************************************************
// Subtract a scalar from each element in the Vector.
    	    	    
Vector& Vector::operator-=(const value_type &s)
{
    for (size_type i=0; i<size(); i++)
	vec[i] -= s;
    return *this;
}    	


//************************** SCALAR (DOT) PRODUCT *****************************
Vector::value_type Vector::operator*(const Vector &v) const
{
    if (v.size() != size())
    	FATAL_ERROR("The input Vector must be the same size as this Vector");

    value_type prod = 0.0;
    const_iterator i,j;
    for (i=vec.begin(), j=v.vec.begin(); i!=vec.end(); i++, j++)
    	prod += *i * *j;

    return prod;
}

//************************** VECTOR CROSS PRODUCT ***************************
Vector Vector::operator%(const Vector &v) const
{
    // This works out the cross product of two vector. These vectors must
    // be of size 3.

    if (size() != 3 || v.size() != 3)
    	FATAL_ERROR("Both Vectors must be of size 3 for the cross product to be"
                    " valid.");

    Vector result;
    result[0] = vec[1]*v[2] - vec[2]*v[1];
    result[1] = vec[2]*v[0] - vec[0]*v[2];
    result[2] = vec[0]*v[1] - vec[1]*v[0];
    return result;
}

//************************** VECTOR MULTIPLICATION BY A SCALAR ****************
Vector Vector::operator*(const value_type& y) const
{ 
    Vector result = *this;
    for (iterator i=result.begin(); i!=result.end(); i++)
    	*i *= y;
    return result;
}

//************************** VECTOR SCALING ***********************************
Vector& Vector::operator*=(const value_type& y)
{
    for (iterator i=vec.begin(); i!=vec.end(); i++)
    	*i *= y;
    return *this;
}

//************************** VECTOR DIVISION BY A SCALAR **********************
Vector Vector::operator/(const value_type& y) const
{
#if defined(DEBUG_VERSION)
    if (y == 0.0) WARNING("Divide by zero");
#endif
    
    Vector result = *this; 
    for (iterator i=result.vec.begin(); i!=result.vec.end(); i++)
    	*i /= y;
    return result;
}

//************************** VECTOR NORMALISATION *****************************
Vector& Vector::operator/=(const value_type& y)
{
#if defined(DEBUG_VERSION)
    if (y == 0.0) WARNING("Divide by zero");
#endif
    
    for (iterator i=vec.begin(); i!=vec.end(); i++)
    	*i /= y;
    return *this;
}

//************************** WRITE OUT A VECTOR e.g. cout << v1; **************
ostream& operator<<(ostream &os, Vector const &v)
{
    if (v.vec.empty()) return os;
    os.setf(ios::showpoint);
    os.setf(ios::fixed, ios::floatfield);

    Vector::const_iterator i=v.vec.begin(); 
    os  << "( " << *i++;
    while ( i!=v.vec.end() )
    	os << ", " << *i++;
    os << " )\n";
    return os;
}

//************************** VECTOR COMPARISON ********************************

unsigned int Vector::IsNotEqualTo(const Vector& v) const
{
    
    if ( size() != v.size() )
    {
    	WARNING("The input Vector should be the same size as this Vector.");
    	return NumericLimits<value_type>::digits10();
    }

    unsigned int error, maxError=0;
    for (size_type i=0; i<size(); i++)
    {
    	error = NumericUtilities<value_type>::AreNotEqual(vec[i], v.vec[i]);
    	if (error > maxError) maxError = error;
    }
    	     
    return maxError;
}