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