void BinaryMX<ScX, ScY>::
  generate(CodeGenerator& g, const std::string& mem,
           const std::vector<int>& arg, const std::vector<int>& res) const {
    // Quick return if nothing to do
    if (nnz()==0) return;

    // Check if inplace
    bool inplace;
    switch (op_) {
    case OP_ADD:
    case OP_SUB:
    case OP_MUL:
    case OP_DIV:
      inplace = res[0]==arg[0];
      break;
    default:
      inplace = false;
      break;
    }

    // Scalar names of arguments (start assuming all scalars)
    string r = g.workel(res[0]);
    string x = g.workel(arg[0]);
    string y = g.workel(arg[1]);

    // Codegen loop, if needed
    if (nnz()>1) {
      // Iterate over result
      g.local("rr", "real_t", "*");
      g.local("i", "int");
      g << "for (i=0, " << "rr=" << g.work(res[0], nnz());
      r = "(*rr++)";

      // Iterate over first argument?
      if (!ScX && !inplace) {
        g.local("cr", "const real_t", "*");
        g << ", cr=" << g.work(arg[0], dep(0).nnz());
        x = "(*cr++)";
      }

      // Iterate over second argument?
      if (!ScY) {
        g.local("cs", "const real_t", "*");
        g << ", cs=" << g.work(arg[1], dep(1).nnz());
        y = "(*cs++)";
      }

      // Close loop
      g << "; i<" << nnz() << "; ++i) ";
    }

    // Perform operation
    g << r << " ";
    if (inplace) {
      g << casadi_math<double>::sep(op_) << "= " << y;
    } else {
      g << " = " << casadi_math<double>::print(op_, x, y);
    }
    g << ";\n";
  }

 interval interval_relation_plugin::meet(interval const& src1, interval const& src2, bool& isempty) {
     isempty = false;
     if (is_empty(0, src1) || is_infinite(src2)) {
         return src1;
     }
     if (is_empty(0, src2) || is_infinite(src1)) {
         return src2;
     }
     bool l_open = src1.is_lower_open();
     bool r_open = src1.is_upper_open();
     ext_numeral low = src1.inf();
     ext_numeral high = src1.sup();
     if (src2.inf() > low || (src2.inf() == low && !l_open)) {
         low = src2.inf();
         l_open = src2.is_lower_open();
     }
     if (src2.sup() < high || (src2.sup() == high && !r_open)) {
         high = src2.sup();
         r_open = src2.is_upper_open();
     }
     if (low > high || (low == high && (l_open || r_open))) {
         isempty = true;
         return interval(dep());
     }
     else {
         return interval(dep(), low, l_open, nullptr, high, r_open, nullptr);
     }
 }

 void Split::generate(CodeGenerator& g, const std::string& mem,
                      const std::vector<int>& arg, const std::vector<int>& res) const {
   int nx = nout();
   for (int i=0; i<nx; ++i) {
     int nz_first = offset_[i];
     int nz_last = offset_[i+1];
     int nz = nz_last-nz_first;
     if (res[i]>=0 && nz>0) { // if anything to assign
       if (nz==1) { // assign scalar
         g.body << "  " << g.workel(res[i]) << " = ";
         if (dep(0).nnz()==1) {
           // rhs is also scalar
           casadi_assert(nz_first==0);
           g.body << g.workel(arg[0]) << ";" << endl;
         } else {
           // rhs is an element in a vector
           g.body << g.work(arg[0], dep(0).nnz()) << "[" << nz_first << "];" << endl;
         }
       } else {
         // assign vector
         std::string r = g.work(arg[0], dep(0).nnz());
         if (nz_first!=0) r = r + "+" + g.to_string(nz_first);
         g.body << "  " << g.copy(r, nz, g.work(res[i], nnz(i))) << endl;
       }
     }
   }
 }

 void InnerProd::evalFwd(const std::vector<std::vector<MX> >& fseed,
                         std::vector<std::vector<MX> >& fsens) {
   for (int d=0; d<fsens.size(); ++d) {
     fsens[d][0] = dep(0)->getInnerProd(fseed[d][1])
       + fseed[d][0]->getInnerProd(dep(1));
   }
 }

 void DenseTranspose::generate(const std::vector<int>& arg, const std::vector<int>& res,
                               CodeGenerator& g) const {
   g.body << "  for (i=0, rr=" << g.work(res[0], nnz()) << ", "
          << "cs=" << g.work(arg[0], nnz()) << "; i<" << dep().size2() << "; ++i) "
          << "for (j=0; j<" << dep().size1() << "; ++j) "
          << "rr[i+j*" << dep().size2() << "] = *cs++;" << endl;
 }

