Kernel.hpp

Go to the documentation of this file.

/*
 * STRUMPACK -- STRUctured Matrices PACKage, Copyright (c) 2014, The
 * Regents of the University of California, through Lawrence Berkeley
 * National Laboratory (subject to receipt of any required approvals
 * from the U.S. Dept. of Energy).  All rights reserved.
 *
 * If you have questions about your rights to use or distribute this
 * software, please contact Berkeley Lab's Technology Transfer
 * Department at TTD@lbl.gov.
 *
 * NOTICE. This software is owned by the U.S. Department of Energy. As
 * such, the U.S. Government has been granted for itself and others
 * acting on its behalf a paid-up, nonexclusive, irrevocable,
 * worldwide license in the Software to reproduce, prepare derivative
 * works, and perform publicly and display publicly.  Beginning five
 * (5) years after the date permission to assert copyright is obtained
 * from the U.S. Department of Energy, and subject to any subsequent
 * five (5) year renewals, the U.S. Government is granted for itself
 * and others acting on its behalf a paid-up, nonexclusive,
 * irrevocable, worldwide license in the Software to reproduce,
 * prepare derivative works, distribute copies to the public, perform
 * publicly and display publicly, and to permit others to do so.
 *
 * Developers: Pieter Ghysels, Francois-Henry Rouet, Xiaoye S. Li.
 *             (Lawrence Berkeley National Lab, Computational Research
 *             Division).
 *
 */
#ifndef STRUMPACK_KERNEL_HPP
#define STRUMPACK_KERNEL_HPP
 
#include "Metrics.hpp"
#include "HSS/HSSOptions.hpp"
#include "dense/DenseMatrix.hpp"
#if defined(STRUMPACK_USE_MPI)
#include "dense/DistributedMatrix.hpp"
#if defined(STRUMPACK_USE_BPACK)
#include "HODLR/HODLROptions.hpp"
#endif
#endif
 
namespace strumpack {
 
  namespace kernel {
 
    template<typename scalar_t> class Kernel {
      using real_t = typename RealType<scalar_t>::value_type;
      using DenseM_t = DenseMatrix<scalar_t>;
      using DenseMW_t = DenseMatrixWrapper<scalar_t>;
#if defined(STRUMPACK_USE_MPI)
      using DistM_t = DistributedMatrix<scalar_t>;
#endif
 
    public:
      Kernel(DenseM_t& data, scalar_t lambda)
        : data_(data), lambda_(lambda) { }
 
      virtual ~Kernel() = default;
 
      std::size_t n() const { return data_.cols(); }
 
      std::size_t d() const { return data_.rows(); }
 
      virtual scalar_t eval(std::size_t i, std::size_t j) const {
        return eval_kernel_function(data_.ptr(0, i), data_.ptr(0, j))
          + ((i == j) ? lambda_ : scalar_t(0.));
      }
 
      void operator()(const std::vector<std::size_t>& I,
                      const std::vector<std::size_t>& J,
                      DenseMatrix<real_t>& B) const {
        assert(B.rows() == I.size() && B.cols() == J.size());
        for (std::size_t j=0; j<J.size(); j++)
          for (std::size_t i=0; i<I.size(); i++) {
            assert(I[i] < n() && J[j] < n());
            B(i, j) = eval(I[i], J[j]);
          }
      }
 
      void operator()(const std::vector<std::size_t>& I,
                      const std::vector<std::size_t>& J,
                      DenseMatrix<std::complex<real_t>>& B) const {
        assert(B.rows() == I.size() && B.cols() == J.size());
        for (std::size_t j=0; j<J.size(); j++)
          for (std::size_t i=0; i<I.size(); i++) {
            assert(I[i] < n() && J[j] < n());
            B(i, j) = eval(I[i], J[j]);
          }
      }
 
      DenseM_t fit_HSS
      (std::vector<scalar_t>& labels, const HSS::HSSOptions<scalar_t>& opts);
 
      std::vector<scalar_t> predict
      (const DenseM_t& test, const DenseM_t& weights) const;
 
#if defined(STRUMPACK_USE_MPI)
      DistM_t fit_HSS
      (const BLACSGrid& grid, std::vector<scalar_t>& labels,
       const HSS::HSSOptions<scalar_t>& opts);
 
      std::vector<scalar_t> predict
      (const DenseM_t& test, const DistM_t& weights) const;
 
#if defined(STRUMPACK_USE_BPACK)
      DenseM_t fit_HODLR
      (const MPIComm& c, std::vector<scalar_t>& labels,
       const HODLR::HODLROptions<scalar_t>& opts);
#endif
#endif
 
      const DenseM_t& data() const { return data_; }
      DenseM_t& data() { return data_; }
 
      std::vector<int>& permutation() { return perm_; }
      const std::vector<int>& permutation() const { return perm_; }
 
      virtual void permute() {}
 
    protected:
      DenseM_t& data_;
      scalar_t lambda_;
      std::vector<int> perm_;
 
      virtual scalar_t eval_kernel_function
      (const scalar_t* x, const scalar_t* y) const = 0;
    };
 
 
    template<typename scalar_t>
    class GaussKernel : public Kernel<scalar_t> {
    public:
      GaussKernel(DenseMatrix<scalar_t>& data, scalar_t h, scalar_t lambda)
        : Kernel<scalar_t>(data, lambda), h_(h) {}
 
