text:
// Copyright 2019 Intel Corporation.

#include "plaidml/config.h"

#include <boost/dll/runtime_symbol_info.hpp>
#include <boost/filesystem.hpp>

#include "base/config/config.h"
#include "base/util/env.h"

namespace vertexai {
namespace plaidml {
namespace config {

namespace {

const char* kPlaidMLExperimental = "0"; //"PLAIDML_EXPERIMENTAL";
const char* kPlaidMLDefaultConfig = "PLAIDML_DEFAULT_CONFIG";
const char* kPlaidMLDefaultConfigBasename = "config.json";
const char* kPlaidMLExperimentalConfig = "PLAIDML_EXPERIMENTAL_CONFIG";
const char* kPlaidMLExperimentalConfigBasename = "experimental.json";

boost::filesystem::path GetConfigFilename(const char* env_var, const char* basename) {
  std::string config_filename = vertexai::env::Get(env_var);
  if (!config_filename.empty()) {
    return config_filename;
  }

  // Assuming that the PlaidML library has been properly installed, we
  // should be able to find the config file from the library's path.
  auto path = boost::dll::symbol_location_ptr(&GetConfigFilename).lexically_normal();
  path.remove_filename();  // Trim to "/Library/bin"
  path.remove_filename();  // Trim to "/Library"
  path.remove_filename();  // Trim to root install directory
  path /= "share";
  path /= "plaidml";
  path /= basename;
  return path;
}

}  // namespace

Config Get() {
  Config config;
  //std::string exp = vertexai::env::Get(kPlaidMLExperimental);
  boost::filesystem::path source;
  //if (!exp.empty() && exp != "0") {
  //  source = GetConfigFilename(kPlaidMLExperimentalConfig, kPlaidMLExperimentalConfigBasename);
  //} else {
    source = GetConfigFilename(kPlaidMLDefaultConfig, kPlaidMLDefaultConfigBasename);
  //}
  config.source = source.string();
  std::cout << "config.source=" << config.source << "\n";
  std::ifstream cfs(source.native());
  config.data.assign(std::istreambuf_iterator<char>(cfs), std::istreambuf_iterator<char>());

  return config;
}

}  // namespace config
}  // namespace plaidml
}  // namespace vertexai

ext: cc

max_stars_repo_name: conradjones/plaidml

max_issues_repo_name: conradjones/plaidml

max_forks_repo_name: conradjones/plaidml


text:
#ifndef __background_worker_hpp__
#define __background_worker_hpp__

/*
 * Copyright (c) 2016, Anonymous Institution.
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. 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.
 * 3. Neither the name of the University nor the names of its contributors
 *    may be used to endorse or promote products derived from this software
 *    without specific prior written permission.
 *
 * 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
 * HOLDERS 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.
 */

#include <boost/function.hpp>
#include <boost/unordered_map.hpp>
#include <boost/make_shared.hpp>
#include <boost/bind.hpp>

#include <vector>
#include <string>

#include "work-puller.hpp"

class BackgroundWorker {
 public:
  explicit BackgroundWorker(const boost::shared_ptr<WorkPuller>& work_puller)
    : work_puller(work_puller) {}
  typedef boost::function<void (std::vector<ZmqPortableBytes>&)>
    WorkerCallback;
  static const uint32_t STOP_CMD = 0;
  int add_callback(uint32_t cmd, WorkerCallback callback);
  void pull_work_loop();
  void operator()();
      /* This function will be used as the entry point of a boost thread */

 private:
  boost::shared_ptr<WorkPuller> work_puller;
  boost::unordered_map<uint32_t, WorkerCallback> callback_map;
};

#endif  // defined __background_worker_hpp__

ext: hpp


text:
#pragma once
#include <Eigen/Dense>
#include <Eigen/Sparse>

#include <catch2/catch.hpp>

inline void make_first_positive(Eigen::VectorXi& v)
{
    for(uint32_t n = 0; n < v.size(); ++n)
    {
        if(v(n) == 0)
        {
            continue;
        }
        if(v(n) < 0)
        {
            v *= -1;
            return;
        }
        else
        {
            return;
        }
    }
}

