global.hpp

#ifndef HAVE_GLOBAL_HPP
#define HAVE_GLOBAL_HPP
// Autogenerated - do not edit by hand !
#include <algorithm>
#include <cmath>
#include <ctime>
#include <iomanip>
#include <iostream>
#include <limits>
#include <set>
#include <sstream>
#include <valarray>
#include <vector>
#include "config.hpp"
#include "radix.hpp"

namespace TMBad {

typedef TMBAD_HASH_TYPE hash_t;
typedef TMBAD_INDEX_TYPE Index;
typedef TMBAD_SCALAR_TYPE Scalar;
typedef std::pair<Index, Index> IndexPair;
typedef TMBAD_INDEX_VECTOR IndexVector;

struct global;
global *get_glob();

template <class T>
std::ostream &operator<<(std::ostream &out, const std::vector<T> &v) {
  out << "{";
  size_t last = v.size() - 1;
  for (size_t i = 0; i < v.size(); ++i) {
    out << v[i];
    if (i != last) out << ", ";
  }
  out << "}";
  return out;
}

template <class T>
struct intervals {
  struct ep : std::pair<T, bool> {
    bool left() const { return !this->second; }
    ep(T x, bool type) : std::pair<T, bool>(x, type) {}
    operator T() { return this->first; }
  };
  std::set<ep> x;
  typedef typename std::set<ep>::iterator iterator;
  bool insert(T a, T b) {
    ep x1(a, false);
    ep x2(b, true);
    iterator it1 = x.upper_bound(x1);
    iterator it2 = x.lower_bound(x2);

    bool insert_x1 = (it1 == x.end()) || it1->left();
    bool insert_x2 = (it2 == x.end()) || it2->left();

    bool change = (it1 != it2) || insert_x1;

    if (it1 != it2) {
      x.erase(it1, it2);
    }

    if (insert_x1) x.insert(x1);
    if (insert_x2) x.insert(x2);
    return change;
  }
  template <class F>
  F &apply(F &f) const {
    for (iterator it = x.begin(); it != x.end();) {
      ep a = *it;
      ++it;
      ep b = *it;
      ++it;
      f(a, b);
    }
    return f;
  }
  struct print_interval {
    void operator()(T a, T b) { Rcout << "[ " << a << " , " << b << " ] "; }
  };
  void print() {
    print_interval f;
    this->apply(f);
    Rcout << "\n";
  }
};

struct Dependencies : std::vector<Index> {
  typedef std::vector<Index> Base;
  std::vector<std::pair<Index, Index> > I;
  Dependencies();
  void clear();
  void add_interval(Index a, Index b);
  void add_segment(Index start, Index size);

  void monotone_transform_inplace(const std::vector<Index> &x);

  template <class F>
  F &apply(F &f) {
    for (size_t i = 0; i < this->size(); i++) f((*this)[i]);
    for (size_t i = 0; i < I.size(); i++) {
      for (Index j = I[i].first; j <= I[i].second; j++) {
        f(j);
      }
    }
    return f;
  }

  template <class F>
  F &apply_if_not_visited(F &f, intervals<Index> &visited) {
    for (size_t i = 0; i < this->size(); i++) f((*this)[i]);
    for (size_t i = 0; i < I.size(); i++) {
      if (visited.insert(I[i].first, I[i].second)) {
        for (Index j = I[i].first; j <= I[i].second; j++) {
          f(j);
        }
      }
    }
    return f;
  }

  bool any(const std::vector<bool> &x) const;
};

enum ArrayAccess { x_read, y_read, y_write, dx_read, dx_write, dy_read };
template <class Args, ArrayAccess What>
struct Accessor {};
template <class Args>
struct Accessor<Args, x_read> {
  typename Args::value_type operator()(const Args &args, Index j) const {
    return args.x(j);
  }
};
template <class Args>
struct Accessor<Args, y_read> {
  typename Args::value_type operator()(const Args &args, Index j) const {
    return args.y(j);
  }
};
template <class Args>
struct Accessor<Args, y_write> {
  typename Args::value_type &operator()(Args &args, Index j) {
    return args.y(j);
  }
};
template <class Args>
struct Accessor<Args, dx_read> {
  typename Args::value_type operator()(const Args &args, Index j) const {
    return args.dx(j);
  }
};
template <class Args>
struct Accessor<Args, dx_write> {
  typename Args::value_type &operator()(Args &args, Index j) {
    return args.dx(j);
  }
};
template <class Args>
struct Accessor<Args, dy_read> {
  typename Args::value_type operator()(const Args &args, Index j) const {
    return args.dy(j);
  }
};

template <class T>
struct IndirectAccessor {
  const std::vector<T> &x;
  const std::vector<Index> &i;
  IndirectAccessor(const std::vector<T> &x, const std::vector<Index> &i)
      : x(x), i(i) {}
  T operator[](size_t j) const { return x[i[j]]; }
  size_t size() const { return i.size(); }
  operator std::vector<T>() const {
    std::vector<T> ans(i.size());
    for (size_t j = 0; j < ans.size(); j++) ans[j] = (*this)[j];
    return ans;
  }
};

template <class Args, ArrayAccess What>
struct segment_ref {
  typedef typename Args::value_type Type;
  Accessor<Args, What> element_access;
  Args args;
  Index from, n;
  segment_ref(const Args &args, Index from, Index n)
      : args(args), from(from), n(n) {}
  template <class Other>
  operator Other() {
    Other ans(n);
    for (size_t i = 0; i < n; i++) {
      ans[i] = element_access(args, from + i);
    }
    return ans;
  }
  Type operator[](Index i) const { return element_access(args, from + i); }
  size_t size() const { return n; }
  template <class Other>
  segment_ref &operator=(const Other &other) {
    for (size_t i = 0; i < n; i++) {
      element_access(args, from + i) = other[i];
    }
    return *this;
  }
  template <class Other>
  segment_ref &operator+=(const Other &other) {
    for (size_t i = 0; i < n; i++) {
      element_access(args, from + i) += other[i];
    }
    return *this;
  }
  template <class Other>
  segment_ref &operator-=(const Other &other) {
    for (size_t i = 0; i < n; i++) {
      element_access(args, from + i) -= other[i];
    }
    return *this;
  }
};

template <class dummy = void>
struct Args {
  const Index *inputs;
  IndexPair ptr;
  Index input(Index j) const { return inputs[ptr.first + j]; }
  Index output(Index j) const { return ptr.second + j; }
  Args(const IndexVector &inputs) : inputs(inputs.data()) {
    ptr.first = 0;
    ptr.second = 0;
  }
};
template <class Type>
struct ForwardArgs : Args<> {
  typedef std::vector<Type> TypeVector;
  typedef Type value_type;
  Type *values;
  global *glob_ptr;
  Type x(Index j) const { return values[input(j)]; }
  Type &y(Index j) { return values[output(j)]; }
  Type *x_ptr(Index j) { return &values[input(j)]; }
  Type *y_ptr(Index j) { return &values[output(j)]; }
  segment_ref<ForwardArgs, x_read> x_segment(Index from, Index size) {
    return segment_ref<ForwardArgs, x_read>(*this, from, size);
  }
  segment_ref<ForwardArgs, y_write> y_segment(Index from, Index size) {
    return segment_ref<ForwardArgs, y_write>(*this, from, size);
  }
  ForwardArgs(const IndexVector &inputs, TypeVector &values,
              global *glob_ptr = NULL)
      : Args<>(inputs), values(values.data()), glob_ptr(glob_ptr) {}
};
template <class Type>
struct ReverseArgs : Args<> {
  typedef std::vector<Type> TypeVector;
  typedef Type value_type;
  Type *values;
  Type *derivs;
  global *glob_ptr;
  Type x(Index j) const { return values[input(j)]; }
  Type y(Index j) const { return values[output(j)]; }
  Type &dx(Index j) { return derivs[input(j)]; }
  Type dy(Index j) const { return derivs[output(j)]; }
  Type *x_ptr(Index j) { return &values[input(j)]; }
  Type *y_ptr(Index j) { return &values[output(j)]; }
  Type *dx_ptr(Index j) { return &derivs[input(j)]; }
  Type *dy_ptr(Index j) { return &derivs[output(j)]; }
  segment_ref<ReverseArgs, x_read> x_segment(Index from, Index size) {
    return segment_ref<ReverseArgs, x_read>(*this, from, size);
  }
  segment_ref<ReverseArgs, y_read> y_segment(Index from, Index size) {
    return segment_ref<ReverseArgs, y_read>(*this, from, size);
  }
  segment_ref<ReverseArgs, dx_write> dx_segment(Index from, Index size) {
    return segment_ref<ReverseArgs, dx_write>(*this, from, size);
  }
  segment_ref<ReverseArgs, dy_read> dy_segment(Index from, Index size) {
    return segment_ref<ReverseArgs, dy_read>(*this, from, size);
  }
  ReverseArgs(const IndexVector &inputs, TypeVector &values, TypeVector &derivs,
              global *glob_ptr = NULL)
      : Args<>(inputs),
        values(values.data()),
        derivs(derivs.data()),
        glob_ptr(glob_ptr) {
    ptr.first = (Index)inputs.size();
    ptr.second = (Index)values.size();
  }
};

template <>
struct ForwardArgs<bool> : Args<> {
  typedef std::vector<bool> BoolVector;
  BoolVector &values;
  intervals<Index> &marked_intervals;
  bool x(Index j) { return values[input(j)]; }
  BoolVector::reference y(Index j) { return values[output(j)]; }
  ForwardArgs(const IndexVector &inputs, BoolVector &values,
              intervals<Index> &marked_intervals)
      : Args<>(inputs), values(values), marked_intervals(marked_intervals) {}
  template <class Operator>
  bool any_marked_input(const Operator &op) {
    if (Operator::implicit_dependencies) {
      Dependencies dep;
      op.dependencies(*this, dep);
      return dep.any(values);
    } else {
      Index ninput = op.input_size();
      for (Index j = 0; j < ninput; j++)
        if (x(j)) return true;
    }
    return false;
  }
  template <class Operator>
  void mark_all_output(const Operator &op) {
    if (Operator::updating && op.output_size() == 0) {
      Dependencies dep;
      op.dependencies_updating(*this, dep);

      for (size_t i = 0; i < dep.size(); i++) values[dep[i]] = true;

      for (size_t i = 0; i < dep.I.size(); i++) {
        Index a = dep.I[i].first;
        Index b = dep.I[i].second;
        bool insert = marked_intervals.insert(a, b);
        if (insert) {
          for (Index j = a; j <= b; j++) {
            values[j] = true;
          }
        }
      }
    } else {
      Index noutput = op.output_size();
      for (Index j = 0; j < noutput; j++) y(j) = true;
    }
  }
  template <class Operator>
  bool mark_dense(const Operator &op) {
    if (any_marked_input(op)) {
      mark_all_output(op);
      return true;
    }
    return false;
  }
};

template <>
struct ReverseArgs<bool> : Args<> {
  typedef std::vector<bool> BoolVector;
  BoolVector &values;
  intervals<Index> &marked_intervals;
  BoolVector::reference x(Index j) { return values[input(j)]; }
  bool y(Index j) { return values[output(j)]; }
  ReverseArgs(IndexVector &inputs, BoolVector &values,
              intervals<Index> &marked_intervals)
      : Args<>(inputs), values(values), marked_intervals(marked_intervals) {
    ptr.first = (Index)inputs.size();
    ptr.second = (Index)values.size();
  }
  template <class Operator>
  bool any_marked_output(const Operator &op) {
    if (Operator::elimination_protected) return true;
    if (Operator::updating && op.output_size() == 0) {
      Dependencies dep;
      op.dependencies_updating(*this, dep);
      return dep.any(values);
    } else {
      Index noutput = op.output_size();
      for (Index j = 0; j < noutput; j++)
        if (y(j)) return true;
    }
    return false;
  }
  template <class Operator>
  void mark_all_input(const Operator &op) {
    if (Operator::implicit_dependencies) {
      Dependencies dep;
      op.dependencies(*this, dep);

      for (size_t i = 0; i < dep.size(); i++) values[dep[i]] = true;

      for (size_t i = 0; i < dep.I.size(); i++) {
        Index a = dep.I[i].first;
        Index b = dep.I[i].second;
        bool insert = marked_intervals.insert(a, b);
        if (insert) {
          for (Index j = a; j <= b; j++) {
            values[j] = true;
          }
        }
      }
    } else {
      Index ninput = op.input_size();
      for (Index j = 0; j < ninput; j++) x(j) = true;
    }
  }
  template <class Operator>
  bool mark_dense(const Operator &op) {
    if (any_marked_output(op)) {
      mark_all_input(op);
      return true;
    }
    return false;
  }
};

std::string tostr(const Index &x);

std::string tostr(const Scalar &x);

struct Writer : std::string {
  static std::ostream *cout;
  Writer(std::string str);
  Writer(Scalar x);
  Writer();

