Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights 1 // Copyright (C) 2010, Guy Barrand. All rights reserved. 2 // See the file tools.license for terms. 2 // See the file tools.license for terms. 3 3 4 #ifndef tools_array 4 #ifndef tools_array 5 #define tools_array 5 #define tools_array 6 6 7 #ifdef TOOLS_MEM 7 #ifdef TOOLS_MEM 8 #include "mem" 8 #include "mem" 9 #endif 9 #endif 10 10 11 #include <vector> 11 #include <vector> 12 #include <string> 12 #include <string> 13 13 14 namespace tools { 14 namespace tools { 15 15 16 // array handles an hyperparallelepiped of cel 16 // array handles an hyperparallelepiped of cells of class T. 17 17 18 template <class T> 18 template <class T> 19 class array { 19 class array { 20 public: 20 public: 21 static const std::string& s_class() { 21 static const std::string& s_class() { 22 static const std::string s_v("tools::array 22 static const std::string s_v("tools::array"); 23 return s_v; 23 return s_v; 24 } 24 } 25 public: 25 public: 26 typedef typename std::vector< std::pair<unsi 26 typedef typename std::vector< std::pair<unsigned int,unsigned int> > cut_t; 27 typedef typename std::vector<unsigned int> u 27 typedef typename std::vector<unsigned int> uints_t; 28 typedef typename std::vector<T>::iterator ve 28 typedef typename std::vector<T>::iterator vec_it_t; 29 typedef typename std::vector<T>::const_itera 29 typedef typename std::vector<T>::const_iterator cons_vec_it_t; 30 public: 30 public: 31 array() { 31 array() { 32 #ifdef TOOLS_MEM 32 #ifdef TOOLS_MEM 33 mem::increment(s_class().c_str()); 33 mem::increment(s_class().c_str()); 34 #endif 34 #endif 35 } 35 } 36 array(const uints_t& a_orders) { 36 array(const uints_t& a_orders) { 37 #ifdef TOOLS_MEM 37 #ifdef TOOLS_MEM 38 mem::increment(s_class().c_str()); 38 mem::increment(s_class().c_str()); 39 #endif 39 #endif 40 configure(a_orders); 40 configure(a_orders); 41 } 41 } 42 array(unsigned int a_dimension,unsigned int 42 array(unsigned int a_dimension,unsigned int a_order) { 43 #ifdef TOOLS_MEM 43 #ifdef TOOLS_MEM 44 mem::increment(s_class().c_str()); 44 mem::increment(s_class().c_str()); 45 #endif 45 #endif 46 // A hypercube of dimension "a_dimension" 46 // A hypercube of dimension "a_dimension" and size "a_order". 47 uints_t _orders(a_dimension); 47 uints_t _orders(a_dimension); 48 for(unsigned int index=0;index<a_dimension 48 for(unsigned int index=0;index<a_dimension;index++) 49 _orders[index] = a_order; 49 _orders[index] = a_order; 50 configure(_orders); 50 configure(_orders); 51 } 51 } 52 virtual ~array() { 52 virtual ~array() { 53 #ifdef TOOLS_MEM 53 #ifdef TOOLS_MEM 54 mem::decrement(s_class().c_str()); 54 mem::decrement(s_class().c_str()); 55 #endif 55 #endif 56 } 56 } 57 public: 57 public: 58 array(const array& a_from) 58 array(const array& a_from) 59 :m_orders(a_from.m_orders) 59 :m_orders(a_from.m_orders) 60 ,m_offsets(a_from.m_offsets) 60 ,m_offsets(a_from.m_offsets) 61 ,m_vector(a_from.m_vector) 61 ,m_vector(a_from.m_vector) 62 ,m_is(a_from.m_is){ 62 ,m_is(a_from.m_is){ 63 #ifdef TOOLS_MEM 63 #ifdef TOOLS_MEM 64 mem::increment(s_class().c_str()); 64 mem::increment(s_class().c_str()); 65 #endif 65 #endif 66 } 66 } 67 array& operator=(const array& a_from) { 67 array& operator=(const array& a_from) { 68 m_orders = a_from.m_orders; 68 m_orders = a_from.m_orders; 69 m_offsets = a_from.m_offsets; 69 m_offsets = a_from.m_offsets; 70 m_vector = a_from.m_vector; 70 m_vector = a_from.m_vector; 71 m_is = a_from.m_is; 71 m_is = a_from.m_is; 72 return *this; 72 return *this; 73 } 73 } 74 public: //operators: 74 public: //operators: 75 array& operator*=(const T& a_T) {multiply(a_ 75 array& operator*=(const T& a_T) {multiply(a_T);return *this;} 76 bool operator==(const array& a_array) const 76 bool operator==(const array& a_array) const { 77 return equal(a_array); 77 return equal(a_array); 78 } 78 } 79 bool operator!=(const array& a_array) const 79 bool operator!=(const array& a_array) const { 80 return !operator==(a_array); 80 return !operator==(a_array); 81 } 81 } 82 array operator*(const T& a_T) const { 82 array operator*(const T& a_T) const { 83 array a(*this); 83 array a(*this); 84 a.multiply(a_T); 84 a.multiply(a_T); 85 return a; 85 return a; 86 } 86 } 87 // the below would need exception to do it p 87 // the below would need exception to do it properly. 88 //array operator+(const array& a_array) cons 88 //array operator+(const array& a_array) const { 89 // array a(*this); 89 // array a(*this); 90 // if(!a.add(a_array)) {} //throw 90 // if(!a.add(a_array)) {} //throw 91 // return a; 91 // return a; 92 //} 92 //} 93 //array& operator+=(const array& a_array) { 93 //array& operator+=(const array& a_array) { 94 // if(!add(a_array)) {} //throw 94 // if(!add(a_array)) {} //throw 95 // return *this; 95 // return *this; 96 //} 96 //} 97 public: 97 public: 98 void copy(const array& a_from) { 98 void copy(const array& a_from) { 99 m_orders = a_from.m_orders; 99 m_orders = a_from.m_orders; 100 m_offsets = a_from.m_offsets; 100 m_offsets = a_from.m_offsets; 101 m_vector = a_from.m_vector; 101 m_vector = a_from.m_vector; 102 m_is = a_from.m_is; 102 m_is = a_from.m_is; 103 } 103 } 104 void clear() { 104 void clear() { 105 m_orders.clear(); 105 m_orders.clear(); 106 m_offsets.clear(); 106 m_offsets.clear(); 107 m_vector.clear(); 107 m_vector.clear(); 108 m_is.clear(); 108 m_is.clear(); 109 } 109 } 110 bool configure(const uints_t& a_orders) { 110 bool configure(const uints_t& a_orders) { 111 m_orders = a_orders; 111 m_orders = a_orders; 112 size_t dim = m_orders.size(); 112 size_t dim = m_orders.size(); 113 if(dim==0) { 113 if(dim==0) { 114 clear(); 114 clear(); 115 return false; 115 return false; 116 } 116 } 117 unsigned int _size = 1; 117 unsigned int _size = 1; 118 for(size_t index=0;index<dim;index++) { 118 for(size_t index=0;index<dim;index++) { 119 //if(m_orders[index]<0) { 119 //if(m_orders[index]<0) { 120 //clear(); 120 //clear(); 121 //return false; 121 //return false; 122 //} 122 //} 123 _size *= m_orders[index]; 123 _size *= m_orders[index]; 124 } 124 } 125 m_vector.resize(_size); 125 m_vector.resize(_size); 126 m_vector.assign(_size,zero()); 126 m_vector.assign(_size,zero()); 127 m_offsets.resize(dim,0); 127 m_offsets.resize(dim,0); 128 m_offsets[0] = 1; 128 m_offsets[0] = 1; 129 for(size_t iaxis=1;iaxis<dim;iaxis++) 129 for(size_t iaxis=1;iaxis<dim;iaxis++) 130 m_offsets[iaxis] = m_offsets[iaxis-1] * 130 m_offsets[iaxis] = m_offsets[iaxis-1] * m_orders[iaxis-1]; 131 m_is.resize(dim); 131 m_is.resize(dim); 132 m_is.assign(dim,0); 132 m_is.assign(dim,0); 133 return true; 133 return true; 134 } 134 } 135 size_t dimension() const { return m_orders.s 135 size_t dimension() const { return m_orders.size();} 136 const uints_t& orders() const { return m_ord 136 const uints_t& orders() const { return m_orders;} 137 size_t size() const {return m_vector.size(); 137 size_t size() const {return m_vector.size();} 138 bool set_value(const uints_t& a_is,const T& 138 bool set_value(const uints_t& a_is,const T& a_value) { 139 unsigned int off = 0; 139 unsigned int off = 0; 140 if(!offset(a_is,off)) return false; 140 if(!offset(a_is,off)) return false; 141 m_vector[off] = a_value; 141 m_vector[off] = a_value; 142 return true; 142 return true; 143 } 143 } 144 bool value(const uints_t& a_is,T& a_value) c 144 bool value(const uints_t& a_is,T& a_value) const { 145 unsigned int off = 0; 145 unsigned int off = 0; 146 if(!offset(a_is,off)) { 146 if(!offset(a_is,off)) { 147 a_value = 0; 147 a_value = 0; 148 return false; 148 return false; 149 } 149 } 150 a_value = m_vector[off]; 150 a_value = m_vector[off]; 151 return true; 151 return true; 152 } 152 } 153 T value_no_check(const uints_t& a_is) const 153 T value_no_check(const uints_t& a_is) const { //TOUCHY 154 unsigned int off = 0; 154 unsigned int off = 0; 155 if(!offset(a_is,off)) {} 155 if(!offset(a_is,off)) {} 156 return m_vector[off]; 156 return m_vector[off]; 157 } 157 } 158 void reset() { 158 void reset() { 159 for(vec_it_t it=m_vector.begin();it!=m_vec 159 for(vec_it_t it=m_vector.begin();it!=m_vector.end();++it) *it = 0; 160 } 160 } 161 const std::vector<T>& vector() const { retur 161 const std::vector<T>& vector() const { return m_vector;} 162 std::vector<T>& vector() { return m_vector;} 162 std::vector<T>& vector() { return m_vector;} 163 bool fill(const std::vector<T>& a_values,cut 163 bool fill(const std::vector<T>& a_values,cut_t* a_cut = 0) { 164 size_t dsize = a_values.size(); 164 size_t dsize = a_values.size(); 165 size_t di = 0; 165 size_t di = 0; 166 unsigned int index = 0; 166 unsigned int index = 0; 167 for(vec_it_t it=m_vector.begin();it!=m_vec 167 for(vec_it_t it=m_vector.begin();it!=m_vector.end();++it,index++) { 168 if(!a_cut || (a_cut && accept(index,*a_c 168 if(!a_cut || (a_cut && accept(index,*a_cut)) ) { 169 if(di>=dsize) return false; //a_values 169 if(di>=dsize) return false; //a_values exhausted too early 170 *it = a_values[di]; 170 *it = a_values[di]; 171 di++; 171 di++; 172 } 172 } 173 } 173 } 174 return true; 174 return true; 175 } 175 } 176 bool fill(unsigned int a_sz,const T* a_data, 176 bool fill(unsigned int a_sz,const T* a_data,cut_t* a_cut = 0) { 177 unsigned int di = 0; 177 unsigned int di = 0; 178 unsigned int index = 0; 178 unsigned int index = 0; 179 for(vec_it_t it=m_vector.begin();it!=m_vec 179 for(vec_it_t it=m_vector.begin();it!=m_vector.end();++it,index++) { 180 if(!a_cut || (a_cut && accept(index,*a_c 180 if(!a_cut || (a_cut && accept(index,*a_cut)) ) { 181 if(di>=a_sz) return false; //a_values 181 if(di>=a_sz) return false; //a_values exhausted too early 182 *it = a_data[di]; 182 *it = a_data[di]; 183 di++; 183 di++; 184 } 184 } 185 } 185 } 186 return true; 186 return true; 187 } 187 } 188 // Related to other array : 188 // Related to other array : 189 bool equal(const array& a_array) const { 189 bool equal(const array& a_array) const { 190 if(m_orders!=a_array.m_orders) return fals 190 if(m_orders!=a_array.m_orders) return false; 191 cons_vec_it_t it = m_vector.begin(); 191 cons_vec_it_t it = m_vector.begin(); 192 cons_vec_it_t ait = a_array.m_vector.begin 192 cons_vec_it_t ait = a_array.m_vector.begin(); 193 for(;it!=m_vector.end();++it,++ait) { 193 for(;it!=m_vector.end();++it,++ait) { 194 if((*it)!=(*ait)) return false; 194 if((*it)!=(*ait)) return false; 195 } 195 } 196 return true; 196 return true; 197 } 197 } 198 bool equal(const array& a_array,T aEpsilon) 198 bool equal(const array& a_array,T aEpsilon) const { 199 if(m_orders!=a_array.m_orders) return fals 199 if(m_orders!=a_array.m_orders) return false; 200 cons_vec_it_t it = m_vector.begin(); 200 cons_vec_it_t it = m_vector.begin(); 201 cons_vec_it_t ait = a_array.m_vector.begin 201 cons_vec_it_t ait = a_array.m_vector.begin(); 202 for(;it!=m_vector.end();++it,++ait) { 202 for(;it!=m_vector.end();++it,++ait) { 203 T diff = (*it) - (*ait); 203 T diff = (*it) - (*ait); 204 if(diff<0) diff *= -1; 204 if(diff<0) diff *= -1; 205 if(diff>=aEpsilon) return false; 205 if(diff>=aEpsilon) return false; 206 } 206 } 207 return true; 207 return true; 208 } 208 } 209 bool is_proportional(const array& a_array,T& 209 bool is_proportional(const array& a_array,T& a_factor) const { 210 // If true, then : a_array = a_factor * th 210 // If true, then : a_array = a_factor * this. 211 a_factor = zero(); 211 a_factor = zero(); 212 if(m_orders!=a_array.m_orders) return fals 212 if(m_orders!=a_array.m_orders) return false; 213 bool first = true; 213 bool first = true; 214 cons_vec_it_t it = m_vector.begin(); 214 cons_vec_it_t it = m_vector.begin(); 215 cons_vec_it_t ait = a_array.m_vector.begin 215 cons_vec_it_t ait = a_array.m_vector.begin(); 216 for(;it!=m_vector.end();++it,++ait) { 216 for(;it!=m_vector.end();++it,++ait) { 217 if( ((*it)==zero()) && ((*ait)==z 217 if( ((*it)==zero()) && ((*ait)==zero())) { 218 continue; 218 continue; 219 } else if( ((*it)!=zero()) && ((*ait)==z 219 } else if( ((*it)!=zero()) && ((*ait)==zero())) { 220 return false; 220 return false; 221 } else if( ((*it)==zero()) && ((*ait)!=z 221 } else if( ((*it)==zero()) && ((*ait)!=zero())) { 222 return false; 222 return false; 223 } else { 223 } else { 224 if(first) { 224 if(first) { 225 a_factor = (*ait)/(*it); 225 a_factor = (*ait)/(*it); 226 first = false; 226 first = false; 227 } else { 227 } else { 228 if((*ait)!=(*it)*a_factor) return fa 228 if((*ait)!=(*it)*a_factor) return false; 229 } 229 } 230 } 230 } 231 } 231 } 232 return true; 232 return true; 233 } 233 } 234 bool add(const array& a_array,cut_t* a_cut = 234 bool add(const array& a_array,cut_t* a_cut = 0) { 235 if(m_orders!=a_array.m_orders) return fals 235 if(m_orders!=a_array.m_orders) return false; 236 vec_it_t it = m_vector.begin(); 236 vec_it_t it = m_vector.begin(); 237 cons_vec_it_t ait = a_array.m_vector.begin 237 cons_vec_it_t ait = a_array.m_vector.begin(); 238 unsigned int index = 0; 238 unsigned int index = 0; 239 for(;it!=m_vector.end();++it,++ait,index++ 239 for(;it!=m_vector.end();++it,++ait,index++) { 240 if(!a_cut || (a_cut && accept(index,*a_c 240 if(!a_cut || (a_cut && accept(index,*a_cut)) ) { 241 (*it) += (*ait); 241 (*it) += (*ait); 242 } 242 } 243 } 243 } 244 return true; 244 return true; 245 } 245 } 246 bool subtract(const array& a_array) { 246 bool subtract(const array& a_array) { 247 if(m_orders!=a_array.m_orders) return fals 247 if(m_orders!=a_array.m_orders) return false; 248 vec_it_t it = m_vector.begin(); 248 vec_it_t it = m_vector.begin(); 249 cons_vec_it_t ait = a_array.m_vector.begin 249 cons_vec_it_t ait = a_array.m_vector.begin(); 250 for(;it!=m_vector.end();++it,++ait) (*it) 250 for(;it!=m_vector.end();++it,++ait) (*it) -= (*ait); 251 return true; 251 return true; 252 } 252 } 253 bool multiply(const array& a_array) { 253 bool multiply(const array& a_array) { 254 if(m_orders!=a_array.m_orders) return fals 254 if(m_orders!=a_array.m_orders) return false; 255 vec_it_t it = m_vector.begin(); 255 vec_it_t it = m_vector.begin(); 256 cons_vec_it_t ait = a_array.m_vector.begin 256 cons_vec_it_t ait = a_array.m_vector.begin(); 257 for(;it!=m_vector.end();++it,++ait) (*it) 257 for(;it!=m_vector.end();++it,++ait) (*it) *= (*ait); 258 return true; 258 return true; 259 } 259 } 260 bool divide(const array& a_array) { 260 bool divide(const array& a_array) { 261 if(m_orders!=a_array.m_orders) return fals 261 if(m_orders!=a_array.m_orders) return false; 262 bool status = true; 262 bool status = true; 263 vec_it_t it = m_vector.begin(); 263 vec_it_t it = m_vector.begin(); 264 cons_vec_it_t ait = a_array.m_vector.begin 264 cons_vec_it_t ait = a_array.m_vector.begin(); 265 for(;it!=m_vector.end();++it,++ait) { 265 for(;it!=m_vector.end();++it,++ait) { 266 if((*ait)==zero()) { 266 if((*ait)==zero()) { 267 (*it) = zero(); //PAW convention. 267 (*it) = zero(); //PAW convention. 268 status = false; 268 status = false; 269 } else { 269 } else { 270 (*it) /= (*ait); 270 (*it) /= (*ait); 271 } 271 } 272 } 272 } 273 return status; 273 return status; 274 } 274 } 275 bool contract(const array& a_array,T& a_valu 275 bool contract(const array& a_array,T& a_value) const { 276 a_value = zero(); 276 a_value = zero(); 277 if(m_orders!=a_array.m_orders) return fals 277 if(m_orders!=a_array.m_orders) return false; 278 cons_vec_it_t it = m_vector.begin(); 278 cons_vec_it_t it = m_vector.begin(); 279 cons_vec_it_t ait = a_array.m_vector.begin 279 cons_vec_it_t ait = a_array.m_vector.begin(); 280 for(;it!=m_vector.end();++it,++ait) a_valu 280 for(;it!=m_vector.end();++it,++ait) a_value += (*it) * (*ait); 281 return true; 281 return true; 282 } 282 } 283 //Else: 283 //Else: 284 void add(const T& a_T,cut_t* a_cut = 0) { 284 void add(const T& a_T,cut_t* a_cut = 0) { 285 unsigned int index = 0; 285 unsigned int index = 0; 286 for(vec_it_t it = m_vector.begin();it!=m_v 286 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it,index++) { 287 if(!a_cut || (a_cut && accept(index,*a_c 287 if(!a_cut || (a_cut && accept(index,*a_cut)) ) { 288 (*it) += a_T; 288 (*it) += a_T; 289 } 289 } 290 } 290 } 291 } 291 } 292 void multiply(const T& a_T) { 292 void multiply(const T& a_T) { 293 for(vec_it_t it = m_vector.begin();it!=m_v 293 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it) (*it) *= a_T; 294 } 294 } 295 bool divide(const T& a_T) { 295 bool divide(const T& a_T) { 296 if(a_T==zero()) return false; 296 if(a_T==zero()) return false; 297 for(vec_it_t it = m_vector.begin();it!=m_v 297 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it) (*it) /= a_T; 298 return true; 298 return true; 299 } 299 } 300 bool invert() { 300 bool invert() { 301 bool status = true; 301 bool status = true; 302 T v_one = one(); 302 T v_one = one(); 303 for(vec_it_t it = m_vector.begin();it!=m_v 303 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it) { 304 if((*it)==zero()) { 304 if((*it)==zero()) { 305 (*it) = zero(); //PAW convention. 305 (*it) = zero(); //PAW convention. 306 status = false; 306 status = false; 307 } else { 307 } else { 308 T _value = (*it); 308 T _value = (*it); 309 (*it) = v_one/_value; 309 (*it) = v_one/_value; 310 } 310 } 311 } 311 } 312 return status; 312 return status; 313 } 313 } 314 bool offset(const uints_t& a_is,unsigned int 314 bool offset(const uints_t& a_is,unsigned int& a_offset) const { 315 size_t dim = m_orders.size(); 315 size_t dim = m_orders.size(); 316 a_offset = 0; 316 a_offset = 0; 317 if(a_is.size()!=dim) return false; 317 if(a_is.size()!=dim) return false; 318 if(dim==0) return false; 318 if(dim==0) return false; 319 for(size_t iaxis=0;iaxis<dim;iaxis++) { 319 for(size_t iaxis=0;iaxis<dim;iaxis++) { 320 unsigned int i = a_is[iaxis]; 320 unsigned int i = a_is[iaxis]; 321 if(i>=m_orders[iaxis]) { 321 if(i>=m_orders[iaxis]) { 322 a_offset = 0; 322 a_offset = 0; 323 return false; 323 return false; 324 } 324 } 325 a_offset += i * m_offsets[iaxis]; 325 a_offset += i * m_offsets[iaxis]; 326 } 326 } 327 return true; 327 return true; 328 } 328 } 329 bool indices(unsigned int a_offset,uints_t& 329 bool indices(unsigned int a_offset,uints_t& a_is) const { 330 if(a_offset>=m_vector.size()) return false 330 if(a_offset>=m_vector.size()) return false; 331 size_t dim = m_orders.size(); 331 size_t dim = m_orders.size(); 332 unsigned int off = a_offset; 332 unsigned int off = a_offset; 333 for(int iaxis=int(dim)-1;iaxis>=0;iaxis--) 333 for(int iaxis=int(dim)-1;iaxis>=0;iaxis--) { 334 a_is[iaxis] = off/m_offsets[iaxis]; 334 a_is[iaxis] = off/m_offsets[iaxis]; 335 off -= a_is[iaxis] * m_offsets[iaxis]; 335 off -= a_is[iaxis] * m_offsets[iaxis]; 336 } 336 } 337 return true; 337 return true; 338 } 338 } 339 bool accept(unsigned int a_index,const cut_t 339 bool accept(unsigned int a_index,const cut_t& a_cut) const { 340 size_t dim = m_orders.size(); 340 size_t dim = m_orders.size(); 341 if(a_cut.size()!=dim) return false; 341 if(a_cut.size()!=dim) return false; 342 if(!indices(a_index,const_cast<uints_t&>(m 342 if(!indices(a_index,const_cast<uints_t&>(m_is))) return false; 343 bool good = true; 343 bool good = true; 344 for(size_t iaxis=0;iaxis<dim;iaxis++) { 344 for(size_t iaxis=0;iaxis<dim;iaxis++) { 345 if(m_is[iaxis]<a_cut[iaxis].first) {good 345 if(m_is[iaxis]<a_cut[iaxis].first) {good = false;break;} 346 if(m_is[iaxis]>a_cut[iaxis].second) {goo 346 if(m_is[iaxis]>a_cut[iaxis].second) {good = false;break;} 347 } 347 } 348 return good; 348 return good; 349 } 349 } 350 void set_constant(const T& a_v) { 350 void set_constant(const T& a_v) { 351 for(vec_it_t it = m_vector.begin();it!=m_v 351 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it) (*it) = a_v; 352 } 352 } 353 void set_zero(){ 353 void set_zero(){ 354 set_constant(zero()); 354 set_constant(zero()); 355 } 355 } 356 public: 356 public: 357 static T zero() { 357 static T zero() { 358 //return (T)0; 358 //return (T)0; 359 //return T(0); 359 //return T(0); 360 return T(); 360 return T(); 361 } 361 } 362 static T one() { 362 static T one() { 363 //return (T)1; 363 //return (T)1; 364 return T(1); 364 return T(1); 365 } 365 } 366 static T minus_one() { 366 static T minus_one() { 367 //return (T)-1; 367 //return (T)-1; 368 return T(-1); 368 return T(-1); 369 } 369 } 370 static T two() { 370 static T two() { 371 //return (T)2; 371 //return (T)2; 372 return T(2); 372 return T(2); 373 } 373 } 374 protected: 374 protected: 375 uints_t m_orders; 375 uints_t m_orders; 376 uints_t m_offsets; 376 uints_t m_offsets; 377 std::vector<T> m_vector; 377 std::vector<T> m_vector; 378 uints_t m_is; 378 uints_t m_is; 379 }; 379 }; 380 380 381 } 381 } 382 382 383 #include <ostream> 383 #include <ostream> 384 384 385 // Helpers array<T> : 385 // Helpers array<T> : 386 namespace tools { 386 namespace tools { 387 387 388 template <class T> 388 template <class T> 389 inline bool contract( 389 inline bool contract( 390 const array<T>& aA 390 const array<T>& aA 391 ,unsigned int aIA 391 ,unsigned int aIA 392 ,const array<T>& aB 392 ,const array<T>& aB 393 ,unsigned int aIB 393 ,unsigned int aIB 394 ,array<T>& aR 394 ,array<T>& aR 395 ){ 395 ){ 396 // Contract a tensor (aA) with a vector (aB) 396 // Contract a tensor (aA) with a vector (aB) 397 // on indice aIA of aA and aIB of aB. 397 // on indice aIA of aA and aIB of aB. 398 // aA.dimension must be > 1. 398 // aA.dimension must be > 1. 399 // aB.dimension must be > 1. 399 // aB.dimension must be > 1. 400 if(&aR==&aA) return false; 400 if(&aR==&aA) return false; 401 if(&aR==&aB) return false; 401 if(&aR==&aB) return false; 402 402 403 if(aA.dimension()==0) return false; 403 if(aA.dimension()==0) return false; 404 if(aB.dimension()==0) return false; 404 if(aB.dimension()==0) return false; 405 if(aIA>=aA.dimension()) return false; 405 if(aIA>=aA.dimension()) return false; 406 if(aIB>=aB.dimension()) return false; 406 if(aIB>=aB.dimension()) return false; 407 if(aA.orders()[aIA]!=aB.orders()[aIB]) retur 407 if(aA.orders()[aIA]!=aB.orders()[aIB]) return false; 408 408 409 unsigned int rdima = aA.dimension()-1; 409 unsigned int rdima = aA.dimension()-1; 410 unsigned int rdimb = aB.dimension()-1; 410 unsigned int rdimb = aB.dimension()-1; 411 411 412 if((rdima+rdimb)==0) { 412 if((rdima+rdimb)==0) { 413 std::vector<unsigned int> rorders(1); 413 std::vector<unsigned int> rorders(1); 414 rorders[0] = 1; 414 rorders[0] = 1; 415 if(!aR.configure(rorders)) return false; 415 if(!aR.configure(rorders)) return false; 416 const std::vector<T>& vA = aA.vector(); 416 const std::vector<T>& vA = aA.vector(); 417 const std::vector<T>& vB = aB.vector(); 417 const std::vector<T>& vB = aB.vector(); 418 T value = array<T>::zero(); 418 T value = array<T>::zero(); 419 for(unsigned int index=0;index<vA.size();i 419 for(unsigned int index=0;index<vA.size();index++) { 420 value += vA[index]*vB[index]; 420 value += vA[index]*vB[index]; 421 } 421 } 422 aR.vector()[0] = value; 422 aR.vector()[0] = value; 423 return true; 423 return true; 424 } 424 } 425 425 426 std::vector<unsigned int> rorders(rdima+rdim 426 std::vector<unsigned int> rorders(rdima+rdimb); 427 unsigned int index; 427 unsigned int index; 428 for(index=0;index<aIA;index++) 428 for(index=0;index<aIA;index++) 429 rorders[index] = aA.orders()[index]; 429 rorders[index] = aA.orders()[index]; 430 for(index=aIA;index<rdima;index++) 430 for(index=aIA;index<rdima;index++) 431 rorders[index] = aA.orders()[index+1]; 431 rorders[index] = aA.orders()[index+1]; 432 432 433 for(index=0;index<aIB;index++) 433 for(index=0;index<aIB;index++) 434 rorders[rdima+index] = aB.orders()[index]; 434 rorders[rdima+index] = aB.orders()[index]; 435 for(index=aIB;index<rdimb;index++) 435 for(index=aIB;index<rdimb;index++) 436 rorders[rdima+index] = aB.orders()[index+1 436 rorders[rdima+index] = aB.orders()[index+1]; 437 437 438 if(!aR.configure(rorders)) return false; 438 if(!aR.configure(rorders)) return false; 439 439 440 std::vector<unsigned int> ais(aA.dimension() 440 std::vector<unsigned int> ais(aA.dimension()); 441 std::vector<unsigned int> bis(aB.dimension() 441 std::vector<unsigned int> bis(aB.dimension()); 442 std::vector<unsigned int> ris(aR.dimension() 442 std::vector<unsigned int> ris(aR.dimension()); 443 443 444 //FIXME : optimize the below. 444 //FIXME : optimize the below. 445 unsigned int order = aA.orders()[aIA]; 445 unsigned int order = aA.orders()[aIA]; 446 unsigned int rsize = aR.size(); 446 unsigned int rsize = aR.size(); 447 447 448 std::vector<T>& rvec = aR.vector(); 448 std::vector<T>& rvec = aR.vector(); 449 449 450 for(unsigned int roffset=0;roffset<rsize;rof 450 for(unsigned int roffset=0;roffset<rsize;roffset++) { 451 if(!aR.indices(roffset,ris)) return false; 451 if(!aR.indices(roffset,ris)) return false; 452 452 453 for(index=0;index<aIA;index++) ais[index] 453 for(index=0;index<aIA;index++) ais[index] = ris[index]; 454 for(index=aIA;index<rdima;index++) ais[ind 454 for(index=aIA;index<rdima;index++) ais[index+1] = ris[index]; 455 455 456 for(index=0;index<aIB;index++) bis[index] 456 for(index=0;index<aIB;index++) bis[index] = ris[rdima+index]; 457 for(index=aIB;index<rdimb;index++) bis[ind 457 for(index=aIB;index<rdimb;index++) bis[index+1] = ris[rdima+index]; 458 458 459 T value = 0; 459 T value = 0; 460 for(index=0;index<order;index++) { 460 for(index=0;index<order;index++) { 461 ais[aIA] = index; 461 ais[aIA] = index; 462 T av; 462 T av; 463 if(!aA.value(ais,av)) return false; 463 if(!aA.value(ais,av)) return false; 464 464 465 bis[aIB] = index; 465 bis[aIB] = index; 466 T bv; 466 T bv; 467 if(!aB.value(bis,bv)) return false; 467 if(!aB.value(bis,bv)) return false; 468 468 469 value += av * bv; 469 value += av * bv; 470 } 470 } 471 471 472 rvec[roffset] = value; 472 rvec[roffset] = value; 473 } 473 } 474 474 475 return true; 475 return true; 476 } 476 } 477 477 478 template <class T> 478 template <class T> 479 inline bool swap( 479 inline bool swap( 480 const array<T>& aV 480 const array<T>& aV 481 ,unsigned int aI1 481 ,unsigned int aI1 482 ,unsigned int aI2 482 ,unsigned int aI2 483 ,array<T>& aR 483 ,array<T>& aR 484 ){ 484 ){ 485 if(&aR==&aV) return false; 485 if(&aR==&aV) return false; 486 486 487 unsigned int dim = aV.dimension(); 487 unsigned int dim = aV.dimension(); 488 if(aI1>=dim) return false; 488 if(aI1>=dim) return false; 489 if(aI2>=dim) return false; 489 if(aI2>=dim) return false; 490 490 491 if(aI1==aI2) { 491 if(aI1==aI2) { 492 aR.copy(aV); 492 aR.copy(aV); 493 return true; 493 return true; 494 } 494 } 495 495 496 if(!aR.configure(aV.orders())) return false; 496 if(!aR.configure(aV.orders())) return false; 497 497 498 std::vector<unsigned int> vis(aV.dimension() 498 std::vector<unsigned int> vis(aV.dimension()); 499 std::vector<unsigned int> ris(aV.dimension() 499 std::vector<unsigned int> ris(aV.dimension()); 500 500 501 const std::vector<T>& vvec = aV.vector(); 501 const std::vector<T>& vvec = aV.vector(); 502 502 503 unsigned int size = aV.size(); 503 unsigned int size = aV.size(); 504 for(unsigned int offset=0;offset<size;offset 504 for(unsigned int offset=0;offset<size;offset++) { 505 T value = vvec[offset]; 505 T value = vvec[offset]; 506 506 507 if(!aV.indices(offset,vis)) return false; 507 if(!aV.indices(offset,vis)) return false; 508 unsigned int i = vis[aI1]; 508 unsigned int i = vis[aI1]; 509 vis[aI1] = vis[aI2]; 509 vis[aI1] = vis[aI2]; 510 vis[aI2] = i; 510 vis[aI2] = i; 511 if(!aR.set_value(vis,value)) return false; 511 if(!aR.set_value(vis,value)) return false; 512 } 512 } 513 513 514 return true; 514 return true; 515 } 515 } 516 516 517 //NOTE : print is a Python keyword. 517 //NOTE : print is a Python keyword. 518 template <class T> 518 template <class T> 519 inline void dump(std::ostream& a_out,const too 519 inline void dump(std::ostream& a_out,const tools::array<T>& a_array,const std::string& a_title){ 520 if(a_title.size()) a_out << a_title << std:: 520 if(a_title.size()) a_out << a_title << std::endl; 521 const std::vector<T>& vec = a_array.vector() 521 const std::vector<T>& vec = a_array.vector(); 522 typedef typename std::vector<T>::const_itera 522 typedef typename std::vector<T>::const_iterator cons_vec_it_t; 523 for(cons_vec_it_t it = vec.begin();it!=vec.e 523 for(cons_vec_it_t it = vec.begin();it!=vec.end();++it) { 524 a_out << (*it) << std::endl; 524 a_out << (*it) << std::endl; 525 } 525 } 526 } 526 } 527 527 528 template <class T> 528 template <class T> 529 inline void diff(std::ostream& a_out,const arr 529 inline void diff(std::ostream& a_out,const array<T>& aA,const array<T>& aB,T a_epsilon){ 530 if(aA.orders()!=aB.orders()) { 530 if(aA.orders()!=aB.orders()) { 531 a_out << "tools::arrays::diff : not same o 531 a_out << "tools::arrays::diff : not same orders" << std::endl; 532 return; 532 return; 533 } 533 } 534 bool header_done = false; 534 bool header_done = false; 535 unsigned int dim = aA.dimension(); 535 unsigned int dim = aA.dimension(); 536 std::vector<unsigned int> is(dim); 536 std::vector<unsigned int> is(dim); 537 unsigned int vsize = aA.vector().size(); 537 unsigned int vsize = aA.vector().size(); 538 for(unsigned int index=0;index<vsize;index++ 538 for(unsigned int index=0;index<vsize;index++) { 539 T diff = aA.vector()[index]-aB.vector()[in 539 T diff = aA.vector()[index]-aB.vector()[index]; 540 if(diff<0) diff *= -1; 540 if(diff<0) diff *= -1; 541 if(diff>=a_epsilon) { 541 if(diff>=a_epsilon) { 542 aA.indices(index,is); 542 aA.indices(index,is); 543 if(!header_done) { 543 if(!header_done) { 544 a_out << "tools::arrays::diff :" << st 544 a_out << "tools::arrays::diff :" << std::endl; 545 header_done = true; 545 header_done = true; 546 } 546 } 547 for(unsigned int i=0;i<dim;i++) a_out << 547 for(unsigned int i=0;i<dim;i++) a_out << is[i] << " "; 548 a_out << aA.vector()[index] << " " << aB 548 a_out << aA.vector()[index] << " " << aB.vector()[index] << std::endl; 549 } 549 } 550 } 550 } 551 } 551 } 552 552 553 ////////////////////////////////////////////// 553 ////////////////////////////////////////////////////////////////////////////// 554 /// common array ///////////////////////////// 554 /// common array ///////////////////////////////////////////////////////////// 555 ////////////////////////////////////////////// 555 ////////////////////////////////////////////////////////////////////////////// 556 556 557 template <class T> 557 template <class T> 558 class kronecker : public array<T> { 558 class kronecker : public array<T> { 559 public: 559 public: 560 kronecker(unsigned int a_order):array<T>(a_o 560 kronecker(unsigned int a_order):array<T>(a_order,a_order){ 561 // epsilon(i1,i2,....in) with n = a_order 561 // epsilon(i1,i2,....in) with n = a_order and i in [0,n[. 562 // Then an Array of a_order * a_order. 562 // Then an Array of a_order * a_order. 563 unsigned int index = 0; 563 unsigned int index = 0; 564 typedef typename std::vector<unsigned int> 564 typedef typename std::vector<unsigned int> _uints_t; 565 _uints_t is(a_order); 565 _uints_t is(a_order); 566 std::vector<T>& vec = array<T>::vector(); 566 std::vector<T>& vec = array<T>::vector(); 567 typedef typename std::vector<T>::iterator 567 typedef typename std::vector<T>::iterator _vec_it_t; 568 _vec_it_t it = vec.begin(); 568 _vec_it_t it = vec.begin(); 569 for(;it!=vec.end();++it,index++) { 569 for(;it!=vec.end();++it,index++) { 570 if(!array<T>::indices(index,is)) return; 570 if(!array<T>::indices(index,is)) return; //FIXME throw. 571 bool good = true; 571 bool good = true; 572 {for(unsigned int iaxis=0;iaxis<a_order;i 572 {for(unsigned int iaxis=0;iaxis<a_order;iaxis++) { 573 unsigned int ival = is[iaxis]; 573 unsigned int ival = is[iaxis]; 574 for(unsigned int iaxis2=iaxis+1;iaxis2 574 for(unsigned int iaxis2=iaxis+1;iaxis2<a_order;iaxis2++) { 575 if(is[iaxis2]==ival) { 575 if(is[iaxis2]==ival) { 576 good = false; 576 good = false; 577 break; 577 break; 578 } 578 } 579 } 579 } 580 if(!good) break; 580 if(!good) break; 581 }} 581 }} 582 if(!good) continue; 582 if(!good) continue; 583 // All indicies are different. 583 // All indicies are different. 584 unsigned int n = 0; 584 unsigned int n = 0; 585 for(unsigned int iaxis=0;iaxis<a_order;) 585 for(unsigned int iaxis=0;iaxis<a_order;) { 586 unsigned int ival = is[iaxis]; 586 unsigned int ival = is[iaxis]; 587 if(ival!=iaxis) { 587 if(ival!=iaxis) { 588 // Swap and add one permutation : 588 // Swap and add one permutation : 589 unsigned int old = is[ival]; 589 unsigned int old = is[ival]; 590 is[ival] = ival; 590 is[ival] = ival; 591 is[iaxis] = old; 591 is[iaxis] = old; 592 n +=1; 592 n +=1; 593 } else { 593 } else { 594 iaxis++; 594 iaxis++; 595 } 595 } 596 } 596 } 597 {unsigned int n_2 = n/2; 597 {unsigned int n_2 = n/2; 598 if(2*n_2==n) (*it) = array<T>::one(); 598 if(2*n_2==n) (*it) = array<T>::one(); 599 else (*it) = array<T>::minus_one();} 599 else (*it) = array<T>::minus_one();} 600 } 600 } 601 } 601 } 602 }; 602 }; 603 603 604 template <class T> 604 template <class T> 605 inline array<T> operator+(const array<T>& a1,c 605 inline array<T> operator+(const array<T>& a1,const array<T>& a2) { 606 array<T> res(a1); 606 array<T> res(a1); 607 if(!res.add(a2)) {} 607 if(!res.add(a2)) {} 608 return res; 608 return res; 609 } 609 } 610 template <class T> 610 template <class T> 611 inline array<T> operator-(const array<T>& a1,c 611 inline array<T> operator-(const array<T>& a1,const array<T>& a2) { 612 array<T> res(a1); 612 array<T> res(a1); 613 if(!res.subtract(a2)) {} 613 if(!res.subtract(a2)) {} 614 return res; 614 return res; 615 } 615 } 616 template <class T> 616 template <class T> 617 inline array<T> operator*(const T& a_fac,const 617 inline array<T> operator*(const T& a_fac,const array<T>& a_m) { 618 array<T> res(a_m); 618 array<T> res(a_m); 619 res *= a_fac; 619 res *= a_fac; 620 return res; 620 return res; 621 } 621 } 622 622 623 } 623 } 624 624 625 #endif 625 #endif