reweighted.h

Go to the documentation of this file.

#ifndef SOPT_REWEIGHTED_H
#define SOPT_REWEIGHTED_H
 
#include "sopt/linear_transform.h"
#include "sopt/types.h"
 
#include <algorithm> // for std::max<>
#include <limits>  // for std::numeric_limits<>
#include <utility> // for std::move<>
 
namespace sopt::algorithm {
 
template <typename ALGORITHM>
class Reweighted;
 
template <typename ALGORITHM>
Reweighted<ALGORITHM> reweighted(ALGORITHM const &algo,
                                 typename Reweighted<ALGORITHM>::t_SetWeights const &set_weights,
                                 typename Reweighted<ALGORITHM>::t_Reweightee const &reweightee);
 
template <typename ALGORITHM>
class Reweighted {
 public:
  using Algorithm = ALGORITHM;
  using Scalar = typename Algorithm::Scalar;
  using Real = typename real_type<Scalar>::type;
  using WeightVector = Vector<Real>;
  using XVector = typename Algorithm::t_Vector;
  using t_IsConverged = ConvergenceFunction<Scalar>;
  using t_Reweightee = std::function<XVector (const Algorithm &, const XVector &)>;
  using t_SetWeights = std::function<void (Algorithm &, const WeightVector &)>;
  using t_DeltaUpdate = std::function<Real (Real)>;
 
  struct ReweightedResult {
    t_uint niters;
    bool good;
    WeightVector weights;
    typename Algorithm::DiagnosticAndResult algo;
    ReweightedResult() : niters(0), good(false), weights(WeightVector::Ones(1)), algo() {}
  };
 
  Reweighted(Algorithm const &algo, t_SetWeights const &setweights, t_Reweightee const &reweightee)
      : algo_(algo),
        setweights_(setweights),
        reweightee_(reweightee),
        itermax_(std::numeric_limits<t_uint>::max()),
        min_delta_(0e0),
        is_converged_(),
        update_delta_([](Real delta) { return 1e-1 * delta; }) {}
 
  Algorithm &algorithm() { return algo_; }
  Algorithm const &algorithm() const { return algo_; }
  Reweighted<Algorithm> &algorithm(Algorithm const &algo) {
    algo_ = algo;
    return *this;
  }
  Reweighted<Algorithm> &algorithm(Algorithm &&algo) {
    algo_ = std::move(algo);
    return *this;
  }
 
  t_SetWeights const &set_weights() const { return setweights_; }
  Reweighted<Algorithm> &set_weights(t_SetWeights const &setweights) const {
    setweights_ = setweights;
    return *this;
  }
  void set_weights(Algorithm &algo, WeightVector const &weights) const {
    return set_weights()(algo, weights);
  }
 
  Reweighted<Algorithm> &reweightee(t_Reweightee const &rw) {
    reweightee_ = rw;
    return *this;
  }
  t_Reweightee const &reweightee() const { return reweightee_; }
  XVector reweightee(XVector const &x) const { return reweightee()(algorithm(), x); }
 
  t_uint itermax() const { return itermax_; }
  Reweighted &itermax(t_uint i) {
    itermax_ = i;
    return *this;
  }
  Real min_delta() const { return min_delta_; }
  Reweighted &min_delta(Real min_delta) {
    min_delta_ = min_delta;
    return *this;
  }
  t_IsConverged const &is_converged() const { return is_converged_; }
  Reweighted &is_converged(t_IsConverged const &convergence) {
    is_converged_ = convergence;
    return *this;
  }
  bool is_converged(XVector const &x) const { return is_converged() ? is_converged()(x) : false; }
 
  template <typename INPUT>
  typename std::enable_if<not(std::is_same<INPUT, typename Algorithm::DiagnosticAndResult>::value or
                              std::is_same<INPUT, ReweightedResult>::value),
                          ReweightedResult>::type
  operator()(INPUT const &input) const;
  ReweightedResult operator()() const;
  ReweightedResult operator()(typename Algorithm::DiagnosticAndResult const &warm) const;
  ReweightedResult operator()(ReweightedResult const &warm) const;
 
  Real update_delta(Real delta) const { return update_delta()(delta); }
  t_DeltaUpdate const &update_delta() const { return update_delta_; }
  Reweighted<Algorithm> update_delta(t_DeltaUpdate const &ud) const { return update_delta_ = ud; }
 