  template <class V>
  std::string vinit(const V &x) {
    std::string y = "{";
    for (size_t i = 0; i < x.size(); i++)
      y = y + (i == 0 ? "" : ",") + tostr(x[i]);
    y = y + "}";
    return y;
  }

  std::string p(std::string x);
  Writer operator+(const Writer &other);
  Writer operator-(const Writer &other);
  Writer operator-();
  Writer operator*(const Writer &other);
  Writer operator/(const Writer &other);

  Writer operator*(const Scalar &other);
  Writer operator+(const Scalar &other);

  void operator=(const Writer &other);
  void operator+=(const Writer &other);
  void operator-=(const Writer &other);
  void operator*=(const Writer &other);
  void operator/=(const Writer &other);

  template <class T>
  friend Writer &operator<<(Writer &w, const T &v) {
    *cout << v;
    return w;
  }
  template <class T>
  friend Writer &operator<<(Writer &w, const std::valarray<T> &x) {
    *cout << w.vinit(x);
    return w;
  }
};

template <>
struct ForwardArgs<Writer> : ForwardArgs<Scalar> {
  typedef std::vector<Scalar> ScalarVector;
  typedef ForwardArgs<Scalar> Base;
  bool const_literals;
  bool indirect;
  void set_indirect() {
    indirect = true;
    ptr.first = 0;
    ptr.second = 0;
  }
  Writer xd(Index j) { return "v[" + tostr(input(j)) + "]"; }
  Writer yd(Index j) { return "v[" + tostr(output(j)) + "]"; }
  Writer xi(Index j) { return "v[i[" + tostr(Index(ptr.first + j)) + "]]"; }
  Writer yi(Index j) { return "v[o[" + tostr(Index(ptr.second + j)) + "]]"; }
  Writer x(Index j) { return (indirect ? xi(j) : xd(j)); }
  Writer y(Index j) { return (indirect ? yi(j) : yd(j)); }
  Writer y_const(Index j) {
    TMBAD_ASSERT2(!indirect, "Attempt to write constants within loop?");
    return tostr(Base::y(j));
  }
  ForwardArgs(IndexVector &inputs, ScalarVector &values)
      : ForwardArgs<Scalar>(inputs, values) {
    const_literals = false;
    indirect = false;
  }
};

template <>
struct ReverseArgs<Writer> : Args<> {
  typedef std::vector<Scalar> ScalarVector;
  bool const_literals;
  bool indirect;
  void set_indirect() {
    indirect = true;
    ptr.first = 0;
    ptr.second = 0;
  }
  Writer dxd(Index j) { return "d[" + tostr(input(j)) + "]"; }
  Writer dyd(Index j) { return "d[" + tostr(output(j)) + "]"; }
  Writer xd(Index j) { return "v[" + tostr(input(j)) + "]"; }
  Writer yd(Index j) { return "v[" + tostr(output(j)) + "]"; }
  Writer dxi(Index j) { return "d[i[" + tostr(Index(ptr.first + j)) + "]]"; }
  Writer dyi(Index j) { return "d[o[" + tostr(Index(ptr.second + j)) + "]]"; }
  Writer xi(Index j) { return "v[i[" + tostr(Index(ptr.first + j)) + "]]"; }
  Writer yi(Index j) { return "v[o[" + tostr(Index(ptr.second + j)) + "]]"; }
  Writer x(Index j) { return (indirect ? xi(j) : xd(j)); }
  Writer y(Index j) { return (indirect ? yi(j) : yd(j)); }
  Writer dx(Index j) { return (indirect ? dxi(j) : dxd(j)); }
  Writer dy(Index j) { return (indirect ? dyi(j) : dyd(j)); }

  ReverseArgs(IndexVector &inputs, ScalarVector &values) : Args<>(inputs) {
    const_literals = false;
    indirect = false;
    ptr.first = (Index)inputs.size();
    ptr.second = (Index)values.size();
  }
};

struct Position {
  Position(Index node, Index first, Index second);
  Position();
  Index node;
  IndexPair ptr;
  bool operator<(const Position &other) const;
};

template <class T>
void sort_inplace(std::vector<T> &x) {
  std::sort(x.begin(), x.end());
}

template <class T>
void sort_unique_inplace(std::vector<T> &x) {
  std::sort(x.begin(), x.end());
  typename std::vector<T>::iterator last = std::unique(x.begin(), x.end());
  x.erase(last, x.end());
}

struct graph {
  std::vector<Index> j;
  std::vector<Index> p;
  graph();
  size_t num_neighbors(Index node);
  Index *neighbors(Index node);
  bool empty();
  size_t num_nodes();
  void print();
  std::vector<bool> mark;
  std::vector<Index> inv2op;
  std::vector<Index> dep2op;
  std::vector<Index> rowcounts();
  std::vector<Index> colcounts();
  void bfs(const std::vector<Index> &start, std::vector<bool> &visited,
           std::vector<Index> &result);
  void search(std::vector<Index> &start, bool sort_input = true,
              bool sort_output = true);
  void search(std::vector<Index> &start, std::vector<bool> &visited,
              bool sort_input = true, bool sort_output = true);
  std::vector<Index> boundary(const std::vector<Index> &subgraph);
  graph(size_t num_nodes, const std::vector<IndexPair> &edges);
};

namespace {
template <class CompleteOperator, bool dynamic>
struct constructOperator {};
template <class CompleteOperator>
struct constructOperator<CompleteOperator, false> {
  CompleteOperator *operator()() {
    static CompleteOperator *pOp = new CompleteOperator();
    return pOp;
  }
};
template <class CompleteOperator>
struct constructOperator<CompleteOperator, true> {
  CompleteOperator *operator()() {
    CompleteOperator *pOp = new CompleteOperator();
    return pOp;
  }

  template <class T1>
  CompleteOperator *operator()(const T1 &x1) {
    CompleteOperator *pOp = new CompleteOperator(x1);
    return pOp;
  }

  template <class T1, class T2>
  CompleteOperator *operator()(const T1 &x1, const T2 &x2) {
    CompleteOperator *pOp = new CompleteOperator(x1, x2);
    return pOp;
  }

  template <class T1, class T2, class T3>
  CompleteOperator *operator()(const T1 &x1, const T2 &x2, const T3 &x3) {
    CompleteOperator *pOp = new CompleteOperator(x1, x2, x3);
    return pOp;
  }

  template <class T1, class T2, class T3, class T4>
  CompleteOperator *operator()(const T1 &x1, const T2 &x2, const T3 &x3,
                               const T4 &x4) {
    CompleteOperator *pOp = new CompleteOperator(x1, x2, x3, x4);
    return pOp;
  }
};
}  // namespace

struct op_info {
  typedef int IntRep;
  IntRep code;
  enum op_flag {
    dynamic,
    smart_pointer,
    is_linear,
    is_constant,
    independent_variable,
    dependent_variable,
    allow_remap,
    elimination_protected,
    updating,
    op_flag_count
  };
  template <class T>
  IntRep get_flags(T op) {
    return

        (op.dynamic * (1 << dynamic)) |
        (op.smart_pointer * (1 << smart_pointer)) |
        (op.is_linear * (1 << is_linear)) |
        (op.is_constant * (1 << is_constant)) |
        (op.independent_variable * (1 << independent_variable)) |
        (op.dependent_variable * (1 << dependent_variable)) |
        (op.allow_remap * (1 << allow_remap)) |
        (op.elimination_protected * (1 << elimination_protected)) |
        (op.updating * (1 << updating));
  }
  op_info();
  op_info(op_flag f);

  template <class T>
  op_info(T op) : code(get_flags(op)) {}
  bool test(op_flag f) const;
  op_info &operator|=(const op_info &other);
  op_info &operator&=(const op_info &other);
};

struct global {
  struct ad_plain;
  struct ad_aug;
  typedef TMBAD_REPLAY_TYPE Replay;
  struct ad_segment;
  struct print_config;
  struct OperatorPure {
    virtual void increment(IndexPair &ptr) = 0;
    virtual void decrement(IndexPair &ptr) = 0;
    virtual void forward(ForwardArgs<Scalar> &args) = 0;
    virtual void reverse(ReverseArgs<Scalar> &args) = 0;
    virtual void forward_incr(ForwardArgs<Scalar> &args) = 0;
    virtual void reverse_decr(ReverseArgs<Scalar> &args) = 0;
    virtual Index input_size() = 0;
    virtual Index output_size() = 0;
    virtual void forward(ForwardArgs<bool> &args) = 0;
    virtual void reverse(ReverseArgs<bool> &args) = 0;
    virtual void forward_incr(ForwardArgs<bool> &args) = 0;
    virtual void reverse_decr(ReverseArgs<bool> &args) = 0;
    virtual void forward_incr_mark_dense(ForwardArgs<bool> &args) = 0;
    virtual void dependencies(Args<> &args, Dependencies &dep) = 0;
    virtual void dependencies_updating(Args<> &args, Dependencies &dep) = 0;
    virtual void forward(ForwardArgs<Replay> &args) = 0;
    virtual void reverse(ReverseArgs<Replay> &args) = 0;
    virtual void forward_incr(ForwardArgs<Replay> &args) = 0;
    virtual void reverse_decr(ReverseArgs<Replay> &args) = 0;
    virtual void forward(ForwardArgs<Writer> &args) = 0;
    virtual void reverse(ReverseArgs<Writer> &args) = 0;
    virtual void forward_incr(ForwardArgs<Writer> &args) = 0;
    virtual void reverse_decr(ReverseArgs<Writer> &args) = 0;
    virtual const char *op_name() { return "NoName"; }
    virtual OperatorPure *self_fuse() = 0;
    virtual OperatorPure *other_fuse(OperatorPure *other) = 0;
    virtual OperatorPure *copy() = 0;
    virtual void deallocate() = 0;
    virtual op_info info() = 0;
    virtual void *operator_data() = 0;
    virtual void *identifier() = 0;
    virtual void print(print_config cfg) = 0;
    virtual void *incomplete() = 0;
    virtual ~OperatorPure() {}
  };

  struct operation_stack : std::vector<OperatorPure *> {
    typedef std::vector<OperatorPure *> Base;
    op_info any;
    operation_stack();
    operation_stack(const operation_stack &other);
    void push_back(OperatorPure *x);
    operation_stack &operator=(const operation_stack &other);
    ~operation_stack();
    void clear();
    void copy_from(const operation_stack &other);
  };

  operation_stack opstack;
  std::vector<Scalar> values;
  std::vector<Scalar> derivs;
  IndexVector inputs;
  std::vector<Index> inv_index;
  std::vector<Index> dep_index;

  mutable std::vector<IndexPair> subgraph_ptr;
  std::vector<Index> subgraph_seq;
  void (*forward_compiled)(Scalar *);
  void (*reverse_compiled)(Scalar *, Scalar *);