void goal::elim_redundancies() {
    if (inconsistent())
        return;
    expr_ref_fast_mark1 neg_lits(m());
    expr_ref_fast_mark2 pos_lits(m());
    unsigned sz = size();
    unsigned j  = 0;
    for (unsigned i = 0; i < sz; i++) {
        expr * f = form(i);
        if (m().is_true(f))
            continue;
        if (m().is_not(f)) {
            expr * atom = to_app(f)->get_arg(0);
            if (neg_lits.is_marked(atom))
                continue;
            if (pos_lits.is_marked(atom)) {
                proof * p = 0;
                if (proofs_enabled()) {
                    proof * prs[2] = { pr(get_idx(atom)), pr(i) };
                    p = m().mk_unit_resolution(2, prs);
                }
                expr_dependency_ref d(m());
                if (unsat_core_enabled())
                    d = m().mk_join(dep(get_idx(atom)), dep(i));
                push_back(m().mk_false(), p, d);                    
                return;
            }
            neg_lits.mark(atom);
        }
        else {
            if (pos_lits.is_marked(f))
                continue;
            if (neg_lits.is_marked(f)) {
                proof * p = 0;
                if (proofs_enabled()) {
                    proof * prs[2] = { pr(get_not_idx(f)), pr(i) };
                    p = m().mk_unit_resolution(2, prs);
                }
                expr_dependency_ref d(m());
                if (unsat_core_enabled())
                    d = m().mk_join(dep(get_not_idx(f)), dep(i));
                push_back(m().mk_false(), p, d);
                return;
            }
            pos_lits.mark(f);
        }
        if (i == j) {
            j++;
            continue;
        }
        m().set(m_forms, j, f);
        if (proofs_enabled())
            m().set(m_proofs, j, pr(i));
        if (unsat_core_enabled())
            m().set(m_dependencies, j, dep(i));
        j++;
    }
    shrink(j);
}

 void HorzRepmat::generate(CodeGenerator& g, const std::string& mem,
                           const std::vector<int>& arg, const std::vector<int>& res) const {
   int nnz = dep(0).nnz();
   g.body << "  for (i=0;i<" << n_ << ";++i) {" << endl;
   g.body << "    " << g.copy(g.work(arg[0], dep(0).nnz()), nnz,
                              g.work(res[0], sparsity().nnz()) + "+ i*" + g.to_string(nnz))
          << endl
   << "  }" << endl;
 }

  void BinaryMX<ScX, ScY>::evalFwd(const std::vector<std::vector<MX> >& fseed,
                                   std::vector<std::vector<MX> >& fsens) {
    // Get partial derivatives
    MX pd[2];
    casadi_math<MX>::der(op_, dep(0), dep(1), shared_from_this<MX>(), pd);

    // Propagate forward seeds
    for (int d=0; d<fsens.size(); ++d) {
      fsens[d][0] = pd[0]*fseed[d][0] + pd[1]*fseed[d][1];
    }
  }

 void Solve<Tr>::print(std::ostream &stream, long& remaining_calls) const{
   if(remaining_calls>0){
     remaining_calls--;
     stream << "(";
     dep(1)->print(stream,remaining_calls);
     if(Tr) stream << "'";
     stream << "\\";
     dep(0)->print(stream,remaining_calls);
     stream << ")";
   } else {
     stream << "...";
   }
 }

