fgt.hpp

Go to the documentation of this file.

// fgt — fast Gauss transforms
// Copyright (C) 2016 Peter J. Gadomski <pete.gadomski@gmail.com>
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 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
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA


#pragma once

#include <cstddef>
#include <memory>

#ifdef _MSC_VER
#pragma warning(push)
#pragma warning(disable : 4365) // unsigned/signed
#pragma warning(disable : 5031) // pragma warning mismatch
#pragma warning(disable : 4820) // padding
#pragma warning(disable : 4626) // assignment operator deleted
#pragma warning(disable : 5027) // move assignment operator deleted
#endif
#include <Eigen/Core>
#ifdef _MSC_VER
#pragma warning(pop)
#endif

namespace fgt {

class fgt_error : public std::runtime_error {
public:
    explicit fgt_error(const std::string what_arg)
        : std::runtime_error(what_arg) {}
};

class ifgt_no_clusters : public fgt_error {
public:
    ifgt_no_clusters()
        : fgt_error("IFGT decided that it didn't need any clusters. These "
                    "parameters cannot be used for IFGT, try another method "
                    "instead.") {}
};

typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
    Matrix;

typedef Eigen::Ref<const Matrix> MatrixRef;

typedef Eigen::VectorXd Vector;

typedef Eigen::Ref<const Vector> VectorRef;

const char* version();

bool with_openmp();

Vector direct(const MatrixRef source, const MatrixRef target, double bandwidth);

Vector direct(const MatrixRef source, const MatrixRef target, double bandwidth,
              const VectorRef weights);

Vector direct_tree(const MatrixRef source, const MatrixRef target,
                   double bandwidth, double epsilon);

Vector direct_tree(const MatrixRef source, const MatrixRef target,
                   double bandwidth, double epsilon, const VectorRef weights);

Vector ifgt(const MatrixRef source, const MatrixRef target, double bandwidth,
            double epsilon);

Vector ifgt(const MatrixRef source, const MatrixRef target, double bandwidth,
            double epsilon, const VectorRef weights);

class Transform {
public:
    Transform(const MatrixRef source, double bandwidth);

    Transform(const Transform&) = delete;
    Transform& operator=(const Transform&) = delete;
    Transform& operator=(Transform&&) = delete;

    virtual ~Transform() {}

    const MatrixRef source() const { return m_source; }
    double bandwidth() const { return m_bandwidth; }

    Vector compute(const MatrixRef target);
    Vector compute(const MatrixRef target, const VectorRef weights);

private:
    virtual Vector compute_impl(const MatrixRef target,
                                const VectorRef weights) const = 0;

    const Matrix m_source;
    double m_bandwidth;
};

class Direct : public Transform {
public:
    Direct(const MatrixRef source, double bandwidth);

    Direct(const Direct&) = delete;
    Direct& operator=(const Direct&) = delete;
    Direct& operator=(Direct&&) = delete;

private:
    virtual Vector compute_impl(const MatrixRef target,
                                const VectorRef weights) const;
};

class DirectTree : public Transform {
public:
    DirectTree(const MatrixRef source, double bandwidth, double epsilon);

    DirectTree(const DirectTree&) = delete;
    DirectTree& operator=(const DirectTree&) = delete;
    DirectTree& operator=(DirectTree&&) = delete;

    virtual ~DirectTree();

    double epsilon() const { return m_epsilon; }

private:
    struct NanoflannTree;

    virtual Vector compute_impl(const MatrixRef target,
                                const VectorRef weights) const;

    double m_epsilon;
    std::unique_ptr<NanoflannTree> m_tree;
};

struct Clustering;

class Ifgt : public Transform {
public:
    Ifgt(const MatrixRef source, double bandwidth, double epsilon);

    Ifgt(const Ifgt&) = delete;
    Ifgt& operator=(const Ifgt&) = delete;
    Ifgt& operator=(Ifgt&&) = delete;

    virtual ~Ifgt();

    double epsilon() const { return m_epsilon; }
    Matrix::Index nclusters() const { return m_nclusters; }
    Matrix::Index truncation_number() const { return m_truncation_number; }
    Matrix::Index p_max_total() const { return m_p_max_total; }

private:
    virtual Vector compute_impl(const MatrixRef target,
                                const VectorRef weights) const;
    Vector compute_monomials(const VectorRef d) const;
    Vector compute_constant_series() const;

    double m_epsilon;
    Matrix::Index m_nclusters;
    std::unique_ptr<Clustering> m_clustering;
    Matrix::Index m_truncation_number;
    Matrix::Index m_p_max_total;
    Vector m_constant_series;
    Vector m_ry_square;
};
}