  global();
  void clear();

  void shrink_to_fit(double tol = .9);

  void clear_deriv(Position start = Position(0, 0, 0));

  Scalar &value_inv(Index i);
  Scalar &deriv_inv(Index i);
  Scalar &value_dep(Index i);
  Scalar &deriv_dep(Index i);
  Position begin();
  Position end();

  struct no_filter {
    CONSTEXPR bool operator[](size_t i) const;
  };
  template <class ForwardArgs, class NodeFilter>
  void forward_loop(ForwardArgs &args, size_t begin,
                    const NodeFilter &node_filter) const {
    for (size_t i = begin; i < opstack.size(); i++) {
      if (node_filter[i])
        opstack[i]->forward_incr(args);
      else
        opstack[i]->increment(args.ptr);
    }
  }
  template <class ForwardArgs>
  void forward_loop(ForwardArgs &args, size_t begin = 0) const {
    forward_loop(args, begin, no_filter());
  }
  template <class ReverseArgs, class NodeFilter>
  void reverse_loop(ReverseArgs &args, size_t begin,
                    const NodeFilter &node_filter) const {
    for (size_t i = opstack.size(); i > begin;) {
      i--;
      if (node_filter[i])
        opstack[i]->reverse_decr(args);
      else
        opstack[i]->decrement(args.ptr);
    }
  }
  template <class ReverseArgs>
  void reverse_loop(ReverseArgs &args, size_t begin = 0) const {
    reverse_loop(args, begin, no_filter());
  }
  template <class ForwardArgs>
  void forward_loop_subgraph(ForwardArgs &args) const {
    subgraph_cache_ptr();
    for (size_t j = 0; j < subgraph_seq.size(); j++) {
      Index i = subgraph_seq[j];
      args.ptr = subgraph_ptr[i];
      opstack[i]->forward(args);
    }
  }
  template <class ReverseArgs>
  void reverse_loop_subgraph(ReverseArgs &args) const {
    subgraph_cache_ptr();
    for (size_t j = subgraph_seq.size(); j > 0;) {
      j--;
      Index i = subgraph_seq[j];
      args.ptr = subgraph_ptr[i];
      opstack[i]->reverse(args);
    }
  }
  template <class Vector>
  void clear_array_subgraph(Vector &array,
                            typename Vector::value_type value =
                                typename Vector::value_type(0)) const {
    if (array.size() != values.size()) {
      array.resize(values.size());
      std::fill(array.begin(), array.end(), value);
      return;
    }
    subgraph_cache_ptr();
    for (size_t j = 0; j < subgraph_seq.size(); j++) {
      Index i = subgraph_seq[j];
      size_t noutput = opstack[i]->output_size();
      for (size_t k = 0; k < noutput; k++)
        array[subgraph_ptr[i].second + k] = value;
    }
  }

  void forward(Position start = Position(0, 0, 0));
  void reverse(Position start = Position(0, 0, 0));
  void forward_sub();
  void reverse_sub();

  void forward(std::vector<bool> &marks);
  void reverse(std::vector<bool> &marks);
  void forward_sub(std::vector<bool> &marks,
                   const std::vector<bool> &node_filter = std::vector<bool>());
  void reverse_sub(std::vector<bool> &marks,
                   const std::vector<bool> &node_filter = std::vector<bool>());
  void forward_dense(std::vector<bool> &marks);

  intervals<Index> updating_intervals() const;

  intervals<Index> updating_intervals_sub() const;

  struct replay {
    std::vector<Replay> values;
    std::vector<Replay> derivs;
    const global &orig;
    global &target;
    global *parent_glob;
    Replay &value_inv(Index i);
    Replay &deriv_inv(Index i);
    Replay &value_dep(Index i);
    Replay &deriv_dep(Index i);
    replay(const global &orig, global &target);
    void start();
    void stop();
    void add_updatable_derivs(const intervals<Index> &I);
    void clear_deriv();
    void forward(bool inv_tags = true, bool dep_tags = true,
                 Position start = Position(0, 0, 0),
                 const std::vector<bool> &node_filter = std::vector<bool>());
    void reverse(bool dep_tags = true, bool inv_tags = false,
                 Position start = Position(0, 0, 0),
                 const std::vector<bool> &node_filter = std::vector<bool>());
    void forward_sub();
    void reverse_sub();
    void clear_deriv_sub();
  };

  void forward_replay(bool inv_tags = true, bool dep_tags = true);

  void subgraph_cache_ptr() const;
  void set_subgraph(const std::vector<bool> &marks, bool append = false);
  void mark_subgraph(std::vector<bool> &marks);
  void unmark_subgraph(std::vector<bool> &marks);
  void subgraph_trivial();
  void clear_deriv_sub();
  global extract_sub(std::vector<Index> &var_remap, global new_glob = global());
  void extract_sub_inplace(std::vector<bool> marks);
  global extract_sub();

  std::vector<Index> var2op();
  std::vector<bool> var2op(const std::vector<bool> &values);
  std::vector<Index> op2var(const std::vector<Index> &seq);
  std::vector<bool> op2var(const std::vector<bool> &seq_mark);
  std::vector<Index> op2idx(const std::vector<Index> &var_subset,
                            Index NA = (Index)-1);
  std::vector<bool> mark_space(size_t n, const std::vector<Index> ind);
  std::vector<bool> inv_marks();
  std::vector<bool> dep_marks();
  std::vector<bool> subgraph_marks();

  struct append_edges {
    size_t &i;
    const std::vector<bool> &keep_var;
    std::vector<Index> &var2op;
    std::vector<IndexPair> &edges;

    std::vector<bool> op_marks;
    size_t pos;
    append_edges(size_t &i, size_t num_nodes, const std::vector<bool> &keep_var,
                 std::vector<Index> &var2op, std::vector<IndexPair> &edges);
    void operator()(Index dep_j);

    void start_iteration();

    void end_iteration();
  };
  graph build_graph(bool transpose, const std::vector<bool> &keep_var);
  graph forward_graph(std::vector<bool> keep_var = std::vector<bool>(0));
  graph reverse_graph(std::vector<bool> keep_var = std::vector<bool>(0));

  bool identical(const global &other) const;

  template <class T>
  void hash(hash_t &h, T x) const {
    static const size_t n =
        (sizeof(T) / sizeof(hash_t)) + (sizeof(T) % sizeof(hash_t) != 0);
    hash_t buffer[n];
    std::fill(buffer, buffer + n, 0);
    for (size_t i = 0; i < sizeof(x); i++)
      ((char *)buffer)[i] = ((char *)&x)[i];
    hash_t A = 54059;
    hash_t B = 76963;
    for (size_t i = 0; i < n; i++) h = (A * h) ^ (B * buffer[i]);
  }

  hash_t hash() const;

  struct hash_config {
    bool strong_inv;
    bool strong_const;
    bool strong_output;
    bool reduce;
    bool deterministic;
    std::vector<Index> inv_seed;
  };

  std::vector<hash_t> hash_sweep(hash_config cfg) const;
  std::vector<hash_t> hash_sweep(bool weak = true) const;

  void eliminate();

  struct print_config {
    std::string prefix, mark;
    int depth;
    print_config();
  };
  void print(print_config cfg);
  void print();

  template <int ninput_, int noutput_ = 1>
  struct Operator {
    static const bool dynamic = false;
    static const int ninput = ninput_;
    static const int noutput = noutput_;
    static const int independent_variable = false;
    static const int dependent_variable = false;
    static const bool have_input_size_output_size = false;
    static const bool have_increment_decrement = false;
    static const bool have_forward_reverse = true;
    static const bool have_forward_incr_reverse_decr = false;
    static const bool have_forward_mark_reverse_mark = false;
    static const bool have_dependencies = false;
    static const bool allow_remap = true;
    static const bool implicit_dependencies = false;
    static const bool add_static_identifier = false;
    static const bool add_forward_replay_copy = false;
    static const bool have_eval = false;
    static const int max_fuse_depth = 2;
    static const bool is_linear = false;
    static const bool is_constant = false;
    static const bool smart_pointer = false;
    static const bool elimination_protected = false;
    static const bool updating = false;
    void dependencies_updating(Args<> &args, Dependencies &dep) const {}
    OperatorPure *other_fuse(OperatorPure *self, OperatorPure *other) {
      return NULL;
    }
    void *operator_data() { return NULL; }
    void print(print_config cfg) {}
  };
  template <int ninput, int noutput>
  struct DynamicOperator : Operator<ninput, noutput> {
    static const bool dynamic = true;
    static const int max_fuse_depth = 0;
  };
  template <int ninput>
  struct DynamicOutputOperator : Operator<ninput, -1> {
    static const bool dynamic = true;
    static const int max_fuse_depth = 0;
    Index noutput;
  };
  template <int noutput = 1>
  struct DynamicInputOperator : Operator<-1, noutput> {
    static const bool dynamic = true;
    static const int max_fuse_depth = 0;
    Index ninput;
  };
  struct DynamicInputOutputOperator : Operator<-1, -1> {
    static const bool dynamic = true;
    static const int max_fuse_depth = 0;
    Index ninput_, noutput_;
    DynamicInputOutputOperator(Index ninput, Index noutput);
    Index input_size() const;
    Index output_size() const;
    static const bool have_input_size_output_size = true;
  };
  struct UniqueDynamicOperator : Operator<-1, -1> {
    static const bool dynamic = true;
    static const int max_fuse_depth = 0;
    static const bool smart_pointer = false;
    static const bool have_input_size_output_size = true;
  };
  struct SharedDynamicOperator : UniqueDynamicOperator {
    static const bool smart_pointer = true;
  };

  template <class OperatorBase>
  struct AddInputSizeOutputSize : OperatorBase {
    INHERIT_CTOR(AddInputSizeOutputSize, OperatorBase)
    Index input_size() const { return this->ninput; }
    Index output_size() const { return this->noutput; }
    static const bool have_input_size_output_size = true;
  };

  template <class OperatorBase>
  struct AddIncrementDecrement : OperatorBase {
    INHERIT_CTOR(AddIncrementDecrement, OperatorBase)
    void increment(IndexPair &ptr) {
      ptr.first += this->input_size();
      ptr.second += this->output_size();
    }
    void decrement(IndexPair &ptr) {
      ptr.first -= this->input_size();
      ptr.second -= this->output_size();
    }
    static const bool have_increment_decrement = true;
  };

  template <class OperatorBase>
  struct AddForwardReverse : OperatorBase {
    INHERIT_CTOR(AddForwardReverse, OperatorBase)

    template <class Type>
    void forward(ForwardArgs<Type> &args) {
      ForwardArgs<Type> args_cpy(args);
      OperatorBase::forward_incr(args_cpy);
    }
    template <class Type>
    void reverse(ReverseArgs<Type> &args) {
      ReverseArgs<Type> args_cpy(args);
      OperatorBase::increment(args_cpy.ptr);
      OperatorBase::reverse_decr(args_cpy);
    }
    static const bool have_forward_reverse = true;
  };

  template <class OperatorBase>
  struct AddForwardIncrReverseDecr : OperatorBase {
    INHERIT_CTOR(AddForwardIncrReverseDecr, OperatorBase)

    template <class Type>
    void forward_incr(ForwardArgs<Type> &args) {
      OperatorBase::forward(args);
      OperatorBase::increment(args.ptr);
    }

    template <class Type>
    void reverse_decr(ReverseArgs<Type> &args) {
      OperatorBase::decrement(args.ptr);
      OperatorBase::reverse(args);
    }
    static const bool have_forward_incr_reverse_decr = true;
  };

  template <class OperatorBase>
  struct AddForwardMarkReverseMark : OperatorBase {
    INHERIT_CTOR(AddForwardMarkReverseMark, OperatorBase)