three_cells_t build23(reduce_t f, cell_t *x, cell_t *y) {
  three_cells_t r;
  r.a = closure_alloc(4);
  r.b = dep(r.a);
  r.c = dep(r.a);
  r.a->func = f;
  r.a->arg[0] = x;
  r.a->arg[1] = y;
  r.a->arg[2] = r.b;
  r.a->arg[3] = r.c;
  r.a->out = 2;
  return r;
}

  void DenseTranspose::evalGen(const T* const* arg, T* const* res,
                               int* iw, T* w) {
    // Get sparsity patterns
    int x_nrow = dep().size1();
    int x_ncol = dep().size2();

    const T* x = arg[0];
    T* xT = res[0];
    for (int i=0; i<x_ncol; ++i) {
      for (int j=0; j<x_nrow; ++j) {
        xT[i+j*x_ncol] = x[j+i*x_nrow];
      }
    }
  }

 /** \brief Check if two nodes are equivalent up to a given depth */
 virtual bool zz_isEqual(const MXNode* node, int depth) const {
   if (op_==node->getOp()) {
     if (isEqual(dep(0), node->dep(0), depth-1) && isEqual(dep(1), node->dep(1), depth-1)) {
       // If arguments are equal
       return true;
     } else {
       // If arguments are flipped
       return operation_checker<CommChecker>(op_) && isEqual(dep(1), node->dep(0), depth-1) &&
         isEqual(dep(0), node->dep(1), depth-1);
     }
   } else {
     return false;
   }
 }

// policy: callback whenever i~k but not j~k
void BinaryRelation::unsafe_merge (Ob i)
{
    Ob j = carrier().find(i);
    POMAGMA_ASSERT4(j < i, "BinaryRelation tried to merge item with self");

    DenseSet diff(item_dim());
    DenseSet rep(item_dim(), nullptr);
    DenseSet dep(item_dim(), nullptr);

    // merge rows (i, _) into (j, _)
    dep.init(m_lines.Lx(i));
    _remove_Rx(i, dep);
    rep.init(m_lines.Lx(j));
    if (rep.merge(dep, diff)) {
        for (auto k = diff.iter(); k.ok(); k.next()) {
            _insert_Rx(j, *k);
        }
    }

    // merge cols (_, i) into (_, j)
    dep.init(m_lines.Rx(i));
    _remove_Lx(dep, i);
    rep.init(m_lines.Rx(j));
    if (rep.merge(dep, diff)) {
        for (auto k = diff.iter(); k.ok(); k.next()) {
            _insert_Lx(*k, j);
        }
    }
}

bool CTrackball::DoTracking(int x, int y)
{
  if (!m_tracking)
    return false;

  Vector2df delta = m_start - Vector2df((float)x, (float)y);
  delta /= 2;

  Vector3df mouse(delta.x, -delta.y, 0.0f);
  Vector3df dep(0.0f, 0.0f, 1.0f);
  Vector3df axis(0.0f, 0.0f, 0.0f);

	// axis will be perpendicular to mouse and dep vectors
	// this is our rotation axis
  axis = mouse ^ dep;
	axis.Normalize();

	// every 100 pixels is 360 degress or rotation
  //GLfloat len = axis.Magnitude() / 100.0f;
  GLfloat len = fabs(delta.Magnitude());
  
	while (len > 360.0f)
		len -= 360.0f;

  ATLTRACE("delta_x : %f, delta_y: %f\n", mouse.x, mouse.y);
  ATLTRACE("axis: x:%f  y:%f  z:%f\n", axis.x, axis.y, axis.z, len);
  ATLTRACE("angle: %f\n", len);
  
  m_rotation_delta.CreateFromAxisAngle(axis.x, axis.y, axis.z, len);
  
  m_quat = m_rotation_delta * m_rotation;
  return true; // still tracking
}

 void Project::evalAdj(const std::vector<std::vector<MX> >& aseed,
                         std::vector<std::vector<MX> >& asens) {
   int nadj = aseed.size();
   for (int d=0; d<nadj; ++d) {
     asens[d][0] += project(aseed[d][0], dep().sparsity(), true);
   }
 }

  std::string SXNode::print_compact(std::map<const SXNode*, casadi_int>& nodeind,
                                   std::vector<std::string>& intermed) const {
    // Get reference to node index
    casadi_int& ind = nodeind[this];

    // If positive, already in intermediate expressions
    if (ind>0) {
      stringstream ss;
      ss << "@" << ind;
      return ss.str();
    }

    // Get expressions for dependencies
    std::string arg[2];
    for (casadi_int i=0; i<n_dep(); ++i) {
      arg[i] = dep(i)->print_compact(nodeind, intermed);
    }

    // Get expression for this
    string s = print(arg[0], arg[1]);

    // Decide what to do with the expression
    if (ind==0) {
      // Inline expression
      return s;
    } else {
      // Add to list of intermediate expressions and return reference
      intermed.push_back(s);
      ind = intermed.size(); // For subsequent references
      stringstream ss;
      ss << "@" << ind;
      return ss.str();
    }
  }

 void GetNonzerosSlice::simplifyMe(MX& ex) {
   // Simplify if identity
   if (isIdentity()) {
     MX t = dep(0);
     ex = t;
   }
 }

