binary_mx_impl.hpp
1 /*
2  * This file is part of CasADi.
3  *
4  * CasADi -- A symbolic framework for dynamic optimization.
5  * Copyright (C) 2010-2023 Joel Andersson, Joris Gillis, Moritz Diehl,
6  * KU Leuven. All rights reserved.
7  * Copyright (C) 2011-2014 Greg Horn
8  *
9  * CasADi is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public
11  * License as published by the Free Software Foundation; either
12  * version 3 of the License, or (at your option) any later version.
13  *
14  * CasADi is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with CasADi; if not, write to the Free Software
21  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22  *
23  */
24 
25 
26 #ifndef CASADI_BINARY_MX_IMPL_HPP
27 #define CASADI_BINARY_MX_IMPL_HPP
28 
29 #include "binary_mx.hpp"
30 #include "casadi_misc.hpp"
31 #include "global_options.hpp"
32 #include "serializing_stream.hpp"
33 #include <sstream>
34 #include <vector>
35 
36 namespace casadi {
37 
38  template<bool ScX, bool ScY>
39  BinaryMX<ScX, ScY>::BinaryMX(Operation op, const MX& x, const MX& y) : op_(op) {
40  set_dep(x, y);
41  if (ScX) {
42  set_sparsity(y.sparsity());
43  } else {
44  set_sparsity(x.sparsity());
45  }
46  }
47 
48  template<bool ScX, bool ScY>
49  BinaryMX<ScX, ScY>::~BinaryMX() {
50  }
51 
52  template<bool ScX, bool ScY>
53  std::string BinaryMX<ScX, ScY>::disp(const std::vector<std::string>& arg) const {
54  return casadi_math<double>::print(op_, arg.at(0), arg.at(1));
55  }
56 
57  template<bool ScX, bool ScY>
58  bool BinaryMX<ScX, ScY>::is_nonnegative() const {
59  return operation_checker<NonnegativeChecker>(op_);
60  }
61 
62  template<bool ScX, bool ScY>
63  void BinaryMX<ScX, ScY>::eval_mx(const std::vector<MX>& arg, std::vector<MX>& res,
64  const std::vector<bool>& unique) const {
65  bool unique_x = !unique.empty() && unique[0];
66  bool unique_y = unique.size() >= 2 && unique[1];
67  res[0] = MX::binary(op_, arg[0], arg[1], unique_x, unique_y);
68  }
69 
70  template<bool ScX, bool ScY>
71  void BinaryMX<ScX, ScY>::eval_linear(const std::vector<std::array<MX, 3> >& arg,
72  std::vector<std::array<MX, 3> >& res) const {
73  casadi_math<MX>::fun_linear(op_, arg[0].data(), arg[1].data(), res[0].data());
74  }
75 
76  template<bool ScX, bool ScY>
77  void BinaryMX<ScX, ScY>::ad_forward(const std::vector<std::vector<MX> >& fseed,
78  std::vector<std::vector<MX> >& fsens) const {
79  // Get partial derivatives
80  MX pd[2];
81  casadi_math<MX>::der(op_, dep(0), dep(1), shared_from_this<MX>(), pd);
82 
83  // Propagate forward seeds
84  for (casadi_int d=0; d<fsens.size(); ++d) {
85  if (op_ == OP_IF_ELSE_ZERO) {
86  fsens[d][0] = if_else_zero(pd[1], fseed[d][1]);
87  } else {
88  fsens[d][0] = pd[0]*fseed[d][0] + pd[1]*fseed[d][1];
89  }
90  }
91  }
92 
93  template<bool ScX, bool ScY>
94  void BinaryMX<ScX, ScY>::ad_reverse(const std::vector<std::vector<MX> >& aseed,
95  std::vector<std::vector<MX> >& asens) const {
96  // Get partial derivatives
97  MX pd[2];
98  casadi_math<MX>::der(op_, dep(0), dep(1), shared_from_this<MX>(), pd);
99 
100  // Propagate adjoint seeds
101  for (casadi_int d=0; d<aseed.size(); ++d) {
102  MX s = aseed[d][0];
103  if (op_ == OP_IF_ELSE_ZERO) {
104  // Special case to avoid NaN propagation
105  if (!s.is_scalar() && dep(1).is_scalar()) {
106  asens[d][1] += dot(dep(0), s);
107  } else {
108  asens[d][1] += if_else_zero(dep(0), s);
109  }
110  } else {
111  // General case
112  for (casadi_int c=0; c<2; ++c) {
113  // Get increment of sensitivity c
114  MX t = pd[c]*s;
115 
116  // If dimension mismatch (i.e. one argument is scalar), then sum all the entries
117  if (!t.is_scalar() && t.size() != dep(c).size()) {
118  if (pd[c].size()!=s.size()) pd[c] = MX(s.sparsity(), pd[c]);
119  t = dot(pd[c], s);
120  }
121 
122  // Propagate the seeds
123  asens[d][c] += t;
124  }
125  }
126  }
127  }
128 
129  template<bool ScX, bool ScY>
130  void BinaryMX<ScX, ScY>::
131  generate(CodeGenerator& g,
132  const std::vector<casadi_int>& arg, const std::vector<casadi_int>& res,
133  const std::vector<bool>& arg_is_ref, std::vector<bool>& res_is_ref) const {
134  // Quick return if nothing to do
135  if (nnz()==0) return;
136 
137  // Check if inplace
138  bool inplace;
139  switch (op_) {
140  case OP_ADD:
141  case OP_SUB:
142  case OP_MUL:
143  case OP_DIV:
144  inplace = res[0]==arg[0] && !arg_is_ref[0];
145  break;
146  default:
147  inplace = false;
148  break;
149  }
150 
151  // Scalar names of arguments (start assuming all scalars)
152  std::string r = g.workel(res[0]);
153  std::string x = g.workel(arg[0]);
154  std::string y = g.workel(arg[1]);
155 
156  // Avoid emitting '/*' which will be mistaken for a comment
157  if (op_==OP_DIV && g.codegen_scalars && dep(1).nnz()==1) {
158  y = "(" + y + ")";
159  }
160 
161  // Codegen loop, if needed
162  if (nnz()>1) {
163  // Iterate over result
164  g.local("rr", "casadi_real", "*");
165  g.local("i", "casadi_int");
166  g << "for (i=0, " << "rr=" << g.work(res[0], nnz(), false);
167  r = "(*rr++)";
168 
169  // Iterate over first argument?
170  if (!ScX && !inplace) {
171  g.local("cr", "const casadi_real", "*");
172  g << ", cr=" << g.work(arg[0], dep(0).nnz(), arg_is_ref[0]);
173  if (op_==OP_OR || op_==OP_AND) {
174  // Avoid short-circuiting with side effects
175  x = "cr[i]";
176  } else {
177  x = "(*cr++)";
178  }
179 
180  }
181 
182  // Iterate over second argument?
183  if (!ScY) {
184  g.local("cs", "const casadi_real", "*");
185  g << ", cs=" << g.work(arg[1], dep(1).nnz(), arg_is_ref[1]);
186  if (op_==OP_OR || op_==OP_AND || op_==OP_IF_ELSE_ZERO) {
187  // Avoid short-circuiting with side effects
188  y = "cs[i]";
189  } else {
190  y = "(*cs++)";
191  }
192  }
193 
194  // Close loop
195  g << "; i<" << nnz() << "; ++i) ";
196  }
197 
198  // Perform operation
199  g << r << " ";
200  if (inplace) {
201  g << casadi_math<double>::sep(op_) << "= " << y;
202  } else {
203  g << " = " << g.print_op(op_, x, y);
204  }
205  g << ";\n";
206  }
207 
208  template<bool ScX, bool ScY>
209  int BinaryMX<ScX, ScY>::
210  eval(const double** arg, double** res, casadi_int* iw, double* w) const {
211  return eval_gen<double>(arg, res, iw, w);
212  }
213 
214  template<bool ScX, bool ScY>
215  int BinaryMX<ScX, ScY>::
216  eval_sx(const SXElem** arg, SXElem** res, casadi_int* iw, SXElem* w) const {
217  return eval_gen<SXElem>(arg, res, iw, w);
218  }
219 
220  template<bool ScX, bool ScY>
221  template<typename T>
222  int BinaryMX<ScX, ScY>::
223  eval_gen(const T* const* arg, T* const* res, casadi_int* iw, T* w) const {
224  // Get data
225  T* output0 = res[0];
226  const T* input0 = arg[0];
227  const T* input1 = arg[1];
228 
229  if (!ScX && !ScY) {
230  casadi_math<T>::fun(op_, input0, input1, output0, nnz());
231  } else if (ScX) {
232  casadi_math<T>::fun(op_, *input0, input1, output0, nnz());
233  } else {
234  casadi_math<T>::fun(op_, input0, *input1, output0, nnz());
235  }
236  return 0;
237  }
238 
239  template<bool ScX, bool ScY>
240  int BinaryMX<ScX, ScY>::
241  sp_forward(const bvec_t** arg, bvec_t** res, casadi_int* iw, bvec_t* w) const {
242  const bvec_t *a0=arg[0], *a1=arg[1];
243  bvec_t *r=res[0];
244  casadi_int n=nnz();
245  for (casadi_int i=0; i<n; ++i) {
246  if (ScX && ScY)
247  *r++ = *a0 | *a1;
248  else if (ScX && !ScY)
249  *r++ = *a0 | *a1++;
250  else if (!ScX && ScY)
251  *r++ = *a0++ | *a1;
252  else
253  *r++ = *a0++ | *a1++;
254  }
255  return 0;
256  }
257 
258  template<bool ScX, bool ScY>
259  int BinaryMX<ScX, ScY>::
260  eval_activity(const bvec_t** arg, bvec_t** res, casadi_int* iw, bvec_t* w) const {
261  // Zero-annihilation per operand, the value-level analog of sp_forward's OR
262  const bool f00 = operation_checker<F00Checker>(op_); // f(0,0)==0
263  const bool f0x = operation_checker<F0XChecker>(op_); // f(0,x)==0 (e.g. mul, div)
264  const bool fx0 = operation_checker<FX0Checker>(op_); // f(x,0)==0 (e.g. mul)
265  const bvec_t nz = ~static_cast<bvec_t>(0);
266  const bvec_t *a0=arg[0], *a1=arg[1];
267  bvec_t *r=res[0];
268  casadi_int n=nnz();
269  for (casadi_int i=0; i<n; ++i) {
270  const bool z0 = (*a0)!=0; // left active
271  const bool z1 = (*a1)!=0; // right active
272  bvec_t out;
273  if (!z0 && !z1) out = f00 ? 0 : nz;
274  else if (!z0 && z1) out = f0x ? 0 : nz;
275  else if ( z0 && !z1) out = fx0 ? 0 : nz;
276  else out = nz;
277  *r++ = out;
278  if (!ScX) a0++;
279  if (!ScY) a1++;
280  }
281  return 0;
282  }
283 
284  template<bool ScX, bool ScY>
285  int BinaryMX<ScX, ScY>::
286  sp_reverse(bvec_t** arg, bvec_t** res, casadi_int* iw, bvec_t* w) const {
287  bvec_t *a0=arg[0], *a1=arg[1], *r = res[0];
288  casadi_int n=nnz();
289  for (casadi_int i=0; i<n; ++i) {
290  bvec_t s = *r;
291  *r++ = 0;
292  if (ScX)
293  *a0 |= s;
294  else
295  *a0++ |= s;
296  if (ScY)
297  *a1 |= s;
298  else
299  *a1++ |= s;
300  }
301  return 0;
302  }
303 
304  template<bool ScX, bool ScY>
305  MX BinaryMX<ScX, ScY>::_get_binary(casadi_int op, const MX& y, bool scX, bool scY,
306  bool unique_x, bool unique_y) const {
307  if (!GlobalOptions::simplification_on_the_fly)
308  return MXNode::_get_binary(op, y, scX, scY, unique_x, unique_y);
309 
310  switch (op_) {
311  case OP_ADD:
312  if (op==OP_SUB && MX::is_equal(y, dep(0), maxDepth())) return dep(1);
313  if (op==OP_SUB && MX::is_equal(y, dep(1), maxDepth())) return dep(0);
314  break;
315  case OP_SUB:
316  if (op==OP_SUB && MX::is_equal(y, dep(0), maxDepth())) return -dep(1);
317  if (op==OP_ADD && MX::is_equal(y, dep(1), maxDepth())) return dep(0);
318  break;
319  default: break; // no rule
320  }
321 
322  // Fallback to default implementation
323  return MXNode::_get_binary(op, y, scX, scY, unique_x, unique_y);
324  }
325 
326  template<bool ScX, bool ScY>
327  void BinaryMX<ScX, ScY>::serialize_body(SerializingStream& s) const {
328  MXNode::serialize_body(s);
329  s.pack("BinaryMX::op", static_cast<int>(op_));
330  }
331 
332  template<bool ScX, bool ScY>
333  void BinaryMX<ScX, ScY>::serialize_type(SerializingStream& s) const {
334  MXNode::serialize_type(s);
335  char type_x = ScX;
336  char type_y = ScY;
337  char type = type_x | (type_y << 1);
338  s.pack("BinaryMX::scalar_flags", type);
339  }
340 
341  template<bool ScX, bool ScY>
342  MXNode* BinaryMX<ScX, ScY>::deserialize(DeserializingStream& s) {
343  char t;
344  s.unpack("BinaryMX::scalar_flags", t);
345  bool scX = t & 1;
346  bool scY = t & 2;
347 
348  if (scX) {
349  if (scY) return new BinaryMX<true, true>(s);
350  return new BinaryMX<true, false>(s);
351  } else {
352  if (scY) return new BinaryMX<false, true>(s);
353  return new BinaryMX<false, false>(s);
354  }
355  }
356 
357  template<bool ScX, bool ScY>
358  BinaryMX<ScX, ScY>::BinaryMX(DeserializingStream& s) : MXNode(s) {
359  int op;
360  s.unpack("BinaryMX::op", op);
361  op_ = Operation(op);
362  }
363 
364  template<bool ScX, bool ScY>
365  MX BinaryMX<ScX, ScY>::get_solve_triu(const MX& r, bool tr) const {
366  // Identify systems with the structure (I - R)
367  if (!ScX && !ScY && op_ == OP_SUB) {
368  // Is the first term a projected unity matrix?
369  if (dep(0).is_op(OP_PROJECT) && dep(0).dep(0).is_eye()) {
370  // Is the second term strictly lower triangular?
371  if (dep(1).is_op(OP_PROJECT) && dep(1).dep(0).sparsity().is_triu(true)) {
372  return dep(1).dep(0)->get_solve_triu_unity(r, tr);
373  }
374  }
375  }
376  // Fall back to default routine
377  return MXNode::get_solve_triu(r, tr);
378  }
379 
380  template<bool ScX, bool ScY>
381  MX BinaryMX<ScX, ScY>::get_solve_tril(const MX& r, bool tr) const {
382  // Identify systems with the structure (I - L)
383  if (!ScX && !ScY && op_ == OP_SUB) {
384  // Is the first term a projected unity matrix?
385  if (dep(0).is_op(OP_PROJECT) && dep(0).dep(0).is_eye()) {
386  // Is the second term strictly lower triangular?
387  if (dep(1).is_op(OP_PROJECT) && dep(1).dep(0).sparsity().is_tril(true)) {
388  return dep(1).dep(0)->get_solve_tril_unity(r, tr);
389  }
390  }
391  }
392  // Fall back to default routine
393  return MXNode::get_solve_tril(r, tr);
394  }
395 
396 } // namespace casadi
397 
398 #endif // CASADI_BINARY_MX_IMPL_HPP
casadi
The casadi namespace.
Definition: archiver.hpp:32