    template <class Type>
    void forward(ForwardArgs<Type> &args) {
      OperatorBase::forward(args);
    }
    template <class Type>
    void reverse(ReverseArgs<Type> &args) {
      OperatorBase::reverse(args);
    }

    void forward(ForwardArgs<bool> &args) { args.mark_dense(*this); }
    void reverse(ReverseArgs<bool> &args) { args.mark_dense(*this); }
    static const bool have_forward_mark_reverse_mark = true;
  };

  template <class OperatorBase>
  struct AddDependencies : OperatorBase {
    INHERIT_CTOR(AddDependencies, OperatorBase)
    void dependencies(Args<> &args, Dependencies &dep) const {
      Index ninput_ = this->input_size();
      for (Index j = 0; j < ninput_; j++) dep.push_back(args.input(j));
    }
    static const bool have_dependencies = true;
  };

  template <class OperatorBase, int ninput>
  struct AddForwardFromEval : OperatorBase {};
  template <class OperatorBase>
  struct AddForwardFromEval<OperatorBase, 1> : OperatorBase {
    INHERIT_CTOR(AddForwardFromEval, OperatorBase)
    template <class Type>
    void forward(ForwardArgs<Type> &args) {
      args.y(0) = this->eval(args.x(0));
    }
  };
  template <class OperatorBase>
  struct AddForwardFromEval<OperatorBase, 2> : OperatorBase {
    INHERIT_CTOR(AddForwardFromEval, OperatorBase)
    template <class Type>
    void forward(ForwardArgs<Type> &args) {
      args.y(0) = this->eval(args.x(0), args.x(1));
    }
  };

  template <bool flag, class dummy>
  struct ReferenceCounter {
    void increment() {}
    void decrement() {}
    size_t operator()() const { return 0; }
  };
  template <class dummy>
  struct ReferenceCounter<true, dummy> {
    size_t counter;
    ReferenceCounter() : counter(0) {}
    void increment() { counter++; }
    void decrement() { counter--; }
    size_t operator()() const { return counter; }
  };

  template <bool flag, class Yes, class No>
  struct if_else {};
  template <class Yes, class No>
  struct if_else<true, Yes, No> {
    typedef Yes type;
  };
  template <class Yes, class No>
  struct if_else<false, Yes, No> {
    typedef No type;
  };

  template <class OperatorBase>
  struct CPL {
    static const bool test1 = !OperatorBase::have_eval;
    typedef typename if_else<
        test1, OperatorBase,
        AddForwardFromEval<OperatorBase, OperatorBase::ninput> >::type Result1;

    static const bool test2 = Result1::have_input_size_output_size;
    typedef
        typename if_else<test2, Result1, AddInputSizeOutputSize<Result1> >::type
            Result2;

    static const bool test3 = !Result2::have_dependencies;
    typedef typename if_else<test3, AddDependencies<Result2>, Result2>::type
        Result3;

    static const bool test4 = Result3::have_increment_decrement;
    typedef
        typename if_else<test4, Result3, AddIncrementDecrement<Result3> >::type
            Result4;

    static const bool test5 = Result4::have_forward_mark_reverse_mark;
    typedef typename if_else<test5, Result4,
                             AddForwardMarkReverseMark<Result4> >::type Result5;

    static const bool test6 = Result5::have_forward_reverse &&
                              !Result5::have_forward_incr_reverse_decr;
    typedef typename if_else<test6, AddForwardIncrReverseDecr<Result5>,
                             Result5>::type Result6;

    static const bool test7 = Result6::have_forward_incr_reverse_decr &&
                              !Result6::have_forward_reverse;
    typedef typename if_else<test7, AddForwardReverse<Result6>, Result6>::type
        Result7;

    typedef Result7 type;
  };

  template <class Operator1, class Operator2>
  struct Fused : Operator<Operator1::ninput + Operator2::ninput,
                          Operator1::noutput + Operator2::noutput> {
    typename CPL<Operator1>::type Op1;
    typename CPL<Operator2>::type Op2;
    static const int independent_variable =
        Operator1::independent_variable && Operator2::independent_variable;
    static const int dependent_variable =
        Operator1::dependent_variable && Operator2::dependent_variable;
    static const int max_fuse_depth =
        (Operator1::max_fuse_depth < Operator2::max_fuse_depth
             ? Operator1::max_fuse_depth - 1
             : Operator2::max_fuse_depth - 1);
    static const bool is_linear = Operator1::is_linear && Operator2::is_linear;
    template <class Type>
    void forward_incr(ForwardArgs<Type> &args) {
      Op1.forward_incr(args);
      Op2.forward_incr(args);
    }
    template <class Type>
    void reverse_decr(ReverseArgs<Type> &args) {
      Op2.reverse_decr(args);
      Op1.reverse_decr(args);
    }
    static const bool have_forward_incr_reverse_decr = true;
    static const bool have_forward_reverse = false;
    const char *op_name() { return "Fused"; }
  };
  template <class Operator1>
  struct Rep : DynamicOperator<-1, -1> {
    typename CPL<Operator1>::type Op;
    static const int independent_variable = Operator1::independent_variable;
    static const int dependent_variable = Operator1::dependent_variable;
    static const bool is_linear = Operator1::is_linear;
    Index n;
    Rep(Index n) : n(n) {}
    Index input_size() const { return Operator1::ninput * n; }
    Index output_size() const { return Operator1::noutput * n; }
    static const bool have_input_size_output_size = true;
    template <class Type>
    void forward_incr(ForwardArgs<Type> &args) {
      for (size_t i = 0; i < (size_t)n; i++) Op.forward_incr(args);
    }
    template <class Type>
    void reverse_decr(ReverseArgs<Type> &args) {
      for (size_t i = 0; i < (size_t)n; i++) Op.reverse_decr(args);
    }
    static const bool have_forward_incr_reverse_decr = true;
    static const bool have_forward_reverse = false;
    OperatorPure *compress() {
      TMBAD_ASSERT(false);
      std::vector<Index> &inputs = get_glob()->inputs;
      size_t k = Op.input_size();
      size_t start = inputs.size() - k * n;
      std::valarray<Index> increment(k);
      if (k > 0) {
        for (size_t i = 0; i < (size_t)n - 1; i++) {
          std::valarray<Index> v1(&inputs[start + i * k], k);
          std::valarray<Index> v2(&inputs[start + (i + 1) * k], k);
          if (i == 0) {
            increment = v2 - v1;
          } else {
            bool ok = (increment == (v2 - v1)).min();
            if (!ok) return NULL;
          }
        }
      }

      size_t reduction = (n - 1) * k;
      inputs.resize(inputs.size() - reduction);
      return get_glob()->getOperator<RepCompress<Operator1> >(n, increment);
    }
    OperatorPure *other_fuse(OperatorPure *self, OperatorPure *other) {
      OperatorPure *op1 = get_glob()->getOperator<Operator1>();
      if (op1 == other) {
        this->n++;
        return self;
      }
      return NULL;
    }
    const char *op_name() { return "Rep"; }
  };
  template <class Operator1>
  struct RepCompress : DynamicOperator<-1, -1> {
    static const int independent_variable = Operator1::independent_variable;
    static const int dependent_variable = Operator1::dependent_variable;
    static const bool is_linear = Operator1::is_linear;
    typename CPL<Operator1>::type Op;
    Index n;

    std::valarray<Index> increment_pattern;
    RepCompress(Index n, std::valarray<Index> v) : n(n), increment_pattern(v) {}
    Index input_size() const { return Operator1::ninput; }
    Index output_size() const { return Operator1::noutput * n; }
    static const bool have_input_size_output_size = true;
    template <class Type>
    void forward(ForwardArgs<Type> &args) {
      std::valarray<Index> inputs(input_size());
      for (size_t i = 0; i < inputs.size(); i++) inputs[i] = args.input(i);
      ForwardArgs<Type> args_cpy = args;
      args_cpy.inputs = &inputs[0];
      args_cpy.ptr.first = 0;
      for (size_t i = 0; i < (size_t)n; i++) {
        Op.forward(args_cpy);
        inputs += this->increment_pattern;
        args_cpy.ptr.second += Op.output_size();
      }
    }
    template <class Type>
    void reverse(ReverseArgs<Type> &args) {
      std::valarray<Index> inputs(input_size());
      for (size_t i = 0; i < inputs.size(); i++) inputs[i] = args.input(i);
      inputs += n * this->increment_pattern;
      ReverseArgs<Type> args_cpy = args;
      args_cpy.inputs = &inputs[0];
      args_cpy.ptr.first = 0;
      args_cpy.ptr.second += n * Op.output_size();
      for (size_t i = 0; i < (size_t)n; i++) {
        inputs -= this->increment_pattern;
        args_cpy.ptr.second -= Op.output_size();
        Op.reverse(args_cpy);
      }
    }
    void dependencies(Args<> &args, Dependencies &dep) const {
      std::valarray<Index> inputs(input_size());
      for (size_t i = 0; i < inputs.size(); i++) inputs[i] = args.input(i);
      for (size_t i = 0; i < (size_t)n; i++) {
        dep.insert(dep.end(), &inputs[0], &inputs[0] + inputs.size());
        inputs += this->increment_pattern;
      }
    }
    static const bool have_dependencies = true;
    void forward(ForwardArgs<Writer> &args) {
      std::valarray<Index> inputs(Op.input_size());
      for (size_t i = 0; i < (size_t)Op.input_size(); i++)
        inputs[i] = args.input(i);
      std::valarray<Index> outputs(Op.output_size());
      for (size_t i = 0; i < (size_t)Op.output_size(); i++)
        outputs[i] = args.output(i);
      Writer w;
      int ninp = Op.input_size();
      int nout = Op.output_size();

      w << "for (int count = 0, "
        << "i[" << ninp << "]=" << inputs << ", "
        << "di[" << ninp << "]=" << increment_pattern << ", "
        << "o[" << nout << "]=" << outputs << "; "
        << "count < " << n << "; count++) {\n";

      w << "    ";
      ForwardArgs<Writer> args_cpy = args;
      args_cpy.set_indirect();
      Op.forward(args_cpy);
      w << "\n";

      w << "    ";
      w << "for (int k=0; k<" << ninp << "; k++) i[k] += di[k];\n";
      w << "    ";
      w << "for (int k=0; k<" << nout << "; k++) o[k] += " << nout << ";\n";

      w << "  ";
      w << "}";
    }
    void reverse(ReverseArgs<Writer> &args) {
      std::valarray<Index> inputs(Op.input_size());
      for (size_t i = 0; i < (size_t)Op.input_size(); i++)
        inputs[i] = args.input(i);
      inputs += n * increment_pattern;
      std::valarray<Index> outputs(Op.output_size());
      for (size_t i = 0; i < (size_t)Op.output_size(); i++)
        outputs[i] = args.output(i);
      outputs += n * Op.output_size();
      Writer w;
      int ninp = Op.input_size();
      int nout = Op.output_size();

      w << "for (int count = 0, "
        << "i[" << ninp << "]=" << inputs << ", "
        << "di[" << ninp << "]=" << increment_pattern << ", "
        << "o[" << nout << "]=" << outputs << "; "
        << "count < " << n << "; count++) {\n";

      w << "    ";
      w << "for (int k=0; k<" << ninp << "; k++) i[k] -= di[k];\n";
      w << "    ";
      w << "for (int k=0; k<" << nout << "; k++) o[k] -= " << nout << ";\n";

      w << "    ";
      ReverseArgs<Writer> args_cpy = args;
      args_cpy.set_indirect();
      Op.reverse(args_cpy);
      w << "\n";

      w << "  ";
      w << "}";
    }
    static const bool have_forward_incr_reverse_decr = false;
    static const bool have_forward_reverse = true;
    static const bool have_forward_mark_reverse_mark = true;
    const char *op_name() { return "CRep"; }

