primal_dual.cc File Reference

#include <catch2/catch_all.hpp>
#include <random>
#include <cassert>
#include <Eigen/Dense>
#include "sopt/imaging_primal_dual.h"
#include "sopt/primal_dual.h"
#include "sopt/proximal.h"
#include "sopt/types.h"

Include dependency graph for primal_dual.cc:

Go to the source code of this file.

Classes
struct	is_primal_dual_ref< T >

Typedefs
using	Scalar = sopt::t_real
using	t_Vector = sopt::Vector< Scalar >
using	t_Matrix = sopt::Matrix< Scalar >

Functions
	TEST_CASE ("Primal Dual Imaging", "[primaldual]")
	TEST_CASE ("Primal Dual with 0.5 * ||x - x0||_2^2 function", "[primaldual]")
	TEST_CASE ("Check type returned on setting variables")

Variables
constexpr auto	N = 5

Typedef Documentation

Scalar

using Scalar = sopt::t_real

Definition at line 13 of file primal_dual.cc.

t_Matrix

using t_Matrix = sopt::Matrix<Scalar>

Definition at line 15 of file primal_dual.cc.

t_Vector

using t_Vector = sopt::Vector<Scalar>

Definition at line 14 of file primal_dual.cc.

Function Documentation

TEST_CASE() [1/3]

TEST_CASE

(

"Check type returned on setting variables"

)

Definition at line 87 of file primal_dual.cc.

                                                      {
  // Yeah, could be static asserts
  using namespace sopt;
  using namespace sopt::algorithm;
  ImagingPrimalDual<double> pd(Vector<double>::Zero(0));
  CHECK(is_primal_dual_ref<decltype(pd.itermax(500))>::value);
  CHECK(is_primal_dual_ref<decltype(pd.sigma(1))>::value);
  CHECK(is_primal_dual_ref<decltype(pd.tau(1))>::value);
  CHECK(is_primal_dual_ref<decltype(pd.rho(1))>::value);
  CHECK(is_primal_dual_ref<decltype(pd.xi(1))>::value);
  CHECK(is_primal_dual_ref<decltype(pd.regulariser_strength(1e0))>::value);
  CHECK(is_primal_dual_ref<decltype(pd.update_scale(1e0))>::value);
  CHECK(is_primal_dual_ref<decltype(pd.positivity_constraint(true))>::value);
  CHECK(is_primal_dual_ref<decltype(pd.real_constraint(true))>::value);
}

TEST_CASE() [2/3]

TEST_CASE	(	"Primal Dual Imaging"	,
		""	[primaldual]
	)

Definition at line 20 of file primal_dual.cc.

                                                 {
  using namespace sopt;
 
  t_Matrix const mId = t_Matrix::Identity(N, N);
 
  t_Vector target = t_Vector::Random(N);
 
  target = sopt::positive_quadrant(target);
 
  auto const epsilon = target.stableNorm() / 2;
 
  auto primaldual = algorithm::ImagingPrimalDual<Scalar>(target)
                        .l1_proximal_weights(t_Vector::Ones(target.size()))
                        .Phi(mId)
                        .Psi(mId)
                        .itermax(5000)
                        .tau(0.1)
                        .regulariser_strength(0.4)
                        .l2ball_proximal_epsilon(epsilon)
                        .relative_variation(1e-4)
                        .residual_convergence(epsilon);
 
  auto const result = primaldual();
  CHECK((result.x - target).stableNorm() <= Approx(epsilon).margin(1e-10));
  CHECK(result.good);
  primaldual
      .l1_proximal([](t_Vector &output, const t_real &regulariser_strength, const t_Vector &input) {
        output = regulariser_strength * input;
      })
      .l1_proximal_weighted(
          [](t_Vector &output, const Vector<t_real> &regulariser_strength, const t_Vector &input) {
            output = 10 * regulariser_strength.array() * input.array();
          });
  CHECK_THROWS(primaldual());
}

References sopt::epsilon(), N, sopt::algorithm::ImagingPrimalDual< SCALAR >::Phi(), sopt::positive_quadrant(), sopt::algorithm::ImagingPrimalDual< SCALAR >::Psi(), sopt::algorithm::ImagingPrimalDual< SCALAR >::residual_convergence(), and sopt::target().

TEST_CASE() [3/3]

TEST_CASE	(	"Primal Dual with 0.5 * ||x - x0||_2^2 function"	,
		""	[primaldual]
	)

Definition at line 55 of file primal_dual.cc.

                                                                            {
  using namespace sopt;
  t_Vector const target0 = t_Vector::Random(N);
  auto const f = [](t_Vector &out, const t_real regulariser_strength, const t_Vector &x) {
    proximal::id(out, regulariser_strength, x);
  };
  auto const g = proximal::L2Norm<Scalar>();
  const t_Vector x_guess = t_Vector::Random(target0.size());
  const t_Vector res = x_guess - target0;
  auto const convergence = [&target0](const t_Vector &x, const t_Vector &res) -> bool {
    return x.isApprox(target0, 1e-9);
  };
  CAPTURE(target0);
  CAPTURE(x_guess);
  CAPTURE(res);
  auto const pd = algorithm::PrimalDual<Scalar>(f, g, target0)
                      .itermax(3000)
                      .regulariser_strength(0.9)
                      .rho(0.5)
                      .update_scale(0.5)
                      .is_converged(convergence);
  auto const result = pd(std::make_tuple(x_guess, res));
  CAPTURE(result.niters);
  CAPTURE(result.x);
  CAPTURE(result.residual);
  CHECK(result.x.isApprox(target0, 1e-9));
  CHECK(result.good);
  CHECK(result.niters < 200);
}

References sopt::proximal::id(), sopt::algorithm::PrimalDual< SCALAR >::is_converged(), and N.

Variable Documentation

N

constexpr auto N = 5

constexpr

Definition at line 18 of file primal_dual.cc.

Referenced by TEST_CASE().