From 17332452caa735264ea5cf5d21eda46f31210804 Mon Sep 17 00:00:00 2001 From: Brad Bell Date: Sat, 2 Mar 2024 09:02:29 -0700 Subject: [PATCH] master: valvector: combine split and join in split_join.hpp and finish out AD version of reverse. --- CMakeLists.txt | 2 +- appendix/whats_new/2024.xrst | 5 + example/valvector/ad_join.cpp | 15 +- example/valvector/ad_split.cpp | 16 +- include/cppad/example/valvector/ad_join.hpp | 315 --------- include/cppad/example/valvector/ad_split.hpp | 313 --------- include/cppad/example/valvector/class.hpp | 3 +- .../cppad/example/valvector/split_join.hpp | 602 ++++++++++++++++++ user_guide.xrst | 2 +- 9 files changed, 639 insertions(+), 634 deletions(-) delete mode 100644 include/cppad/example/valvector/ad_join.hpp delete mode 100644 include/cppad/example/valvector/ad_split.hpp create mode 100644 include/cppad/example/valvector/split_join.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 797b0ff27..516d03b93 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,7 +16,7 @@ IF( POLICY CMP0054 ) ENDIF( POLICY CMP0054 ) # # cppad_version is used by version.sh to get the version number. -SET(cppad_version "20240301") +SET(cppad_version "20240302") SET(cppad_url "https://coin-or.github.io/CppAD" ) SET(cppad_description "Differentiation of C++ Algorithms" ) IF( NOT DEFINED CMAKE_BUILD_TYPE) diff --git a/appendix/whats_new/2024.xrst b/appendix/whats_new/2024.xrst index 520ca8fb4..b985a5d44 100644 --- a/appendix/whats_new/2024.xrst +++ b/appendix/whats_new/2024.xrst @@ -28,6 +28,11 @@ Release Notes for 2024 mm-dd ***** +03-02 +===== +Add the :ref:`valvector_ad_split-name` and :ref:`valvector_ad_join-name` +examples. + 03-01 ===== Fix bug in the optimizer :ref:`optimize@options@val_graph` option. diff --git a/example/valvector/ad_join.cpp b/example/valvector/ad_join.cpp index 7a10e8ffa..ac5cea4ba 100644 --- a/example/valvector/ad_join.cpp +++ b/example/valvector/ad_join.cpp @@ -2,7 +2,7 @@ // SPDX-FileCopyrightText: Bradley M. Bell // SPDX-FileContributor: 2024 Bradley M. Bell // --------------------------------------------------------------------------- -# include +# include # include /* {xrst_begin valvector_ad_join.cpp} @@ -153,6 +153,19 @@ bool ad_join(void) for(size_t j = 0; j < n; ++j) ok &= dz[0][j] == 2.0 * x[j][0]; // + // h + CPPAD_TESTVECTOR( ad_valvector ) aw(m), adw(n); + aw[0] = valvector( 1.0 ); + CppAD::Independent(ax); + af.Forward(0, ax); + adw = af.Reverse(1, aw); + CppAD::ADFun h(ax, adw); + // + // ok + dw = h.Forward(0, x); + for(size_t j = 0; j < n; ++j) + ok &= dw[j][0] == 2.0 * x[j][0]; + // // ok f.optimize(); z = f.Forward(0, x); diff --git a/example/valvector/ad_split.cpp b/example/valvector/ad_split.cpp index 8a925f970..36a64fb9c 100644 --- a/example/valvector/ad_split.cpp +++ b/example/valvector/ad_split.cpp @@ -2,7 +2,7 @@ // SPDX-FileCopyrightText: Bradley M. Bell // SPDX-FileContributor: 2024 Bradley M. Bell // --------------------------------------------------------------------------- -# include +# include # include /* {xrst_begin valvector_ad_split.cpp} @@ -135,6 +135,20 @@ bool ad_split(void) for(size_t i = 0; i < m; ++i) ok &= dy[i][0] == 2.0 * x[0][i]; // + // h + CPPAD_TESTVECTOR( ad_valvector ) aw(m), adw(n); + for(size_t i = 0; i < m; ++i) + aw[i] = valvector( 1.0 ); + CppAD::Independent(ax); + af.Forward(0, ax); + adw = af.Reverse(1, aw); + CppAD::ADFun h(ax, adw); + // + // ok + dw = h.Forward(0, x); + for(size_t i = 0; i < m; ++i) + ok &= dw[0][i] == 2.0 * x[0][i]; + // // ok f.optimize(); y = f.Forward(0, x); diff --git a/include/cppad/example/valvector/ad_join.hpp b/include/cppad/example/valvector/ad_join.hpp deleted file mode 100644 index 5982d5159..000000000 --- a/include/cppad/example/valvector/ad_join.hpp +++ /dev/null @@ -1,315 +0,0 @@ -# ifndef CPPAD_EXAMPLE_VALVECTOR_AD_JOIN_HPP -# define CPPAD_EXAMPLE_VALVECTOR_AD_JOIN_HPP -// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later -// SPDX-FileCopyrightText: Bradley M. Bell -// SPDX-FileContributor: 2024 Bradley M. Bell -// ---------------------------------------------------------------------------- -# include -# include -/* -{xrst_begin valvector_ad_join} -{xrst_spell - ajoin - valvectors -} - -Join a Vector of valvectors Into One valvector -############################################## -Join a vector of valvectors each with size one, into one valvector. - -Under Construction -****************** - -Syntax -****** -| valvector_ad_join *ajoin* -| *ajoin( *ax_vec* , *ay* ) - - -m -* -We use *m* to denote *ax_vec* .size() . - -ax_vec -****** -This is a :ref:`SimpleVector-name` with elements of type CppAD::AD . -It is ``const`` and is passed by reference. -For *i* = 0 , ... , *m* - 1 , -the size *ax* [ *i* ].size() is one. - -ay -** -This CppAD::AD is passed by reference and its input value -does not matter. -Upon return, its size in *m* and -for *i* = 0 , ... , *m* - 1 :: - - ay[i] = ax_vec[i][0] - - -{xrst_toc_hidden - example/valvector/ad_join.cpp -} -Example -******* -The file :ref:`valvector_ad_join.cpp-name` is an example and test -of this operation. - -{xrst_end valvector_ad_join} -*/ -class valvector_ad_join_atom : public CppAD::atomic_four { -public: - // - // ctor - valvector_ad_join_atom(const std::string& name) : - CppAD::atomic_four(name) - { } -private: - // - // ad_valvector - typedef CppAD::AD ad_valvector; - // - // scalar_type - typedef valvector::scalar_type scalar_type; - // ------------------------------------------------------------------------ - // for_type - bool for_type( - size_t call_id , - const CppAD::vector& type_x , - CppAD::vector& type_y ) override - { - assert( call_id == 0 ); // default value - assert( type_y.size() == 1 ); // m - // - // type_y - size_t n = type_x.size(); - type_y[0] = CppAD::identical_zero_enum; - for(size_t j = 0; j < n; ++j) - type_y[0] = std::max( type_y[0], type_x[j] ); - return true; - } - // ------------------------------------------------------------------------ - // forward - bool forward( - size_t call_id , - const CppAD::vector& select_y , - size_t order_low , - size_t order_up , - const CppAD::vector& taylor_x , - CppAD::vector& taylor_y ) override - { // - // ok - bool ok = true; - // - // q, n - size_t q = order_up + 1; - size_t n = taylor_x.size() / q; - // -# ifndef NDEBUG - size_t m = taylor_y.size() / q; - assert( call_id == 0 ); - assert( m == 1 ); - assert( m == select_y.size() ); - for(size_t k = 0; k < q; ++k) - assert( taylor_x[k].size() == 1 ); -# endif - // - // taylor_y - if( select_y[0] ) - { for(size_t k = order_low; k < q; ++k) - { // - // taylor_y[k] - taylor_y[k].resize(n); - for(size_t j = 0; j < n; ++j) - taylor_y[k][j] = taylor_x[j * q + k][0]; - } - } - // - return ok; - } - // ------------------------------------------------------------------------ - // forward - bool forward( - size_t call_id , - const CppAD::vector& select_y , - size_t order_low , - size_t order_up , - const CppAD::vector& ataylor_x , - CppAD::vector& ataylor_y ) override - { // - // ok - bool ok = true; - // - // q, n - size_t q = order_up + 1; - size_t n = ataylor_x.size() / q; - // -# ifndef NDEBUG - size_t m = ataylor_y.size() / q; - assert( call_id == 0 ); - assert( m == 1 ); - assert( m == select_y.size() ); -# endif - if( select_y[0] ) - { // - // ay_k - CppAD::vector ay_k(1); - CppAD::vector ax_k(n); - // - // ataylor_y - for(size_t k = order_low; k < q; ++k) - { for(size_t j = 0; j < n; ++j) - ax_k[j] = ataylor_x[j * q + k]; - (*this)(call_id, ax_k, ay_k); - ataylor_y[k] = ay_k[0]; - } - } - // - return ok; - } - // ------------------------------------------------------------------------ - // reverse - bool reverse( - size_t call_id , - const CppAD::vector& select_x , - size_t order_up , - const CppAD::vector& taylor_x , - const CppAD::vector& taylor_y , - CppAD::vector& partial_x , - const CppAD::vector& partial_y ) override - { // - // ok - bool ok = true; - // - // q, m - size_t q = order_up + 1; - size_t n = taylor_x.size() / q; - // -# ifndef NDEBUG - size_t m = taylor_y.size() / q; - assert( call_id == 0 ); - assert( m == 1 ); - assert( n == select_x.size() ); -# endif - // - // partial_x - for(size_t k = 0; k < q; ++k) - { assert( taylor_y[k].size() == n ); - for(size_t j = 0; j < n; ++j) - { partial_x[j * q + k].resize(1); - partial_x[j * q + k][0] = partial_y[k][j]; - } - } - return ok; - } - // ------------------------------------------------------------------------ - // jac_sparsity - bool jac_sparsity( - size_t call_id , - bool dependency , - const CppAD::vector& ident_zero_x , - const CppAD::vector& select_x , - const CppAD::vector& select_y , - CppAD::sparse_rc< CppAD::vector >& pattern_out ) override - { // - // ok - bool ok = true; - // - // m, n - size_t m = select_y.size(); - size_t n = select_x.size(); - // - assert( call_id == 0 ); - assert( m == 1 ); - // - // nnz - size_t nnz = 0; - if( select_y[0] ) - { for(size_t j = 0; j < n; ++j) - { if( select_x[j] ) - ++nnz; - } - } - // - // pattern_out - pattern_out.resize(m, n, nnz); - size_t k = 0; - if( select_y[0] ) - { for(size_t j = 0; j < n; ++j) - { if( select_x[j] ) - pattern_out.set(k++, 0, j); - } - } - return ok; - } - // ------------------------------------------------------------------------ - // hes_sparsity - bool hes_sparsity( - size_t call_id , - const CppAD::vector& ident_zero_x , - const CppAD::vector& select_x , - const CppAD::vector& select_y , - CppAD::sparse_rc< CppAD::vector >& pattern_out ) override - { // - // ok - bool ok = true; - // - // n - size_t n = select_x.size(); - // - assert( call_id == 0 ); - assert( n == select_x.size() ); - // - // pattern_out - size_t nnz = 0; - pattern_out.resize(n, n, nnz); - // - return ok; - } - // ------------------------------------------------------------------------ - // rev_depend - bool rev_depend( - size_t call_id , - const CppAD::vector& ident_zero_x , - CppAD::vector& depend_x , - const CppAD::vector& depend_y ) override - { // - // ok - bool ok = true; - // - // n - size_t n = depend_x.size(); - // -# ifndef NDEBUG - size_t m = depend_y.size(); - assert( call_id == 0 ); - assert( m == 1 ); -# endif - // - // depend_x - for(size_t j = 0; j < n; ++j) - depend_x[j] = depend_y[0]; - // - return ok; - } -}; - -class valvector_ad_join { -private: - typedef CppAD::AD ad_valvector; - valvector_ad_join_atom atomic_fun_; -public: - valvector_ad_join(void) : atomic_fun_("valvector_ad_join_atom") - { } - template - void operator()(const ADVector& ax, ad_valvector& ay) - { - size_t call_id = 0; - ADVector ay_vec(1); - atomic_fun_(call_id, ax, ay_vec); - ay = ay_vec[0]; - return; - } -}; - -# endif diff --git a/include/cppad/example/valvector/ad_split.hpp b/include/cppad/example/valvector/ad_split.hpp deleted file mode 100644 index 8f31e6471..000000000 --- a/include/cppad/example/valvector/ad_split.hpp +++ /dev/null @@ -1,313 +0,0 @@ -# ifndef CPPAD_EXAMPLE_VALVECTOR_AD_SPLIT_HPP -# define CPPAD_EXAMPLE_VALVECTOR_AD_SPLIT_HPP -// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later -// SPDX-FileCopyrightText: Bradley M. Bell -// SPDX-FileContributor: 2024 Bradley M. Bell -// ---------------------------------------------------------------------------- -# include -# include -/* -{xrst_begin valvector_ad_split} -{xrst_spell - asplit - valvectors -} - -Split A valvector Into a Vector of valvector -############################################ -Split a valvector into a vector of valvectors each with size one. - -Syntax -****** -| valvector_ad_split *asplit* -| *asplit( *ax* , *ay_vec* ) - -m -* -We use *m* to denote *ay_vec* .size() - -ax -** -This CppAD::AD is ``const`` and passed by reference. -The size *ax*.size() must be equal to one or *m* - -ay_vec -****** -This is a :ref:`SimpleVector-name` with elements of type CppAD::AD -and is passed by reference. -The size of *ay_vec* is not changed. -Upon return, for *i* = 0 , ... , *m* - 1 , -the size *ay* [ *i* ].size() is one and:: - - ay[i][0] = ax[i] - -{xrst_toc_hidden - example/valvector/ad_split.cpp -} -Example -******* -The file :ref:`valvector_ad_split.cpp-name` is an example and test -of this operation. - -{xrst_end valvector_ad_split} -*/ -class valvector_ad_split_atom : public CppAD::atomic_four { -public: - // - // ctor - valvector_ad_split_atom(const std::string& name) : - CppAD::atomic_four(name) - { } -private: - // - // ad_valvector - typedef CppAD::AD ad_valvector; - // - // scalar_type - typedef valvector::scalar_type scalar_type; - // ------------------------------------------------------------------------ - // for_type - bool for_type( - size_t call_id , - const CppAD::vector& type_x , - CppAD::vector& type_y ) override - { - assert( call_id == 0 ); // default value - assert( type_x.size() == 1 ); // n - // - // type_y - size_t m = type_y.size(); - for(size_t i = 0; i < m; ++i) - type_y[i] = type_x[0]; - return true; - } - // ------------------------------------------------------------------------ - // forward - bool forward( - size_t call_id , - const CppAD::vector& select_y , - size_t order_low , - size_t order_up , - const CppAD::vector& taylor_x , - CppAD::vector& taylor_y ) override - { // - // ok - bool ok = true; - // - // q, m - size_t q = order_up + 1; - size_t m = taylor_y.size() / q; - // -# ifndef NDEBUG - size_t n = taylor_x.size() / q; - assert( call_id == 0 ); - assert( n == 1 ); - assert( m == select_y.size() ); - for(size_t k = 0; k < q; ++k) - assert( taylor_x[k].size() == 1 || taylor_x[k].size() == m ); -# endif - // - // taylor_y - for(size_t i = 0; i < m; ++i) if( select_y[i] ) - { for(size_t k = order_low; k < q; ++k) - { // - // taylor_y[i * q + k] - taylor_y[i * q + k].resize(1); - taylor_y[i * q + k][0] = taylor_x[k][i]; - } - } - // - return ok; - } - // ------------------------------------------------------------------------ - // forward - bool forward( - size_t call_id , - const CppAD::vector& select_y , - size_t order_low , - size_t order_up , - const CppAD::vector& ataylor_x , - CppAD::vector& ataylor_y ) override - { // - // ok - bool ok = true; - // - // q, m - size_t q = order_up + 1; - size_t m = ataylor_y.size() / q; - // -# ifndef NDEBUG - size_t n = ataylor_x.size() / q; - assert( call_id == 0 ); - assert( n == 1 ); - assert( m == select_y.size() ); -# endif - // - // ay_k - CppAD::vector ay_k(m); - CppAD::vector ax_k(1); - // - // ataylor_y - for(size_t k = order_low; k < q; ++k) - { ax_k[0] = ataylor_x[k]; - (*this)(call_id, ax_k, ay_k); - { for(size_t i = 0; i < m; ++i) - ataylor_y[i * q + k] = ay_k[i]; - } - } - // - return ok; - } - // ------------------------------------------------------------------------ - // reverse - bool reverse( - size_t call_id , - const CppAD::vector& select_x , - size_t order_up , - const CppAD::vector& taylor_x , - const CppAD::vector& taylor_y , - CppAD::vector& partial_x , - const CppAD::vector& partial_y ) override - { // - // ok - bool ok = true; - // - // q, m - size_t q = order_up + 1; - size_t m = taylor_y.size() / q; - // -# ifndef NDEBUG - size_t n = taylor_x.size() / q; - assert( call_id == 0 ); - assert( n == 1 ); - assert( 1 == select_x.size() ); -# endif - if( ! select_x[0] ) - return ok; - // - // partial_x - for(size_t k = 0; k < q; ++k) - { partial_x[k].resize(m); - for(size_t i = 0; i < m; ++i) - { // - // partial_x[k][i] - assert( taylor_y[i * q + k].size() == 1 ); - partial_x[k][i] = partial_y[i * q + k][0]; - } - } - // - return ok; - } - // ------------------------------------------------------------------------ - // jac_sparsity - bool jac_sparsity( - size_t call_id , - bool dependency , - const CppAD::vector& ident_zero_x , - const CppAD::vector& select_x , - const CppAD::vector& select_y , - CppAD::sparse_rc< CppAD::vector >& pattern_out ) override - { // - // ok - bool ok = true; - // - // m, n - size_t m = select_y.size(); - size_t n = select_x.size(); - // - assert( call_id == 0 ); - assert( n == 1 ); - // - // nnz - size_t nnz = 0; - if( select_x[0] ) - { for(size_t i = 0; i < m; ++i) - { if( select_y[i] ) - ++nnz; - } - } - // - // pattern_out - pattern_out.resize(m, n, nnz); - size_t k = 0; - if( select_x[0] ) - { for(size_t i = 0; i < m; ++i) - { if( select_y[i] ) - pattern_out.set(k++, i, 0); - } - } - return ok; - } - // ------------------------------------------------------------------------ - // hes_sparsity - bool hes_sparsity( - size_t call_id , - const CppAD::vector& ident_zero_x , - const CppAD::vector& select_x , - const CppAD::vector& select_y , - CppAD::sparse_rc< CppAD::vector >& pattern_out ) override - { // - // ok - bool ok = true; - // - // n - size_t n = select_x.size(); - // - assert( call_id == 0 ); - assert( n == 1 ); - assert( n == select_x.size() ); - // - // pattern_out - size_t nnz = 0; - pattern_out.resize(n, n, nnz); - // - return ok; - } - // ------------------------------------------------------------------------ - // rev_depend - bool rev_depend( - size_t call_id , - const CppAD::vector& ident_zero_x , - CppAD::vector& depend_x , - const CppAD::vector& depend_y ) override - { // - // ok - bool ok = true; - // - // m - size_t m = depend_y.size(); - // -# ifndef NDEBUG - size_t n = depend_x.size(); - assert( call_id == 0 ); - assert( n == 1 ); -# endif - // - // depend_x - depend_x[0] = false; - for(size_t i = 0; i < m; ++i) - depend_x[0] |= depend_y[i]; - // - return ok; - } -}; - -class valvector_ad_split { -private: - typedef CppAD::AD ad_valvector; - valvector_ad_split_atom atomic_fun_; -public: - valvector_ad_split(void) : atomic_fun_("valvector_ad_split_atom") - { } - template - void operator()(const ad_valvector& ax, ADVector& ay) - { - size_t call_id = 0; - ADVector ax_vec(1); - ax_vec[0] = ax; - atomic_fun_(call_id, ax_vec, ay); - return; - } -}; - -# endif diff --git a/include/cppad/example/valvector/class.hpp b/include/cppad/example/valvector/class.hpp index a3d7771f5..ea37cd668 100644 --- a/include/cppad/example/valvector/class.hpp +++ b/include/cppad/example/valvector/class.hpp @@ -37,8 +37,7 @@ git the expected results. Operations ********** {xrst_toc_table after - include/cppad/example/valvector/ad_split.hpp - include/cppad/example/valvector/ad_join.hpp + include/cppad/example/valvector/split_join.hpp example/valvector/get_started.cpp } diff --git a/include/cppad/example/valvector/split_join.hpp b/include/cppad/example/valvector/split_join.hpp new file mode 100644 index 000000000..73141ad8f --- /dev/null +++ b/include/cppad/example/valvector/split_join.hpp @@ -0,0 +1,602 @@ +# ifndef CPPAD_EXAMPLE_VALVECTOR_SPLIT_JOIN_HPP +# define CPPAD_EXAMPLE_VALVECTOR_SPLIT_JOIN_HPP +// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later +// SPDX-FileCopyrightText: Bradley M. Bell +// SPDX-FileContributor: 2024 Bradley M. Bell +/* +------------------------------------------------------------------------------- +{xrst_begin valvector_ad_split} +{xrst_spell + asplit + valvectors +} + +Split A valvector Into a Vector of valvector +############################################ +Split a valvector into a vector of valvectors each with size one. + +Syntax +****** +| valvector_ad_split *asplit* +| *asplit( *ax* , *ay_vec* ) + +m +* +We use *m* to denote *ay_vec* .size() + +ax +** +This CppAD::AD is ``const`` and passed by reference. +The size *ax*.size() must be equal to one or *m* + +ay_vec +****** +This is a :ref:`SimpleVector-name` with elements of type CppAD::AD +and is passed by reference. +The size of *ay_vec* is not changed. +Upon return, for *i* = 0 , ... , *m* - 1 , +the size *ay* [ *i* ].size() is one and:: + + ay[i][0] = ax[i] + +{xrst_toc_hidden + example/valvector/ad_split.cpp +} +Example +******* +The file :ref:`valvector_ad_split.cpp-name` is an example and test +of this operation. + +{xrst_end valvector_ad_split} +------------------------------------------------------------------------------- +{xrst_begin valvector_ad_join} +{xrst_spell + ajoin + valvectors +} + +Join a Vector of valvectors Into One valvector +############################################## +Join a vector of valvectors each with size one, into one valvector. + +Under Construction +****************** + +Syntax +****** +| valvector_ad_join *ajoin* +| *ajoin( *ax_vec* , *ay* ) + + +m +* +We use *m* to denote *ax_vec* .size() . + +ax_vec +****** +This is a :ref:`SimpleVector-name` with elements of type CppAD::AD . +It is ``const`` and is passed by reference. +For *i* = 0 , ... , *m* - 1 , +the size *ax* [ *i* ].size() is one. + +ay +** +This CppAD::AD is passed by reference and its input value +does not matter. +Upon return, its size in *m* and +for *i* = 0 , ... , *m* - 1 :: + + ay[i] = ax_vec[i][0] + + +{xrst_toc_hidden + example/valvector/ad_join.cpp +} +Example +******* +The file :ref:`valvector_ad_join.cpp-name` is an example and test +of this operation. + +{xrst_end valvector_ad_join} +------------------------------------------------------------------------------- +*/ +# include +# include +// +// valvector_split_join +// split and join are in same atomic function so that in AD version of reverse +// split can call join and join can call split. +class valvector_split_join : public CppAD::atomic_four { +public: + // + // ctor + valvector_split_join(const std::string& name) : + CppAD::atomic_four(name) + { } + // + // is_split + static bool is_split(size_t call_id) + { return call_id == 1; } + // + // is_join + static bool is_join(size_t call_id) + { return call_id == 2; } +private: + // + // ad_valvector + typedef CppAD::AD ad_valvector; + // + // ------------------------------------------------------------------------ + // for_type + bool for_type( + size_t call_id , + const CppAD::vector& type_x , + CppAD::vector& type_y ) override + { + if( is_split(call_id) ) + { // + assert( type_x.size() == 1 ); // n + // + // type_y + size_t m = type_y.size(); + for(size_t i = 0; i < m; ++i) + type_y[i] = type_x[0]; + } + else + { assert( is_join(call_id) ); + // + assert( type_y.size() == 1 ); // m + // + // type_y + size_t n = type_x.size(); + type_y[0] = CppAD::identical_zero_enum; + for(size_t j = 0; j < n; ++j) + type_y[0] = std::max( type_y[0], type_x[j] ); + } + return true; + } + // ------------------------------------------------------------------------ + // forward + bool forward( + size_t call_id , + const CppAD::vector& select_y , + size_t order_low , + size_t order_up , + const CppAD::vector& taylor_x , + CppAD::vector& taylor_y ) override + { // + if( is_split(call_id) ) + { // + // q, m + size_t q = order_up + 1; + size_t m = taylor_y.size() / q; + // +# ifndef NDEBUG + size_t n = taylor_x.size() / q; + assert( n == 1 ); + assert( m == select_y.size() ); + for(size_t k = 0; k < q; ++k) + assert( taylor_x[k].size() == 1 || taylor_x[k].size() == m ); +# endif + // + // taylor_y + for(size_t i = 0; i < m; ++i) if( select_y[i] ) + { for(size_t k = order_low; k < q; ++k) + { // + // taylor_y[i * q + k] + taylor_y[i * q + k].resize(1); + taylor_y[i * q + k][0] = taylor_x[k][i]; + } + } + } + else + { assert( is_join(call_id) ); + // + // q, n + size_t q = order_up + 1; + size_t n = taylor_x.size() / q; + // +# ifndef NDEBUG + size_t m = taylor_y.size() / q; + assert( m == 1 ); + assert( m == select_y.size() ); + for(size_t k = 0; k < q; ++k) + assert( taylor_x[k].size() == 1 ); +# endif + // + // taylor_y + if( select_y[0] ) + { for(size_t k = order_low; k < q; ++k) + { // + // taylor_y[k] + taylor_y[k].resize(n); + for(size_t j = 0; j < n; ++j) + taylor_y[k][j] = taylor_x[j * q + k][0]; + } + } + // + } + return true; + } + // ------------------------------------------------------------------------ + // forward + bool forward( + size_t call_id , + const CppAD::vector& select_y , + size_t order_low , + size_t order_up , + const CppAD::vector& ataylor_x , + CppAD::vector& ataylor_y ) override + { // + if( is_split(call_id) ) + { // + // q, m + size_t q = order_up + 1; + size_t m = ataylor_y.size() / q; + // +# ifndef NDEBUG + size_t n = ataylor_x.size() / q; + assert( n == 1 ); + assert( m == select_y.size() ); +# endif + // + // ax_k, ay_k + CppAD::vector ax_k(1); + CppAD::vector ay_k(m); + // + // ataylor_y + for(size_t k = order_low; k < q; ++k) + { ax_k[0] = ataylor_x[k]; + (*this)(call_id, ax_k, ay_k); + { for(size_t i = 0; i < m; ++i) + ataylor_y[i * q + k] = ay_k[i]; + } + } + } + else + { assert( is_join(call_id) ); + // + // q, n + size_t q = order_up + 1; + size_t n = ataylor_x.size() / q; + // +# ifndef NDEBUG + size_t m = ataylor_y.size() / q; + assert( m == 1 ); + assert( m == select_y.size() ); +# endif + if( select_y[0] ) + { // + // ax_k, ay_k + CppAD::vector ax_k(n); + CppAD::vector ay_k(1); + // + // ataylor_y + for(size_t k = order_low; k < q; ++k) + { for(size_t j = 0; j < n; ++j) + ax_k[j] = ataylor_x[j * q + k]; + (*this)(call_id, ax_k, ay_k); + ataylor_y[k] = ay_k[0]; + } + } + } + return true; + } + // ------------------------------------------------------------------------ + // reverse + bool reverse( + size_t call_id , + const CppAD::vector& select_x , + size_t order_up , + const CppAD::vector& taylor_x , + const CppAD::vector& taylor_y , + CppAD::vector& partial_x , + const CppAD::vector& partial_y ) override + { // + if( is_split(call_id) ) + { // + // q, m + size_t q = order_up + 1; + size_t m = taylor_y.size() / q; + // +# ifndef NDEBUG + size_t n = taylor_x.size() / q; + assert( n == 1 ); + assert( 1 == select_x.size() ); +# endif + if( ! select_x[0] ) + return true; + // + // partial_x + for(size_t k = 0; k < q; ++k) + { partial_x[k].resize(m); + for(size_t i = 0; i < m; ++i) + { // + // partial_x[k][i] + assert( taylor_y[i * q + k].size() == 1 ); + partial_x[k][i] = partial_y[i * q + k][0]; + } + } + } + else + { assert( is_join(call_id) ); + // + // q, m + size_t q = order_up + 1; + size_t n = taylor_x.size() / q; + // +# ifndef NDEBUG + size_t m = taylor_y.size() / q; + assert( m == 1 ); + assert( n == select_x.size() ); +# endif + // + // partial_x + for(size_t k = 0; k < q; ++k) + { assert( taylor_y[k].size() == n ); + for(size_t j = 0; j < n; ++j) + { partial_x[j * q + k].resize(1); + partial_x[j * q + k][0] = partial_y[k][j]; + } + } + } + return true; + } + // ------------------------------------------------------------------------ + // reverse + bool reverse( + size_t call_id , + const CppAD::vector& select_x , + size_t order_up , + const CppAD::vector& ataylor_x , + const CppAD::vector& ataylor_y , + CppAD::vector& apartial_x , + const CppAD::vector& apartial_y ) override + { // + if( is_split(call_id) ) + { // + // q, m + size_t q = order_up + 1; + size_t m = ataylor_y.size() / q; + // +# ifndef NDEBUG + size_t n = ataylor_x.size() / q; + assert( n == 1 ); + assert( 1 == select_x.size() ); +# endif + if( ! select_x[0] ) + return true; + // + // ay_k, ax_k + CppAD::vector ay_k(m); + CppAD::vector ax_k(1); + // + // join_call_id + size_t join_call_id = 2; + assert( is_join(join_call_id) ); + // + // apartial_x + for(size_t k = 0; k < q; ++k) + { for(size_t i = 0; i < m; ++i) + ay_k[i] = apartial_y[i * q + k]; + // split calls join here + (*this)(join_call_id, ay_k, ax_k); + apartial_x[k] = ax_k[0]; + } + } + else + { assert( is_join(call_id) ); + // + // q, m + size_t q = order_up + 1; + size_t n = ataylor_x.size() / q; + // +# ifndef NDEBUG + size_t m = ataylor_y.size() / q; + assert( m == 1 ); + assert( n == select_x.size() ); +# endif + // + // ay_k, ax_k + CppAD::vector ay_k(1); + CppAD::vector ax_k(n); + // + // split_call_id + size_t split_call_id = 1; + assert( is_split(split_call_id) ); + // + // apartial_x + for(size_t k = 0; k < q; ++k) + { ay_k[0] = apartial_y[k]; + // Join calls split here + (*this)(split_call_id, ay_k, ax_k); + for(size_t j = 0; j < n; ++j) + apartial_x[j * q + k] = ax_k[j]; + } + } + return true; + } + // ------------------------------------------------------------------------ + // jac_sparsity + bool jac_sparsity( + size_t call_id , + bool dependency , + const CppAD::vector& ident_zero_x , + const CppAD::vector& select_x , + const CppAD::vector& select_y , + CppAD::sparse_rc< CppAD::vector >& pattern_out ) override + { // + if( is_split(call_id) ) + { // + // m, n + size_t m = select_y.size(); + size_t n = select_x.size(); + // + assert( n == 1 ); + // + // nnz + size_t nnz = 0; + if( select_x[0] ) + { for(size_t i = 0; i < m; ++i) + { if( select_y[i] ) + ++nnz; + } + } + // + // pattern_out + pattern_out.resize(m, n, nnz); + size_t k = 0; + if( select_x[0] ) + { for(size_t i = 0; i < m; ++i) + { if( select_y[i] ) + pattern_out.set(k++, i, 0); + } + } + } + else + { assert( is_join(call_id) ); + // + // m, n + size_t m = select_y.size(); + size_t n = select_x.size(); + // + assert( m == 1 ); + // + // nnz + size_t nnz = 0; + if( select_y[0] ) + { for(size_t j = 0; j < n; ++j) + { if( select_x[j] ) + ++nnz; + } + } + // + // pattern_out + pattern_out.resize(m, n, nnz); + size_t k = 0; + if( select_y[0] ) + { for(size_t j = 0; j < n; ++j) + { if( select_x[j] ) + pattern_out.set(k++, 0, j); + } + } + } + return true; + } + // ------------------------------------------------------------------------ + // hes_sparsity + bool hes_sparsity( + size_t call_id , + const CppAD::vector& ident_zero_x , + const CppAD::vector& select_x , + const CppAD::vector& select_y , + CppAD::sparse_rc< CppAD::vector >& pattern_out ) override + { // + if( is_split(call_id) ) + { // + // n + size_t n = select_x.size(); + // + assert( n == 1 ); + assert( n == select_x.size() ); + // + // pattern_out + size_t nnz = 0; + pattern_out.resize(n, n, nnz); + } + else + { assert( is_join(call_id) ); + // + // n + size_t n = select_x.size(); + // + assert( n == select_x.size() ); + // + // pattern_out + size_t nnz = 0; + pattern_out.resize(n, n, nnz); + } + return true; + } + // ------------------------------------------------------------------------ + // rev_depend + bool rev_depend( + size_t call_id , + const CppAD::vector& ident_zero_x , + CppAD::vector& depend_x , + const CppAD::vector& depend_y ) override + { // + if( is_split(call_id) ) + { // + // m + size_t m = depend_y.size(); + // +# ifndef NDEBUG + size_t n = depend_x.size(); + assert( n == 1 ); +# endif + // + // depend_x + depend_x[0] = false; + for(size_t i = 0; i < m; ++i) + depend_x[0] |= depend_y[i]; + } + else + { assert( is_join(call_id) ); + // + // n + size_t n = depend_x.size(); + // +# ifndef NDEBUG + size_t m = depend_y.size(); + assert( m == 1 ); +# endif + // + // depend_x + for(size_t j = 0; j < n; ++j) + depend_x[j] = depend_y[0]; + } + return true; + } +}; +// --------------------------------------------------------------------------- +// valvector_ad_split +class valvector_ad_split { +private: + typedef CppAD::AD ad_valvector; + valvector_split_join atomic_fun_; +public: + valvector_ad_split(void) : atomic_fun_("valvector_split_join") + { } + template + void operator()(const ad_valvector& ax, ADVector& ay) + { + size_t call_id = 1; + assert( atomic_fun_.is_split(call_id) ); + ADVector ax_vec(1); + ax_vec[0] = ax; + atomic_fun_(call_id, ax_vec, ay); + return; + } +}; +// --------------------------------------------------------------------------- +// valvector_ad_join +class valvector_ad_join { +private: + typedef CppAD::AD ad_valvector; + valvector_split_join atomic_fun_; +public: + valvector_ad_join(void) : atomic_fun_("valvector_split_join") + { } + template + void operator()(const ADVector& ax, ad_valvector& ay) + { + size_t call_id = 2; + ADVector ay_vec(1); + assert( atomic_fun_.is_join(call_id) ); + atomic_fun_(call_id, ax, ay_vec); + ay = ay_vec[0]; + return; + } +}; +# endif diff --git a/user_guide.xrst b/user_guide.xrst index af910f98a..4b6844ade 100644 --- a/user_guide.xrst +++ b/user_guide.xrst @@ -13,7 +13,7 @@ {xrst_comment BEGIN: Before changing see new_release.sh and check_version.sh} -cppad-20240301: CppAD User's Manual +cppad-20240302: CppAD User's Manual ################################### .. image:: {xrst_dir coin.png}