libStatGen/general/MathVector.cpp

/*
 *  Copyright (C) 2010  Regents of the University of Michigan
 *
 *   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/>.
 */

#include "MathVector.h"
#include "MathMatrix.h"
#include "MathConstant.h"
#include "Sort.h"
#include "Error.h"

#ifdef  _MSC_VER
#define _USE_MATH_DEFINES
#endif

#include <string.h>
#include <math.h>

int Vector::alloc = 32;

void Vector::Init()
{
    dim = size = 0;
    label = "Unknown";
    data = NULL;
}

Vector::~Vector()
{
    // printf(" Deleting vector %s ...\n", (const char *) label);
    if (data != NULL) delete [] data;
}

void Vector::Dimension(int d)
{
    if (d > size)
    {
        if (size < 1024)
        {
            size = (d + alloc) / alloc * alloc;
            double * newData = new double [size];
            if (data != NULL)
            {
                for (int i = 0; i < dim; i++)
                    newData[i] = data[i];
                delete [] data;
            }
            data = newData;
        }
        else
        {
            while (size <= d)
                size *= 2;

            double * newData = new double [size];
            if (data != NULL)
            {
                for (int i = 0; i < dim; i++)
                    newData[i] = data[i];
                delete [] data;
            }
            data = newData;
        }
    }
    dim = d;
}

void Vector::Dimension(int d, double value)
{
    int original = dim;

    Dimension(d);

    for (int i = original; i < dim; i++)
        data[i] = value;
}

void Vector::Negate()
{
    for (int i = 0; i < dim; i++)
        data[i] = -data[i];
}

void Vector::Add(double n)
{
    for (int i = 0; i< dim; i++)
        data[i] += n;
}

void Vector::Multiply(double k)
{
    for (int i = 0; i < dim; i++)
        data[i] *= k;
}

void Vector::Copy(const Vector & v)
{
    Dimension(v.dim);

    if (v.data != NULL)
        for (int i=0; i < dim; i++)
            data[i] = v.data[i];
}

Vector & Vector::operator = (const Vector & rhs)
{
    Copy(rhs);
    return *this;
}

void Vector::Add(Vector & v)
{
    if (dim != v.dim)
        error("Vector::Add - vectors have different dimensions\n"
              "Vectors     - %s [%d] + %s [%d] ",
              (const char *) label, dim, (const char  *) v.label, v.dim);

    for (int i = 0; i < dim; i++)
        data[i] += v.data[i];
}

void Vector::AddMultiple(double k, Vector & v)
{
    if (dim != v.dim)
        error("Vector::AddMultiple - vectors are incompatible\n"
              "Vectors             - %s [%d] + %s [%d] ",
              (const char  *) label, dim, (const char  *) v.label, v.dim);

    for (int i = 0; i < dim; i++)
        data[i] += k * v.data[i];
}


void Vector::Subtract(Vector & v)
{
    if (dim != v.dim)
        error("Vector::Subtract - vectors have different dimensions\n"
              "Vectors          - %s [%d] + %s [%d] ",
              (const char  *) label, dim, (const char  *) v.label, v.dim);

    for (int i = 0; i < dim; i++)
        data[i] -= v.data[i];
}


void Vector::Zero()
{
    for (int i = 0; i < dim; i++)
        data[i] = 0.0;
}

void Vector::Set(double k)
{
    for (int i = 0; i < dim; i++)
        data[i] = k;
}

void Vector::SetMultiple(double k, Vector & v)
{
    Dimension(v.dim);

    for (int i = 0; i < dim; i++)
        data[i] = k * v[i];
}

double Vector::InnerProduct(Vector & v)
{
    if (dim != v.dim)
        error("Vector::InnerProduct - vectors have different dimensions\n"
              "Vectors              - %s[%d] * %s[%d] ",
              (const char  *) label, dim, (const char  *) v.label, v.dim);

    double sum = 0.0;
    for (int i = 0; i < dim; i++)
        sum += data[i] * v.data[i];

    return sum;
}

void Vector::Insert(int n, double value)
{
    Dimension(dim + 1);

    for (int i = dim - 1; i > n; i--)
        data[i] = data[i - 1];
    data[n] = value;
}

void Vector::DeleteDimension(int n)
{
    for (int i = n; i < dim - 1; i++)
        data[i] = data[i + 1];
    dim --;
}

void Vector::Product(Matrix & m, Vector & v)
{
    if (m.cols != v.dim)
        error("Vector::Product - Cannot Multiply Matrix by Vector\n"
              "Vectors         - %s [%d, %d] * %s [%d]\n",
              (const char  *) m.label, m.rows, m.cols,
              (const char  *) v.label, v.dim);

    Dimension(m.rows);
    Zero();

    for (int i = 0; i < m.rows; i++)
        for (int j = 0; j < m.cols; j++)
            data[i] += m[i][j] * v[j];
}

