Kernel Quantum Probability Library
The KQP library aims at providing tools for working with quantums probabilities
unary_kernel.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_UNARY_KERNEL_SPACE_H__
18 #define __KQP_UNARY_KERNEL_SPACE_H__
19 
20 #include <boost/lexical_cast.hpp>
21 #include <boost/unordered_map.hpp>
22 #include <kqp/feature_matrix.hpp>
23 #include <kqp/space_factory.hpp>
24 
25 namespace kqp {
26 
27  template<typename Scalar>
29  Scalar m_value;
30  Scalar m_lower;
31  Scalar m_upper;
32 
33  public:
34  BoundedParameter(Scalar value, Scalar lower = -std::numeric_limits<Scalar>::infinity(), Scalar upper = std::numeric_limits<Scalar>::infinity())
35  : m_value(value), m_lower(lower), m_upper(upper) {
36 
37  }
38 
39  BoundedParameter(const picojson::object &o, const std::string &name, Scalar defaultValue, Scalar defaultLower, Scalar defaultUpper) {
40  auto p = o.find(name);
41  if (p != o.end() && p->second.is<picojson::object>()) {
42  m_value = getNumeric<Scalar>(name, p->second, "value", defaultValue);
43  m_lower = getNumeric<Scalar>(name, p->second, "lower", defaultLower);
44  m_upper = getNumeric<Scalar>(name, p->second, "upper", defaultUpper);
45  } else {
46  m_value = getNumeric<Scalar>("", o, "name", defaultValue);
47  }
48  }
49 
50  inline Scalar upper() const { return m_upper; }
51  inline Scalar lower() const { return m_lower; }
52  inline operator Scalar() const { return m_value; }
53 
54  inline BoundedParameter &operator=(Scalar value) {
55  m_value = value;
56  return *this;
57  }
58 
59  picojson::value json() const {
60  picojson::value::object object;
61  object["value"] = picojson::value(m_value);
62  object["lower"] = picojson::value(m_lower);
63  object["upper"] = picojson::value(m_upper);
64  return picojson::value(object);
65  }
66 
67  };
68 
69  template<typename Scalar>
70  class UnaryKernelSpace : public SpaceBase< Scalar > {
71  public:
72  KQP_SCALAR_TYPEDEFS(Scalar);
73 
74  UnaryKernelSpace(const FSpacePtr &base) : m_base(base) {}
75 
76  ~UnaryKernelSpace() {}
77 
78  virtual FMatrixBasePtr newMatrix() const override { return m_base->newMatrix(); }
79  // REMOVE
80  // virtual FMatrixBasePtr newMatrix(const FeatureMatrixBase<Scalar> &other) const override { return m_base->newMatrix(other); }
81 
82  virtual const ScalarMatrix &k(const FMatrixBase &mX) const override {
83  // FIXME: Big memory leak here (the feature space should be notified when a feature matrix is deleted) !!!
84  // So we suppose we are not multithreaded... and compute it each time.
85  // ScalarMatrix &gram = m_gramCache[static_cast<const void*>(&mX)];
86 
87  static ScalarMatrix gram;
88 
89  gram.resize(0,0);
90 
91  if (mX.size() == gram.rows()) return gram;
92 
93  // We lose space here, could be used otherwise???
94  Index current = gram.rows();
95  if (current < mX.size())
96  gram.conservativeResize(mX.size(), mX.size());
97 
98  Index tofill = mX.size() - current;
99 
100  fillGram(gram, tofill, mX);
101 
102  gram.bottomLeftCorner(tofill, current).noalias() = gram.topRightCorner(current, tofill).adjoint().eval();
103 
104  return gram;
105  }
106 
107  virtual void load(const picojson::object &json) {
108  auto p = json.find("base");
109  if (p == json.end())
110  KQP_THROW_EXCEPTION(exception, "A unary kernel element should contain a base kernel");
111 
112  m_base = kqp::our_dynamic_cast< SpaceBase<Scalar> >(SpaceFactory::load(p->second.get<picojson::object>()));
113  }
114 
115  virtual picojson::object save() const override {
116  auto json = SpaceBase<Scalar>::save();
117  if (m_base)
118  json["base"] = picojson::value(m_base->save());
119  return json;
120  }
121 
122 
123  FSpaceCPtr base() const {
124  return m_base;
125  }
126 
127  virtual int numberOfKernelValues() const override {
128  return 1 + m_base->numberOfKernelValues();
129  }
130 
131 
132 
133  protected:
139  virtual void fillGram(ScalarMatrix &gram, Index tofill, const FMatrixBase &mX) const = 0;
140 
142  FSpacePtr m_base;
143 
144  private:
146  mutable boost::unordered_map<const void *, ScalarMatrix> m_gramCache;
147 
148  };
149 
150 
151 
153  template<typename Scalar> class GaussianSpace : public UnaryKernelSpace<Scalar> {
154  public:
155  typedef GaussianSpace<Scalar> Self;
156  KQP_SPACE_TYPEDEFS("gaussian", Scalar);
157 #ifndef SWIG
160 #endif
161  GaussianSpace(Real sigma, const FSpacePtr &base) : UnaryKernelSpace<Scalar>(base), m_sigma(sigma) {}
162  GaussianSpace() : UnaryKernelSpace<Scalar>(FSpacePtr()), m_sigma(1) {}
163 
164  virtual FSpacePtr copy() const override { return FSpacePtr(new GaussianSpace<Scalar>(m_sigma, m_base)); }
165 
166  virtual Index dimension() const override { return -1; }
167  virtual bool _canLinearlyCombine() const override { return false; }
168 
169 
170  virtual ScalarMatrix k(const FMatrixBase &mX1, const ScalarAltMatrix &mY1, const RealAltVector &mD1,
171  const FMatrixBase &mX2, const ScalarAltMatrix &mY2, const RealAltVector &mD2) const override {
172  return mD1.asDiagonal() * mY1.adjoint()
173  * this->f(m_base->k(mX1, mX2), m_base->k(mX1).diagonal(), m_base->k(mX2).diagonal())
174  * mY2 * mD2.asDiagonal();
175  }
176 
177 
178  virtual void update(std::vector< KernelValues<Scalar> > &values, int kOffset = 0) const override {
179  auto &self = values[kOffset];
180  auto &child = values[kOffset+1];
181  m_base->update(values, kOffset+1);
182  Scalar sigma_2 = m_sigma * m_sigma;
183  Scalar re_inner = Eigen::internal::real(child._inner);
184  self._inner = Eigen::internal::exp((2. * re_inner - child._innerX - child._innerY) / sigma_2);
185  self._innerX = 1;
186  self._innerY = 1;
187  }
188 
189  virtual void updatePartials(Real alpha, std::vector<Real> &partials, int offset, const std::vector< KernelValues<Scalar> > &values, int kOffset, int mode) const {
190  auto &child = values[kOffset+1];
191  Scalar exp_v = 1;
192  Scalar sigma_2 = m_sigma * m_sigma;
193  if (mode == 0) {
194  Scalar re_inner = Eigen::internal::real(child.inner(0));
195  Scalar v = (2. * re_inner - child.innerX(0) - child.innerY(0)) / sigma_2;
196  exp_v = Eigen::internal::exp(v);
197  partials[offset] += alpha * -2. / m_sigma * v * Eigen::internal::exp(v);
198  }
199 
200 
201  Scalar beta = alpha * (exp_v / sigma_2);
202  m_base->updatePartials(2 * beta, partials, offset+1, values, kOffset + 1, 0);
203  m_base->updatePartials(- beta, partials, offset+1, values, kOffset + 1, -1);
204  m_base->updatePartials(- beta, partials, offset+1, values, kOffset + 1, 1);
205  }
206 
207  virtual int numberOfParameters(bool onlyFreeParameters) const {
208  return 1 + m_base->numberOfParameters(onlyFreeParameters);
209  }
210 
211  virtual void getParameters(bool onlyFreeParameters, std::vector<Real> & parameters, int offset = 0) const {
212  parameters[offset] = m_sigma;
213  m_base->getParameters(onlyFreeParameters, parameters, offset + 1);
214  }
215 
216  virtual void getBounds(bool onlyFreeParameters, std::vector<Real> &lower, std::vector<Real> &upper, int offset = 0) const
217  {
218  lower[offset] = m_sigma.lower();
219  upper[offset] = m_sigma.upper();
220  m_base->getBounds(onlyFreeParameters, lower, upper, offset + 1);
221  }
222 
223  virtual void setParameters(bool onlyFreeParameters, const std::vector<Real> & parameters, int offset = 0) {
224  m_sigma = parameters[offset];
225  if (m_sigma < 0) m_sigma = -m_sigma;
226  if (m_sigma < epsilon()) m_sigma = kqp::epsilon();
227  m_base->setParameters(onlyFreeParameters, parameters, offset + 1);
228  }
229 
230  virtual void load(const picojson::object &json) override {
232  m_sigma = BoundedParameter<Real>(json, "sigma", 1., epsilon(), std::numeric_limits<Real>::infinity());
233  }
234 
235  virtual picojson::object save() const override {
236  picojson::object json = UnaryKernelSpace<Scalar>::save();
237  json["sigma"] = m_sigma.json();
238  return json;
239  }
240  protected:
241  virtual void fillGram(ScalarMatrix &gram, Index tofill, const FMatrixBase &mX) const override {
242  // Compute the remaining inner products
243  const ScalarMatrix &baseGram = m_base->k(mX);
244  gram.rightCols(tofill) = this->f(baseGram.rightCols(tofill), baseGram.diagonal(), baseGram.diagonal().tail(tofill));
245  }
246 
247 
248  private:
249 
250 
251  template<typename Derived, typename DerivedRow, typename DerivedCol>
252  inline ScalarMatrix f(const Eigen::MatrixBase<Derived> &k,
253  const Eigen::MatrixBase<DerivedRow>& rowNorms,
254  const Eigen::MatrixBase<DerivedCol>& colNorms) const {
255  return (-(rowNorms.derived().rowwise().replicate(k.cols()) + colNorms.derived().adjoint().colwise().replicate(k.rows()) - 2 * k.derived().real()) / (m_sigma*m_sigma)).array().exp();
256  }
257 
258  BoundedParameter<Real> m_sigma;
259  };
260 
262  template<typename Scalar> class PolynomialSpace : public UnaryKernelSpace<Scalar> {
263  public:
264  typedef GaussianSpace<Scalar> Self;
265  KQP_SPACE_TYPEDEFS("polynomial", Scalar);
266 #ifndef SWIG
269 #endif
270  virtual FSpacePtr copy() const override { return FSpacePtr(new PolynomialSpace<Scalar>(m_bias, m_degree, m_base)); }
271 
272  PolynomialSpace(Real bias, int degree, const FSpacePtr &base) : UnaryKernelSpace<Scalar>(base), m_bias(bias, 0), m_degree(degree) {}
273  PolynomialSpace() : UnaryKernelSpace<Scalar>(FSpacePtr()), m_bias(0), m_degree(1) {}
274 
275  virtual Index dimension() const override { return -1; }
276  virtual bool _canLinearlyCombine() const override { return false; }
277 
278  virtual ScalarMatrix k(const FMatrixBase &mX1, const ScalarAltMatrix &mY1, const RealAltVector &mD1,
279  const FMatrixBase &mX2, const ScalarAltMatrix &mY2, const RealAltVector &mD2) const override {
280  return mD1.asDiagonal() * mY1.adjoint() * this->f(m_base->k(mX1,mX2)) * mY2 * mD2.asDiagonal();
281  }
282 
283  virtual void updatePartials(Real alpha, std::vector<Real> &partials, int offset,
284  const std::vector< KernelValues<Scalar> > &values, int kOffset, int mode) const override {
285  Scalar v = (Scalar)m_degree * Eigen::internal::pow(values[kOffset+1].inner(mode) + m_bias, (Scalar)m_degree - 1);
286  partials[offset] += alpha * v;
287 
288  m_base->updatePartials(alpha * v, partials, offset+1, values, kOffset + 1, mode);
289  }
290 
291  virtual void update(std::vector< KernelValues<Scalar> > &values, int kOffset = 0) const override {
292  m_base->update(values, kOffset+1);
293  auto &self = values[kOffset];
294  auto &child = values[kOffset+1];
295  self._inner = Eigen::internal::pow(child._inner + m_bias, (Scalar)m_degree);
296  self._innerX = Eigen::internal::pow(child._innerX + m_bias, (Scalar)m_degree);
297  self._innerY = Eigen::internal::pow(child._innerY + m_bias, (Scalar)m_degree);
298  }
299 
300  virtual int numberOfParameters(bool onlyFreeParameters) const {
301  return 1 + m_base->numberOfParameters(onlyFreeParameters);
302  }
303 
304  virtual void getParameters(bool onlyFreeParameters, std::vector<Real> & parameters, int offset = 0) const override {
305  parameters[offset] = m_bias;
306  m_base->getParameters(onlyFreeParameters, parameters, offset + 1);
307  }
308 
309  virtual void getBounds(bool onlyFreeParameters, std::vector<Real> &lower, std::vector<Real> &upper, int offset = 0) const override
310  {
311  lower[offset] = m_bias.lower();
312  upper[offset] = m_bias.upper();
313  m_base->getBounds(onlyFreeParameters, lower, upper, offset + 1);
314  }
315 
316  virtual void setParameters(bool onlyFreeParameters, const std::vector<Real> & parameters, int offset = 0) override {
317  m_bias = parameters[offset];
318  m_base->setParameters(onlyFreeParameters, parameters, offset + 1);
319  }
320 
321 
322  virtual void load(const picojson::object &json) {
324  m_bias = BoundedParameter<Real>(json, "bias", 1., 0., std::numeric_limits<Real>::infinity());
325  m_degree = getNumeric<int>("", json, "degree", 2);
326  }
327 
328  virtual picojson::object save() const override {
329  picojson::object json = UnaryKernelSpace<Scalar>::save();
330  json["degree"] = picojson::value((double)m_degree);
331  json["bias"] = m_bias.json();
332  return json;
333  }
334 
335  protected:
336  virtual void fillGram(ScalarMatrix &gram, Index tofill, const FMatrixBase &mX) const override {
337  gram.rightCols(tofill) = this->f(m_base->k(mX).rightCols(tofill));
338  }
339 
340 
341  private:
342  template<typename Derived>
343  inline ScalarMatrix f(const Eigen::MatrixBase<Derived> &k) const {
344  return (k.derived() + ScalarMatrix::Constant(k.rows(), k.cols(), m_bias)).array().pow(m_degree);
345  }
346 
347  BoundedParameter<Real> m_bias;
348  int m_degree;
349  };
350 
351 }
352 
353 #endif