covariance.hpp

Go to the documentation of this file.

//============================================================================
// MCKL/include/mckl/utility/covariance.hpp
//----------------------------------------------------------------------------
// MCKL: Monte Carlo Kernel Library
//----------------------------------------------------------------------------
// Copyright (c) 2013-2018, Yan Zhou
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
//   Redistributions of source code must retain the above copyright notice,
//   this list of conditions and the following disclaimer.
//
//   Redistributions in binary form must reproduce the above copyright notice,
//   this list of conditions and the following disclaimer in the documentation
//   and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS AS IS
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//============================================================================

#ifndef MCKL_UTILITY_COVARIANCE_HPP
#define MCKL_UTILITY_COVARIANCE_HPP

#include <mckl/internal/common.hpp>
#include <mckl/internal/cblas.hpp>

namespace mckl {

template <typename RealType = double>
class Covariance
{
    static_assert(internal::is_blas_floating_point<RealType>::value,
        "**Covariance** used with RealType other than float or double");

  public:
    using result_type = RealType;

    void operator()(MatrixLayout layout, std::size_t n, std::size_t p,
        const result_type *x, const result_type *w, result_type *mean,
        result_type *cov, MatrixLayout cov_layout = RowMajor,
        bool cov_upper = false, bool cov_packed = false)
    {
        if (n * p == 0) {
            return;
        }

        if (x == nullptr) {
            return;
        }

        if (mean == nullptr && cov == nullptr) {
            return;
        }

        internal::size_check<MCKL_BLAS_INT>(p, "Covariance::operator()");
        internal::size_check<MCKL_BLAS_INT>(n, "Covariance::operator()");

        result_type sw = w == nullptr ?
            static_cast<result_type>(n) :
            std::accumulate(w, w + n, const_zero<result_type>());
        mean_.resize(p);
        if (w == nullptr) {
            if (layout == RowMajor) {
                std::fill(mean_.begin(), mean_.end(), 0);
                for (std::size_t i = 0; i != n; ++i) {
                    add(p, x + i * p, mean_.data(), mean_.data());
                }
            } else {
                for (std::size_t i = 0; i != p; ++i) {
                    mean_[i] = std::accumulate(
                        x + i * n, x + (i + 1) * n, const_zero<result_type>());
                }
            }
        } else {
            mean_init(layout, n, p, x, w);
        }
        div(p, mean_.data(), sw, mean_.data());
        if (mean != nullptr) {
            std::memcpy(
                mean, mean_.data(), sizeof(result_type) * mean_.size());
        }
        if (cov == nullptr) {
            return;
        }

        result_type sw2 =
            w == nullptr ? static_cast<result_type>(n) : swsqr(n, w);
        result_type B = sw / (sw * sw - sw2);
        result_type BW = B * sw;
        cov_.resize(p * p);
        std::fill(cov_.begin(), cov_.end(), 0);
        cov_init(layout, p, static_cast<result_type *>(nullptr));
        if (w == nullptr) {
            cov_update(layout, n, p, x, B, BW);
        } else {
            wsqrt_.resize(n);
            x_.resize(n * p);
            sqrt(n, w, wsqrt_.data());
            if (layout == RowMajor) {
                for (std::size_t i = 0; i != n; ++i) {
                    mul(p, x + i * p, wsqrt_[i], x_.data() + i * p);
                }
            } else {
                for (std::size_t i = 0; i != p; ++i) {
                    mul(n, x + i * n, wsqrt_.data(), x_.data() + i * n);
                }
            }
            cov_update(layout, n, p, x_.data(), B, BW);
        }
        cov_pack(p, cov, layout, cov_layout, cov_upper, cov_packed);
    }

  private:
    Vector<result_type> mean_;
    Vector<result_type> cov_;
    Vector<result_type> wsqrt_;
    Vector<result_type> x_;

    void mean_init(MatrixLayout layout, std::size_t n, std::size_t p,
        const float *x, const float *w)
    {
        internal::cblas_sgemv(layout == RowMajor ? internal::CblasRowMajor :
                                                   internal::CblasColMajor,
            internal::CblasTrans, static_cast<MCKL_BLAS_INT>(n),
            static_cast<MCKL_BLAS_INT>(p), 1.0, x,
            static_cast<MCKL_BLAS_INT>(layout == RowMajor ? p : n), w, 1, 0.0,
            mean_.data(), 1);
    }

    void mean_init(MatrixLayout layout, std::size_t n, std::size_t p,
        const double *x, const double *w)
    {
        internal::cblas_dgemv(layout == RowMajor ? internal::CblasRowMajor :
                                                   internal::CblasColMajor,
            internal::CblasTrans, static_cast<MCKL_BLAS_INT>(n),
            static_cast<MCKL_BLAS_INT>(p), 1.0, x,
            static_cast<MCKL_BLAS_INT>(layout == RowMajor ? p : n), w, 1, 0.0,
            mean_.data(), 1);
    }

