vector_sparse.hpp
Go to the documentation of this file.
1 /*!
2  * \brief Sparse vector class
3  *
4  * \author O. Krause
5  * \date 2013
6  *
7  *
8  * \par Copyright 1995-2015 Shark Development Team
9  *
10  * <BR><HR>
11  * This file is part of Shark.
12  * <http://image.diku.dk/shark/>
13  *
14  * Shark is free software: you can redistribute it and/or modify
15  * it under the terms of the GNU Lesser General Public License as published
16  * by the Free Software Foundation, either version 3 of the License, or
17  * (at your option) any later version.
18  *
19  * Shark is distributed in the hope that it will be useful,
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22  * GNU Lesser General Public License for more details.
23  *
24  * You should have received a copy of the GNU Lesser General Public License
25  * along with Shark. If not, see <http://www.gnu.org/licenses/>.
26  *
27  */
28 #ifndef SHARK_LINALG_BLAS_VECTOR_SPARSE_HPP
29 #define SHARK_LINALG_BLAS_VECTOR_SPARSE_HPP
30 
31 #include "vector_proxy.hpp"
32 #include <vector>
33 
34 namespace shark {
35 namespace blas {
36 
37 /** \brief Compressed array based sparse vector
38  *
39  * a sparse vector of values of type T of variable size. The non zero values are stored as
40  * two seperate arrays: an index array and a value array. The index array is always sorted
41  * and there is at most one entry for each index. Inserting an element can be time consuming.
42  * If the vector contains a few zero entries, then it is better to have a normal vector.
43  * If the vector has a very high dimension with a few non-zero values, then this vector is
44  * very memory efficient (at the cost of a few more computations).
45  *
46  * For a \f$n\f$-dimensional compressed vector and \f$0 \leq i < n\f$ the non-zero elements
47  * \f$v_i\f$ are mapped to consecutive elements of the index and value container, i.e. for
48  * elements \f$k = v_{i_1}\f$ and \f$k + 1 = v_{i_2}\f$ of these containers holds \f$i_1 < i_2\f$.
49  *
50  * Supported parameters for the adapted array (indices and values) are \c unbounded_array<> ,
51  * \c bounded_array<> and \c std::vector<>.
52  *
53  * \tparam T the type of object stored in the vector (like double, float, complex, etc...)
54  * \tparam I the indices stored in the vector
55  */
56 template<class T, class I>
57 class compressed_vector:public vector_container<compressed_vector<T, I> > {
58 
59  typedef T& true_reference;
60  typedef compressed_vector<T, I> self_type;
61 public:
62 
63  typedef std::size_t size_type;
64  typedef std::ptrdiff_t difference_type;
65  typedef T value_type;
66  typedef value_type scalar_type;
67  typedef T *pointer;
68  typedef const T *const_pointer;
69  typedef const T& const_reference;
70 
71  typedef I index_type;
72  typedef index_type const* const_index_pointer;
73  typedef index_type index_pointer;
74 
75  class reference {
76  private:
77 
78  const_reference value()const {
79  return const_cast<self_type const&>(m_vector)(m_i);
80  }
81  value_type& ref() const {
82  //find position of the index in the array
83  index_type const* start = m_vector.indices();
84  index_type const* end = start + m_vector.nnz();
85  index_type const *pos = std::lower_bound(start,end,m_i);
86 
87  if (pos != end&& *pos == m_i)
88  return m_vector.m_values[pos-start];
89  else {
90  //create iterator to the insertion position and insert new element
91  iterator posIter(m_vector.values(),m_vector.indices(),pos-start);
92  return *m_vector.set_element(posIter, m_i, m_vector.m_zero);
93  }
94  }
95 
96  public:
97  // Construction and destruction
98  reference(self_type& m, size_type i):
99  m_vector(m), m_i(i) {}
100 
101  // Assignment
102  value_type& operator = (value_type d)const {
103  return ref()=d;
104  }
105 
106  value_type& operator=(reference const& v ){
107  return ref() = v.value();
108  }
109 
110  value_type& operator += (value_type d)const {
111  return ref()+=d;
112  }
113  value_type& operator -= (value_type d)const {
114  return ref()-=d;
115  }
116  value_type& operator *= (value_type d)const {
117  return ref()*=d;
118  }
119  value_type& operator /= (value_type d)const {
120  return ref()/=d;
121  }
122 
123  // Comparison
124  bool operator == (value_type d) const {
125  return value() == d;
126  }
127  bool operator != (value_type d) const {
128  return value() != d;
129  }
130 
131  operator const_reference() const{
132  return value();
133  }
134  private:
135  self_type& m_vector;
136  size_type m_i;
137  };
138 
141  typedef sparse_tag storage_category;
142  typedef elementwise_tag evaluation_category;
143 
144  // Construction and destruction
145  compressed_vector():m_size(0), m_nnz(0),m_indices(1,0),m_zero(0){}
146  explicit compressed_vector(size_type size, value_type value = value_type(), size_type non_zeros = 0)
147  :m_size(size), m_nnz(0), m_indices(non_zeros,0), m_values(non_zeros),m_zero(0){}
148  template<class AE>
149  compressed_vector(vector_expression<AE> const& ae, size_type non_zeros = 0)
150  :m_size(ae().size()), m_nnz(0), m_indices(non_zeros,0), m_values(non_zeros),m_zero(0)
151  {
152  assign(*this, ae);
153  }
154 
155  // Accessors
156  size_type size() const {
157  return m_size;
158  }
159  size_type nnz_capacity() const {
160  return m_indices.size();
161  }
162  size_type nnz() const {
163  return m_nnz;
164  }
165 
166  // Storage accessors
167  void set_filled(size_type filled) {
168  SIZE_CHECK(filled <= nnz_capacity());
169  m_nnz = filled;
170  }
171 
172  index_type const* indices() const{
173  if(nnz_capacity() == 0)
174  return 0;
175  return& m_indices[0];
176  }
177  index_type* indices(){
178  if(nnz_capacity() == 0)
179  return 0;
180  return& m_indices[0];
181  }
182  value_type const* values() const {
183  if(nnz_capacity() == 0)
184  return 0;
185  return& m_values[0];
186  }
187  value_type* values(){
188  if(nnz_capacity() == 0)
189  return 0;
190  return& m_values[0];
191  }
192 
193  void resize(size_type size) {
194  m_size = size;
195  m_nnz = 0;
196  }
197  void reserve(size_type non_zeros) {
198  if(non_zeros <= nnz_capacity()) return;
199  non_zeros = std::min(size(),non_zeros);
200  m_indices.resize(non_zeros);
201  m_values.resize(non_zeros);
202  }
203 
204  // Element access
205  const_reference operator()(size_type i) const {
206  SIZE_CHECK(i < m_size);
207  std::size_t pos = lower_bound(i);
208  if (pos == nnz() || m_indices[pos] != i)
209  return m_zero;
210  return m_values [pos];
211  }
212  reference operator()(size_type i) {
213  return reference(*this,i);
214  }
215 
216 
217  const_reference operator [](size_type i) const {
218  return (*this)(i);
219  }
220  reference operator [](size_type i) {
221  return (*this)(i);
222  }
223 
224  // Zeroing
225  void clear() {
226  m_nnz = 0;
227  }
228 
229  // Assignment
231  m_size = v.m_size;
232  m_nnz = v.m_nnz;
233  m_indices = v.m_indices;
234  m_values = v.m_values;
235  return *this;
236  }
237  template<class C> // Container assignment without temporary
239  resize(v().size(), false);
240  assign(*this, v);
241  return *this;
242  }
243  template<class AE>
245  self_type temporary(ae, nnz_capacity());
246  swap(temporary);
247  return *this;
248  }
249 
250  // Swapping
252  std::swap(m_size, v.m_size);
253  std::swap(m_nnz, v.m_nnz);
254  m_indices.swap(v.m_indices);
255  m_values.swap(v.m_values);
256  }
257 
258  friend void swap(compressed_vector& v1, compressed_vector& v2){
259  v1.swap(v2);
260  }
261 
262  // Iterator types
263  typedef compressed_storage_iterator<value_type const, index_type const> const_iterator;
264  typedef compressed_storage_iterator<value_type, index_type const> iterator;
265 
266  const_iterator begin() const {
267  return const_iterator(values(),indices(),0);
268  }
269 
270  const_iterator end() const {
271  return const_iterator(values(),indices(),nnz());
272  }
273 
274  iterator begin() {
275  return iterator(values(),indices(),0);
276  }
277 
278  iterator end() {
279  return iterator(values(),indices(),nnz());
280  }
281 
282  // Element assignment
283  iterator set_element(iterator pos, size_type index, value_type value) {
284  RANGE_CHECK(size_type(pos - begin()) <=m_size);
285 
286  if(pos != end() && pos.index() == index){
287  *pos = value;
288  return pos;
289  }
290  //get position of the new element in the array.
291  difference_type arrayPos = pos - begin();
292  if (m_nnz <= nnz_capacity())//reserve more space if needed, this invalidates pos.
293  reserve(std::max<std::size_t>(2 * nnz_capacity(),1));
294 
295  //copy the remaining elements to make space for the new ones
296  std::copy_backward(
297  m_values.begin()+arrayPos,m_values.begin() + m_nnz , m_values.begin() + m_nnz +1
298  );
299  std::copy_backward(
300  m_indices.begin()+arrayPos,m_indices.begin() + m_nnz , m_indices.begin() + m_nnz +1
301  );
302  //insert new element
303  m_values[arrayPos] = value;
304  m_indices[arrayPos] = index;
305  ++m_nnz;
306 
307 
308  //return new iterator to the inserted element.
309  return iterator(values(),indices(),arrayPos);
310  }
311 
312  iterator clear_range(iterator start, iterator end) {
313  //get position of the elements in the array.
314  difference_type startPos = start - begin();
315  difference_type endPos = end - begin();
316 
317  //remove the elements in the range
318  std::copy(
319  m_values.begin()+endPos,m_values.begin() + m_nnz, m_values.begin() + startPos
320  );
321  std::copy(
322  m_indices.begin()+endPos,m_indices.begin() + m_nnz , m_indices.begin() + startPos
323  );
324  m_nnz -= endPos - startPos;
325  //return new iterator to the next element
326  return iterator(values(),indices(), startPos);
327  }
328 
329  iterator clear_element(iterator pos){
330  //get position of the element in the array.
331  difference_type arrayPos = pos - begin();
332  if(arrayPos == m_nnz-1){//last element
333  --m_nnz;
334  return end();
335  }
336 
337  std::copy(
338  m_values.begin()+arrayPos+1,m_values.begin() + m_nnz , m_values.begin() + arrayPos
339  );
340  std::copy(
341  m_indices.begin()+arrayPos+1,m_indices.begin() + m_nnz , m_indices.begin() + arrayPos
342  );
343  //return new iterator to the next element
344  return iterator(values(),indices(),arrayPos);
345  }
346 
347  // Serialization
348  template<class Archive>
349  void serialize(Archive& ar, const unsigned int /* file_version */) {
350  boost::serialization::collection_size_type s(m_size);
351  ar & boost::serialization::make_nvp("size",s);
352  if (Archive::is_loading::value) {
353  m_size = s;
354  }
355  // ISSUE: m_indices and m_values are undefined between m_nnz and capacity (trouble with 'nan'-values)
356  ar & boost::serialization::make_nvp("nnz", m_nnz);
357  ar & boost::serialization::make_nvp("indices", m_indices);
358  ar & boost::serialization::make_nvp("values", m_values);
359  }
360 
361 private:
362  std::size_t lower_bound( index_type t)const{
363  index_type const* begin = indices();
364  index_type const* end = indices()+nnz();
365  return std::lower_bound(begin, end, t)-begin;
366  }
367 
368  size_type m_size;
369  size_type m_nnz;
370  std::vector<index_type> m_indices;
371  std::vector<value_type> m_values;
372  value_type m_zero;
373 };
374 
375 template<class T>
376 struct vector_temporary_type<T,sparse_bidirectional_iterator_tag>{
377  typedef compressed_vector<T> type;
378 };
379 
380 template<class T,class I>
381 struct const_expression<compressed_vector<T,I> >{
382  typedef compressed_vector<T,I> const type;
383 };
384 template<class T,class I>
385 struct const_expression<compressed_vector<T,I> const>{
386  typedef compressed_vector<T,I> const type;
387 };
388 
389 }}
390 
391 #endif