double Vector::Average() const
{
    if (dim == 0)
        error("Average undefined for null vector %s",
              (const char  *) label);

    return Sum() / dim;
}

double Vector::Product() const
{
    double product = 1.0;

    for (int j = 0; j < dim; j++)
        product *= data[j];

    return product;
}

double Vector::Sum() const
{
    double sum = 0.0;

    for (int j=0; j<dim; j++)
        sum += data[j];

    return sum;
}

double Vector::SumSquares() const
{
    double sum = 0.0;

    for (int j=0; j<dim; j++)
        sum += data[j] * data[j];

    return sum;
}

void Vector::AveVar(double & ave, double & var) const
{
    // uses a two pass method to correct for
    // round-off errors

    if (dim == 0)
        error("Average and Variance undefined for null vector %s",
              (const char  *) label);

    double s, ep;

    ave = var = ep = 0.0;

    for (int j=0; j<dim; j++)
        ave += data[j];

    ave /= dim;

    for (int j=0; j<dim; j++)
    {
        s = data[j] - ave;
        ep += s;
        var += s*s;
    }

    if (dim > 1)
        var = (var - ep*ep/dim)/(dim-1);
}

double Vector::Var() const
{
    double mean, var;
    AveVar(mean, var);
    return var;
}

double Vector::StandardDeviation() const
{
    double var = Var();

    if (var < 0.0) var = 0.0;

    return sqrt(var);
}

void Vector::Print(FILE * f, int d)
{
    if (d == -1 || d > dim) d = dim;

    fprintf(f, "%.15s : ", (const char  *) label);
    for (int i = 0; i < d; i++)
        fprintf(f, "%7.3f ", data[i]);
    fprintf(f, "\n");
}

int Vector::CompareDouble(const double * a, const double * b)
{
    if (*a < *b) return -1;
    if (*a > *b) return 1;
    return 0;
}

void Vector::Sort()
{
    QuickSort(data, dim, sizeof(double), COMPAREFUNC CompareDouble);
}

void Vector::Sort(Vector & freeRider)
{
    QuickSort2(data, freeRider.data, dim, sizeof(double),
               COMPAREFUNC CompareDouble);
}

int Vector::BinarySearch(double element)
{
    void * pointer = ::BinarySearch
                     (&element, data, dim, sizeof(double), COMPAREFUNC CompareDouble);

    if (pointer == NULL)
        return -1;

    return ((double *) pointer) - data;
}

void Vector::RemoveDuplicates()
{
    int out = 0;

    for (int in = 1; in < Length(); in++)
        if (data[in] != data[out])
            data[++out] = data[in];

    Dimension(out + 1);
}

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

    for (int i = 0; i < dim; i++)
        if (data[i] != rhs[i])
            return false;
    return true;
}

// These functions are useful for simulation
//

int Vector::CountIfGreater(double threshold) const
{
    int count = 0;

    for (int i = 0; i < dim; i++)
        if (data[i] > threshold)
            count++;

    return count;
}

int Vector::CountIfGreaterOrEqual(double treshold) const
{
    int count = 0;

    for (int i = 0; i < dim; i++)
        if (data[i] >= treshold)
            count++;

    return count;
}

// Min and max functions
//

double Vector::Min() const
{
    if (dim == 0)
        return 0.0;

    double min = data[0];

    for (int i = 1; i < dim; i++)
        if (data[i] < min)
            min = data[i];

    return min;
}

double Vector::Max() const
{
    if (dim == 0)
        return 0.0;

    double max = data[0];

    for (int i = 1; i < dim; i++)
        if (data[i] > max)
            max = data[i];

    return max;
}

// Push and Pop functions for using vector as a stack
//

void Vector::Push(double value)
{
    Dimension(dim + 1);
    data[dim - 1] = value;
}

void Vector::Stack(const Vector & v)
{
    int end = dim;

    Dimension(dim + v.dim);

    for (int i = 0; i < v.dim; i++)
        data[i + end] = v[i];
}

// Check if all values are in ascending or descending order
//

bool Vector::isAscending()
{
    for (int i = 1; i < dim; i++)
        if (data[i] < data[i - 1])
            return false;
    return true;
}

bool Vector::isDescending()
{
    for (int i = 1; i < dim; i++)
        if (data[i] > data[i - 1])
            return false;
    return true;
}

// VectorFunc class
//

VectorFunc::VectorFunc()
{
    f = NULL;
}

VectorFunc::VectorFunc(double(*func)(Vector &))
{
    f = func;
}