    static float swsqr(std::size_t n, const float *w)
    {
        return internal::cblas_sdot(static_cast<MCKL_BLAS_INT>(n), w, 1, w, 1);
    }

    static double swsqr(std::size_t n, const double *w)
    {
        return internal::cblas_ddot(static_cast<MCKL_BLAS_INT>(n), w, 1, w, 1);
    }

    void cov_init(MatrixLayout layout, std::size_t p, float *)
    {
        internal::cblas_ssyr(layout == RowMajor ? internal::CblasRowMajor :
                                                  internal::CblasColMajor,
            internal::CblasLower, static_cast<MCKL_BLAS_INT>(p), 1,
            mean_.data(), 1, cov_.data(), static_cast<MCKL_BLAS_INT>(p));
    }

    void cov_init(MatrixLayout layout, std::size_t p, double *)
    {
        internal::cblas_dsyr(layout == RowMajor ? internal::CblasRowMajor :
                                                  internal::CblasColMajor,
            internal::CblasLower, static_cast<MCKL_BLAS_INT>(p), 1,
            mean_.data(), 1, cov_.data(), static_cast<MCKL_BLAS_INT>(p));
    }

    void cov_update(MatrixLayout layout, std::size_t n, std::size_t p,
        const float *x, float B, float BW)
    {
        internal::cblas_ssyrk(layout == RowMajor ? internal::CblasRowMajor :
                                                   internal::CblasColMajor,
            internal::CblasLower, internal::CblasTrans,
            static_cast<MCKL_BLAS_INT>(p), static_cast<MCKL_BLAS_INT>(n), B, x,
            static_cast<MCKL_BLAS_INT>(layout == RowMajor ? p : n), -BW,
            cov_.data(), static_cast<MCKL_BLAS_INT>(p));
    }

    void cov_update(MatrixLayout layout, std::size_t n, std::size_t p,
        const double *x, double B, double BW)
    {
        internal::cblas_dsyrk(layout == RowMajor ? internal::CblasRowMajor :
                                                   internal::CblasColMajor,
            internal::CblasLower, internal::CblasTrans,
            static_cast<MCKL_BLAS_INT>(p), static_cast<MCKL_BLAS_INT>(n), B, x,
            static_cast<MCKL_BLAS_INT>(layout == RowMajor ? p : n), -BW,
            cov_.data(), static_cast<MCKL_BLAS_INT>(p));
    }

    void cov_pack(std::size_t p, result_type *cov, MatrixLayout layout,
        MatrixLayout cov_layout, bool cov_upper, bool cov_packed)
    {
        if (layout == RowMajor) {
            for (std::size_t i = 0; i != p; ++i) {
                for (std::size_t j = 0; j != i; ++j) {
                    cov_[j * p + i] = cov_[i * p + j];
                }
            }
        }

        if (layout == ColMajor) {
            for (std::size_t i = 0; i != p; ++i) {
                for (std::size_t j = 0; j != i; ++j) {
                    cov_[i * p + j] = cov_[j * p + i];
                }
            }
        }

        if (!cov_packed) {
            std::memcpy(cov, cov_.data(), sizeof(result_type) * cov_.size());
            return;
        }

        unsigned l = cov_layout == RowMajor ? 0 : 1;
        unsigned u = cov_upper ? 1 : 0;
        unsigned c = (l << 1) + u;
        switch (c) {
            case 0: // Row, Lower, Pack
                for (size_t i = 0; i != p; ++i) {
                    for (std::size_t j = 0; j <= i; ++j) {
                        *cov++ = cov_[i * p + j];
                    }
                }
                break;
            case 1: // Row, Upper, Pack
                for (std::size_t i = 0; i != p; ++i) {
                    for (std::size_t j = i; j != p; ++j) {
                        *cov++ = cov_[i * p + j];
                    }
                }
                break;
            case 2: // Col, Lower, Pack
                for (std::size_t j = 0; j != p; ++j) {
                    for (std::size_t i = j; i != p; ++i) {
                        *cov++ = cov_[j * p + i];
                    }
                }
                break;
            case 3: // Col, Upper, Pack
                for (std::size_t j = 0; j != p; ++j) {
                    for (std::size_t i = 0; i <= j; ++i) {
                        *cov++ = cov_[j * p + i];
                    }
                }
                break;
            default:
                break;
        }
    }
}; // class Covariance

} // namespace mckl

#endif // MCKL_UTILITY_COVARIANCE_HPP