1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
|
// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
// SPDX-FileContributor: 2003-22 Bradley M. Bell
// ----------------------------------------------------------------------------
/*
{xrst_begin json_sparse.cpp}
Json Representation of a Sparse Matrix: Example and Test
########################################################
Discussion
**********
The example using a CppAD Json to represent the sparse matrix
.. math::
\partial_x f(x, p) = \left( \begin{array}{ccc}
p_0 & 0 & 0 \\
0 & p_1 & 0 \\
0 & 0 & c_0
\end{array} \right)
where :math:`c_0` is the constant 3.
Source Code
***********
{xrst_literal
// BEGIN C++
// END C++
}
{xrst_end json_sparse.cpp}
*/
// BEGIN C++
# include <cppad/cppad.hpp>
bool sparse(void)
{ bool ok = true;
using CppAD::vector;
typedef vector<size_t> s_vector;
typedef vector<double> d_vector;
//
// An AD graph example
// node_1 : p[0]
// node_2 : p[1]
// node_3 : x[0]
// node_4 : x[1]
// node_5 : x[2]
// node_6 : c[0]
// node_7 : p[0] * x[0]
// node_8 : p[1] * x[1]
// node_9 : c[0] * x[2]
// y[0] = p[0] * x[0]
// y[1] = p[1] * x[1]
// y[2] = c[0] * x[0]
// use single quote to avoid having to escape double quote
std::string json =
"{\n"
" 'function_name' : 'sparse example',\n"
" 'op_define_vec' : [ 1, [\n"
" { 'op_code':1, 'name':'mul', 'n_arg':2 } ]\n"
" ],\n"
" 'n_dynamic_ind' : 2,\n"
" 'n_variable_ind' : 3,\n"
" 'constant_vec' : [ 1, [ 3.0 ] ],\n"
" 'op_usage_vec' : [ 3, [\n"
" [ 1, 1, 3 ] , \n"
" [ 1, 2, 4 ] , \n"
" [ 1, 6, 5 ] ] \n"
" ],\n"
" 'dependent_vec' : [ 3, [7, 8, 9] ] \n"
"}\n";
// Convert the single quote to double quote
for(size_t i = 0; i < json.size(); ++i)
if( json[i] == '\'' ) json[i] = '"';
//
CppAD::ADFun<double> fun;
fun.from_json(json);
ok &= fun.Domain() == 3;
ok &= fun.Range() == 3;
ok &= fun.size_dyn_ind() == 2;
//
// Point at which we are going to evaluate the Jacobian of f(x, p).
// The Jacobian is constant w.r.t. x, so the value of x does not matter.
vector<double> p(2), x(3);
p[0] = 1.0;
p[1] = 2.0;
for(size_t j = 0; j < x.size(); ++j)
x[j] = 0.0;
//
// set the dynamic parameters
fun.new_dynamic(p);
//
// compute the sparsity pattern for f_x(x, p)
vector<bool> select_domain(3);
vector<bool> select_range(3);
for(size_t i = 0; i < 3; ++i)
{ select_domain[i] = true;
select_range[i] = true;
}
bool transpose = false;
CppAD::sparse_rc<s_vector> pattern;
fun.subgraph_sparsity(select_domain, select_range, transpose, pattern);
//
// compute the entire Jacobian
CppAD::sparse_rcv<s_vector, d_vector> subset(pattern);
fun.subgraph_jac_rev(x, subset);
//
// information in the sparse Jacobian
const s_vector& row( subset.row() );
const s_vector& col( subset.col() );
const d_vector& val( subset.val() );
size_t nnz = subset.nnz();
s_vector row_major = subset.row_major();
//
// check number of non-zero elements in sparse matrix
ok &= nnz == 3;
//
// check first element of matrix (in row major order)
size_t k = row_major[0];
ok &= row[k] == 0;
ok &= col[k] == 0;
ok &= val[k] == p[0];
//
// check second element of matrix
k = row_major[1];
ok &= row[k] == 1;
ok &= col[k] == 1;
ok &= val[k] == p[1];
//
// check third element of matrix
k = row_major[2];
ok &= row[k] == 2;
ok &= col[k] == 2;
ok &= val[k] == 3.0;
//
return ok;
}
// END C++
|
