1 /*
2 * Copyright © 2004 Ondra Kamenik
3 * Copyright © 2019 Dynare Team
4 *
5 * This file is part of Dynare.
6 *
7 * Dynare is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * Dynare is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with Dynare. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21 #include "sparse_tensor.hh"
22 #include "fs_tensor.hh"
23 #include "tl_exception.hh"
24
25 #include <iostream>
26 #include <iomanip>
27 #include <cmath>
28
29 /* This is straightforward. Before we insert anything, we do a few
30 checks. Then we reset ‘first_nz_row’ and ‘last_nz_row’ if necessary. */
31
32 void
insert(IntSequence key,int r,double c)33 SparseTensor::insert(IntSequence key, int r, double c)
34 {
35 TL_RAISE_IF(r < 0 || r >= nr,
36 "Row number out of dimension of tensor in SparseTensor::insert");
37 TL_RAISE_IF(key.size() != dimen(),
38 "Wrong length of key in SparseTensor::insert");
39 TL_RAISE_IF(!std::isfinite(c),
40 "Insertion of non-finite value in SparseTensor::insert");
41
42 auto first_pos = m.lower_bound(key);
43
44 // check that pair ‘key’ and ‘r’ is unique
45 auto last_pos = m.upper_bound(key);
46 for (auto it = first_pos; it != last_pos; ++it)
47 TL_RAISE_IF(it->second.first == r, "Duplicate <key, r> insertion in SparseTensor::insert");
48
49 m.emplace_hint(first_pos, std::move(key), std::make_pair(r, c));
50 if (first_nz_row > r)
51 first_nz_row = r;
52 if (last_nz_row < r)
53 last_nz_row = r;
54 }
55
56 /* This returns true if all items are finite (not NaN nor ∞). */
57
58 bool
isFinite() const59 SparseTensor::isFinite() const
60 {
61 bool res = true;
62 auto run = m.begin();
63 while (res && run != m.end())
64 {
65 if (!std::isfinite(run->second.second))
66 res = false;
67 ++run;
68 }
69 return res;
70 }
71
72 /* This returns a ratio of a number of non-zero columns in folded
73 tensor to the total number of columns. */
74
75 double
getFoldIndexFillFactor() const76 SparseTensor::getFoldIndexFillFactor() const
77 {
78 int cnt = 0;
79 auto start_col = m.begin();
80 while (start_col != m.end())
81 {
82 cnt++;
83 const IntSequence &key = start_col->first;
84 start_col = m.upper_bound(key);
85 }
86
87 return static_cast<double>(cnt)/ncols();
88 }
89
90 /* This returns a ratio of a number of non-zero columns in unfolded
91 tensor to the total number of columns. */
92
93 double
getUnfoldIndexFillFactor() const94 SparseTensor::getUnfoldIndexFillFactor() const
95 {
96 int cnt = 0;
97 auto start_col = m.begin();
98 while (start_col != m.end())
99 {
100 const IntSequence &key = start_col->first;
101 cnt += key.getSymmetry().noverseq();
102 start_col = m.upper_bound(key);
103 }
104
105 return static_cast<double>(cnt)/ncols();
106 }
107
108 /* This prints the fill factor and all items. */
109
110 void
print() const111 SparseTensor::print() const
112 {
113 std::cout << "Fill: "
114 << std::fixed << std::setprecision(2) << 100*getFillFactor()
115 << std::setprecision(6) << std::defaultfloat << " %\n";
116 auto start_col = m.begin();
117 while (start_col != m.end())
118 {
119 const IntSequence &key = start_col->first;
120 std::cout << "Column: ";
121 key.print();
122 auto end_col = m.upper_bound(key);
123 int cnt = 1;
124 for (auto run = start_col; run != end_col; ++run, cnt++)
125 {
126 if (cnt % 7 == 0)
127 std::cout << "\n";
128 std::cout << run->second.first << '(' << run->second.second << ") ";
129 }
130 std::cout << "\n";
131 start_col = end_col;
132 }
133 }
134
FSSparseTensor(int d,int nvar,int r)135 FSSparseTensor::FSSparseTensor(int d, int nvar, int r)
136 : SparseTensor(d, r, FFSTensor::calcMaxOffset(nvar, d)),
137 nv(nvar), sym{d}
138 {
139 }
140
141 void
insert(IntSequence key,int r,double c)142 FSSparseTensor::insert(IntSequence key, int r, double c)
143 {
144 TL_RAISE_IF(!key.isSorted(),
145 "Key is not sorted in FSSparseTensor::insert");
146 TL_RAISE_IF(key[key.size()-1] >= nv || key[0] < 0,
147 "Wrong value of the key in FSSparseTensor::insert");
148 SparseTensor::insert(std::move(key), r, c);
149 }
150
151 /* We go through the tensor ‘t’ which is supposed to have single
152 column. If the item of ‘t’ is nonzero, we make a key by sorting the
153 index, and then we go through all items having the same key (it is its
154 column), obtain the row number and the element, and do the
155 multiplication.
156
157 The test for non-zero is ‘a != 0.0’, since there will be items which
158 are exact zeros.
159
160 I have also tried to make the loop through the sparse tensor outer, and
161 find index of tensor ‘t’ within the loop. Surprisingly, it is little
162 slower (for monomial tests with probability of zeros equal 0.3). But
163 everything depends how filled is the sparse tensor. */
164
165 void
multColumnAndAdd(const Tensor & t,Vector & v) const166 FSSparseTensor::multColumnAndAdd(const Tensor &t, Vector &v) const
167 {
168 // check compatibility of input parameters
169 TL_RAISE_IF(v.length() != nrows(),
170 "Wrong size of output vector in FSSparseTensor::multColumnAndAdd");
171 TL_RAISE_IF(t.dimen() != dimen(),
172 "Wrong dimension of tensor in FSSparseTensor::multColumnAndAdd");
173 TL_RAISE_IF(t.ncols() != 1,
174 "The input tensor is not single-column in FSSparseTensor::multColumnAndAdd");
175
176 for (Tensor::index it = t.begin(); it != t.end(); ++it)
177 {
178 int ind = *it;
179 double a = t.get(ind, 0);
180 if (a != 0.0)
181 {
182 IntSequence key(it.getCoor());
183 key.sort();
184
185 // check that ‘key’ is within the range
186 TL_RAISE_IF(key[0] < 0 || key[key.size()-1] >= nv,
187 "Wrong coordinates of index in FSSparseTensor::multColumnAndAdd");
188
189 auto first_pos = m.lower_bound(key);
190 auto last_pos = m.upper_bound(key);
191 for (auto cit = first_pos; cit != last_pos; ++cit)
192 {
193 int r = cit->second.first;
194 double c = cit->second.second;
195 v[r] += c * a;
196 }
197 }
198 }
199 }
200
201 void
print() const202 FSSparseTensor::print() const
203 {
204 std::cout << "FS Sparse tensor: dim=" << dim << ", nv=" << nv << ", (" << nr << 'x' << nc << ")\n";
205 SparseTensor::print();
206 }
207
208 // GSSparseTensor slicing constructor
209 /* This is the same as FGSTensor slicing constructor from FSSparseTensor. */
GSSparseTensor(const FSSparseTensor & t,const IntSequence & ss,const IntSequence & coor,TensorDimens td)210 GSSparseTensor::GSSparseTensor(const FSSparseTensor &t, const IntSequence &ss,
211 const IntSequence &coor, TensorDimens td)
212 : SparseTensor(td.dimen(), t.nrows(), td.calcFoldMaxOffset()),
213 tdims(std::move(td))
214 {
215 // set ‘lb’ and ‘ub’ to lower and upper bounds of slice indices
216 /* The same code is present in FGSTensor slicing constructor, see it for
217 details. */
218 IntSequence s_offsets(ss.size(), 0);
219 for (int i = 1; i < ss.size(); i++)
220 s_offsets[i] = s_offsets[i-1] + ss[i-1];
221
222 IntSequence lb(coor.size());
223 IntSequence ub(coor.size());
224 for (int i = 0; i < coor.size(); i++)
225 {
226 lb[i] = s_offsets[coor[i]];
227 ub[i] = s_offsets[coor[i]] + ss[coor[i]] - 1;
228 }
229
230 auto lbi = t.getMap().lower_bound(lb);
231 auto ubi = t.getMap().upper_bound(ub);
232 for (auto run = lbi; run != ubi; ++run)
233 {
234 if (lb.lessEq(run->first) && run->first.lessEq(ub))
235 {
236 IntSequence c(run->first);
237 c.add(-1, lb);
238 insert(c, run->second.first, run->second.second);
239 }
240 }
241
242 }
243
244 void
insert(IntSequence s,int r,double c)245 GSSparseTensor::insert(IntSequence s, int r, double c)
246 {
247 TL_RAISE_IF(!s.less(tdims.getNVX()),
248 "Wrong coordinates of index in GSSparseTensor::insert");
249 SparseTensor::insert(std::move(s), r, c);
250 }
251
252 void
print() const253 GSSparseTensor::print() const
254 {
255 std::cout << "GS Sparse tensor: (" << nr << 'x' << nc << ")\nSymmetry: ";
256 tdims.getSym().print();
257 std::cout << "NVS: ";
258 tdims.getNVS().print();
259 SparseTensor::print();
260 }
261