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  void BinaryMX<ScX, ScY>::eval_mx(const std::vector<MX>& arg, std::vector<MX>& res,
59  const std::vector<bool>& unique) const {
60  bool unique_x = !unique.empty() && unique[0];
61  bool unique_y = unique.size() >= 2 && unique[1];
62  res[0] = MX::binary(op_, arg[0], arg[1], unique_x, unique_y);
63  }
64 
65  template<bool ScX, bool ScY>
66  void BinaryMX<ScX, ScY>::eval_linear(const std::vector<std::array<MX, 3> >& arg,
67  std::vector<std::array<MX, 3> >& res) const {
68  casadi_math<MX>::fun_linear(op_, arg[0].data(), arg[1].data(), res[0].data());
69  }
70 
71  template<bool ScX, bool ScY>
72  void BinaryMX<ScX, ScY>::ad_forward(const std::vector<std::vector<MX> >& fseed,
73  std::vector<std::vector<MX> >& fsens) const {
74  // Get partial derivatives
75  MX pd[2];
76  casadi_math<MX>::der(op_, dep(0), dep(1), shared_from_this<MX>(), pd);
77 
78  // Propagate forward seeds
79  for (casadi_int d=0; d<fsens.size(); ++d) {
80  if (op_ == OP_IF_ELSE_ZERO) {
81  fsens[d][0] = if_else_zero(pd[1], fseed[d][1]);
82  } else {
83  fsens[d][0] = pd[0]*fseed[d][0] + pd[1]*fseed[d][1];
84  }
85  }
86  }
87 
88  template<bool ScX, bool ScY>
89  void BinaryMX<ScX, ScY>::ad_reverse(const std::vector<std::vector<MX> >& aseed,
90  std::vector<std::vector<MX> >& asens) const {
91  // Get partial derivatives
92  MX pd[2];
93  casadi_math<MX>::der(op_, dep(0), dep(1), shared_from_this<MX>(), pd);
94 
95  // Propagate adjoint seeds
96  for (casadi_int d=0; d<aseed.size(); ++d) {
97  MX s = aseed[d][0];
98  if (op_ == OP_IF_ELSE_ZERO) {
99  // Special case to avoid NaN propagation
100  if (!s.is_scalar() && dep(1).is_scalar()) {
101  asens[d][1] += dot(dep(0), s);
102  } else {
103  asens[d][1] += if_else_zero(dep(0), s);
104  }
105  } else {
106  // General case
107  for (casadi_int c=0; c<2; ++c) {
108  // Get increment of sensitivity c
109  MX t = pd[c]*s;
110 
111  // If dimension mismatch (i.e. one argument is scalar), then sum all the entries
112  if (!t.is_scalar() && t.size() != dep(c).size()) {
113  if (pd[c].size()!=s.size()) pd[c] = MX(s.sparsity(), pd[c]);
114  t = dot(pd[c], s);
115  }
116 
117  // Propagate the seeds
118  asens[d][c] += t;
119  }
120  }
121  }
122  }
123 
124  template<bool ScX, bool ScY>
125  void BinaryMX<ScX, ScY>::
126  generate(CodeGenerator& g,
127  const std::vector<casadi_int>& arg, const std::vector<casadi_int>& res,
128  const std::vector<bool>& arg_is_ref, std::vector<bool>& res_is_ref) const {
129  // Quick return if nothing to do
130  if (nnz()==0) return;
131 
132  // Check if inplace
133  bool inplace;
134  switch (op_) {
135  case OP_ADD:
136  case OP_SUB:
137  case OP_MUL:
138  case OP_DIV:
139  inplace = res[0]==arg[0] && !arg_is_ref[0];
140  break;
141  default:
142  inplace = false;
143  break;
144  }
145 
146  // Scalar names of arguments (start assuming all scalars)
147  std::string r = g.workel(res[0]);
148  std::string x = g.workel(arg[0]);
149  std::string y = g.workel(arg[1]);
150 
151  // Avoid emitting '/*' which will be mistaken for a comment
152  if (op_==OP_DIV && g.codegen_scalars && dep(1).nnz()==1) {
153  y = "(" + y + ")";
154  }
155 
156  // Codegen loop, if needed
157  if (nnz()>1) {
158  // Iterate over result
159  g.local("rr", "casadi_real", "*");
160  g.local("i", "casadi_int");
161  g << "for (i=0, " << "rr=" << g.work(res[0], nnz(), false);
162  r = "(*rr++)";
163 
164  // Iterate over first argument?
165  if (!ScX && !inplace) {
166  g.local("cr", "const casadi_real", "*");
167  g << ", cr=" << g.work(arg[0], dep(0).nnz(), arg_is_ref[0]);
168  if (op_==OP_OR || op_==OP_AND) {
169  // Avoid short-circuiting with side effects
170  x = "cr[i]";
171  } else {
172  x = "(*cr++)";
173  }
174 
175  }
176 
177  // Iterate over second argument?
178  if (!ScY) {
179  g.local("cs", "const casadi_real", "*");
180  g << ", cs=" << g.work(arg[1], dep(1).nnz(), arg_is_ref[1]);
181  if (op_==OP_OR || op_==OP_AND || op_==OP_IF_ELSE_ZERO) {
182  // Avoid short-circuiting with side effects
183  y = "cs[i]";
184  } else {
185  y = "(*cs++)";
186  }
187  }
188 
189  // Close loop
190  g << "; i<" << nnz() << "; ++i) ";
191  }
192 
193  // Perform operation
194  g << r << " ";
195  if (inplace) {
196  g << casadi_math<double>::sep(op_) << "= " << y;
197  } else {
198  g << " = " << g.print_op(op_, x, y);
199  }
200  g << ";\n";
201  }
202 
203  template<bool ScX, bool ScY>
204  int BinaryMX<ScX, ScY>::
205  eval(const double** arg, double** res, casadi_int* iw, double* w) const {
206  return eval_gen<double>(arg, res, iw, w);
207  }
208 
209  template<bool ScX, bool ScY>
210  int BinaryMX<ScX, ScY>::
211  eval_sx(const SXElem** arg, SXElem** res, casadi_int* iw, SXElem* w) const {
212  return eval_gen<SXElem>(arg, res, iw, w);
213  }
214 
215  template<bool ScX, bool ScY>
216  template<typename T>
217  int BinaryMX<ScX, ScY>::
218  eval_gen(const T* const* arg, T* const* res, casadi_int* iw, T* w) const {
219  // Get data
220  T* output0 = res[0];
221  const T* input0 = arg[0];
222  const T* input1 = arg[1];
223 
224  if (!ScX && !ScY) {
225  casadi_math<T>::fun(op_, input0, input1, output0, nnz());
226  } else if (ScX) {
227  casadi_math<T>::fun(op_, *input0, input1, output0, nnz());
228  } else {
229  casadi_math<T>::fun(op_, input0, *input1, output0, nnz());
230  }
231  return 0;
232  }
233 
234  template<bool ScX, bool ScY>
235  int BinaryMX<ScX, ScY>::
236  sp_forward(const bvec_t** arg, bvec_t** res, casadi_int* iw, bvec_t* w) const {
237  const bvec_t *a0=arg[0], *a1=arg[1];
238  bvec_t *r=res[0];
239  casadi_int n=nnz();
240  for (casadi_int i=0; i<n; ++i) {
241  if (ScX && ScY)
242  *r++ = *a0 | *a1;
243  else if (ScX && !ScY)
244  *r++ = *a0 | *a1++;
245  else if (!ScX && ScY)
246  *r++ = *a0++ | *a1;
247  else
248  *r++ = *a0++ | *a1++;
249  }
250  return 0;
251  }
252 
253  template<bool ScX, bool ScY>
254  int BinaryMX<ScX, ScY>::
255  eval_activity(const bvec_t** arg, bvec_t** res, casadi_int* iw, bvec_t* w) const {
256  // Zero-annihilation per operand, the value-level analog of sp_forward's OR
257  const bool f00 = operation_checker<F00Checker>(op_); // f(0,0)==0
258  const bool f0x = operation_checker<F0XChecker>(op_); // f(0,x)==0 (e.g. mul, div)
259  const bool fx0 = operation_checker<FX0Checker>(op_); // f(x,0)==0 (e.g. mul)
260  const bvec_t nz = ~static_cast<bvec_t>(0);
261  const bvec_t *a0=arg[0], *a1=arg[1];
262  bvec_t *r=res[0];
263  casadi_int n=nnz();
264  for (casadi_int i=0; i<n; ++i) {
265  const bool z0 = (*a0)!=0; // left active
266  const bool z1 = (*a1)!=0; // right active
267  bvec_t out;
268  if (!z0 && !z1) out = f00 ? 0 : nz;
269  else if (!z0 && z1) out = f0x ? 0 : nz;
270  else if ( z0 && !z1) out = fx0 ? 0 : nz;
271  else out = nz;
272  *r++ = out;
273  if (!ScX) a0++;
274  if (!ScY) a1++;
275  }
276  return 0;
277  }
278 
279  template<bool ScX, bool ScY>
280  int BinaryMX<ScX, ScY>::
281  sp_reverse(bvec_t** arg, bvec_t** res, casadi_int* iw, bvec_t* w) const {
282  bvec_t *a0=arg[0], *a1=arg[1], *r = res[0];
283  casadi_int n=nnz();
284  for (casadi_int i=0; i<n; ++i) {
285  bvec_t s = *r;
286  *r++ = 0;
287  if (ScX)
288  *a0 |= s;
289  else
290  *a0++ |= s;
291  if (ScY)
292  *a1 |= s;
293  else
294  *a1++ |= s;
295  }
296  return 0;
297  }
298 
299  template<bool ScX, bool ScY>
300  MX BinaryMX<ScX, ScY>::_get_binary(casadi_int op, const MX& y, bool scX, bool scY,
301  bool unique_x, bool unique_y) const {
302  if (!GlobalOptions::simplification_on_the_fly)
303  return MXNode::_get_binary(op, y, scX, scY, unique_x, unique_y);
304 
305  switch (op_) {
306  case OP_ADD:
307  if (op==OP_SUB && MX::is_equal(y, dep(0), maxDepth())) return dep(1);
308  if (op==OP_SUB && MX::is_equal(y, dep(1), maxDepth())) return dep(0);
309  break;
310  case OP_SUB:
311  if (op==OP_SUB && MX::is_equal(y, dep(0), maxDepth())) return -dep(1);
312  if (op==OP_ADD && MX::is_equal(y, dep(1), maxDepth())) return dep(0);
313  break;
314  default: break; // no rule
315  }
316 
317  // Fallback to default implementation
318  return MXNode::_get_binary(op, y, scX, scY, unique_x, unique_y);
319  }
320 
321  template<bool ScX, bool ScY>
322  void BinaryMX<ScX, ScY>::serialize_body(SerializingStream& s) const {
323  MXNode::serialize_body(s);
324  s.pack("BinaryMX::op", static_cast<int>(op_));
325  }
326 
327  template<bool ScX, bool ScY>
328  void BinaryMX<ScX, ScY>::serialize_type(SerializingStream& s) const {
329  MXNode::serialize_type(s);
330  char type_x = ScX;
331  char type_y = ScY;
332  char type = type_x | (type_y << 1);
333  s.pack("BinaryMX::scalar_flags", type);
334  }
335 
336  template<bool ScX, bool ScY>
337  MXNode* BinaryMX<ScX, ScY>::deserialize(DeserializingStream& s) {
338  char t;
339  s.unpack("BinaryMX::scalar_flags", t);
340  bool scX = t & 1;
341  bool scY = t & 2;
342 
343  if (scX) {
344  if (scY) return new BinaryMX<true, true>(s);
345  return new BinaryMX<true, false>(s);
346  } else {
347  if (scY) return new BinaryMX<false, true>(s);
348  return new BinaryMX<false, false>(s);
349  }
350  }
351 
352  template<bool ScX, bool ScY>
353  BinaryMX<ScX, ScY>::BinaryMX(DeserializingStream& s) : MXNode(s) {
354  int op;
355  s.unpack("BinaryMX::op", op);
356  op_ = Operation(op);
357  }
358 
359  template<bool ScX, bool ScY>
360  MX BinaryMX<ScX, ScY>::get_solve_triu(const MX& r, bool tr) const {
361  // Identify systems with the structure (I - R)
362  if (!ScX && !ScY && op_ == OP_SUB) {
363  // Is the first term a projected unity matrix?
364  if (dep(0).is_op(OP_PROJECT) && dep(0).dep(0).is_eye()) {
365  // Is the second term strictly lower triangular?
366  if (dep(1).is_op(OP_PROJECT) && dep(1).dep(0).sparsity().is_triu(true)) {
367  return dep(1).dep(0)->get_solve_triu_unity(r, tr);
368  }
369  }
370  }
371  // Fall back to default routine
372  return MXNode::get_solve_triu(r, tr);
373  }
374 
375  template<bool ScX, bool ScY>
376  MX BinaryMX<ScX, ScY>::get_solve_tril(const MX& r, bool tr) const {
377  // Identify systems with the structure (I - L)
378  if (!ScX && !ScY && op_ == OP_SUB) {
379  // Is the first term a projected unity matrix?
380  if (dep(0).is_op(OP_PROJECT) && dep(0).dep(0).is_eye()) {
381  // Is the second term strictly lower triangular?
382  if (dep(1).is_op(OP_PROJECT) && dep(1).dep(0).sparsity().is_tril(true)) {
383  return dep(1).dep(0)->get_solve_tril_unity(r, tr);
384  }
385  }
386  }
387  // Fall back to default routine
388  return MXNode::get_solve_tril(r, tr);
389  }
390 
391 } // namespace casadi
392 
393 #endif // CASADI_BINARY_MX_IMPL_HPP
casadi
The casadi namespace.
Definition: archiver.hpp:32