diff --git a/external/apktool/apktool b/external/apktool/apktool
new file mode 100755
index 0000000000..e99086b911
--- /dev/null
+++ b/external/apktool/apktool
@@ -0,0 +1,77 @@
+#!/bin/bash
+#
+# Copyright (C) 2007 The Android Open Source Project
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# This script is a wrapper for smali.jar, so you can simply call "smali",
+# instead of java -jar smali.jar. It is heavily based on the "dx" script
+# from the Android SDK
+
+# Set up prog to be the path of this script, including following symlinks,
+# and set up progdir to be the fully-qualified pathname of its directory.
+prog="$0"
+while [ -h "${prog}" ]; do
+ newProg=`/bin/ls -ld "${prog}"`
+
+ newProg=`expr "${newProg}" : ".* -> \(.*\)$"`
+ if expr "x${newProg}" : 'x/' >/dev/null; then
+ prog="${newProg}"
+ else
+ progdir=`dirname "${prog}"`
+ prog="${progdir}/${newProg}"
+ fi
+done
+oldwd=`pwd`
+progdir=`dirname "${prog}"`
+cd "${progdir}"
+progdir=`pwd`
+prog="${progdir}"/`basename "${prog}"`
+cd "${oldwd}"
+
+jarfile=apktool.jar
+libdir="$progdir"
+if [ ! -r "$libdir/$jarfile" ]
+then
+ echo `basename "$prog"`": can't find $jarfile"
+ exit 1
+fi
+
+javaOpts=""
+
+# If you want DX to have more memory when executing, uncomment the following
+# line and adjust the value accordingly. Use "java -X" for a list of options
+# you can pass here.
+#
+javaOpts="-Xmx512M -Dfile.encoding=utf-8"
+
+# Alternatively, this will extract any parameter "-Jxxx" from the command line
+# and pass them to Java (instead of to dx). This makes it possible for you to
+# add a command-line parameter such as "-JXmx256M" in your ant scripts, for
+# example.
+while expr "x$1" : 'x-J' >/dev/null; do
+ opt=`expr "$1" : '-J\(.*\)'`
+ javaOpts="${javaOpts} -${opt}"
+ shift
+done
+
+if [ "$OSTYPE" = "cygwin" ] ; then
+ jarpath=`cygpath -w "$libdir/$jarfile"`
+else
+ jarpath="$libdir/$jarfile"
+fi
+
+# add current location to path for aapt
+PATH=$PATH:`pwd`;
+export PATH;
+exec java $javaOpts -jar "$jarpath" "$@"
\ No newline at end of file
diff --git a/external/apktool/apktool.jar b/external/apktool/apktool.jar
new file mode 120000
index 0000000000..566a437afc
--- /dev/null
+++ b/external/apktool/apktool.jar
@@ -0,0 +1 @@
+apktool_2.3.0.jar
\ No newline at end of file
diff --git a/external/apktool/apktool_2.3.0.jar b/external/apktool/apktool_2.3.0.jar
new file mode 100644
index 0000000000..bb9003ec40
Binary files /dev/null and b/external/apktool/apktool_2.3.0.jar differ
diff --git a/external/azcopy/azcopy b/external/azcopy/azcopy
new file mode 100755
index 0000000000..d78ea5e2fe
--- /dev/null
+++ b/external/azcopy/azcopy
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b7d08710ae4c0bbdc73184453b4d40da979fa3b0abbdb4e974980f9250730127
+size 18914367
diff --git a/external/bin/MP4Box b/external/bin/MP4Box
new file mode 120000
index 0000000000..014bd86c6a
--- /dev/null
+++ b/external/bin/MP4Box
@@ -0,0 +1 @@
+../gpac/bin/MP4Box
\ No newline at end of file
diff --git a/external/bin/apktool b/external/bin/apktool
new file mode 120000
index 0000000000..88f53d4ad8
--- /dev/null
+++ b/external/bin/apktool
@@ -0,0 +1 @@
+../apktool/apktool
\ No newline at end of file
diff --git a/external/bin/azcopy b/external/bin/azcopy
new file mode 120000
index 0000000000..b0c92398bf
--- /dev/null
+++ b/external/bin/azcopy
@@ -0,0 +1 @@
+../azcopy/azcopy
\ No newline at end of file
diff --git a/external/bin/capnp b/external/bin/capnp
new file mode 120000
index 0000000000..4fa0c5d67d
--- /dev/null
+++ b/external/bin/capnp
@@ -0,0 +1 @@
+../capnp/bin/capnp
\ No newline at end of file
diff --git a/external/bin/capnpc b/external/bin/capnpc
new file mode 120000
index 0000000000..833204686c
--- /dev/null
+++ b/external/bin/capnpc
@@ -0,0 +1 @@
+../capnp/bin/capnpc
\ No newline at end of file
diff --git a/external/bin/capnpc-c b/external/bin/capnpc-c
new file mode 120000
index 0000000000..5f36acd2aa
--- /dev/null
+++ b/external/bin/capnpc-c
@@ -0,0 +1 @@
+../capnp/bin/capnpc-c
\ No newline at end of file
diff --git a/external/bin/capnpc-c++ b/external/bin/capnpc-c++
new file mode 120000
index 0000000000..c63802a798
--- /dev/null
+++ b/external/bin/capnpc-c++
@@ -0,0 +1 @@
+../capnp/bin/capnpc-c++
\ No newline at end of file
diff --git a/external/bin/capnpc-capnp b/external/bin/capnpc-capnp
new file mode 120000
index 0000000000..6ce8504c1f
--- /dev/null
+++ b/external/bin/capnpc-capnp
@@ -0,0 +1 @@
+../capnp/bin/capnpc-capnp
\ No newline at end of file
diff --git a/external/bin/capnpc-java b/external/bin/capnpc-java
new file mode 120000
index 0000000000..62e7133c48
--- /dev/null
+++ b/external/bin/capnpc-java
@@ -0,0 +1 @@
+../capnp/bin/capnpc-java
\ No newline at end of file
diff --git a/external/bin/ffmpeg b/external/bin/ffmpeg
new file mode 120000
index 0000000000..68be8c0682
--- /dev/null
+++ b/external/bin/ffmpeg
@@ -0,0 +1 @@
+../ffmpeg/bin/ffmpeg
\ No newline at end of file
diff --git a/external/bin/ffprobe b/external/bin/ffprobe
new file mode 120000
index 0000000000..3948a822b4
--- /dev/null
+++ b/external/bin/ffprobe
@@ -0,0 +1 @@
+../ffmpeg/bin/ffprobe
\ No newline at end of file
diff --git a/external/bin/flame b/external/bin/flame
new file mode 120000
index 0000000000..890f902508
--- /dev/null
+++ b/external/bin/flame
@@ -0,0 +1 @@
+../pyflame/flame.sh
\ No newline at end of file
diff --git a/external/bin/img2simg b/external/bin/img2simg
new file mode 120000
index 0000000000..14f3994a9f
--- /dev/null
+++ b/external/bin/img2simg
@@ -0,0 +1 @@
+../simg2img/img2simg
\ No newline at end of file
diff --git a/external/bin/ipfs b/external/bin/ipfs
new file mode 120000
index 0000000000..223f8fb14c
--- /dev/null
+++ b/external/bin/ipfs
@@ -0,0 +1 @@
+../ipfs/ipfs
\ No newline at end of file
diff --git a/external/bin/simg2img b/external/bin/simg2img
new file mode 120000
index 0000000000..db45bc1221
--- /dev/null
+++ b/external/bin/simg2img
@@ -0,0 +1 @@
+../simg2img/simg2img
\ No newline at end of file
diff --git a/external/capnp/bin/capnp b/external/capnp/bin/capnp
new file mode 100755
index 0000000000..5d12ba79d6
Binary files /dev/null and b/external/capnp/bin/capnp differ
diff --git a/external/capnp/bin/capnpc b/external/capnp/bin/capnpc
new file mode 120000
index 0000000000..5668473f09
--- /dev/null
+++ b/external/capnp/bin/capnpc
@@ -0,0 +1 @@
+capnp
\ No newline at end of file
diff --git a/external/capnp/bin/capnpc-c b/external/capnp/bin/capnpc-c
new file mode 100755
index 0000000000..657d51de88
Binary files /dev/null and b/external/capnp/bin/capnpc-c differ
diff --git a/external/capnp/bin/capnpc-c++ b/external/capnp/bin/capnpc-c++
new file mode 100755
index 0000000000..bb18064839
Binary files /dev/null and b/external/capnp/bin/capnpc-c++ differ
diff --git a/external/capnp/bin/capnpc-capnp b/external/capnp/bin/capnpc-capnp
new file mode 100755
index 0000000000..a645fe78f3
Binary files /dev/null and b/external/capnp/bin/capnpc-capnp differ
diff --git a/external/capnp/bin/capnpc-java b/external/capnp/bin/capnpc-java
new file mode 100755
index 0000000000..55af5b5763
Binary files /dev/null and b/external/capnp/bin/capnpc-java differ
diff --git a/external/capnp/build.sh b/external/capnp/build.sh
new file mode 100755
index 0000000000..deb7723ef6
--- /dev/null
+++ b/external/capnp/build.sh
@@ -0,0 +1,47 @@
+set -e
+echo "Installing capnp"
+
+ONE=${HOME}/one
+VERSION=0.6.1
+wget https://capnproto.org/capnproto-c++-${VERSION}.tar.gz
+tar xvf capnproto-c++-${VERSION}.tar.gz
+cd capnproto-c++-${VERSION}
+CXXFLAGS="-fPIC" ./configure --prefix=${ONE}/external/capnp
+make -j9
+
+# manually build binaries statically
+g++ -std=gnu++11 -I./src -I./src -DKJ_HEADER_WARNINGS -DCAPNP_HEADER_WARNINGS -DCAPNP_INCLUDE_DIR=\"/usr/local/include\" -pthread -O2 -DNDEBUG -pthread -pthread -o .libs/capnp src/capnp/compiler/module-loader.o src/capnp/compiler/capnp.o ./.libs/libcapnpc.a ./.libs/libcapnp.a ./.libs/libkj.a -lpthread -pthread
+
+g++ -std=gnu++11 -I./src -I./src -DKJ_HEADER_WARNINGS -DCAPNP_HEADER_WARNINGS -DCAPNP_INCLUDE_DIR=\"/usr/local/include\" -pthread -O2 -DNDEBUG -pthread -pthread -o .libs/capnpc-c++ src/capnp/compiler/capnpc-c++.o ./.libs/libcapnp.a ./.libs/libkj.a -lpthread -pthread
+
+g++ -std=gnu++11 -I./src -I./src -DKJ_HEADER_WARNINGS -DCAPNP_HEADER_WARNINGS -DCAPNP_INCLUDE_DIR=\"/usr/local/include\" -pthread -O2 -DNDEBUG -pthread -pthread -o .libs/capnpc-capnp src/capnp/compiler/capnpc-capnp.o ./.libs/libcapnp.a ./.libs/libkj.a -lpthread -pthread
+
+
+make -j4 install
+
+# --------
+echo "Installing c-capnp"
+
+git clone git@github.com:commaai/c-capnproto.git
+cd c-capnproto
+git submodule update --init --recursive
+autoreconf -f -i -s
+CFLAGS="-fPIC" ./configure --prefix=${ONE}/external/capnp
+make -j4
+
+# manually build binaries statically
+gcc -fPIC -o .libs/capnpc-c compiler/capnpc-c.o compiler/schema.capnp.o compiler/str.o ./.libs/libcapnp_c.a -Wl,-rpath -Wl,${ONE}/external/capnp/lib
+
+make install
+
+# --------
+echo "Installing java-capnp"
+
+git clone https://github.com/dwrensha/capnproto-java.git
+cd capnproto-java
+git reset --hard 2c43bd712fb218da0eabdf241a750b9c05903e8e
+g++ compiler/src/main/cpp/capnpc-java.c++ -std=c++11 -pthread -I${ONE}/external/capnp/include -L${ONE}/external/capnp/lib -l:libcapnp.a -l:libkj.a -pthread -lpthread -o capnpc-java
+cp capnpc-java ${ONE}/external/capnp/bin/
+
+rm -rf capnproto-c++-${VERSION}.tar.gz
+rm -rf capnproto-c++-${VERSION}
diff --git a/external/cppad/build.txt b/external/cppad/build.txt
new file mode 100644
index 0000000000..5e9181b4ad
--- /dev/null
+++ b/external/cppad/build.txt
@@ -0,0 +1,7 @@
+wget 'http://www.coin-or.org/download/source/CppAD/cppad-20170816.epl.tgz'
+tar xvf cppad-20170816.epl.tgz
+cd cppad-20170816/
+mkdir build
+cd build
+cmake -D cppad_prefix="$HOME/one/external/cppad" ..
+make install
diff --git a/external/cppad/include/cppad/base_require.hpp b/external/cppad/include/cppad/base_require.hpp
new file mode 100644
index 0000000000..1c02a96995
--- /dev/null
+++ b/external/cppad/include/cppad/base_require.hpp
@@ -0,0 +1,185 @@
+// $Id: base_require.hpp 3845 2016-11-19 01:50:47Z bradbell $
+# ifndef CPPAD_BASE_REQUIRE_HPP
+# define CPPAD_BASE_REQUIRE_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+$begin base_require$$
+$spell
+ azmul
+ ostream
+ alloc
+ eps
+ std
+ Lt
+ Le
+ Eq
+ Ge
+ Gt
+ cppad.hpp
+ namespace
+ optimizations
+ bool
+ const
+ CppAD
+ enum
+ Lt
+ Le
+ Eq
+ Ge
+ Gt
+ inline
+ Op
+ std
+ CondExp
+$$
+
+$section AD Requirements for a CppAD Base Type$$
+
+$head Syntax$$
+$code # include $$
+
+$head Purpose$$
+This section lists the requirements for the type
+$icode Base$$ so that the type $codei%AD<%Base%>%$$ can be used.
+
+$head API Warning$$
+Defining a CppAD $icode Base$$ type is an advanced use of CppAD.
+This part of the CppAD API changes with time. The most common change
+is adding more requirements.
+Search for $code base_require$$ in the
+current $cref whats_new$$ section for these changes.
+
+$head Standard Base Types$$
+In the case where $icode Base$$ is
+$code float$$,
+$code double$$,
+$code std::complex$$,
+$code std::complex$$,
+or $codei%AD<%Other%>%$$,
+these requirements are provided by including the file
+$code cppad/cppad.hpp$$.
+
+$head Include Order$$
+If you are linking a non-standard base type to CppAD,
+you must first include the file $code cppad/base_require.hpp$$,
+then provide the specifications below,
+and then include the file $code cppad/cppad.hpp$$.
+
+$head Numeric Type$$
+The type $icode Base$$ must support all the operations for a
+$cref NumericType$$.
+
+$head Output Operator$$
+The type $icode Base$$ must support the syntax
+$codei%
+ %os% << %x%
+%$$
+where $icode os$$ is an $code std::ostream&$$
+and $icode x$$ is a $code const base_alloc&$$.
+For example, see
+$cref/base_alloc/base_alloc.hpp/Output Operator/$$.
+
+$head Integer$$
+The type $icode Base$$ must support the syntax
+$codei%
+ %i% = CppAD::Integer(%x%)
+%$$
+which converts $icode x$$ to an $code int$$.
+The argument $icode x$$ has prototype
+$codei%
+ const %Base%& %x%
+%$$
+and the return value $icode i$$ has prototype
+$codei%
+ int %i%
+%$$
+
+$subhead Suggestion$$
+In many cases, the $icode Base$$ version of the $code Integer$$ function
+can be defined by
+$codei%
+namespace CppAD {
+ inline int Integer(const %Base%& x)
+ { return static_cast(x); }
+}
+%$$
+For example, see
+$cref/base_float/base_float.hpp/Integer/$$ and
+$cref/base_alloc/base_alloc.hpp/Integer/$$.
+
+$head Absolute Zero, azmul$$
+The type $icode Base$$ must support the syntax
+$codei%
+ %z% = azmul(%x%, %y%)
+%$$
+see; $cref azmul$$.
+The following preprocessor macro invocation suffices
+(for most $icode Base$$ types):
+$codei%
+namespace CppAD {
+ CPPAD_AZMUL(%Base%)
+}
+%$$
+where the macro is defined by
+$srccode%cpp% */
+# define CPPAD_AZMUL(Base) \
+ inline Base azmul(const Base& x, const Base& y) \
+ { Base zero(0.0); \
+ if( x == zero ) \
+ return zero; \
+ return x * y; \
+ }
+/* %$$
+
+$childtable%
+ omh/base_require/base_member.omh%
+ cppad/core/base_cond_exp.hpp%
+ omh/base_require/base_identical.omh%
+ omh/base_require/base_ordered.omh%
+ cppad/core/base_std_math.hpp%
+ cppad/core/base_limits.hpp%
+ cppad/core/base_to_string.hpp%
+ cppad/core/base_hash.hpp%
+ omh/base_require/base_example.omh
+%$$
+
+$end
+*/
+
+// definitions that must come before base implementations
+# include
+# include
+# include
+# include
+
+// grouping documentation by feature
+# include
+# include
+# include
+# include
+# include
+
+// must define template class numeric_limits before the base cases
+# include
+# include // deprecated
+
+// base cases that come with CppAD
+# include
+# include
+# include
+
+// deprecated base type
+# include
+
+# endif
diff --git a/external/cppad/include/cppad/configure.hpp b/external/cppad/include/cppad/configure.hpp
new file mode 100644
index 0000000000..26cf15a63a
--- /dev/null
+++ b/external/cppad/include/cppad/configure.hpp
@@ -0,0 +1,216 @@
+# ifndef CPPAD_CONFIGURE_HPP
+# define CPPAD_CONFIGURE_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*!
+ \file configure.hpp
+Replacement for config.h so that all preprocessor symbols begin with CPPAD_
+*/
+
+/*!
+\def CPPAD_COMPILER_IS_GNUCXX
+is the compiler a variant of g++
+*/
+# define CPPAD_COMPILER_IS_GNUCXX 1
+
+/*!
+\def CPPAD_DISABLE_SOME_MICROSOFT_COMPILER_WARNINGS
+This macro is only used to document the pragmas that disables the
+follow warnings:
+
+\li C4100
+unreferenced formal parameter.
+\li C4127
+conditional expression is constant.
+*/
+# define CPPAD_DISABLE_SOME_MICROSOFT_COMPILER_WARNINGS 1
+# if _MSC_VER
+# pragma warning( disable : 4100 )
+# pragma warning( disable : 4127 )
+# endif
+# undef CPPAD_DISABLE_SOME_MICROSOFT_COMPILER_WARNINGS
+
+/*!
+\def CPPAD_USE_CPLUSPLUS_2011
+Should CppAD use C++11 features. This will be true if the current
+compiler flags request C++11 features and the install procedure
+determined that all the necessary features are avaiable.
+*/
+# if _MSC_VER
+# define CPPAD_USE_CPLUSPLUS_2011 0
+# else //
+# if __cplusplus >= 201100
+# define CPPAD_USE_CPLUSPLUS_2011 0
+# else //
+# define CPPAD_USE_CPLUSPLUS_2011 0
+# endif //
+# endif //
+
+/*!
+\def CPPAD_PACKAGE_STRING
+cppad-yyyymmdd as a C string where yyyy is year, mm is month, and dd is day.
+*/
+# define CPPAD_PACKAGE_STRING "cppad-20170816"
+
+/*!
+def CPPAD_HAS_ADOLC
+Was a adolc_prefix specified on the cmake command line.
+*/
+# define CPPAD_HAS_ADOLC 0
+
+/*!
+def CPPAD_HAS_COLPACK
+Was a colpack_prefix specified on the cmake command line.
+*/
+# define CPPAD_HAS_COLPACK 0
+
+/*!
+def CPPAD_HAS_EIGEN
+Was a eigen_prefix specified on the cmake command line.
+*/
+# define CPPAD_HAS_EIGEN 0
+
+/*!
+def CPPAD_HAS_IPOPT
+Was a ipopt_prefix specified on the cmake command line.
+*/
+# define CPPAD_HAS_IPOPT 0
+
+/*!
+\def CPPAD_DEPRECATED
+This symbol is not currently being used.
+*/
+# define CPPAD_DEPRECATED 0
+
+/*!
+\def CPPAD_BOOSTVECTOR
+If this symbol is one, and _MSC_VER is not defined,
+we are using boost vector for CPPAD_TESTVECTOR.
+It this symbol is zero,
+we are not using boost vector for CPPAD_TESTVECTOR.
+*/
+# define CPPAD_BOOSTVECTOR 0
+
+/*!
+\def CPPAD_CPPADVECTOR
+If this symbol is one,
+we are using CppAD vector for CPPAD_TESTVECTOR.
+It this symbol is zero,
+we are not using CppAD vector for CPPAD_TESTVECTOR.
+*/
+# define CPPAD_CPPADVECTOR 1
+
+/*!
+\def CPPAD_STDVECTOR
+If this symbol is one,
+we are using standard vector for CPPAD_TESTVECTOR.
+It this symbol is zero,
+we are not using standard vector for CPPAD_TESTVECTOR.
+*/
+# define CPPAD_STDVECTOR 0
+
+/*!
+\def CPPAD_EIGENVECTOR
+If this symbol is one,
+we are using Eigen vector for CPPAD_TESTVECTOR.
+If this symbol is zero,
+we are not using Eigen vector for CPPAD_TESTVECTOR.
+*/
+# define CPPAD_EIGENVECTOR 0
+
+/*!
+\def CPPAD_HAS_GETTIMEOFDAY
+If this symbol is one, and _MSC_VER is not defined,
+this system supports the gettimeofday funcgtion.
+Otherwise, this smybol should be zero.
+*/
+# define CPPAD_HAS_GETTIMEOFDAY 1
+
+/*!
+\def CPPAD_SIZE_T_NOT_UNSIGNED_INT
+If this symbol is zero, the type size_t is the same as the type unsigned int,
+otherwise this symbol is one.
+*/
+# define CPPAD_SIZE_T_NOT_UNSIGNED_INT 1
+
+/*!
+\def CPPAD_TAPE_ADDR_TYPE
+Is the type used to store address on the tape. If not size_t, then
+sizeof(CPPAD_TAPE_ADDR_TYPE) <= sizeof( size_t )
+to conserve memory.
+This type must support \c std::numeric_limits,
+the \c <= operator,
+and conversion to \c size_t.
+Make sure that the type chosen returns true for is_pod
+in pod_vector.hpp.
+This type is later defined as \c addr_t in the CppAD namespace.
+*/
+# define CPPAD_TAPE_ADDR_TYPE unsigned int
+
+/*!
+\def CPPAD_TAPE_ID_TYPE
+Is the type used to store tape identifiers. If not size_t, then
+sizeof(CPPAD_TAPE_ID_TYPE) <= sizeof( size_t )
+to conserve memory.
+This type must support \c std::numeric_limits,
+the \c <= operator,
+and conversion to \c size_t.
+Make sure that the type chosen returns true for is_pod
+in pod_vector.hpp.
+This type is later defined as \c tape_id_t in the CppAD namespace.
+*/
+# define CPPAD_TAPE_ID_TYPE unsigned int
+
+/*!
+\def CPPAD_MAX_NUM_THREADS
+Specifies the maximum number of threads that CppAD can support
+(must be greater than or equal four).
+
+The user may define CPPAD_MAX_NUM_THREADS before including any of the CppAD
+header files. If it is not yet defined,
+*/
+# ifndef CPPAD_MAX_NUM_THREADS
+# define CPPAD_MAX_NUM_THREADS 48
+# endif
+
+/*!
+\def CPPAD_HAS_MKSTEMP
+It true, mkstemp works in C++ on this system.
+*/
+# define CPPAD_HAS_MKSTEMP 1
+
+/*!
+\def CPPAD_HAS_TMPNAM_S
+It true, tmpnam_s works in C++ on this system.
+*/
+# define CPPAD_HAS_TMPNAM_S 0
+
+// ---------------------------------------------------------------------------
+// defines that only depend on values above
+// ---------------------------------------------------------------------------
+/*!
+\def CPPAD_NULL
+This preprocessor symbol is used for a null pointer.
+
+If it is not yet defined,
+it is defined when cppad/core/define.hpp is included.
+*/
+# ifndef CPPAD_NULL
+# if CPPAD_USE_CPLUSPLUS_2011
+# define CPPAD_NULL nullptr
+# else
+# define CPPAD_NULL 0
+# endif
+# endif
+
+# endif
diff --git a/external/cppad/include/cppad/core/abort_recording.hpp b/external/cppad/include/cppad/core/abort_recording.hpp
new file mode 100644
index 0000000000..19938b0c44
--- /dev/null
+++ b/external/cppad/include/cppad/core/abort_recording.hpp
@@ -0,0 +1,60 @@
+# ifndef CPPAD_CORE_ABORT_RECORDING_HPP
+# define CPPAD_CORE_ABORT_RECORDING_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+/*
+$begin abort_recording$$
+$spell
+$$
+
+$section Abort Recording of an Operation Sequence$$
+$mindex tape$$
+
+
+$head Syntax$$
+$codei%AD<%Base%>::abort_recording()%$$
+
+$head Purpose$$
+Sometimes it is necessary to abort the recording of an operation sequence
+that started with a call of the form
+$codei%
+ Independent(%x%)
+%$$
+If such a recording is currently in progress,
+$code abort_recording$$ will stop the recording and delete the
+corresponding information.
+Otherwise, $code abort_recording$$ has no effect.
+
+$children%
+ example/general/abort_recording.cpp
+%$$
+$head Example$$
+The file
+$cref abort_recording.cpp$$
+contains an example and test of this operation.
+It returns true if it succeeds and false otherwise.
+
+$end
+----------------------------------------------------------------------------
+*/
+
+
+namespace CppAD {
+ template
+ void AD::abort_recording(void)
+ { local::ADTape* tape = AD::tape_ptr();
+ if( tape != CPPAD_NULL )
+ AD::tape_manage(tape_manage_delete);
+ }
+}
+
+# endif
diff --git a/external/cppad/include/cppad/core/abs.hpp b/external/cppad/include/cppad/core/abs.hpp
new file mode 100644
index 0000000000..11c2e17c20
--- /dev/null
+++ b/external/cppad/include/cppad/core/abs.hpp
@@ -0,0 +1,114 @@
+# ifndef CPPAD_CORE_ABS_HPP
+# define CPPAD_CORE_ABS_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+-------------------------------------------------------------------------------
+$begin abs$$
+$spell
+ fabs
+ Vec
+ std
+ faq
+ Taylor
+ Cpp
+ namespace
+ const
+ abs
+$$
+
+$section AD Absolute Value Functions: abs, fabs$$
+
+$head Syntax$$
+$icode%y% = abs(%x%)
+%$$
+$icode%y% = fabs(%x%)%$$
+
+$head x, y$$
+See the $cref/possible types/unary_standard_math/Possible Types/$$
+for a unary standard math function.
+
+$head Atomic$$
+In the case where $icode x$$ is an AD type,
+this is an $cref/atomic operation/glossary/Operation/Atomic/$$.
+
+$head Complex Types$$
+The functions $code abs$$ and $icode fabs$$
+are not defined for the base types
+$code std::complex$$ or $code std::complex$$
+because the complex $code abs$$ function is not complex differentiable
+(see $cref/complex types faq/Faq/Complex Types/$$).
+
+$head Derivative$$
+CppAD defines the derivative of the $code abs$$ function is
+the $cref sign$$ function; i.e.,
+$latex \[
+{\rm abs}^{(1)} ( x ) = {\rm sign} (x ) =
+\left\{ \begin{array}{rl}
+ +1 & {\rm if} \; x > 0 \\
+ 0 & {\rm if} \; x = 0 \\
+ -1 & {\rm if} \; x < 0
+\end{array} \right.
+\] $$
+The result for $icode%x% == 0%$$ used to be a directional derivative.
+
+$head Example$$
+$children%
+ example/general/fabs.cpp
+%$$
+The file
+$cref fabs.cpp$$
+contains an example and test of this function.
+It returns true if it succeeds and false otherwise.
+
+$end
+-------------------------------------------------------------------------------
+*/
+
+// BEGIN CppAD namespace
+namespace CppAD {
+
+template
+AD AD::abs_me (void) const
+{
+ AD result;
+ result.value_ = abs(value_);
+ CPPAD_ASSERT_UNKNOWN( Parameter(result) );
+
+ if( Variable(*this) )
+ { // add this operation to the tape
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::AbsOp) == 1 );
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::AbsOp) == 1 );
+ local::ADTape *tape = tape_this();
+
+ // corresponding operand address
+ tape->Rec_.PutArg(taddr_);
+ // put operator in the tape
+ result.taddr_ = tape->Rec_.PutOp(local::AbsOp);
+ // make result a variable
+ result.tape_id_ = tape->id_;
+ }
+ return result;
+}
+
+template
+inline AD abs(const AD &x)
+{ return x.abs_me(); }
+
+template
+inline AD abs(const VecAD_reference &x)
+{ return x.ADBase().abs_me(); }
+
+} // END CppAD namespace
+
+# endif
diff --git a/external/cppad/include/cppad/core/abs_normal_fun.hpp b/external/cppad/include/cppad/core/abs_normal_fun.hpp
new file mode 100644
index 0000000000..f974046ca8
--- /dev/null
+++ b/external/cppad/include/cppad/core/abs_normal_fun.hpp
@@ -0,0 +1,867 @@
+# ifndef CPPAD_CORE_ABS_NORMAL_FUN_HPP
+# define CPPAD_CORE_ABS_NORMAL_FUN_HPP
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+/*
+$begin abs_normal_fun$$
+$spell
+ const
+$$
+
+
+$section Create An Abs-normal Representation of a Function$$
+
+$head Syntax$$
+$icode%f%.abs_normal_fun(%g%, %a%)%$$
+
+$head f$$
+The object $icode f$$ has prototype
+$codei%
+ ADFun<%Base%>& %f%
+%$$
+It represents a function $latex f : \B{R}^n \rightarrow \B{R}^m$$.
+We assume that the only non-smooth terms in the representation are
+absolute value functions and use $latex s \in \B{Z}_+$$
+to represent the number of these terms.
+It is effectively $code const$$, except that some internal state
+that is not relevant to the user; see
+$cref/const ADFun/wish_list/const ADFun/$$.
+
+$subhead n$$
+We use $icode n$$ to denote the dimension of the domain space for $icode f$$.
+
+$subhead m$$
+We use $icode m$$ to denote the dimension of the range space for $icode f$$.
+
+$subhead s$$
+We use $icode s$$ to denote the number of absolute value terms in $icode f$$.
+
+
+$head a$$
+The object $icode a$$ has prototype
+$codei%
+ ADFun<%Base%> %a%
+%$$
+The initial function representation in $icode a$$ is lost.
+Upon return it represents the result of the absolute terms
+$latex a : \B{R}^n \rightarrow \B{R}^s$$; see $latex a(x)$$ defined below.
+Note that $icode a$$ is constructed by copying $icode f$$
+and then changing the dependent variables. There may
+be many calculations in this representation that are not necessary
+and can be removed using
+$codei%
+ %a%.optimize()
+%$$
+This optimization is not done automatically by $code abs_normal_fun$$
+because it may take a significant amount of time.
+
+$subhead zeta$$
+Let $latex \zeta_0 ( x )$$
+denote the argument for the first absolute value term in $latex f(x)$$,
+$latex \zeta_1 ( x , |\zeta_0 (x)| )$$ for the second term, and so on.
+
+$subhead a(x)$$
+For $latex i = 0 , \ldots , {s-1}$$ define
+$latex \[
+a_i (x)
+=
+| \zeta_i ( x , a_0 (x) , \ldots , a_{i-1} (x ) ) |
+\] $$
+This defines $latex a : \B{R}^n \rightarrow \B{R}^s$$.
+
+$head g$$
+The object $icode g$$ has prototype
+$codei%
+ ADFun<%Base%> %g%
+%$$
+The initial function representation in $icode g$$ is lost.
+Upon return it represents the smooth function
+$latex g : \B{R}^{n + s} \rightarrow \B{R}^{m + s}$$ is defined by
+$latex \[
+g( x , u )
+=
+\left[ \begin{array}{c} y(x, u) \\ z(x, u) \end{array} \right]
+\] $$
+were $latex y(x, u)$$ and $latex z(x, u)$$ are defined below.
+
+$subhead z(x, u)$$
+Define the smooth function
+$latex z : \B{R}^{n + s} \rightarrow \B{R}^s$$ by
+$latex \[
+z_i ( x , u ) = \zeta_i ( x , u_0 , \ldots , u_{i-1} )
+\] $$
+Note that the partial of $latex z_i$$ with respect to $latex u_j$$ is zero
+for $latex j \geq i$$.
+
+$subhead y(x, u)$$
+There is a smooth function
+$latex y : \B{R}^{n + s} \rightarrow \B{R}^m$$
+such that $latex y( x , u ) = f(x)$$ whenever $latex u = a(x)$$.
+
+$head Affine Approximation$$
+We define the affine approximations
+$latex \[
+\begin{array}{rcl}
+y[ \hat{x} ]( x , u )
+& = &
+y ( \hat{x}, a( \hat{x} ) )
+ + \partial_x y ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} )
+ + \partial_u y ( \hat{x}, a( \hat{x} ) ) ( u - a( \hat{x} ) )
+\\
+z[ \hat{x} ]( x , u )
+& = &
+z ( \hat{x}, a( \hat{x} ) )
+ + \partial_x z ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} )
+ + \partial_u z ( \hat{x}, a( \hat{x} ) ) ( u - a( \hat{x} ) )
+\end{array}
+\] $$
+It follows that
+$latex \[
+\begin{array}{rcl}
+y( x , u )
+& = &
+y[ \hat{x} ]( x , u ) + o ( x - \hat{x}, u - a( \hat{x} ) )
+\\
+z( x , u )
+& = &
+z[ \hat{x} ]( x , u ) + o ( x - \hat{x}, u - a( \hat{x} ) )
+\end{array}
+\] $$
+
+$head Abs-normal Approximation$$
+
+$subhead Approximating a(x)$$
+The function $latex a(x)$$ is not smooth, but it is equal to
+$latex | z(x, u) |$$ when $latex u = a(x)$$.
+Furthermore
+$latex \[
+z[ \hat{x} ]( x , u )
+=
+z ( \hat{x}, a( \hat{x} ) )
+ + \partial_x z ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} )
+ + \partial_u z ( \hat{x}, a( \hat{x} ) ) ( u - a( \hat{x} ) )
+\] $$
+The partial of $latex z_i$$ with respect to $latex u_j$$ is zero
+for $latex j \geq i$$. It follows that
+$latex \[
+z_i [ \hat{x} ]( x , u )
+=
+z_i ( \hat{x}, a( \hat{x} ) )
+ + \partial_x z_i ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} )
+ + \sum_{j < i} \partial_{u(j)}
+ z_i ( \hat{x}, a( \hat{x} ) ) ( u_j - a_j ( \hat{x} ) )
+\] $$
+Considering the case $latex i = 0$$ we define
+$latex \[
+a_0 [ \hat{x} ]( x )
+=
+| z_0 [ \hat{x} ]( x , u ) |
+=
+\left|
+ z_0 ( \hat{x}, a( \hat{x} ) )
+ + \partial_x z_0 ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} )
+\right|
+\] $$
+It follows that
+$latex \[
+ a_0 (x) = a_0 [ \hat{x} ]( x ) + o ( x - \hat{x} )
+\] $$
+In general, we define $latex a_i [ \hat{x} ]$$ using
+$latex a_j [ \hat{x} ]$$ for $latex j < i$$ as follows:
+$latex \[
+a_i [ \hat{x} ]( x )
+=
+\left |
+ z_i ( \hat{x}, a( \hat{x} ) )
+ + \partial_x z_i ( \hat{x}, a( \hat{x} ) ) ( x - \hat{x} )
+ + \sum_{j < i} \partial_{u(j)}
+ z_i ( \hat{x}, a( \hat{x} ) )
+ ( a_j [ \hat{x} ] ( x ) - a_j ( \hat{x} ) )
+\right|
+\] $$
+It follows that
+$latex \[
+ a (x) = a[ \hat{x} ]( x ) + o ( x - \hat{x} )
+\] $$
+Note that in the case where $latex z(x, u)$$ and $latex y(x, u)$$ are
+affine,
+$latex \[
+ a[ \hat{x} ]( x ) = a( x )
+\] $$
+
+
+$subhead Approximating f(x)$$
+$latex \[
+f(x)
+=
+y ( x , a(x ) )
+=
+y [ \hat{x} ] ( x , a[ \hat{x} ] ( x ) )
++ o( \Delta x )
+\] $$
+
+$head Correspondence to Literature$$
+Using the notation
+$latex Z = \partial_x z(\hat{x}, \hat{u})$$,
+$latex L = \partial_u z(\hat{x}, \hat{u})$$,
+$latex J = \partial_x y(\hat{x}, \hat{u})$$,
+$latex Y = \partial_u y(\hat{x}, \hat{u})$$,
+the approximation for $latex z$$ and $latex y$$ are
+$latex \[
+\begin{array}{rcl}
+z[ \hat{x} ]( x , u )
+& = &
+z ( \hat{x}, a( \hat{x} ) ) + Z ( x - \hat{x} ) + L ( u - a( \hat{x} ) )
+\\
+y[ \hat{x} ]( x , u )
+& = &
+y ( \hat{x}, a( \hat{x} ) ) + J ( x - \hat{x} ) + Y ( u - a( \hat{x} ) )
+\end{array}
+\] $$
+Moving the terms with $latex \hat{x}$$ together, we have
+$latex \[
+\begin{array}{rcl}
+z[ \hat{x} ]( x , u )
+& = &
+z ( \hat{x}, a( \hat{x} ) ) - Z \hat{x} - L a( \hat{x} ) + Z x + L u
+\\
+y[ \hat{x} ]( x , u )
+& = &
+y ( \hat{x}, a( \hat{x} ) ) - J \hat{x} - Y a( \hat{x} ) + J x + Y u
+\end{array}
+\] $$
+Using the notation
+$latex c = z ( \hat{x}, \hat{u} ) - Z \hat{x} - L \hat{u}$$,
+$latex b = y ( \hat{x}, \hat{u} ) - J \hat{x} - Y \hat{u}$$,
+we have
+$latex \[
+\begin{array}{rcl}
+z[ \hat{x} ]( x , u ) & = & c + Z x + L u
+\\
+y[ \hat{x} ]( x , u ) & = & b + J x + Y u
+\end{array}
+\] $$
+Considering the affine case, where the approximations are exact,
+and choosing $latex u = a(x) = |z(x, u)|$$, we obtain
+$latex \[
+\begin{array}{rcl}
+z( x , a(x ) ) & = & c + Z x + L |z( x , a(x ) )|
+\\
+y( x , a(x ) ) & = & b + J x + Y |z( x , a(x ) )|
+\end{array}
+\] $$
+This is Equation (2) of the
+$cref/reference/abs_normal/Reference/$$.
+
+$children%example/abs_normal/get_started.cpp
+%$$
+$head Example$$
+The file $cref abs_get_started.cpp$$ contains
+an example and test using this operation.
+
+$end
+-------------------------------------------------------------------------------
+*/
+/*!
+file abs_normal_fun.hpp
+Create an abs-normal representation of a function
+*/
+
+namespace CppAD { // BEGIN_CPPAD_NAMESPACE
+/*!
+Create an abs-normal representation of an ADFun object.
+
+\tparam Base
+base type for this abs-normal form and for the function beging represented;
+i.e., f.
+
+\param f
+is the function that this object will represent in abs-normal form.
+This is effectively const except that the play back state play_
+is used.
+*/
+
+# define NOT_YET_COMPILING 0
+
+template
+void ADFun::abs_normal_fun(ADFun& g, ADFun& a)
+{ using namespace local;
+
+ // -----------------------------------------------------------------------
+ // Forward sweep to determine number of absolute value operations in f
+ // -----------------------------------------------------------------------
+ // The argument and result index in f for each absolute value operator
+ CppAD::vector f_abs_arg;
+ CppAD::vector f_abs_res;
+ //
+ OpCode op; // this operator
+ const addr_t* arg = CPPAD_NULL; // arguments for this operator
+ size_t i_op; // index of this operator
+ size_t i_var; // variable index for this operator
+ play_.forward_start(op, arg, i_op, i_var);
+ CPPAD_ASSERT_UNKNOWN( op == BeginOp );
+ //
+ bool more_operators = true;
+ while( more_operators )
+ {
+ // next op
+ play_.forward_next(op, arg, i_op, i_var);
+ switch( op )
+ { // absolute value operator
+ case AbsOp:
+ CPPAD_ASSERT_NARG_NRES(op, 1, 1);
+ f_abs_arg.push_back( arg[0] );
+ f_abs_res.push_back( i_var );
+ break;
+
+ case CSumOp:
+ // CSumOp has a variable number of arguments
+ play_.forward_csum(op, arg, i_op, i_var);
+ break;
+
+ case CSkipOp:
+ // CSkip has a variable number of arguments
+ play_.forward_cskip(op, arg, i_op, i_var);
+ break;
+
+ case EndOp:
+ more_operators = false;
+ break;
+
+ default:
+ break;
+ }
+ }
+ // ------------------------------------------------------------------------
+ // Forward sweep to create new recording
+ // ------------------------------------------------------------------------
+ // recorder for new operation sequence
+ recorder rec;
+ //
+ // number of variables in both operation sequences
+ // (the AbsOp operators are replace by InvOp operators)
+ const size_t num_var = play_.num_var_rec();
+ //
+ // mapping from old variable index to new variable index
+ CPPAD_ASSERT_UNKNOWN(
+ std::numeric_limits::max() >= num_var
+ );
+ CppAD::vector f2g_var(num_var);
+ for(i_var = 0; i_var < num_var; i_var++)
+ f2g_var[i_var] = addr_t( num_var ); // invalid (should not be used)
+ //
+ // record the independent variables in f
+ play_.forward_start(op, arg, i_op, i_var);
+ CPPAD_ASSERT_UNKNOWN( op == BeginOp );
+ more_operators = true;
+ while( more_operators )
+ { switch( op )
+ {
+ // phantom variable
+ case BeginOp:
+ CPPAD_ASSERT_NARG_NRES(op, 1, 1);
+ CPPAD_ASSERT_UNKNOWN( arg[0] == 0 );
+ rec.PutArg(0);
+ f2g_var[i_var] = rec.PutOp(op);
+ break;
+
+ // independent variables
+ case InvOp:
+ CPPAD_ASSERT_NARG_NRES(op, 0, 1);
+ f2g_var[i_var] = rec.PutOp(op);
+ break;
+
+ // end of independent variables
+ default:
+ more_operators = false;
+ break;
+ }
+ if( more_operators )
+ play_.forward_next(op, arg, i_op, i_var);
+ }
+ // add one for the phantom variable
+ CPPAD_ASSERT_UNKNOWN( 1 + Domain() == i_var );
+ //
+ // record the independent variables corresponding AbsOp results
+ size_t index_abs;
+ for(index_abs = 0; index_abs < f_abs_res.size(); index_abs++)
+ f2g_var[ f_abs_res[index_abs] ] = rec.PutOp(InvOp);
+ //
+ // used to hold new argument vector
+ addr_t new_arg[6];
+ //
+ // Parameters in recording of f
+ const Base* f_parameter = play_.GetPar();
+ //
+ // now loop through the rest of the
+ more_operators = true;
+ index_abs = 0;
+ while( more_operators )
+ { addr_t mask; // temporary used in some switch cases
+ switch( op )
+ {
+ // check setting of f_abs_arg and f_abs_res;
+ case AbsOp:
+ CPPAD_ASSERT_NARG_NRES(op, 1, 1);
+ CPPAD_ASSERT_UNKNOWN( f_abs_arg[index_abs] == arg[0] );
+ CPPAD_ASSERT_UNKNOWN( f_abs_res[index_abs] == i_var );
+ CPPAD_ASSERT_UNKNOWN( f2g_var[i_var] > 0 );
+ ++index_abs;
+ break;
+
+ // These operators come at beginning of take and are handled above
+ case InvOp:
+ CPPAD_ASSERT_UNKNOWN(false);
+ break;
+
+ // ---------------------------------------------------------------
+ // Unary operators, argument a parameter, one result
+ case ParOp:
+ CPPAD_ASSERT_NARG_NRES(op, 1, 1);
+ new_arg[0] = rec.PutPar( f_parameter[ arg[0] ] );
+ rec.PutArg( new_arg[0] );
+ f2g_var[i_var] = rec.PutOp(op);
+ break;
+
+ // --------------------------------------------------------------
+ // Unary operators, argument a variable, one result
+ // (excluding the absolute value operator AbsOp)
+ case AcosOp:
+ case AcoshOp:
+ case AsinOp:
+ case AsinhOp:
+ case AtanOp:
+ case AtanhOp:
+ case CosOp:
+ case CoshOp:
+ case ExpOp:
+ case Expm1Op:
+ case LogOp:
+ case Log1pOp:
+ case SignOp:
+ case SinOp:
+ case SinhOp:
+ case SqrtOp:
+ case TanOp:
+ case TanhOp:
+ // some of these operators have an auxillary result; e.g.,
+ // sine and cosine are computed togeather.
+ CPPAD_ASSERT_UNKNOWN( NumArg(op) == 1 );
+ CPPAD_ASSERT_UNKNOWN( NumRes(op) == 1 || NumRes(op) == 2 );
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[0] ] ) < num_var );
+ new_arg[0] = f2g_var[ arg[0] ];
+ rec.PutArg( new_arg[0] );
+ f2g_var[i_var] = rec.PutOp( op );
+ break;
+
+ case ErfOp:
+ CPPAD_ASSERT_NARG_NRES(op, 3, 5);
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[0] ] ) < num_var );
+ // Error function is a special case
+ // second argument is always the parameter 0
+ // third argument is always the parameter 2 / sqrt(pi)
+ rec.PutArg( rec.PutPar( Base(0.0) ) );
+ rec.PutArg( rec.PutPar(
+ Base( 1.0 / std::sqrt( std::atan(1.0) ) )
+ ) );
+ f2g_var[i_var] = rec.PutOp(op);
+ break;
+ // --------------------------------------------------------------
+ // Binary operators, left variable, right parameter, one result
+ case SubvpOp:
+ case DivvpOp:
+ case PowvpOp:
+ case ZmulvpOp:
+ CPPAD_ASSERT_NARG_NRES(op, 2, 1);
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[0] ] ) < num_var );
+ new_arg[0] = f2g_var[ arg[0] ];
+ new_arg[1] = rec.PutPar( f_parameter[ arg[1] ] );
+ rec.PutArg( new_arg[0], new_arg[1] );
+ f2g_var[i_var] = rec.PutOp(op);
+ break;
+ // ---------------------------------------------------
+ // Binary operators, left index, right variable, one result
+ case DisOp:
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
+ new_arg[0] = arg[0];
+ new_arg[1] = f2g_var[ arg[1] ];
+ rec.PutArg( new_arg[0], new_arg[1] );
+ f2g_var[i_var] = rec.PutOp(op);
+ break;
+
+ // --------------------------------------------------------------
+ // Binary operators, left parameter, right variable, one result
+ case AddpvOp:
+ case SubpvOp:
+ case MulpvOp:
+ case DivpvOp:
+ case PowpvOp:
+ case ZmulpvOp:
+ CPPAD_ASSERT_NARG_NRES(op, 2, 1);
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
+ new_arg[0] = rec.PutPar( f_parameter[ arg[0] ] );
+ new_arg[1] = f2g_var[ arg[1] ];
+ rec.PutArg( new_arg[0], new_arg[1] );
+ f2g_var[i_var] = rec.PutOp(op);
+ break;
+ // --------------------------------------------------------------
+ // Binary operators, left and right variables, one result
+ case AddvvOp:
+ case SubvvOp:
+ case MulvvOp:
+ case DivvvOp:
+ case ZmulvvOp:
+ CPPAD_ASSERT_NARG_NRES(op, 2, 1);
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[0] ] ) < num_var );
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
+ new_arg[0] = f2g_var[ arg[0] ];
+ new_arg[1] = f2g_var[ arg[1] ];
+ rec.PutArg( new_arg[0], new_arg[1] );
+ f2g_var[i_var] = rec.PutOp(op);
+ break;
+ // ---------------------------------------------------
+ // Conditional expression operators
+ case CExpOp:
+ CPPAD_ASSERT_NARG_NRES(op, 6, 1);
+ new_arg[0] = arg[0];
+ new_arg[1] = arg[1];
+ mask = 1;
+ for(size_t i = 2; i < 6; i++)
+ { if( arg[1] & mask )
+ { CPPAD_ASSERT_UNKNOWN( size_t(f2g_var[arg[i]]) < num_var );
+ new_arg[i] = f2g_var[ arg[i] ];
+ }
+ else
+ new_arg[i] = rec.PutPar( f_parameter[ arg[i] ] );
+ mask = mask << 1;
+ }
+ rec.PutArg(
+ new_arg[0] ,
+ new_arg[1] ,
+ new_arg[2] ,
+ new_arg[3] ,
+ new_arg[4] ,
+ new_arg[5]
+ );
+ f2g_var[i_var] = rec.PutOp(op);
+ break;
+
+ // --------------------------------------------------
+ // Operators with no arguments and no results
+ case EndOp:
+ CPPAD_ASSERT_NARG_NRES(op, 0, 0);
+ rec.PutOp(op);
+ more_operators = false;
+ break;
+
+ // ---------------------------------------------------
+ // Operations with two arguments and no results
+ case LepvOp:
+ case LtpvOp:
+ case EqpvOp:
+ case NepvOp:
+ CPPAD_ASSERT_NARG_NRES(op, 2, 0);
+ new_arg[0] = rec.PutPar( f_parameter[ arg[0] ] );
+ new_arg[1] = f2g_var[ arg[1] ];
+ rec.PutArg(new_arg[0], new_arg[1]);
+ rec.PutOp(op);
+ break;
+ //
+ case LevpOp:
+ case LtvpOp:
+ CPPAD_ASSERT_NARG_NRES(op, 2, 0);
+ new_arg[0] = f2g_var[ arg[0] ];
+ new_arg[1] = rec.PutPar( f_parameter[ arg[1] ] );
+ rec.PutArg(new_arg[0], new_arg[1]);
+ rec.PutOp(op);
+ break;
+ //
+ case LevvOp:
+ case LtvvOp:
+ case EqvvOp:
+ case NevvOp:
+ CPPAD_ASSERT_NARG_NRES(op, 2, 0);
+ new_arg[0] = f2g_var[ arg[0] ];
+ new_arg[1] = f2g_var[ arg[1] ];
+ rec.PutArg(new_arg[0], new_arg[1]);
+ rec.PutOp(op);
+ break;
+
+ // ---------------------------------------------------
+ // print forward operator
+ case PriOp:
+ CPPAD_ASSERT_NARG_NRES(op, 5, 0);
+ //
+ // arg[0]
+ new_arg[0] = arg[0];
+ //
+ // arg[1]
+ if( arg[0] & 1 )
+ {
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
+ new_arg[1] = f2g_var[ arg[1] ];
+ }
+ else
+ { new_arg[1] = rec.PutPar( f_parameter[ arg[1] ] );
+ }
+ //
+ // arg[3]
+ if( arg[0] & 2 )
+ {
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[3] ] ) < num_var );
+ new_arg[3] = f2g_var[ arg[3] ];
+ }
+ else
+ { new_arg[3] = rec.PutPar( f_parameter[ arg[3] ] );
+ }
+ new_arg[2] = rec.PutTxt( play_.GetTxt( arg[2] ) );
+ new_arg[4] = rec.PutTxt( play_.GetTxt( arg[4] ) );
+ //
+ rec.PutArg(
+ new_arg[0] ,
+ new_arg[1] ,
+ new_arg[2] ,
+ new_arg[3] ,
+ new_arg[4]
+ );
+ // no result
+ rec.PutOp(op);
+ break;
+
+ // ---------------------------------------------------
+ // VecAD operators
+
+ // Load using a parameter index
+ case LdpOp:
+ CPPAD_ASSERT_NARG_NRES(op, 3, 1);
+ new_arg[0] = arg[0];
+ new_arg[1] = arg[1];
+ new_arg[2] = arg[2];
+ rec.PutArg(
+ new_arg[0],
+ new_arg[1],
+ new_arg[2]
+ );
+ f2g_var[i_var] = rec.PutLoadOp(op);
+ break;
+
+ // Load using a variable index
+ case LdvOp:
+ CPPAD_ASSERT_NARG_NRES(op, 3, 1);
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
+ new_arg[0] = arg[0];
+ new_arg[1] = f2g_var[ arg[1] ];
+ new_arg[2] = arg[2];
+ rec.PutArg(
+ new_arg[0],
+ new_arg[1],
+ new_arg[2]
+ );
+ f2g_var[i_var] = rec.PutLoadOp(op);
+ break;
+
+ // Store a parameter using a parameter index
+ case StppOp:
+ CPPAD_ASSERT_NARG_NRES(op, 3, 0);
+ new_arg[0] = arg[0];
+ new_arg[1] = rec.PutPar( f_parameter[ arg[1] ] );
+ new_arg[2] = rec.PutPar( f_parameter[ arg[2] ] );
+ rec.PutArg(
+ new_arg[0],
+ new_arg[1],
+ new_arg[2]
+ );
+ rec.PutOp(op);
+ break;
+
+ // Store a parameter using a variable index
+ case StvpOp:
+ CPPAD_ASSERT_NARG_NRES(op, 3, 0);
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
+ new_arg[0] = arg[0];
+ new_arg[1] = f2g_var[ arg[1] ];
+ new_arg[2] = rec.PutPar( f_parameter[ arg[2] ] );
+ rec.PutArg(
+ new_arg[0],
+ new_arg[1],
+ new_arg[2]
+ );
+ rec.PutOp(op);
+ break;
+
+ // Store a variable using a parameter index
+ case StpvOp:
+ CPPAD_ASSERT_NARG_NRES(op, 3, 0);
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[2] ] ) < num_var );
+ new_arg[0] = arg[0];
+ new_arg[1] = rec.PutPar( f_parameter[ arg[1] ] );
+ new_arg[2] = f2g_var[ arg[2] ];
+ rec.PutArg(
+ new_arg[0],
+ new_arg[1],
+ new_arg[2]
+ );
+ rec.PutOp(op);
+ break;
+
+ // Store a variable using a variable index
+ case StvvOp:
+ CPPAD_ASSERT_NARG_NRES(op, 3, 0);
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[1] ] ) < num_var );
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[ arg[2] ] ) < num_var );
+ new_arg[0] = arg[0];
+ new_arg[1] = f2g_var[ arg[1] ];
+ new_arg[2] = f2g_var[ arg[2] ];
+ rec.PutArg(
+ new_arg[0],
+ new_arg[1],
+ new_arg[2]
+ );
+ break;
+
+ // -----------------------------------------------------------
+ // user atomic function call operators
+
+ case UserOp:
+ CPPAD_ASSERT_NARG_NRES(op, 4, 0);
+ // atomic_index, user_old, user_n, user_m
+ rec.PutArg(arg[0], arg[1], arg[2], arg[3]);
+ rec.PutOp(UserOp);
+ break;
+
+ case UsrapOp:
+ CPPAD_ASSERT_NARG_NRES(op, 1, 0);
+ new_arg[0] = rec.PutPar( f_parameter[ arg[0] ] );
+ rec.PutArg(new_arg[0]);
+ rec.PutOp(UsrapOp);
+ break;
+
+ case UsravOp:
+ CPPAD_ASSERT_NARG_NRES(op, 1, 0);
+ CPPAD_ASSERT_UNKNOWN( size_t( f2g_var[arg[0]] ) < num_var );
+ new_arg[0] = f2g_var[ arg[0] ];
+ rec.PutArg(new_arg[0]);
+ rec.PutOp(UsravOp);
+ break;
+
+ case UsrrpOp:
+ CPPAD_ASSERT_NARG_NRES(op, 1, 0);
+ new_arg[0] = rec.PutPar( f_parameter[ arg[0] ] );
+ rec.PutArg(new_arg[0]);
+ rec.PutOp(UsrrpOp);
+ break;
+
+ case UsrrvOp:
+ CPPAD_ASSERT_NARG_NRES(op, 0, 1);
+ f2g_var[i_var] = rec.PutOp(UsrrvOp);
+ break;
+ // ---------------------------------------------------
+
+ // all cases should be handled above
+ default:
+ CPPAD_ASSERT_UNKNOWN(false);
+ }
+ if( more_operators )
+ play_.forward_next(op, arg, i_op, i_var);
+ }
+ // Check a few expected results
+ CPPAD_ASSERT_UNKNOWN( rec.num_op_rec() == play_.num_op_rec() );
+ CPPAD_ASSERT_UNKNOWN( rec.num_var_rec() == play_.num_var_rec() );
+ CPPAD_ASSERT_UNKNOWN( rec.num_load_op_rec() == play_.num_load_op_rec() );
+
+ // -----------------------------------------------------------------------
+ // Use rec to create the function g
+ // -----------------------------------------------------------------------
+
+ // number of variables in the recording
+ g.num_var_tape_ = rec.num_var_rec();
+
+ // dimension cskip_op vector to number of operators
+ g.cskip_op_.erase();
+ g.cskip_op_.extend( rec.num_op_rec() );
+
+ // independent variables in g: (x, u)
+ size_t s = f_abs_res.size();
+ size_t n = Domain();
+ g.ind_taddr_.resize(n + s);
+ // (x, u)
+ for(size_t j = 0; j < n; j++)
+ { g.ind_taddr_[j] = f2g_var[ ind_taddr_[j] ];
+ CPPAD_ASSERT_UNKNOWN( g.ind_taddr_[j] == j + 1 );
+ }
+ for(size_t j = 0; j < s; j++)
+ { g.ind_taddr_[n + j] = f2g_var[ f_abs_res[j] ];
+ CPPAD_ASSERT_UNKNOWN( g.ind_taddr_[n + j] == n + j + 1 );
+ }
+
+ // dependent variable in g: (y, z)
+ CPPAD_ASSERT_UNKNOWN( s == f_abs_arg.size() );
+ size_t m = Range();
+ g.dep_taddr_.resize(m + s);
+ for(size_t i = 0; i < m; i++)
+ { g.dep_taddr_[i] = f2g_var[ dep_taddr_[i] ];
+ CPPAD_ASSERT_UNKNOWN( g.dep_taddr_[i] < num_var );
+ }
+ for(size_t i = 0; i < s; i++)
+ { g.dep_taddr_[m + i] = f2g_var[ f_abs_arg[i] ];
+ CPPAD_ASSERT_UNKNOWN( g.dep_taddr_[m + i] < num_var );
+ }
+
+ // which dependent variables are parameters
+ g.dep_parameter_.resize(m + s);
+ for(size_t i = 0; i < m; i++)
+ g.dep_parameter_[i] = dep_parameter_[i];
+ for(size_t i = 0; i < s; i++)
+ g.dep_parameter_[m + i] = false;
+
+ // free memory allocated for sparse Jacobian calculation
+ // (the resutls are no longer valid)
+ g.for_jac_sparse_pack_.resize(0, 0);
+ g.for_jac_sparse_set_.resize(0, 0);
+
+ // free taylor coefficient memory
+ g.taylor_.free();
+ g.num_order_taylor_ = 0;
+ g.cap_order_taylor_ = 0;
+
+ // Transferring the recording swaps its vectors so do this last
+ // replace the recording in g (this ADFun object)
+ g.play_.get(rec);
+
+ // ------------------------------------------------------------------------
+ // Create the function a
+ // ------------------------------------------------------------------------
+
+ // start with a copy of f
+ a = *this;
+
+ // dependent variables in a(x)
+ CPPAD_ASSERT_UNKNOWN( s == f_abs_arg.size() );
+ a.dep_taddr_.resize(s);
+ for(size_t i = 0; i < s; i++)
+ { a.dep_taddr_[i] = f_abs_res[i];
+ CPPAD_ASSERT_UNKNOWN( a.dep_taddr_[i] < num_var );
+ }
+
+ // free memory allocated for sparse Jacobian calculation
+ // (the resutls are no longer valid)
+ a.for_jac_sparse_pack_.resize(0, 0);
+ a.for_jac_sparse_set_.resize(0, 0);
+
+ // free taylor coefficient memory
+ a.taylor_.free();
+ a.num_order_taylor_ = 0;
+ a.cap_order_taylor_ = 0;
+}
+
+} // END_CPPAD_NAMESPACE
+
+# endif
diff --git a/external/cppad/include/cppad/core/acosh.hpp b/external/cppad/include/cppad/core/acosh.hpp
new file mode 100644
index 0000000000..7763d8a713
--- /dev/null
+++ b/external/cppad/include/cppad/core/acosh.hpp
@@ -0,0 +1,95 @@
+# ifndef CPPAD_CORE_ACOSH_HPP
+# define CPPAD_CORE_ACOSH_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+-------------------------------------------------------------------------------
+$begin acosh$$
+$spell
+ acosh
+ const
+ Vec
+ std
+ cmath
+ CppAD
+$$
+$section The Inverse Hyperbolic Cosine Function: acosh$$
+
+$head Syntax$$
+$icode%y% = acosh(%x%)%$$
+
+$head Description$$
+The inverse hyperbolic cosine function is defined by
+$icode%x% == cosh(%y%)%$$.
+
+$head x, y$$
+See the $cref/possible types/unary_standard_math/Possible Types/$$
+for a unary standard math function.
+
+$head CPPAD_USE_CPLUSPLUS_2011$$
+
+$subhead true$$
+If this preprocessor symbol is true ($code 1$$),
+and $icode x$$ is an AD type,
+this is an $cref/atomic operation/glossary/Operation/Atomic/$$.
+
+$subhead false$$
+If this preprocessor symbol is false ($code 0$$),
+CppAD uses the representation
+$latex \[
+\R{acosh} (x) = \log \left( x + \sqrt{ x^2 - 1 } \right)
+\] $$
+to compute this function.
+
+$head Example$$
+$children%
+ example/general/acosh.cpp
+%$$
+The file
+$cref acosh.cpp$$
+contains an example and test of this function.
+It returns true if it succeeds and false otherwise.
+
+$end
+-------------------------------------------------------------------------------
+*/
+# include
+# if ! CPPAD_USE_CPLUSPLUS_2011
+
+// BEGIN CppAD namespace
+namespace CppAD {
+
+template
+Type acosh_template(const Type &x)
+{ return CppAD::log( x + CppAD::sqrt( x * x - Type(1) ) );
+}
+
+inline float acosh(const float &x)
+{ return acosh_template(x); }
+
+inline double acosh(const double &x)
+{ return acosh_template(x); }
+
+template
+inline AD acosh(const AD &x)
+{ return acosh_template(x); }
+
+template
+inline AD acosh(const VecAD_reference &x)
+{ return acosh_template( x.ADBase() ); }
+
+
+} // END CppAD namespace
+
+# endif // CPPAD_USE_CPLUSPLUS_2011
+# endif // CPPAD_ACOSH_INCLUDED
diff --git a/external/cppad/include/cppad/core/ad.hpp b/external/cppad/include/cppad/core/ad.hpp
new file mode 100644
index 0000000000..937b7dbb29
--- /dev/null
+++ b/external/cppad/include/cppad/core/ad.hpp
@@ -0,0 +1,291 @@
+# ifndef CPPAD_CORE_AD_HPP
+# define CPPAD_CORE_AD_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+// simple AD operations that must be defined for AD as well as base class
+# include
+# include
+
+// define the template classes that are used by the AD template class
+# include
+# include
+# include
+# include
+
+namespace CppAD { // BEGIN_CPPAD_NAMESPACE
+
+typedef enum {
+ tape_manage_new,
+ tape_manage_delete,
+ tape_manage_clear
+} tape_manage_job;
+
+template
+class AD {
+private :
+ // -----------------------------------------------------------------------
+ // value_ corresponding to this object
+ Base value_;
+
+ // Tape identifier corresponding to taddr
+ tape_id_t tape_id_;
+
+ // taddr_ in tape for this variable
+ addr_t taddr_;
+ // -----------------------------------------------------------------------
+
+ // enable use of AD in parallel mode
+ template
+ friend void parallel_ad(void);
+
+ // template friend functions where template parameter is not bound
+ template
+ friend void Independent(VectorAD &x, size_t abort_op_index);
+
+ // one argument functions
+ friend bool Parameter
+ (const AD &u);
+ friend bool Parameter
+ (const VecAD &u);
+ friend bool Variable
+ (const AD &u);
+ friend bool Variable
+ (const VecAD &u);
+ friend int Integer
+ (const AD &u);
+ friend AD Var2Par
+ (const AD &u);
+
+ // power function
+ friend AD pow
+ (const AD &x, const AD &y);
+
+ // azmul function
+ friend AD azmul
+ (const AD &x, const AD &y);
+
+ // order determining functions, see ordered.hpp
+ friend bool GreaterThanZero (const AD &x);
+ friend bool GreaterThanOrZero (const AD &x);
+ friend bool LessThanZero (const AD &x);
+ friend bool LessThanOrZero (const AD &x);
+ friend bool abs_geq
+ (const AD& x, const AD& y);
+
+ // The identical property functions, see identical.hpp
+ friend bool IdenticalPar (const AD &x);
+ friend bool IdenticalZero (const AD &x);
+ friend bool IdenticalOne (const AD &x);
+ friend bool IdenticalEqualPar
+ (const AD &x, const AD &y);
+
+ // EqualOpSeq function
+ friend bool EqualOpSeq
+ (const AD &u, const AD &v);
+
+ // NearEqual function
+ friend bool NearEqual (
+ const AD &x, const AD &y, const Base &r, const Base &a);
+
+ friend bool NearEqual (
+ const Base &x, const AD &y, const Base &r, const Base &a);
+
+ friend bool NearEqual (
+ const AD &x, const Base &y, const Base &r, const Base &a);
+
+ // CondExp function
+ friend AD CondExpOp (
+ enum CompareOp cop ,
+ const AD &left ,
+ const AD &right ,
+ const AD &trueCase ,
+ const AD &falseCase
+ );
+
+ // classes
+ friend class local::ADTape;
+ friend class ADFun;
+ friend class atomic_base;
+ friend class discrete;
+ friend class VecAD;
+ friend class VecAD_reference;
+
+ // arithematic binary operators
+ friend AD operator +
+ (const AD &left, const AD &right);
+ friend AD operator -
+ (const AD &left, const AD &right);
+ friend AD operator *
+ (const AD &left, const AD &right);
+ friend AD operator /
+ (const AD &left, const AD &right);
+
+ // comparison operators
+ friend bool operator <
+ (const AD &left, const AD &right);
+ friend bool operator <=
+ (const AD &left, const AD &right);
+ friend bool operator >
+ (const AD &left, const AD &right);
+ friend bool operator >=
+ (const AD &left, const AD &right);
+ friend bool operator ==
+ (const AD &left, const AD &right);
+ friend bool operator !=
+ (const AD &left, const AD &right);
+
+ // input operator
+ friend std::istream& operator >>
+ (std::istream &is, AD &x);
+
+ // output operations
+ friend std::ostream& operator <<
+ (std::ostream &os, const AD &x);
+ friend void PrintFor (
+ const AD& flag ,
+ const char* before ,
+ const AD& var ,
+ const char* after
+ );
+public:
+ // type of value
+ typedef Base value_type;
+
+ // implicit default constructor
+ inline AD(void);
+
+ // use default implicit copy constructor and assignment operator
+ // inline AD(const AD &x);
+ // inline AD& operator=(const AD &x);
+
+ // implicit construction and assingment from base type
+ inline AD(const Base &b);
+ inline AD& operator=(const Base &b);
+
+ // implicit contructor and assignment from VecAD::reference
+ inline AD(const VecAD_reference &x);
+ inline AD& operator=(const VecAD_reference &x);
+
+ // explicit construction from some other type (depricated)
+ template inline explicit AD(const T &t);
+
+ // assignment from some other type
+ template inline AD& operator=(const T &right);
+
+ // base type corresponding to an AD object
+ friend Base Value (const AD &x);
+
+ // compound assignment operators
+ inline AD& operator += (const AD &right);
+ inline AD& operator -= (const AD &right);
+ inline AD& operator *= (const AD &right);
+ inline AD& operator /= (const AD &right);
+
+ // unary operators
+ inline AD operator +(void) const;
+ inline AD operator -(void) const;
+
+ // destructor
+ ~AD(void)
+ { }
+
+ // interface so these functions need not be friends
+ inline AD abs_me(void) const;
+ inline AD acos_me(void) const;
+ inline AD asin_me(void) const;
+ inline AD atan_me(void) const;
+ inline AD cos_me(void) const;
+ inline AD cosh_me(void) const;
+ inline AD exp_me(void) const;
+ inline AD fabs_me(void) const;
+ inline AD log_me(void) const;
+ inline AD sin_me(void) const;
+ inline AD sign_me(void) const;
+ inline AD sinh_me(void) const;
+ inline AD sqrt_me(void) const;
+ inline AD tan_me(void) const;
+ inline AD tanh_me(void) const;
+# if CPPAD_USE_CPLUSPLUS_2011
+ inline AD erf_me(void) const;
+ inline AD asinh_me(void) const;
+ inline AD acosh_me(void) const;
+ inline AD atanh_me(void) const;
+ inline AD expm1_me(void) const;
+ inline AD log1p_me(void) const;
+# endif
+
+ // ----------------------------------------------------------
+ // static public member functions
+
+ // abort current AD recording
+ static void abort_recording(void);
+
+ // set the maximum number of OpenMP threads (deprecated)
+ static void omp_max_thread(size_t number);
+
+ // These functions declared public so can be accessed by user through
+ // a macro interface and are not intended for direct use.
+ // The macro interface is documented in bool_fun.hpp.
+ // Developer documentation for these fucntions is in bool_fun.hpp
+ static inline bool UnaryBool(
+ bool FunName(const Base &x),
+ const AD &x
+ );
+ static inline bool BinaryBool(
+ bool FunName(const Base &x, const Base &y),
+ const AD &x , const AD &y
+ );
+
+private:
+ //
+ // Make this variable a parameter
+ //
+ void make_parameter(void)
+ { CPPAD_ASSERT_UNKNOWN( Variable(*this) ); // currently a var
+ tape_id_ = 0;
+ }
+ //
+ // Make this parameter a new variable
+ //
+ void make_variable(tape_id_t id, addr_t taddr)
+ { CPPAD_ASSERT_UNKNOWN( Parameter(*this) ); // currently a par
+ CPPAD_ASSERT_UNKNOWN( taddr > 0 ); // sure valid taddr
+
+ taddr_ = taddr;
+ tape_id_ = id;
+ }
+ // ---------------------------------------------------------------
+ // tape linking functions
+ //
+ // not static
+ inline local::ADTape* tape_this(void) const;
+ //
+ // static
+ inline static tape_id_t** tape_id_handle(size_t thread);
+ inline static tape_id_t* tape_id_ptr(size_t thread);
+ inline static local::ADTape** tape_handle(size_t thread);
+ static local::ADTape* tape_manage(tape_manage_job job);
+ inline static local::ADTape* tape_ptr(void);
+ inline static local::ADTape* tape_ptr(tape_id_t tape_id);
+};
+// ---------------------------------------------------------------------------
+
+} // END_CPPAD_NAMESPACE
+
+// tape linking private functions
+# include
+
+// operations that expect the AD template class to be defined
+
+
+# endif
diff --git a/external/cppad/include/cppad/core/ad_assign.hpp b/external/cppad/include/cppad/core/ad_assign.hpp
new file mode 100644
index 0000000000..735706310d
--- /dev/null
+++ b/external/cppad/include/cppad/core/ad_assign.hpp
@@ -0,0 +1,140 @@
+# ifndef CPPAD_CORE_AD_ASSIGN_HPP
+# define CPPAD_CORE_AD_ASSIGN_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+------------------------------------------------------------------------------
+
+$begin ad_assign$$
+$spell
+ Vec
+ const
+$$
+
+
+$section AD Assignment Operator$$
+$mindex assign Base VecAD$$
+
+$head Syntax$$
+$icode%y% = %x%$$
+
+$head Purpose$$
+Assigns the value in $icode x$$ to the object $icode y$$.
+In either case,
+
+$head x$$
+The argument $icode x$$ has prototype
+$codei%
+ const %Type% &%x%
+%$$
+where $icode Type$$ is
+$codei%VecAD<%Base%>::reference%$$,
+$codei%AD<%Base%>%$$,
+$icode Base$$,
+or any type that has an implicit constructor of the form
+$icode%Base%(%x%)%$$.
+
+$head y$$
+The target $icode y$$ has prototype
+$codei%
+ AD<%Base%> %y%
+%$$
+
+$head Example$$
+$children%
+ example/general/ad_assign.cpp
+%$$
+The file $cref ad_assign.cpp$$ contain examples and tests of these operations.
+It test returns true if it succeeds and false otherwise.
+
+$end
+------------------------------------------------------------------------------
+*/
+
+namespace CppAD { // BEGIN_CPPAD_NAMESPACE
+
+/*!
+\file ad_assign.hpp
+AD constructors and and copy operations.
+*/
+
+/*!
+\page AD_default_assign
+Use default assignment operator
+because they may be optimized better than the code below:
+\code
+template
+inline AD& AD::operator=(const AD &right)
+{ value_ = right.value_;
+ tape_id_ = right.tape_id_;
+ taddr_ = right.taddr_;
+
+ return *this;
+}
+\endcode
+*/
+
+/*!
+Assignment to Base type value.
+
+\tparam Base
+Base type for this AD object.
+
+\param b
+is the Base type value being assignment to this AD object.
+The tape identifier will be an invalid tape identifier,
+so this object is initially a parameter.
+*/
+template
+inline AD& AD::operator=(const Base &b)
+{ value_ = b;
+ tape_id_ = 0;
+
+ // check that this is a parameter
+ CPPAD_ASSERT_UNKNOWN( Parameter(*this) );
+
+ return *this;
+}
+
+/*!
+Assignment to an ADVec element drops the vector information.
+
+\tparam Base
+Base type for this AD object.
+*/
+template
+inline AD& AD::operator=(const VecAD_reference &x)
+{ return *this = x.ADBase(); }
+
+/*!
+Assignment from any other type, converts to Base type, and then uses assignment
+from Base type.
+
+\tparam Base
+Base type for this AD object.
+
+\tparam T
+is the the type that is being assigned to AD.
+There must be an assignment for Base from Type.
+
+\param t
+is the object that is being assigned to an AD object.
+*/
+template
+template
+inline AD& AD::operator=(const T &t)
+{ return *this = Base(t); }
+
+
+} // END_CPPAD_NAMESPACE
+# endif
diff --git a/external/cppad/include/cppad/core/ad_binary.hpp b/external/cppad/include/cppad/core/ad_binary.hpp
new file mode 100644
index 0000000000..03bbeb8091
--- /dev/null
+++ b/external/cppad/include/cppad/core/ad_binary.hpp
@@ -0,0 +1,144 @@
+# ifndef CPPAD_CORE_AD_BINARY_HPP
+# define CPPAD_CORE_AD_BINARY_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+-------------------------------------------------------------------------------
+$begin ad_binary$$
+$spell
+ Op
+ VecAD
+ const
+$$
+
+$section AD Binary Arithmetic Operators$$
+$mindex + add plus - subtract minus * multiply times / divide$$
+
+
+
+
+
+
+
+$head Syntax$$
+$icode%z% = %x% %Op% %y%$$
+
+$head Purpose$$
+Performs arithmetic operations where either $icode x$$ or $icode y$$
+has type
+$codei%AD<%Base%>%$$ or
+$cref%VecAD::reference%VecAD%VecAD::reference%$$.
+
+$head Op$$
+The operator $icode Op$$ is one of the following
+$table
+$bold Op$$ $cnext $bold Meaning$$ $rnext
+$code +$$ $cnext $icode z$$ is $icode x$$ plus $icode y$$ $rnext
+$code -$$ $cnext $icode z$$ is $icode x$$ minus $icode y$$ $rnext
+$code *$$ $cnext $icode z$$ is $icode x$$ times $icode y$$ $rnext
+$code /$$ $cnext $icode z$$ is $icode x$$ divided by $icode y$$
+$tend
+
+$head Base$$
+The type $icode Base$$ is determined by the operand that
+has type $codei%AD<%Base%>%$$ or $codei%VecAD<%Base%>::reference%$$.
+
+$head x$$
+The operand $icode x$$ has the following prototype
+$codei%
+ const %Type% &%x%
+%$$
+where $icode Type$$ is
+$codei%VecAD<%Base%>::reference%$$,
+$codei%AD<%Base%>%$$,
+$icode Base$$, or
+$code double$$.
+
+$head y$$
+The operand $icode y$$ has the following prototype
+$codei%
+ const %Type% &%y%
+%$$
+where $icode Type$$ is
+$codei%VecAD<%Base%>::reference%$$,
+$codei%AD<%Base%>%$$,
+$icode Base$$, or
+$code double$$.
+
+
+$head z$$
+The result $icode z$$ has the following prototype
+$codei%
+ %Type% %z%
+%$$
+where $icode Type$$ is
+$codei%AD<%Base%>%$$.
+
+$head Operation Sequence$$
+This is an $cref/atomic/glossary/Operation/Atomic/$$
+$cref/AD of Base/glossary/AD of Base/$$ operation
+and hence it is part of the current
+AD of $icode Base$$
+$cref/operation sequence/glossary/Operation/Sequence/$$.
+
+$children%
+ example/general/add.cpp%
+ example/general/sub.cpp%
+ example/general/mul.cpp%
+ example/general/div.cpp
+%$$
+
+$head Example$$
+The following files contain examples and tests of these functions.
+Each test returns true if it succeeds and false otherwise.
+$table
+$rref add.cpp$$
+$rref sub.cpp$$
+$rref mul.cpp$$
+$rref div.cpp$$
+$tend
+
+$head Derivative$$
+If $latex f$$ and $latex g$$ are
+$cref/Base functions/glossary/Base Function/$$
+
+$subhead Addition$$
+$latex \[
+ \D{[ f(x) + g(x) ]}{x} = \D{f(x)}{x} + \D{g(x)}{x}
+\] $$
+
+$subhead Subtraction$$
+$latex \[
+ \D{[ f(x) - g(x) ]}{x} = \D{f(x)}{x} - \D{g(x)}{x}
+\] $$
+
+$subhead Multiplication$$
+$latex \[
+ \D{[ f(x) * g(x) ]}{x} = g(x) * \D{f(x)}{x} + f(x) * \D{g(x)}{x}
+\] $$
+
+$subhead Division$$
+$latex \[
+ \D{[ f(x) / g(x) ]}{x} =
+ [1/g(x)] * \D{f(x)}{x} - [f(x)/g(x)^2] * \D{g(x)}{x}
+\] $$
+
+$end
+-----------------------------------------------------------------------------
+*/
+# include
+# include
+# include
+# include
+
+# endif
diff --git a/external/cppad/include/cppad/core/ad_ctor.hpp b/external/cppad/include/cppad/core/ad_ctor.hpp
new file mode 100644
index 0000000000..ab10718d1a
--- /dev/null
+++ b/external/cppad/include/cppad/core/ad_ctor.hpp
@@ -0,0 +1,167 @@
+# ifndef CPPAD_CORE_AD_CTOR_HPP
+# define CPPAD_CORE_AD_CTOR_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+------------------------------------------------------------------------------
+
+$begin ad_ctor$$
+$spell
+ cppad
+ ctor
+ initializes
+ Vec
+ const
+$$
+
+
+$section AD Constructors $$
+$mindex convert Base VecAD$$
+
+$head Syntax$$
+$codei%AD<%Base%> %y%()
+%$$
+$codei%AD<%Base%> %y%(%x%)
+%$$
+
+$head Purpose$$
+creates a new $codei%AD<%Base%>%$$ object $icode y$$
+and initializes its value as equal to $icode x$$.
+
+$head x$$
+
+$subhead implicit$$
+There is an implicit constructor where $icode x$$ has one of the following
+prototypes:
+$codei%
+ const %Base%& %x%
+ const VecAD<%Base%>& %x%
+%$$
+
+$subhead explicit$$
+There is an explicit constructor where $icode x$$ has prototype
+$codei%
+ const %Type%& %x%
+%$$
+for any type that has an explicit constructor of the form
+$icode%Base%(%x%)%$$.
+
+$head y$$
+The target $icode y$$ has prototype
+$codei%
+ AD<%Base%> %y%
+%$$
+
+$head Example$$
+$children%
+ example/general/ad_ctor.cpp
+%$$
+The files $cref ad_ctor.cpp$$ contain examples and tests of these operations.
+It test returns true if it succeeds and false otherwise.
+
+$end
+------------------------------------------------------------------------------
+*/
+
+namespace CppAD { // BEGIN_CPPAD_NAMESPACE
+
+/*!
+\file ad_ctor.hpp
+AD constructors and and copy operations.
+*/
+
+/*!
+\page AD_default_ctor
+Use default copy constructor
+because they may be optimized better than the code below:
+\code
+template
+inline AD::AD(const AD &x)
+{
+ value_ = x.value_;
+ tape_id_ = x.tape_id_;
+ taddr_ = x.taddr_;
+
+ return;
+}
+\endcode
+*/
+
+/*!
+Default Constructor.
+
+\tparam Base
+Base type for this AD object.
+*/
+template
+inline AD::AD(void)
+: value_()
+, tape_id_(0)
+, taddr_(0)
+{ }
+
+
+/*!
+Constructor from Base type.
+
+\tparam Base
+Base type for this AD object.
+
+\param b
+is the Base type value corresponding to this AD object.
+The tape identifier will be an invalid tape identifier,
+so this object is initially a parameter.
+*/
+template
+inline AD::AD(const Base &b)
+: value_(b)
+, tape_id_(0)
+, taddr_(0)
+{ // check that this is a parameter
+ CPPAD_ASSERT_UNKNOWN( Parameter(*this) );
+}
+
+/*!
+Constructor from an ADVec element drops the vector information.
+
+\tparam Base
+Base type for this AD object.
+*/
+template
+inline AD::AD(const VecAD_reference &x)
+{ *this = x.ADBase(); }
+
+/*!
+Constructor from any other type, converts to Base type, and uses constructor
+from Base type.
+
+\tparam Base
+Base type for this AD object.
+
+\tparam T
+is the the type that is being converted to AD.
+There must be a constructor for Base from Type.
+
+\param t
+is the object that is being converted from T to AD.
+*/
+template
+template
+inline AD::AD(const T &t)
+: value_(Base(t))
+, tape_id_(0)
+, taddr_(0)
+{ }
+
+} // END_CPPAD_NAMESPACE
+# endif
diff --git a/external/cppad/include/cppad/core/ad_fun.hpp b/external/cppad/include/cppad/core/ad_fun.hpp
new file mode 100644
index 0000000000..0388beb5fc
--- /dev/null
+++ b/external/cppad/include/cppad/core/ad_fun.hpp
@@ -0,0 +1,712 @@
+# ifndef CPPAD_CORE_AD_FUN_HPP
+# define CPPAD_CORE_AD_FUN_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+/*
+$begin ADFun$$
+$spell
+ xk
+ Ind
+ bool
+ taylor_
+ sizeof
+ const
+ std
+ ind_taddr_
+ dep_taddr_
+$$
+
+$spell
+$$
+
+$section ADFun Objects$$
+
+
+$head Purpose$$
+An AD of $icode Base$$
+$cref/operation sequence/glossary/Operation/Sequence/$$
+is stored in an $code ADFun$$ object by its $cref FunConstruct$$.
+The $code ADFun$$ object can then be used to calculate function values,
+derivative values, and other values related to the corresponding function.
+
+$childtable%
+ omh/adfun.omh%
+ cppad/core/optimize.hpp%
+ example/abs_normal/abs_normal.omh%
+ cppad/core/fun_check.hpp%
+ cppad/core/check_for_nan.hpp
+%$$
+
+$end
+*/
+
+namespace CppAD { // BEGIN_CPPAD_NAMESPACE
+/*!
+\file ad_fun.hpp
+File used to define the ADFun class.
+*/
+
+/*!
+Class used to hold function objects
+
+\tparam Base
+A function object has a recording of AD operations.
+It does it calculations using \c Base operations.
+*/
+
+template
+class ADFun {
+// ------------------------------------------------------------
+// Private member variables
+private:
+ /// Has this ADFun object been optmized
+ bool has_been_optimized_;
+
+ /// Check for nan's and report message to user (default value is true).
+ bool check_for_nan_;
+
+ /// If zero, ignoring comparison operators. Otherwise is the
+ /// compare change count at which to store the operator index.
+ size_t compare_change_count_;
+
+ /// If compare_change_count_ is zero, compare_change_number_ is also zero.
+ /// Otherwise, it is set to the number of comparison operations that had a
+ /// different result during the subsequent zero order forward.
+ size_t compare_change_number_;
+
+ /// If compare_change_count is zero, compare_change_op_index_ is also
+ /// zero. Otherwise it is the operator index for the comparison operator
+ //// that corresponded to the number changing from count-1 to count.
+ size_t compare_change_op_index_;
+
+ /// number of orders stored in taylor_
+ size_t num_order_taylor_;
+
+ /// maximum number of orders that will fit in taylor_
+ size_t cap_order_taylor_;
+
+ /// number of directions stored in taylor_
+ size_t num_direction_taylor_;
+
+ /// number of variables in the recording (play_)
+ size_t num_var_tape_;
+
+ /// tape address for the independent variables
+ CppAD::vector ind_taddr_;
+
+ /// tape address and parameter flag for the dependent variables
+ CppAD::vector dep_taddr_;
+
+ /// which dependent variables are actually parameters
+ CppAD::vector dep_parameter_;
+
+ /// results of the forward mode calculations
+ local::pod_vector taylor_;
+
+ /// which operations can be conditionally skipped
+ /// Set during forward pass of order zero
+ local::pod_vector cskip_op_;
+
+ /// Variable on the tape corresponding to each vecad load operation
+ /// (if zero, the operation corresponds to a parameter).
+ local::pod_vector load_op_;
+
+ /// the operation sequence corresponding to this object
+ local::player play_;
+
+ /// Packed results of the forward mode Jacobian sparsity calculations.
+ /// for_jac_sparse_pack_.n_set() != 0 implies other sparsity results
+ /// are empty
+ local::sparse_pack for_jac_sparse_pack_;
+
+ /// Set results of the forward mode Jacobian sparsity calculations
+ /// for_jac_sparse_set_.n_set() != 0 implies for_sparse_pack_ is empty.
+ local::sparse_list for_jac_sparse_set_;
+
+// ------------------------------------------------------------
+// Private member functions
+
+ /// change the operation sequence corresponding to this object
+ template
+ void Dependent(local::ADTape *tape, const ADvector &y);
+
+ // ------------------------------------------------------------
+ // vector of bool version of ForSparseJac
+ // (see doxygen in for_sparse_jac.hpp)
+ template
+ void ForSparseJacCase(
+ bool set_type ,
+ bool transpose ,
+ bool dependency,
+ size_t q ,
+ const VectorSet& r ,
+ VectorSet& s
+ );
+ // vector of std::set version of ForSparseJac
+ // (see doxygen in for_sparse_jac.hpp)
+ template
+ void ForSparseJacCase(
+ const std::set& set_type ,
+ bool transpose ,
+ bool dependency,
+ size_t q ,
+ const VectorSet& r ,
+ VectorSet& s
+ );
+ // ------------------------------------------------------------
+ // vector of bool version of RevSparseJac
+ // (see doxygen in rev_sparse_jac.hpp)
+ template
+ void RevSparseJacCase(
+ bool set_type ,
+ bool transpose ,
+ bool dependency,
+ size_t p ,
+ const VectorSet& s ,
+ VectorSet& r
+ );
+ // vector of std::set version of RevSparseJac
+ // (see doxygen in rev_sparse_jac.hpp)
+ template
+ void RevSparseJacCase(
+ const std::set& set_type ,
+ bool transpose ,
+ bool dependency,
+ size_t p ,
+ const VectorSet& s ,
+ VectorSet& r
+ );
+ // ------------------------------------------------------------
+ // vector of bool version of ForSparseHes
+ // (see doxygen in rev_sparse_hes.hpp)
+ template
+ void ForSparseHesCase(
+ bool set_type ,
+ const VectorSet& r ,
+ const VectorSet& s ,
+ VectorSet& h
+ );
+ // vector of std::set version of ForSparseHes
+ // (see doxygen in rev_sparse_hes.hpp)
+ template
+ void ForSparseHesCase(
+ const std::set& set_type ,
+ const VectorSet& r ,
+ const VectorSet& s ,
+ VectorSet& h
+ );
+ // ------------------------------------------------------------
+ // vector of bool version of RevSparseHes
+ // (see doxygen in rev_sparse_hes.hpp)
+ template
+ void RevSparseHesCase(
+ bool set_type ,
+ bool transpose ,
+ size_t q ,
+ const VectorSet& s ,
+ VectorSet& h
+ );
+ // vector of std::set version of RevSparseHes
+ // (see doxygen in rev_sparse_hes.hpp)
+ template
+ void RevSparseHesCase(
+ const std::set& set_type ,
+ bool transpose ,
+ size_t q ,
+ const VectorSet& s ,
+ VectorSet& h
+ );
+ // ------------------------------------------------------------
+ // Forward mode version of SparseJacobian
+ // (see doxygen in sparse_jacobian.hpp)
+ template
+ size_t SparseJacobianFor(
+ const VectorBase& x ,
+ VectorSet& p_transpose ,
+ const VectorSize& row ,
+ const VectorSize& col ,
+ VectorBase& jac ,
+ sparse_jacobian_work& work
+ );
+ // Reverse mode version of SparseJacobian
+ // (see doxygen in sparse_jacobian.hpp)
+ template
+ size_t SparseJacobianRev(
+ const VectorBase& x ,
+ VectorSet& p ,
+ const VectorSize& row ,
+ const VectorSize& col ,
+ VectorBase& jac ,
+ sparse_jacobian_work& work
+ );
+ // ------------------------------------------------------------
+ // combined sparse_list and sparse_pack version of
+ // SparseHessian (see doxygen in sparse_hessian.hpp)
+ template
+ size_t SparseHessianCompute(
+ const VectorBase& x ,
+ const VectorBase& w ,
+ VectorSet& sparsity ,
+ const VectorSize& row ,
+ const VectorSize& col ,
+ VectorBase& hes ,
+ sparse_hessian_work& work
+ );
+// ------------------------------------------------------------
+public:
+ /// copy constructor
+ ADFun(const ADFun& g)
+ : num_var_tape_(0)
+ { CppAD::ErrorHandler::Call(
+ true,
+ __LINE__,
+ __FILE__,
+ "ADFun(const ADFun& g)",
+ "Attempting to use the ADFun copy constructor.\n"
+ "Perhaps you are passing an ADFun object "
+ "by value instead of by reference."
+ );
+ }
+
+ /// default constructor
+ ADFun(void);
+
+ // assignment operator
+ // (see doxygen in fun_construct.hpp)
+ void operator=(const ADFun& f);
+
+ /// sequence constructor
+ template
+ ADFun(const ADvector &x, const ADvector &y);
+
+ /// destructor
+ ~ADFun(void)
+ { }
+
+ /// set value of check_for_nan_
+ void check_for_nan(bool value)
+ { check_for_nan_ = value; }
+ bool check_for_nan(void) const
+ { return check_for_nan_; }
+
+ /// assign a new operation sequence
+ template
+ void Dependent(const ADvector &x, const ADvector &y);
+
+ /// forward mode user API, one order multiple directions.
+ template
+ VectorBase Forward(size_t q, size_t r, const VectorBase& x);
+
+ /// forward mode user API, multiple directions one order.
+ template
+ VectorBase Forward(size_t q,
+ const VectorBase& x, std::ostream& s = std::cout
+ );
+
+ /// reverse mode sweep
+ template
+ VectorBase Reverse(size_t p, const VectorBase &v);
+
+ // ---------------------------------------------------------------------
+ // Jacobian sparsity
+ template
+ VectorSet ForSparseJac(
+ size_t q, const VectorSet &r, bool transpose = false,
+ bool dependency = false
+ );
+ template
+ VectorSet RevSparseJac(
+ size_t q, const VectorSet &s, bool transpose = false,
+ bool dependency = false
+ );
+ // ---------------------------------------------------------------------
+ template
+ size_t sparse_jac_for(
+ size_t group_max ,
+ const BaseVector& x ,
+ sparse_rcv& subset ,
+ const sparse_rc& pattern ,
+ const std::string& coloring ,
+ sparse_jac_work& work
+ );
+ template
+ size_t sparse_jac_rev(
+ const BaseVector& x ,
+ sparse_rcv& subset ,
+ const sparse_rc& pattern ,
+ const std::string& coloring ,
+ sparse_jac_work& work
+ );
+ template
+ size_t sparse_hes(
+ const BaseVector& x ,
+ const BaseVector& w ,
+ sparse_rcv& subset ,
+ const sparse_rc& pattern ,
+ const std::string& coloring ,
+ sparse_hes_work& work
+ );
+ // ---------------------------------------------------------------------
+ template
+ void for_jac_sparsity(
+ const sparse_rc& pattern_in ,
+ bool transpose ,
+ bool dependency ,
+ bool internal_bool ,
+ sparse_rc& pattern_out
+ );
+ template
+ void rev_jac_sparsity(
+ const sparse_rc& pattern_in ,
+ bool transpose ,
+ bool dependency ,
+ bool internal_bool ,
+ sparse_rc& pattern_out
+ );
+ template
+ void rev_hes_sparsity(
+ const BoolVector& select_range ,
+ bool transpose ,
+ bool internal_bool ,
+ sparse_rc& pattern_out
+ );
+ template
+ void for_hes_sparsity(
+ const BoolVector& select_domain ,
+ const BoolVector& select_range ,
+ bool internal_bool ,
+ sparse_rc& pattern_out
+ );
+ // ---------------------------------------------------------------------
+ // forward mode Hessian sparsity
+ // (see doxygen documentation in rev_sparse_hes.hpp)
+ template
+ VectorSet ForSparseHes(
+ const VectorSet &r, const VectorSet &s
+ );
+ // internal set sparsity version of ForSparseHes
+ // (used by checkpoint functions only)
+ void ForSparseHesCheckpoint(
+ vector& r ,
+ vector& s ,
+ local::sparse_list& h
+ );
+ // reverse mode Hessian sparsity
+ // (see doxygen documentation in rev_sparse_hes.hpp)
+ template
+ VectorSet RevSparseHes(
+ size_t q, const VectorSet &s, bool transpose = false
+ );
+ // internal set sparsity version of RevSparseHes
+ // (used by checkpoint functions only)
+ void RevSparseHesCheckpoint(
+ size_t q ,
+ vector& s ,
+ bool transpose ,
+ local::sparse_list& h
+ );
+ // internal set sparsity version of RevSparseJac
+ // (used by checkpoint functions only)
+ void RevSparseJacCheckpoint(
+ size_t q ,
+ const local::sparse_list& r ,
+ bool transpose ,
+ bool dependency ,
+ local::sparse_list& s
+ );
+ // internal set sparsity version of RevSparseJac
+ // (used by checkpoint functions only)
+ void ForSparseJacCheckpoint(
+ size_t q ,
+ const local::sparse_list& r ,
+ bool transpose ,
+ bool dependency ,
+ local::sparse_list& s
+ );
+
+ /// amount of memory used for boolean Jacobain sparsity pattern
+ size_t size_forward_bool(void) const
+ { return for_jac_sparse_pack_.memory(); }
+
+ /// free memory used for Jacobain sparsity pattern
+ void size_forward_bool(size_t zero)
+ { CPPAD_ASSERT_KNOWN(
+ zero == 0,
+ "size_forward_bool: argument not equal to zero"
+ );
+ for_jac_sparse_pack_.resize(0, 0);
+ }
+
+ /// amount of memory used for vector of set Jacobain sparsity pattern
+ size_t size_forward_set(void) const
+ { return for_jac_sparse_set_.memory(); }
+
+ /// free memory used for Jacobain sparsity pattern
+ void size_forward_set(size_t zero)
+ { CPPAD_ASSERT_KNOWN(
+ zero == 0,
+ "size_forward_bool: argument not equal to zero"
+ );
+ for_jac_sparse_set_.resize(0, 0);
+ }
+
+ /// number of operators in the operation sequence
+ size_t size_op(void) const
+ { return play_.num_op_rec(); }
+
+ /// number of operator arguments in the operation sequence
+ size_t size_op_arg(void) const
+ { return play_.num_op_arg_rec(); }
+
+ /// amount of memory required for the operation sequence
+ size_t size_op_seq(void) const
+ { return play_.Memory(); }
+
+ /// number of parameters in the operation sequence
+ size_t size_par(void) const
+ { return play_.num_par_rec(); }
+
+ /// number taylor coefficient orders calculated
+ size_t size_order(void) const
+ { return num_order_taylor_; }
+
+ /// number taylor coefficient directions calculated
+ size_t size_direction(void) const
+ { return num_direction_taylor_; }
+
+ /// number of characters in the operation sequence
+ size_t size_text(void) const
+ { return play_.num_text_rec(); }
+
+ /// number of variables in opertion sequence
+ size_t size_var(void) const
+ { return num_var_tape_; }
+
+ /// number of VecAD indices in the operation sequence
+ size_t size_VecAD(void) const
+ { return play_.num_vec_ind_rec(); }
+
+ /// set number of orders currently allocated (user API)
+ void capacity_order(size_t c);
+
+ /// set number of orders and directions currently allocated
+ void capacity_order(size_t c, size_t r);
+
+ /// number of variables in conditional expressions that can be skipped
+ size_t number_skip(void);
+
+ /// number of independent variables
+ size_t Domain(void) const
+ { return ind_taddr_.size(); }
+
+ /// number of dependent variables
+ size_t Range(void) const
+ { return dep_taddr_.size(); }
+
+ /// is variable a parameter
+ bool Parameter(size_t i)
+ { CPPAD_ASSERT_KNOWN(
+ i < dep_taddr_.size(),
+ "Argument to Parameter is >= dimension of range space"
+ );
+ return dep_parameter_[i];
+ }
+
+ /// Deprecated: number of comparison operations that changed
+ /// for the previous zero order forward (than when function was recorded)
+ size_t CompareChange(void) const
+ { return compare_change_number_; }
+
+ /// count as which to store operator index
+ void compare_change_count(size_t count)
+ { compare_change_count_ = count;
+ compare_change_number_ = 0;
+ compare_change_op_index_ = 0;
+ }
+
+ /// number of comparison operations that changed
+ size_t compare_change_number(void) const
+ { return compare_change_number_; }
+
+ /// operator index for the count-th comparison change
+ size_t compare_change_op_index(void) const
+ { if( has_been_optimized_ )
+ return 0;
+ return compare_change_op_index_;
+ }
+
+ /// calculate entire Jacobian
+ template
+ VectorBase Jacobian(const VectorBase &x);
+
+ /// calculate Hessian for one component of f
+ template
+ VectorBase Hessian(const VectorBase &x, const VectorBase &w);
+ template
+ VectorBase Hessian(const VectorBase &x, size_t i);
+
+ /// forward mode calculation of partial w.r.t one domain component
+ template
+ VectorBase ForOne(
+ const VectorBase &x ,
+ size_t j );
+
+ /// reverse mode calculation of derivative of one range component
+ template
+ VectorBase RevOne(
+ const VectorBase &x ,
+ size_t i );
+
+ /// forward mode calculation of a subset of second order partials
+ template
+ VectorBase ForTwo(
+ const VectorBase &x ,
+ const VectorSize_t &J ,
+ const VectorSize_t &K );
+
+ /// reverse mode calculation of a subset of second order partials
+ template
+ VectorBase RevTwo(
+ const VectorBase &x ,
+ const VectorSize_t &I ,
+ const VectorSize_t &J );
+
+ /// calculate sparse Jacobians
+ template
+ VectorBase SparseJacobian(
+ const VectorBase &x
+ );
+ template
+ VectorBase SparseJacobian(
+ const VectorBase &x ,
+ const VectorSet &p
+ );
+ template
+ size_t SparseJacobianForward(
+ const VectorBase& x ,
+ const VectorSet& p ,
+ const VectorSize& r ,
+ const VectorSize& c ,
+ VectorBase& jac ,
+ sparse_jacobian_work& work
+ );
+ template
+ size_t SparseJacobianReverse(
+ const VectorBase& x ,
+ const VectorSet& p ,
+ const VectorSize& r ,
+ const VectorSize& c ,
+ VectorBase& jac ,
+ sparse_jacobian_work& work
+ );
+
+ /// calculate sparse Hessians
+ template
+ VectorBase SparseHessian(
+ const VectorBase& x ,
+ const VectorBase& w
+ );
+ template
+ VectorBase SparseHessian(
+ const VectorBase& x ,
+ const VectorBase& w ,
+ const VectorBool& p
+ );
+ template
+ size_t SparseHessian(
+ const VectorBase& x ,
+ const VectorBase& w ,
+ const VectorSet& p ,
+ const VectorSize& r ,
+ const VectorSize& c ,
+ VectorBase& hes ,
+ sparse_hessian_work& work
+ );
+
+ // Optimize the tape
+ // (see doxygen documentation in optimize.hpp)
+ void optimize( const std::string& options = "" );
+
+ // create abs-normal representation of the function f(x)
+ void abs_normal_fun( ADFun& g, ADFun& a );
+ // ------------------- Deprecated -----------------------------
+
+ /// deprecated: assign a new operation sequence
+ template
+ void Dependent(const ADvector &y);
+
+ /// Deprecated: number of variables in opertion sequence
+ size_t Size(void) const
+ { return num_var_tape_; }
+
+ /// Deprecated: # taylor_ coefficients currently stored
+ /// (per variable,direction)
+ size_t Order(void) const
+ { return num_order_taylor_ - 1; }
+
+ /// Deprecated: amount of memory for this object
+ /// Note that an approximation is used for the std::set memory
+ size_t Memory(void) const
+ { size_t pervar = cap_order_taylor_ * sizeof(Base)
+ + for_jac_sparse_pack_.memory()
+ + for_jac_sparse_set_.memory();
+ size_t total = num_var_tape_ * pervar + play_.Memory();
+ return total;
+ }
+
+ /// Deprecated: # taylor_ coefficient orderss stored
+ /// (per variable,direction)
+ size_t taylor_size(void) const
+ { return num_order_taylor_; }
+
+ /// Deprecated: Does this AD operation sequence use
+ /// VecAD::reference operands
+ bool use_VecAD(void) const
+ { return play_.num_vec_ind_rec() > 0; }
+
+ /// Deprecated: # taylor_ coefficient orders calculated
+ /// (per variable,direction)
+ size_t size_taylor(void) const
+ { return num_order_taylor_; }
+
+ /// Deprecated: set number of orders currently allocated
+ /// (per variable,direction)
+ void capacity_taylor(size_t per_var);
+};
+// ---------------------------------------------------------------------------
+
+} // END_CPPAD_NAMESPACE
+
+// non-user interfaces
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+
+// user interfaces
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+
+# endif
diff --git a/external/cppad/include/cppad/core/ad_io.hpp b/external/cppad/include/cppad/core/ad_io.hpp
new file mode 100644
index 0000000000..7b8a71b6a2
--- /dev/null
+++ b/external/cppad/include/cppad/core/ad_io.hpp
@@ -0,0 +1,223 @@
+# ifndef CPPAD_CORE_AD_IO_HPP
+# define CPPAD_CORE_AD_IO_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+$begin ad_input$$
+$spell
+ VecAD
+ std
+ istream
+ const
+$$
+
+
+$section AD Output Stream Operator$$
+$mindex >> input write$$
+
+$head Syntax$$
+$icode%is% >> %x%$$
+
+
+$head Purpose$$
+Sets $icode x$$ to a $cref/parameter/glossary/Parameter/$$
+with value $icode b$$ corresponding to
+$codei%
+ %is% >> %b%
+%$$
+where $icode b$$ is a $icode Base$$ object.
+It is assumed that this $icode Base$$ input operation returns
+a reference to $icode is$$.
+
+$head is$$
+The operand $icode is$$ has prototype
+$codei%
+ std::istream& %is%
+%$$
+
+$head x$$
+The operand $icode x$$ has one of the following prototypes
+$codei%
+ AD<%Base%>& %x%
+%$$
+
+$head Result$$
+The result of this operation can be used as a reference to $icode is$$.
+For example, if the operand $icode y$$ has prototype
+$codei%
+ AD<%Base%> %y%
+%$$
+then the syntax
+$codei%
+ %is% >> %x% >> %y%
+%$$
+will first read the $icode Base$$ value of $icode x$$ from $icode is$$,
+and then read the $icode Base$$ value to $icode y$$.
+
+$head Operation Sequence$$
+The result of this operation is not an
+$cref/AD of Base/glossary/AD of Base/$$ object.
+Thus it will not be recorded as part of an
+AD of $icode Base$$
+$cref/operation sequence/glossary/Operation/Sequence/$$.
+
+$head Example$$
+$children%
+ example/general/ad_input.cpp
+%$$
+The file
+$cref ad_input.cpp$$
+contains an example and test of this operation.
+It returns true if it succeeds and false otherwise.
+
+$end
+------------------------------------------------------------------------------
+$begin ad_output$$
+$spell
+ VecAD
+ std
+ ostream
+ const
+$$
+
+
+$section AD Output Stream Operator$$
+$mindex <<$$
+
+$head Syntax$$
+$icode%os% << %x%$$
+
+
+$head Purpose$$
+Writes the $icode Base$$ value, corresponding to $icode x$$,
+to the output stream $icode os$$.
+
+$head Assumption$$
+If $icode b$$ is a $icode Base$$ object,
+$codei%
+ %os% << %b%
+%$$
+returns a reference to $icode os$$.
+
+$head os$$
+The operand $icode os$$ has prototype
+$codei%
+ std::ostream& %os%
+%$$
+
+$head x$$
+The operand $icode x$$ has one of the following prototypes
+$codei%
+ const AD<%Base%>& %x%
+ const VecAD<%Base%>::reference& %x%
+%$$
+
+$head Result$$
+The result of this operation can be used as a reference to $icode os$$.
+For example, if the operand $icode y$$ has prototype
+$codei%
+ AD<%Base%> %y%
+%$$
+then the syntax
+$codei%
+ %os% << %x% << %y%
+%$$
+will output the value corresponding to $icode x$$
+followed by the value corresponding to $icode y$$.
+
+$head Operation Sequence$$
+The result of this operation is not an
+$cref/AD of Base/glossary/AD of Base/$$ object.
+Thus it will not be recorded as part of an
+AD of $icode Base$$
+$cref/operation sequence/glossary/Operation/Sequence/$$.
+
+$head Example$$
+$children%
+ example/general/ad_output.cpp
+%$$
+The file
+$cref ad_output.cpp$$
+contains an example and test of this operation.
+It returns true if it succeeds and false otherwise.
+
+$end
+------------------------------------------------------------------------------
+*/
+namespace CppAD { // BEGIN_CPPAD_NAMESPACE
+/*!
+\file ad_io.hpp
+AD input and ouput stream operators.
+*/
+// ---------------------------------------------------------------------------
+/*!
+Read an AD object from an input stream.
+
+\tparam Base
+Base type for the AD object.
+
+\param is [in,out]
+Is the input stream from which that value is read.
+
+\param x [out]
+is the object that is being set to a value.
+Upone return, x.value_ is read from the input stream
+and x.tape_is_ is zero; i.e., x is a parameter.
+*/
+template
+CPPAD_INLINE_FRIEND_TEMPLATE_FUNCTION
+std::istream& operator >> (std::istream& is, AD& x)
+{ // like assignment to a base type value
+ x.tape_id_ = 0;
+ CPPAD_ASSERT_UNKNOWN( Parameter(x) );
+ return (is >> x.value_);
+}
+// ---------------------------------------------------------------------------
+/*!
+Write an AD object to an output stream.
+
+\tparam Base
+Base type for the AD object.
+
+\param os [in,out]
+Is the output stream to which that value is written.
+
+\param x
+is the object that is being written to the output stream.
+This is equivalent to writing x.value_ to the output stream.
+*/
+template
+CPPAD_INLINE_FRIEND_TEMPLATE_FUNCTION
+std::ostream& operator << (std::ostream &os, const AD &x)
+{ return (os << x.value_); }
+// ---------------------------------------------------------------------------
+/*!
+Write a VecAD_reference object to an output stream.
+
+\tparam Base
+Base type for the VecAD_reference object.
+
+\param os [in,out]
+Is the output stream to which that value is written.
+
+\param x
+is the element of the VecAD object that is being written to the output stream.
+This is equivalent to writing the corresponing Base value to the stream.
+*/
+template
+CPPAD_INLINE_FRIEND_TEMPLATE_FUNCTION
+std::ostream& operator << (std::ostream &os, const VecAD_reference &x)
+{ return (os << x.ADBase()); }
+
+} // END_CPPAD_NAMESPACE
+# endif
diff --git a/external/cppad/include/cppad/core/ad_to_string.hpp b/external/cppad/include/cppad/core/ad_to_string.hpp
new file mode 100644
index 0000000000..dc431b1c0e
--- /dev/null
+++ b/external/cppad/include/cppad/core/ad_to_string.hpp
@@ -0,0 +1,71 @@
+// $Id$
+# ifndef CPPAD_CORE_AD_TO_STRING_HPP
+# define CPPAD_CORE_AD_TO_STRING_HPP
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+$begin ad_to_string$$
+$spell
+ const
+ std
+$$
+
+$section Convert An AD or Base Type to String$$
+
+$head Syntax$$
+$icode%s% = to_string(%value%)%$$.
+
+$head See Also$$
+$cref to_string$$, $cref base_to_string$$
+
+$head value$$
+The argument $icode value$$ has prototype
+$codei%
+ const AD<%Base%>& %value%
+ const %Base%& %value%
+%$$
+where $icode Base$$ is a type that supports the
+$cref base_to_string$$ type requirement.
+
+$head s$$
+The return value has prototype
+$codei%
+ std::string %s%
+%$$
+and contains a representation of the specified $icode value$$.
+If $icode value$$ is an AD type,
+the result has the same precision as for the $icode Base$$ type.
+
+$head Example$$
+The file $cref to_string.cpp$$
+includes an example and test of $code to_string$$ with AD types.
+It returns true if it succeeds and false otherwise.
+
+$end
+*/
+# include
+# include
+
+namespace CppAD {
+
+ // Template definition is in cppad/utility/to_string.hpp.
+ // Partial specialzation for AD types
+ template
+ struct to_string_struct< CppAD::AD >
+ { std::string operator()(const CppAD::AD& value)
+ { to_string_struct ts;
+ return ts( Value( Var2Par( value ) ) ); }
+ };
+
+}
+
+# endif
diff --git a/external/cppad/include/cppad/core/ad_valued.hpp b/external/cppad/include/cppad/core/ad_valued.hpp
new file mode 100644
index 0000000000..810703822b
--- /dev/null
+++ b/external/cppad/include/cppad/core/ad_valued.hpp
@@ -0,0 +1,49 @@
+// $Id$
+# ifndef CPPAD_CORE_AD_VALUED_HPP
+# define CPPAD_CORE_AD_VALUED_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+$begin ADValued$$
+$spell
+$$
+
+
+$section AD Valued Operations and Functions$$
+
+$comment atomic.omh includes atomic_base.omh which atomic_base.hpp$$
+$childtable%
+ cppad/core/arithmetic.hpp%
+ cppad/core/standard_math.hpp%
+ cppad/core/cond_exp.hpp%
+ cppad/core/discrete.hpp%
+ cppad/core/numeric_limits.hpp%
+ omh/atomic.omh
+%$$
+
+$end
+*/
+
+// include MathOther.h after CondExp.h because some MathOther.h routines use
+// CondExp.h and CondExp.h is not sufficently declared in Declare.h
+
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+# include
+
+# endif
diff --git a/external/cppad/include/cppad/core/add.hpp b/external/cppad/include/cppad/core/add.hpp
new file mode 100644
index 0000000000..ad15e36b16
--- /dev/null
+++ b/external/cppad/include/cppad/core/add.hpp
@@ -0,0 +1,97 @@
+// $Id$
+# ifndef CPPAD_CORE_ADD_HPP
+# define CPPAD_CORE_ADD_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+// BEGIN CppAD namespace
+namespace CppAD {
+
+template
+AD operator + (const AD &left , const AD &right)
+{
+ // compute the Base part of this AD object
+ AD result;
+ result.value_ = left.value_ + right.value_;
+ CPPAD_ASSERT_UNKNOWN( Parameter(result) );
+
+ // check if there is a recording in progress
+ local::ADTape* tape = AD::tape_ptr();
+ if( tape == CPPAD_NULL )
+ return result;
+ tape_id_t tape_id = tape->id_;
+
+ // tape_id cannot match the default value for tape_id_; i.e., 0
+ CPPAD_ASSERT_UNKNOWN( tape_id > 0 );
+ bool var_left = left.tape_id_ == tape_id;
+ bool var_right = right.tape_id_ == tape_id;
+
+ if( var_left )
+ { if( var_right )
+ { // result = variable + variable
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::AddvvOp) == 1 );
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::AddvvOp) == 2 );
+
+ // put operand addresses in tape
+ tape->Rec_.PutArg(left.taddr_, right.taddr_);
+ // put operator in the tape
+ result.taddr_ = tape->Rec_.PutOp(local::AddvvOp);
+ // make result a variable
+ result.tape_id_ = tape_id;
+ }
+ else if( IdenticalZero(right.value_) )
+ { // result = variable + 0
+ result.make_variable(left.tape_id_, left.taddr_);
+ }
+ else
+ { // result = variable + parameter
+ // = parameter + variable
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::AddpvOp) == 1 );
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::AddpvOp) == 2 );
+
+ // put operand addresses in tape
+ addr_t p = tape->Rec_.PutPar(right.value_);
+ tape->Rec_.PutArg(p, left.taddr_);
+ // put operator in the tape
+ result.taddr_ = tape->Rec_.PutOp(local::AddpvOp);
+ // make result a variable
+ result.tape_id_ = tape_id;
+ }
+ }
+ else if( var_right )
+ { if( IdenticalZero(left.value_) )
+ { // result = 0 + variable
+ result.make_variable(right.tape_id_, right.taddr_);
+ }
+ else
+ { // result = parameter + variable
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::AddpvOp) == 1 );
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::AddpvOp) == 2 );
+
+ // put operand addresses in tape
+ addr_t p = tape->Rec_.PutPar(left.value_);
+ tape->Rec_.PutArg(p, right.taddr_);
+ // put operator in the tape
+ result.taddr_ = tape->Rec_.PutOp(local::AddpvOp);
+ // make result a variable
+ result.tape_id_ = tape_id;
+ }
+ }
+ return result;
+}
+
+// convert other cases into the case above
+CPPAD_FOLD_AD_VALUED_BINARY_OPERATOR(+)
+
+} // END CppAD namespace
+
+# endif
diff --git a/external/cppad/include/cppad/core/add_eq.hpp b/external/cppad/include/cppad/core/add_eq.hpp
new file mode 100644
index 0000000000..dc22c0969f
--- /dev/null
+++ b/external/cppad/include/cppad/core/add_eq.hpp
@@ -0,0 +1,92 @@
+// $Id$
+# ifndef CPPAD_CORE_ADD_EQ_HPP
+# define CPPAD_CORE_ADD_EQ_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-16 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+// BEGIN CppAD namespace
+namespace CppAD {
+
+template
+AD& AD::operator += (const AD &right)
+{
+ // compute the Base part
+ Base left;
+ left = value_;
+ value_ += right.value_;
+
+ // check if there is a recording in progress
+ local::ADTape* tape = AD::tape_ptr();
+ if( tape == CPPAD_NULL )
+ return *this;
+ tape_id_t tape_id = tape->id_;
+
+ // tape_id cannot match the default value for tape_id_; i.e., 0
+ CPPAD_ASSERT_UNKNOWN( tape_id > 0 );
+ bool var_left = tape_id_ == tape_id;
+ bool var_right = right.tape_id_ == tape_id;
+
+ if( var_left )
+ { if( var_right )
+ { // this = variable + variable
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::AddvvOp) == 1 );
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::AddvvOp) == 2 );
+
+ // put operand addresses in tape
+ tape->Rec_.PutArg(taddr_, right.taddr_);
+ // put operator in the tape
+ taddr_ = tape->Rec_.PutOp(local::AddvvOp);
+ // make this a variable
+ CPPAD_ASSERT_UNKNOWN( tape_id_ == tape_id );
+ }
+ else if( ! IdenticalZero( right.value_ ) )
+ { // this = variable + parameter
+ // = parameter + variable
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::AddpvOp) == 1 );
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::AddpvOp) == 2 );
+
+ // put operand addresses in tape
+ addr_t p = tape->Rec_.PutPar(right.value_);
+ tape->Rec_.PutArg(p, taddr_);
+ // put operator in the tape
+ taddr_ = tape->Rec_.PutOp(local::AddpvOp);
+ // make this a variable
+ CPPAD_ASSERT_UNKNOWN( tape_id_ == tape_id );
+ }
+ }
+ else if( var_right )
+ { if( IdenticalZero(left) )
+ { // this = 0 + right
+ make_variable(right.tape_id_, right.taddr_);
+ }
+ else
+ { // this = parameter + variable
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::AddpvOp) == 1 );
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::AddpvOp) == 2 );
+
+ // put operand addresses in tape
+ addr_t p = tape->Rec_.PutPar(left);
+ tape->Rec_.PutArg(p, right.taddr_);
+ // put operator in the tape
+ taddr_ = tape->Rec_.PutOp(local::AddpvOp);
+ // make this a variable
+ tape_id_ = tape_id;
+ }
+ }
+ return *this;
+}
+
+CPPAD_FOLD_ASSIGNMENT_OPERATOR(+=)
+
+} // END CppAD namespace
+
+# endif
diff --git a/external/cppad/include/cppad/core/arithmetic.hpp b/external/cppad/include/cppad/core/arithmetic.hpp
new file mode 100644
index 0000000000..e1b3fe3dee
--- /dev/null
+++ b/external/cppad/include/cppad/core/arithmetic.hpp
@@ -0,0 +1,42 @@
+# ifndef CPPAD_CORE_ARITHMETIC_HPP
+# define CPPAD_CORE_ARITHMETIC_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+-------------------------------------------------------------------------------
+$begin Arithmetic$$
+$spell
+ Op
+ const
+$$
+
+
+
+$section AD Arithmetic Operators and Compound Assignments$$
+
+$childtable%
+ cppad/core/unary_plus.hpp%
+ cppad/core/unary_minus.hpp%
+ cppad/core/ad_binary.hpp%
+ cppad/core/compound_assign.hpp
+%$$
+
+$end
+-------------------------------------------------------------------------------
+*/
+# include
+# include
+# include
+# include
+
+# endif
diff --git a/external/cppad/include/cppad/core/asinh.hpp b/external/cppad/include/cppad/core/asinh.hpp
new file mode 100644
index 0000000000..fdebe63bda
--- /dev/null
+++ b/external/cppad/include/cppad/core/asinh.hpp
@@ -0,0 +1,96 @@
+# ifndef CPPAD_CORE_ASINH_HPP
+# define CPPAD_CORE_ASINH_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+-------------------------------------------------------------------------------
+
+$begin asinh$$
+$spell
+ asinh
+ const
+ Vec
+ std
+ cmath
+ CppAD
+$$
+$section The Inverse Hyperbolic Sine Function: asinh$$
+
+$head Syntax$$
+$icode%y% = asinh(%x%)%$$
+
+$head Description$$
+The inverse hyperbolic sine function is defined by
+$icode%x% == sinh(%y%)%$$.
+
+$head x, y$$
+See the $cref/possible types/unary_standard_math/Possible Types/$$
+for a unary standard math function.
+
+$head CPPAD_USE_CPLUSPLUS_2011$$
+
+$subhead true$$
+If this preprocessor symbol is true ($code 1$$),
+and $icode x$$ is an AD type,
+this is an $cref/atomic operation/glossary/Operation/Atomic/$$.
+
+$subhead false$$
+If this preprocessor symbol is false ($code 0$$),
+CppAD uses the representation
+$latex \[
+\R{asinh} (x) = \log \left( x + \sqrt{ 1 + x^2 } \right)
+\] $$
+to compute this function.
+
+$head Example$$
+$children%
+ example/general/asinh.cpp
+%$$
+The file
+$cref asinh.cpp$$
+contains an example and test of this function.
+It returns true if it succeeds and false otherwise.
+
+$end
+-------------------------------------------------------------------------------
+*/
+# include
+# if ! CPPAD_USE_CPLUSPLUS_2011
+
+// BEGIN CppAD namespace
+namespace CppAD {
+
+template
+Type asinh_template(const Type &x)
+{ return CppAD::log( x + CppAD::sqrt( Type(1) + x * x ) );
+}
+
+inline float asinh(const float &x)
+{ return asinh_template(x); }
+
+inline double asinh(const double &x)
+{ return asinh_template(x); }
+
+template
+inline AD asinh(const AD &x)
+{ return asinh_template(x); }
+
+template
+inline AD asinh(const VecAD_reference &x)
+{ return asinh_template( x.ADBase() ); }
+
+
+} // END CppAD namespace
+
+# endif // CPPAD_USE_CPLUSPLUS_2011
+# endif // CPPAD_ASINH_INCLUDED
diff --git a/external/cppad/include/cppad/core/atan2.hpp b/external/cppad/include/cppad/core/atan2.hpp
new file mode 100644
index 0000000000..103e7313cb
--- /dev/null
+++ b/external/cppad/include/cppad/core/atan2.hpp
@@ -0,0 +1,141 @@
+# ifndef CPPAD_CORE_ATAN2_HPP
+# define CPPAD_CORE_ATAN2_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+-------------------------------------------------------------------------------
+$begin atan2$$
+$spell
+ Vec
+ CppAD
+ namespace
+ std
+ atan
+ const
+$$
+
+
+$section AD Two Argument Inverse Tangent Function$$
+$mindex tan atan2$$
+
+$head Syntax$$
+$icode%theta% = atan2(%y%, %x%)%$$
+
+
+$head Purpose$$
+Determines an angle $latex \theta \in [ - \pi , + \pi ]$$
+such that
+$latex \[
+\begin{array}{rcl}
+ \sin ( \theta ) & = & y / \sqrt{ x^2 + y^2 } \\
+ \cos ( \theta ) & = & x / \sqrt{ x^2 + y^2 }
+\end{array}
+\] $$
+
+$head y$$
+The argument $icode y$$ has one of the following prototypes
+$codei%
+ const AD<%Base%> &%y%
+ const VecAD<%Base%>::reference &%y%
+%$$
+
+$head x$$
+The argument $icode x$$ has one of the following prototypes
+$codei%
+ const AD<%Base%> &%x%
+ const VecAD<%Base%>::reference &%x%
+%$$
+
+$head theta$$
+The result $icode theta$$ has prototype
+$codei%
+ AD<%Base%> %theta%
+%$$
+
+$head Operation Sequence$$
+The AD of $icode Base$$
+operation sequence used to calculate $icode theta$$ is
+$cref/independent/glossary/Operation/Independent/$$
+of $icode x$$ and $icode y$$.
+
+$head Example$$
+$children%
+ example/general/atan2.cpp
+%$$
+The file
+$cref atan2.cpp$$
+contains an example and test of this function.
+It returns true if it succeeds and false otherwise.
+
+$end
+-------------------------------------------------------------------------------
+*/
+
+namespace CppAD { // BEGIN CppAD namespace
+
+inline float atan2(float x, float y)
+{ return std::atan2(x, y); }
+
+inline double atan2(double x, double y)
+{ return std::atan2(x, y); }
+
+// The code below is used as an example by the CondExp documentation.
+// BEGIN CondExp
+template
+AD atan2 (const AD &y, const AD &x)
+{ AD alpha;
+ AD beta;
+ AD theta;
+
+ AD zero(0.);
+ AD pi2(2. * atan(1.));
+ AD pi(2. * pi2);
+
+ AD ax = fabs(x);
+ AD ay = fabs(y);
+
+ // if( ax > ay )
+ // theta = atan(ay / ax);
+ // else theta = pi2 - atan(ax / ay);
+ alpha = atan(ay / ax);
+ beta = pi2 - atan(ax / ay);
+ theta = CondExpGt(ax, ay, alpha, beta); // use of CondExp
+
+ // if( x <= 0 )
+ // theta = pi - theta;
+ theta = CondExpLe(x, zero, pi - theta, theta); // use of CondExp
+
+ // if( y <= 0 )
+ // theta = - theta;
+ theta = CondExpLe(y, zero, -theta, theta); // use of CondExp
+
+ return theta;
+}
+// END CondExp
+
+template
+inline AD atan2 (const VecAD_reference &y, const AD &x)
+{ return atan2( y.ADBase() , x ); }
+
+template
+inline AD atan2 (const AD &y, const VecAD_reference &x)
+{ return atan2( y , x.ADBase() ); }
+
+template
+inline AD atan2
+(const VecAD_reference &y, const VecAD_reference &x)
+{ return atan2( y.ADBase() , x.ADBase() ); }
+
+} // END CppAD namespace
+
+# endif
diff --git a/external/cppad/include/cppad/core/atanh.hpp b/external/cppad/include/cppad/core/atanh.hpp
new file mode 100644
index 0000000000..9880c1f23b
--- /dev/null
+++ b/external/cppad/include/cppad/core/atanh.hpp
@@ -0,0 +1,96 @@
+# ifndef CPPAD_CORE_ATANH_HPP
+# define CPPAD_CORE_ATANH_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+/*
+-------------------------------------------------------------------------------
+$begin atanh$$
+$spell
+ atanh
+ const
+ Vec
+ std
+ cmath
+ CppAD
+ tanh
+$$
+$section The Inverse Hyperbolic Tangent Function: atanh$$
+
+$head Syntax$$
+$icode%y% = atanh(%x%)%$$
+
+$head Description$$
+The inverse hyperbolic tangent function is defined by
+$icode%x% == tanh(%y%)%$$.
+
+$head x, y$$
+See the $cref/possible types/unary_standard_math/Possible Types/$$
+for a unary standard math function.
+
+$head CPPAD_USE_CPLUSPLUS_2011$$
+
+$subhead true$$
+If this preprocessor symbol is true ($code 1$$),
+and $icode x$$ is an AD type,
+this is an $cref/atomic operation/glossary/Operation/Atomic/$$.
+
+$subhead false$$
+If this preprocessor symbol is false ($code 0$$),
+CppAD uses the representation
+$latex \[
+\R{atanh} (x) = \frac{1}{2} \log \left( \frac{1 + x}{1 - x} \right)
+\] $$
+to compute this function.
+
+$head Example$$
+$children%
+ example/general/atanh.cpp
+%$$
+The file
+$cref atanh.cpp$$
+contains an example and test of this function.
+It returns true if it succeeds and false otherwise.
+
+$end
+-------------------------------------------------------------------------------
+*/
+# include
+# if ! CPPAD_USE_CPLUSPLUS_2011
+
+// BEGIN CppAD namespace
+namespace CppAD {
+
+template
+Type atanh_template(const Type &x)
+{ return CppAD::log( (Type(1) + x) / (Type(1) - x) ) / Type(2);
+}
+
+inline float atanh(const float &x)
+{ return atanh_template(x); }
+
+inline double atanh(const double &x)
+{ return atanh_template(x); }
+
+template
+inline AD atanh(const AD &x)
+{ return atanh_template(x); }
+
+template
+inline AD atanh(const VecAD_reference &x)
+{ return atanh_template( x.ADBase() ); }
+
+
+} // END CppAD namespace
+
+# endif // CPPAD_USE_CPLUSPLUS_2011
+# endif // CPPAD_ATANH_INCLUDED
diff --git a/external/cppad/include/cppad/core/atomic_base.hpp b/external/cppad/include/cppad/core/atomic_base.hpp
new file mode 100644
index 0000000000..4aa8b7021c
--- /dev/null
+++ b/external/cppad/include/cppad/core/atomic_base.hpp
@@ -0,0 +1,2423 @@
+# ifndef CPPAD_CORE_ATOMIC_BASE_HPP
+# define CPPAD_CORE_ATOMIC_BASE_HPP
+
+/* --------------------------------------------------------------------------
+CppAD: C++ Algorithmic Differentiation: Copyright (C) 2003-17 Bradley M. Bell
+
+CppAD is distributed under multiple licenses. This distribution is under
+the terms of the
+ Eclipse Public License Version 1.0.
+
+A copy of this license is included in the COPYING file of this distribution.
+Please visit http://www.coin-or.org/CppAD/ for information on other licenses.
+-------------------------------------------------------------------------- */
+
+# include
+# include
+# include
+// needed before one can use CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL
+# include
+
+namespace CppAD { // BEGIN_CPPAD_NAMESPACE
+/*!
+\file atomic_base.hpp
+Base class for atomic user operations.
+*/
+
+template
+class atomic_base {
+// ===================================================================
+public:
+ enum option_enum {
+ pack_sparsity_enum ,
+ bool_sparsity_enum ,
+ set_sparsity_enum
+ };
+private:
+ // ------------------------------------------------------
+ // constants
+ //
+ /// index of this object in class_object
+ const size_t index_;
+
+ // -----------------------------------------------------
+ // variables
+ //
+ /// sparsity pattern this object is currently using
+ /// (set by constructor and option member functions)
+ option_enum sparsity_;
+
+ /// temporary work space used by member functions, declared here to avoid
+ // memory allocation/deallocation for each usage
+ struct work_struct {
+ vector vx;
+ vector vy;
+ vector tx;
+ vector ty;
+ //
+ vector bool_t;
+ //
+ vectorBool pack_h;
+ vectorBool pack_r;
+ vectorBool pack_s;
+ vectorBool pack_u;
+ //
+ vector bool_h;
+ vector bool_r;
+ vector bool_s;
+ vector bool_u;
+ //
+ vector< std::set > set_h;
+ vector< std::set > set_r;
+ vector< std::set > set_s;
+ vector< std::set > set_u;
+ };
+ // Use pointers, to avoid false sharing between threads.
+ // Not using: vector work_;
+ // so that deprecated atomic examples do not result in a memory leak.
+ work_struct* work_[CPPAD_MAX_NUM_THREADS];
+ // -----------------------------------------------------
+ // static member functions
+ //
+ /// List of all the object in this class
+ static std::vector& class_object(void)
+ { CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
+ static std::vector list_;
+ return list_;
+ }
+ /// List of names for each object in this class
+ static std::vector& class_name(void)
+ { CPPAD_ASSERT_FIRST_CALL_NOT_PARALLEL;
+ static std::vector list_;
+ return list_;
+ }
+ // =====================================================================
+public:
+ // -----------------------------------------------------
+ // member functions not in user API
+ //
+ /// current sparsity setting
+ option_enum sparsity(void) const
+ { return sparsity_; }
+
+ /// Name corresponding to a base_atomic object
+ const std::string& afun_name(void) const
+ { return class_name()[index_]; }
+/*
+$begin atomic_ctor$$
+$spell
+ enum
+ sq
+ std
+ afun
+ arg
+ CppAD
+ bool
+ ctor
+ const
+ mat_mul_xam.cpp
+ hpp
+$$
+
+$section Atomic Function Constructor$$
+
+$head Syntax$$
+$icode%atomic_user afun%(%ctor_arg_list%)
+%$$
+$codei%atomic_base<%Base%>(%name%, %sparsity%)
+%$$
+
+$head atomic_user$$
+
+$subhead ctor_arg_list$$
+Is a list of arguments for the $icode atomic_user$$ constructor.
+
+$subhead afun$$
+The object $icode afun$$ must stay in scope for as long
+as the corresponding atomic function is used.
+This includes use by any $cref/ADFun/ADFun/$$ that
+has this $icode atomic_user$$ operation in its
+$cref/operation sequence/glossary/Operation/Sequence/$$.
+
+$subhead Implementation$$
+The user defined $icode atomic_user$$ class is a publicly derived class of
+$codei%atomic_base<%Base%>%$$.
+It should be declared as follows:
+$codei%
+ class %atomic_user% : public CppAD::atomic_base<%Base%> {
+ public:
+ %atomic_user%(%ctor_arg_list%) : atomic_base<%Base%>(%name%, %sparsity%)
+ %...%
+ };
+%$$
+where $icode ...$$
+denotes the rest of the implementation of the derived class.
+This includes completing the constructor and
+all the virtual functions that have their
+$code atomic_base$$ implementations replaced by
+$icode atomic_user$$ implementations.
+
+$head atomic_base$$
+
+$subhead Restrictions$$
+The $code atomic_base$$ constructor cannot be called in
+$cref/parallel/ta_in_parallel/$$ mode.
+
+$subhead Base$$
+The template parameter determines the
+$icode Base$$ type for this $codei%AD<%Base%>%$$ atomic operation.
+
+$subhead name$$
+This $code atomic_base$$ constructor argument has the following prototype
+$codei%
+ const std::string& %name%
+%$$
+It is the name for this atomic function and is used for error reporting.
+The suggested value for $icode name$$ is $icode afun$$ or $icode atomic_user$$,
+i.e., the name of the corresponding atomic object or class.
+
+$subhead sparsity$$
+This $code atomic_base$$ constructor argument has prototype
+$codei%
+ atomic_base<%Base%>::option_enum %sparsity%
+%$$
+The current $icode sparsity$$ for an $code atomic_base$$ object
+determines which type of sparsity patterns it uses
+and its value is one of the following:
+$table
+$icode sparsity$$ $cnext sparsity patterns $rnext
+$codei%atomic_base<%Base%>::pack_sparsity_enum%$$ $pre $$ $cnext
+ $cref/vectorBool/CppAD_vector/vectorBool/$$
+$rnext
+$codei%atomic_base<%Base%>::bool_sparsity_enum%$$ $pre $$ $cnext
+ $cref/vector/CppAD_vector/$$$code $$
+$rnext
+$codei%atomic_base<%Base%>::set_sparsity_enum%$$ $pre $$ $cnext
+ $cref/vector/CppAD_vector/$$$code >$$
+$tend
+There is a default value for $icode sparsity$$ if it is not
+included in the constructor (which may be either the bool or set option).
+
+$head Example$$
+
+$subhead Define Constructor$$
+The following is an example of a user atomic function constructor definitions:
+$cref%get_started.cpp%atomic_get_started.cpp%Constructor%$$.
+
+$subhead Use Constructor$$
+The following is an example using a user atomic function constructor:
+$cref%get_started.cpp%atomic_get_started.cpp%Use Atomic Function%Constructor%$$.
+
+$end
+*/
+/*!
+Base class for atomic_user functions.
+
+\tparam Base
+This class is used for defining an AD atomic operation y = f(x).
+*/
+/// make sure user does not invoke the default constructor
+atomic_base(void)
+{ CPPAD_ASSERT_KNOWN(false,
+ "Attempt to use the atomic_base default constructor"
+ );
+}
+/*!
+Constructor
+
+\param name
+name used for error reporting
+
+\param sparsity [in]
+what type of sparsity patterns are computed by this function,
+bool_sparsity_enum or set_sparsity_enum. Default value is
+bool sparsity patterns.
+*/
+atomic_base(
+ const std::string& name,
+ option_enum sparsity = bool_sparsity_enum
+) :
+index_ ( class_object().size() ) ,
+sparsity_( sparsity )
+{ CPPAD_ASSERT_KNOWN(
+ ! thread_alloc::in_parallel() ,
+ "atomic_base: constructor cannot be called in parallel mode."
+ );
+ class_object().push_back(this);
+ class_name().push_back(name);
+ CPPAD_ASSERT_UNKNOWN( class_object().size() == class_name().size() );
+ //
+ // initialize work pointers as null;
+ for(size_t thread = 0; thread < CPPAD_MAX_NUM_THREADS; thread++)
+ work_[thread] = CPPAD_NULL;
+}
+/// destructor informs CppAD that this atomic function with this index
+/// has dropped out of scope by setting its pointer to null
+virtual ~atomic_base(void)
+{ CPPAD_ASSERT_UNKNOWN( class_object().size() > index_ );
+ // change object pointer to null, but leave name for error reporting
+ class_object()[index_] = CPPAD_NULL;
+ //
+ // free temporary work memory
+ for(size_t thread = 0; thread < CPPAD_MAX_NUM_THREADS; thread++)
+ free_work(thread);
+}
+/// allocates work_ for a specified thread
+void allocate_work(size_t thread)
+{ if( work_[thread] == CPPAD_NULL )
+ { // allocate the raw memory
+ size_t min_bytes = sizeof(work_struct);
+ size_t num_bytes;
+ void* v_ptr = thread_alloc::get_memory(min_bytes, num_bytes);
+ // save in work_
+ work_[thread] = reinterpret_cast( v_ptr );
+ // call constructor
+ new( work_[thread] ) work_struct;
+ }
+ return;
+}
+/// frees work_ for a specified thread
+void free_work(size_t thread)
+{ if( work_[thread] != CPPAD_NULL )
+ { // call destructor
+ work_[thread]->~work_struct();
+ // return memory to avialable pool for this thread
+ thread_alloc::return_memory( reinterpret_cast(work_[thread]) );
+ // mark this thread as not allocated
+ work_[thread] = CPPAD_NULL;
+ }
+ return;
+}
+/// atomic_base function object corresponding to a certain index
+static atomic_base* class_object(size_t index)
+{ CPPAD_ASSERT_UNKNOWN( class_object().size() > index );
+ return class_object()[index];
+}
+/// atomic_base function name corresponding to a certain index
+static const std::string& class_name(size_t index)
+{ CPPAD_ASSERT_UNKNOWN( class_name().size() > index );
+ return class_name()[index];
+}
+/*
+$begin atomic_option$$
+$spell
+ sq
+ enum
+ afun
+ bool
+ CppAD
+ std
+ typedef
+$$
+
+$section Set Atomic Function Options$$
+
+$head Syntax$$
+$icode%afun%.option(%option_value%)%$$
+These settings do not apply to individual $icode afun$$ calls,
+but rather all subsequent uses of the corresponding atomic operation
+in an $cref ADFun$$ object.
+
+$head atomic_sparsity$$
+Note that, if you use $cref optimize$$, these sparsity patterns are used
+to determine the $cref/dependency/dependency.cpp/$$ relationship between
+argument and result variables.
+
+$subhead pack_sparsity_enum$$
+If $icode option_value$$ is $codei%atomic_base<%Base%>::pack_sparsity_enum%$$,
+then the type used by $icode afun$$ for
+$cref/sparsity patterns/glossary/Sparsity Pattern/$$,
+(after the option is set) will be
+$codei%
+ typedef CppAD::vectorBool %atomic_sparsity%
+%$$
+If $icode r$$ is a sparsity pattern
+for a matrix $latex R \in B^{p \times q}$$:
+$icode%r%.size() == %p% * %q%$$.
+
+$subhead bool_sparsity_enum$$
+If $icode option_value$$ is $codei%atomic_base<%Base%>::bool_sparsity_enum%$$,
+then the type used by $icode afun$$ for
+$cref/sparsity patterns/glossary/Sparsity Pattern/$$,
+(after the option is set) will be
+$codei%
+ typedef CppAD::vector %atomic_sparsity%
+%$$
+If $icode r$$ is a sparsity pattern
+for a matrix $latex R \in B^{p \times q}$$:
+$icode%r%.size() == %p% * %q%$$.
+
+$subhead set_sparsity_enum$$
+If $icode option_value$$ is $icode%atomic_base<%Base%>::set_sparsity_enum%$$,
+then the type used by $icode afun$$ for
+$cref/sparsity patterns/glossary/Sparsity Pattern/$$,
+(after the option is set) will be
+$codei%
+ typedef CppAD::vector< std::set > %atomic_sparsity%
+%$$
+If $icode r$$ is a sparsity pattern
+for a matrix $latex R \in B^{p \times q}$$:
+$icode%r%.size() == %p%$$, and for $latex i = 0 , \ldots , p-1$$,
+the elements of $icode%r%[%i%]%$$ are between zero and $latex q-1$$ inclusive.
+
+$end
+*/
+void option(enum option_enum option_value)
+{ switch( option_value )
+ { case pack_sparsity_enum:
+ case bool_sparsity_enum:
+ case set_sparsity_enum:
+ sparsity_ = option_value;
+ break;
+
+ default:
+ CPPAD_ASSERT_KNOWN(
+ false,
+ "atoic_base::option: option_value is not valid"
+ );
+ }
+ return;
+}
+/*
+-----------------------------------------------------------------------------
+$begin atomic_afun$$
+
+$spell
+ sq
+ mul
+ afun
+ const
+ CppAD
+ mat_mul.cpp
+$$
+
+$section Using AD Version of Atomic Function$$
+
+$head Syntax$$
+$icode%afun%(%ax%, %ay%)%$$
+
+$head Purpose$$
+Given $icode ax$$,
+this call computes the corresponding value of $icode ay$$.
+If $codei%AD<%Base%>%$$ operations are being recorded,
+it enters the computation as an atomic operation in the recording;
+see $cref/start recording/Independent/Start Recording/$$.
+
+$head ADVector$$
+The type $icode ADVector$$ must be a
+$cref/simple vector class/SimpleVector/$$ with elements of type
+$codei%AD<%Base%>%$$; see $cref/Base/atomic_ctor/atomic_base/Base/$$.
+
+$head afun$$
+is a $cref/atomic_user/atomic_ctor/atomic_user/$$ object
+and this $icode afun$$ function call is implemented by the
+$cref/atomic_base/atomic_ctor/atomic_base/$$ class.
+
+$head ax$$
+This argument has prototype
+$codei%
+ const %ADVector%& %ax%
+%$$
+and size must be equal to $icode n$$.
+It specifies vector $latex x \in B^n$$
+at which an $codei%AD<%Base%>%$$ version of
+$latex y = f(x)$$ is to be evaluated; see
+$cref/Base/atomic_ctor/atomic_base/Base/$$.
+
+$head ay$$
+This argument has prototype
+$codei%
+ %ADVector%& %ay%
+%$$
+and size must be equal to $icode m$$.
+The input values of its elements
+are not specified (must not matter).
+Upon return, it is an $codei%AD<%Base%>%$$ version of
+$latex y = f(x)$$.
+
+$head Examples$$
+The following files contain example uses of
+the AD version of atomic functions during recording:
+$cref%get_started.cpp%atomic_get_started.cpp%Use Atomic Function%Recording%$$,
+$cref%norm_sq.cpp%atomic_norm_sq.cpp%Use Atomic Function%Recording%$$,
+$cref%reciprocal.cpp%atomic_reciprocal.cpp%Use Atomic Function%Recording%$$,
+$cref%tangent.cpp%atomic_tangent.cpp%Use Atomic Function%Recording%$$,
+$cref%mat_mul.cpp%atomic_mat_mul.cpp%Use Atomic Function%Recording%$$.
+
+$end
+-----------------------------------------------------------------------------
+*/
+/*!
+Implement the user call to afun(ax, ay) and old_atomic call to
+afun(ax, ay, id).
+
+\tparam ADVector
+A simple vector class with elements of type AD.
+
+\param id
+optional extra information vector that is just passed through by CppAD,
+and used by old_atomic derived class (not other derived classes).
+This is an extra parameter to the virtual callbacks for old_atomic;
+see the set_old member function.
+
+\param ax
+is the argument vector for this call,
+ax.size() determines the number of arguments.
+
+\param ay
+is the result vector for this call,
+ay.size() determines the number of results.
+*/
+template
+void operator()(
+ const ADVector& ax ,
+ ADVector& ay ,
+ size_t id = 0 )
+{ size_t i, j;
+ size_t n = ax.size();
+ size_t m = ay.size();
+# ifndef NDEBUG
+ bool ok;
+ std::string msg = "atomic_base: " + afun_name() + ".eval: ";
+ if( (n == 0) | (m == 0) )
+ { msg += "ax.size() or ay.size() is zero";
+ CPPAD_ASSERT_KNOWN(false, msg.c_str() );
+ }
+# endif
+ size_t thread = thread_alloc::thread_num();
+ allocate_work(thread);
+ vector & tx = work_[thread]->tx;
+ vector & ty = work_[thread]->ty;
+ vector & vx = work_[thread]->vx;
+ vector & vy = work_[thread]->vy;
+ //
+ if( vx.size() != n )
+ { vx.resize(n);
+ tx.resize(n);
+ }
+ if( vy.size() != m )
+ { vy.resize(m);
+ ty.resize(m);
+ }
+ //
+ // Determine tape corresponding to variables in ax
+ tape_id_t tape_id = 0;
+ local::ADTape* tape = CPPAD_NULL;
+ for(j = 0; j < n; j++)
+ { tx[j] = ax[j].value_;
+ vx[j] = Variable( ax[j] );
+ if( vx[j] )
+ {
+ if( tape_id == 0 )
+ { tape = ax[j].tape_this();
+ tape_id = ax[j].tape_id_;
+ CPPAD_ASSERT_UNKNOWN( tape != CPPAD_NULL );
+ }
+# ifndef NDEBUG
+ if( tape_id != ax[j].tape_id_ )
+ { msg += afun_name() +
+ ": ax contains variables from different threads.";
+ CPPAD_ASSERT_KNOWN(false, msg.c_str());
+ }
+# endif
+ }
+ }
+ // Use zero order forward mode to compute values
+ size_t p = 0, q = 0;
+ set_old(id);
+# ifdef NDEBUG
+ forward(p, q, vx, vy, tx, ty);
+# else
+ ok = forward(p, q, vx, vy, tx, ty);
+ if( ! ok )
+ { msg += afun_name() + ": ok is false for "
+ "zero order forward mode calculation.";
+ CPPAD_ASSERT_KNOWN(false, msg.c_str());
+ }
+# endif
+ bool record_operation = false;
+ for(i = 0; i < m; i++)
+ {
+ // pass back values
+ ay[i].value_ = ty[i];
+
+ // initialize entire vector parameters (not in tape)
+ ay[i].tape_id_ = 0;
+ ay[i].taddr_ = 0;
+
+ // we need to record this operation if
+ // any of the eleemnts of ay are variables,
+ record_operation |= vy[i];
+ }
+# ifndef NDEBUG
+ if( record_operation & (tape == CPPAD_NULL) )
+ { msg +=
+ "all elements of vx are false but vy contains a true element";
+ CPPAD_ASSERT_KNOWN(false, msg.c_str() );
+ }
+# endif
+ // if tape is not null, ay is on the tape
+ if( record_operation )
+ {
+ // Operator that marks beginning of this atomic operation
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::UserOp) == 0 );
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::UserOp) == 4 );
+ CPPAD_ASSERT_KNOWN( std::numeric_limits::max() >=
+ std::max( std::max( std::max(index_, id), n), m ),
+ "atomic_base: cppad_tape_addr_type maximum not large enough"
+ );
+ tape->Rec_.PutArg(addr_t(index_), addr_t(id), addr_t(n), addr_t(m));
+ tape->Rec_.PutOp(local::UserOp);
+
+ // Now put n operators, one for each element of argument vector
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::UsravOp) == 0 );
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::UsrapOp) == 0 );
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::UsravOp) == 1 );
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::UsrapOp) == 1 );
+ for(j = 0; j < n; j++)
+ { if( vx[j] )
+ { // information for an argument that is a variable
+ tape->Rec_.PutArg(ax[j].taddr_);
+ tape->Rec_.PutOp(local::UsravOp);
+ }
+ else
+ { // information for an argument that is parameter
+ addr_t par = tape->Rec_.PutPar(ax[j].value_);
+ tape->Rec_.PutArg(par);
+ tape->Rec_.PutOp(local::UsrapOp);
+ }
+ }
+
+ // Now put m operators, one for each element of result vector
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::UsrrpOp) == 1 );
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::UsrrpOp) == 0 );
+ CPPAD_ASSERT_UNKNOWN( local::NumArg(local::UsrrvOp) == 0 );
+ CPPAD_ASSERT_UNKNOWN( local::NumRes(local::UsrrvOp) == 1 );
+ for(i = 0; i < m; i++)
+ { if( vy[i] )
+ { ay[i].taddr_ = tape->Rec_.PutOp(local::UsrrvOp);
+ ay[i].tape_id_ = tape_id;
+ }
+ else
+ { addr_t par = tape->Rec_.PutPar(ay[i].value_);
+ tape->Rec_.PutArg(par);
+ tape->Rec_.PutOp(local::UsrrpOp);
+ }
+ }
+
+ // Put a duplicate UserOp at end of UserOp sequence
+ CPPAD_ASSERT_KNOWN( std::numeric_limits::max() >=
+ std::max( std::max( std::max(index_, id), n), m ),
+ "atomic_base: cppad_tape_addr_type maximum not large enough"
+ );
+ tape->Rec_.PutArg(addr_t(index_), addr_t(id), addr_t(n), addr_t(m));
+ tape->Rec_.PutOp(local::UserOp);
+ }
+ return;
+}
+/*
+-----------------------------------------------------------------------------
+$begin atomic_forward$$
+$spell
+ sq
+ mul.hpp
+ hes
+ afun
+ vx
+ vy
+ ty
+ Taylor
+ const
+ CppAD
+ bool
+$$
+
+$section Atomic Forward Mode$$
+$mindex callback virtual$$
+
+
+$head Syntax$$
+$icode%ok% = %afun%.forward(%p%, %q%, %vx%, %vy%, %tx%, %ty%)%$$
+
+$head Purpose$$
+This virtual function is used by $cref atomic_afun$$
+to evaluate function values.
+It is also used buy $cref/forward/Forward/$$
+to compute function vales and derivatives.
+
+$head Implementation$$
+This virtual function must be defined by the
+$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
+It can just return $icode%ok% == false%$$
+(and not compute anything) for values
+of $icode%q% > 0%$$ that are greater than those used by your
+$cref/forward/Forward/$$ mode calculations.
+
+$head p$$
+The argument $icode p$$ has prototype
+$codei%
+ size_t %p%
+%$$
+It specifies the lowest order Taylor coefficient that we are evaluating.
+During calls to $cref atomic_afun$$, $icode%p% == 0%$$.
+
+$head q$$
+The argument $icode q$$ has prototype
+$codei%
+ size_t %q%
+%$$
+It specifies the highest order Taylor coefficient that we are evaluating.
+During calls to $cref atomic_afun$$, $icode%q% == 0%$$.
+
+$head vx$$
+The $code forward$$ argument $icode vx$$ has prototype
+$codei%
+ const CppAD::vector& %vx%
+%$$
+The case $icode%vx%.size() > 0%$$ only occurs while evaluating a call to
+$cref atomic_afun$$.
+In this case,
+$icode%p% == %q% == 0%$$,
+$icode%vx%.size() == %n%$$, and
+for $latex j = 0 , \ldots , n-1$$,
+$icode%vx%[%j%]%$$ is true if and only if
+$icode%ax%[%j%]%$$ is a $cref/variable/glossary/Variable/$$
+in the corresponding call to
+$codei%
+ %afun%(%ax%, %ay%)
+%$$
+If $icode%vx%.size() == 0%$$,
+then $icode%vy%.size() == 0%$$ and neither of these vectors
+should be used.
+
+$head vy$$
+The $code forward$$ argument $icode vy$$ has prototype
+$codei%
+ CppAD::vector& %vy%
+%$$
+If $icode%vy%.size() == 0%$$, it should not be used.
+Otherwise,
+$icode%q% == 0%$$ and $icode%vy%.size() == %m%$$.
+The input values of the elements of $icode vy$$
+are not specified (must not matter).
+Upon return, for $latex j = 0 , \ldots , m-1$$,
+$icode%vy%[%i%]%$$ is true if and only if
+$icode%ay%[%i%]%$$ is a variable
+(CppAD uses $icode vy$$ to reduce the necessary computations).
+
+$head tx$$
+The argument $icode tx$$ has prototype
+$codei%
+ const CppAD::vector<%Base%>& %tx%
+%$$
+and $icode%tx%.size() == (%q%+1)*%n%$$.
+For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , q$$,
+we use the Taylor coefficient notation
+$latex \[
+\begin{array}{rcl}
+ x_j^k & = & tx [ j * ( q + 1 ) + k ]
+ \\
+ X_j (t) & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^q t^q
+\end{array}
+\] $$
+Note that superscripts represent an index for $latex x_j^k$$
+and an exponent for $latex t^k$$.
+Also note that the Taylor coefficients for $latex X(t)$$ correspond
+to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way:
+$latex \[
+ x_j^k = \frac{1}{ k ! } X_j^{(k)} (0)
+\] $$
+
+$head ty$$
+The argument $icode ty$$ has prototype
+$codei%
+ CppAD::vector<%Base%>& %ty%
+%$$
+and $icode%tx%.size() == (%q%+1)*%m%$$.
+Upon return,
+For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , q$$,
+$latex \[
+\begin{array}{rcl}
+ Y_i (t) & = & f_i [ X(t) ]
+ \\
+ Y_i (t) & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^q t^q + o ( t^q )
+ \\
+ ty [ i * ( q + 1 ) + k ] & = & y_i^k
+\end{array}
+\] $$
+where $latex o( t^q ) / t^q \rightarrow 0$$ as $latex t \rightarrow 0$$.
+Note that superscripts represent an index for $latex y_j^k$$
+and an exponent for $latex t^k$$.
+Also note that the Taylor coefficients for $latex Y(t)$$ correspond
+to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way:
+$latex \[
+ y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0)
+\] $$
+If $latex p > 0$$,
+for $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , p-1$$,
+the input of $icode ty$$ satisfies
+$latex \[
+ ty [ i * ( q + 1 ) + k ] = y_i^k
+\]$$
+and hence the corresponding elements need not be recalculated.
+
+$head ok$$
+If the required results are calculated, $icode ok$$ should be true.
+Otherwise, it should be false.
+
+$head Discussion$$
+For example, suppose that $icode%q% == 2%$$,
+and you know how to compute the function $latex f(x)$$,
+its first derivative $latex f^{(1)} (x)$$,
+and it component wise Hessian $latex f_i^{(2)} (x)$$.
+Then you can compute $icode ty$$ using the following formulas:
+$latex \[
+\begin{array}{rcl}
+y_i^0 & = & Y(0)
+ = f_i ( x^0 )
+\\
+y_i^1 & = & Y^{(1)} ( 0 )
+ = f_i^{(1)} ( x^0 ) X^{(1)} ( 0 )
+ = f_i^{(1)} ( x^0 ) x^1
+\\
+y_i^2
+& = & \frac{1}{2 !} Y^{(2)} (0)
+\\
+& = & \frac{1}{2} X^{(1)} (0)^\R{T} f_i^{(2)} ( x^0 ) X^{(1)} ( 0 )
+ + \frac{1}{2} f_i^{(1)} ( x^0 ) X^{(2)} ( 0 )
+\\
+& = & \frac{1}{2} (x^1)^\R{T} f_i^{(2)} ( x^0 ) x^1
+ + f_i^{(1)} ( x^0 ) x^2
+\end{array}
+\] $$
+For $latex i = 0 , \ldots , m-1$$, and $latex k = 0 , 1 , 2$$,
+$latex \[
+ ty [ i * (q + 1) + k ] = y_i^k
+\] $$
+
+$children%
+ example/atomic/forward.cpp
+%$$
+$head Examples$$
+The file $cref atomic_forward.cpp$$ contains an example and test
+that uses this routine.
+It returns true if the test passes and false if it fails.
+
+$end
+-----------------------------------------------------------------------------
+*/
+/*!
+Link from atomic_base to forward mode
+
+\param p [in]
+lowerest order for this forward mode calculation.
+
+\param q [in]
+highest order for this forward mode calculation.
+
+\param vx [in]
+if size not zero, which components of \c x are variables
+
+\param vy [out]
+if size not zero, which components of \c y are variables
+
+\param tx [in]
+Taylor coefficients corresponding to \c x for this calculation.
+
+\param ty [out]
+Taylor coefficient corresponding to \c y for this calculation
+
+See the forward mode in user's documentation for base_atomic
+*/
+virtual bool forward(
+ size_t p ,
+ size_t q ,
+ const vector& vx ,
+ vector& vy ,
+ const vector& tx ,
+ vector& ty )
+{ return false; }
+/*
+-----------------------------------------------------------------------------
+$begin atomic_reverse$$
+$spell
+ sq
+ mul.hpp
+ afun
+ ty
+ px
+ py
+ Taylor
+ const
+ CppAD
+$$
+
+$section Atomic Reverse Mode$$
+$spell
+ bool
+$$
+
+$head Syntax$$
+$icode%ok% = %afun%.reverse(%q%, %tx%, %ty%, %px%, %py%)%$$
+
+$head Purpose$$
+This function is used by $cref/reverse/Reverse/$$
+to compute derivatives.
+
+$head Implementation$$
+If you are using
+$cref/reverse/Reverse/$$ mode,
+this virtual function must be defined by the
+$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
+It can just return $icode%ok% == false%$$
+(and not compute anything) for values
+of $icode q$$ that are greater than those used by your
+$cref/reverse/Reverse/$$ mode calculations.
+
+$head q$$
+The argument $icode q$$ has prototype
+$codei%
+ size_t %q%
+%$$
+It specifies the highest order Taylor coefficient that
+computing the derivative of.
+
+$head tx$$
+The argument $icode tx$$ has prototype
+$codei%
+ const CppAD::vector<%Base%>& %tx%
+%$$
+and $icode%tx%.size() == (%q%+1)*%n%$$.
+For $latex j = 0 , \ldots , n-1$$ and $latex k = 0 , \ldots , q$$,
+we use the Taylor coefficient notation
+$latex \[
+\begin{array}{rcl}
+ x_j^k & = & tx [ j * ( q + 1 ) + k ]
+ \\
+ X_j (t) & = & x_j^0 + x_j^1 t^1 + \cdots + x_j^q t^q
+\end{array}
+\] $$
+Note that superscripts represent an index for $latex x_j^k$$
+and an exponent for $latex t^k$$.
+Also note that the Taylor coefficients for $latex X(t)$$ correspond
+to the derivatives of $latex X(t)$$ at $latex t = 0$$ in the following way:
+$latex \[
+ x_j^k = \frac{1}{ k ! } X_j^{(k)} (0)
+\] $$
+
+$head ty$$
+The argument $icode ty$$ has prototype
+$codei%
+ const CppAD::vector<%Base%>& %ty%
+%$$
+and $icode%tx%.size() == (%q%+1)*%m%$$.
+For $latex i = 0 , \ldots , m-1$$ and $latex k = 0 , \ldots , q$$,
+we use the Taylor coefficient notation
+$latex \[
+\begin{array}{rcl}
+ Y_i (t) & = & f_i [ X(t) ]
+ \\
+ Y_i (t) & = & y_i^0 + y_i^1 t^1 + \cdots + y_i^q t^q + o ( t^q )
+ \\
+ y_i^k & = & ty [ i * ( q + 1 ) + k ]
+\end{array}
+\] $$
+where $latex o( t^q ) / t^q \rightarrow 0$$ as $latex t \rightarrow 0$$.
+Note that superscripts represent an index for $latex y_j^k$$
+and an exponent for $latex t^k$$.
+Also note that the Taylor coefficients for $latex Y(t)$$ correspond
+to the derivatives of $latex Y(t)$$ at $latex t = 0$$ in the following way:
+$latex \[
+ y_j^k = \frac{1}{ k ! } Y_j^{(k)} (0)
+\] $$
+
+
+$head F$$
+We use the notation $latex \{ x_j^k \} \in B^{n \times (q+1)}$$ for
+$latex \[
+ \{ x_j^k \W{:} j = 0 , \ldots , n-1, k = 0 , \ldots , q \}
+\]$$
+We use the notation $latex \{ y_i^k \} \in B^{m \times (q+1)}$$ for
+$latex \[
+ \{ y_i^k \W{:} i = 0 , \ldots , m-1, k = 0 , \ldots , q \}
+\]$$
+We define the function
+$latex F : B^{n \times (q+1)} \rightarrow B^{m \times (q+1)}$$ by
+$latex \[
+ y_i^k = F_i^k [ \{ x_j^k \} ]
+\] $$
+Note that
+$latex \[
+ F_i^0 ( \{ x_j^k \} ) = f_i ( X(0) ) = f_i ( x^0 )
+\] $$
+We also note that
+$latex F_i^\ell ( \{ x_j^k \} )$$ is a function of
+$latex x^0 , \ldots , x^\ell$$
+and is determined by the derivatives of $latex f_i (x)$$
+up to order $latex \ell$$.
+
+
+$head G, H$$
+We use $latex G : B^{m \times (q+1)} \rightarrow B$$
+to denote an arbitrary scalar valued function of $latex \{ y_i^k \}$$.
+We use $latex H : B^{n \times (q+1)} \rightarrow B$$
+defined by
+$latex \[
+ H ( \{ x_j^k \} ) = G[ F( \{ x_j^k \} ) ]
+\] $$
+
+$head py$$
+The argument $icode py$$ has prototype
+$codei%
+ const CppAD::vector<%Base%>& %py%
+%$$
+and $icode%py%.size() == m * (%q%+1)%$$.
+For $latex i = 0 , \ldots , m-1$$, $latex k = 0 , \ldots , q$$,
+$latex \[
+ py[ i * (q + 1 ) + k ] = \partial G / \partial y_i^k
+\] $$
+
+$subhead px$$
+The $icode px$$ has prototype
+$codei%
+ CppAD::vector<%Base%>& %px%
+%$$
+and $icode%px%.size() == n * (%q%+1)%$$.
+The input values of the elements of $icode px$$
+are not specified (must not matter).
+Upon return,
+for $latex j = 0 , \ldots , n-1$$ and $latex \ell = 0 , \ldots , q$$,
+$latex \[
+\begin{array}{rcl}
+px [ j * (q + 1) + \ell ] & = & \partial H / \partial x_j^\ell
+\\
+& = &
+( \partial G / \partial \{ y_i^k \} ) \cdot
+ ( \partial \{ y_i^k \} / \partial x_j^\ell )
+\\
+& = &
+\sum_{k=0}^q
+\sum_{i=0}^{m-1}
+( \partial G / \partial y_i^k ) ( \partial y_i^k / \partial x_j^\ell )
+\\
+& = &
+\sum_{k=\ell}^q
+\sum_{i=0}^{m-1}
+py[ i * (q + 1 ) + k ] ( \partial F_i^k / \partial x_j^\ell )
+\end{array}
+\] $$
+Note that we have used the fact that for $latex k < \ell$$,
+$latex \partial F_i^k / \partial x_j^\ell = 0$$.
+
+$head ok$$
+The return value $icode ok$$ has prototype
+$codei%
+ bool %ok%
+%$$
+If it is $code true$$, the corresponding evaluation succeeded,
+otherwise it failed.
+
+$children%
+ example/atomic/reverse.cpp
+%$$
+$head Examples$$
+The file $cref atomic_forward.cpp$$ contains an example and test
+that uses this routine.
+It returns true if the test passes and false if it fails.
+
+$end
+-----------------------------------------------------------------------------
+*/
+/*!
+Link from reverse mode sweep to users routine.
+
+\param q [in]
+highest order for this reverse mode calculation.
+
+\param tx [in]
+Taylor coefficients corresponding to \c x for this calculation.
+
+\param ty [in]
+Taylor coefficient corresponding to \c y for this calculation
+
+\param px [out]
+Partials w.r.t. the \c x Taylor coefficients.
+
+\param py [in]
+Partials w.r.t. the \c y Taylor coefficients.
+
+See atomic_reverse mode use documentation
+*/
+virtual bool reverse(
+ size_t q ,
+ const vector& tx ,
+ const vector& ty ,
+ vector& px ,
+ const vector& py )
+{ return false; }
+/*
+-------------------------------------- ---------------------------------------
+$begin atomic_for_sparse_jac$$
+$spell
+ sq
+ mul.hpp
+ afun
+ Jacobian
+ jac
+ const
+ CppAD
+ std
+ bool
+ std
+$$
+
+$section Atomic Forward Jacobian Sparsity Patterns$$
+
+$head Syntax$$
+$icode%ok% = %afun%.for_sparse_jac(%q%, %r%, %s%, %x%)
+%$$
+
+$head Deprecated 2016-06-27$$
+$icode%ok% = %afun%.for_sparse_jac(%q%, %r%, %s%)
+%$$
+
+$head Purpose$$
+This function is used by $cref ForSparseJac$$ to compute
+Jacobian sparsity patterns.
+For a fixed matrix $latex R \in B^{n \times q}$$,
+the Jacobian of $latex f( x + R * u)$$ with respect to $latex u \in B^q$$ is
+$latex \[
+ S(x) = f^{(1)} (x) * R
+\] $$
+Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
+$code for_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$.
+
+$head Implementation$$
+If you are using
+$cref ForSparseJac$$,
+$cref ForSparseHes$$, or
+$cref RevSparseHes$$,
+one of the versions of this
+virtual function must be defined by the
+$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
+
+$subhead q$$
+The argument $icode q$$ has prototype
+$codei%
+ size_t %q%
+%$$
+It specifies the number of columns in
+$latex R \in B^{n \times q}$$ and the Jacobian
+$latex S(x) \in B^{m \times q}$$.
+
+$subhead r$$
+This argument has prototype
+$codei%
+ const %atomic_sparsity%& %r%
+%$$
+and is a $cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
+$latex R \in B^{n \times q}$$.
+
+$subhead s$$
+This argument has prototype
+$codei%
+ %atomic_sparsity%& %s%
+%$$
+The input values of its elements
+are not specified (must not matter).
+Upon return, $icode s$$ is a
+$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
+$latex S(x) \in B^{m \times q}$$.
+
+$subhead x$$
+$index deprecated$$
+The argument has prototype
+$codei%
+ const CppAD::vector<%Base%>& %x%
+%$$
+and size is equal to the $icode n$$.
+This is the $cref Value$$ value corresponding to the parameters in the
+vector $cref/ax/atomic_afun/ax/$$ (when the atomic function was called).
+To be specific, if
+$codei%
+ if( Parameter(%ax%[%i%]) == true )
+ %x%[%i%] = Value( %ax%[%i%] );
+ else
+ %x%[%i%] = CppAD::numeric_limits<%Base%>::quiet_NaN();
+%$$
+The version of this function with out the $icode x$$ argument is deprecated;
+i.e., you should include the argument even if you do not use it.
+
+$head ok$$
+The return value $icode ok$$ has prototype
+$codei%
+ bool %ok%
+%$$
+If it is $code true$$, the corresponding evaluation succeeded,
+otherwise it failed.
+
+$children%
+ example/atomic/for_sparse_jac.cpp
+%$$
+$head Examples$$
+The file $cref atomic_for_sparse_jac.cpp$$ contains an example and test
+that uses this routine.
+It returns true if the test passes and false if it fails.
+
+$end
+-----------------------------------------------------------------------------
+*/
+/*!
+Link, after case split, from for_jac_sweep to atomic_base.
+
+\param q
+is the column dimension for the Jacobian sparsity partterns.
+
+\param r
+is the Jacobian sparsity pattern for the argument vector x
+
+\param s
+is the Jacobian sparsity pattern for the result vector y
+
+\param x
+is the integer value for x arguments that are parameters.
+*/
+virtual bool for_sparse_jac(
+ size_t q ,
+ const vector< std::set >& r ,
+ vector< std::set >& s ,
+ const vector& x )
+{ return false; }
+virtual bool for_sparse_jac(
+ size_t q ,
+ const vector& r ,
+ vector& s ,
+ const vector& x )
+{ return false; }
+virtual bool for_sparse_jac(
+ size_t q ,
+ const vectorBool& r ,
+ vectorBool& s ,
+ const vector& x )
+{ return false; }
+// deprecated versions
+virtual bool for_sparse_jac(
+ size_t q ,
+ const vector< std::set >& r ,
+ vector< std::set >& s )
+{ return false; }
+virtual bool for_sparse_jac(
+ size_t q ,
+ const vector& r ,
+ vector& s )
+{ return false; }
+virtual bool for_sparse_jac(
+ size_t q ,
+ const vectorBool& r ,
+ vectorBool& s )
+{ return false; }
+
+/*!
+Link, before case split, from for_jac_sweep to atomic_base.
+
+\tparam InternalSparsity
+Is the used internaly for sparsity calculations; i.e.,
+sparse_pack or sparse_list.
+
+\param x
+is parameter arguments to the function, other components are nan.
+
+\param x_index
+is the variable index, on the tape, for the arguments to this function.
+This size of x_index is n, the number of arguments to this function.
+
+\param y_index
+is the variable index, on the tape, for the results for this function.
+This size of y_index is m, the number of results for this function.
+
+\param var_sparsity
+On input, for j = 0, ... , n-1, the sparsity pattern with index x_index[j],
+is the sparsity for the j-th argument to this atomic function.
+On input, for i = 0, ... , m-1, the sparsity pattern with index y_index[i],
+is empty. On output, it is the sparsity
+for the j-th result for this atomic function.
+*/
+template
+void for_sparse_jac(
+ const vector& x ,
+ const vector& x_index ,
+ const vector& y_index ,
+ InternalSparsity& var_sparsity )
+{
+ // intial results are empty during forward mode
+ size_t q = var_sparsity.end();
+ bool input_empty = true;
+ bool zero_empty = true;
+ bool transpose = false;
+ size_t m = y_index.size();
+ bool ok = false;
+ size_t thread = thread_alloc::thread_num();
+ allocate_work(thread);
+ //
+ std::string msg = ": atomic_base.for_sparse_jac: returned false";
+ if( sparsity_ == pack_sparsity_enum )
+ { vectorBool& pack_r ( work_[thread]->pack_r );
+ vectorBool& pack_s ( work_[thread]->pack_s );
+ local::get_internal_sparsity(
+ transpose, x_index, var_sparsity, pack_r
+ );
+ //
+ pack_s.resize(m * q );
+ ok = for_sparse_jac(q, pack_r, pack_s, x);
+ if( ! ok )
+ ok = for_sparse_jac(q, pack_r, pack_s);
+ if( ! ok )
+ { msg = afun_name() + msg + " sparsity = pack_sparsity_enum";
+ CPPAD_ASSERT_KNOWN(false, msg.c_str());
+ }
+ local::set_internal_sparsity(zero_empty, input_empty,
+ transpose, y_index, var_sparsity, pack_s
+ );
+ }
+ else if( sparsity_ == bool_sparsity_enum )
+ { vector& bool_r ( work_[thread]->bool_r );
+ vector& bool_s ( work_[thread]->bool_s );
+ local::get_internal_sparsity(
+ transpose, x_index, var_sparsity, bool_r
+ );
+ bool_s.resize(m * q );
+ ok = for_sparse_jac(q, bool_r, bool_s, x);
+ if( ! ok )
+ ok = for_sparse_jac(q, bool_r, bool_s);
+ if( ! ok )
+ { msg = afun_name() + msg + " sparsity = bool_sparsity_enum";
+ CPPAD_ASSERT_KNOWN(false, msg.c_str());
+ }
+ local::set_internal_sparsity(zero_empty, input_empty,
+ transpose, y_index, var_sparsity, bool_s
+ );
+ }
+ else
+ { CPPAD_ASSERT_UNKNOWN( sparsity_ == set_sparsity_enum );
+ vector< std::set >& set_r ( work_[thread]->set_r );
+ vector< std::set >& set_s ( work_[thread]->set_s );
+ local::get_internal_sparsity(
+ transpose, x_index, var_sparsity, set_r
+ );
+ //
+ set_s.resize(m);
+ ok = for_sparse_jac(q, set_r, set_s, x);
+ if( ! ok )
+ ok = for_sparse_jac(q, set_r, set_s);
+ if( ! ok )
+ { msg = afun_name() + msg + " sparsity = set_sparsity_enum";
+ CPPAD_ASSERT_KNOWN(false, msg.c_str());
+ }
+ local::set_internal_sparsity(zero_empty, input_empty,
+ transpose, y_index, var_sparsity, set_s
+ );
+ }
+ return;
+}
+/*
+-------------------------------------- ---------------------------------------
+$begin atomic_rev_sparse_jac$$
+$spell
+ sq
+ mul.hpp
+ rt
+ afun
+ Jacobian
+ jac
+ CppAD
+ std
+ bool
+ const
+ hes
+$$
+
+$section Atomic Reverse Jacobian Sparsity Patterns$$
+
+$head Syntax$$
+$icode%ok% = %afun%.rev_sparse_jac(%q%, %rt%, %st%, %x%)
+%$$
+
+$head Deprecated 2016-06-27$$
+$icode%ok% = %afun%.rev_sparse_jac(%q%, %rt%, %st%)
+%$$
+
+$head Purpose$$
+This function is used by
+$cref RevSparseJac$$ to compute
+Jacobian sparsity patterns.
+If you are using $cref RevSparseJac$$,
+one of the versions of this
+virtual function must be defined by the
+$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
+$pre
+
+$$
+For a fixed matrix $latex R \in B^{q \times m}$$,
+the Jacobian of $latex R * f( x )$$ with respect to $latex x \in B^n$$ is
+$latex \[
+ S(x) = R * f^{(1)} (x)
+\] $$
+Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for $latex R$$,
+$code rev_sparse_jac$$ computes a sparsity pattern for $latex S(x)$$.
+
+$head Implementation$$
+If you are using
+$cref RevSparseJac$$ or $cref ForSparseHes$$,
+this virtual function must be defined by the
+$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
+
+$subhead q$$
+The argument $icode q$$ has prototype
+$codei%
+ size_t %q%
+%$$
+It specifies the number of rows in
+$latex R \in B^{q \times m}$$ and the Jacobian
+$latex S(x) \in B^{q \times n}$$.
+
+$subhead rt$$
+This argument has prototype
+$codei%
+ const %atomic_sparsity%& %rt%
+%$$
+and is a
+$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
+$latex R^\R{T} \in B^{m \times q}$$.
+
+$subhead st$$
+This argument has prototype
+$codei%
+ %atomic_sparsity%& %st%
+%$$
+The input value of its elements
+are not specified (must not matter).
+Upon return, $icode s$$ is a
+$cref/atomic_sparsity/atomic_option/atomic_sparsity/$$ pattern for
+$latex S(x)^\R{T} \in B^{n \times q}$$.
+
+$subhead x$$
+$index deprecated$$
+The argument has prototype
+$codei%
+ const CppAD::vector<%Base%>& %x%
+%$$
+and size is equal to the $icode n$$.
+This is the $cref Value$$ corresponding to the parameters in the
+vector $cref/ax/atomic_afun/ax/$$ (when the atomic function was called).
+To be specific, if
+$codei%
+ if( Parameter(%ax%[%i%]) == true )
+ %x%[%i%] = Value( %ax%[%i%] );
+ else
+ %x%[%i%] = CppAD::numeric_limits<%Base%>::quiet_NaN();
+%$$
+The version of this function with out the $icode x$$ argument is deprecated;
+i.e., you should include the argument even if you do not use it.
+
+$head ok$$
+The return value $icode ok$$ has prototype
+$codei%
+ bool %ok%
+%$$
+If it is $code true$$, the corresponding evaluation succeeded,
+otherwise it failed.
+
+$children%
+ example/atomic/rev_sparse_jac.cpp
+%$$
+$head Examples$$
+The file $cref atomic_rev_sparse_jac.cpp$$ contains an example and test
+that uses this routine.
+It returns true if the test passes and false if it fails.
+
+$end
+-----------------------------------------------------------------------------
+*/
+/*!
+Link, after case split, from rev_jac_sweep to atomic_base
+
+\param q [in]
+is the row dimension for the Jacobian sparsity partterns
+
+\param rt [out]
+is the tansposed Jacobian sparsity pattern w.r.t to range variables y
+
+\param st [in]
+is the tansposed Jacobian sparsity pattern for the argument variables x
+
+\param x
+is the integer value for x arguments that are parameters.
+*/
+virtual bool rev_sparse_jac(
+ size_t q ,
+ const vector< std::set >& rt ,
+ vector< std::set >& st ,
+ const vector& x )
+{ return false; }
+virtual bool rev_sparse_jac(
+ size_t q ,
+ const vector& rt ,
+ vector& st ,
+ const vector& x )
+{ return false; }
+virtual bool rev_sparse_jac(
+ size_t q ,
+ const vectorBool& rt ,
+ vectorBool& st ,
+ const vector& x )
+{ return false; }
+// deprecated versions
+virtual bool rev_sparse_jac(
+ size_t q ,
+ const vector< std::set >& rt ,
+ vector< std::set >& st )
+{ return false; }
+virtual bool rev_sparse_jac(
+ size_t q ,
+ const vector& rt ,
+ vector& st )
+{ return false; }
+virtual bool rev_sparse_jac(
+ size_t q ,
+ const vectorBool& rt ,
+ vectorBool& st )
+{ return false; }
+
+/*!
+Link, before case split, from rev_jac_sweep to atomic_base.
+
+\tparam InternalSparsity
+Is the used internaly for sparsity calculations; i.e.,
+sparse_pack or sparse_list.
+
+\param x
+is parameter arguments to the function, other components are nan.
+
+\param x_index
+is the variable index, on the tape, for the arguments to this function.
+This size of x_index is n, the number of arguments to this function.
+
+\param y_index
+is the variable index, on the tape, for the results for this function.
+This size of y_index is m, the number of results for this function.
+
+\param var_sparsity
+On input, for i = 0, ... , m-1, the sparsity pattern with index y_index[i],
+is the sparsity for the i-th argument to this atomic function.
+On output, for j = 0, ... , n-1, the sparsity pattern with index x_index[j],
+the sparsity has been updated to remove y as a function of x.
+*/
+template
+void rev_sparse_jac(
+ const vector& x ,
+ const vector& x_index ,
+ const vector& y_index ,
+ InternalSparsity& var_sparsity )
+{
+ // initial results may be non-empty during reverse mode
+ size_t q = var_sparsity.end();
+ bool input_empty = false;
+ bool zero_empty = true;
+ bool transpose = false;
+ size_t n = x_index.size();
+ bool ok = false;
+ size_t thread = thread_alloc::thread_num();
+ allocate_work(thread);
+ //
+ std::string msg = ": atomic_base.rev_sparse_jac: returned false";
+ if( sparsity_ == pack_sparsity_enum )
+ { vectorBool& pack_rt ( work_[thread]->pack_r );
+ vectorBool& pack_st ( work_[thread]->pack_s );
+ local::get_internal_sparsity(
+ transpose, y_index, var_sparsity, pack_rt
+ );
+ //
+ pack_st.resize(n * q );
+ ok = rev_sparse_jac(q, pack_rt, pack_st, x);
+ if( ! ok )
+ ok = rev_sparse_jac(q, pack_rt, pack_st);
+ if( ! ok )
+ { msg = afun_name() + msg + " sparsity = pack_sparsity_enum";
+ CPPAD_ASSERT_KNOWN(false, msg.c_str());
+ }
+ local::set_internal_sparsity(zero_empty, input_empty,
+ transpose, x_index, var_sparsity, pack_st
+ );
+ }
+ else if( sparsity_ == bool_sparsity_enum )
+ { vector& bool_rt ( work_[thread]->bool_r );
+ vector& bool_st ( work_[thread]->bool_s );
+ local::get_internal_sparsity(
+ transpose, y_index, var_sparsity, bool_rt
+ );
+ bool_st.resize(n * q );
+ ok = rev_sparse_jac(q, bool_rt, bool_st, x);
+ if( ! ok )
+ ok = rev_sparse_jac(q, bool_rt, bool_st);
+ if( ! ok )
+ { msg = afun_name() + msg + " sparsity = bool_sparsity_enum";
+ CPPAD_ASSERT_KNOWN(false, msg.c_str());
+ }
+ local::set_internal_sparsity(zero_empty, input_empty,
+ transpose, x_index, var_sparsity, bool_st
+ );
+ }
+ else
+ { CPPAD_ASSERT_UNKNOWN( sparsity_ == set_sparsity_enum );
+ vector< std::set >& set_rt ( work_[thread]->set_r );
+ vector< std::set >& set_st ( work_[thread]->set_s );
+ local::get_internal_sparsity(
+ transpose, y_index, var_sparsity, set_rt
+ );
+ set_st.resize(n);
+ ok = rev_sparse_jac(q, set_rt, set_st, x);
+ if( ! ok )
+ ok = rev_sparse_jac(q, set_rt, set_st);
+ if( ! ok )
+ { msg = afun_name() + msg + " sparsity = set_sparsity_enum";
+ CPPAD_ASSERT_KNOWN(false, msg.c_str());
+ }
+ local::set_internal_sparsity(zero_empty, input_empty,
+ transpose, x_index, var_sparsity, set_st
+ );
+ }
+ return;
+}
+/*
+-------------------------------------- ---------------------------------------
+$begin atomic_for_sparse_hes$$
+$spell
+ sq
+ mul.hpp
+ vx
+ afun
+ Jacobian
+ jac
+ CppAD
+ std
+ bool
+ hes
+ const
+$$
+
+$section Atomic Forward Hessian Sparsity Patterns$$
+
+$head Syntax$$
+$icode%ok% = %afun%.for_sparse_hes(%vx%, %r%, %s%, %h%, %x%)%$$
+
+$head Deprecated 2016-06-27$$
+$icode%ok% = %afun%.for_sparse_hes(%vx%, %r%, %s%, %h%)%$$
+
+$head Purpose$$
+This function is used by $cref ForSparseHes$$ to compute
+Hessian sparsity patterns.
+If you are using $cref ForSparseHes$$,
+one of the versions of this
+virtual function must be defined by the
+$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
+$pre
+
+$$
+Given a $cref/sparsity pattern/glossary/Sparsity Pattern/$$ for
+a diagonal matrix $latex R \in B^{n \times n}$$, and
+a row vector $latex S \in B^{1 \times m}$$,
+this routine computes the sparsity pattern for
+$latex \[
+ H(x) = R^\R{T} \cdot (S \cdot f)^{(2)}( x ) \cdot R
+\] $$
+
+$head Implementation$$
+If you are using and $cref ForSparseHes$$,
+this virtual function must be defined by the
+$cref/atomic_user/atomic_ctor/atomic_user/$$ class.
+
+$subhead vx$$
+The argument $icode vx$$ has prototype
+$codei%
+ const CppAD:vector