inlib
1.2.0
|
00001 // Copyright (C) 2010, Guy Barrand. All rights reserved. 00002 // See the file inlib.license for terms. 00003 00004 #ifndef inlib_array 00005 #define inlib_array 00006 00007 #ifdef INLIB_MEM 00008 #include "mem" 00009 #endif 00010 00011 #include <vector> 00012 00013 namespace inlib { 00014 00015 // array handles an hyperparallelepiped of cells of class T. 00016 00017 template <class T> 00018 class array { 00019 #ifdef INLIB_MEM 00020 static const std::string& s_class() { 00021 static const std::string s_v("inlib::array"); 00022 return s_v; 00023 } 00024 #endif 00025 public: 00026 typedef typename std::vector< std::pair<unsigned int,unsigned int> > cut_t; 00027 typedef typename std::vector<unsigned int> uints_t; 00028 typedef typename std::vector<T>::iterator vec_it_t; 00029 typedef typename std::vector<T>::const_iterator cons_vec_it_t; 00030 public: 00031 array() { 00032 #ifdef INLIB_MEM 00033 mem::increment(s_class().c_str()); 00034 #endif 00035 } 00036 array(const uints_t& a_orders) { 00037 #ifdef INLIB_MEM 00038 mem::increment(s_class().c_str()); 00039 #endif 00040 configure(a_orders); 00041 } 00042 array(unsigned int a_dimension,unsigned int a_order) { 00043 #ifdef INLIB_MEM 00044 mem::increment(s_class().c_str()); 00045 #endif 00046 // A hypercube of dimension "a_dimension" and size "a_order". 00047 uints_t orders(a_dimension); 00048 for(unsigned int index=0;index<a_dimension;index++) 00049 orders[index] = a_order; 00050 configure(orders); 00051 } 00052 virtual ~array() { 00053 #ifdef INLIB_MEM 00054 mem::decrement(s_class().c_str()); 00055 #endif 00056 } 00057 public: 00058 array(const array& a_from) 00059 :m_orders(a_from.m_orders) 00060 ,m_offsets(a_from.m_offsets) 00061 ,m_vector(a_from.m_vector) 00062 ,m_is(a_from.m_is){ 00063 #ifdef INLIB_MEM 00064 mem::increment(s_class().c_str()); 00065 #endif 00066 } 00067 array& operator=(const array& a_from) { 00068 m_orders = a_from.m_orders; 00069 m_offsets = a_from.m_offsets; 00070 m_vector = a_from.m_vector; 00071 m_is = a_from.m_is; 00072 return *this; 00073 } 00074 public: //operators: 00075 array& operator*=(const T& a_T) {multiply(a_T);return *this;} 00076 bool operator==(const array& a_array) const { 00077 return equal(a_array); 00078 } 00079 bool operator!=(const array& a_array) const { 00080 return !operator==(a_array); 00081 } 00082 array operator*(const T& a_T) const { 00083 array a(*this); 00084 a.multiply(a_T); 00085 return a; 00086 } 00087 // the below would need exception to do it properly. 00088 //array operator+(const array& a_array) const { 00089 // array a(*this); 00090 // if(!a.add(a_array)) {} //throw 00091 // return a; 00092 //} 00093 //array& operator+=(const array& a_array) { 00094 // if(!add(a_array)) {} //throw 00095 // return *this; 00096 //} 00097 public: 00098 void copy(const array& a_from) { 00099 m_orders = a_from.m_orders; 00100 m_offsets = a_from.m_offsets; 00101 m_vector = a_from.m_vector; 00102 m_is = a_from.m_is; 00103 } 00104 void clear() { 00105 m_orders.clear(); 00106 m_offsets.clear(); 00107 m_vector.clear(); 00108 m_is.clear(); 00109 } 00110 bool configure(const uints_t& a_orders) { 00111 m_orders = a_orders; 00112 unsigned int dim = m_orders.size(); 00113 if(dim==0) { 00114 clear(); 00115 return false; 00116 } 00117 unsigned int size = 1; 00118 for(unsigned int index=0;index<dim;index++) { 00119 //if(m_orders[index]<0) { 00120 //clear(); 00121 //return false; 00122 //} 00123 size *= m_orders[index]; 00124 } 00125 m_vector.resize(size,zero()); 00126 m_offsets.resize(dim,0); 00127 m_offsets[0] = 1; 00128 for(unsigned int iaxis=1;iaxis<dim;iaxis++) 00129 m_offsets[iaxis] = m_offsets[iaxis-1] * m_orders[iaxis-1]; 00130 m_is.resize(dim,0); 00131 return true; 00132 } 00133 unsigned int dimension() const { return m_orders.size();} 00134 const uints_t& orders() const { return m_orders;} 00135 unsigned int size() const {return m_vector.size();} 00136 bool set_value(const uints_t& a_is,const T& a_value) { 00137 unsigned int off = 0; 00138 if(!offset(a_is,off)) return false; 00139 m_vector[off] = a_value; 00140 return true; 00141 } 00142 bool value(const uints_t& a_is,T& a_value) const { 00143 unsigned int off = 0; 00144 if(!offset(a_is,off)) { 00145 a_value = 0; 00146 return false; 00147 } 00148 a_value = m_vector[off]; 00149 return true; 00150 } 00151 void reset() { 00152 for(vec_it_t it=m_vector.begin();it!=m_vector.end();++it) *it = 0; 00153 } 00154 const std::vector<T>& vector() const { return m_vector;} 00155 std::vector<T>& vector() { return m_vector;} 00156 bool fill(const std::vector<T>& a_values,cut_t* a_cut = 0) { 00157 unsigned int dsize = a_values.size(); 00158 unsigned di = 0; 00159 unsigned int index = 0; 00160 for(vec_it_t it=m_vector.begin();it!=m_vector.end();++it,index++) { 00161 if(!a_cut || (a_cut && accept(index,*a_cut)) ) { 00162 if(di>=dsize) return false; //a_values exhausted too early 00163 *it = a_values[di]; 00164 di++; 00165 } 00166 } 00167 return true; 00168 } 00169 // Related to other array : 00170 bool equal(const array& a_array) const { 00171 if(m_orders!=a_array.m_orders) return false; 00172 cons_vec_it_t it = m_vector.begin(); 00173 cons_vec_it_t ait = a_array.m_vector.begin(); 00174 for(;it!=m_vector.end();++it,++ait) { 00175 if((*it)!=(*ait)) return false; 00176 } 00177 return true; 00178 } 00179 bool equal(const array& a_array,T aEpsilon) const { 00180 if(m_orders!=a_array.m_orders) return false; 00181 cons_vec_it_t it = m_vector.begin(); 00182 cons_vec_it_t ait = a_array.m_vector.begin(); 00183 for(;it!=m_vector.end();++it,++ait) { 00184 T diff = (*it) - (*ait); 00185 if(diff<0) diff *= -1; 00186 if(diff>=aEpsilon) return false; 00187 } 00188 return true; 00189 } 00190 bool is_proportional(const array& a_array,T& a_factor) const { 00191 // If true, then : a_array = a_factor * this. 00192 a_factor = zero(); 00193 if(m_orders!=a_array.m_orders) return false; 00194 bool first = true; 00195 cons_vec_it_t it = m_vector.begin(); 00196 cons_vec_it_t ait = a_array.m_vector.begin(); 00197 for(;it!=m_vector.end();++it,++ait) { 00198 if( ((*it)==zero()) && ((*ait)==zero())) { 00199 continue; 00200 } else if( ((*it)!=zero()) && ((*ait)==zero())) { 00201 return false; 00202 } else if( ((*it)==zero()) && ((*ait)!=zero())) { 00203 return false; 00204 } else { 00205 if(first) { 00206 a_factor = (*ait)/(*it); 00207 first = false; 00208 } else { 00209 if((*ait)!=(*it)*a_factor) return false; 00210 } 00211 } 00212 } 00213 return true; 00214 } 00215 bool add(const array& a_array,cut_t* a_cut = 0) { 00216 if(m_orders!=a_array.m_orders) return false; 00217 vec_it_t it = m_vector.begin(); 00218 cons_vec_it_t ait = a_array.m_vector.begin(); 00219 unsigned int index = 0; 00220 for(;it!=m_vector.end();++it,++ait,index++) { 00221 if(!a_cut || (a_cut && accept(index,*a_cut)) ) { 00222 (*it) += (*ait); 00223 } 00224 } 00225 return true; 00226 } 00227 bool subtract(const array& a_array) { 00228 if(m_orders!=a_array.m_orders) return false; 00229 vec_it_t it = m_vector.begin(); 00230 cons_vec_it_t ait = a_array.m_vector.begin(); 00231 for(;it!=m_vector.end();++it,++ait) (*it) -= (*ait); 00232 return true; 00233 } 00234 bool multiply(const array& a_array) { 00235 if(m_orders!=a_array.m_orders) return false; 00236 vec_it_t it = m_vector.begin(); 00237 cons_vec_it_t ait = a_array.m_vector.begin(); 00238 for(;it!=m_vector.end();++it,++ait) (*it) *= (*ait); 00239 return true; 00240 } 00241 bool divide(const array& a_array) { 00242 if(m_orders!=a_array.m_orders) return false; 00243 bool status = true; 00244 vec_it_t it = m_vector.begin(); 00245 cons_vec_it_t ait = a_array.m_vector.begin(); 00246 for(;it!=m_vector.end();++it,++ait) { 00247 if((*ait)==zero()) { 00248 (*it) = zero(); //PAW convention. 00249 status = false; 00250 } else { 00251 (*it) /= (*ait); 00252 } 00253 } 00254 return status; 00255 } 00256 bool contract(const array& a_array,T& a_value) const { 00257 a_value = zero(); 00258 if(m_orders!=a_array.m_orders) return false; 00259 cons_vec_it_t it = m_vector.begin(); 00260 cons_vec_it_t ait = a_array.m_vector.begin(); 00261 for(;it!=m_vector.end();++it,++ait) a_value += (*it) * (*ait); 00262 return true; 00263 } 00264 //Else: 00265 void add(const T& a_T,cut_t* a_cut = 0) { 00266 unsigned int index = 0; 00267 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it,index++) { 00268 if(!a_cut || (a_cut && accept(index,*a_cut)) ) { 00269 (*it) += a_T; 00270 } 00271 } 00272 } 00273 void multiply(const T& a_T) { 00274 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it) (*it) *= a_T; 00275 } 00276 bool divide(const T& a_T) { 00277 if(a_T==zero()) return false; 00278 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it) (*it) /= a_T; 00279 return true; 00280 } 00281 bool invert() { 00282 bool status = true; 00283 T v_one = one(); 00284 for(vec_it_t it = m_vector.begin();it!=m_vector.end();++it) { 00285 if((*it)==zero()) { 00286 (*it) = zero(); //PAW convention. 00287 status = false; 00288 } else { 00289 T value = (*it); 00290 (*it) = v_one/value; 00291 } 00292 } 00293 return status; 00294 } 00295 bool offset(const uints_t& a_is,unsigned int& a_offset) const { 00296 unsigned int dim = m_orders.size(); 00297 a_offset = 0; 00298 if(a_is.size()!=dim) return false; 00299 if(dim==0) return false; 00300 for(unsigned int iaxis=0;iaxis<dim;iaxis++) { 00301 unsigned int i = a_is[iaxis]; 00302 if(i>=m_orders[iaxis]) { 00303 a_offset = 0; 00304 return false; 00305 } 00306 a_offset += i * m_offsets[iaxis]; 00307 } 00308 return true; 00309 } 00310 bool indices(unsigned int a_offset,uints_t& a_is) const { 00311 if(a_offset>=m_vector.size()) return false; 00312 unsigned int dim = m_orders.size(); 00313 unsigned int off = a_offset; 00314 for(int iaxis=dim-1;iaxis>=0;iaxis--) { 00315 a_is[iaxis] = off/m_offsets[iaxis]; 00316 off -= a_is[iaxis] * m_offsets[iaxis]; 00317 } 00318 return true; 00319 } 00320 bool accept(unsigned int a_index,const cut_t& a_cut) const { 00321 unsigned int dim = m_orders.size(); 00322 if(a_cut.size()!=dim) return false; 00323 if(!indices(a_index,const_cast<uints_t&>(m_is))) return false; 00324 bool good = true; 00325 for(unsigned iaxis=0;iaxis<dim;iaxis++) { 00326 if(m_is[iaxis]<a_cut[iaxis].first) {good = false;break;} 00327 if(m_is[iaxis]>a_cut[iaxis].second) {good = false;break;} 00328 } 00329 return good; 00330 } 00331 static T zero() { 00332 //return (T)0; 00333 //return T(0); 00334 return T(); 00335 } 00336 static T one() { 00337 //return (T)1; 00338 return T(1); 00339 } 00340 static T minus_one() { 00341 //return (T)-1; 00342 return T(-1); 00343 } 00344 static T two() { 00345 //return (T)2; 00346 return T(2); 00347 } 00348 protected: 00349 uints_t m_orders; 00350 uints_t m_offsets; 00351 std::vector<T> m_vector; 00352 uints_t m_is; 00353 }; 00354 00355 } 00356 00357 #include <ostream> 00358 00359 // Helpers array<T> : 00360 namespace inlib { 00361 00362 template <class T> 00363 inline bool contract( 00364 const array<T>& aA 00365 ,unsigned int aIA 00366 ,const array<T>& aB 00367 ,unsigned int aIB 00368 ,array<T>& aR 00369 ){ 00370 // Contract a tensor (aA) with a vector (aB) 00371 // on indice aIA of aA and aIB of aB. 00372 // aA.dimension must be > 1. 00373 // aB.dimension must be > 1. 00374 if(&aR==&aA) return false; 00375 if(&aR==&aB) return false; 00376 00377 if(aA.dimension()==0) return false; 00378 if(aB.dimension()==0) return false; 00379 if(aIA>=aA.dimension()) return false; 00380 if(aIB>=aB.dimension()) return false; 00381 if(aA.orders()[aIA]!=aB.orders()[aIB]) return false; 00382 00383 unsigned int rdima = aA.dimension()-1; 00384 unsigned int rdimb = aB.dimension()-1; 00385 00386 if((rdima+rdimb)==0) return false; 00387 00388 std::vector<unsigned int> rorders(rdima+rdimb); 00389 unsigned int index; 00390 for(index=0;index<aIA;index++) 00391 rorders[index] = aA.orders()[index]; 00392 for(index=aIA;index<rdima;index++) 00393 rorders[index] = aA.orders()[index+1]; 00394 00395 for(index=0;index<aIB;index++) 00396 rorders[rdima+index] = aB.orders()[index]; 00397 for(index=aIB;index<rdimb;index++) 00398 rorders[rdima+index] = aB.orders()[index+1]; 00399 00400 if(!aR.configure(rorders)) return false; 00401 00402 std::vector<unsigned int> ais(aA.dimension()); 00403 std::vector<unsigned int> bis(aB.dimension()); 00404 std::vector<unsigned int> ris(aR.dimension()); 00405 00406 //FIXME : optimize the below. 00407 unsigned int order = aA.orders()[aIA]; 00408 unsigned int rsize = aR.size(); 00409 00410 std::vector<T>& rvec = aR.vector(); 00411 00412 for(unsigned int roffset=0;roffset<rsize;roffset++) { 00413 if(!aR.indices(roffset,ris)) return false; 00414 00415 for(index=0;index<aIA;index++) ais[index] = ris[index]; 00416 for(index=aIA;index<rdima;index++) ais[index+1] = ris[index]; 00417 00418 for(index=0;index<aIB;index++) bis[index] = ris[rdima+index]; 00419 for(index=aIB;index<rdimb;index++) bis[index+1] = ris[rdima+index]; 00420 00421 T value = 0; 00422 for(index=0;index<order;index++) { 00423 ais[aIA] = index; 00424 T av; 00425 if(!aA.value(ais,av)) return false; 00426 00427 bis[aIB] = index; 00428 T bv; 00429 if(!aB.value(bis,bv)) return false; 00430 00431 value += av * bv; 00432 } 00433 00434 rvec[roffset] = value; 00435 } 00436 00437 return true; 00438 } 00439 00440 template <class T> 00441 inline bool swap( 00442 const array<T>& aV 00443 ,unsigned int aI1 00444 ,unsigned int aI2 00445 ,array<T>& aR 00446 ){ 00447 if(&aR==&aV) return false; 00448 00449 unsigned int dim = aV.dimension(); 00450 if(aI1>=dim) return false; 00451 if(aI2>=dim) return false; 00452 00453 if(aI1==aI2) { 00454 aR.copy(aV); 00455 return true; 00456 } 00457 00458 if(!aR.configure(aV.orders())) return false; 00459 00460 std::vector<unsigned int> vis(aV.dimension()); 00461 std::vector<unsigned int> ris(aV.dimension()); 00462 00463 const std::vector<T>& vvec = aV.vector(); 00464 00465 unsigned int size = aV.size(); 00466 for(unsigned int offset=0;offset<size;offset++) { 00467 T value = vvec[offset]; 00468 00469 if(!aV.indices(offset,vis)) return false; 00470 unsigned int i = vis[aI1]; 00471 vis[aI1] = vis[aI2]; 00472 vis[aI2] = i; 00473 if(!aR.set_value(vis,value)) return false; 00474 } 00475 00476 return true; 00477 } 00478 00479 //NOTE : print is a Python keyword. 00480 template <class T> 00481 inline void dump(std::ostream& a_out,const inlib::array<T>& a_array,const std::string& a_title){ 00482 if(a_title.size()) a_out << a_title << std::endl; 00483 const std::vector<T>& vec = a_array.vector(); 00484 typedef typename std::vector<T>::const_iterator cons_vec_it_t; 00485 for(cons_vec_it_t it = vec.begin();it!=vec.end();++it) { 00486 a_out << (*it) << std::endl; 00487 } 00488 } 00489 00490 template <class T> 00491 inline void diff(std::ostream& a_out,const array<T>& aA,const array<T>& aB,T a_epsilon){ 00492 if(aA.orders()!=aB.orders()) { 00493 a_out << "inlib::arrays::diff : not same orders" << std::endl; 00494 return; 00495 } 00496 bool header_done = false; 00497 unsigned int dim = aA.dimension(); 00498 std::vector<unsigned int> is(dim); 00499 unsigned int vsize = aA.vector().size(); 00500 for(unsigned int index=0;index<vsize;index++) { 00501 T diff = aA.vector()[index]-aB.vector()[index]; 00502 if(diff<0) diff *= -1; 00503 if(diff>=a_epsilon) { 00504 aA.indices(index,is); 00505 if(!header_done) { 00506 a_out << "inlib::arrays::diff :" << std::endl; 00507 header_done = true; 00508 } 00509 for(unsigned int i=0;i<dim;i++) a_out << is[i] << " "; 00510 a_out << aA.vector()[index] << " " << aB.vector()[index] << std::endl; 00511 } 00512 } 00513 } 00514 00518 00519 template <class T> 00520 class kronecker : public array<T> { 00521 public: 00522 kronecker(unsigned int a_order):array<T>(a_order,a_order){ 00523 // epsilon(i1,i2,....in) with n = a_order and i in [0,n[. 00524 // Then an Array of a_order * a_order. 00525 unsigned int index = 0; 00526 typedef typename std::vector<unsigned int> uints_t; 00527 uints_t is(a_order); 00528 std::vector<T>& vec = array<T>::vector(); 00529 typedef typename std::vector<T>::iterator vec_it_t; 00530 vec_it_t it = vec.begin(); 00531 for(;it!=vec.end();++it,index++) { 00532 if(!array<T>::indices(index,is)) return; //FIXME throw. 00533 bool good = true; 00534 {for(unsigned iaxis=0;iaxis<a_order;iaxis++) { 00535 unsigned int ival = is[iaxis]; 00536 for(unsigned iaxis2=iaxis+1;iaxis2<a_order;iaxis2++) { 00537 if(is[iaxis2]==ival) { 00538 good = false; 00539 break; 00540 } 00541 } 00542 if(!good) break; 00543 }} 00544 if(!good) continue; 00545 // All indicies are different. 00546 unsigned int n = 0; 00547 for(unsigned iaxis=0;iaxis<a_order;) { 00548 unsigned int ival = is[iaxis]; 00549 if(ival!=iaxis) { 00550 // Swap and add one permutation : 00551 unsigned int old = is[ival]; 00552 is[ival] = ival; 00553 is[iaxis] = old; 00554 n +=1; 00555 } else { 00556 iaxis++; 00557 } 00558 } 00559 {unsigned int n_2 = n/2; 00560 if(2*n_2==n) (*it) = array<T>::one(); 00561 else (*it) = array<T>::minus_one();} 00562 } 00563 } 00564 }; 00565 00566 } 00567 00568 #endif