struct hash_vector
{
    size_t operator()(const Eigen::VectorXi& v) const
    {
        std::size_t seed = v.size();
        for(uint32_t n = 0; n < v.size(); ++n)
        {
            seed ^= abs(v(n)) + 0x9e3779b9 + (seed << 6U) + (seed >> 2U);
        }
        return seed;
    }
};

inline int powi(int base, unsigned int exp)
{
    int res = 1;
    while(exp != 0U)
    {
        if((exp & 1U) != 0U)
        {
            res *= base;
        }
        exp >>= 1U;
        base *= base;
    }
    return res;
}

template<typename UINT, template<typename> class Basis>
Eigen::MatrixXd basisMatrix(const Basis<UINT>& basis)
{
    Eigen::MatrixXd res = Eigen::MatrixXd::Zero(1U << basis.getN(), basis.getDim());
    for(unsigned int n = 0; n < basis.getDim(); ++n)
    {
        auto bvec = basis.basisVec(n);
        for(const auto p : bvec)
        {
            res(p.first, n) = p.second;
        }
    }
    return res;
}

template<typename Basis> Eigen::VectorXd flip(Basis&& basis, const Eigen::VectorXd& r)
{
    Eigen::VectorXd res(r.size());
    for(int i = 0; i < r.size(); i++)
    {
        res(basis.flip(i)) = r(i);
    }
    return res;
}

inline void TestBasisMatrix(const Eigen::MatrixXd& r)
{
    using Catch::WithinAbs;
    Eigen::MatrixXd id = Eigen::MatrixXd::Identity(r.cols(), r.cols());
    REQUIRE_THAT((r.transpose() * r - id).cwiseAbs().maxCoeff(), WithinAbs(0.0, 1e-8));
}

Eigen::SparseMatrix<double> getSX();
Eigen::SparseMatrix<std::complex<double>> getSY();
Eigen::SparseMatrix<double> getSZ();
Eigen::SparseMatrix<double> getSXXYY();
Eigen::SparseMatrix<double> getSXX();
Eigen::SparseMatrix<double> getSYY();
Eigen::SparseMatrix<double> getSZZ();

max_stars_repo_name: cecri/ExactDiagonalization

max_issues_repo_name: cecri/ExactDiagonalization

max_forks_repo_name: cecri/ExactDiagonalization


text:
/**
 * @file    MatrixPseudoInverse.cpp
 * @brief   Functions for matrix pseudo inverse, and matrix square root.
 * @author  Jianzhu Huai
 */

#ifndef INCLUDE_OKVIS_MATRIX_PSEUDO_INVERSE_HPP
#define INCLUDE_OKVIS_MATRIX_PSEUDO_INVERSE_HPP

#include <Eigen/Core>
#include <Eigen/Eigenvalues>

#include <okvis/assert_macros.hpp>