    protected:
      scalar_t h_; // kernel width parameter
 
      scalar_t eval_kernel_function
      (const scalar_t* x, const scalar_t* y) const override {
        return std::exp
          (-Euclidean_distance_squared(this->d(), x, y)
           / (scalar_t(2.) * h_ * h_));
      }
    };
 
 
    template<typename scalar_t>
    class LaplaceKernel : public Kernel<scalar_t> {
    public:
      LaplaceKernel(DenseMatrix<scalar_t>& data, scalar_t h, scalar_t lambda)
        : Kernel<scalar_t>(data, lambda), h_(h) {}
 
    protected:
      scalar_t h_; // kernel width parameter
 
      scalar_t eval_kernel_function
      (const scalar_t* x, const scalar_t* y) const override {
        return std::exp(-norm1_distance(this->d(), x, y) / h_);
      }
    };
 
    template<typename scalar_t>
    class ANOVAKernel : public Kernel<scalar_t> {
    public:
      ANOVAKernel
      (DenseMatrix<scalar_t>& data, scalar_t h, scalar_t lambda, int p=1)
        : Kernel<scalar_t>(data, lambda), h_(h), p_(p) {
        assert(p >= 1 && p <= int(this->d()));
      }
 
    protected:
      scalar_t h_; // kernel width parameter
      int p_;      // kernel degree parameter 1 <= p_ <= this->d()
 
      scalar_t eval_kernel_function
      (const scalar_t* x, const scalar_t* y) const override {
        std::vector<scalar_t> Ks(p_), Kss(p_), Kpp(p_+1);
        Kpp[0] = 1;
        for (int j=0; j<p_; j++) Kss[j] = 0;
        for (std::size_t i=0; i<this->d(); i++) {
          scalar_t tmp = std::exp
            (-Euclidean_distance_squared(1, x+i, y+i)
             / (scalar_t(2.) * h_ * h_));
          Ks[0] = tmp;
          Kss[0] += Ks[0];
          for (int j=1; j<p_; j++) {
            Ks[j] = Ks[j-1]*tmp;
            Kss[j] += Ks[j];
          }
        }
        for (int i=1; i<=p_; i++) {
          Kpp[i] = 0;
          for (int s=1; s<=i; s++)
            Kpp[i] += std::pow(-1,s+1)*Kpp[i-s]*Kss[s-1];
          Kpp[i] /= i;
        }
        return Kpp[p_];
      }
    };
 
 
    template<typename scalar_t>
    class DenseKernel : public Kernel<scalar_t> {
    public:
      DenseKernel(DenseMatrix<scalar_t>& data,
                  DenseMatrix<scalar_t>& A, scalar_t lambda)
        : Kernel<scalar_t>(data, lambda), A_(A) {}
 
      scalar_t eval(std::size_t i, std::size_t j) const override {
        return A_(i, j) + ((i == j) ? this->lambda_ : scalar_t(0.));
      }
 
      void permute() override {
        A_.lapmt(this->perm_, true);
        A_.lapmr(this->perm_, true);
      }
 
    protected:
      DenseMatrix<scalar_t>& A_; // kernel matrix
 
      scalar_t eval_kernel_function
      (const scalar_t* x, const scalar_t* y) const override {
        assert(false);
      }
    };
 
 
    enum class KernelType {
      DENSE,    
      GAUSS,    
      LAPLACE,  
      ANOVA     
    };
 
    inline std::string get_name(KernelType k) {
      switch (k) {
      case KernelType::DENSE: return "dense";
      case KernelType::GAUSS: return "Gauss";
      case KernelType::LAPLACE: return "Laplace";
      case KernelType::ANOVA: return "ANOVA";
      default: return "UNKNOWN";
      }
    }
 
    inline KernelType kernel_type(const std::string& k) {
      if (k == "dense") return KernelType::DENSE;
      else if (k == "Gauss") return KernelType::GAUSS;
      else if (k == "Laplace") return KernelType::LAPLACE;
      else if (k == "ANOVA") return KernelType::ANOVA;
      std::cerr << "ERROR: Kernel type not recogonized, "
                << " setting kernel type to Gauss."
                << std::endl;
      return KernelType::GAUSS;
    }
 
    template<typename scalar_t>
    std::unique_ptr<Kernel<scalar_t>> create_kernel
    (KernelType k, DenseMatrix<scalar_t>& data,
     scalar_t h, scalar_t lambda, int p=1) {
      switch (k) {
      // case KernelType::DENSE:
      //   return std::unique_ptr<Kernel<scalar_t>>
      //     (new DenseKernel<scalar_t>(args ...));
      case KernelType::GAUSS:
        return std::unique_ptr<Kernel<scalar_t>>
          (new GaussKernel<scalar_t>(data, h, lambda));
      case KernelType::LAPLACE:
        return std::unique_ptr<Kernel<scalar_t>>
          (new LaplaceKernel<scalar_t>(data, h, lambda));
      case KernelType::ANOVA:
        return std::unique_ptr<Kernel<scalar_t>>
          (new ANOVAKernel<scalar_t>(data, h, lambda, p));
      default:
        return std::unique_ptr<Kernel<scalar_t>>
          (new GaussKernel<scalar_t>(data, h, lambda));
      }
    }
 
  } // end namespace kernel
 
} // end namespace strumpack
 
#endif // STRUMPACK_KERNEL_HPP