covariance.hpp

Go to the documentation of this file.

/*
 * Copyright (C) 1998-2018 ALPS Collaboration. See COPYRIGHT.TXT
 * All rights reserved. Use is subject to license terms. See LICENSE.TXT
 * For use in publications, see ACKNOWLEDGE.TXT
 */
#pragma once

#include <alps/alea/core.hpp>
#include <alps/alea/util.hpp>
#include <alps/alea/bundle.hpp>
#include <alps/alea/complex_op.hpp>
#include <alps/alea/computed.hpp>
#include <alps/alea/var_strategy.hpp>

#include <memory>

// Forward declarations

namespace alps { namespace alea {
    template <typename T, typename Str> class cov_data;
    template <typename T, typename Str> class cov_acc;
    template <typename T, typename Str> class cov_result;

    template <typename T> class batch_result;

    template <typename T, typename Str>
    void serialize(serializer &, const std::string &, const cov_result<T,Str> &);

    template <typename T, typename Str>
    void deserialize(deserializer &, const std::string &, cov_result<T,Str> &);

    template <typename T, typename Str>
    std::ostream &operator<<(std::ostream &, const cov_result<T,Str> &);
}}

// Actual declarations

namespace alps { namespace alea {

template <typename T, typename Strategy=circular_var>
class cov_data
{
public:
    typedef typename bind<Strategy, T>::value_type value_type;
    typedef typename bind<Strategy, T>::cov_type cov_type;
    typedef typename eigen<cov_type>::matrix cov_matrix_type;

public:
    cov_data(size_t size);

    void reset();

    size_t size() const { return data_.rows(); }

    size_t count() const { return count_; }

    size_t &count() { return count_; }

    double count2() const { return count2_; }

    double &count2() { return count2_; }

    const column<value_type> &data() const { return data_; }

    column<value_type> &data() { return data_; }

    const cov_matrix_type &data2() const { return data2_; }

    cov_matrix_type &data2() { return data2_; }

    void convert_to_mean();

    void convert_to_sum();

private:
    column<T> data_;
    cov_matrix_type data2_;
    size_t count_;
    double count2_;

    friend class cov_acc<T, Strategy>;
    friend class cov_result<T, Strategy>;
    friend void serialize<>(serializer &, const std::string &, const cov_result<T,Strategy> &);
    friend void deserialize<>(deserializer &, const std::string &, cov_result<T,Strategy> &);
};

template <typename T, typename Strategy>
struct traits< cov_data<T,Strategy> >
{
    typedef Strategy strategy_type;
    typedef typename bind<Strategy, T>::value_type value_type;
    typedef typename bind<Strategy, T>::var_type var_type;
    typedef typename bind<Strategy, T>::cov_type cov_type;
};

extern template class cov_data<double>;
extern template class cov_data<std::complex<double>, circular_var>;
extern template class cov_data<std::complex<double>, elliptic_var>;


template <typename T, typename Strategy=circular_var>
class cov_acc
{
public:
    using value_type = T;
    using var_type = typename bind<Strategy, T>::var_type;
    using cov_type =  typename bind<Strategy, T>::cov_type;
    using cov_matrix_type = typename eigen<cov_type>::matrix;

public:
    cov_acc(size_t size=1, size_t batch_size=1);

    cov_acc(const cov_acc &other);

    cov_acc &operator=(const cov_acc &other);

    void reset();

    void set_size(size_t size);

    void set_batch_size(size_t batch_size);

    bool valid() const { return (bool)store_; }

    size_t size() const { return current_.size(); }

    size_t batch_size() const { return current_.target(); }

    cov_acc& operator<<(const computed<T>& src){ add(src, 1); return *this; }

    cov_acc &operator<<(const cov_result<T,Strategy> &result);

    size_t count() const { return store_->count(); }

    cov_result<T,Strategy> result() const;

    cov_result<T,Strategy> finalize();

    const bundle<value_type> &current() const { return current_; }

    const cov_data<T,Strategy> &store() const { return *store_; }

protected:
    void add(const computed<T> &source, size_t count);

    void add_bundle();

