| 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