IMP logo
IMP Reference Guide  2.13.0
The Integrative Modeling Platform
UnitSimplexD.h
Go to the documentation of this file.
1 /**
2  * \file IMP/algebra/UnitSimplexD.h
3  * \brief Simple unit simplex class.
4  *
5  * Copyright 2007-2020 IMP Inventors. All rights reserved.
6  *
7  */
8 
9 #ifndef IMPALGEBRA_UNIT_SIMPLEX_D_H
10 #define IMPALGEBRA_UNIT_SIMPLEX_D_H
11 
12 #include <IMP/algebra/algebra_config.h>
13 #include <IMP/algebra/Triangle3D.h>
15 #include <IMP/algebra/VectorD.h>
16 #include <IMP/showable_macros.h>
17 #include <boost/math/special_functions/gamma.hpp>
18 #include <algorithm>
19 #include <cmath>
20 #include <limits>
21 #include <numeric>
22 
23 IMPALGEBRA_BEGIN_NAMESPACE
24 
25 //! Base class for a unit simplex embedded in D-dimensional real space.
26 /** See UnitSimplexD for more information. */
27 template <int D>
29  public:
30  //! Get dimension D of embedded real space.
31  virtual int get_dimension() const = 0;
32 
33  //! Get center of mass on simplex in embedded coordinates.
34  /** This is equivalent to the ones vector divided by D. */
36  double ci = 1.0 / static_cast<double>(get_dimension());
37  return get_ones_vector_kd<D>(get_dimension(), ci);
38  }
39 
40  //! Get whether the point is on the unit simplex
41  /**
42  \param[in] p Point
43  \param[in] atol Absolute tolerance of error for each element. Total
44  tolerance is atol*d.
45  */
47  const VectorD<D> &p,
48  double atol = std::numeric_limits<double>::epsilon()) const {
49  int d = p.get_dimension();
50  if (d != get_dimension()) return false;
51  double nrm = 0;
52  for (int i = 0; i < d; ++i) {
53  if (p[i] < 0) return false;
54  nrm += p[i];
55  }
56  if (std::abs(nrm - 1) > d * atol) return false;
57  return true;
58  }
59 };
60 
61 //! Represent a unit simplex embedded in D-dimensional real space.
62 /** The unit simplex (also known as standard simplex or probability simplex)
63  is the manifold of vectors of probabilities that sum to 1:
64  \f$ \Delta^{D-1} = \{ (t_1,\dots,t_D) \in \mathbb{R}^D \mid \sum_{d=1}^D t_d = 1, t_d \ge 0 \} \f$
65 
66  \geometry
67  */
68 template <int D>
69 class UnitSimplexD : public UnitSimplexBaseD<D> {
70  public:
71  UnitSimplexD() { IMP_USAGE_CHECK(D > 0, "Dimension must be positive."); }
72 
73  //! Get dimension D of embedded real space.
74  int get_dimension() const IMP_OVERRIDE { return D; }
75 
76  IMP_SHOWABLE_INLINE(UnitSimplexD<D>, { out << "UnitSimplex" << D << "D"; });
77 };
78 
79 //! Represent a unit simplex embedded in d-dimensional real space.
80 /** This version's dimension is not known at compile time and can operate with
81  variable-dimension vectors.
82  */
83 template <>
84 class UnitSimplexD<-1> : public UnitSimplexBaseD<-1> {
85  public:
86  //! Construct from dimension d of embedded real space.
87  UnitSimplexD(int d = 1) : d_(d) {
88  IMP_USAGE_CHECK(d > 0, "Dimension must be positive.");
89  }
90 
91  //! Get dimension of embedded real space.
92  int get_dimension() const IMP_OVERRIDE { return d_; }
93 
94  // FIXME: SWIG doesn't seem to use this.
96  { out << "UnitSimplexKD(" << d_ << ")"; });
97 
98  private:
99  int d_;
100 };
101 
102 typedef UnitSimplexD<1> UnitSimplex1D;
103 IMP_VALUES(UnitSimplex1D, UnitSimplex1Ds);
105 IMP_VALUES(UnitSimplex2D, UnitSimplex2Ds);
107 IMP_VALUES(UnitSimplex3D, UnitSimplex3Ds);
109 IMP_VALUES(UnitSimplex4D, UnitSimplex4Ds);
111 IMP_VALUES(UnitSimplex5D, UnitSimplex5Ds);
113 IMP_VALUES(UnitSimplex6D, UnitSimplex6Ds);
114 typedef UnitSimplexD<-1> UnitSimplexKD;
115 IMP_VALUES(UnitSimplexKD, UnitSimplexKDs);
116 
117 //! Return a list of the vertices (bases) of the unit simplex.
118 template <int D>
120  Vector<VectorD<D> > ret;
121  int d = s.get_dimension();
122  for (int i = 0; i < d; ++i)
123  ret.push_back(get_basis_vector_kd<D>(d, i));
124  return ret;
125 }
126 
127 inline algebra::Triangle3D get_triangle_3d(const UnitSimplex3D &s) {
128  Vector3Ds ps = get_vertices(s);
129  return Triangle3D(ps[0], ps[1], ps[2]);
130 }
131 
132 //! Convert point on simplex from embedded to increasing coordinates.
133 /** Increasing coordinates are defined as the vector whose elements contain the
134  cumulative sum of all embedded coordinates of lower indices.
135  */
136 template <int D>
138  const VectorD<D> &p) {
139  int d = s.get_dimension();
140  IMP_USAGE_CHECK(static_cast<int>(p.get_dimension()) == d,
141  "Dimension of point must match dimension of simplex.");
142  IMP_INTERNAL_CHECK(s.get_contains(p), "Input point is not on simplex");
144  std::partial_sum(p.begin(), p.end(), q.begin());
145  return q;
146 }
147 
148 //! Convert point on simplex from increasing to embedded coordinates.
149 /** \see get_increasing_from_embedded */
150 template <int D>
152  const VectorD<D> &p) {
153  int d = s.get_dimension();
154  IMP_USAGE_CHECK(static_cast<int>(p.get_dimension()) == d,
155  "Dimension of point must match dimension of simplex.");
157  std::adjacent_difference(p.begin(), p.end(), q.begin());
158  IMP_INTERNAL_CHECK(s.get_contains(q), "Output point is not on simplex");
159  return q;
160 }
161 
162 // Perform Euclidean projection of p onto the unit simplex s.
163 /**
164  Any negative weights are set to 0, and the unit l1-norm is enforced. The
165  algorithm used to project to the simplex is described in
166  arXiv:1309.1541. It finds a threshold below which all weights are set to
167  0 and above which all weights shifted by the threshold sum to 1.
168  */
169 template <int D>
171  int d = s.get_dimension();
172  IMP_USAGE_CHECK(d == static_cast<int>(p.get_dimension()),
173  "Dimension of point must match dimension of simplex.");
174 
175  if (s.get_contains(p)) return p;
176 
177  VectorD<D> u(p);
178  std::sort(u.begin(), u.end(), std::greater<double>());
179 
180  Floats u_cumsum(d);
181  std::partial_sum(u.begin(), u.end(), u_cumsum.begin());
182 
183  int rho = 1;
184  while (rho < d) {
185  if (u[rho] + (1 - u_cumsum[rho]) / (rho + 1) < 0) break;
186  rho += 1;
187  }
188  double lam = (1 - u_cumsum[rho - 1]) / rho;
189 
190  for (int i = 0; i < d; ++i) {
191  double ui = p[i] + lam;
192  u[i] = ui > 0 ? ui : 0.0;
193  }
194 
195  return u;
196 }
197 
198 IMPALGEBRA_END_NAMESPACE
199 
200 #endif /* IMPALGEBRA_UNIT_SIMPLEX_D_H */
Base class for geometric types.
#define IMP_SHOWABLE_INLINE(Name, how_to_show)
Declare the methods needed by an object that can be printed.
VectorD< D > get_zero_vector_kd(int Di)
Return a dynamically sized vector of zeros.
Definition: VectorD.h:276
Vector< VectorD< 3 > > Vector3Ds
Definition: VectorD.h:423
VectorD< D > get_barycenter() const
Get center of mass on simplex in embedded coordinates.
Definition: UnitSimplexD.h:35
Represent a unit simplex embedded in D-dimensional real space.
Definition: UnitSimplexD.h:69
Represent a triangle in 3D.
Definition: Triangle3D.h:23
A more IMP-like version of the std::vector.
Definition: Vector.h:39
int get_dimension() const
Get dimension of embedded real space.
Definition: UnitSimplexD.h:92
#define IMP_INTERNAL_CHECK(expr, message)
An assertion to check for internal errors in IMP. An IMP::ErrorException will be thrown.
Definition: check_macros.h:139
#define IMP_VALUES(Name, PluralName)
Define the type for storing sets of values.
Definition: value_macros.h:23
Base class for geometric types.
bool get_contains(const VectorD< D > &p, double atol=std::numeric_limits< double >::epsilon()) const
Get whether the point is on the unit simplex.
Definition: UnitSimplexD.h:46
VectorD< D > get_increasing_from_embedded(const UnitSimplexD< D > &s, const VectorD< D > &p)
Convert point on simplex from embedded to increasing coordinates.
Definition: UnitSimplexD.h:137
A Cartesian vector in D-dimensions.
Definition: VectorD.h:52
VectorD< D > get_embedded_from_increasing(const UnitSimplexD< D > &s, const VectorD< D > &p)
Convert point on simplex from increasing to embedded coordinates.
Definition: UnitSimplexD.h:151
Simple D vector class.
Represent a triangle in 3D.
UnitSimplexD(int d=1)
Construct from dimension d of embedded real space.
Definition: UnitSimplexD.h:87
Base class for a unit simplex embedded in D-dimensional real space.
Definition: UnitSimplexD.h:28
Represent a unit simplex embedded in d-dimensional real space.
Definition: UnitSimplexD.h:84
Vector< VectorD< D > > get_vertices(const UnitSimplexD< D > &s)
Return a list of the vertices (bases) of the unit simplex.
Definition: UnitSimplexD.h:119
VectorD< D > get_projected(const UnitSimplexD< D > &s, const VectorD< D > &p)
Definition: UnitSimplexD.h:170
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
int get_dimension() const
Get dimension D of embedded real space.
Definition: UnitSimplexD.h:74
Various general useful macros for IMP.
#define IMP_OVERRIDE
Cause a compile error if this method does not override a parent method.