-
|
I want to compute Jacobian of a function The source code: #include <iostream> // standard input/output
#include <cppad/cppad.hpp> // the CppAD package
// main program
int main(void)
{
using CppAD::AD;
using CppAD::ADFun;
using CppAD::Independent;
size_t n = 1;
size_t m = 1;
size_t k = 1;
CPPAD_TESTVECTOR(AD<double>) ax(n), ay(m), ap(k);
ap[0] = AD<double>(4.0);
Independent(ax, ap);
ay[0] = ap[0] * ax[0] * ax[0];
ADFun<double> f(ax, ay);
ADFun<AD<double>, double> af = f.base2ad();
Independent(ax, ap);
CPPAD_TESTVECTOR(AD<double>) az(m * n);
az = af.Jacobian(ax);
ADFun<double> g;
g.Dependent(ax, az);
g.function_name_set("g");
g.optimize("no_compare_op");
std::string c_type = "double";
g.to_csrc(std::cout, c_type);
}The output: // includes
# include <stddef.h>
# include <math.h>
// typedefs
typedef double float_point_t;
// externals
// azmul
static float_point_t azmul(float_point_t x, float_point_t y)
{ if( x == 0.0 ) return 0.0;
return x * y;
}
// sign
static float_point_t sign(float_point_t x)
{ if( x > 0.0 ) return 1.0;
if( x == 0.0 ) return 0.0;
return -1.0;
}
// This JIT function
__declspec(dllexport) int __cdecl cppad_jit_g(
size_t nx ,
const float_point_t* x ,
size_t ny ,
float_point_t* y ,
size_t* compare_change )
{ // begin function body
// declare variables
float_point_t v[7];
size_t i;
// check nx, ny
if( nx != 2) return 1;
if( ny != 1) return 2;
// initialize
v[0] = NAN; // const
// independent variables
// set v[1+i] for i = 0, ..., nx-1
for(i = 0; i < nx; ++i)
v[1+i] = x[i];
// constants
// set v[1+nx+i] for i = 0, ..., nc-1
// nc = 2
v[1+nx+0] = 4;
v[1+nx+1] = 0;
// result nodes
// set v[1+nx+nc+i] for i = 0, ..., n_result_node-1
// n_result_node = 2
v[5] = v[3] * v[2];
v[6] = v[5] + v[5];
// dependent variables
// set y[i] for i = 0, ny-1
y[0] = v[6];
return 0;
}The value of Could you give some advice? |
Beta Was this translation helpful? Give feedback.
Replies: 4 comments 5 replies
-
|
Could you post a fully complete working example please? |
Beta Was this translation helpful? Give feedback.
-
|
What is the problem you encounter when compiling the code above? |
Beta Was this translation helpful? Give feedback.
-
|
I think this is what you want, please check the results and see that it works for you. |
Beta Was this translation helpful? Give feedback.
-
|
Thanks for the elegant solution! We need to declare This approach can be used to compute the Hessian as well. #include <iostream> // standard input/output
#include <cppad/cppad.hpp> // the CppAD package
// main program
int main(void)
{
using CppAD::AD;
using CppAD::ADFun;
using CppAD::Independent;
std::string c_type = "double";
size_t n = 1;
size_t m = 1;
size_t k = 1;
CPPAD_TESTVECTOR(AD<double>) ax(n), ay(m), ap(k);
Independent(ax, ap);
ay[0] = ap[0] * ax[0] * ax[0];
ADFun<double> f;
f.Dependent(ax, ay);
f.function_name_set("f");
f.optimize("no_compare_op");
f.to_csrc(std::cout, c_type);
ADFun<AD<double>, double> af = f.base2ad();
CPPAD_TESTVECTOR(AD<double>) axp(n + k);
{
Independent(axp);
for (size_t j = 0; j < n; ++j)
ax[j] = axp[j];
for (size_t j = 0; j < k; ++j)
ap[j] = axp[n + j];
CPPAD_TESTVECTOR(AD<double>) az(m * n);
af.new_dynamic(ap);
az = af.Jacobian(ax);
ADFun<double> g;
g.Dependent(axp, az);
g.function_name_set("jacobi");
g.optimize("no_compare_op");
g.to_csrc(std::cout, c_type);
}
{
Independent(axp);
for (size_t j = 0; j < n; ++j)
ax[j] = axp[j];
for (size_t j = 0; j < k; ++j)
ap[j] = axp[n + j];
CPPAD_TESTVECTOR(AD<double>) az(n * n);
af.new_dynamic(ap);
az = af.Hessian(ax, 0);
ADFun<double> h;
h.Dependent(axp, az);
h.function_name_set("hessian");
h.optimize("no_compare_op");
h.to_csrc(std::cout, c_type);
}
} |
Beta Was this translation helpful? Give feedback.
-
|
The sparse version works as well. I post my code here for reference: #include <iostream> // standard input/output
#include <cppad/cppad.hpp> // the CppAD package
// main program
int main(void)
{
using CppAD::AD;
using CppAD::ADFun;
using CppAD::Independent;
using CppAD::sparse_rc;
using CppAD::sparse_rcv;
typedef CPPAD_TESTVECTOR(double) d_vector;
typedef CPPAD_TESTVECTOR(AD<double>) ad_vector;
typedef CPPAD_TESTVECTOR(size_t) s_vector;
std::string c_type = "double";
size_t n = 2;
size_t m = 2;
size_t k = 1;
CPPAD_TESTVECTOR(AD<double>) ax(n), ay(m), ap(k);
Independent(ax, ap);
ay[0] = ap[0] * ax[0] * ax[0];
ay[1] = ax[1] * ax[1] + ap[0] * ax[1];
ADFun<double> f;
f.Dependent(ax, ay);
f.function_name_set("f");
f.optimize("no_compare_op");
f.to_csrc(std::cout, c_type);
// n by n identity matrix sparsity
sparse_rc<s_vector> pattern_in;
pattern_in.resize(n, n, n);
for (size_t k = 0; k < n; k++)
pattern_in.set(k, k, k);
// sparsity for J(x)
bool transpose = false;
bool dependency = false;
bool internal_bool = true;
sparse_rc<s_vector> pattern_jac;
f.for_jac_sparsity(pattern_in, transpose, dependency, internal_bool, pattern_jac);
// include all x components in sparsity pattern
CPPAD_TESTVECTOR(bool) select_domain(n);
for (size_t j = 0; j < n; j++)
select_domain[j] = true;
// compute sparsity pattern for H(x) = F_1''(x) + F_2''(x)
CPPAD_TESTVECTOR(bool) select_range(m);
select_range[0] = true;
select_range[1] = true;
internal_bool = true;
sparse_rc<s_vector> pattern_hes;
f.for_hes_sparsity(select_domain, select_range, internal_bool, pattern_hes);
ADFun<AD<double>, double> af = f.base2ad();
CPPAD_TESTVECTOR(AD<double>) axp(n + k);
{
Independent(axp);
for (size_t j = 0; j < n; ++j)
ax[j] = axp[j];
for (size_t j = 0; j < k; ++j)
ap[j] = axp[n + j];
af.new_dynamic(ap);
// compute entire forward mode Jacobian
sparse_rcv<s_vector, ad_vector> subset(pattern_jac);
CppAD::sparse_jac_work work;
std::string coloring = "cppad";
size_t group_max = 10;
size_t n_color = af.sparse_jac_for(group_max, ax, subset, pattern_jac, coloring, work);
ADFun<double> g;
g.Dependent(axp, subset.val());
g.function_name_set("jacobi");
g.optimize("no_compare_op");
g.to_csrc(std::cout, c_type);
}
{
Independent(axp);
for (size_t j = 0; j < n; ++j)
ax[j] = axp[j];
for (size_t j = 0; j < k; ++j)
ap[j] = axp[n + j];
af.new_dynamic(ap);
// compute entire sparse Hessian (really only need lower triangle)
CppAD::sparse_rcv<s_vector, ad_vector> subset(pattern_hes);
CppAD::sparse_hes_work work;
std::string coloring = "cppad.symmetric";
ad_vector w(n);
for (size_t j = 0; j < n; j++)
w[j] = 1.0;
size_t n_sweep = af.sparse_hes(ax, w, subset, pattern_hes, coloring, work);
ADFun<double> h;
h.Dependent(axp, subset.val());
h.function_name_set("hessian");
h.optimize("no_compare_op");
h.to_csrc(std::cout, c_type);
}
} |
Beta Was this translation helpful? Give feedback.
-
Can I just check if these code snippets are correct, are you sure these assignments, and the others, should not be axp[i] = ax[i]? Also you never set the value of ap to anything? Trying to learn form these posted examples so just wanted to make sure. I'm also struggling with why the arguments to Independent and Dependent are always different, could you annotate the examples with a few comments saying what's going on? Finally, really thrown by new_dynamic, could you explain what you're doing in the context of these examples and the use of new_dynamic with Dependant, Independent? Thanks |
Beta Was this translation helpful? Give feedback.
-
Does this make the code above clear to you ? |
Beta Was this translation helpful? Give feedback.
-
|
Hi, thanks for the explanation - I'll be honest, I'm kind of struggling (my fault). Would you be able to try to answer my other questions above? Thanks |
Beta Was this translation helpful? Give feedback.
-
|
@metab0t @a-jp I forgot that both the dynamic parameters and the variables are in the jit function. The dynamic parameters come first. I suggest you look at the file dynamic.c (created by this example). Also, see the heading 03-06 on |
Beta Was this translation helpful? Give feedback.
I think this is what you want, please check the results and see that it works for you.