Skip to content

Commit cf409dc

Browse files
AlePaluAlePalu
committed
virtually dispatchable distributions (simd_distribution), gsr sets its distribution at runtime, gsr gcv with external caching support, core submodule update
1 parent d04d771 commit cf409dc

File tree

6 files changed

+206
-148
lines changed

6 files changed

+206
-148
lines changed

‎fdaPDE/core‎

‎fdaPDE/src/distributions.h‎

Lines changed: 81 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -93,8 +93,27 @@ inline double inverse_erf(double x, double eps = 1e-10) {
9393
} while (std::fabs(delta) > eps);
9494
return y;
9595
}
96+
97+
} // namespace internals
98+
99+
class simd_distribution {
100+
#ifdef __FDAPDE_HAS_EIGEN__ // SIMD vectorized
101+
public:
102+
using matrix_t = Eigen::Matrix<double, Dynamic, 1>;
103+
constexpr simd_distribution() = default;
104+
105+
virtual matrix_t variance(const matrix_t& x) const = 0;
106+
virtual matrix_t link (const matrix_t& x) const = 0;
107+
virtual matrix_t inv_link(const matrix_t& x) const = 0;
108+
virtual matrix_t der_link(const matrix_t& x) const = 0;
109+
virtual double deviance(const matrix_t& x, const matrix_t& y) const = 0;
110+
#endif
111+
virtual ~simd_distribution() = default;
112+
};
96113

