/usr/include/sopt/maths.h is in libsopt-dev 2.0.0-4.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 | #ifndef SOPT_MATHS_H
#define SOPT_MATHS_H
#include "sopt/config.h"
#include <algorithm>
#include <complex>
#include <type_traits>
#include <Eigen/Core>
#include "sopt/exception.h"
#include "sopt/types.h"
namespace sopt {
//! Computes the standard deviation of a vector
template <class T>
typename real_type<typename T::Scalar>::type standard_deviation(Eigen::ArrayBase<T> const &x) {
return (x - x.mean()).matrix().stableNorm() / std::sqrt(x.size());
}
//! Computes the standard deviation of a vector
template <class T>
typename real_type<typename T::Scalar>::type standard_deviation(Eigen::MatrixBase<T> const &x) {
return standard_deviation(x.array());
}
//! abs(x) < threshhold ? 0: x - sgn(x) * threshhold
template <class SCALAR>
typename std::enable_if<std::is_arithmetic<SCALAR>::value or is_complex<SCALAR>::value,
SCALAR>::type
soft_threshhold(SCALAR const &x, typename real_type<SCALAR>::type const &threshhold) {
auto const normalized = std::abs(x);
return normalized < threshhold ? SCALAR(0) : (x * (SCALAR(1) - threshhold / normalized));
}
namespace details {
//! Expression to create projection onto positive orthant
template <class SCALAR> class ProjectPositiveQuadrant {
public:
SCALAR operator()(const SCALAR &value) const { return std::max(value, SCALAR(0)); }
};
//! Specialization for complex numbers
template <class SCALAR> class ProjectPositiveQuadrant<std::complex<SCALAR>> {
public:
SCALAR operator()(SCALAR const &value) const { return std::max(value, SCALAR(0)); }
std::complex<SCALAR> operator()(std::complex<SCALAR> const &value) const {
return {operator()(value.real()), SCALAR(0)};
}
};
//! Helper template typedef to instantiate soft_threshhold that takes an Eigen object
template <class SCALAR>
using SoftThreshhold = decltype(
std::bind(soft_threshhold<SCALAR>, std::placeholders::_1, typename real_type<SCALAR>::type(1)));
}
//! Expression to create projection onto positive quadrant
template <class T>
Eigen::CwiseUnaryOp<const details::ProjectPositiveQuadrant<typename T::Scalar>, const T>
positive_quadrant(Eigen::DenseBase<T> const &input) {
typedef details::ProjectPositiveQuadrant<typename T::Scalar> Projector;
typedef Eigen::CwiseUnaryOp<const Projector, const T> UnaryOp;
return UnaryOp(input.derived(), Projector());
}
//! Expression to create soft-threshhold
template <class T>
Eigen::CwiseUnaryOp<const details::SoftThreshhold<typename T::Scalar>, const T>
soft_threshhold(Eigen::DenseBase<T> const &input,
typename real_type<typename T::Scalar>::type const &threshhold) {
typedef typename T::Scalar Scalar;
typedef typename real_type<Scalar>::type Real;
return Eigen::CwiseUnaryOp<const details::SoftThreshhold<typename T::Scalar>, const T>
{input.derived(),
std::bind(soft_threshhold<Scalar>, std::placeholders::_1, Real(threshhold))};
}
//! \brief Expression to create soft-threshhold with multiple parameters
//! \details Operates over a vector of threshholds: ``out(i) = soft_threshhold(x(i), h(i))``
//! Threshhold and input vectors must have the same size and type. The latter condition is enforced
//! by CwiseBinaryOp, unfortunately.
template <class T0, class T1>
typename std::enable_if<std::is_arithmetic<typename T0::Scalar>::value
and std::is_arithmetic<typename T1::Scalar>::value,
Eigen::CwiseBinaryOp<typename T0::Scalar (*)(typename T0::Scalar const &,
typename T1::Scalar const &),
const T0, const T1>>::type
soft_threshhold(Eigen::DenseBase<T0> const &input, Eigen::DenseBase<T1> const &threshhold) {
if(input.size() != threshhold.size())
SOPT_THROW("Threshhold and input should have the same size");
return {input.derived(), threshhold.derived(), soft_threshhold<typename T0::Scalar>};
}
//! \brief Expression to create soft-threshhold with multiple parameters
//! \details Operates over a vector of threshholds: ``out(i) = soft_threshhold(x(i), h(i))``
//! Threshhold and input vectors must have the same size and type. The latter condition is enforced
//! by CwiseBinaryOp, unfortunately. So we cast threshhold from real to complex and back.
template <class T0, class T1>
typename std::
enable_if<is_complex<typename T0::Scalar>::value
and std::is_arithmetic<typename T1::Scalar>::value,
Eigen::CwiseBinaryOp<
typename T0::Scalar (*)(typename T0::Scalar const &, typename T0::Scalar const &),
const T0,
decltype(std::declval<const T1>().template cast<typename T0::Scalar>())>>::type
soft_threshhold(Eigen::DenseBase<T0> const &input, Eigen::DenseBase<T1> const &threshhold) {
if(input.size() != threshhold.size())
SOPT_THROW("Threshhold and input should have the same size: ") << threshhold.size() << " vs "
<< input.size();
typedef typename T0::Scalar Complex;
auto const func
= [](Complex const &x, Complex const &t) -> Complex { return soft_threshhold(x, t.real()); };
return {input.derived(), threshhold.derived().template cast<Complex>(), func};
}
//! Computes weighted L1 norm
template <class T0, class T1>
typename real_type<typename T0::Scalar>::type
l1_norm(Eigen::ArrayBase<T0> const &input, Eigen::ArrayBase<T1> const &weights) {
if(weights.size() == 1)
return input.cwiseAbs().sum() * std::abs(weights(0));
return (input.cwiseAbs() * weights).real().sum();
}
//! Computes weighted L1 norm
template <class T0, class T1>
typename real_type<typename T0::Scalar>::type
l1_norm(Eigen::MatrixBase<T0> const &input, Eigen::MatrixBase<T1> const &weight) {
return l1_norm(input.array(), weight.array());
}
//! Computes L1 norm
template <class T0>
typename real_type<typename T0::Scalar>::type l1_norm(Eigen::ArrayBase<T0> const &input) {
return input.cwiseAbs().sum();
}
//! Computes L1 norm
template <class T0>
typename real_type<typename T0::Scalar>::type l1_norm(Eigen::MatrixBase<T0> const &input) {
return l1_norm(input.array());
}
//! Computes weighted L2 norm
template <class T0, class T1>
typename real_type<typename T0::Scalar>::type
l2_norm(Eigen::ArrayBase<T0> const &input, Eigen::ArrayBase<T1> const &weights) {
if(weights.size() == 1)
return input.matrix().stableNorm() * std::abs(weights(0));
return (input * weights).matrix().stableNorm();
}
//! Computes weighted L2 norm
template <class T0, class T1>
typename real_type<typename T0::Scalar>::type
l2_norm(Eigen::MatrixBase<T0> const &input, Eigen::MatrixBase<T1> const &weights) {
return l2_norm(input.derived().array(), weights.derived().array());
}
namespace details {
//! Greatest common divisor
inline t_int gcd(t_int a, t_int b) { return b == 0 ? a : gcd(b, a % b); }
}
} /* sopt */
#endif
|