    struct operator_data_t {
      OperatorPure *Op;
      Index n;
      std::valarray<Index> ip;
      operator_data_t(const RepCompress &x)
          : Op(get_glob()->getOperator<Operator1>()),
            n(x.n),
            ip(x.increment_pattern) {}
      ~operator_data_t() { Op->deallocate(); }
      bool operator==(const operator_data_t &other) {
        return (Op == other.Op) && (ip.size() == other.ip.size()) &&
               ((ip - other.ip).min() == 0);
      }
    };
    void *operator_data() { return new operator_data_t(*this); }
    OperatorPure *other_fuse(OperatorPure *self, OperatorPure *other) {
      if (this->op_name() == other->op_name()) {
        operator_data_t *p1 =
            static_cast<operator_data_t *>(self->operator_data());
        operator_data_t *p2 =
            static_cast<operator_data_t *>(other->operator_data());
        bool match = (*p1 == *p2);
        int other_n = p2->n;
        delete p1;
        delete p2;
        if (match) {
          std::vector<Index> &inputs = get_glob()->inputs;
          size_t reduction = increment_pattern.size();
          inputs.resize(inputs.size() - reduction);
          this->n += other_n;
          other->deallocate();
          return self;
        }
      }
      return NULL;
    }
  };

  template <class OperatorBase>
  struct Complete : OperatorPure {
    typename CPL<OperatorBase>::type Op;
    INHERIT_CTOR(Complete, Op)
    ~Complete() {}
    void forward(ForwardArgs<Scalar> &args) { Op.forward(args); }
    void reverse(ReverseArgs<Scalar> &args) { Op.reverse(args); }
    void forward_incr(ForwardArgs<Scalar> &args) { Op.forward_incr(args); }
    void reverse_decr(ReverseArgs<Scalar> &args) { Op.reverse_decr(args); }

    void forward(ForwardArgs<Replay> &args) {
      if (Op.add_forward_replay_copy)
        forward_replay_copy(args);
      else
        Op.forward(args);
    }
    void reverse(ReverseArgs<Replay> &args) { Op.reverse(args); }
    void forward_incr(ForwardArgs<Replay> &args) {
      if (Op.add_forward_replay_copy) {
        forward_replay_copy(args);
        increment(args.ptr);
      } else
        Op.forward_incr(args);
    }
    void reverse_decr(ReverseArgs<Replay> &args) { Op.reverse_decr(args); }

    void forward(ForwardArgs<bool> &args) { Op.forward(args); }
    void reverse(ReverseArgs<bool> &args) { Op.reverse(args); }
    void forward_incr(ForwardArgs<bool> &args) { Op.forward_incr(args); }
    void reverse_decr(ReverseArgs<bool> &args) { Op.reverse_decr(args); }
    void forward_incr_mark_dense(ForwardArgs<bool> &args) {
      args.mark_dense(Op);
      Op.increment(args.ptr);
    };

    void forward(ForwardArgs<Writer> &args) { Op.forward(args); }
    void reverse(ReverseArgs<Writer> &args) { Op.reverse(args); }
    void forward_incr(ForwardArgs<Writer> &args) { Op.forward_incr(args); }
    void reverse_decr(ReverseArgs<Writer> &args) { Op.reverse_decr(args); }
    std::vector<ad_plain> operator()(const std::vector<ad_plain> &x) {
      TMBAD_ASSERT2(OperatorBase::dynamic,
                    "Stack to heap copy only allowed for dynamic operators");
      Complete *pOp = new Complete(*this);
      TMBAD_ASSERT2(pOp->ref_count() == 0, "Operator already on the heap");
      pOp->ref_count.increment();
      return get_glob()->add_to_stack<OperatorBase>(pOp, x);
    }
    ad_segment operator()(const ad_segment &x) {
      TMBAD_ASSERT2(OperatorBase::dynamic,
                    "Stack to heap copy only allowed for dynamic operators");
      Complete *pOp = new Complete(*this);
      TMBAD_ASSERT2(pOp->ref_count() == 0, "Operator already on the heap");
      pOp->ref_count.increment();
      return get_glob()->add_to_stack<OperatorBase>(pOp, x);
    }
    ad_segment operator()(const ad_segment &x, const ad_segment &y) {
      TMBAD_ASSERT2(OperatorBase::dynamic,
                    "Stack to heap copy only allowed for dynamic operators");
      Complete *pOp = new Complete(*this);
      TMBAD_ASSERT2(pOp->ref_count() == 0, "Operator already on the heap");
      pOp->ref_count.increment();
      return get_glob()->add_to_stack<OperatorBase>(pOp, x, y);
    }
    template <class T>
    std::vector<T> operator()(const std::vector<T> &x) {
      std::vector<ad_plain> x_(x.begin(), x.end());
      std::vector<ad_plain> y_ = (*this)(x_);
      std::vector<T> y(y_.begin(), y_.end());
      return y;
    }
    void forward_replay_copy(ForwardArgs<Replay> &args) {
      std::vector<ad_plain> x(Op.input_size());
      for (size_t i = 0; i < x.size(); i++) x[i] = args.x(i);
      std::vector<ad_plain> y =
          get_glob()->add_to_stack<OperatorBase>(this->copy(), x);
      for (size_t i = 0; i < y.size(); i++) args.y(i) = y[i];
    }
    void dependencies(Args<> &args, Dependencies &dep) {
      Op.dependencies(args, dep);
    }
    void dependencies_updating(Args<> &args, Dependencies &dep) {
      Op.dependencies_updating(args, dep);
    }
    void increment(IndexPair &ptr) { Op.increment(ptr); }
    void decrement(IndexPair &ptr) { Op.decrement(ptr); }
    Index input_size() { return Op.input_size(); }
    Index output_size() { return Op.output_size(); }
    const char *op_name() { return Op.op_name(); }
    void print(print_config cfg) { Op.print(cfg); }

    template <class Operator_, int depth>
    struct SelfFuse {
      typedef Rep<Operator_> type;
      OperatorPure *operator()() {
        return get_glob()->template getOperator<type>(2);
      }
    };
    template <class Operator_>
    struct SelfFuse<Operator_, 0> {
      OperatorPure *operator()() { return NULL; }
    };
    OperatorPure *self_fuse() {
      return SelfFuse<OperatorBase, OperatorBase::max_fuse_depth>()();
    }
    OperatorPure *other_fuse(OperatorPure *other) {
      return Op.other_fuse(this, other);
    }
    ReferenceCounter<OperatorBase::smart_pointer, void> ref_count;
    OperatorPure *copy() {
      if (Op.smart_pointer) {
        ref_count.increment();
        return this;
      } else if (Op.dynamic)
        return new Complete(*this);
      else
        return this;
    }
    void deallocate() {
      if (!Op.dynamic) return;
      if (Op.smart_pointer) {
        if (ref_count() > 1) {
          ref_count.decrement();
          return;
        }
      }
      delete this;
    }
    op_info info() {
      op_info info(Op);
      return info;
    }
    void *identifier() {
      if (Op.add_static_identifier) {
        static void *id = new char();
        return id;
      } else
        return (void *)this;
    }
    void *operator_data() { return Op.operator_data(); }
    void *incomplete() { return &Op; }
  };

  template <class OperatorBase>
  Complete<OperatorBase> *getOperator() const {
    return constructOperator<Complete<OperatorBase>, OperatorBase::dynamic>()();
  }
  template <class OperatorBase, class T1>
  Complete<OperatorBase> *getOperator(const T1 &x1) const {
    return constructOperator<Complete<OperatorBase>, OperatorBase::dynamic>()(
        x1);
  }
  template <class OperatorBase, class T1, class T2>
  Complete<OperatorBase> *getOperator(const T1 &x1, const T2 &x2) const {
    return constructOperator<Complete<OperatorBase>, OperatorBase::dynamic>()(
        x1, x2);
  }
  template <class OperatorBase, class T1, class T2, class T3>
  Complete<OperatorBase> *getOperator(const T1 &x1, const T2 &x2,
                                      const T3 &x3) const {
    return constructOperator<Complete<OperatorBase>, OperatorBase::dynamic>()(
        x1, x2, x3);
  }
  template <class OperatorBase, class T1, class T2, class T3, class T4>
  Complete<OperatorBase> *getOperator(const T1 &x1, const T2 &x2, const T3 &x3,
                                      const T4 &x4) const {
    return constructOperator<Complete<OperatorBase>, OperatorBase::dynamic>()(
        x1, x2, x3, x4);
  }
  struct InvOp : Operator<0> {
    static const int independent_variable = true;
    template <class Type>
    void forward(ForwardArgs<Type> &args) {}
    template <class Type>
    void reverse(ReverseArgs<Type> &args) {}
    const char *op_name();
  };

  struct DepOp : Operator<1> {
    static const bool is_linear = true;
    static const int dependent_variable = true;
    static const bool have_eval = true;
    template <class Type>
    Type eval(Type x0) {
      return x0;
    }
    template <class Type>
    void reverse(ReverseArgs<Type> &args) {
      args.dx(0) += args.dy(0);
    }
    const char *op_name();
  };

  struct ConstOp : Operator<0, 1> {
    static const bool is_linear = true;
    static const bool is_constant = true;
    template <class Type>
    void forward(ForwardArgs<Type> &args) {}
    void forward(ForwardArgs<Replay> &args);
    template <class Type>
    void reverse(ReverseArgs<Type> &args) {}
    const char *op_name();
    void forward(ForwardArgs<Writer> &args);
  };
  struct DataOp : DynamicOutputOperator<0> {
    typedef DynamicOutputOperator<0> Base;
    static const bool is_linear = true;
    DataOp(Index n);
    template <class Type>
    void forward(ForwardArgs<Type> &args) {}
    template <class Type>
    void reverse(ReverseArgs<Type> &args) {}
    const char *op_name();
    void forward(ForwardArgs<Writer> &args);
  };
  struct ZeroOp : DynamicOutputOperator<0> {
    typedef DynamicOutputOperator<0> Base;
    static const bool add_forward_replay_copy = true;
    ZeroOp(Index n);
    template <class Type>
    void forward(ForwardArgs<Type> &args) {
      for (Index i = 0; i < Base::noutput; i++) args.y(i) = Type(0);
    }
    template <class Type>
    void reverse(ReverseArgs<Type> &args) {}
    const char *op_name();
    void forward(ForwardArgs<Writer> &args);
    void operator()(Replay *x, Index n);
  };
  struct NullOp : Operator<0, 0> {
    NullOp();
    const char *op_name();
    template <class T>
    void forward(ForwardArgs<T> &args) {}
    template <class T>
    void reverse(ReverseArgs<T> &args) {}
  };
  struct NullOp2 : DynamicInputOutputOperator {
    NullOp2(Index ninput, Index noutput);
    const char *op_name();
    template <class T>
    void forward(ForwardArgs<T> &args) {}
    template <class T>
    void reverse(ReverseArgs<T> &args) {}
  };
  struct RefOp : DynamicOperator<0, 1> {
    static const bool dynamic = true;
    global *glob;
    Index i;
    RefOp(global *glob, Index i);
    void forward(ForwardArgs<Scalar> &args);
    void forward(ForwardArgs<Replay> &args);
    template <class Type>
    void reverse(ReverseArgs<Type> &args) {
      TMBAD_ASSERT2(false,
                    "Reverse mode updates are forbidden until all references "
                    "are resolved");
    }
    void reverse(ReverseArgs<Replay> &args);
    const char *op_name();
  };

  typedef Operator<1> UnaryOperator;
  typedef Operator<2> BinaryOperator;

  OperatorPure *Fuse(OperatorPure *Op1, OperatorPure *Op2);

  static bool fuse;

  void set_fuse(bool flag);

