Warning, /include/Geant4/tools/array is written in an unsupported language. File is not indexed.
0001 // Copyright (C) 2010, Guy Barrand. All rights reserved.
0002 // See the file tools.license for terms.
0003
0004 #ifndef tools_array
0005 #define tools_array
0006
0007 #ifdef TOOLS_MEM
0008 #include "mem"
0009 #endif
0010
0011 #include <vector>
0012 #include <string>
0013
0014 namespace tools {
0015
0016 // array handles an hyperparallelepiped of cells of class T.
0017
0018 template <class T>
0019 class array {
0020 public:
0021 static const std::string& s_class() {
0022 static const std::string s_v("tools::array");
0023 return s_v;
0024 }
0025 public:
0026 typedef typename std::vector< std::pair<unsigned int,unsigned int> > cut_t;
0027 typedef typename std::vector<unsigned int> uints_t;
0028 typedef typename std::vector<T>::iterator vec_it_t;
0029 typedef typename std::vector<T>::const_iterator cons_vec_it_t;
0030 public:
0031 array() {
0032 #ifdef TOOLS_MEM
0033 mem::increment(s_class().c_str());
0034 #endif
0035 }
0036 array(const uints_t& a_orders) {
0037 #ifdef TOOLS_MEM
0038 mem::increment(s_class().c_str());
0039 #endif
0040 configure(a_orders);
0041 }
0042 array(unsigned int a_dimension,unsigned int a_order) {
0043 #ifdef TOOLS_MEM
0044 mem::increment(s_class().c_str());
0045 #endif
0046 // A hypercube of dimension "a_dimension" and size "a_order".
0047 uints_t _orders(a_dimension);
0048 for(unsigned int index=0;index<a_dimension;index++)
0049 _orders[index] = a_order;
0050 configure(_orders);
0051 }
0052 virtual ~array() {
0053 #ifdef TOOLS_MEM
0054 mem::decrement(s_class().c_str());
0055 #endif
0056 }
0057 public:
0058 array(const array& a_from)
0059 :m_orders(a_from.m_orders)
0060 ,m_offsets(a_from.m_offsets)
0061 ,m_vector(a_from.m_vector)
0062 ,m_is(a_from.m_is){
0063 #ifdef TOOLS_MEM
0064 mem::increment(s_class().c_str());
0065 #endif
0066 }
0067 array& operator=(const array& a_from) {
0068 m_orders = a_from.m_orders;
0069 m_offsets = a_from.m_offsets;
0070 m_vector = a_from.m_vector;
0071 m_is = a_from.m_is;
0072 return *this;
0073 }
0074 public: //operators:
0075 array& operator*=(const T& a_T) {multiply(a_T);return *this;}
0076 bool operator==(const array& a_array) const {
0077 return equal(a_array);
0078 }
0079 bool operator!=(const array& a_array) const {
0080 return !operator==(a_array);
0081 }
0082 array operator*(const T& a_T) const {
0083 array a(*this);
0084 a.multiply(a_T);
0085 return a;
0086 }
0087 // the below would need exception to do it properly.
0088 //array operator+(const array& a_array) const {
0089 // array a(*this);
0090 // if(!a.add(a_array)) {} //throw
0091 // return a;
0092 //}
0093 //array& operator+=(const array& a_array) {
0094 // if(!add(a_array)) {} //throw
0095 // return *this;
0096 //}
0097 public:
0098 void copy(const array& a_from) {
0099 m_orders = a_from.m_orders;
0100 m_offsets = a_from.m_offsets;
0101 m_vector = a_from.m_vector;
0102 m_is = a_from.m_is;
0103 }
0104 void clear() {
0105 m_orders.clear();
0106 m_offsets.clear();
0107 m_vector.clear();
0108 m_is.clear();
0109 }
0110 bool configure(const uints_t& a_orders) {
0111 m_orders = a_orders;
0112 size_t dim = m_orders.size();
0113 if(dim==0) {
0114 clear();
0115 return false;
0116 }
0117 unsigned int _size = 1;
0118 for(size_t index=0;index<dim;index++) {
0119 //if(m_orders[index]<0) {
0120 //clear();
0121 //return false;
0122 //}
0123 _size *= m_orders[index];
0124 }
0125 m_vector.resize(_size);
0126 m_vector.assign(_size,zero());
0127 m_offsets.resize(dim,0);
0128 m_offsets[0] = 1;
0129 for(size_t iaxis=1;iaxis<dim;iaxis++)
0130 m_offsets[iaxis] = m_offsets[iaxis-1] * m_orders[iaxis-1];
0131 m_is.resize(dim);
0132 m_is.assign(dim,0);
0133 return true;
0134 }
0135 size_t dimension() const { return m_orders.size();}
0136 const uints_t& orders() const { return m_orders;}
0137 size_t size() const {return m_vector.size();}
0138 bool set_value(const uints_t& a_is,const T& a_value) {
0139 unsigned int off = 0;
0140 if(!offset(a_is,off)) return false;
0141 m_vector[off] = a_value;
0142 return true;
0143 }
0144 bool value(const uints_t& a_is,T& a_value) const {
0145 unsigned int off = 0;
0146 if(!offset(a_is,off)) {
0147 a_value = 0;
0148 return false;
0149 }
0150 a_value = m_vector[off];
0151 return true;
0152 }
0153 T value_no_check(const uints_t& a_is) const { //TOUCHY
0154 unsigned int off = 0;
0155 if(!offset(a_is,off)) {}
0156 return m_vector[off];
0157 }
0158 void reset() {
0159 for(vec_it_t it=m_vector.begin();it!=m_vector.end();++it) *it = 0;
0160 }
0161 const std::vector<T>& vector() const { return m_vector;}
0162 std::vector<T>& vector() { return m_vector;}
0163 bool fill(const std::vector<T>& a_values,cut_t* a_cut = 0) {
0164 size_t dsize = a_values.size();
0165 size_t di = 0;
0166 unsigned int index = 0;
0167 for(vec_it_t it=m_vector.begin();it!=m_vector.end();++it,index++) {
0168 if(!a_cut || (a_cut && accept(index,*a_cut)) ) {
0169 if(di>=dsize) return false; //a_values exhausted too early
0170 *it = a_values[di];
0171 di++;
0172 }
0173 }
0174 return true;
0175 }
0176 bool fill(unsigned int a_sz,const T* a_data,cut_t* a_cut = 0) {
0177 unsigned int di = 0;
0178 unsigned int index = 0;
0179 for(vec_it_t it=m_vector.begin();it!=m_vector.end();++it,index++) {
0180 if(!a_cut || (a_cut && accept(index,*a_cut)) ) {
0181 if(di>=a_sz) return false; //a_values exhausted too early
0182 *it = a_data[di];
0183 di++;
0184 }
0185 }
0186 return true;
0187 }
0188 // Related to other array :
0189 bool equal(const array& a_array) const {
0190 if(m_orders!=a_array.m_orders) return false;
0191 cons_vec_it_t it = m_vector.begin();
0192 cons_vec_it_t ait = a_array.m_vector.begin();
0193 for(;it!=m_vector.end();++it,++ait) {
0194 if((*it)!=(*ait)) return false;
0195 }
0196 return true;
0197 }
0198 bool equal(const array& a_array,T aEpsilon) const {
0199 if(m_orders!=a_array.m_orders) return false;
0200 cons_vec_it_t it = m_vector.begin();
0201 cons_vec_it_t ait = a_array.m_vector.begin();
0202 for(;it!=m_vector.end();++it,++ait) {
0203 T diff = (*it) - (*ait);
0204 if(diff<0) diff *= -1;
0205 if(diff>=aEpsilon) return false;
0206 }
0207 return true;
0208 }
0209 bool is_proportional(const array& a_array,T& a_factor) const {
0210 // If true, then : a_array = a_factor * this.
0211 a_factor = zero();
0212 if(m_orders!=a_array.m_orders) return false;
0213 bool first = true;
0214 cons_vec_it_t it = m_vector.begin();
0215 cons_vec_it_t ait = a_array.m_vector.begin();
0216 for(;it!=m_vector.end();++it,++ait) {
0217 if( ((*it)==zero()) && ((*ait)==zero())) {
0218 continue;
0219 } else if( ((*it)!=zero()) && ((*ait)==zero())) {
0220 return false;
0221 } else if( ((*it)==zero()) && ((*ait)!=zero())) {
0222 return false;
0223 } else {
0224 if(first) {
0225 a_factor = (*ait)/(*it);
0226 first = false;
0227 } else {
0228 if((*ait)!=(*it)*a_factor) return false;
0229 }
0230 }
0231 }
0232 return true;
0233 }
0234 bool add(const array& a_array,cut_t* a_cut = 0) {
0235 if(m_orders!=a_array.m_orders) return false;
0236 vec_it_t it = m_vector.begin();
0237 cons_vec_it_t ait = a_array.m_vector.begin();
0238 unsigned int index = 0;
0239 for(;it!=m_vector.end();++it,++ait,index++) {
0240 if(!a_cut || (a_cut && accept(index,*a_cut)) ) {
0241 (*it) += (*ait);
0242 }
0243 }
0244 return true;
0245 }
0246 bool subtract(const array& a_array) {
0247 if(m_orders!=a_array.m_orders) return false;
0248 vec_it_t it = m_vector.begin();
0249 cons_vec_it_t ait = a_array.m_vector.begin();
0250 for(;it!=m_vector.end();++it,++ait) (*it) -= (*ait);
0251 return true;
0252 }
0253 bool multiply(const array& a_array) {
0254 if(m_orders!=a_array.m_orders) return false;
0255 vec_it_t it = m_vector.begin();
0256 cons_vec_it_t ait = a_array.m_vector.begin();
0257 for(;it!=m_vector.end();++it,++ait) (*it) *= (*ait);
0258 return true;
0259 }
0260 bool divide(const array& a_array) {
0261 if(m_orders!=a_array.m_orders) return false;
0262 bool status = true;
0263 vec_it_t it = m_vector.begin();
0264 cons_vec_it_t ait = a_array.m_vector.begin();
0265 for(;it!=m_vector.end();++it,++ait) {
0266 if((*ait)==zero()) {
0267 (*it) = zero(); //PAW convention.
0268 status = false;
0269 } else {
0270 (*it) /= (*ait);
0271 }
0272 }
0273 return status;
0274 }
0275 bool contract(const array& a_array,T& a_value) const {
0276 a_value = zero();
0277 if(m_orders!=a_array.m_orders) return false;
0278 cons_vec_it_t it = m_vector.begin();
0279 cons_vec_it_t ait = a_array.m_vector.begin();
0280 for(;it!=m_vector.end();++it,++ait) a_value += (*it) * (*ait);
0281 return true;
0282 }
0283 //Else:
0284 void add(const T& a_T,cut_t* a_cut = 0) {
0285 unsigned int index = 0;
0286 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it,index++) {
0287 if(!a_cut || (a_cut && accept(index,*a_cut)) ) {
0288 (*it) += a_T;
0289 }
0290 }
0291 }
0292 void multiply(const T& a_T) {
0293 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it) (*it) *= a_T;
0294 }
0295 bool divide(const T& a_T) {
0296 if(a_T==zero()) return false;
0297 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it) (*it) /= a_T;
0298 return true;
0299 }
0300 bool invert() {
0301 bool status = true;
0302 T v_one = one();
0303 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it) {
0304 if((*it)==zero()) {
0305 (*it) = zero(); //PAW convention.
0306 status = false;
0307 } else {
0308 T _value = (*it);
0309 (*it) = v_one/_value;
0310 }
0311 }
0312 return status;
0313 }
0314 bool offset(const uints_t& a_is,unsigned int& a_offset) const {
0315 size_t dim = m_orders.size();
0316 a_offset = 0;
0317 if(a_is.size()!=dim) return false;
0318 if(dim==0) return false;
0319 for(size_t iaxis=0;iaxis<dim;iaxis++) {
0320 unsigned int i = a_is[iaxis];
0321 if(i>=m_orders[iaxis]) {
0322 a_offset = 0;
0323 return false;
0324 }
0325 a_offset += i * m_offsets[iaxis];
0326 }
0327 return true;
0328 }
0329 bool indices(unsigned int a_offset,uints_t& a_is) const {
0330 if(a_offset>=m_vector.size()) return false;
0331 size_t dim = m_orders.size();
0332 unsigned int off = a_offset;
0333 for(int iaxis=int(dim)-1;iaxis>=0;iaxis--) {
0334 a_is[iaxis] = off/m_offsets[iaxis];
0335 off -= a_is[iaxis] * m_offsets[iaxis];
0336 }
0337 return true;
0338 }
0339 bool accept(unsigned int a_index,const cut_t& a_cut) const {
0340 size_t dim = m_orders.size();
0341 if(a_cut.size()!=dim) return false;
0342 if(!indices(a_index,const_cast<uints_t&>(m_is))) return false;
0343 bool good = true;
0344 for(size_t iaxis=0;iaxis<dim;iaxis++) {
0345 if(m_is[iaxis]<a_cut[iaxis].first) {good = false;break;}
0346 if(m_is[iaxis]>a_cut[iaxis].second) {good = false;break;}
0347 }
0348 return good;
0349 }
0350 void set_constant(const T& a_v) {
0351 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it) (*it) = a_v;
0352 }
0353 void set_zero(){
0354 set_constant(zero());
0355 }
0356 public:
0357 static T zero() {
0358 //return (T)0;
0359 //return T(0);
0360 return T();
0361 }
0362 static T one() {
0363 //return (T)1;
0364 return T(1);
0365 }
0366 static T minus_one() {
0367 //return (T)-1;
0368 return T(-1);
0369 }
0370 static T two() {
0371 //return (T)2;
0372 return T(2);
0373 }
0374 protected:
0375 uints_t m_orders;
0376 uints_t m_offsets;
0377 std::vector<T> m_vector;
0378 uints_t m_is;
0379 };
0380
0381 }
0382
0383 #include <ostream>
0384
0385 // Helpers array<T> :
0386 namespace tools {
0387
0388 template <class T>
0389 inline bool contract(
0390 const array<T>& aA
0391 ,unsigned int aIA
0392 ,const array<T>& aB
0393 ,unsigned int aIB
0394 ,array<T>& aR
0395 ){
0396 // Contract a tensor (aA) with a vector (aB)
0397 // on indice aIA of aA and aIB of aB.
0398 // aA.dimension must be > 1.
0399 // aB.dimension must be > 1.
0400 if(&aR==&aA) return false;
0401 if(&aR==&aB) return false;
0402
0403 if(aA.dimension()==0) return false;
0404 if(aB.dimension()==0) return false;
0405 if(aIA>=aA.dimension()) return false;
0406 if(aIB>=aB.dimension()) return false;
0407 if(aA.orders()[aIA]!=aB.orders()[aIB]) return false;
0408
0409 unsigned int rdima = aA.dimension()-1;
0410 unsigned int rdimb = aB.dimension()-1;
0411
0412 if((rdima+rdimb)==0) {
0413 std::vector<unsigned int> rorders(1);
0414 rorders[0] = 1;
0415 if(!aR.configure(rorders)) return false;
0416 const std::vector<T>& vA = aA.vector();
0417 const std::vector<T>& vB = aB.vector();
0418 T value = array<T>::zero();
0419 for(unsigned int index=0;index<vA.size();index++) {
0420 value += vA[index]*vB[index];
0421 }
0422 aR.vector()[0] = value;
0423 return true;
0424 }
0425
0426 std::vector<unsigned int> rorders(rdima+rdimb);
0427 unsigned int index;
0428 for(index=0;index<aIA;index++)
0429 rorders[index] = aA.orders()[index];
0430 for(index=aIA;index<rdima;index++)
0431 rorders[index] = aA.orders()[index+1];
0432
0433 for(index=0;index<aIB;index++)
0434 rorders[rdima+index] = aB.orders()[index];
0435 for(index=aIB;index<rdimb;index++)
0436 rorders[rdima+index] = aB.orders()[index+1];
0437
0438 if(!aR.configure(rorders)) return false;
0439
0440 std::vector<unsigned int> ais(aA.dimension());
0441 std::vector<unsigned int> bis(aB.dimension());
0442 std::vector<unsigned int> ris(aR.dimension());
0443
0444 //FIXME : optimize the below.
0445 unsigned int order = aA.orders()[aIA];
0446 unsigned int rsize = aR.size();
0447
0448 std::vector<T>& rvec = aR.vector();
0449
0450 for(unsigned int roffset=0;roffset<rsize;roffset++) {
0451 if(!aR.indices(roffset,ris)) return false;
0452
0453 for(index=0;index<aIA;index++) ais[index] = ris[index];
0454 for(index=aIA;index<rdima;index++) ais[index+1] = ris[index];
0455
0456 for(index=0;index<aIB;index++) bis[index] = ris[rdima+index];
0457 for(index=aIB;index<rdimb;index++) bis[index+1] = ris[rdima+index];
0458
0459 T value = 0;
0460 for(index=0;index<order;index++) {
0461 ais[aIA] = index;
0462 T av;
0463 if(!aA.value(ais,av)) return false;
0464
0465 bis[aIB] = index;
0466 T bv;
0467 if(!aB.value(bis,bv)) return false;
0468
0469 value += av * bv;
0470 }
0471
0472 rvec[roffset] = value;
0473 }
0474
0475 return true;
0476 }
0477
0478 template <class T>
0479 inline bool swap(
0480 const array<T>& aV
0481 ,unsigned int aI1
0482 ,unsigned int aI2
0483 ,array<T>& aR
0484 ){
0485 if(&aR==&aV) return false;
0486
0487 unsigned int dim = aV.dimension();
0488 if(aI1>=dim) return false;
0489 if(aI2>=dim) return false;
0490
0491 if(aI1==aI2) {
0492 aR.copy(aV);
0493 return true;
0494 }
0495
0496 if(!aR.configure(aV.orders())) return false;
0497
0498 std::vector<unsigned int> vis(aV.dimension());
0499 std::vector<unsigned int> ris(aV.dimension());
0500
0501 const std::vector<T>& vvec = aV.vector();
0502
0503 unsigned int size = aV.size();
0504 for(unsigned int offset=0;offset<size;offset++) {
0505 T value = vvec[offset];
0506
0507 if(!aV.indices(offset,vis)) return false;
0508 unsigned int i = vis[aI1];
0509 vis[aI1] = vis[aI2];
0510 vis[aI2] = i;
0511 if(!aR.set_value(vis,value)) return false;
0512 }
0513
0514 return true;
0515 }
0516
0517 //NOTE : print is a Python keyword.
0518 template <class T>
0519 inline void dump(std::ostream& a_out,const tools::array<T>& a_array,const std::string& a_title){
0520 if(a_title.size()) a_out << a_title << std::endl;
0521 const std::vector<T>& vec = a_array.vector();
0522 typedef typename std::vector<T>::const_iterator cons_vec_it_t;
0523 for(cons_vec_it_t it = vec.begin();it!=vec.end();++it) {
0524 a_out << (*it) << std::endl;
0525 }
0526 }
0527
0528 template <class T>
0529 inline void diff(std::ostream& a_out,const array<T>& aA,const array<T>& aB,T a_epsilon){
0530 if(aA.orders()!=aB.orders()) {
0531 a_out << "tools::arrays::diff : not same orders" << std::endl;
0532 return;
0533 }
0534 bool header_done = false;
0535 unsigned int dim = aA.dimension();
0536 std::vector<unsigned int> is(dim);
0537 unsigned int vsize = aA.vector().size();
0538 for(unsigned int index=0;index<vsize;index++) {
0539 T diff = aA.vector()[index]-aB.vector()[index];
0540 if(diff<0) diff *= -1;
0541 if(diff>=a_epsilon) {
0542 aA.indices(index,is);
0543 if(!header_done) {
0544 a_out << "tools::arrays::diff :" << std::endl;
0545 header_done = true;
0546 }
0547 for(unsigned int i=0;i<dim;i++) a_out << is[i] << " ";
0548 a_out << aA.vector()[index] << " " << aB.vector()[index] << std::endl;
0549 }
0550 }
0551 }
0552
0553 //////////////////////////////////////////////////////////////////////////////
0554 /// common array /////////////////////////////////////////////////////////////
0555 //////////////////////////////////////////////////////////////////////////////
0556
0557 template <class T>
0558 class kronecker : public array<T> {
0559 public:
0560 kronecker(unsigned int a_order):array<T>(a_order,a_order){
0561 // epsilon(i1,i2,....in) with n = a_order and i in [0,n[.
0562 // Then an Array of a_order * a_order.
0563 unsigned int index = 0;
0564 typedef typename std::vector<unsigned int> _uints_t;
0565 _uints_t is(a_order);
0566 std::vector<T>& vec = array<T>::vector();
0567 typedef typename std::vector<T>::iterator _vec_it_t;
0568 _vec_it_t it = vec.begin();
0569 for(;it!=vec.end();++it,index++) {
0570 if(!array<T>::indices(index,is)) return; //FIXME throw.
0571 bool good = true;
0572 {for(unsigned int iaxis=0;iaxis<a_order;iaxis++) {
0573 unsigned int ival = is[iaxis];
0574 for(unsigned int iaxis2=iaxis+1;iaxis2<a_order;iaxis2++) {
0575 if(is[iaxis2]==ival) {
0576 good = false;
0577 break;
0578 }
0579 }
0580 if(!good) break;
0581 }}
0582 if(!good) continue;
0583 // All indicies are different.
0584 unsigned int n = 0;
0585 for(unsigned int iaxis=0;iaxis<a_order;) {
0586 unsigned int ival = is[iaxis];
0587 if(ival!=iaxis) {
0588 // Swap and add one permutation :
0589 unsigned int old = is[ival];
0590 is[ival] = ival;
0591 is[iaxis] = old;
0592 n +=1;
0593 } else {
0594 iaxis++;
0595 }
0596 }
0597 {unsigned int n_2 = n/2;
0598 if(2*n_2==n) (*it) = array<T>::one();
0599 else (*it) = array<T>::minus_one();}
0600 }
0601 }
0602 };
0603
0604 template <class T>
0605 inline array<T> operator+(const array<T>& a1,const array<T>& a2) {
0606 array<T> res(a1);
0607 if(!res.add(a2)) {}
0608 return res;
0609 }
0610 template <class T>
0611 inline array<T> operator-(const array<T>& a1,const array<T>& a2) {
0612 array<T> res(a1);
0613 if(!res.subtract(a2)) {}
0614 return res;
0615 }
0616 template <class T>
0617 inline array<T> operator*(const T& a_fac,const array<T>& a_m) {
0618 array<T> res(a_m);
0619 res *= a_fac;
0620 return res;
0621 }
0622
0623 }
0624
0625 #endif