double VectorFunc::Evaluate(Vector & v)
{
    return f(v);
}

#ifndef   M_SQRT2
#define   M_SQRT2        1.41421356
#endif

#define   MAXROUNDS      10
#define   SQRT_HALF      (1.0/M_SQRT2)
#define   TWO            (M_SQRT2 * M_SQRT2)

void VectorFunc::Derivative(Vector & x, Vector & d, double h_start)
{
    double a[MAXROUNDS][MAXROUNDS];

    // Calculate derivatives along each direction ...
    for (int k = 0; k < x.dim; k++)
    {
        double left, right;
        double save_x = x[k];
        double h = h_start;

        // Evaluate function to the left of x along direction k
        x[k] = save_x - h;
        left  = Evaluate(x);

        // Initialize or update dfmin if appropriate...
        if (k == 0 || left < dfmin)
            dfmin = left, dpmin = x;

        // Evaluate function to the right of x along direction k
        x[k] = save_x + h;
        right = Evaluate(x);

        // Update dfmin
        if (right < dfmin)
            dfmin = left, dpmin = x;

        // Initial crude estimate
        a[0][0] = (right - left) / (2.0 * h);

        // Initial guess of error is large
        double err = 1e30;

        // At each round, update Neville tableau with smaller stepsize and higher
        // order extrapolation ...
        for (int i = 1; i < MAXROUNDS; i++)
        {
            // Decrease h
            h *= SQRT_HALF;

            // Re-evaluate function and update dfmin as required
            x[k]  = save_x - h;
            left  = Evaluate(x);
            if (left < dfmin) dfmin = left, dpmin = x;
            x[k]  = save_x + h;
            right = Evaluate(x);
            if (right < dfmin) dfmin = right, dpmin = x;

            // Improved estimate of derivative
            a[0][i] = (right - left) / (2.0 * h);

            // Calculate extrapolations of various orders ...
            double factor = TWO;

            for (int j = 1; j <= i; j++)
            {
                a[j][i] = (a[j-1][i] * factor - a[j-1][i-1])/(factor - 1.0);

                factor *= TWO;

                double error = max(fabs(a[j][i] - a[j-1][i]), fabs(a[j][i] - a[j-1][i-1]));

                // Did we improve solution?
                if (error < err)
                {
                    err = error;
                    d[k] = a[j][i];
                }
            }

            // Stop if solution is deteriorating ...
            if (fabs(a[i][i] - a[i-1][i-1]) >= 2.0 * err)
            {
                x[k] = save_x;
                break;
            }
        }

        x[k] = save_x;
    }
}

int Vector::SafeCount() const
{
    int nonMissing = dim;

    for (int i = 0; i < dim; i++)
        if (data[i] == _NAN_)
            nonMissing--;

    return nonMissing;
}

double Vector::SafeMin() const
{
    double min = _NAN_;
    int i;

    for (i = 0; i < dim; i++)
        if (data[i] != _NAN_)
        {
            min = data[i];
            break;
        }

    for (; i < dim; i++)
        if (data[i] != _NAN_ && data[i] < min)
            min = data[i];

    return min;
}

double Vector::SafeMax() const
{
    double max = _NAN_;
    int i;

    for (i = 0; i < dim; i++)
        if (data[i] != _NAN_)
        {
            max = data[i];
            break;
        }

    for (; i < dim; i++)
        if (data[i] != _NAN_ && data[i] > max)
            max = data[i];

    return max;
}

void Vector::Reverse()
{
    for (int i = 0, j = dim - 1; i < j; i++, j--)
        Swap(i, j);
}

void Vector::InsertInSortedList(int value)
{
    // Skip through large elements
    int pos = dim - 1;

    while (pos >= 0 && data[pos] > value)
        pos--;

    // If the value is already in the list, we are done
    if (pos >= 0 && data[pos] == value)
        return;

    // Otherwise we need to grow array
    Dimension(dim + 1);

    // And then shift larger elements to the right
    pos++;
    for (int i = dim - 1; i > pos; i--)
        data[i] = data[i - 1];

    data[pos] = value;
}

void Vector::Swap(Vector & rhs)
{
    double * temp = rhs.data;
    rhs.data = data;
    data = temp;

    int swap = rhs.dim;
    rhs.dim = dim;
    dim = swap;

    swap = rhs.size;
    rhs.size = size;
    size = swap;
}

double Vector::Average(double returnIfNull)
{
    if (Length() == 0)
        return returnIfNull;

    return Average();
}

double Vector::Var(double returnIfNull)
{
    if (Length() == 0)
        return returnIfNull;

    return Var();
}

double Vector::StandardDeviation(double returnIfNull)
{
    if (Length() == 0)
        return returnIfNull;

    return StandardDeviation();
}