    void finalize_to(cov_result<T,Strategy> &result);

private:
    std::unique_ptr<cov_data<T,Strategy> > store_;
    bundle<value_type> current_;

    friend class batch_result<T>;
};

template <typename T, typename Strategy>
struct traits< cov_acc<T,Strategy> >
{
    typedef Strategy strategy_type;
    typedef typename bind<Strategy, T>::value_type value_type;
    typedef typename bind<Strategy, T>::var_type var_type;
    typedef typename bind<Strategy, T>::cov_type cov_type;
    typedef cov_result<T,Strategy> result_type;
    typedef cov_data<T, Strategy> store_type;
};

extern template class cov_acc<double>;
extern template class cov_acc<std::complex<double>, circular_var>;
extern template class cov_acc<std::complex<double>, elliptic_var>;


template <typename T, typename Strategy=circular_var>
class cov_result
{
public:
    typedef typename bind<Strategy, T>::value_type value_type;
    typedef typename bind<Strategy, T>::var_type var_type;
    typedef typename bind<Strategy, T>::cov_type cov_type;

public:
    cov_result() { }

    cov_result(const cov_data<T,Strategy> &acc_data)
        : store_(new cov_data<T,Strategy>(acc_data))
    { }

    cov_result(const cov_result &other);

    cov_result &operator=(const cov_result &other);

    bool valid() const { return (bool)store_; }

    size_t size() const { return store_->size(); }

    size_t count() const { return store_->count(); }

    double count2() const { return store_->count2(); }

    double batch_size() const { return store_->count2() / store_->count(); }

    double observations() const { return count() / batch_size(); }

    const column<T> &mean() const { return store_->data(); }

    // TODO: this is essentially a weighted variance thing.  The weighted
    // variance differs from the pooled on by a factor.  We should probably
    // split the two things.

    column<var_type> var() const { return batch_size() * store_->data2().diagonal().real(); }

    typename eigen<cov_type>::matrix cov() const { return batch_size() * store_->data2(); }

    column<var_type> stderror() const;

    const cov_data<T,Strategy> &store() const { return *store_; }

    cov_data<T,Strategy> &store() { return *store_; }

    void reduce(const reducer &r) { reduce(r, true, true); }

    friend void serialize<>(serializer &, const std::string &, const cov_result &);

    friend void deserialize<>(deserializer &, const std::string &, cov_result &);

    friend std::ostream &operator<< <>(std::ostream &, const cov_result &);

protected:
    void reduce(const reducer &, bool do_pre_commit, bool do_post_commit);

private:
    std::unique_ptr<cov_data<T,Strategy> > store_;

    friend class cov_acc<T,Strategy>;
};

template <typename T, typename Strategy>
bool operator==(const cov_result<T, Strategy> &r1, const cov_result<T, Strategy> &r2);
template <typename T, typename Strategy>
bool operator!=(const cov_result<T, Strategy> &r1, const cov_result<T, Strategy> &r2)
{
    return !operator==(r1, r2);
}

template<typename T> struct is_alea_acc<cov_acc<T, circular_var>> :
    std::true_type {};
template<typename T> struct is_alea_acc<cov_acc<T, elliptic_var>> :
    std::true_type {};
template<typename T> struct is_alea_result<cov_result<T, circular_var>> :
    std::true_type {};
template<typename T> struct is_alea_result<cov_result<T, elliptic_var>> :
    std::true_type {};

template <typename T, typename Strategy>
struct traits< cov_result<T,Strategy> >
{
    typedef Strategy strategy_type;
    typedef typename bind<Strategy, T>::value_type value_type;
    typedef typename bind<Strategy, T>::var_type var_type;
    typedef typename bind<Strategy, T>::cov_type cov_type;

    const static bool HAVE_MEAN  = true;
    const static bool HAVE_VAR   = true;
    const static bool HAVE_COV   = true;
    const static bool HAVE_TAU   = false;
    const static bool HAVE_BATCH = false;
};

extern template class cov_result<double>;
extern template class cov_result<std::complex<double>, circular_var>;
extern template class cov_result<std::complex<double>, elliptic_var>;

}} /* namespace alps::alea */