CppADCodeGen  HEAD
A C++ Algorithmic Differentiation Package with Source Code Generation
abstract_atomic_fun.hpp
1 #ifndef CPPAD_CG_ABSTRACT_ATOMIC_FUN_INCLUDED
2 #define CPPAD_CG_ABSTRACT_ATOMIC_FUN_INCLUDED
3 /* --------------------------------------------------------------------------
4  * CppADCodeGen: C++ Algorithmic Differentiation with Source Code Generation:
5  * Copyright (C) 2013 Ciengis
6  * Copyright (C) 2018 Joao Leal
7  *
8  * CppADCodeGen is distributed under multiple licenses:
9  *
10  * - Eclipse Public License Version 1.0 (EPL1), and
11  * - GNU General Public License Version 3 (GPL3).
12  *
13  * EPL1 terms and conditions can be found in the file "epl-v10.txt", while
14  * terms and conditions for the GPL3 can be found in the file "gpl3.txt".
15  * ----------------------------------------------------------------------------
16  * Author: Joao Leal
17  */
18 
19 namespace CppAD {
20 namespace cg {
21 
27 template <class Base>
28 class CGAbstractAtomicFun : public BaseAbstractAtomicFun<Base> {
29 public:
30  using Super = BaseAbstractAtomicFun<Base>;
31  using CGB = CppAD::cg::CG<Base>;
32  using Arg = Argument<Base>;
33 protected:
37  const size_t id_;
43  bool standAlone_;
44 
45 protected:
46 
57  explicit CGAbstractAtomicFun(const std::string& name,
58  bool standAlone = false) :
59  Super(name),
60  id_(createNewAtomicFunctionID()),
61  standAlone_(standAlone) {
62  CPPADCG_ASSERT_KNOWN(!name.empty(), "The atomic function name cannot be empty")
63  this->option(CppAD::atomic_base<CGB>::set_sparsity_enum);
64  }
65 
66 public:
67  virtual ~CGAbstractAtomicFun() = default;
68 
69  template <class ADVector>
70  void operator()(const ADVector& ax,
71  ADVector& ay,
72  size_t id = 0) {
73  this->BaseAbstractAtomicFun<Base>::operator()(ax, ay, id);
74  }
75 
81  inline size_t getId() const {
82  return id_;
83  }
84 
90  inline bool isStandAlone() const {
91  return standAlone_;
92  }
93 
94  bool forward(size_t q,
95  size_t p,
96  const CppAD::vector<bool>& vx,
97  CppAD::vector<bool>& vy,
98  const CppAD::vector<CGB>& tx,
99  CppAD::vector<CGB>& ty) override {
100  using CppAD::vector;
101 
102  CppAD::vector<CGB> x;
103 
104  bool valuesDefined = BaseAbstractAtomicFun<Base>::isValuesDefined(tx);
105  if (vx.size() > 0) {
106  size_t n = vx.size();
107  x.resize(n);
108  for (size_t j = 0; j < n; j++) {
109  x[j] = tx[j * (p + 1)];
110  }
111 
112  zeroOrderDependency(vx, vy, x);
113  }
114 
115  bool allParameters = BaseAbstractAtomicFun<Base>::isParameters(tx);
116  if (allParameters) {
117  vector<Base> tyb;
118  if (!evalForwardValues(q, p, tx, tyb, ty.size()))
119  return false;
120 
121  CPPADCG_ASSERT_UNKNOWN(tyb.size() == ty.size())
122  for (size_t i = 0; i < ty.size(); i++) {
123  ty[i] = tyb[i];
124  }
125  return true;
126  }
127 
128  size_t m = ty.size() / (p + 1);
129 
130  vector<bool> vyLocal;
131  if (p == 0) {
132  vyLocal = vy;
133  } else if (p >= 1) {
139  size_t n = tx.size() / (p + 1);
140 
141  vector<std::set<size_t> > r(n);
142  for (size_t j = 0; j < n; j++) {
143  if (!tx[j * (p + 1) + 1].isIdenticalZero())
144  r[j].insert(0);
145  }
146  vector<std::set<size_t> > s(m);
147 
148  if(x.size() == 0) {
149  x.resize(n);
150  for (size_t j = 0; j < n; j++) {
151  x[j] = tx[j * (p + 1)];
152  }
153  }
154 
155  bool good = this->for_sparse_jac(1, r, s, x);
156  if (!good)
157  return false;
158 
159  vyLocal.resize(ty.size());
160  for (size_t i = 0; i < vyLocal.size(); i++) {
161  vyLocal[i] = true;
162  }
163 
164  for (size_t i = 0; i < m; i++) {
165  vyLocal[i * (p + 1) + 1] = !s[i].empty();
166  }
167 
168  if (p == 1) {
169  bool allZero = true;
170  for (size_t i = 0; i < vyLocal.size(); i++) {
171  if (vyLocal[i]) {
172  allZero = false;
173  break;
174  }
175  }
176 
177  if (allZero) {
178  for (size_t i = 0; i < ty.size(); i++) {
179  ty[i] = Base(0.0);
180  }
181  return true;
182  }
183  }
184  }
185 
186  vector<Base> tyb;
187  if (valuesDefined) {
188  if (!evalForwardValues(q, p, tx, tyb, ty.size()))
189  return false;
190  }
191 
192  CodeHandler<Base>* handler = findHandler(tx);
193  CPPADCG_ASSERT_UNKNOWN(handler != nullptr)
194 
195  size_t p1 = p + 1;
196 
197  std::vector<OperationNode<Base>*> txArray(p1), tyArray(p1);
198  for (size_t k = 0; k < p1; k++) {
199  if (k == 0)
200  txArray[k] = BaseAbstractAtomicFun<Base>::makeArray(*handler, tx, p, k);
201  else
202  txArray[k] = BaseAbstractAtomicFun<Base>::makeSparseArray(*handler, tx, p, k);
203  tyArray[k] = BaseAbstractAtomicFun<Base>::makeZeroArray(*handler, m);
204  }
205 
206  std::vector<Argument<Base> > args(2 * p1);
207  for (size_t k = 0; k < p1; k++) {
208  args[0 * p1 + k] = *txArray[k];
209  args[1 * p1 + k] = *tyArray[k];
210  }
211 
212  OperationNode<Base>* atomicOp = handler->makeNode(CGOpCode::AtomicForward,{id_, q, p}, args);
213  handler->registerAtomicFunction(*this);
214 
215  for (size_t k = 0; k < p1; k++) {
216  for (size_t i = 0; i < m; i++) {
217  size_t pos = i * p1 + k;
218  if (vyLocal.size() == 0 || vyLocal[pos]) {
219  ty[pos] = CGB(*handler->makeNode(CGOpCode::ArrayElement, {i}, {*tyArray[k], *atomicOp}));
220  if (valuesDefined) {
221  ty[pos].setValue(tyb[pos]);
222  }
223  } else {
224  CPPADCG_ASSERT_KNOWN(tyb.size() == 0 || IdenticalZero(tyb[pos]), "Invalid value")
225  ty[pos] = 0; // not a variable (zero)
226  }
227  }
228  }
229 
230  return true;
231  }
232 
233  bool reverse(size_t p,
234  const CppAD::vector<CGB>& tx,
235  const CppAD::vector<CGB>& ty,
236  CppAD::vector<CGB>& px,
237  const CppAD::vector<CGB>& py) override {
238  using CppAD::vector;
239 
240  bool allParameters = BaseAbstractAtomicFun<Base>::isParameters(tx);
241  if (allParameters) {
242  allParameters = BaseAbstractAtomicFun<Base>::isParameters(ty);
243  if (allParameters) {
244  allParameters = BaseAbstractAtomicFun<Base>::isParameters(py);
245  }
246  }
247 
248  if (allParameters) {
249  vector<Base> pxb;
250 
251  if (!evalReverseValues(p, tx, ty, pxb, py))
252  return false;
253 
254  CPPADCG_ASSERT_UNKNOWN(pxb.size() == px.size())
255 
256  for (size_t i = 0; i < px.size(); i++) {
257  px[i] = pxb[i];
258  }
259  return true;
260  }
261 
266  vector<bool> vxLocal(px.size());
267  for (size_t j = 0; j < vxLocal.size(); j++) {
268  vxLocal[j] = true;
269  }
270 
271  size_t p1 = p + 1;
272  // k == 0
273  size_t m = ty.size() / p1;
274  size_t n = tx.size() / p1;
275 
276  vector<std::set<size_t> > rt(m);
277  for (size_t i = 0; i < m; i++) {
278  if (!py[i * p1].isIdenticalZero()) {
279  rt[i].insert(0);
280  }
281  }
282 
283  CppAD::vector<CGB> x(n);
284  for (size_t j = 0; j < n; j++) {
285  x[j] = tx[j * p1];
286  }
287 
288  vector<std::set<size_t> > st(n);
289  bool good = this->rev_sparse_jac(1, rt, st, x);
290  if (!good) {
291  return false;
292  }
293 
294  for (size_t j = 0; j < n; j++) {
295  vxLocal[j * p1 + p] = !st[j].empty();
296  }
297 
298  if (p >= 1) {
303  vector<bool> vx(n);
304  vector<bool> s(m);
305  vector<bool> t(n);
306  vector<std::set<size_t> > r(n);
307  vector<std::set<size_t> > u(m);
308  vector<std::set<size_t> > v(n);
309 
310  for (size_t j = 0; j < n; j++) {
311  vx[j] = !tx[j * p1].isParameter();
312  if (!tx[j * p1 + 1].isIdenticalZero()) {
313  r[j].insert(0);
314  }
315  }
316  for (size_t i = 0; i < m; i++) {
317  s[i] = !py[i * p1 + 1].isIdenticalZero();
318  }
319 
320  this->rev_sparse_hes(vx, s, t, 1, r, u, v, x);
321 
322  for (size_t j = 0; j < n; j++) {
323  vxLocal[j * p1 + p - 1] = !v[j].empty();
324  }
325  }
326 
327  bool allZero = true;
328  for (size_t j = 0; j < vxLocal.size(); j++) {
329  if (vxLocal[j]) {
330  allZero = false;
331  break;
332  }
333  }
334 
335  if (allZero) {
336  for (size_t j = 0; j < px.size(); j++) {
337  px[j] = Base(0.0);
338  }
339  return true;
340  }
341 
342  bool valuesDefined = BaseAbstractAtomicFun<Base>::isValuesDefined(tx);
343  if (valuesDefined) {
344  valuesDefined = BaseAbstractAtomicFun<Base>::isValuesDefined(ty);
345  if (valuesDefined) {
346  valuesDefined = BaseAbstractAtomicFun<Base>::isValuesDefined(py);
347  }
348  }
349 
350  vector<Base> pxb;
351  if (valuesDefined) {
352  if (!evalReverseValues(p, tx, ty, pxb, py))
353  return false;
354  }
355 
356  CodeHandler<Base>* handler = findHandler(tx);
357  if (handler == nullptr) {
358  handler = findHandler(ty);
359  if (handler == nullptr) {
360  handler = findHandler(py);
361  }
362  }
363  CPPADCG_ASSERT_UNKNOWN(handler != nullptr)
364 
365  std::vector<OperationNode<Base>*> txArray(p1), tyArray(p1), pxArray(p1), pyArray(p1);
366  for (size_t k = 0; k <= p; k++) {
367  if (k == 0)
368  txArray[k] = BaseAbstractAtomicFun<Base>::makeArray(*handler, tx, p, k);
369  else
370  txArray[k] = BaseAbstractAtomicFun<Base>::makeSparseArray(*handler, tx, p, k);
371 
372  if (standAlone_) {
373  tyArray[k] = BaseAbstractAtomicFun<Base>::makeEmptySparseArray(*handler, m);
374  } else {
375  tyArray[k] = BaseAbstractAtomicFun<Base>::makeSparseArray(*handler, ty, p, k);
376  }
377 
378  if (k == 0)
379  pxArray[k] = BaseAbstractAtomicFun<Base>::makeZeroArray(*handler, n);
380  else
381  pxArray[k] = BaseAbstractAtomicFun<Base>::makeEmptySparseArray(*handler, n);
382 
383  if (k == 0)
384  pyArray[k] = BaseAbstractAtomicFun<Base>::makeSparseArray(*handler, py, p, k);
385  else
386  pyArray[k] = BaseAbstractAtomicFun<Base>::makeArray(*handler, py, p, k);
387  }
388 
389  std::vector<Argument<Base> > args(4 * p1);
390  for (size_t k = 0; k <= p; k++) {
391  args[0 * p1 + k] = *txArray[k];
392  args[1 * p1 + k] = *tyArray[k];
393  args[2 * p1 + k] = *pxArray[k];
394  args[3 * p1 + k] = *pyArray[k];
395  }
396 
397  OperationNode<Base>* atomicOp = handler->makeNode(CGOpCode::AtomicReverse,{id_, p}, args);
398  handler->registerAtomicFunction(*this);
399 
400  for (size_t k = 0; k < p1; k++) {
401  for (size_t j = 0; j < n; j++) {
402  size_t pos = j * p1 + k;
403  if (vxLocal[pos]) {
404  px[pos] = CGB(*handler->makeNode(CGOpCode::ArrayElement, {j}, {*pxArray[k], *atomicOp}));
405  if (valuesDefined) {
406  px[pos].setValue(pxb[pos]);
407  }
408  } else {
409  // CPPADCG_ASSERT_KNOWN(pxb.size() == 0 || IdenticalZero(pxb[j]), "Invalid value");
410  // pxb[j] might be non-zero but it is not required (it might have been used to determine other pxbs)
411  px[pos] = Base(0); // not a variable (zero)
412  }
413  }
414  }
415 
416  return true;
417  }
418 
419  inline virtual CppAD::vector<std::set<size_t>> jacobianForwardSparsitySet(size_t m,
420  const CppAD::vector<CGB>& x) {
421  size_t n = x.size();
422 
423  CppAD::vector<std::set<size_t>> r(n); // identity matrix
424  for (size_t i = 0; i < n; i++)
425  r[i].insert(i);
426 
427  CppAD::vector<std::set<size_t> > s(m);
428  bool good = this->for_sparse_jac(n, r, s, x);
429  if (!good)
430  throw CGException("Failed to compute jacobian sparsity pattern for atomic function '", this->atomic_name(), "'");
431 
432  return s;
433  }
434 
435  inline virtual CppAD::vector<std::set<size_t>> jacobianReverseSparsitySet(size_t m,
436  const CppAD::vector<CGB>& x) {
437  size_t n = x.size();
438 
439  CppAD::vector<std::set<size_t> > rt(m); // identity matrix
440  for (size_t i = 0; i < m; i++)
441  rt[i].insert(i);
442 
443  CppAD::vector<std::set<size_t> > st(n);
444  bool good = this->rev_sparse_jac(m, rt, st, x);
445  if (!good)
446  throw CGException("Failed to compute jacobian sparsity pattern for atomic function '", this->atomic_name(), "'");
447 
448  CppAD::vector<std::set<size_t>> s = transposePattern(st, n, m);
449 
450  return s;
451  }
452 
453  inline virtual CppAD::vector<std::set<size_t>> hessianSparsitySet(size_t m,
454  const CppAD::vector<CGB>& x) {
455  CppAD::vector<bool> s(m);
456  for (size_t i = 0; i < m; ++i)
457  s[i] = true;
458 
459  return hessianSparsitySet(s, x);
460  }
461 
462  inline virtual CppAD::vector<std::set<size_t>> hessianSparsitySet(const CppAD::vector<bool>& s,
463  const CppAD::vector<CGB>& x) {
464  size_t n = x.size();
465  size_t m = s.size();
466 
470  CppAD::vector<std::set<size_t>> r(n); // identity matrix
471  for (size_t j = 0; j < n; j++)
472  r[j].insert(j);
473 
474  CppAD::vector<bool> vx(n); // which x's are variables
475  for (size_t i = 0; i < n; ++i)
476  vx[i] = true;
477 
478  CppAD::vector<bool> t(n);
479  for (size_t i = 0; i < n; ++i)
480  t[i] = false;
481 
482  const CppAD::vector<std::set<size_t> > u(m); // empty
483  CppAD::vector<std::set<size_t> > v(n);
484 
485  bool good = this->rev_sparse_hes(vx, s, t, n, r, u, v, x);
486  if (!good)
487  throw CGException("Failed to compute Hessian sparsity pattern for atomic function '", this->atomic_name(), "'");
488 
489  return v;
490  }
491 
495  static size_t createNewAtomicFunctionID() {
496  CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
497  static size_t count = 0;
498  count++;
499  return count;
500  }
501 
502 protected:
503 
504  virtual void zeroOrderDependency(const CppAD::vector<bool>& vx,
505  CppAD::vector<bool>& vy,
506  const CppAD::vector<CGB>& x) = 0;
507 
522  virtual bool atomicForward(size_t q,
523  size_t p,
524  const CppAD::vector<Base>& tx,
525  CppAD::vector<Base>& ty) = 0;
538  virtual bool atomicReverse(size_t p,
539  const CppAD::vector<Base>& tx,
540  const CppAD::vector<Base>& ty,
541  CppAD::vector<Base>& px,
542  const CppAD::vector<Base>& py) = 0;
543 
544 private:
545 
546  inline bool evalForwardValues(size_t q,
547  size_t p,
548  const CppAD::vector<CGB>& tx,
549  CppAD::vector<Base>& tyb,
550  size_t ty_size) {
551  CppAD::vector<Base> txb(tx.size());
552  tyb.resize(ty_size);
553 
554  for (size_t i = 0; i < txb.size(); i++) {
555  txb[i] = tx[i].getValue();
556  }
557 
558  return atomicForward(q, p, txb, tyb);
559  }
560 
561  inline bool evalReverseValues(size_t p,
562  const CppAD::vector<CGB>& tx,
563  const CppAD::vector<CGB>& ty,
564  CppAD::vector<Base>& pxb,
565  const CppAD::vector<CGB>& py) {
566  using CppAD::vector;
567 
568  vector<Base> txb(tx.size());
569  vector<Base> tyb(ty.size());
570  pxb.resize(tx.size());
571  vector<Base> pyb(py.size());
572 
573  for (size_t i = 0; i < txb.size(); i++) {
574  txb[i] = tx[i].getValue();
575  }
576  for (size_t i = 0; i < tyb.size(); i++) {
577  tyb[i] = ty[i].getValue();
578  }
579  for (size_t i = 0; i < pyb.size(); i++) {
580  pyb[i] = py[i].getValue();
581  }
582 
583  return atomicReverse(p, txb, tyb, pxb, pyb);
584  }
585 
586 };
587 
588 } // END cg namespace
589 } // END CppAD namespace
590 
591 #endif
CppAD::cg::CGAbstractAtomicFun::hessianSparsitySet
virtual CppAD::vector< std::set< size_t > > hessianSparsitySet(const CppAD::vector< bool > &s, const CppAD::vector< CGB > &x)
Definition: abstract_atomic_fun.hpp:462
CppAD::cg::CGAbstractAtomicFun::atomicReverse
virtual bool atomicReverse(size_t p, const CppAD::vector< Base > &tx, const CppAD::vector< Base > &ty, CppAD::vector< Base > &px, const CppAD::vector< Base > &py)=0
CppAD::cg::CGAbstractAtomicFun::atomicForward
virtual bool atomicForward(size_t q, size_t p, const CppAD::vector< Base > &tx, CppAD::vector< Base > &ty)=0
CppAD::cg::CGAbstractAtomicFun::forward
bool forward(size_t q, size_t p, const CppAD::vector< bool > &vx, CppAD::vector< bool > &vy, const CppAD::vector< CGB > &tx, CppAD::vector< CGB > &ty) override
Definition: abstract_atomic_fun.hpp:94
CppAD::IdenticalZero
bool IdenticalZero(const cg::CG< Base > &x)
Definition: identical.hpp:37
CppAD::cg::CGAbstractAtomicFun::reverse
bool reverse(size_t p, const CppAD::vector< CGB > &tx, const CppAD::vector< CGB > &ty, CppAD::vector< CGB > &px, const CppAD::vector< CGB > &py) override
Definition: abstract_atomic_fun.hpp:233
CppAD::cg::CGAbstractAtomicFun::CGAbstractAtomicFun
CGAbstractAtomicFun(const std::string &name, bool standAlone=false)
Definition: abstract_atomic_fun.hpp:57