IMP logo
IMP Reference Guide  2.20.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-2022 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 <cereal/access.hpp>
19 #include <algorithm>
20 #include <cmath>
21 #include <limits>
22 #include <numeric>
23 
24 IMPALGEBRA_BEGIN_NAMESPACE
25 
26 //! Base class for a unit simplex embedded in D-dimensional real space.
27 /** See UnitSimplexD for more information. */
28 template <int D>
30  public:
31  //! Get dimension D of embedded real space.
32  virtual int get_dimension() const = 0;
33 
34  //! Get center of mass on simplex in embedded coordinates.
35  /** This is equivalent to the ones vector divided by D. */
37  double ci = 1.0 / static_cast<double>(get_dimension());
38  return get_ones_vector_kd<D>(get_dimension(), ci);
39  }
40 
41  //! Get whether the point is on the unit simplex
42  /**
43  \param[in] p Point
44  \param[in] atol Absolute tolerance of error for each element. Total
45  tolerance is atol*d.
46  */
48  const VectorD<D> &p,
49  double atol = std::numeric_limits<double>::epsilon()) const {
50  int d = p.get_dimension();
51  if (d != get_dimension()) return false;
52  double nrm = 0;
53  for (int i = 0; i < d; ++i) {
54  if (p[i] < 0) return false;
55  nrm += p[i];
56  }
57  if (std::abs(nrm - 1) > d * atol) return false;
58  return true;
59  }
60 };
61 
62 //! Represent a unit simplex embedded in D-dimensional real space.
63 /** The unit simplex (also known as standard simplex or probability simplex)
64  is the manifold of vectors of probabilities that sum to 1:
65  \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$
66 
67  \geometry
68  */
69 template <int D>
70 class UnitSimplexD : public UnitSimplexBaseD<D> {
71  friend class cereal::access;
72 
73  template<class Archive> void serialize(Archive &ar) {}
74 
75  public:
76  UnitSimplexD() { IMP_USAGE_CHECK(D > 0, "Dimension must be positive."); }
77 
78  //! Get dimension D of embedded real space.
79  int get_dimension() const override { return D; }
80 
81  IMP_SHOWABLE_INLINE(UnitSimplexD<D>, { out << "UnitSimplex" << D << "D"; });
82 };
83 
84 //! Represent a unit simplex embedded in d-dimensional real space.
85 /** This version's dimension is not known at compile time and can operate with
86  variable-dimension vectors.
87  */
88 template <>
89 class UnitSimplexD<-1> : public UnitSimplexBaseD<-1> {
90  friend class cereal::access;
91 
92  template<class Archive> void serialize(Archive &ar) {
93  ar(d_);
94  }
95 
96  public:
97  //! Construct from dimension d of embedded real space.
98  UnitSimplexD(int d = 1) : d_(d) {
99  IMP_USAGE_CHECK(d > 0, "Dimension must be positive.");
100  }
101 
102  //! Get dimension of embedded real space.
103  int get_dimension() const override { return d_; }
104 
105  // FIXME: SWIG doesn't seem to use this.
107  { out << "UnitSimplexKD(" << d_ << ")"; });
108 
109  private:
110  int d_;
111 };
112 
113 typedef UnitSimplexD<1> UnitSimplex1D;
114 IMP_VALUES(UnitSimplex1D, UnitSimplex1Ds);
116 IMP_VALUES(UnitSimplex2D, UnitSimplex2Ds);
118 IMP_VALUES(UnitSimplex3D, UnitSimplex3Ds);
120 IMP_VALUES(UnitSimplex4D, UnitSimplex4Ds);
122 IMP_VALUES(UnitSimplex5D, UnitSimplex5Ds);
124 IMP_VALUES(UnitSimplex6D, UnitSimplex6Ds);
125 typedef UnitSimplexD<-1> UnitSimplexKD;
126 IMP_VALUES(UnitSimplexKD, UnitSimplexKDs);
127 
128 //! Return a list of the vertices (bases) of the unit simplex.
129 template <int D>
131  Vector<VectorD<D> > ret;
132  int d = s.get_dimension();
133  for (int i = 0; i < d; ++i)
134  ret.push_back(get_basis_vector_kd<D>(d, i));
135  return ret;
136 }
137 
138 inline algebra::Triangle3D get_triangle_3d(const UnitSimplex3D &s) {
139  Vector3Ds ps = get_vertices(s);
140  return Triangle3D(ps[0], ps[1], ps[2]);
141 }
142 
143 //! Convert point on simplex from embedded to increasing coordinates.
144 /** Increasing coordinates are defined as the vector whose elements contain the
145  cumulative sum of all embedded coordinates of lower indices.
146  */
147 template <int D>
149  const VectorD<D> &p) {
150  int d = s.get_dimension();
151  IMP_USAGE_CHECK(static_cast<int>(p.get_dimension()) == d,
152  "Dimension of point must match dimension of simplex.");
153  IMP_INTERNAL_CHECK(s.get_contains(p), "Input point is not on simplex");
155  std::partial_sum(p.begin(), p.end(), q.begin());
156  return q;
157 }
158 
159 //! Convert point on simplex from increasing to embedded coordinates.
160 /** \see get_increasing_from_embedded */
161 template <int D>
163  const VectorD<D> &p) {
164  int d = s.get_dimension();
165  IMP_USAGE_CHECK(static_cast<int>(p.get_dimension()) == d,
166  "Dimension of point must match dimension of simplex.");
168  std::adjacent_difference(p.begin(), p.end(), q.begin());
169  IMP_INTERNAL_CHECK(s.get_contains(q), "Output point is not on simplex");
170  return q;
171 }
172 
173 // Perform Euclidean projection of p onto the unit simplex s.
174 /**
175  Any negative weights are set to 0, and the unit l1-norm is enforced. The
176  algorithm used to project to the simplex is described in
177  arXiv:1309.1541. It finds a threshold below which all weights are set to
178  0 and above which all weights shifted by the threshold sum to 1.
179  */
180 template <int D>
182  int d = s.get_dimension();
183  IMP_USAGE_CHECK(d == static_cast<int>(p.get_dimension()),
184  "Dimension of point must match dimension of simplex.");
185 
186  if (s.get_contains(p)) return p;
187 
188  VectorD<D> u(p);
189  std::sort(u.begin(), u.end(), std::greater<double>());
190 
191  Floats u_cumsum(d);
192  std::partial_sum(u.begin(), u.end(), u_cumsum.begin());
193 
194  int rho = 1;
195  while (rho < d) {
196  if (u[rho] + (1 - u_cumsum[rho]) / (rho + 1) < 0) break;
197  rho += 1;
198  }
199  double lam = (1 - u_cumsum[rho - 1]) / rho;
200 
201  for (int i = 0; i < d; ++i) {
202  double ui = p[i] + lam;
203  u[i] = ui > 0 ? ui : 0.0;
204  }
205 
206  return u;
207 }
208 
209 IMPALGEBRA_END_NAMESPACE
210 
211 #endif /* IMPALGEBRA_UNIT_SIMPLEX_D_H */
int get_dimension() const override
Get dimension D of embedded real space.
Definition: UnitSimplexD.h:79
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:280
Vector< VectorD< 3 > > Vector3Ds
Definition: VectorD.h:427
VectorD< D > get_barycenter() const
Get center of mass on simplex in embedded coordinates.
Definition: UnitSimplexD.h:36
Represent a unit simplex embedded in D-dimensional real space.
Definition: UnitSimplexD.h:70
Represent a triangle in 3D.
Definition: Triangle3D.h:24
A more IMP-like version of the std::vector.
Definition: Vector.h:42
#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:47
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:148
A Cartesian vector in D-dimensions.
Definition: VectorD.h:56
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:162
int get_dimension() const override
Get dimension of embedded real space.
Definition: UnitSimplexD.h:103
Simple D vector class.
Represent a triangle in 3D.
UnitSimplexD(int d=1)
Construct from dimension d of embedded real space.
Definition: UnitSimplexD.h:98
Base class for a unit simplex embedded in D-dimensional real space.
Definition: UnitSimplexD.h:29
Represent a unit simplex embedded in d-dimensional real space.
Definition: UnitSimplexD.h:89
Vector< VectorD< D > > get_vertices(const UnitSimplexD< D > &s)
Return a list of the vertices (bases) of the unit simplex.
Definition: UnitSimplexD.h:130
VectorD< D > get_projected(const UnitSimplexD< D > &s, const VectorD< D > &p)
Definition: UnitSimplexD.h:181
#define IMP_USAGE_CHECK(expr, message)
A runtime test for incorrect usage of a class or method.
Definition: check_macros.h:168
Macros to help with objects that can be printed to a stream.