src/modules/tsa/arima.cpp

/* ----------------------------------------------------------------------- *//**
 *
 * @file arima.cpp
 *
 * @brief ARIMA
 *
 * @date Aug 21, 2013
 *//* ----------------------------------------------------------------------- */


#include <dbconnector/dbconnector.hpp>
#include <math.h>
#include <iostream>
#include <algorithm>
#include <functional>
#include <numeric>
#include <sstream>
#include "arima.hpp"

namespace madlib {
namespace modules {
namespace tsa {

using namespace dbal::eigen_integration;
using madlib::dbconnector::postgres::madlib_construct_array;
using madlib::dbconnector::postgres::madlib_construct_md_array;
using namespace std;


AnyType arima_residual::run (AnyType & args)
{
    int32_t distid = args[0].getAs<int32_t>();
    ArrayHandle<double> tvals = args[1].getAs<ArrayHandle<double> >();
    int p = args[2].getAs<int>();
    int d = args[3].getAs<int>();
    int q = args[4].getAs<int>();
    ArrayHandle<double> phi(NULL);
    if(p > 0)
        phi  = args[5].getAs<ArrayHandle<double> >();
    ArrayHandle<double> theta(NULL);
    if(q > 0)
        theta = args[6].getAs<ArrayHandle<double> >();
    double mean = 0.0;
    if(!args[7].isNull())
        mean = args[7].getAs<double>();
    ArrayHandle<double> prez(NULL);
    if(q > 0)
        prez = args[8].getAs<ArrayHandle<double> >();

    int ret_size = static_cast<int>((distid == 1) ? \
            (tvals.size()+d) : (tvals.size()-p));
    // FIXME: construct_array functions circumvent the abstraction layer. These
    // should be replaced with appropriate Allocator:: calls.
    MutableArrayHandle<double> res(
        madlib_construct_array(
            NULL, ret_size, FLOAT8OID, sizeof(double), true, 'd'));

    if (q < 1) {
        // compute the errors
        for(size_t i = p; i < tvals.size(); i++){
            double err = tvals[i] - mean;
            for(int j = 0; j < p; j++)
                err -= phi[j] * (tvals[i - j - 1] - mean);
            // note that for distid = 1, the first p residuals
            // will always be 0
            res[(distid == 1) ? i+d : (i - p)] = err;
        }
        return res;
    } else {
        // in this way we don't need to update prez explicitly
        double * errs = new double[ret_size + q];
        memset(errs, 0, sizeof(double) * (ret_size+q));
        // how to judge ArrayHandle is Null?
        if(distid != 1)
            for(int i = 0; i < q; i++) errs[i] = prez[i];

        // compute the errors
        for(size_t t = p; t < tvals.size(); t++){
            double err = tvals[t] - mean;

            for(int j = 0; j < p; j++)
                err -= phi[j] * (tvals[t - j - 1] - mean);

            for(int j = 0; j < q; j++)
                if(distid == 1)
                    err -= theta[j] * errs[t+q+d-j-1];
                else
                    err -= theta[j] * errs[t-p+q-j-1];
            errs[(distid == 1) ? (t+q+d) : (t - p + q)] = err;
        }
        memcpy(res.ptr(), errs + q, ret_size * sizeof(double));
        delete[] errs;
    }

    return res;
}

// ----------------------------------------------------------------------

static int * diff_coef (int d)
{
    int * coef = new int[d + 1];
    coef[0] = 1, coef[1] = -1;
    for(int i = 2; i <= d; i++)
        coef[i] = 0;

    for (int i = 1; i < d; i++) {
        coef[1 + i] = - coef[i];
        for(int j = i; j >= 1; j--)
            coef[j] = coef[j] - coef[j - 1];
    }

    return coef;
}

// ----------------------------------------------------------------------

AnyType arima_diff::run (AnyType & args)
{
    ArrayHandle<double> tvals = args[0].getAs<ArrayHandle<double> >();
    uint32_t d  = args[1].getAs<uint32_t>();
    int sz = static_cast<int>(tvals.size() - d);
    // FIXME: construct_array functions circumvent the abstraction layer. These
    // should be replaced with appropriate Allocator:: calls.
    MutableArrayHandle<double> diffs(
        madlib_construct_array(
            NULL, sz, FLOAT8OID, sizeof(double), true, 'd'));
    // get diff coef
    int* coef = diff_coef(d);

    // in-place diff
    for(size_t i = tvals.size() - 1; i >= d; i--){
        diffs[i-d] = 0;
        for(size_t j = 0; j <=d; j++)
            diffs[i-d] += coef[j] * tvals[i - j];
        // tvals[i] = diff;
    }

    // // set the first d elements to zero
    // for(int i = 0; i < d; i++)
    //     tvals[i] = 0;

    delete [] coef;
    return diffs;
}

// ----------------------------------------------------------------------

AnyType arima_adjust::run (AnyType & args)
{
    int distid = args[0].getAs<int>();

    if (distid == 1) return args[1];

    ArrayHandle<double> cur_tvals = args[1].getAs<ArrayHandle<double> >();
    ArrayHandle<double> pre_tvals = args[2].getAs<ArrayHandle<double> >();
    int32_t p  = args[3].getAs<int32_t>();

    // FIXME: construct_array functions circumvent the abstraction layer. These
    // should be replaced with appropriate Allocator:: calls.
    // note that curr_tvals.size() could be different from prez_tvals.size()
    MutableArrayHandle<double> res(
        madlib_construct_array(NULL,
                               static_cast<int>(cur_tvals.size() + p),
                               FLOAT8OID, sizeof(double), true, 'd'));

    // fill in the last p values from the previous tvals
    for(int i = 0; i < p; i++)
        res[i] = pre_tvals[pre_tvals.size() - p + i];
    // copy the current tvals
    for(size_t i = 0; i < cur_tvals.size(); i++)
        res[p + i] = cur_tvals[i];

    return res;
}

// ----------------------------------------------------------------------

AnyType arima_lm_delta::run (AnyType & args)
{
    MappedColumnVector jj = args[0].getAs<MappedColumnVector>();
    MappedColumnVector g = args[1].getAs<MappedColumnVector>();
    double u = args[2].getAs<double>();

    size_t l = g.size();
    Matrix m_jj(jj);
    m_jj.resize(l, l);

    Matrix a = m_jj.diagonal().asDiagonal();
    a = m_jj.transpose() + u * a;

    ColumnVector x = a.lu().solve(g);
    return x;
}

// ----------------------------------------------------------------------

AnyType arima_lm::run (AnyType & args)
{
    int distid = args[0].getAs<int>();
    MutableArrayHandle<double> tvals = args[1].getAs<MutableArrayHandle<double> >();
    int p = args[2].getAs<int>();
    int q = args[3].getAs<int>();
    ArrayHandle<double> phi(NULL);
    if (p > 0)
        phi = args[4].getAs<ArrayHandle<double> >();
    ArrayHandle<double> theta(NULL);
    if (q > 0)
        theta = args[5].getAs<ArrayHandle<double> >();

    bool include_mean = true;
    double mean = 0.0;
    if (args[6].isNull())
        include_mean = false;
    else
        mean = args[6].getAs<double>();

    int l = p + q;
    if (include_mean) l++;

    // FIXME: construct_array functions circumvent the abstraction layer. These
    // should be replaced with appropriate Allocator:: calls.
    MutableArrayHandle<double> prez(
        madlib_construct_array(
            NULL, 0, FLOAT8OID, sizeof(double), true, 'd'));
    MutableArrayHandle<double> prej(
        madlib_construct_array(
            NULL, 0, FLOAT8OID, sizeof(double), true, 'd'));
    if (q > 0) {
        if (distid == 1) {
            prez = madlib_construct_array(
                NULL, q, FLOAT8OID, sizeof(double), true, 'd');
            prej = madlib_construct_array(
                NULL, q * l,  FLOAT8OID, sizeof(double), true, 'd');
        } else {
            prez = args[7].getAs<MutableArrayHandle<double> >();
            prej = args[8].getAs<MutableArrayHandle<double> >();
        }
    }

    // minus the mean
    if (include_mean)
        for(size_t i = 0; i < tvals.size(); i++)
            tvals[i] -= mean;

    double z2 = 0;
    // FIXME: construct_array functions circumvent the abstraction layer. These
    // should be replaced with appropriate Allocator:: calls.
    MutableArrayHandle<double> jj(
        madlib_construct_array(
            NULL, l * l, FLOAT8OID, sizeof(double), true, 'd'));
    MutableArrayHandle<double> jz(
        madlib_construct_array(
            NULL, l, FLOAT8OID, sizeof(double), true, 'd'));

    for (size_t t = p; t < tvals.size(); t++) {
        // compute the error and Jacob
        double * jacob = new double[l];
        for (int i = 0; i < l; i++) jacob[i] = 0;
        double err = 0;

        // compute error
        err = tvals[t];
        for (int i = 0; i < p; i++)
            err -= phi[i] * tvals[t - 1 - i];
        for (int i = 0; i < q; i++)
            err -= theta[i] * prez[q - i - 1];

        // compute J
        // compute the partial derivatives over phi
        for (int i = 0; i < p; i++) {
            jacob[i] = tvals[t - i - 1];
            // recursive part
            for(int j = 0; j < q; j++)
                jacob[i] -= theta[j] * prej[(q - j - 1) * l + i];
        }

        // compute the partial derivatives over theta
        for (int i = 0; i < q; i++) {
            jacob[p + i] = prez[q - i - 1];
            // recursive part
            for(int j = 0; j < q; j++)
                jacob[p + i] -= theta[j] * prej[(q - j - 1) * l + p + i];
        }

        // compute the partial derivatives over mean
        if (include_mean) {
            jacob[p + q] = 1;
            for(int i = 0; i < p; i++)
                jacob[p + q] -= phi[i];
            for(int i = 0; i < q; i++)
                jacob[p + q] -= theta[i] * prej[(q - i - 1) * l + p + q];
        }

        // update Z^2
        z2 += err * err;

        if(q != 0){
            // update PreZ
            for(int i = 0; i < q - 1; i++)
                prez[i] = prez[i + 1];
            prez[q-1] = err;

            // update PreJ
            for(int i = 0; i < l * (q - 1); i++)
                prej[i] = prej[i + l];
            for(int i = 0; i < l; i++)
                prej[l * (q - 1) + i] = jacob[i];
        }

        // update J^T * J
        for(int i = 0; i < l; i++)
            for(int j = 0; j < l; j++)
                jj[i * l + j] += jacob[i] * jacob[j];

        // update jz
        for(int i = 0; i < l; i++)
            jz[i] += jacob[i] * err;

        // delete jacob
        if(jacob) delete[] jacob;
    }

    AnyType tuple;
    tuple << z2 << jj << jz << prez << prej;
    return tuple;
}

// ----------------------------------------------------------------------

AnyType arima_lm_result_sfunc::run (AnyType& args)
{
    ArrayHandle<double> jj = args[1].getAs<ArrayHandle<double> >();
    ArrayHandle<double> jz = args[2].getAs<ArrayHandle<double> >();
    double z2 = args[3].getAs<double>();
    int l = static_cast<int>(jz.size());
    int l2 = l*l;

    MutableArrayHandle<double> state(NULL);

    // FIXME: construct_array functions circumvent the abstraction layer. These
    // should be replaced with appropriate Allocator:: calls.
    if (args[0].isNull()) {
        // state[0:(l*l-1)]         - jj
        // state[(l*l):(l*l+l-1)]   - jz
        // state[l*l+l]             - z2
        // state[l*l+l+1]           - l
        state = madlib_construct_array(NULL, l2 + l + 2,
                                       FLOAT8OID, sizeof(double), true, 'd');

        for (int i = 0; i < l2; i++) state[i] = jj[i];
        for (int i = 0; i < l; i++) state[l2 + i] = jz[i];
        state[l2 + l] = z2;
        state[l2 + l + 1] = l;
    } else {
        state = args[0].getAs<MutableArrayHandle<double> >();
        for (int i = 0; i < l2; i++) state[i] += jj[i];
        for (int i = 0; i < l; i++) state[l2 + i] += jz[i];
        state[l2 + l] += z2;
    }

    return state;
}

// ----------------------------------------------------------------------

AnyType arima_lm_result_pfunc::run (AnyType& args)
{
    if (args[0].isNull() && args[1].isNull()) return args[0];
    if (args[0].isNull()) return args[1].getAs<ArrayHandle<double> >();
    if (args[1].isNull()) return args[0].getAs<ArrayHandle<double> >();

    MutableArrayHandle<double> state1 = args[0].getAs<MutableArrayHandle<double> >();
    ArrayHandle<double> state2 = args[1].getAs<ArrayHandle<double> >();
    for (size_t i = 0; i < state1.size() - 1; i++)
        state1[i] += state2[i];

    return state1;
}

// ----------------------------------------------------------------------

AnyType arima_lm_result_ffunc::run (AnyType& args)
{
    ArrayHandle<double> state = args[0].getAs<ArrayHandle<double> >();
    int l = static_cast<int>(state[state.size() - 1]);

    double z2 = state[state.size() - 2];
    const double* jj = state.ptr();
    const double* jz = state.ptr() + l * l;

    double mx = 0;
    for (int i = 0; i < l; i++) {
        int ll = (l + 1) * i;
        if (jj[ll] > mx)
            mx = jj[ll];
    }

    // FIXME: construct_array functions circumvent the abstraction layer. These
    // should be replaced with appropriate Allocator:: calls.
    MutableArrayHandle<double> arr_jj(
        madlib_construct_array(
            NULL, l*l, FLOAT8OID, sizeof(double), true, 'd'));
    MutableArrayHandle<double> arr_jz(
        madlib_construct_array(
            NULL, l, FLOAT8OID, sizeof(double), true, 'd'));

    memcpy(arr_jj.ptr(), jj, l*l*sizeof(double));
    memcpy(arr_jz.ptr(), jz, l*sizeof(double));

    AnyType tuple;
    tuple << arr_jj << arr_jz << z2 << mx;
    return tuple;
}

// ----------------------------------------------------------------------

static double error(int tid, const double * tvals, int p, int q,
                    const double * phi, const double * theta, const double * prez)
{
    double err = 0.0;
    if(tid > p){
        err = tvals[p];
        for(int i = 0; i < p; i++)
            err -= phi[i] * tvals[p - i - 1];
        for(int i = 0; i < q; i++)
            err -= theta[i] * prez[q - i - 1];
    }

    return err;
}

// ----------------------------------------------------------------------

static double error(int tid, const double * tvals, int p, int q,
                    double * phi, double * theta, const double * prez,
                    double delta, int pos1, int pos2, int sign1, int sign2)
{
    double dmean = 0.0; // change of mean
    // Add delta
    if (pos1 == pos2){
        if (pos1 < p)
            phi[pos1] += sign1 * delta;
        else if (pos1 < p + q)
            theta[pos1-p] += sign1 * delta;
        else
            dmean = sign1 * delta;
    } else {
        if (pos1 < p)
            phi[pos1] += sign1 * delta / 2;
        else if (pos1 < p + q)
            theta[pos1-p] += sign1 * delta / 2;
        else
            dmean = sign1 * delta / 2;

        if (pos2 < p)
            phi[pos2] += sign2 * delta / 2;
        else if (pos2 < p + q)
            theta[pos2-p] += sign2 * delta / 2;
        else
            dmean = sign2 * delta / 2;
    }

    double err = 0.0;
    if (tid > p) {
        err = tvals[p] - dmean;
        for (int i = 0; i < p; i++)
            err -= phi[i] * (tvals[p - i - 1] - dmean);
        for (int i = 0; i < q; i++)
            err -= theta[i] * prez[q - i - 1];
    }

    // Restore the original coefficients
    if (pos1 == pos2) {
        if (pos1 < p)
            phi[pos1] -= sign1 * delta;
        else if (pos1 < p + q)
            theta[pos1-p] -= sign1 * delta;
    } else {
        if (pos1 < p)
            phi[pos1] -= sign1 * delta / 2;
        else if (pos1 < p + q)
            theta[pos1-p] -= sign1 * delta / 2;

        if (pos2 < p)
            phi[pos2] -= sign2 * delta / 2;
        else if (pos2 < p + q)
            theta[pos2-p] -= sign2 * delta / 2;
    }
    return err;
}

// ----------------------------------------------------------------------

static double * update_prez(double * prez, int q, double z)
{
    if (q != 0) {
        for (int i = 1; i < q; i++)
            prez[i - 1] = prez[i];
        prez[q-1] = z;
    }
    return prez;
}

// ----------------------------------------------------------------------

AnyType arima_lm_stat_sfunc::run (AnyType& args)
{
    int distid = args[1].getAs<int>();
    MutableArrayHandle<double> tvals = args[2].getAs<MutableArrayHandle<double> >();
    int p = args[3].getAs<int>();
    int q = args[4].getAs<int>();
    MutableArrayHandle<double> phi = args[5].getAs<MutableArrayHandle<double> >();
    MutableArrayHandle<double> theta = args[6].getAs<MutableArrayHandle<double> >();
    bool include_mean = true;
    double mean = 0.0;
    if (args[7].isNull())
        include_mean = false;
    else
        mean = args[7].getAs<double>();
    double delta = args[8].getAs<double>();

    int l = p + q;
    if(include_mean)
        l++;

    // minus the mean
    if (include_mean)
        for(size_t i = 0; i < tvals.size(); i++)
            tvals[i] -= mean;

    MutableArrayHandle<double> state(NULL);

    // FIXME: construct_array functions circumvent the abstraction layer. These
    // should be replaced with appropriate Allocator:: calls.
    if (args[0].isNull()) {
        // Eqs. (1.2.21, 1.2.22) tells how many Z^2 are needed
        // Also Hessian is symmetric
        int sz = (2 * l * l + 1) * (1 + q) + 4;
        // state[0]                     -- l
        // state[1]                     -- delta
        // state[2]                     -- N
        // state[3:(2*l*l+3)]           -- all the Z^2 needed for numerically differentiation
        // state[(2*l*l+4):((2*l*l+1)*q+(2*l*l+3))]  -- all the PreZ
        // last one   --- p
        state = madlib_construct_array(
            NULL, sz, FLOAT8OID, sizeof(double), true, 'd');
        for (size_t i = 0; i < state.size(); i++) state[i] = 0.;
        state[0] = l;
        state[1] = delta;
        state[2] = 0;
        state[sz-1] = p;
    } else {
        state = args[0].getAs<MutableArrayHandle<double> >();
    }

    // Referring to Eqs. (1.2.21, 1.2.22):
    // If a != b in Eqs. (1.2.21, 1.2.22), we need 4 sum(z^2) for each pair
    // of coef values. [4 * l(l-1)/2 = 2l(l-1)]
    // If a == b, we only need three. One of the three sum(z^2) is for
    // coef without any changes. [ 2l + 1]
    // Also H_{ab} = H_{ba}, so totally we need to measure 2l^2 + 1
    // differernt sum(z^2)'s

    // Comute Z^2 and update the state
    // The one without delta
    int prez_offset = 4 + 2 * l * l;
    for (size_t t = p; t < tvals.size(); t++) {
        int dtid = static_cast<int>((distid == 1) ? (1+t) : (p+1));
        int dtv = static_cast<int>(t - p);
        double * prez = state.ptr() + prez_offset;

        double err = error(dtid, tvals.ptr()+dtv, p, q, phi.ptr(), theta.ptr(), prez);
        state[3] += err * err;
        update_prez(prez, q, err);

        // The others with delta
        int count = 0;
        for (int i = 0; i < l; i++) {
            for (int j = i; j < l; j++) {
                if (i == j) {
                    int sign[][2] = {{1, 1}, {-1, -1}};
                    for(int k = 0; k < 2; k++){
                        prez = state.ptr() + prez_offset + q + (i * 2 + k) * q;
                        err = error(dtid, tvals.ptr()+dtv, p, q, phi.ptr(), theta.ptr(), prez, delta, i, j, sign[k][0], sign[k][1]);
                        state[4 + i * 2 + k] += err * err;
                        update_prez(prez, q, err);
                    }
                } else {
                    int sign[][2] = {{1, 1}, {-1, 1}, {1, -1}, {-1, -1}};
                    for (int k = 0; k < 4; k++) {
                        prez = state.ptr() + prez_offset + (1 + 2 * l) * q + q * (count + k);
                        err = error(dtid, tvals.ptr()+dtv, p, q, phi.ptr(), theta.ptr(), prez, delta, i, j, sign[k][0], sign[k][1]);
                        state[4 + 2 * l + count + k] += err * err;
                        update_prez(prez, q, err);
                    }
                    count += 4;
                }
            }
        }
    }

    if (distid == 1) { state[2] += static_cast<double>(tvals.size()); }
    else { state[2] += static_cast<double>(tvals.size() - p); }

    return state;
}

// ----------------------------------------------------------------------

AnyType arima_lm_stat_ffunc::run (AnyType& args)
{
    // extract residual variance, compute log-likelihood
    ArrayHandle<double> state = args[0].getAs<ArrayHandle<double> >();
    int l = static_cast<int>(state[0]);
    double delta = state[1];
    int N = static_cast<int>(state[2]);
    int p = static_cast<int>(state[state.size() - 1]);

    double z2 = state[3];
    double sigma2 = z2 / (N-p);
    double loglik = -(1 + log(2 * M_PI * sigma2)) * N / 2;
    double delta2 = delta * delta * 2. * z2 / N;

    // compute the Hessian matrix
    double * hessian = new double[l * l];
    int count = 0;
    int offset = 4 + 2 * l;
    for (int i = 0; i < l; i++) {
        for (int j = i; j < l; j++) {
            if(j == i) {
                hessian[l * i + i] = (state[4 + i * 2] - 2 * z2 + state[4 + i * 2 + 1]) / delta2;
            } else {
                hessian[l * i + j] = (state[offset + count] - state[offset + count + 1]
                    - state[offset + count + 2] + state[offset + count + 3]) / delta2;
                hessian[l * j + i] = hessian[l * i + j];
                count += 4;
            }
        }
    }

    // inverse the Hessian matrixand to obtain the std errors
    Eigen::Map<Matrix> m(hessian, l, l);
    SymmetricPositiveDefiniteEigenDecomposition<Matrix> decomposition(
         m.transpose(), EigenvaluesOnly, ComputePseudoInverse);
    ColumnVector diag = decomposition.pseudoInverse().diagonal();

    // FIXME: construct_array functions circumvent the abstraction layer. These
    // should be replaced with appropriate Allocator:: calls.
    MutableArrayHandle<double> std_err(
        madlib_construct_array(
            NULL, l, FLOAT8OID, sizeof(double), true, 'd'));
    for (int i = 0; i < l; i++)
        std_err[i] = sqrt(diag[i]);

    delete [] hessian;

    AnyType tuple;
    tuple << std_err << sigma2 << loglik;
    return tuple;
}

}
}
}