  void add_to_opstack(OperatorPure *pOp);
  template <class OperatorBase>
  ad_plain add_to_stack(Scalar result = 0) {
    ad_plain ans;
    ans.index = this->values.size();

    this->values.push_back(result);

    Complete<OperatorBase> *pOp = this->template getOperator<OperatorBase>();
    add_to_opstack(pOp);

    TMBAD_ASSERT(!TMBAD_INDEX_OVERFLOW(values.size()));
    return ans;
  }
  template <class OperatorBase>
  ad_plain add_to_stack(const ad_plain &x) {
    ad_plain ans;
    ans.index = this->values.size();

    this->values.push_back(OperatorBase().eval(x.Value()));

    this->inputs.push_back(x.index);

    Complete<OperatorBase> *pOp = this->template getOperator<OperatorBase>();
    add_to_opstack(pOp);

    TMBAD_ASSERT(!TMBAD_INDEX_OVERFLOW(values.size()));
    TMBAD_ASSERT(!TMBAD_INDEX_OVERFLOW(inputs.size()));
    return ans;
  }
  template <class OperatorBase>
  ad_plain add_to_stack(const ad_plain &x, const ad_plain &y) {
    ad_plain ans;
    ans.index = this->values.size();

    this->values.push_back(OperatorBase().eval(x.Value(), y.Value()));

    this->inputs.push_back(x.index);
    this->inputs.push_back(y.index);

    Complete<OperatorBase> *pOp = this->template getOperator<OperatorBase>();
    add_to_opstack(pOp);

    TMBAD_ASSERT(!TMBAD_INDEX_OVERFLOW(values.size()));
    TMBAD_ASSERT(!TMBAD_INDEX_OVERFLOW(inputs.size()));
    return ans;
  }
  template <class OperatorBase>
  ad_segment add_to_stack(ad_segment lhs, ad_segment rhs,
                          ad_segment more = ad_segment()) {
    IndexPair ptr((Index)inputs.size(), (Index)values.size());
    Complete<OperatorBase> *pOp =
        this->template getOperator<OperatorBase>(lhs, rhs);
    size_t n = pOp->output_size();
    ad_segment ans(values.size(), n);
    inputs.push_back(lhs.index());
    inputs.push_back(rhs.index());
    if (more.size() > 0) inputs.push_back(more.index());
    opstack.push_back(pOp);
    values.resize(values.size() + n);
    ForwardArgs<Scalar> args(inputs, values, this);
    args.ptr = ptr;
    pOp->forward(args);

    TMBAD_ASSERT(!TMBAD_INDEX_OVERFLOW(values.size()));
    TMBAD_ASSERT(!TMBAD_INDEX_OVERFLOW(inputs.size()));
    return ans;
  }

  template <class OperatorBase>
  ad_segment add_to_stack(Complete<OperatorBase> *pOp, ad_segment lhs,
                          ad_segment rhs = ad_segment()) {
    static_assert(
        OperatorBase::dynamic,
        "Unlikely that you want to use this method for static operators?");
    static_assert(
        OperatorBase::ninput == 0 || OperatorBase::implicit_dependencies,
        "Operators with pointer inputs should always implement "
        "'implicit_dependencies'");

    IndexPair ptr((Index)inputs.size(), (Index)values.size());
    size_t n = pOp->output_size();
    ad_segment ans(values.size(), n);
    TMBAD_ASSERT((Index)(lhs.size() > 0) + (Index)(rhs.size() > 0) ==
                 pOp->input_size());
    if (lhs.size() > 0) inputs.push_back(lhs.index());
    if (rhs.size() > 0) inputs.push_back(rhs.index());
    opstack.push_back(pOp);
    values.resize(values.size() + n);
    ForwardArgs<Scalar> args(inputs, values, this);
    args.ptr = ptr;
    pOp->forward(args);

    TMBAD_ASSERT(!TMBAD_INDEX_OVERFLOW(values.size()));
    TMBAD_ASSERT(!TMBAD_INDEX_OVERFLOW(inputs.size()));
    return ans;
  }
  template <class OperatorBase>
  std::vector<ad_plain> add_to_stack(OperatorPure *pOp,
                                     const std::vector<ad_plain> &x) {
    IndexPair ptr((Index)inputs.size(), (Index)values.size());
    size_t m = pOp->input_size();
    size_t n = pOp->output_size();
    ad_segment ans(values.size(), n);
    for (size_t i = 0; i < m; i++) inputs.push_back(x[i].index);
    opstack.push_back(pOp);
    values.resize(values.size() + n);
    ForwardArgs<Scalar> args(inputs, values, this);
    args.ptr = ptr;
    pOp->forward(args);

    TMBAD_ASSERT(!TMBAD_INDEX_OVERFLOW(values.size()));
    TMBAD_ASSERT(!TMBAD_INDEX_OVERFLOW(inputs.size()));
    std::vector<ad_plain> out(n);
    for (size_t i = 0; i < n; i++) out[i].index = ans.index() + i;
    return out;
  }

  struct ad_plain {
    Index index;
    static const Index NA = (Index)-1;
    bool initialized() const;
    bool on_some_tape() const;
    void addToTape() const;
    global *glob() const;
    void override_by(const ad_plain &x) const;

    ad_plain();

    ad_plain(Scalar x);
    ad_plain(ad_aug x);

    struct CopyOp : Operator<1> {
      static const bool have_eval = true;
      template <class Type>
      Type eval(Type x0) {
        return x0;
      }
      Replay eval(Replay x0);
      template <class Type>
      void reverse(ReverseArgs<Type> &args) {
        args.dx(0) += args.dy(0);
      }
      const char *op_name();
    };
    ad_plain copy() const;
    struct ValOp : Operator<1> {
      static const bool have_dependencies = true;
      static const bool have_eval = true;
      template <class Type>
      Type eval(Type x0) {
        return x0;
      }
      Replay eval(Replay x0);
      template <class Type>
      void reverse(ReverseArgs<Type> &args) {}
      void dependencies(Args<> &args, Dependencies &dep) const;
      const char *op_name();
    };
    ad_plain copy0() const;

    template <bool left_var, bool right_var>
    struct AddOp_ : BinaryOperator {
      static const bool is_linear = true;
      static const bool have_eval = true;
      template <class Type>
      Type eval(Type x0, Type x1) {
        return x0 + x1;
      }
      template <class Type>
      void reverse(ReverseArgs<Type> &args) {
        if (left_var) args.dx(0) += args.dy(0);
        if (right_var) args.dx(1) += args.dy(0);
      }
      const char *op_name() { return "AddOp"; }
      OperatorPure *other_fuse(OperatorPure *self, OperatorPure *other) {
        if (other == get_glob()->getOperator<MulOp>()) {
          return get_glob()->getOperator<Fused<AddOp_, MulOp> >();
        }
        return NULL;
      }
    };
    typedef AddOp_<true, true> AddOp;
    ad_plain operator+(const ad_plain &other) const;

    template <bool left_var, bool right_var>
    struct SubOp_ : BinaryOperator {
      static const bool is_linear = true;
      static const bool have_eval = true;
      template <class Type>
      Type eval(Type x0, Type x1) {
        return x0 - x1;
      }
      template <class Type>
      void reverse(ReverseArgs<Type> &args) {
        if (left_var) args.dx(0) += args.dy(0);
        if (right_var) args.dx(1) -= args.dy(0);
      }
      const char *op_name() { return "SubOp"; }
    };
    typedef SubOp_<true, true> SubOp;
    ad_plain operator-(const ad_plain &other) const;

    template <bool left_var, bool right_var>
    struct MulOp_ : BinaryOperator {
      static const bool have_eval = true;
      static const bool is_linear = !left_var || !right_var;
      template <class Type>
      Type eval(Type x0, Type x1) {
        return x0 * x1;
      }
      template <class Type>
      void reverse(ReverseArgs<Type> &args) {
        if (left_var) args.dx(0) += args.x(1) * args.dy(0);
        if (right_var) args.dx(1) += args.x(0) * args.dy(0);
      }
      const char *op_name() { return "MulOp"; }
    };
    typedef MulOp_<true, true> MulOp;
    ad_plain operator*(const ad_plain &other) const;
    ad_plain operator*(const Scalar &other) const;

    template <bool left_var, bool right_var>
    struct DivOp_ : BinaryOperator {
      static const bool have_eval = true;
      template <class Type>
      Type eval(Type x0, Type x1) {
        return x0 / x1;
      }
      template <class Type>
      void reverse(ReverseArgs<Type> &args) {
        Type tmp0 = args.dy(0) / args.x(1);
        if (left_var) args.dx(0) += tmp0;
        if (right_var) args.dx(1) -= args.y(0) * tmp0;
      }
      const char *op_name() { return "DivOp"; }
    };
    typedef DivOp_<true, true> DivOp;
    ad_plain operator/(const ad_plain &other) const;

    struct NegOp : UnaryOperator {
      static const bool is_linear = true;
      static const bool have_eval = true;
      template <class Type>
      Type eval(Type x0) {
        return -x0;
      }
      template <class Type>
      void reverse(ReverseArgs<Type> &args) {
        args.dx(0) -= args.dy(0);
      }
      const char *op_name();
    };
    ad_plain operator-() const;

    ad_plain &operator+=(const ad_plain &other);
    ad_plain &operator-=(const ad_plain &other);
    ad_plain &operator*=(const ad_plain &other);
    ad_plain &operator/=(const ad_plain &other);

    void Dependent();

    void Independent();
    Scalar &Value();
    Scalar Value() const;
    Scalar Value(global *glob) const;
    Scalar &Deriv();
  };
  global *parent_glob;
  bool in_use;
  void ad_start();
  void ad_stop();
  void Independent(std::vector<ad_plain> &x);
  struct ad_segment {
    ad_plain x;
    size_t n;
    size_t c;
    ad_segment();
    ad_segment(ad_plain x, size_t n);
    ad_segment(ad_aug x);
    ad_segment(Scalar x);
    ad_segment(Index idx, size_t n);
    ad_segment(ad_plain x, size_t r, size_t c);
    ad_segment(Replay *x, size_t n, bool zero_check = false);
    bool identicalZero();
    bool all_on_active_tape(Replay *x, size_t n);
    bool is_contiguous(Replay *x, size_t n);
    bool all_zero(Replay *x, size_t n);
    bool all_constant(Replay *x, size_t n);
    size_t size() const;
    size_t rows() const;
    size_t cols() const;

    ad_plain operator[](size_t i) const;
    ad_plain offset() const;
    Index index() const;
  };
  struct ad_aug {
    mutable ad_plain taped_value;
    TMBAD_UNION_OR_STRUCT {
      Scalar value;
      mutable global *glob;
    }
    data;
    bool on_some_tape() const;
    bool on_active_tape() const;
    bool ontape() const;
    bool constant() const;
    Index index() const;
    global *glob() const;
    Scalar Value() const;
    ad_aug();
    ad_aug(Scalar x);
    ad_aug(ad_plain x);
    void addToTape() const;
    void override_by(const ad_plain &x) const;
    bool in_context_stack(global *glob) const;
    ad_aug copy() const;
    ad_aug copy0() const;
    bool identicalZero() const;
    bool identicalOne() const;
    bool bothConstant(const ad_aug &other) const;
    bool identical(const ad_aug &other) const;
    ad_aug operator+(const ad_aug &other) const;
    ad_aug operator-(const ad_aug &other) const;
    ad_aug operator-() const;
    ad_aug operator*(const ad_aug &other) const;
    ad_aug operator/(const ad_aug &other) const;
    ad_aug &operator+=(const ad_aug &other);
    ad_aug &operator-=(const ad_aug &other);
    ad_aug &operator*=(const ad_aug &other);
    ad_aug &operator/=(const ad_aug &other);
    void Dependent();
    void Independent();
    Scalar &Value();
    Scalar &Deriv();
  };
  void Independent(std::vector<ad_aug> &x);
};

template <class S, class T>
std::ostream &operator<<(std::ostream &os, const std::pair<S, T> &x) {
  os << "(" << x.first << ", " << x.second << ")";
  return os;
}

std::ostream &operator<<(std::ostream &os, const global::ad_plain &x);
std::ostream &operator<<(std::ostream &os, const global::ad_aug &x);

template <class T>
struct adaptive : T {
  INHERIT_CTOR(adaptive, T)
  bool operator==(const T &other) const {
    return this->Value() == other.Value();
  }
  bool operator!=(const T &other) const {
    return this->Value() != other.Value();
  }
  bool operator>=(const T &other) const {
    return this->Value() >= other.Value();
  }
  bool operator<=(const T &other) const {
    return this->Value() <= other.Value();
  }
  bool operator<(const T &other) const { return this->Value() < other.Value(); }
  bool operator>(const T &other) const { return this->Value() > other.Value(); }