 protected:
  Algorithm algo_;
  t_SetWeights setweights_;
  t_Reweightee reweightee_;
  t_uint itermax_;
  Real min_delta_;
  t_IsConverged is_converged_;
  t_DeltaUpdate update_delta_;
};
 
template <typename ALGORITHM>
template <typename INPUT>
typename std::enable_if<
    not(std::is_same<INPUT, typename ALGORITHM::DiagnosticAndResult>::value or
        std::is_same<INPUT, typename Reweighted<ALGORITHM>::ReweightedResult>::value),
    typename Reweighted<ALGORITHM>::ReweightedResult>::type
Reweighted<ALGORITHM>::operator()(INPUT const &input) const {
  Algorithm algo = algorithm();
  set_weights(algo, WeightVector::Ones(1));
  return operator()(algo(input));
}
 
template <typename ALGORITHM>
typename Reweighted<ALGORITHM>::ReweightedResult Reweighted<ALGORITHM>::operator()() const {
  Algorithm algo = algorithm();
  set_weights(algo, WeightVector::Ones(1));
  return operator()(algo());
}
 
template <typename ALGORITHM>
typename Reweighted<ALGORITHM>::ReweightedResult Reweighted<ALGORITHM>::operator()(
    typename Algorithm::DiagnosticAndResult const &warm) const {
  ReweightedResult result;
  result.algo = warm;
  result.weights = WeightVector::Ones(1);
  return operator()(result);
}
 
template <typename ALGORITHM>
typename Reweighted<ALGORITHM>::ReweightedResult Reweighted<ALGORITHM>::operator()(
    ReweightedResult const &warm) const {
  SOPT_HIGH_LOG("Starting reweighted scheme");
  // Copies inner algorithm, so that operator() can be constant
  Algorithm algo(algorithm());
  ReweightedResult result(warm);
 
  auto delta = std::max(standard_deviation(reweightee(warm.algo.x)), min_delta());
  SOPT_LOW_LOG("-   Initial delta: {}", delta);
  for (result.niters = 0; result.niters < itermax(); ++result.niters) {
    SOPT_LOW_LOG("Reweigting iteration {}/{} ", result.niters, itermax());
    SOPT_LOW_LOG("  - delta: {}", delta);
    result.weights = delta / (delta + reweightee(result.algo.x).array().abs());
    set_weights(algo, result.weights);
    result.algo = algo(result.algo);
    if (is_converged(result.algo.x)) {
      SOPT_MEDIUM_LOG("Reweighting scheme did converge in {} iterations", result.niters);
      result.good = true;
      break;
    }
    delta = std::max(min_delta(), update_delta(delta));
  }
  // result is always good if no convergence function is defined
  if (not is_converged())
    result.good = true;
  else if (not result.good)
    SOPT_ERROR("Reweighting scheme did *not* converge in {} iterations", itermax());
  return result;
}
 
template <typename ALGORITHM>
Reweighted<ALGORITHM> reweighted(ALGORITHM const &algo,
                                 typename Reweighted<ALGORITHM>::t_SetWeights const &set_weights,
                                 typename Reweighted<ALGORITHM>::t_Reweightee const &reweightee) {
  return {algo, set_weights, reweightee};
}
 
template <typename SCALAR>
class ImagingProximalADMM;
template <typename ALGORITHM>
class PositiveQuadrant;
template <typename T>
Eigen::CwiseUnaryOp<const details::ProjectPositiveQuadrant<typename T::Scalar>, const T>
positive_quadrant(Eigen::DenseBase<T> const &input);
 
template <typename SCALAR>
Reweighted<PositiveQuadrant<ImagingProximalADMM<SCALAR>>> reweighted(
    ImagingProximalADMM<SCALAR> const &algo) {
  auto const posq = positive_quadrant(algo);
  using Algorithm = typename std::remove_const<decltype(posq)>::type;
  using RW = Reweighted<Algorithm>;
  auto const reweightee =
      [](Algorithm const &posq, typename RW::XVector const &x) -> typename RW::XVector {
    return posq.algorithm().Psi().adjoint() * x;
  };
  auto const set_weights = [](Algorithm &posq, typename RW::WeightVector const &weights) -> void {
    posq.algorithm().l1_proximal_weights(weights);
  };
  return {posq, set_weights, reweightee};
}
 
template <typename SCALAR>
class PrimalDual;
template <typename ALGORITHM>
class PositiveQuadrant;
template <typename T>
Eigen::CwiseUnaryOp<const details::ProjectPositiveQuadrant<typename T::Scalar>, const T>
positive_quadrant(Eigen::DenseBase<T> const &input);
 
template <typename SCALAR>
Reweighted<PositiveQuadrant<PrimalDual<SCALAR>>> reweighted(PrimalDual<SCALAR> const &algo) {
  auto const posq = positive_quadrant(algo);
  using Algorithm = typename std::remove_const<decltype(posq)>::type;
  using RW = Reweighted<Algorithm>;
  auto const reweightee =
      [](Algorithm const &posq, typename RW::XVector const &x) -> typename RW::XVector {
    return posq.algorithm().Psi().adjoint() * x;
  };
  auto const set_weights = [](Algorithm &posq, typename RW::WeightVector const &weights) -> void {
    posq.algorithm().l1_proximal_weights(weights);
  };
  return {posq, set_weights, reweightee};
}
 
} // namespace sopt::algorithm
#endif