Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
kernel_sum.hpp
1 /*
2  This file is part of the Kernel Quantum Probability library (KQP).
3 
4  KQP is free software: you can redistribute it and/or modify
5  it under the terms of the GNU General Public License as published by
6  the Free Software Foundation, either version 3 of the License, or
7  (at your option) any later version.
8 
9  KQP is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU General Public License for more details.
13 
14  You should have received a copy of the GNU General Public License
15  along with KQP. If not, see <http://www.gnu.org/licenses/>.
16  */
17 #ifndef __KQP_KERNEL_SUM_H__
18 #define __KQP_KERNEL_SUM_H__
19 
20 #include <vector>
21 #include <boost/lexical_cast.hpp>
22 #include <kqp/feature_matrix.hpp>
23 #include <kqp/space_factory.hpp>
24 
25 namespace kqp {
26  template<typename Scalar> class KernelSumSpace;
27 
34  template <typename Scalar>
35  class KernelSumMatrix : public FeatureMatrixBase<Scalar> {
36  public:
37  KQP_SCALAR_TYPEDEFS(Scalar);
39 
41  for(size_t i = 0; i < space.m_spaces.size(); i++) {
42  m_matrices.push_back(space.m_spaces[i]->newMatrix());
43  }
44  }
45 
46  KernelSumMatrix(const Self &other) : m_matrices(other.m_matrices) {
47  }
48 
49  KernelSumMatrix(const std::vector<FMatrixPtr> &matrices) : m_matrices(matrices) {
50  }
51 
52 # ifndef SWIG
53  KernelSumMatrix(Self &&other) {
54  m_matrices.swap(other.m_matrices);
55  }
56  Self &operator=(Self &&other) {
57  m_matrices.swap(other.m_matrices);
58  return *this;
59  }
60 # endif
61 
62  const FMatrixCPtr operator[](size_t i) const {
63  return m_matrices[i];
64  }
65 
66  virtual Index size() const override {
67  return m_matrices.empty() ? 0 : m_matrices[0]->size();
68  }
69 
70  virtual void add(const FMatrixBase &f, const std::vector<bool> *which = NULL) override {
71  const Self &other = f.template as<Self>();
72  for(size_t i = 0; i < m_matrices.size(); i++)
73  m_matrices[i]->add(other[i], which);
74  }
75 
76  virtual FMatrixBasePtr subset(const std::vector<bool>::const_iterator &begin,
77  const std::vector<bool>::const_iterator &end) const override {
78  boost::shared_ptr<Self> matrix = boost::shared_ptr<Self>(new Self());
79 
80  for(size_t i = 0; i < m_matrices.size(); i++) {
81  matrix->m_matrices.push_back(FMatrixPtr(m_matrices[i]->subset(begin, end)));
82  }
83  return matrix;
84  }
85 
86 
88  virtual FMatrixBase &operator=(const FeatureMatrixBase<Scalar> &other) override {
89  this->m_matrices = other.template as<Self>().m_matrices;
90  return *this;
91  }
92 
94  virtual FMatrixBasePtr copy() const override {
95  return FMatrixBasePtr(new Self(*this));
96  }
97 
98 
99  private:
100  KernelSumMatrix() {}
101  std::vector<FMatrixPtr> m_matrices;
102  };
103 
111  template <typename Scalar>
112  class KernelSumSpace : public SpaceBase<Scalar> {
113  public:
115  KQP_SPACE_TYPEDEFS("sum", Scalar);
116 
118  friend class KernelSumMatrix<Scalar>;
119 
121  void addSpace(Real weight, const FSpacePtr &space) {
122  m_weights.push_back(weight);
123  m_sum += weight * weight;
124  m_spaces.push_back(space);
125  }
126 
127  KernelSumSpace() : m_sum(0) {}
128  KernelSumSpace(const KernelSumSpace &other) : m_spaces(other.m_spaces), m_weights(other.m_weights), m_sum(other.m_sum) {}
129  ~KernelSumSpace() {}
130 
131  static FSpacePtr create() { return FSpacePtr(new KernelSumSpace()); }
132 
133  virtual Index dimension() const override {
134  Index n = 1;
135  for(size_t i = 0; i < m_spaces.size(); i++) {
136  n *= m_spaces[i]->dimension();
137  if (n < 0) return -1;
138  }
139 
140  return n;
141  }
142 
143  virtual boost::shared_ptr< SpaceBase<Scalar> > copy() const {
144  return FSpacePtr(new KernelSumSpace());
145  }
146 
147  virtual bool _canLinearlyCombine() const override { return false; }
148 
149  virtual FMatrixBasePtr newMatrix() const override {
150  return FMatrixBasePtr(new TMatrix(*this));
151  }
152  virtual FMatrixBasePtr newMatrix(const FeatureMatrixBase<Scalar> &other) const override {
153  return FMatrixBasePtr(new TMatrix(other.template as<TMatrix>()));
154  }
155 
156 
157  virtual const ScalarMatrix &k(const FMatrixBase &mX) const {
158  // FIXME: Big memory leak here (the feature space should be notified when a feature matrix is deleted) !!!
159  // So we suppose we are not multithreaded... and compute it each time.
160  // ScalarMatrix &gram = m_gramCache[static_cast<const void*>(&mX)];
161 
162  static ScalarMatrix gram;
163 
164  gram = ScalarMatrix::Zero(mX.size(), mX.size());
165 
166  for(size_t i = 0; i < m_spaces.size(); i++) {
167  gram += getNormalizedWeight(i) * m_spaces[i]->k(mX.template as<TMatrix>()[i]);
168  }
169  return gram;
170  }
171 
172  virtual ScalarMatrix k(const FMatrixBase &mX1, const ScalarAltMatrix &mY1, const RealAltVector &mD1,
173  const FMatrixBase &mX2, const ScalarAltMatrix &mY2, const RealAltVector &mD2) const override {
174  ScalarMatrix k = ScalarMatrix::Zero(mD1.rows(), mD2.rows());
175 
176  for(size_t i = 0; i < m_spaces.size(); i++) {
177  k += getNormalizedWeight(i) * m_spaces[i]->k(mX1.template as<TMatrix>()[i], mY1, mD1, mX2.template as<TMatrix>()[i], mY2, mD2);
178  }
179 
180  return k;
181  }
182 
183 
184 
185  virtual void update(std::vector< KernelValues<Scalar> > &values, int kOffset = 0) const override {
186  auto &self = values[kOffset];
187  kOffset++;
188 
189  self._inner = self._innerX = self._innerY = 0;
190  for(size_t i = 0; i < m_spaces.size(); i++) {
191  m_spaces[i]->update(values, kOffset);
192  const auto &child = values[kOffset];
193  self._inner += getNormalizedWeight(i) * child._inner;
194  self._innerX += getNormalizedWeight(i) * child._innerX;
195  self._innerY += getNormalizedWeight(i) * child._innerY;
196  kOffset += m_spaces[i]->numberOfKernelValues();
197  }
198 
199  }
200 
201  virtual void updatePartials(Real alpha, std::vector<Real> &partials, int offset, const std::vector< KernelValues<Scalar> > &values, int kOffset, int mode) const override {
202  Scalar k = values[kOffset].inner(mode);
203  kOffset += 1;
204 
205  for(size_t i = 0; i < m_spaces.size(); i++) {
206  partials[offset] += alpha * 2. * m_weights[i] * (values[kOffset].inner(mode) - k) / m_sum;
207  m_spaces[i]->updatePartials(alpha * getNormalizedWeight(i), partials, offset + 1, values, kOffset, mode);
208 
209  offset += m_spaces[i]->numberOfParameters(false) + 1;
210  kOffset += m_spaces[i]->numberOfKernelValues();
211  }
212 
213  }
214 
215  virtual int numberOfKernelValues() const {
216  int n = 1;
217  for(size_t i = 0; i < m_spaces.size(); i++)
218  n += m_spaces[i]->numberOfKernelValues();
219  return n;
220  }
221 
222  virtual int numberOfParameters(bool onlyFreeParameters) const {
223  int n = m_spaces.size();
224  for(size_t i = 0; i < m_spaces.size(); i++)
225  n += m_spaces[i]->numberOfParameters(onlyFreeParameters);
226  return onlyFreeParameters ? n - 1 : n;
227  }
228 
229  virtual void getParameters(bool onlyFreeParameters, std::vector<Real> & parameters, int offset = 0) const override {
230  for(size_t i = 0; i < m_spaces.size(); i++) {
231  if (!onlyFreeParameters || i < m_spaces.size() - 1) {
232  parameters[offset] = m_weights[i];
233  offset += 1;
234  }
235  m_spaces[i]->getParameters(onlyFreeParameters, parameters, offset);
236  offset += m_spaces[i]->numberOfParameters(onlyFreeParameters);
237  }
238  }
239 
240  virtual void setParameters(bool onlyFreeParameters, const std::vector<Real> & parameters, int offset = 0) override {
241  m_sum = 0;
242  for(size_t i = 0; i < m_spaces.size(); i++) {
243  if (!onlyFreeParameters || i < m_spaces.size() - 1) {
244  m_weights[i] = parameters[offset];
245  offset += 1;
246  }
247  m_sum += m_weights[i] * m_weights[i];
248  m_spaces[i]->setParameters(onlyFreeParameters, parameters, offset);
249  offset += m_spaces[i]->numberOfParameters(onlyFreeParameters);
250  }
251 
252  // If in only-free parameter mode, set the last parameter to be 1 - m_sum
253  if (onlyFreeParameters) {
254  if (m_sum > 1) {
255  m_weights[m_spaces.size()-1] = 0;
256  } else {
257  m_weights[m_spaces.size()-1] = std::sqrt(1 - m_sum);
258  m_sum = 1;
259  }
260  }
261  }
262 
263  virtual void getBounds(bool onlyFreeParameters, std::vector<Real> &lower, std::vector<Real> &upper, int offset = 0) const override
264  {
265  for(size_t i = 0; i < m_spaces.size(); i++) {
266  if (!onlyFreeParameters || i < m_spaces.size() - 1) {
267  lower[offset] = 0;
268  upper[offset] = 1;
269  offset += 1;
270  }
271  m_spaces[i]->getBounds(onlyFreeParameters, lower, upper, offset);
272  offset += m_spaces[i]->numberOfParameters(onlyFreeParameters);
273  }
274  }
275 
276 
277  virtual int getNumberOfConstraints(bool onlyFreeParameters) const
278  {
279  int n = onlyFreeParameters ? 1 : 0;
280  for(size_t i = 0; i < m_spaces.size(); i++) {
281  n += m_spaces[i]->getNumberOfConstraints(onlyFreeParameters);
282  }
283  return n;
284  }
285 
286  virtual void getConstraints(bool onlyFreeParameters, std::vector<Real> &constraintValues, int offset = 0) const
287  {
288  if (onlyFreeParameters) {
289  // Computes sum(free weights) - 1: (greater than 0 == constraint not respected)
290  Real c = -1.;
291  for(size_t i = 0; i < m_spaces.size() - 1; i++) {
292  c += m_weights[i] * m_weights[i];
293  }
294  constraintValues[offset] = c;
295  offset += 1;
296  }
297  for(size_t i = 0; i < m_spaces.size(); i++) {
298  m_spaces[i]->getConstraints(onlyFreeParameters, constraintValues, offset);
299  offset += m_spaces[i]->getNumberOfConstraints(onlyFreeParameters);
300  }
301  }
302 
303  FSpacePtr space(size_t i) { return m_spaces[i]; }
304  FSpaceCPtr space(size_t i) const { return m_spaces[i]; }
305  Real weight(size_t i) const { return m_weights[i]; }
306  size_t size() const { return m_spaces.size(); }
307 
308  virtual void load(const picojson::object &json) {
309  m_sum = 0;
310  m_spaces.clear();
311  m_weights.clear();
312 
313  auto p = json.find("list");
314  if (p == json.end())
315  KQP_THROW_EXCEPTION(exception, "No 'list' field in kernel sum");
316  if (!p->second.is<picojson::array>())
317  KQP_THROW_EXCEPTION(exception, "'list' field in kernel sum is not an array");
318 
319  for(auto child: p->second.get<picojson::array>()) {
320  if (!child.is<picojson::object>())
321  KQP_THROW_EXCEPTION(exception, "'list' field in kernel sum is not an array of objets");
322 
323  auto list = child.get<picojson::object>();
324 
325  Real weight = getNumeric<Real>("", list, "weight", 1.);
326  auto spaceJson = list.find("space");
327  if (spaceJson == list.end())
328  KQP_THROW_EXCEPTION(exception, "objet in 'list' does not contain a space");
329 
330 
331  auto space = kqp::our_dynamic_cast< SpaceBase<Scalar> >(SpaceFactory::load(spaceJson->second.get<picojson::object>()));
332 
333  addSpace(weight, space);
334  }
335  }
336 
337  virtual picojson::object save() const override {
338  picojson::object json = SpaceBase<Scalar>::save();
339  picojson::array array;
340  for(size_t i = 0; i < m_spaces.size(); i++) {
341  picojson::object listSpace;
342  listSpace["weight"] = picojson::value(m_weights[i]);
343  listSpace["space"] = picojson::value(m_spaces[i]->save());
344  array.push_back(picojson::value(listSpace));
345  }
346  json["list"] = picojson::value(array);
347  return json;
348  }
349 
350  inline Real getNormalizedWeight(size_t i) const {
351  return m_weights[i] * m_weights[i] / m_sum;
352  }
353  private:
354  std::vector<FSpacePtr> m_spaces;
355  std::vector<Real> m_weights;
356  mutable Real m_sum;
357  };
358 
359 }
360 
361 
362 #endif