auto parse_format_string() {
    constexpr auto placeholder_head = here::find_placeholder(String+I, String+E)-String;
    return here::static_if([&] (auto dep, std::enable_if_t<dep(placeholder_head) == E>* = nullptr) {
        return text<String, I, dep(placeholder_head)> {};
    }, [&] (auto dep) {
        constexpr auto head = dep(placeholder_head);
        constexpr auto index_specifier_pos = here::find_index_specifier(String+head+1, String+E)-String;
        constexpr auto has_index = index_specifier_pos != head+1;
        constexpr auto index = has_index ? here::toi<unsigned>(String+head+1, String+index_specifier_pos-1)-1 : PlaceholderIndex;
        using placeholder = here::placeholder<String, index_specifier_pos, E, typename at_c<Seq, index>::type, index>;
        constexpr auto next_placeholder_index = PlaceholderIndex+!has_index+placeholder::used_indexes_count;
        return here::compose(text<String, I, head> {},
                             here::compose(placeholder{},
                                           here::parse_format_string<String, placeholder::end_pos, E, Seq, next_placeholder_index>()));
    });
}

  void Multiplication<TrX,TrY>::generateOperation(std::ostream &stream, const std::vector<std::string>& arg, const std::vector<std::string>& res, CodeGenerator& gen) const{
    // Check if inplace
    bool inplace = arg.at(0).compare(res.front())==0;

    // Copy first argument if not inplace
    if(!inplace){      
      stream << "  for(i=0; i<" << this->size() << "; ++i) " << res.front() << "[i]=" << arg.at(0) << "[i];" << endl;
    }

    // Perform sparse matrix multiplication
    gen.addAuxiliary(CodeGenerator::AUX_MM_TN_SPARSE);
    stream << "  casadi_mm_tn_sparse(";
    stream << arg.at(1) << ",s" << gen.getSparsity(dep(1).sparsity()) << ",";
    stream << arg.at(2) << ",s" << gen.getSparsity(dep(2).sparsity()) << ",";
    stream << res.front() << ",s" << gen.getSparsity(sparsity()) << ");" << endl;
  }

  void Reshape::evaluateMX(const MXPtrV& input, MXPtrV& output, const MXPtrVV& fwdSeed,
                           MXPtrVV& fwdSens, const MXPtrVV& adjSeed, MXPtrVV& adjSens,
                           bool output_given) {
    // Quick return if inplace
    if (input[0]==output[0]) return;

    if (!output_given) {
      *output[0] = reshape(*input[0], shape());
    }

    // Forward sensitivities
    int nfwd = fwdSens.size();
    for (int d = 0; d<nfwd; ++d) {
      *fwdSens[d][0] = reshape(*fwdSeed[d][0], shape());
    }

    // Adjoint sensitivities
    int nadj = adjSeed.size();
    for (int d=0; d<nadj; ++d) {
      MX& aseed = *adjSeed[d][0];
      MX& asens = *adjSens[d][0];
      asens.addToSum(reshape(aseed, dep().shape()));
      aseed = MX();
    }
  }