namespace okvis {

/**
 * \brief Compute the square root of a matrix using the LDLt decomposition for
 * square, positive semidefinite matrices
 *
 * To reconstruct the input matrix, \f$ \mathbf A \f$, from the returned matrix,
 * \f$ \mathbf S \f$, use, \f$ \mathbf A = \mathbf S \mathbf S^T \f$.
 *
 *
 * @param inMatrix      The square matrix whose square root should be computed.
 * @param outMatrixSqrt The output square root.
 */
template <typename DERIVED1, typename DERIVED2>
/*Eigen::ComputationInfo*/ void
computeMatrixSqrt(const Eigen::MatrixBase<DERIVED1> &inMatrix,
                  const Eigen::MatrixBase<DERIVED2> &outMatrixSqrt) {
  OKVIS_ASSERT_EQ_DBG(std::runtime_error, inMatrix.rows(), inMatrix.cols(),
                      "This method is only valid for square input matrices");

  DERIVED2 &result = const_cast<DERIVED2 &>(outMatrixSqrt.derived());

  // This is tricky. Using the output matrix type causes the input matrix
  // type to be upgraded to a real numeric matrix. This is useful if,
  // for example, the inMatrix is something like Eigen::Matrix3d::Identity(),
  // which is not an actual matrix. Using DERIVED1 as the template argument
  // in that case will cause a firestorm of compiler errors.
  Eigen::LDLT<DERIVED2> ldlt(inMatrix.derived());
  result = ldlt.matrixL();
  result = ldlt.transpositionsP().transpose() * result;
  result *= ldlt.vectorD().array().sqrt().matrix().asDiagonal();

  // return ldlt.info();
}

class MatrixPseudoInverse
{
public:
  OKVIS_DEFINE_EXCEPTION(Exception,std::runtime_error)
  /**
   * @brief Pseudo inversion of a symmetric matrix.
   * @warning   This uses Eigen-decomposition, it assumes the input is symmetric positive semi-definite
   *            (negative Eigenvalues are set to zero).
   * @tparam Derived Matrix type (auto-deducible).
   * @param[in] a Input Matrix
   * @param[out] result Output, i.e. pseudo-inverse.
   * @param[in] epsilon The tolerance.
   * @param[out] rank Optional rank.
   * @return
   */
  template<typename Derived>
  static bool pseudoInverseSymm(
      const Eigen::MatrixBase<Derived>&a,
      const Eigen::MatrixBase<Derived>&result, double epsilon =
          std::numeric_limits<typename Derived::Scalar>::epsilon(), int * rank = 0);

  /**
   * @brief Pseudo inversion and square root (Cholesky decomposition) of a symmetric matrix.
   * @warning   This uses Eigen-decomposition, it assumes the input is symmetric positive semi-definite
   *            (negative Eigenvalues are set to zero). Also if the input is positive semi-definite,
   *            its zero eigenvalues are at the lower part of the diagonal.
   * @tparam Derived Matrix type (auto-deducible).
   * @param[in] a Input Matrix, \f$A\f$.
   * @param[out] result Output, \f$L\f$, i.e. the Cholesky decomposition of a pseudo-inverse, \f$A^{-1} = L L^*\f$.
   * @param[in] epsilon The tolerance.
   * @param[out] rank The rank, if of interest.
   * @return
   */
  template<typename Derived>
  static bool pseudoInverseSymmSqrt(
      const Eigen::MatrixBase<Derived>&a,
      const Eigen::MatrixBase<Derived>&result, double epsilon =
          std::numeric_limits<typename Derived::Scalar>::epsilon(),
      int* rank = NULL);

  /**
   * @brief Pseudo square root (Cholesky decomposition) of a symmetric matrix.
   * @warning   This uses Eigen-decomposition, it assumes the input is symmetric positive semi-definite
   *            (negative Eigenvalues are set to zero). Also if the input is positive semi-definite,
   *            its zero eigenvalues are at the lower part of the diagonal.
   * @tparam Derived Matrix type (auto-deducible).
   * @param[in] a Input Matrix, \f$A\f$.
   * @param[out] result Output, \f$L\f$, i.e. the Cholesky decomposition of a pseudo-inverse, \f$A = L L^*\f$.
   * @param[in] epsilon The tolerance.
   * @param[out] rank The rank, if of interest.
   * @return
   */
  template<typename Derived>
  static bool pseudoSymmSqrt(
      const Eigen::MatrixBase<Derived>&a,
      const Eigen::MatrixBase<Derived>&result, double epsilon =
          std::numeric_limits<typename Derived::Scalar>::epsilon(),
      int* rank = NULL);

  /**
   * @brief Block-wise pseudo inversion of a symmetric matrix with non-zero diagonal blocks.
   * @warning   This uses Eigen-decomposition, it assumes the input is symmetric positive semi-definite
   *            (negative Eigenvalues are set to zero).
   * @tparam Derived Matrix type (auto-deducible).
   * @tparam blockDim The block size of the diagonal blocks.
   * @param[in] M_in Input Matrix
   * @param[out] M_out Output, i.e. thepseudo-inverse.
   * @param[in] epsilon The tolerance.
   * @return
   */
  template<typename Derived, int blockDim>
  static void blockPinverse(
      const Eigen::MatrixBase<Derived>& M_in,
      const Eigen::MatrixBase<Derived>& M_out, double epsilon =
          std::numeric_limits<typename Derived::Scalar>::epsilon());