97-
template <typename Distribution_> class distribution_base {
114+
namespace internals {
115+
116+
template <typename Distribution_> class distribution_base : public simd_distribution {
98117
public:
99118
using Distribution = Distribution_;
100119
constexpr distribution_base() noexcept { }
@@ -104,8 +123,7 @@ template <typename Distribution_> class distribution_base {
104123

105124
template <typename T, typename F>
106125
requires(
107-
internals::is_vector_like_v<T> &&
108-
std::is_convertible_v<internals::subscript_result_of_t<T, int>, double>)
126+
internals::is_vector_like_v<T> && std::is_convertible_v<internals::subscript_result_of_t<T, int>, double>)
109127
constexpr std::vector<double> apply_(const T& x, F&& f) const {
110128
std::vector<double> res(x.size());
111129
for (std::size_t i = 0; i < x.size(); ++i) { res[i] = f(x[i]); }
@@ -118,7 +136,7 @@ template <typename Distribution_> class distribution_base {
118136
struct bernoulli_distribution : public internals::distribution_base<std::bernoulli_distribution> {
119137
using result_type = double;
120138
using param_type = double;
121-
private:
139+
protected:
122140
using Base = internals::distribution_base<std::bernoulli_distribution>;
123141
using Base::distr_;
124142
param_type p_ = 0;
@@ -159,13 +177,19 @@ struct bernoulli_distribution : public internals::distribution_base<std::bernoul
159177
constexpr result_type deviance(T x, T y) const {
160178
return almost_zero(y) ? 2 * std::log(1.0 / (1.0 - x)) : 2.0 * std::log(1.0 / x);
161179
}
162-
// SIMD vectorized
163-
using matrix_t = Eigen::Matrix<double, Dynamic, Dynamic>;
164-
matrix_t variance(const matrix_t& x) const { return x.array() * (1 - x.array()); }
165-
matrix_t link (const matrix_t& x) const { return ((1 - x.array()).inverse() * x.array()).log(); }
166-
matrix_t inv_link(const matrix_t& x) const { return (1 + ((-x).array().exp())).inverse(); }
167-
matrix_t der_link(const matrix_t& x) const { return (x.array() * (1 - x.array())).inverse(); }
168-
180+
#ifdef __FDAPDE_HAS_EIGEN__ // SIMD vectorized
181+
using matrix_t = Eigen::Matrix<double, Dynamic, 1>;
182+
matrix_t variance(const matrix_t& x) const override { return x.array() * (1 - x.array()); }
183+
matrix_t link (const matrix_t& x) const override { return ((1 - x.array()).inverse() * x.array()).log(); }
184+
matrix_t inv_link(const matrix_t& x) const override { return (1 + ((-x).array().exp())).inverse(); }
185+
matrix_t der_link(const matrix_t& x) const override { return (x.array() * (1 - x.array())).inverse(); }
186+
double deviance(const matrix_t& x, const matrix_t& y) const override {
187+
fdapde_assert(x.cols() == 1 && y.cols() == 1 && x.rows() == y.rows());
188+
double dev_ = 0;
189+
for (int i = 0, n = x.rows(); i < n; ++i) { dev_ += deviance(x(i, 0), y(i, 0)); }
190+
return dev_;
191+
}
192+
#endif
169193
template <typename T> constexpr auto transform(const T& data) const {
170194
if constexpr (internals::is_eigen_dense_xpr_v<T>) {
171195
return 0.5 * (data.array() + 0.5);
@@ -176,11 +200,11 @@ struct bernoulli_distribution : public internals::distribution_base<std::bernoul
176200
void set_param(param_type p) { p_ = p; }
177201
};
178202

179-
struct rademacher_distribution : public internals::distribution_base<std::bernoulli_distribution> {
203+
struct rademacher_distribution : public bernoulli_distribution {
180204
using result_type = double;
181-
using param_type = double;
205+
using param_type = double;
182206
private:
183-
using Base = internals::distribution_base<std::bernoulli_distribution>;
207+
using Base = bernoulli_distribution;
184208
using Base::distr_;
185209
public:
186210
constexpr rademacher_distribution() noexcept : Base(0.5) { }
@@ -249,13 +273,18 @@ struct poisson_distribution : public internals::distribution_base<std::poisson_d
249273
constexpr result_type deviance(T x, T y) const {
250274
return y > 0 ? y * std::log(y / x) - (y - x) : x;
251275
}
252-
// SIMD vectorized
276+
#ifdef __FDAPDE_HAS_EIGEN__ // SIMD vectorized
253277
using matrix_t = Eigen::Matrix<double, Dynamic, 1>;
254-
matrix_t variance(const matrix_t& x) const { return x; }
255-
matrix_t link (const matrix_t& x) const { return x.array().log(); }
256-
matrix_t inv_link(const matrix_t& x) const { return x.array().exp(); }
257-
matrix_t der_link(const matrix_t& x) const { return x.array().inverse(); }
258-
278+
matrix_t variance(const matrix_t& x) const override { return x; }
279+
matrix_t link (const matrix_t& x) const override { return x.array().log(); }
280+
matrix_t inv_link(const matrix_t& x) const override { return x.array().exp(); }
281+
matrix_t der_link(const matrix_t& x) const override { return x.array().inverse(); }
282+
double deviance(const matrix_t& x, const matrix_t& y) const override {
283+
fdapde_assert(x.cols() == 1 && y.cols() == 1 && x.rows() == y.rows());
284+
return ((y.array() > 0).select(y.array() * ((y.array() / x.array()).log() - 1) + x.array(), x.array()))
285+
.sum();
286+
}
287+
#endif
259288
template <typename T> constexpr auto transform(const T& data) const {
260289
if constexpr (internals::is_eigen_dense_xpr_v<T>) {
261290
return matrix_t((data.array() <= 0).select(1.0, data));
@@ -311,12 +340,16 @@ struct exponential_distribution : public internals::distribution_base<std::expon
311340
constexpr result_type deviance(T x, T y) const {
312341
return 2 * ((y - x) / x - std::log(y / x));
313342
}
314-
#ifdef __FDAPDE_HAS_EIGEN__ // SIMD vectorized
343+
#ifdef __FDAPDE_HAS_EIGEN__ // SIMD vectorized
315344
using matrix_t = Eigen::Matrix<double, Dynamic, 1>;
316-
matrix_t variance(const matrix_t& x) const { return x.array().pow(2); }
317-
matrix_t link (const matrix_t& x) const { return (-x).array().inverse(); }
318-
matrix_t inv_link(const matrix_t& x) const { return (-x).array().inverse(); }
319-
matrix_t der_link(const matrix_t& x) const { return x.array().pow(2).inverse(); }
345+
matrix_t variance(const matrix_t& x) const override { return x.array().pow(2); }
346+
matrix_t link (const matrix_t& x) const override { return (-x).array().inverse(); }
347+
matrix_t inv_link(const matrix_t& x) const override { return (-x).array().inverse(); }
348+
matrix_t der_link(const matrix_t& x) const override { return x.array().pow(2).inverse(); }
349+
double deviance(const matrix_t& x, const matrix_t& y) const override {
350+
fdapde_assert(x.cols() == 1 && y.cols() == 1 && x.rows() == y.rows());
351+
return (2 * ((y.array() - x.array()) / x.array() - (y.array() / x.array()).log())).sum();
352+
}
320353
#endif
321354
void set_param(param_type l) { l_ = l; }
322355
};
@@ -364,12 +397,16 @@ class gamma_distribution : public internals::distribution_base<std::gamma_distri
364397
constexpr result_type deviance(T x, T y) const {
365398
return 2 * ((y - x) / x - std::log(y / x));
366399
}
367-
#ifdef __FDAPDE_HAS_EIGEN__ // SIMD vectorized
400+
#ifdef __FDAPDE_HAS_EIGEN__ // SIMD vectorized
368401
using matrix_t = Eigen::Matrix<double, Dynamic, 1>;
369-
matrix_t variance(const matrix_t& x) const { return x.array().pow(2); }
370-
matrix_t link (const matrix_t& x) const { return (-x).array().inverse(); }
371-
matrix_t inv_link(const matrix_t& x) const { return (-x).array().inverse(); }
372-
matrix_t der_link(const matrix_t& x) const { return x.array().pow(2).inverse(); }
402+
matrix_t variance(const matrix_t& x) const override { return x.array().pow(2); }
403+
matrix_t link (const matrix_t& x) const override { return (-x).array().inverse(); }
404+
matrix_t inv_link(const matrix_t& x) const override { return (-x).array().inverse(); }
405+
matrix_t der_link(const matrix_t& x) const override { return x.array().pow(2).inverse(); }
406+
double deviance(const matrix_t& x, const matrix_t& y) const override {
407+
fdapde_assert(x.cols() == 1 && y.cols() == 1 && x.rows() == y.rows());
408+
return (2 * ((y.array() - x.array()) / x.array() - (y.array() / x.array()).log())).sum();
409+
}
373410
#endif
374411
void set_param(param_type k, param_type theta) {
375412
k_ = k;
@@ -420,12 +457,16 @@ class normal_distribution : public internals::distribution_base<std::normal_dist
420457
constexpr result_type deviance(T x, T y) const {
421458
return (x - y) * (x - y);
422459
}
423-
#ifdef __FDAPDE_HAS_EIGEN__ // SIMD vectorized
460+
#ifdef __FDAPDE_HAS_EIGEN__ // SIMD vectorized
424461
using matrix_t = Eigen::Matrix<double, Dynamic, 1>;
425-
matrix_t variance(const matrix_t& x) const { return matrix_t::Ones(x.rows()); }
426-
const matrix_t& link(const matrix_t& x) const { return x; }
427-
const matrix_t& inv_link(const matrix_t& x) const { return x; }
428-
matrix_t der_link(const matrix_t& x) const { return matrix_t::Ones(x.rows()); }
462+
matrix_t variance(const matrix_t& x) const override { return matrix_t::Ones(x.rows()); }
463+
matrix_t link (const matrix_t& x) const override { return x; }
464+
matrix_t inv_link(const matrix_t& x) const override { return x; }
465+
matrix_t der_link(const matrix_t& x) const override { return matrix_t::Ones(x.rows()); }
466+
double deviance(const matrix_t& x, const matrix_t& y) const override {
467+
fdapde_assert(x.cols() == 1 && y.cols() == 1 && x.rows() == y.rows());
468+
return (x - y).squaredNorm();
469+
}
429470
#endif
430471
void set_param(param_type mu, param_type sigma) {
431472
mu_ = mu;
@@ -484,10 +525,11 @@ class chi_squared_distribution : public internals::distribution_base<std::chi_sq
484525
}
485526
#ifdef __FDAPDE_HAS_EIGEN__ // SIMD vectorized
486527
using matrix_t = Eigen::Matrix<double, Dynamic, 1>;
487-
matrix_t variance(const matrix_t& x) const { return gamma_.variance(x); }
488-
matrix_t link (const matrix_t& x) const { return gamma_.link(x); }
489-
matrix_t inv_link(const matrix_t& x) const { return gamma_.inv_link(x); }
490-
matrix_t der_link(const matrix_t& x) const { return gamma_.der_link(x); }
528+
matrix_t variance(const matrix_t& x) const override { return gamma_.variance(x); }
529+
matrix_t link (const matrix_t& x) const override { return gamma_.link(x); }
530+
matrix_t inv_link(const matrix_t& x) const override { return gamma_.inv_link(x); }
531+
matrix_t der_link(const matrix_t& x) const override { return gamma_.der_link(x); }
532+
double deviance(const matrix_t& x, const matrix_t& y) const override { return gamma_.deviance(x, y); }
491533
#endif
492534
void set_param(param_type n) { n_ = n; }
493535
};

0 commit comments

Comments
 (0)