  adaptive operator+(const T &other) const {
    return adaptive(T(*this) + other);
  }
  adaptive operator-(const T &other) const {
    return adaptive(T(*this) - other);
  }
  adaptive operator*(const T &other) const {
    return adaptive(T(*this) * other);
  }
  adaptive operator/(const T &other) const {
    return adaptive(T(*this) / other);
  }

  adaptive operator-() const { return adaptive(-(T(*this))); }
};

typedef global::ad_plain ad_plain;
typedef global::ad_aug ad_aug;
typedef global::Replay Replay;
typedef adaptive<ad_aug> ad_adapt;
typedef global::OperatorPure OperatorPure;
typedef global::operation_stack operation_stack;
typedef global::print_config print_config;
struct ad_plain_index : ad_plain {
  ad_plain_index(const Index &i);
  ad_plain_index(const ad_plain &x);
};
struct ad_aug_index : ad_aug {
  ad_aug_index(const Index &i);
  ad_aug_index(const ad_aug &x);
  ad_aug_index(const ad_plain &x);
};

template <class T>
void Independent(std::vector<T> &x) {
  for (size_t i = 0; i < x.size(); i++) x[i].Independent();
}
template <class T>
void Dependent(std::vector<T> &x) {
  for (size_t i = 0; i < x.size(); i++) x[i].Dependent();
}
template <class T>
Scalar Value(T x) {
  return x.Value();
}
Scalar Value(Scalar x);

template <class V>
bool isContiguous(V &x) {
  bool ok = true;
  Index j_previous;
  for (size_t i = 0; i < (size_t)x.size(); i++) {
    if (!x[i].on_some_tape()) {
      ok = false;
      break;
    }
    Index j = ad_plain(x[i]).index;
    if (i > 0) {
      if (j != j_previous + 1) {
        ok = false;
        break;
      }
    }
    j_previous = j;
  }
  return ok;
}
template <class V>
V getContiguous(const V &x) {
  V y(x.size());
  for (size_t i = 0; i < (size_t)x.size(); i++) y[i] = x[i].copy();
  return y;
}
template <class V>
void forceContiguous(V &x) {
  if (!isContiguous(x)) x = getContiguous(x);
}
ad_aug operator+(const double &x, const ad_aug &y);
ad_aug operator-(const double &x, const ad_aug &y);
ad_aug operator*(const double &x, const ad_aug &y);
ad_aug operator/(const double &x, const ad_aug &y);

bool operator<(const double &x, const ad_adapt &y);
bool operator<=(const double &x, const ad_adapt &y);
bool operator>(const double &x, const ad_adapt &y);
bool operator>=(const double &x, const ad_adapt &y);
bool operator==(const double &x, const ad_adapt &y);
bool operator!=(const double &x, const ad_adapt &y);
using ::round;
using ::trunc;
using std::ceil;
using std::floor;
Writer floor(const Writer &x);
struct FloorOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return floor(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {}
  const char *op_name();
};
ad_plain floor(const ad_plain &x);
ad_aug floor(const ad_aug &x);
Writer ceil(const Writer &x);
struct CeilOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return ceil(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {}
  const char *op_name();
};
ad_plain ceil(const ad_plain &x);
ad_aug ceil(const ad_aug &x);
Writer trunc(const Writer &x);
struct TruncOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return trunc(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {}
  const char *op_name();
};
ad_plain trunc(const ad_plain &x);
ad_aug trunc(const ad_aug &x);
Writer round(const Writer &x);
struct RoundOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return round(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {}
  const char *op_name();
};
ad_plain round(const ad_plain &x);
ad_aug round(const ad_aug &x);

double sign(const double &x);
Writer sign(const Writer &x);
struct SignOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return sign(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {}
  const char *op_name();
};
ad_plain sign(const ad_plain &x);
ad_aug sign(const ad_aug &x);

double ge0(const double &x);
double lt0(const double &x);
Writer ge0(const Writer &x);
struct Ge0Op : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return ge0(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {}
  const char *op_name();
};
ad_plain ge0(const ad_plain &x);
ad_aug ge0(const ad_aug &x);
Writer lt0(const Writer &x);
struct Lt0Op : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return lt0(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {}
  const char *op_name();
};
ad_plain lt0(const ad_plain &x);
ad_aug lt0(const ad_aug &x);
using ::expm1;
using ::fabs;
using ::log1p;
using std::acos;
using std::acosh;
using std::asin;
using std::asinh;
using std::atan;
using std::atanh;
using std::cos;
using std::cosh;
using std::exp;
using std::log;
using std::sin;
using std::sinh;
using std::sqrt;
using std::tan;
using std::tanh;

Writer fabs(const Writer &x);
struct AbsOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return fabs(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * sign(args.x(0));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain fabs(const ad_plain &x);
ad_aug fabs(const ad_aug &x);
ad_adapt fabs(const ad_adapt &x);
Writer cos(const Writer &x);
ad_aug cos(const ad_aug &x);
Writer sin(const Writer &x);
struct SinOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return sin(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * cos(args.x(0));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain sin(const ad_plain &x);
ad_aug sin(const ad_aug &x);
ad_adapt sin(const ad_adapt &x);
Writer cos(const Writer &x);
struct CosOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return cos(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * -sin(args.x(0));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain cos(const ad_plain &x);
ad_aug cos(const ad_aug &x);
ad_adapt cos(const ad_adapt &x);
Writer exp(const Writer &x);
struct ExpOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return exp(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * args.y(0);
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain exp(const ad_plain &x);
ad_aug exp(const ad_aug &x);
ad_adapt exp(const ad_adapt &x);
Writer log(const Writer &x);
struct LogOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return log(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * Type(1.) / args.x(0);
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain log(const ad_plain &x);
ad_aug log(const ad_aug &x);
ad_adapt log(const ad_adapt &x);
Writer sqrt(const Writer &x);
struct SqrtOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return sqrt(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * Type(0.5) / args.y(0);
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain sqrt(const ad_plain &x);
ad_aug sqrt(const ad_aug &x);
ad_adapt sqrt(const ad_adapt &x);
Writer tan(const Writer &x);
struct TanOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return tan(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * Type(1.) / (cos(args.x(0)) * cos(args.x(0)));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain tan(const ad_plain &x);
ad_aug tan(const ad_aug &x);
ad_adapt tan(const ad_adapt &x);
Writer cosh(const Writer &x);
ad_aug cosh(const ad_aug &x);
Writer sinh(const Writer &x);
struct SinhOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return sinh(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * cosh(args.x(0));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain sinh(const ad_plain &x);
ad_aug sinh(const ad_aug &x);
ad_adapt sinh(const ad_adapt &x);
Writer cosh(const Writer &x);
struct CoshOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return cosh(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * sinh(args.x(0));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain cosh(const ad_plain &x);
ad_aug cosh(const ad_aug &x);
ad_adapt cosh(const ad_adapt &x);
Writer tanh(const Writer &x);
struct TanhOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return tanh(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * Type(1.) / (cosh(args.x(0)) * cosh(args.x(0)));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain tanh(const ad_plain &x);
ad_aug tanh(const ad_aug &x);
ad_adapt tanh(const ad_adapt &x);
Writer expm1(const Writer &x);
struct Expm1 : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return expm1(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * args.y(0) + Type(1.);
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain expm1(const ad_plain &x);
ad_aug expm1(const ad_aug &x);
ad_adapt expm1(const ad_adapt &x);
Writer log1p(const Writer &x);
struct Log1p : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return log1p(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * Type(1.) / (args.x(0) + Type(1.));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain log1p(const ad_plain &x);
ad_aug log1p(const ad_aug &x);
ad_adapt log1p(const ad_adapt &x);
Writer asin(const Writer &x);
struct AsinOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return asin(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) +=
        args.dy(0) * Type(1.) / sqrt(Type(1.) - args.x(0) * args.x(0));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain asin(const ad_plain &x);
ad_aug asin(const ad_aug &x);
ad_adapt asin(const ad_adapt &x);
Writer acos(const Writer &x);
struct AcosOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return acos(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) +=
        args.dy(0) * Type(-1.) / sqrt(Type(1.) - args.x(0) * args.x(0));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain acos(const ad_plain &x);
ad_aug acos(const ad_aug &x);
ad_adapt acos(const ad_adapt &x);
Writer atan(const Writer &x);
struct AtanOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return atan(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * Type(1.) / (Type(1.) + args.x(0) * args.x(0));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain atan(const ad_plain &x);
ad_aug atan(const ad_aug &x);
ad_adapt atan(const ad_adapt &x);
Writer asinh(const Writer &x);
struct AsinhOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return asinh(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) +=
        args.dy(0) * Type(1.) / sqrt(args.x(0) * args.x(0) + Type(1.));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain asinh(const ad_plain &x);
ad_aug asinh(const ad_aug &x);
ad_adapt asinh(const ad_adapt &x);
Writer acosh(const Writer &x);
struct AcoshOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return acosh(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) +=
        args.dy(0) * Type(1.) / sqrt(args.x(0) * args.x(0) - Type(1.));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain acosh(const ad_plain &x);
ad_aug acosh(const ad_aug &x);
ad_adapt acosh(const ad_adapt &x);
Writer atanh(const Writer &x);
struct AtanhOp : global::UnaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x) {
    return atanh(x);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * Type(1.) / (Type(1) - args.x(0) * args.x(0));
  }
  void reverse(ReverseArgs<Scalar> &args);
  const char *op_name();
};
ad_plain atanh(const ad_plain &x);
ad_aug atanh(const ad_aug &x);
ad_adapt atanh(const ad_adapt &x);

template <class T>
T abs(const T &x) {
  return fabs(x);
}
using std::pow;
Writer pow(const Writer &x1, const Writer &x2);
struct PowOp : global::BinaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x1, Type x2) {
    return pow(x1, x2);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * args.x(1) * pow(args.x(0), args.x(1) - Type(1.));
    args.dx(1) += args.dy(0) * args.y(0) * log(args.x(0));
  }
  const char *op_name();
};
ad_plain pow(const ad_plain &x1, const ad_plain &x2);
ad_aug pow(const ad_aug &x1, const ad_aug &x2);
ad_adapt pow(const ad_adapt &x1, const ad_adapt &x2);
using std::atan2;
Writer atan2(const Writer &x1, const Writer &x2);
struct Atan2 : global::BinaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x1, Type x2) {
    return atan2(x1, x2);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * args.x(1) /
                  (args.x(0) * args.x(0) + args.x(1) * args.x(1));
    args.dx(1) += args.dy(0) * -args.x(0) /
                  (args.x(0) * args.x(0) + args.x(1) * args.x(1));
  }
  const char *op_name();
};
ad_plain atan2(const ad_plain &x1, const ad_plain &x2);
ad_aug atan2(const ad_aug &x1, const ad_aug &x2);
ad_adapt atan2(const ad_adapt &x1, const ad_adapt &x2);
using std::max;
Writer max(const Writer &x1, const Writer &x2);
struct MaxOp : global::BinaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x1, Type x2) {
    return max(x1, x2);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * ge0(args.x(0) - args.x(1));
    args.dx(1) += args.dy(0) * lt0(args.x(0) - args.x(1));
  }
  const char *op_name();
};
ad_plain max(const ad_plain &x1, const ad_plain &x2);
ad_aug max(const ad_aug &x1, const ad_aug &x2);
ad_adapt max(const ad_adapt &x1, const ad_adapt &x2);

using std::min;
Writer min(const Writer &x1, const Writer &x2);
struct MinOp : global::BinaryOperator {
  static const bool have_eval = true;
  template <class Type>
  Type eval(Type x1, Type x2) {
    return min(x1, x2);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    args.dx(0) += args.dy(0) * ge0(args.x(1) - args.x(0));
    args.dx(1) += args.dy(0) * lt0(args.x(1) - args.x(0));
  }
  const char *op_name();
};
ad_plain min(const ad_plain &x1, const ad_plain &x2);
ad_aug min(const ad_aug &x1, const ad_aug &x2);
ad_adapt min(const ad_adapt &x1, const ad_adapt &x2);
Replay CondExpEq(const Replay &x0, const Replay &x1, const Replay &x2,
                 const Replay &x3);
struct CondExpEqOp : global::Operator<4, 1> {
  void forward(ForwardArgs<Scalar> &args);
  void reverse(ReverseArgs<Scalar> &args);
  void forward(ForwardArgs<Replay> &args);
  void reverse(ReverseArgs<Replay> &args);
  void forward(ForwardArgs<Writer> &args);
  void reverse(ReverseArgs<Writer> &args);
  template <class Type>
  void forward(ForwardArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  const char *op_name();
};
Scalar CondExpEq(const Scalar &x0, const Scalar &x1, const Scalar &x2,
                 const Scalar &x3);
ad_plain CondExpEq(const ad_plain &x0, const ad_plain &x1, const ad_plain &x2,
                   const ad_plain &x3);
ad_aug CondExpEq(const ad_aug &x0, const ad_aug &x1, const ad_aug &x2,
                 const ad_aug &x3);
Replay CondExpNe(const Replay &x0, const Replay &x1, const Replay &x2,
                 const Replay &x3);
struct CondExpNeOp : global::Operator<4, 1> {
  void forward(ForwardArgs<Scalar> &args);
  void reverse(ReverseArgs<Scalar> &args);
  void forward(ForwardArgs<Replay> &args);
  void reverse(ReverseArgs<Replay> &args);
  void forward(ForwardArgs<Writer> &args);
  void reverse(ReverseArgs<Writer> &args);
  template <class Type>
  void forward(ForwardArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  const char *op_name();
};
Scalar CondExpNe(const Scalar &x0, const Scalar &x1, const Scalar &x2,
                 const Scalar &x3);
ad_plain CondExpNe(const ad_plain &x0, const ad_plain &x1, const ad_plain &x2,
                   const ad_plain &x3);
ad_aug CondExpNe(const ad_aug &x0, const ad_aug &x1, const ad_aug &x2,
                 const ad_aug &x3);
Replay CondExpGt(const Replay &x0, const Replay &x1, const Replay &x2,
                 const Replay &x3);
struct CondExpGtOp : global::Operator<4, 1> {
  void forward(ForwardArgs<Scalar> &args);
  void reverse(ReverseArgs<Scalar> &args);
  void forward(ForwardArgs<Replay> &args);
  void reverse(ReverseArgs<Replay> &args);
  void forward(ForwardArgs<Writer> &args);
  void reverse(ReverseArgs<Writer> &args);
  template <class Type>
  void forward(ForwardArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  const char *op_name();
};
Scalar CondExpGt(const Scalar &x0, const Scalar &x1, const Scalar &x2,
                 const Scalar &x3);
ad_plain CondExpGt(const ad_plain &x0, const ad_plain &x1, const ad_plain &x2,
                   const ad_plain &x3);
ad_aug CondExpGt(const ad_aug &x0, const ad_aug &x1, const ad_aug &x2,
                 const ad_aug &x3);
Replay CondExpLt(const Replay &x0, const Replay &x1, const Replay &x2,
                 const Replay &x3);
struct CondExpLtOp : global::Operator<4, 1> {
  void forward(ForwardArgs<Scalar> &args);
  void reverse(ReverseArgs<Scalar> &args);
  void forward(ForwardArgs<Replay> &args);
  void reverse(ReverseArgs<Replay> &args);
  void forward(ForwardArgs<Writer> &args);
  void reverse(ReverseArgs<Writer> &args);
  template <class Type>
  void forward(ForwardArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  const char *op_name();
};
Scalar CondExpLt(const Scalar &x0, const Scalar &x1, const Scalar &x2,
                 const Scalar &x3);
ad_plain CondExpLt(const ad_plain &x0, const ad_plain &x1, const ad_plain &x2,
                   const ad_plain &x3);
ad_aug CondExpLt(const ad_aug &x0, const ad_aug &x1, const ad_aug &x2,
                 const ad_aug &x3);
Replay CondExpGe(const Replay &x0, const Replay &x1, const Replay &x2,
                 const Replay &x3);
struct CondExpGeOp : global::Operator<4, 1> {
  void forward(ForwardArgs<Scalar> &args);
  void reverse(ReverseArgs<Scalar> &args);
  void forward(ForwardArgs<Replay> &args);
  void reverse(ReverseArgs<Replay> &args);
  void forward(ForwardArgs<Writer> &args);
  void reverse(ReverseArgs<Writer> &args);
  template <class Type>
  void forward(ForwardArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  const char *op_name();
};
Scalar CondExpGe(const Scalar &x0, const Scalar &x1, const Scalar &x2,
                 const Scalar &x3);
ad_plain CondExpGe(const ad_plain &x0, const ad_plain &x1, const ad_plain &x2,
                   const ad_plain &x3);
ad_aug CondExpGe(const ad_aug &x0, const ad_aug &x1, const ad_aug &x2,
                 const ad_aug &x3);
Replay CondExpLe(const Replay &x0, const Replay &x1, const Replay &x2,
                 const Replay &x3);
struct CondExpLeOp : global::Operator<4, 1> {
  void forward(ForwardArgs<Scalar> &args);
  void reverse(ReverseArgs<Scalar> &args);
  void forward(ForwardArgs<Replay> &args);
  void reverse(ReverseArgs<Replay> &args);
  void forward(ForwardArgs<Writer> &args);
  void reverse(ReverseArgs<Writer> &args);
  template <class Type>
  void forward(ForwardArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    TMBAD_ASSERT(false);
  }
  const char *op_name();
};
Scalar CondExpLe(const Scalar &x0, const Scalar &x1, const Scalar &x2,
                 const Scalar &x3);
ad_plain CondExpLe(const ad_plain &x0, const ad_plain &x1, const ad_plain &x2,
                   const ad_plain &x3);
ad_aug CondExpLe(const ad_aug &x0, const ad_aug &x1, const ad_aug &x2,
                 const ad_aug &x3);

template <class Info>
struct InfoOp : global::DynamicOperator<-1, 0> {
  Index n;
  Info info;
  InfoOp(Index n, Info info) : n(n), info(info) {}
  static const bool elimination_protected = true;
  static const bool add_forward_replay_copy = true;
  static const bool have_input_size_output_size = true;
  template <class Type>
  void forward(ForwardArgs<Type> &args) {}
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {}
  Index input_size() const { return n; }
  Index output_size() const { return 0; }
  const char *op_name() { return "InfoOp"; }
  void print(global::print_config cfg) {
    Rcout << cfg.prefix << info << std::endl;
  }
  void *operator_data() { return &info; }
};
template <class Info>
void addInfo(const std::vector<ad_aug> &x, const Info &info) {
  global::Complete<InfoOp<Info> >(x.size(), info)(x);
}
template <class Info>
void addInfo(const std::vector<double> &x, const Info &info) {}

struct SumOp : global::DynamicOperator<-1, 1> {
  static const bool is_linear = true;
  static const bool have_input_size_output_size = true;
  static const bool add_forward_replay_copy = true;
  size_t n;
  Index input_size() const;
  Index output_size() const;
  SumOp(size_t n);
  template <class Type>
  void forward(ForwardArgs<Type> &args) {
    args.y(0) = 0;
    for (size_t i = 0; i < n; i++) {
      args.y(0) += args.x(i);
    }
  }
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    for (size_t i = 0; i < n; i++) {
      args.dx(i) += args.dy(0);
    }
  }
  const char *op_name();
};
template <class T>
T sum(const std::vector<T> &x) {
  return global::Complete<SumOp>(x.size())(x)[0];
}

ad_plain logspace_sum(const std::vector<ad_plain> &x);
struct LogSpaceSumOp : global::DynamicOperator<-1, 1> {
  size_t n;
  static const bool have_input_size_output_size = true;
  Index input_size() const;
  Index output_size() const;
  LogSpaceSumOp(size_t n);
  void forward(ForwardArgs<Scalar> &args);
  void forward(ForwardArgs<Replay> &args);
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    for (size_t i = 0; i < n; i++) {
      args.dx(i) += exp(args.x(i) - args.y(0)) * args.dy(0);
    }
  }
  const char *op_name();
};
ad_plain logspace_sum(const std::vector<ad_plain> &x);
template <class T>
T logspace_sum(const std::vector<T> &x_) {
  std::vector<ad_plain> x(x_.begin(), x_.end());
  return logspace_sum(x);
}

ad_plain logspace_sum_stride(const std::vector<ad_plain> &x,
                             const std::vector<Index> &stride, size_t n);
struct LogSpaceSumStrideOp : global::DynamicOperator<-1, 1> {
  std::vector<Index> stride;
  size_t n;
  static const bool have_input_size_output_size = true;

  Index number_of_terms() const;
  template <class Type>
  Type &entry(Type **px, size_t i, size_t j) const {
    return px[j][0 + i * stride[j]];
  }
  template <class Type>
  Type rowsum(Type **px, size_t i) const {
    size_t m = stride.size();
    Type s = (Scalar)(0);
    for (size_t j = 0; j < m; j++) {
      s += entry(px, i, j);
    }
    return s;
  }
  Index input_size() const;
  Index output_size() const;
  LogSpaceSumStrideOp(std::vector<Index> stride, size_t n);
  void forward(ForwardArgs<Scalar> &args);
  void forward(ForwardArgs<Replay> &args);
  template <class Type>
  void reverse(ReverseArgs<Type> &args) {
    size_t m = stride.size();
    std::vector<Type *> wrk1(m);
    std::vector<Type *> wrk2(m);
    Type **px = &(wrk1[0]);
    Type **pdx = &(wrk2[0]);
    for (size_t i = 0; i < m; i++) {
      px[i] = args.x_ptr(i);
      pdx[i] = args.dx_ptr(i);
    }
    for (size_t i = 0; i < n; i++) {
      Type s = rowsum(px, i);
      Type tmp = exp(s - args.y(0)) * args.dy(0);
      for (size_t j = 0; j < m; j++) {
        entry(pdx, i, j) += tmp;
      }
    }
  }
  void dependencies(Args<> &args, Dependencies &dep) const;
  static const bool have_dependencies = true;
  static const bool implicit_dependencies = true;
  static const bool allow_remap = false;
  const char *op_name();

  void forward(ForwardArgs<Writer> &args);
  void reverse(ReverseArgs<Writer> &args);
};
ad_plain logspace_sum_stride(const std::vector<ad_plain> &x,
                             const std::vector<Index> &stride, size_t n);
template <class T>
T logspace_sum_stride(const std::vector<T> &x_,
                      const std::vector<Index> &stride, size_t n) {
  std::vector<ad_plain> x(x_.begin(), x_.end());
  return logspace_sum_stride(x, stride, n);
}
}  // namespace TMBad
#endif  // HAVE_GLOBAL_HPP