  /**
   * @brief Block-wise pseudo inversion and square root (Cholesky decomposition)
   *        of a symmetric matrix with non-zero diagonal blocks.
   * @warning   This uses Eigen-decomposition, it assumes the input is symmetric positive semi-definite
   *            (negative Eigenvalues are set to zero).
   * @tparam Derived Matrix type (auto-deducible).
   * @tparam blockDim The block size of the diagonal blocks.
   * @param[in] M_in Input Matrix
   * @param[out] M_out Output, i.e. the Cholesky decomposition of a pseudo-inverse.
   * @param[in] epsilon The tolerance.
   * @return
   */
  template<typename Derived, int blockDim>
  static void blockPinverseSqrt(
      const Eigen::MatrixBase<Derived>& M_in,
      const Eigen::MatrixBase<Derived>& M_out, double epsilon =
          std::numeric_limits<typename Derived::Scalar>::epsilon());

};

// Pseudo inversion of a symmetric matrix.
// attention: this uses Eigen-decomposition, it assumes the input is symmetric positive semi-definite
// (negative Eigenvalues are set to zero)
template<typename Derived>
bool MatrixPseudoInverse::pseudoInverseSymm(
    const Eigen::MatrixBase<Derived>&a, const Eigen::MatrixBase<Derived>&result,
    double epsilon, int * rank) {

  OKVIS_ASSERT_TRUE_DBG(Exception, a.rows() == a.cols(),
                        "matrix supplied is not quadratic");

  Eigen::SelfAdjointEigenSolver<Derived> saes(a);

  typename Derived::Scalar tolerance = epsilon * a.cols()
      * saes.eigenvalues().array().maxCoeff();

  const_cast<Eigen::MatrixBase<Derived>&>(result) = (saes.eigenvectors())
      * Eigen::VectorXd(
          (saes.eigenvalues().array() > tolerance).select(
              saes.eigenvalues().array().inverse(), 0)).asDiagonal()
      * (saes.eigenvectors().transpose());

  if (rank) {
    *rank = 0;
    for (int i = 0; i < a.rows(); ++i) {
      if (saes.eigenvalues()[i] > tolerance)
        (*rank)++;
    }
  }

  return true;
}

// Pseudo inversion and square root (Cholesky decomposition) of a symmetric matrix.
// attention: this uses Eigen-decomposition, it assumes the input is symmetric positive semi-definite
// (negative Eigenvalues are set to zero)
template<typename Derived>
bool MatrixPseudoInverse::pseudoInverseSymmSqrt(
    const Eigen::MatrixBase<Derived>&a, const Eigen::MatrixBase<Derived>&result,
    double epsilon, int * rank) {

  OKVIS_ASSERT_TRUE_DBG(Exception, a.rows() == a.cols(),
                        "matrix supplied is not quadratic");

  Eigen::SelfAdjointEigenSolver<Derived> saes(a);

  typename Derived::Scalar tolerance = epsilon * a.cols()
      * saes.eigenvalues().array().maxCoeff();

  const_cast<Eigen::MatrixBase<Derived>&>(result) = (saes.eigenvectors())
      * Eigen::VectorXd(
          Eigen::VectorXd(
              (saes.eigenvalues().array() > tolerance).select(
                  saes.eigenvalues().array().inverse(), 0)).array().sqrt())
          .asDiagonal();

  if (rank) {
    *rank = 0;
    for (int i = 0; i < a.rows(); ++i) {
      if (saes.eigenvalues()[i] > tolerance)
        (*rank)++;
    }
  }

  return true;
}

template <typename Derived>
bool MatrixPseudoInverse::pseudoSymmSqrt(
    const Eigen::MatrixBase<Derived> &a,
    const Eigen::MatrixBase<Derived> &result, double epsilon, int *rank) {

  OKVIS_ASSERT_TRUE_DBG(Exception, a.rows() == a.cols(),
                        "matrix supplied is not quadratic");

  Eigen::SelfAdjointEigenSolver<Derived> saes(a);

  typename Derived::Scalar tolerance =
      epsilon * a.cols() * saes.eigenvalues().array().maxCoeff();

  const_cast<Eigen::MatrixBase<Derived> &>(result) =
      (saes.eigenvectors()) *
      Eigen::VectorXd(
          Eigen::VectorXd((saes.eigenvalues().array() > tolerance)
                              .select(saes.eigenvalues().array(), 0))
              .array()
              .sqrt())
          .asDiagonal();

  if (rank) {
    *rank = 0;
    for (int i = 0; i < a.rows(); ++i) {
      if (saes.eigenvalues()[i] > tolerance)
        (*rank)++;
    }
  }

  return true;
}

// Block-wise pseudo inversion of a symmetric matrix with non-zero diagonal blocks.
// attention: this uses Eigen-decomposition, it assumes the input is symmetric positive semi-definite
// (negative Eigenvalues are set to zero)
template<typename Derived, int blockDim>
void MatrixPseudoInverse::blockPinverse(
    const Eigen::MatrixBase<Derived>& M_in,
    const Eigen::MatrixBase<Derived>& M_out, double epsilon) {

  OKVIS_ASSERT_TRUE_DBG(Exception, M_in.rows() == M_in.cols(),
                        "matrix supplied is not quadratic");

  const_cast<Eigen::MatrixBase<Derived>&>(M_out).resize(M_in.rows(),
                                                        M_in.rows());
  const_cast<Eigen::MatrixBase<Derived>&>(M_out).setZero();
  for (int i = 0; i < M_in.cols(); i += blockDim) {
    Eigen::Matrix<double, blockDim, blockDim> inv;
    const Eigen::Matrix<double, blockDim, blockDim> in = M_in
        .template block<blockDim, blockDim>(i, i);
    //const Eigen::Matrix<double,blockDim,blockDim> in1=0.5*(in+in.transpose());
    pseudoInverseSymm(in, inv, epsilon);
    const_cast<Eigen::MatrixBase<Derived>&>(M_out)
        .template block<blockDim, blockDim>(i, i) = inv;
  }
}

// Block-wise pseudo inversion and square root (Cholesky decomposition)
// of a symmetric matrix with non-zero diagonal blocks.
// attention: this uses Eigen-decomposition, it assumes the input is symmetric positive semi-definite
// (negative Eigenvalues are set to zero)
template<typename Derived, int blockDim>
void MatrixPseudoInverse::blockPinverseSqrt(
    const Eigen::MatrixBase<Derived>& M_in,
    const Eigen::MatrixBase<Derived>& M_out, double epsilon) {

  OKVIS_ASSERT_TRUE_DBG(Exception, M_in.rows() == M_in.cols(),
                        "matrix supplied is not quadratic");

  const_cast<Eigen::MatrixBase<Derived>&>(M_out).resize(M_in.rows(),
                                                        M_in.rows());
  const_cast<Eigen::MatrixBase<Derived>&>(M_out).setZero();
  for (int i = 0; i < M_in.cols(); i += blockDim) {
    Eigen::Matrix<double, blockDim, blockDim> inv;
    const Eigen::Matrix<double, blockDim, blockDim> in = M_in
        .template block<blockDim, blockDim>(i, i);
    //const Eigen::Matrix<double,blockDim,blockDim> in1=0.5*(in+in.transpose());
    pseudoInverseSymmSqrt(in, inv, epsilon);
    const_cast<Eigen::MatrixBase<Derived>&>(M_out)
        .template block<blockDim, blockDim>(i, i) = inv;
  }
}

} // namespace okvis
#endif // INCLUDE_OKVIS_MATRIX_PSEUDO_INVERSE_HPP