introduction  —  namespaces  —  modules  —  classes  —  files  —  globals  —  members  —  examples  —  Marc Toussaint

array.h

Go to the documentation of this file.
00001 /*  Copyright (C) 2000, 2006  Marc Toussaint (mtoussai@inf.ed.ac.uk)
00002     under the terms of the GNU LGPL (http://www.gnu.org/copyleft/lesser.html)
00003     see the `std.h' file for a full copyright statement  */
00004 
00007 #ifndef MT_array_h
00008 #define MT_array_h
00009 
00010 //#define MT_EXPRESSIONS
00011 
00012 #include"std.h"
00013 
00014 #define FOR1D(x,i) for(i=0;i<x.d0;i++)
00015 #define FOR2D(x,i,j) for(i=0;i<x.d0;i++) for(j=0;j<x.d1;j++)
00016 
00017 #define forAll(i,A) for(i=A.p;i!=A.pstop;i++)
00018 
00019 //===========================================================================
00020 //
00021 // handle expressions of arrays more efficiently -- PRELIMINARY
00022 //
00023 
00024 #ifdef MT_EXPRESSIONS
00025 namespace MT{
00026 enum ExOp { UNI, PLUS, MINUS, MUL, Div, PROD, OUT };
00027 struct Ex{
00028   ExOp op;
00029   double add,mul;
00030   bool trans;
00031   void *A,*B;
00032   Ex(MT::ExOp OP,void *AA,void *BB){ op=OP; A=AA; B=BB; }
00033   Ex(void *_A,double _add=0,double _mul=1,bool _trans=false){ op=UNI; A=_A; add=_add; mul=_mul; trans=_trans; }
00034 };
00035 }
00036 #  define IFEX(x) x
00037 #else
00038 #  define IFEX(x)
00039 #endif
00040 
00041 
00042 //===========================================================================
00043 //
00044 // Array class
00045 //
00046 
00047 namespace MT{
00054 template<class T>
00055 class Array{
00056 public:
00057   T *p;     
00058   T *pstop; 
00059   uint N;   
00060   uint nd;  
00061   uint d0;  
00062   uint d1;  
00063   uint d2;  
00064   uint memSize;
00065   bool reference;
00066   T **pp;   
00067   IFEX( Ex *ex; )
00068   bool temp;
00069 
00070 public:
00071 
00076   bool memMove;
00077 
00083   bool flexiMem;
00085 private:
00086   static char memMoveInit;
00087   static int sizeT;
00088   void init(){
00089     reference=false;
00090     memMove=false;
00091     if(sizeT==-1) sizeT=sizeof(T);
00092     if(memMoveInit==-1){
00093       memMoveInit=0;
00094       if(typeid(T)==typeid(bool) ||
00095   typeid(T)==typeid(char) ||
00096   typeid(T)==typeid(unsigned char) ||
00097   typeid(T)==typeid(int) ||
00098   typeid(T)==typeid(unsigned int) ||
00099   typeid(T)==typeid(short) ||
00100   typeid(T)==typeid(unsigned short) ||
00101   typeid(T)==typeid(long) ||
00102   typeid(T)==typeid(unsigned long) ||
00103   typeid(T)==typeid(float) ||
00104   typeid(T)==typeid(double)) memMoveInit=1;
00105     }
00106     memMove=(memMoveInit==1);
00107     flexiMem=true;
00108     p=pstop=0;
00109     memSize=N=nd=d0=d1=d2=0;
00110     pp=0;
00111     IFEX( ex=0; );
00112     temp=false;
00113   }
00114 
00115 public:
00116 
00118   Array(){ init(); }
00119 
00121   Array(const Array<T>& a){ init(); operator=(a); }
00122 
00124   Array(uint i){ init(); resize(i); }
00125 
00127   Array(uint i,uint j){ init(); resize(i,j); }
00128 
00130   Array(uint i,uint j,uint k){ init(); resize(i,j,k); }
00131 
00133   Array(const Array<T>& a,uint dim){ init(); referToSubDim(a,dim); }
00134 
00136   Array(const Array<T>& a,uint i,uint j){ init(); referToSubDim(a,i,j); }
00137 
00139   Array(T* p,uint size){ init(); referTo(p,size); }
00140 
00142   Array(const char* str){ init(); set(str); }
00143 
00144   IFEX( Array(Ex *_ex){ ex=_ex; } )
00145 
00146   ~Array(){
00147     IFEX( if(ex){ delete ex; return; } )
00148     freeMEM();
00149   }
00150   
00151 
00152 public:
00153 
00155   void clear(){ resize(0); }
00156 
00158   void resize(uint D0){ nd=1; d0=D0; resizeMEM(d0,false); }
00159   void resizeCopy(uint D0){ nd=1; d0=D0; resizeMEM(d0,true); }
00160   void resizeZero(uint D0){ nd=1; d0=D0; resizeMEM(d0,false); setZero(); }
00161   void reshape(uint D0){ CHECK(N==D0,"reshape must preserve total memory size"); nd=1; d0=D0; d1=d2=0; }
00162 
00164   void resize(uint D0,uint D1){ nd=2; d0=D0; d1=D1; resizeMEM(d0*d1,false); }
00165   void resizeCopy(uint D0,uint D1){ nd=2; d0=D0; d1=D1; resizeMEM(d0*d1,true); }
00166   void reshape(uint D0,uint D1){ CHECK(N==D0*D1,"reshape must preserve total memory size"); nd=2; d0=D0; d1=D1; d2=0; }
00167 
00169   void resize(uint D0,uint D1,uint D2){ nd=3; d0=D0; d1=D1; d2=D2; resizeMEM(d0*d1*d2,false); }
00170   void resizeCopy(uint D0,uint D1,uint D2){ nd=3; d0=D0; d1=D1; d2=D2; resizeMEM(d0*d1*d2,true); }
00171   void reshape(uint D0,uint D1,uint D2){ CHECK(N==D0*D1*D2,"reshape must preserve total memory size"); nd=3; d0=D0; d1=D1; d2=D2; }
00172 
00173 #ifndef MT_MSVC
00174 
00175   void resize(const Array<T>& a){ nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2; resizeMEM(a.N,false); }
00176   void resizeCopy(const Array<T>& a){ nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2; resizeMEM(a.N,true); }
00177 #endif
00178 
00180   template<class S>
00181   void resize(const Array<S>& a){ nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2; resizeMEM(a.N,false); }
00182   template<class S>
00183   void resizeCopy(const Array<S>& a){ nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2; resizeMEM(a.N,true); }
00185 
00186 private: //core rountine to allocate memory
00187   void resizeMEM(uint size,bool copy){
00188     if(reference){
00189       CHECK(size==N,"real resize of subarray is not allowed! (only a resize without changing memory size)");
00190       return;
00191     }
00192     if(N==size) return;
00193     uint i;
00194     T *pold=p;
00195     if(!flexiMem){
00196       if(memSize!=size){
00197   if(size) p=new T [size];
00198   if(!p) HALT("MT::Array failed memory allocation of "<<size*sizeT<<"bytes");
00199   if(copy && !memMove) for(i=N<size?N:size;i--;) p[i]=pold[i];
00200   if(copy && memMove) memmove(p,pold,sizeT*(N<size?N:size));
00201   if(memSize) delete[] pold;
00202   N=memSize=size;
00203   pstop=p+N;
00204       }
00205     }else{
00206       if(size>0 && memSize==0){ //first time
00207   p=new T [size];
00208   if(!p) HALT("MT::Array failed memory allocation of "<<size*sizeT<<"bytes");
00209   memSize=size;
00210       }else if(size>memSize || 10+2*size<memSize/2){
00211   uint oversize=10+2*size;
00212   p=new T [oversize];
00213   if(!p) HALT("MT::Array failed memory allocation of "<<size*sizeT<<"bytes");
00214   if(copy && !memMove) for(i=N<size?N:size;i--;) p[i]=pold[i];
00215   if(copy && memMove) memmove(p,pold,sizeT*(N<size?N:size));
00216   if(memSize) delete[] pold;
00217   memSize=oversize;
00218       }
00219       N=size;
00220       pstop=p+N;
00221     }
00222     if(pp){ delete[] pp; pp=0; }
00223   }
00224   void freeMEM(){
00225     if(memSize) delete[] p;
00226     if(pp) delete[] pp;
00227     p=pstop=0;
00228     memSize=N=nd=d0=d1=d2=0;
00229   }
00230 
00231 
00232 public:
00233 
00235   void append(const T& x){
00236     if(nd==2 && d1==1)
00237       resizeCopy(d0+1,d1);
00238     else
00239       resizeCopy(N+1);
00240     p[N-1]=x;
00241   }
00243   void append(const Array<T>& x){
00244     uint oldN=N,i;
00245     if(nd==2 && d1==x.N) 
00246       resizeCopy(d0+1,d1);
00247     else
00248       resizeCopy(N+x.N);
00249     if(!memMove){
00250       for(i=0;i<x.N;i++) p[oldN+i]=x.p[i];
00251     }else memmove(p+oldN,x.p,sizeT*x.N);
00252   }
00254   void append(const T *p,uint n){
00255     uint oldN=N,i;
00256     if(nd==2 && d1==n) 
00257       resizeCopy(d0+1,d1);
00258     else
00259       resizeCopy(N+n);
00260     if(!memMove){
00261       for(i=0;i<n;i++) p[n+i]=p[i];
00262     }else memmove(p+oldN,p,sizeT*n);
00263   }
00264 
00265   /* inserts x at the position i -- the array becomes 1D! [only with memMove!] */
00266   void insert(uint i,const T& x){
00267     CHECK(memMove,"only with memMove");
00268     uint Nold=N;
00269     resizeCopy(Nold+1);
00270     if(i<Nold) memmove(p+i+1,p+i,sizeT*(Nold-i));
00271     p[i]=x;
00272   }
00273 
00275   void remove(uint i,uint n=1){
00276     CHECK(memMove,"only with memMove");
00277     if(N>i+n) memmove(p+i,p+i+n,sizeT*(N-i-n));
00278     resizeCopy(N-n);
00279   }
00280 
00282   void replace(uint i,uint n,const Array<T>& x){
00283     CHECK(memMove,"only with memMove");
00284     uint Nold=N;
00285     if(n==x.N){
00286       memmove(p+i,x.p,sizeT*(x.N));
00287     }else if(n>x.N){
00288       memmove(p+i+x.N,p+i+n,sizeT*(Nold-i-n));
00289       if(i+n<Nold) memmove(p+i,x.p,sizeT*(x.N));
00290       resizeCopy(Nold-n+x.N);
00291     }else{
00292       resizeCopy(Nold+x.N-n);
00293       if(i+n<Nold) memmove(p+i+x.N,p+i+n,sizeT*(Nold-i-n));
00294       memmove(p+i,x.p,sizeT*(x.N));
00295     }
00296   }
00297 
00299   void delRow(uint i){
00300     CHECK(memMove,"only with memMove");
00301     CHECK(i<d0,"range check error");
00302     uint n=d1;
00303     if(i+1<d0) memmove(p+i*n,p+(i+1)*n,sizeT*(d0-i-1)*n);
00304     resizeCopy(d0-1,n);
00305   }
00306 
00308   void delColumn(uint i){
00309     CHECK(memMove,"only with memMove");
00310     CHECK(i<d1,"range check error");
00311     uint n=d1;
00312     for(uint j=0;j<d0;j++){
00313       memmove(p+j*(n-1)  ,p+j*n  ,sizeT*i);
00314       memmove(p+j*(n-1)+i,p+j*n+(i+1),sizeT*(n-i-1));
00315     }
00316     resizeCopy(d0,n-1);
00317   }
00318   
00319 
00320 
00321 public:
00322 
00323   T& elem(uint i) const{
00324     CHECK(i<N,"range error ("<<i<<">="<<N<<")");
00325     return p[i];
00326   }
00327 
00328   /* scalar reference (legal iff N==1) */
00329   /*operator T&() const{ 
00330     CHECK(N==1,"scalar reference ("<<N<<"!=1)");
00331     return *e; }*/
00332 
00334   T& operator()(uint i) const{ 
00335     CHECK(nd==1 && i<d0,"1D range error ("<<nd<<"=1,"<<i<<"<"<<d0<<")");
00336     return p[i]; }
00337 
00339   T& operator()(uint i,uint j) const{
00340     CHECK(nd==2 && i<d0 && j<d1,
00341     "2D range error ("<<nd<<"=2,"<<i<<"<"<<d0<<","<<j<<"<"<<d1<<")");
00342     return p[i*d1+j];
00343   }
00344 
00346   T& operator()(uint i,uint j,uint k) const{
00347     CHECK(nd==3 && i<d0 && j<d1 && k<d2,
00348     "3D range error ("<<
00349     nd<<"=3,"<<i<<"<"<<d0<<","<<j<<"<"<<d1<<","<<k<<"<"<<d2<<")");
00350     return p[(i*d1+j)*d2+k];
00351   }
00352 
00354   Array<T> operator[](uint i) const{ return Array(*this,i); }
00355 
00357   Array<T>& operator()(){ return (*this); }
00358 
00361   uint maxIndex(){ uint i,m=0; for(i=0;i<N;i++) if(p[i]>p[m]) m=i; return m; }
00362 
00365   void maxIndex(uint& i,uint& j){ CHECK(nd==2,"needs 2D array"); j=maxIndex(); i=j/d1; j=j%d1; }
00366 
00369   void maxIndex(uint& i,uint& j,uint& k){ CHECK(nd==3,"needs 3D array"); k=maxIndex(); i=k/(d1*d2); k=k%(d1*d2); j=k/d2; k=k%d2; }
00370 
00372   uint minIndex(){ uint i,m=0; for(i=0;i<N;i++) if(p[i]<p[m]) m=i; return m; }
00373 
00374   bool contains(T& x){ uint i; for(i=0;i<N;i++) if(p[i]==x) return true; return false; }
00375 
00379   Array<T> sub(uint i,int I) const{
00380     CHECK(nd==1,"1D range error ");
00381     Array<T> x;
00382     if(I==-1) I=d0-1;
00383     CHECK(i<=I,"lower limit higher than upper!");
00384     x.resize(I-i+1);
00385     uint k;
00386     for(k=i;k<=I;k++) x(k-i)=operator()(k);
00387     return x;
00388   }
00389 
00393   Array<T> sub(uint i,int I,uint j,int J) const{
00394     CHECK(nd==2,"2D range error ");
00395     Array<T> x;
00396     if(I==-1) I=d0-1;
00397     if(J==-1) J=d1-1;
00398     CHECK(i<(uint)I && j<(uint)J,"lower limit must be higher than upper!");
00399     x.resize(I-i+1,J-j+1);
00400     uint k,l;
00401     for(k=i;k<=(uint)I;k++) for(l=j;l<=(uint)J;l++) x(k-i,l-j)=operator()(k,l);
00402     return x;
00403   }
00404 
00409   Array<T> sub(uint i,int I,Array<uint> cols) const{
00410     CHECK(nd==2,"2D range error ");
00411     Array<T> x;
00412     if(I==-1) I=d0-1;
00413     CHECK(i<=(uint)I,"lower limit higher than upper!");
00414     x.resize(I-i+1,cols.N);
00415     uint k,l;
00416     for(k=i;k<=(uint)I;k++) for(l=0;l<cols.N;l++) x(k-i,l)=operator()(k,cols(l));
00417     return x;
00418   }
00419   
00420 public:
00421 
00422   // allocates, sets and returns the pp-pointer (of type T**)
00423   T** getCarray(){
00424     CHECK(nd==2,"only 2D array gives C-array of type T**");
00425     if(pp) return pp;
00426     pp=new T* [d0];
00427     for(uint i=0;i<d0;i++) pp[i]=p+i*d1;
00428     return pp;
00429   }
00430 
00433   T** getPointers(Array<T*>& array2d) const{
00434     CHECK(nd==2,"only 2D array gives C-array of type T**");
00435     array2d.resize(d0);
00436     for(uint i=0;i<d0;i++) array2d(i)=p+i*d1;
00437     return array2d.p;
00438   }
00439 
00441   T*** getPointers(Array<T**>& array3d,Array<T*>& array2d) const{
00442     CHECK(nd==3,"only 3D array gives C-array of type T*** ");
00443     array2d.resize(d0,d1);
00444     for(uint i=0;i<d0;i++){
00445       for(uint j=0;j<d1;j++) array2d(i,j)=&operator()(i,j,0);
00446       array3d(i)=&array2d(i,0);
00447     }
00448     return array3d.p;
00449   }
00450 
00451 
00452 
00453 public:
00454 
00456   Array<T>& operator=(const T v){
00457     uint i;
00458     //if(memMove && typeid(T)==typeid(T)) memset(p,*((int*)&v),N); else
00459     for(i=0;i<N;i++) p[i]=v;
00460     return *this;
00461   }
00462 
00464   Array<T>& operator=(const Array<T>& a){
00465     IFEX( if(a.ex){ assignEx(a); return *this; } );
00466     //if(a.temp){ takeOver(*((Array<T>*)&a)); return *this; }
00467     resize(a);
00468     uint i;
00469     if(memMove) memmove(p,a.p,sizeT*N);
00470     else for(i=0;i<N;i++) p[i]=a.p[i];
00471     return *this;
00472   }
00473 
00475   void set(const char* str){
00476     std::istringstream s(str);
00477     readRaw(s);
00478   }
00479 
00480 #ifdef MT_EXPRESSIONS
00481   void checkEx() const{ if(ex) assignEx(); }
00482   void assignEx(const Array<T>& a) const;
00483   void assignEx() const;
00484   //Array<T>& operator=(const Ex& ex);
00485 #else
00486   void checkEx() const{}
00487 #endif
00488 
00490   template<class S>
00491   Array<T>& copy(const Array<S>& a){
00492     resize(a);
00493     uint i;
00494     if(memMove && typeid(T)==typeid(S)) memmove(p,a.p,sizeT*N);
00495     else for(i=0;i<N;i++) p[i]=(T)a.p[i];
00496     return *this;
00497   }
00498 
00501   void setZero(byte zero=0){
00502     CHECK(memMove,"can set array's memory to zero only if memMove option is true");
00503     memset(p,zero,sizeT*N);
00504   }
00505 
00507   template<class S>
00508   void copy(S *buffer,uint D0){
00509     resize(D0);
00510     uint i;
00511     if(memMove && typeid(T)==typeid(S))
00512       memmove(p,buffer,sizeT*d0);
00513     else for(i=0;i<d0;i++) operator()(i)=(T)buffer[i];
00514   }
00515 
00517   template<class S>
00518   void copy(S **buffer,uint D0,uint D1){
00519     resize(D0,D1);
00520     uint i,j;
00521     for(i=0;i<d0;i++){
00522       if(memMove && typeid(T)==typeid(S))
00523   memmove(p+i*d1,buffer[i],sizeT*d1);
00524       else for(j=0;j<d1;j++) operator()(i,j)=(T)buffer[i][j];
00525     }
00526   }
00527 
00530   template<class S>
00531   void copyInto(S *buffer){
00532     CHECK(nd==1,"can only copy 1D Array into 1D C-array");
00533     uint i;
00534     if(memMove && typeid(T)==typeid(S)) memmove(buffer,p,sizeT*d0);
00535     else for(i=0;i<d0;i++) buffer[i]=(S)operator()(i);
00536   }
00537 
00540   template<class S>
00541   void copyInto2D(S **buffer){
00542     CHECK(nd==2,"can only copy 2D Array into 2D C-array");
00543     uint i,j;
00544     for(i=0;i<d0;i++){
00545       if(memMove && typeid(T)==typeid(S)) memmove(buffer[i],p+i*d1,sizeT*d1);
00546       else for(j=0;j<d1;j++) buffer[i][j]=(S)operator()(i,j);
00547     }
00548   }
00549 
00551   void referTo(const T *buffer,uint n){
00552     freeMEM();
00553     reference=true;
00554     nd=1; d0=n; d1=d2=0; N=n;
00555     p=(T*)buffer;
00556     pstop=p+N;
00557   }
00558 
00560   void referTo(const Array<T>& a){
00561     freeMEM();
00562     reference=true; memMove=a.memMove; flexiMem=a.flexiMem;
00563     N=a.N; nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2;
00564     p=a.p;
00565     pstop=a.pstop;
00566   }
00567 
00569   void referToSubRange(const Array<T>& a,uint i,int I){
00570     CHECK(a.nd==1,"not implemented yet");
00571     freeMEM();
00572     reference=true; memMove=a.memMove; flexiMem=a.flexiMem;
00573     if(I==-1) I=a.N-1;
00574     N=I+1-i; nd=1; d0=N; d1=0; d2=0;
00575     p=a.p+i;
00576     pstop=p+N;
00577   }
00578 
00580   void referToSubDim(const Array<T>& a,uint dim){
00581     CHECK(a.nd>1,"can't create subarray of array less than 2 dimensions");
00582     freeMEM();
00583     reference=true; memMove=a.memMove; flexiMem=a.flexiMem;
00584     if(a.nd==2){
00585       nd=1; d0=a.d1; d1=d2=0; N=d0;
00586       p=&a(dim,0);
00587       pstop=p+N;
00588     }
00589     if(a.nd==3){
00590       nd=2; d0=a.d1; d1=a.d2; d2=0; N=d0*d1;
00591       p=&a(dim,0,0);
00592       pstop=p+N;
00593     }
00594   }
00595 
00597   void referToSubDim(const Array<T>& a,uint i,uint j){
00598     CHECK(a.nd>2,"can't create subsubarray of array less than 3 dimensions");
00599     freeMEM();
00600     reference=true; memMove=a.memMove; flexiMem=a.flexiMem;
00601     if(a.nd==3){
00602       nd=1; d0=a.d2; d1=0; d2=0; N=d0;
00603       p=&a(i,j,0);
00604       pstop=p+N;
00605     }
00606   }
00607 
00609   void takeOver(Array<T>& a){
00610     freeMEM();
00611     memMove=a.memMove; flexiMem=a.flexiMem;
00612     N=a.N; nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2;
00613     p=a.p;
00614     pstop=a.pstop;
00615     memSize=a.memSize;
00616     a.reference=true;
00617     a.memSize=0;
00618   }
00619 
00620 
00621 public:
00622   void First(T*& i) const{ i=p; }
00623   void Last(T*& i) const{ i=(p+N)-1; }
00624   void Next(T*& i) const{ i++; }
00625   bool CheckNext(T*& i) const{ return i!=pstop; }
00626   void Prev(T*& i) const{ i--; }
00627   bool CheckPrev(T*& i) const{ return i!=p-1; }
00628 
00629 
00630 public:
00631 
00632   void permute(uint i,uint j){ T x=p[i]; p[i]=p[j]; p[j]=x; }
00633 
00635   void randomPermute(){
00636     int j,r;
00637     for(j=N-1;j>=1;j--){
00638       r=rnd(j+1);
00639       permute(r,j);
00640     }
00641   }
00642 
00643 
00644 public:
00645 
00647   void write(std::ostream& os,char *ARRAYSEP=" ",bool binary=false) const{
00648     CHECK(!binary || memMove,"binary write works only for memMoveable data");
00649     uint i,j,k;
00650     if(!IOraw){//write tag
00651       if(nd==0){ os <<"<0>"; return; }
00652       if(nd==1)  os <<"<1:" <<d0;
00653       if(nd==2)  os <<"<2:" <<d0 <<',' <<d1;
00654       if(nd==3)  os <<"<3:" <<d0 <<',' <<d1 <<',' <<d2;
00655     }
00656     //write data
00657     if(binary){
00658       os.write((char*)p,sizeT*N);
00659     }else{
00660       if(nd==1){
00661   //os <<' ';
00662   for(i=0;i<d0;i++) os <<ARRAYSEP <<operator()(i);
00663       }
00664       if(nd==2) for(j=0;j<d0;j++){
00665   if(j || !IOraw) os <<"\n";
00666   for(i=0;i<d1;i++) os <<ARRAYSEP <<operator()(j,i);
00667       }
00668       if(nd==3) for(k=0;k<d0;k++){
00669         if(k || !IOraw) os <<"\n";
00670         for(j=0;j<d1;j++){
00671           for(i=0;i<d2;i++) os <<ARRAYSEP <<operator()(k,j,i);
00672           os <<"\n";
00673         }
00674       }
00675     }
00676     if(!IOraw){//write end tag
00677       os <<">";
00678       if(nd>1) os <<std::endl;
00679     }
00680   }
00681 
00683   void read(std::istream& is,bool binary=false){
00684     CHECK(!binary || memMove,"binary write works only for memMoveable data");
00685     uint d,i,j,k;
00686     if(!IOraw){//read tag
00687       is >>"<" >>d;
00688       if(is.fail()) HALT ("could not read array tag");
00689       if(d==0){ is >>">"; return; }
00690       if(d==1){
00691         is >>":" >>i;
00692         if(is.fail()) HALT ("could not read array's dimensions");
00693   resize(i);
00694       }
00695       if(d==2){
00696         is >>":" >>i >>"," >>j;
00697   if(is.fail()) HALT ("could not read array's dimensions");
00698         resize(i,j);
00699       }
00700       if(d==3){
00701   is >>":" >>i >>"," >>j >>"," >>k;
00702   if(is.fail()) HALT ("could not read array's dimensions");
00703         resize(i,j,k);
00704       }
00705     }
00706     //read data
00707     if(binary){
00708       is.read((char*)p,sizeT*N);
00709     }else{
00710       for(i=0;i<N;i++){
00711   if(is.fail()) HALT("could not read "<<i<<"-th element of an array");
00712   is >>p[i];
00713       }
00714     }
00715     //read end tag
00716     if(!IOraw){
00717       is >>">";
00718       if(is.fail()) HALT ("could not read array end tag");
00719     }
00720   }
00721 
00723   const Array<T>& ioraw() const{ IOraw=true; return *this; }
00724 
00726   const Array<T>& ionoraw() const{ IOraw=false; return *this; }
00727 
00729   void readRaw(std::istream& is,uint d1=0,uint d2=0){
00730     resize(0);
00731     T x;
00732     for(;;){
00733       MT::skip(is);
00734       if(!is.good()) break;
00735       is >>x;
00736       if(is.fail()) HALT("error when reading "<<N<<"th element in readRaw");
00737       if(d1 && d2 && N>=d1*d2) break;
00738       append(x);
00739     }
00740     if(d1 && N>d1){
00741       CHECK(!(N%d1),"read "<<N<<"elements - not dividable by "<<d1);
00742       resizeCopy(N/d1,d1);
00743     }
00744   }
00745   void readRaw(const char* filename,uint d1=0,uint d2=0){
00746     std::ifstream is; MT::open(is,filename); readRaw(is,d1,d2);
00747   }
00748   void writeRaw(const char* filename){
00749     std::ofstream os; MT::open(os,filename); os <<(*this).ioraw();
00750   }
00752 
00753 #ifndef MT_EXPRESSIONS
00754 
00755   friend inline Array<T> operator~(const Array<T>& y){ Array<T> x; transpose(x,y); x.temp=true; return x; }
00757   friend inline Array<T> operator*(const Array<T>& y,const Array<T>& z){ Array<T> x; innerProduct(x,y,z); x.temp=true; return x; }
00759   friend inline Array<T> operator%(const Array<T>& y,const Array<T>& z){ Array<T> x; mult(x,y,z); x.temp=true; return x; }
00761   friend inline Array<T> operator^(const Array<T>& y,const Array<T>& z){ Array<T> x; outerProduct(x,y,z); x.temp=true; return x; }
00763   friend inline Array<T> operator*(const Array<T>& y,T a){ Array<T> x; scalarMultiplication(x,y,a); x.temp=true; return x; }
00765   friend inline Array<T> operator*(T a,const Array<T>& y){ Array<T> x; scalarMultiplication(x,y,a); x.temp=true; return x; }
00766 #endif
00767 };
00768 }
00769 
00770 
00771 //===========================================================================
00772 //
00774 //
00775 
00777 template<class T> inline std::istream& operator>>(std::istream& is,MT::Array<T>& x){ x.read(is);return is; }
00778 
00780 template<class T> inline std::ostream& operator<<(std::ostream& os,const MT::Array<T>& x){
00781   IFEX( if(x.ex) x.assignEx(); )
00782   x.write(os); return os;
00783 }
00784 
00786 template<class T,class S>
00787 inline bool samedim(const MT::Array<T>& a,const MT::Array<S>& b){
00788   return (b.nd==a.nd && b.d0==a.d0 && b.d1==a.d1 && b.d2==a.d2);
00789 }
00790 
00792 template<class T>
00793 inline bool operator==(const MT::Array<T>& v,const MT::Array<T>& w){
00794   if(!samedim(v,w)) return false;
00795   const T *iv,*iw;
00796   for(iv=v.p,iw=w.p; iv!=v.pstop; iv++,iw++)
00797     if (*iv != *iw) return false;
00798   return true;
00799 }
00800 
00802 template<class T>
00803 inline bool operator!=(const MT::Array<T>& v,const MT::Array<T>& w){
00804   return !(v==w);
00805 }
00806 
00808 template<class T>
00809 inline bool operator<(const MT::Array<T>& v,const MT::Array<T>& w){
00810   if(v.N==w.N){
00811     for(uint i=0;i<v.N;i++){
00812       if(v.p[i]>w.p[i]) return false;
00813       if(v.p[i]<w.p[i]) return true;
00814     }
00815     return false; //they are equal
00816   }
00817   return v.N<w.N;
00818 }
00819 
00820 
00821 //===========================================================================
00822 //
00824 //
00825 
00826 #ifndef MT_doxy // exclude these macros when generating the documentation
00827                 // (doxygen can't handle them...)
00828 
00829 //---------- unary operators
00830 
00831 #define UnaryOperator( op )                                               \
00832 template<class T>                                                         \
00833 inline MT::Array<T> operator op (const MT::Array<T>& y){                  \
00834   MT::Array<T> x;                                                         \
00835   x.resize(y);                  \
00836   T *ix=x.p;                                                              \
00837   const T *iy=y.p;                                                        \
00838   for(; iy!=y.pstop; iy++,ix++) *ix= op *iy;            \
00839   x.temp=true; return x;                                                  \
00840 }
00841 
00842 //UnaryOperator( ! )
00843 //UnaryOperator( ~ )
00844 //UnaryOperator( + )
00845 UnaryOperator( - )
00846 #undef UnaryOperator
00847 
00848 
00849 //---------- binary operators
00850 
00851 #ifndef MT_EXPRESSIONS
00852 #define BinaryOperator( op )                                              \
00853 template<class T>                                                         \
00854 inline MT::Array<T> operator op(const MT::Array<T>& v,const MT::Array<T>& w){ \
00855   v.checkEx(); w.checkEx(); \
00856   CHECK(v.N==w.N,               \
00857   "binary operator on different array dimensions ("<<v.N<<", "<<w.N<<")");\
00858   MT::Array<T> z;                                                         \
00859   z.resize(v);                  \
00860   T *iz=z.p;                  \
00861   const T *iw=w.p,*iv=v.p;              \
00862   for(; iv!=v.pstop; iv++,iw++,iz++) *iz = *iv op *iw;                    \
00863   z.temp=true; return z;                                                               \
00864 }                                                                         \
00865                                                                           \
00866 template<class T>                                                         \
00867 inline MT::Array<T> operator op (const MT::Array<T>& v,T w){              \
00868   v.checkEx(); \
00869   MT::Array<T> z;                                                         \
00870   z.resize(v);                  \
00871   T *iz=z.p;                  \
00872   const T *iv=v.p;                \
00873   for(; iv!=v.pstop; iv++,iz++) *iz = *iv op w;                           \
00874   z.temp=true; return z;                                                               \
00875 }                                                                         \
00876                                                                           \
00877 template<class T>                                                         \
00878 inline MT::Array<T> operator op ( T v, const MT::Array<T>& w ){           \
00879   w.checkEx(); \
00880   MT::Array<T> z;                                                         \
00881   z.resize(w);                  \
00882   T *iz=z.p;                  \
00883   const T *iw=w.p;                \
00884   for(; iw!=w.pstop; iw++,iz++) *iz = v op *iw;                           \
00885   z.temp=true; return z;                                                               \
00886 }
00887 
00888 BinaryOperator( || )
00889 BinaryOperator( && )
00890 BinaryOperator( | )
00891   //BinaryOperator( ^ ) I reserve this for the outer product!
00892 BinaryOperator( & )
00893 #ifndef MT_EXPRESSIONS
00894 BinaryOperator( + )
00895 BinaryOperator( - )
00896 #endif
00897   //BinaryOperator( * ) I reserve this operator for the inner product!
00898 BinaryOperator( / )
00899 BinaryOperator( % )
00900 #undef BinaryOperator
00901 #endif
00902 
00903 
00904 //---------- compound assignment operators
00905 
00906 #define CompoundAssignmentOperator( op )                                  \
00907 template<class T,class S>                                                 \
00908 inline MT::Array<T>& operator op (MT::Array<T>& x,const MT::Array<S>& y){ \
00909   y.checkEx(); \
00910   CHECK(x.N==y.N,               \
00911   "binary operator on different array dimensions ("<<x.N<<", "<<y.N<<")");\
00912   T *ix=x.p;                  \
00913   const S *iy=y.p;                \
00914   for(; ix!=x.pstop; ix++,iy++) *ix op *iy;                               \
00915   x.temp=true; return x;                                                               \
00916 }                                                                         \
00917                                                                           \
00918 template<class T>                                                         \
00919 inline MT::Array<T>& operator op ( MT::Array<T>& x, T y ){                \
00920   T *ix=x.p;                  \
00921   for(; ix!=x.pstop; ix++) *ix op y;                                      \
00922   x.temp=true; return x;                                                               \
00923 }                                                                         \
00924 
00925 CompoundAssignmentOperator( |= )
00926 CompoundAssignmentOperator( ^= )
00927 CompoundAssignmentOperator( &= )
00928 CompoundAssignmentOperator( += )
00929 CompoundAssignmentOperator( -= )
00930 CompoundAssignmentOperator( *= )
00931 CompoundAssignmentOperator( /= )
00932 CompoundAssignmentOperator( %= )
00933 #undef CompoundAssignmentOperator
00934 
00935 
00936 //---------- unary functions
00937 
00938 #define UnaryFunction( func )                                             \
00939 template<class T>                                                         \
00940 inline MT::Array<T> func (const MT::Array<T>& v){                         \
00941   v.checkEx(); \
00942   MT::Array<T> z;                                                         \
00943   z.resize(v);                                                            \
00944   T *iz=z.p;                  \
00945   const T *iv=v.p;                \
00946   for(; iv!=v.pstop; iv++, iz++) *iz = ::func( *iv );                     \
00947   z.temp=true; return z;                                                               \
00948 }                   \
00949 
00950 // trigonometric functions
00951 UnaryFunction( acos )
00952 UnaryFunction( asin )
00953 UnaryFunction( atan )
00954 UnaryFunction( cos )
00955 UnaryFunction( sin )
00956 UnaryFunction( tan )
00957 
00958 // hyperbolic functions
00959 UnaryFunction( cosh )
00960 UnaryFunction( sinh )
00961 UnaryFunction( tanh )
00962 UnaryFunction( acosh )
00963 UnaryFunction( asinh )
00964 UnaryFunction( atanh )
00965 
00966 // exponential and logarithmic functions
00967 UnaryFunction( exp )
00968 UnaryFunction( log )
00969 UnaryFunction( log10 )
00970 
00971 // power functions
00972 UnaryFunction( sqrt )
00973 UnaryFunction( cbrt )
00974 
00975 // nearest integer and absolute value
00976 UnaryFunction( ceil )
00977 UnaryFunction( fabs )
00978 UnaryFunction( floor )
00979 #undef UnaryFunction
00980 
00981 
00982 //---------- binary functions
00983 
00984 #define BinaryFunction( func )                                            \
00985 template<class T>                                                         \
00986 inline MT::Array<T> func(const MT::Array<T>& v,const MT::Array<T>& w){    \
00987   v.checkEx(); w.checkEx(); \
00988   CHECK(v.N==w.N,               \
00989   "binary operator on different array dimensions ("<<v.N<<", "<<w.N<<")");\
00990   MT::Array<T> z;                                                         \
00991   z.resize(v);                                                            \
00992   for(uint i=z.N;i--; ) z.p[i]= func(v.p[i],w.p[i]);        \
00993   z.temp=true; return z;                                                               \
00994 }                                                                         \
00995                                                                           \
00996 template<class T>                                                         \
00997 inline MT::Array<T> func(const MT::Array<T>& v,T w){                      \
00998   v.checkEx(); \
00999   MT::Array<T> z;                                                         \
01000   z.resize(v);                                                            \
01001   for(uint i=z.N;i--; ) z.p[i]= func(v.p[i],w);       \
01002   z.temp=true; return z;                                                               \
01003 }                                                                         \
01004                                                                           \
01005 template<class T>                                                         \
01006 inline MT::Array<T> func(T v,const MT::Array<T>& w){                      \
01007   w.checkEx(); \
01008   MT::Array<T> z;                                                         \
01009   z.resize(w);                                                            \
01010   for(uint i=z.N;i--; ) z.p[i]= func(v,w.p[i]);       \
01011   z.temp=true; return z;                                                               \
01012 }
01013 
01014 BinaryFunction( atan2 )
01015 BinaryFunction( pow )
01016 BinaryFunction( fmod )
01017 #undef BinaryFunction
01018 
01019 #endif //(doxygen exclusion)
01020 
01021 
01022 
01023 //===========================================================================
01024 //
01026 //
01027 
01028 #ifdef MT_MSVC
01029 #  ifdef plus
01030 #    undef plus
01031 #  endif
01032 #endif
01033 
01035 template<class T>
01036 void plus(MT::Array<T>& x,const MT::Array<T>& y,const MT::Array<T>& z){
01037   CHECK(y.N==z.N,"must have same size for adding!");
01038   uint i,n=y.N;
01039   x.resize(y);
01040   for(i=0;i<n;i++) x.p[i]=y.p[i]+z.p[i];
01041 }
01042 
01044 template<class T>
01045 void scalarPlus(MT::Array<T>& x,const MT::Array<T>& y,T a){
01046   uint i,n=y.N;
01047   x.resize(y);
01048   for(i=0;i<n;i++) x.p[i]=y.p[i]+a;
01049 }
01050 
01052 template<class T>
01053 void plus(MT::Array<T>& x,T a,const MT::Array<T>& y,T b,const MT::Array<T>& z){
01054   CHECK(y.N==z.N,"must have same size for adding!");
01055   uint i,n=y.N;
01056   x.resize(y);
01057   for(i=0;i<n;i++) x.p[i]=a*y.p[i]+b*z.p[i];
01058 }
01059 
01061 template<class T>
01062 void minus(MT::Array<T>& x,const MT::Array<T>& y,const MT::Array<T>& z){
01063   CHECK(y.N==z.N,"must have same size for subtracting!");
01064   uint i,n=y.N;
01065   x.resize(y);
01066   for(i=0;i<n;i++) x.p[i]=y.p[i]-z.p[i];
01067 }
01068 
01070 template<class T>
01071 void mult(MT::Array<T>& x,const MT::Array<T>& y,const MT::Array<T>& z){
01072   CHECK(y.N==z.N,"must have same size for element-wise multiplication!");
01073   uint i,n=y.N;
01074   x.resize(y);
01075   for(i=0;i<n;i++) x.p[i]=y.p[i]*z.p[i];
01076 }
01077 
01079 template<class T>
01080 void div(MT::Array<T>& x,const MT::Array<T>& y,const MT::Array<T>& z){
01081   CHECK(y.N==z.N,"must have same size for element-wise multiplication!");
01082   uint i,n=y.N;
01083   x.resize(y);
01084   for(i=0;i<n;i++) x.p[i]=y.p[i]/z.p[i];
01085 }
01086 
01087 
01088 //===========================================================================
01089 //
01091 //
01092 
01094 template<class T>
01095 void setIdentity(MT::Array<T>& mat,uint dim){
01096   mat.resize(dim,dim); mat=0.;
01097   for(uint i=0;i<dim;i++) mat(i,i)=1.;
01098 }
01099 
01101 inline MT::Array<double> Identity(uint dim){
01102   MT::Array<double> z;
01103   setIdentity(z,dim);
01104   return z;
01105 }
01106 
01108 inline void makeSymmetric(MT::Array<double>& A){
01109   CHECK(A.nd==2 && A.d0==A.d1,"not symmetric");
01110   uint n=A.d0,i,j;
01111   for(i=1;i<n;i++) for(j=0;j<i;j++) A(j,i) = A(i,j) = .5 * (A(i,j) + A(j,i));
01112 }
01113 
01115 template<class T>
01116 void transpose(MT::Array<T>& x,const MT::Array<T>& y){
01117   CHECK(y.nd<=2,"can only transpose up to 2D arrays");
01118   if(y.nd==2){
01119     uint i,j,d0=y.d1,d1=y.d0;
01120     x.resize(d0,d1);
01121     for(i=0;i<d0;i++) for(j=0;j<d1;j++) x.p[i*d1+j]=y.p[j*d0+i];
01122     return;
01123   }
01124   if(y.nd==1){
01125     x=y;
01126     x.resizeCopy(1,y.N);
01127     return;
01128   }
01129   HALT("transpose not implemented for this dims");
01130 }
01131 
01133 template<class T>
01134 inline void diag(MT::Array<T>& x,const MT::Array<T>& y){
01135   CHECK(y.nd==2 && y.d0==y.d1,"can only give diagonal of symmetric 2D matrix");
01136   x.resize(y.d0);
01137   uint i;
01138   for(i=0;i<x.d0;i++) x(i)=y(i,i);
01139 }
01140 
01142 template<class T>
01143 inline void setDiagonal(MT::Array<T>& x,const MT::Array<T>& v){
01144   CHECK(v.nd==1,"can only give diagonal of 1D array");
01145   x.resize(v.d0,v.d0);
01146   x.setZero();
01147   uint i;
01148   for(i=0;i<v.d0;i++) x(i,i)=v(i);
01149 }
01150 
01152 template<class T>
01153 inline void setDiagonal(MT::Array<T>& x,uint d,double v){
01154   x.resize(d,d);
01155   x.setZero();
01156   uint i;
01157   for(i=0;i<d;i++) x(i,i)=v;
01158 }
01159 
01161 template<class T>
01162 void inverse2d(MT::Array<T>& invA,const MT::Array<T>& A){
01163   invA.resize(2,2);
01164   invA(0,0)=A(1,1); invA(1,1)=A(0,0); invA(0,1)=-A(0,1); invA(1,0)=-A(1,0);
01165   invA/=A(0,0)*A(1,1)-A(0,1)*A(1,0);
01166 }
01167 
01169 template<class T>
01170 void blockMatrix(MT::Array<T>& X,const MT::Array<T>& A,const MT::Array<T>& B,const MT::Array<T>& C,const MT::Array<T>& D){
01171   CHECK(A.nd==2 && B.nd==2 && C.nd==2 && D.nd==2,"");
01172   CHECK(A.d0==B.d0 && A.d1==C.d1 && B.d1==D.d1 && C.d0==D.d0,"");
01173   uint i,j,a=A.d0,b=A.d1;
01174   X.resize(A.d0+C.d0,A.d1+B.d1);
01175   for(i=0;i<A.d0;i++) for(j=0;j<A.d1;j++) X(i  ,j  )=A(i,j);
01176   for(i=0;i<B.d0;i++) for(j=0;j<B.d1;j++) X(i  ,j+b)=B(i,j);
01177   for(i=0;i<C.d0;i++) for(j=0;j<C.d1;j++) X(i+a,j  )=C(i,j);
01178   for(i=0;i<D.d0;i++) for(j=0;j<D.d1;j++) X(i+a,j+b)=D(i,j);
01179 }
01180 
01181 
01182 //===========================================================================
01183 //
01185 //
01186 
01188 template<class T>
01189 void scalarMultiplication(MT::Array<T>& x,const MT::Array<T>& y,const T& a){
01190   x.resize(y);
01191   T *ix=x.p;
01192   const T *iy=y.p;
01193   for(; ix!=x.pstop; ix++,iy++) *ix = *iy * a;
01194 }
01195 
01199 template<class T> 
01200 void innerProduct(MT::Array<T>& x,const MT::Array<T>& y, const MT::Array<T>& z){
01201   /*
01202   if(y.nd==2 && z.nd==2 && y.N==z.N && y.d1==1 && z.d1==1){  //elem-wise
01203     HALT("make element-wise multiplication explicite!");
01204     mult(x,y,z);
01205     return;
01206   }
01207   */
01208   if(y.nd==2 && z.nd==2){ //plain matrix multiplication
01209     CHECK(y.d1==z.d0,"wrong dimensions for inner product");
01210     uint i,j,d0=y.d0,d1=z.d1,dk=y.d1;
01211     T *a,*astop,*b,*c;
01212     x.resize(d0,d1); x.setZero();
01213     c=x.p;
01214     for(i=0;i<d0;i++) for(j=0;j<d1;j++){
01215       //for(s=0.,k=0;k<dk;k++) s+=y.p[i*dk+k]*z.p[k*d1+j];
01216       //this is faster:
01217       a=y.p+i*dk; astop=a+dk; b=z.p+j;
01218       for(;a<astop; a++, b+=d1) (*c)+=(*a) * (*b);
01219       c++;
01220     }
01221     return;
01222   }
01223   if(y.nd==2 && z.nd==1){ //matrix x vector -> vector
01224     CHECK(y.d1==z.d0,"wrong dimensions for inner product");
01225     uint i,k,d0=y.d0,dk=y.d1;
01226     T *a,*b;
01227     T s;
01228     x.resize(d0);
01229     for(i=0;i<d0;i++){
01230       //for(s=0.,k=0;k<dk;k++) s+=y.p[i*dk+k]*z.p[k];
01231       //this is faster:
01232       a=y.p+i*dk; b=z.p;
01233       for(s=0,k=0;k<dk;k++){ s+=(*a) * (*b); a++; b++; }
01234       x.p[i]=s;
01235     }
01236     return;
01237   }
01238   if(y.nd==1 && z.nd==2 && z.d0==1){ //vector x vector^T -> matrix (outer product)
01239     CHECK(y.d0==z.d1,"wrong dimensions for inner product");
01240     uint i,j,d0=y.d0,d1=z.d1;
01241     x.resize(d0,d1);
01242     for(i=0;i<d0;i++) for(j=0;j<d1;j++) x(i,j)=y(i)*z(0,j);
01243     return;
01244   }
01245   if(y.nd==1 && z.nd==2){ //vector^T x matrix -> vector^T
01246     CHECK(y.d0==z.d0,"wrong dimensions for inner product");
01247     uint i,k,d0=z.d1,dk=y.d0;
01248     x.resize(d0);
01249     T s;
01250     for(i=0;i<d0;i++){
01251       for(s=0,k=0;k<dk;k++) s+=y.p[k]*z.p[k*d0+i];
01252       x.p[i]=s;
01253     }
01254     return;
01255   }
01256   if(y.nd==1 && z.nd==1 && z.N==1){ //vector multiplied with scalar (disguised as 1D vector)
01257     uint k,dk=y.N;
01258     x.resize(y.N);
01259     for(k=0;k<dk;k++) x.p[k]=y.p[k]*z.p[0];
01260     return;
01261   }
01262   if(y.nd==1 && z.nd==1){ //should be scalar product, but be careful
01263     HALT("what do you want? scalar product or element wise multiplication?");
01264     CHECK(y.d0==z.d0,"wrong dimensions for inner product");
01265     uint k,dk=y.d0;
01266     x.resize(1);
01267     T s;
01268     for(s=0,k=0;k<dk;k++) s+=y.p[k]*z.p[k];
01269     x.p[0]=s;
01270     return;
01271   }
01272   HALT("inner product - not yet implemented for these dimensions");
01273 }
01274 
01277 template<class T> 
01278 void outerProduct(MT::Array<T>& x,const MT::Array<T>& y, const MT::Array<T>& z){
01279   if(y.nd==1 && z.nd==1){
01280     uint i,j,d0=y.d0,d1=z.d0;
01281     x.resize(d0,d1);
01282     for(i=0;i<d0;i++) for(j=0;j<d1;j++) x.p[i*d1+j]=y.p[i]*z.p[j];
01283     return;
01284   }
01285   HALT("outer product - not yet implemented for these dimensions");
01286 }
01287 
01289 template<class T> 
01290 inline T scalarProduct(const MT::Array<T>& v, const MT::Array<T>& w){
01291   CHECK(v.N==w.N, 
01292     "scalar product on different array dimensions ("<<v.N<<", "<<w.N<<")");
01293   T t(0);
01294   for(uint i=v.N; i--; t+=v.p[i]*w.p[i]);
01295   return t;
01296 }
01297 
01299 template<class T> 
01300 inline T scalarProduct(const MT::Array<T>& g,const MT::Array<T>& v, const MT::Array<T>& w){
01301   CHECK(v.N==w.N && g.nd==2 && g.d0==v.N && g.d1==w.N,
01302     "scalar product on different array dimensions ("<<v.N<<", "<<w.N<<")");
01303   T t(0);
01304   for(uint i=0;i<g.d0;i++) for(uint j=0;j<g.d1;j++) t+=g(i,j)*v.p[i]*w.p[j];
01305   return t;
01306 }
01307 
01308 
01309 //===========================================================================
01310 //
01312 //
01313 
01315 template<class T> 
01316 inline T sqrDistance(const MT::Array<T>& v, const MT::Array<T>& w){
01317   CHECK(v.N==w.N, 
01318     "sqrDistance on different array dimensions ("<<v.N<<", "<<w.N<<")");
01319   T d,t(0);
01320   for(uint i=v.N;i--;){ d=v.p[i]-w.p[i]; t+=d*d; }
01321   return t;
01322 }
01323 
01325 template<class T> 
01326 inline T sqrDistance(const MT::Array<T>& v, const MT::Array<T>& w,const MT::Array<bool>& mask){
01327   CHECK(v.N==w.N, 
01328     "sqrDistance on different array dimensions ("<<v.N<<", "<<w.N<<")");
01329   T d,t(0);
01330   for(uint i=v.N;i--;) if(mask(i)){ d=v.p[i]-w.p[i]; t+=d*d; }
01331   return t;
01332 }
01333 
01335 template<class T>
01336 inline T sqrDistance(const MT::Array<T>& g,const MT::Array<T>& v,const MT::Array<T>& w){
01337   MT::Array<T> d;
01338   minus(d,v,w);
01339   return scalarProduct(g,d,d);
01340 }
01341 
01343 template<class T> 
01344 inline T euclideanDistance(const MT::Array<T>& v, const MT::Array<T>& w){
01345   return ::sqrt(sqrDistance(v,w));
01346 }
01347 
01348 
01349 
01350 //===========================================================================
01351 //
01353 //
01354 
01356 template<class T>
01357 inline T sum(const MT::Array<T>& v){
01358   T t(0);
01359   for(uint i=v.N; i--; t+=v.p[i]);
01360   return t;
01361 }
01362 
01364 template<class T>
01365 inline T sumOfAbs(const MT::Array<T>& v){
01366   T t(0);
01367   for(uint i=v.N; i--; t+=fabs(v.p[i]));
01368   return t;
01369 }
01370 
01372 template<class T>
01373 inline T sumOfSqr(const MT::Array<T>& v){
01374   T t(0);
01375   for(uint i=v.N; i--; t+=v.p[i]*v.p[i]);
01376   return t;
01377 }
01378 
01380 template<class T>
01381 inline T norm(const MT::Array<T>& v){ return sqrt(sumOfSqr(v)); }
01382 
01384 template<class T>
01385 inline T trace(const MT::Array<T>& v){
01386   CHECK(v.nd==2 && v.d0==v.d1,"only for squared matrix");
01387   T t(0);
01388   for(uint i=0;i<v.d0;i++) t+=v(i,i);
01389   return t;
01390 }
01391 
01392 
01393 //===========================================================================
01394 //
01396 //
01397 
01400 template<class T>
01401 inline T entropy(const MT::Array<T>& v){
01402   T t(0);
01403   for(uint i=v.N; i--;) if(v.p[i]) t-=v.p[i]*::log(v.p[i]);
01404   return t/MT_LN2;
01405 }
01406 
01408 template<class T>
01409 inline T normalizeDist(MT::Array<T>& v){
01410   double Z=sum(v);
01411   if(Z) v/=Z; else v=1./(double)v.N;
01412   return Z;
01413 }
01414 
01415 template<class T>
01416 inline void eliminate(MT::Array<T>& x,const MT::Array<T>& y,uint d){
01417   CHECK(y.nd==2,"only implemented for 2D yet");
01418   uint i,j;
01419   if(d==1){
01420     x.resize(y.d0); x=0.;
01421     for(i=0;i<y.d0;i++) for(j=0;j<y.d1;j++) x(i)+=y(i,j);
01422   }
01423   if(d==0){
01424     x.resize(y.d1); x=0.;
01425     for(i=0;i<y.d0;i++) for(j=0;j<y.d1;j++) x(j)+=y(i,j);
01426   }
01427 }
01428 
01429 template<class T>
01430 inline void eliminate(MT::Array<T>& x,const MT::Array<T>& y,uint d,uint e){
01431   CHECK(y.nd==3,"only implemented for 3D yet");
01432   uint i,j,k;
01433   if(d==1 && e==2){
01434     x.resize(y.d0); x=0.;
01435     for(i=0;i<y.d0;i++) for(j=0;j<y.d1;j++) for(k=0;k<y.d2;k++) x(i)+=y(i,j,k);
01436   }
01437   if(d==0 && e==2){
01438     x.resize(y.d1); x=0.;
01439     for(i=0;i<y.d0;i++) for(j=0;j<y.d1;j++) for(k=0;k<y.d2;k++) x(j)+=y(i,j,k);
01440   }
01441   if(d==0 && e==1){
01442     x.resize(y.d2); x=0.;
01443     for(i=0;i<y.d0;i++) for(j=0;j<y.d1;j++) for(k=0;k<y.d2;k++) x(k)+=y(i,j,k);
01444   }
01445 }
01446 
01447 
01448 //===========================================================================
01449 //
01451 //
01452 
01454 template<class T>
01455 inline T product(const MT::Array<T>& v){
01456   T t(1);
01457   for(uint i=v.N; i--; t *= v.p[i]);
01458   return t;
01459 }
01460 
01461 
01462 //===========================================================================
01463 //
01465 //
01466 
01468 template<class T>
01469 inline T minA(const MT::Array<T>& v){
01470   CHECK(v.N>0,"");
01471   T t(v.p[0]);
01472   for(uint i=1; i<v.N; ++i) if(v.p[i]<t) t=v.p[i];
01473   return t;
01474 }
01475 
01477 template<class T>
01478 inline T minA(const MT::Array<T>& v, uint & ind, uint start=0, uint end=0){
01479   CHECK(v.N>0,"");
01480   CHECK(v.N>start,"");  
01481   CHECK(v.N>=end,"");
01482   CHECK(end>=start,"");
01483   T t(v(start));
01484   ind=start;
01485   if (end==0) end=v.N;
01486   for(uint i=start; i<end; ++i) if(v.p[i]<t){
01487     t  =v.p[i];
01488     ind=i;
01489   }
01490   return t;
01491 }
01492 
01494 template<class T>
01495 inline T maxA(const MT::Array<T>& v){
01496   CHECK(v.N>0,"");
01497   T t(v.p[0]);
01498   for(uint i=1; i<v.N; ++i) if(v.p[i]>t) t=v.p[i];
01499   return t;
01500 }
01501 
01503 template<class T>
01504 inline T maxA(const MT::Array<T>& v, uint & ind, uint start=0, uint end=0){
01505   CHECK(v.N>0,"");
01506   CHECK(v.N>start,"");
01507   CHECK(v.N>=end,"");
01508   CHECK(end>=start,"");
01509   T t(v(start));
01510   ind=start;
01511   if(end==0){
01512     end=v.N;
01513   }
01514   for(uint i=start; i<end; ++i) if(v.p[i]>t){
01515     t  =v.p[i];
01516     ind=long(i);
01517   }
01518   return t;
01519 }
01520 
01522 template<class T>
01523 inline void minmaxA(const MT::Array<T>& v, T& minVal, T& maxVal){
01524   if(v.N>0){
01525     minVal=maxVal=v.p[0];
01526     for(uint i=1; i<v.N; ++i){
01527       if(v.p[i]<minVal) minVal=v.p[i];
01528       else if(v.p[i]>maxVal) maxVal=v.p[i];
01529     }
01530   }
01531 }
01532 
01533 template<class T>
01534 inline T absMax(const MT::Array<T>& x){
01535   CHECK(x.N>0,"");
01536   T t(x.p[0]);
01537   for(uint i=1; i<x.N; ++i) if(fabs(x.p[i])>t) t=fabs(x.p[i]);
01538   return t;
01539 }
01540 
01541 
01542 //===========================================================================
01543 //
01545 //
01546 
01548 template<class T>
01549 void rndInt(MT::Array<T>& a,int low=0.,int high=1.){
01550   for(uint i=0;i<a.N;i++) a.p[i]=low+(int)rnd.num(1+high-low);
01551 }
01552 
01554 template<class T>
01555 void rndUni(MT::Array<T>& a,double low=0.,double high=1.){
01556   for(uint i=0;i<a.N;i++) a.p[i]=rnd.uni(low,high);
01557 }
01558 
01560 template<class T>
01561 void rndGauss(MT::Array<T>& a,double stdDev,bool add=false){
01562   if(a.nd==2 && a.d0!=a.d1)  //assumes that it is a batch of data
01563     stdDev/=::sqrt((double)a.d1);
01564   else
01565     stdDev/=::sqrt((double)a.N);
01566   if(!add) for(uint i=0;i<a.N;i++) a.p[i]=rnd.gauss(stdDev);
01567   else     for(uint i=0;i<a.N;i++) a.p[i]+=rnd.gauss(stdDev);
01568 }
01569 
01571 template<class T>
01572 MT::Array<T>& rndGauss(double stdDev,uint dim){
01573   static MT::Array<T> z;
01574   stdDev/=::sqrt(dim);
01575   z.resize(dim);
01576   rndGauss(z,stdDev);
01577   return z;
01578 }
01579 
01583 template<class T>
01584 uint softMax(const MT::Array<T>& a,MT::Array<double>& soft,double beta){
01585   double norm=0.,r;
01586   uint i; int sel=-1;
01587   soft.resize(a);
01588   for(i=0;i<a.N;i++){
01589     soft(i)=exp(beta*a(i));
01590     norm+=soft(i);
01591   }
01592   r=rnd.uni();
01593   for(i=0;i<a.N;i++){
01594     soft(i)/=norm;
01595     r-=soft(i);
01596     if(sel==-1 && r<0.) sel=i;
01597   }
01598   return sel;
01599 }
01600 
01601 
01602 //===========================================================================
01603 //
01605 //
01606 
01607 template<class T>
01608 MT::Array<T>& grid(MT::Array<T>& a,uint dim,T lo,T hi,uint steps){
01609   uint i,j;
01610   if(dim==1){
01611     a.resize(steps+1,1);
01612     for(i=0;i<a.d0;i++) a(i,0)=lo+(hi-lo)*i/steps;
01613     return a;
01614   }
01615   if(dim==2){
01616     a.resize(steps+1,steps+1,2);
01617     for(i=0;i<a.d0;i++) for(j=0;j<a.d1;j++){
01618       a(i,j,0)=lo+(hi-lo)*j/steps;
01619       a(i,j,1)=lo+(hi-lo)*i/steps;
01620     }
01621     a.reshape(a.d0*a.d1,2);
01622     return a;
01623   }
01624   HALT("not implemented yet");
01625 }
01626 
01627 
01628 //===========================================================================
01629 //
01631 //
01632 
01633 #ifdef MT_EXPRESSIONS
01634 
01635 inline MT::Array<double> operator*(double a,const MT::Array<double>& b){ return MT::Array<double>(new MT::Ex((void*)&b,0,a,false)); }
01636 inline MT::Array<double> operator*(const MT::Array<double>& a,double b){ return MT::Array<double>(new MT::Ex((void*)&a,0,b,false)); }
01637 inline MT::Array<double> operator/(const MT::Array<double>& a,double b){ return MT::Array<double>(new MT::Ex((void*)&a,0,1./b,false)); }
01638 inline MT::Array<double> operator+(double a,const MT::Array<double>& b){ return MT::Array<double>(new MT::Ex((void*)&b,a,1,false)); }
01639 inline MT::Array<double> operator+(const MT::Array<double>& a,double b){ return MT::Array<double>(new MT::Ex((void*)&a,b,1,false)); }
01640 inline MT::Array<double> operator-(double a,const MT::Array<double>& b){ return MT::Array<double>(new MT::Ex((void*)&b,-a,1,false)); }
01641 inline MT::Array<double> operator-(const MT::Array<double>& a,double b){ return MT::Array<double>(new MT::Ex((void*)&a,-b,1,false)); }
01642 inline MT::Array<double> operator~(const MT::Array<double>& a){ return MT::Array<double>(new MT::Ex((void*)&a,0,1,true)); }
01643 
01644 inline MT::Array<double> operator*(const MT::Array<double>& y,const MT::Array<double>& z){ return MT::Array<double>(new MT::Ex(MT::PROD,(void*)&y,(void*)&z)); }
01645 inline MT::Array<double> operator^(const MT::Array<double>& y,const MT::Array<double>& z){ return MT::Array<double>(new MT::Ex(MT::OUT ,(void*)&y,(void*)&z)); }
01646 inline MT::Array<double> operator%(const MT::Array<double>& y,const MT::Array<double>& z){ return MT::Array<double>(new MT::Ex(MT::MUL ,(void*)&y,(void*)&z)); }
01647 
01648 inline MT::Array<double> operator/(const MT::Array<double>& y,const MT::Array<double>& z){ return MT::Array<double>(new MT::Ex(MT::Div,(void*)&y,(void*)&z)); }
01649 inline MT::Array<double> operator+(const MT::Array<double>& y,const MT::Array<double>& z){ return MT::Array<double>(new MT::Ex(MT::PLUS,(void*)&y,(void*)&z)); }
01650 inline MT::Array<double> operator-(const MT::Array<double>& y,const MT::Array<double>& z){ return MT::Array<double>(new MT::Ex(MT::MINUS,(void*)&y,(void*)&z)); }
01651 
01652 void assign(MT::Array<double>& x,const MT::Array<double>& a);
01653 void assign(MT::Array<double>& x);
01654 template<class T> inline void MT::Array<T>::assignEx(const MT::Array<T>& a) const{ ::assign(*((MT::Array<double>*)this),*((MT::Array<double>*)&a)); }
01655 template<class T> inline void MT::Array<T>::assignEx() const{ ::assign(*((MT::Array<double>*)this)); }
01656 
01657 #endif //MT_EXPRESSIONS
01658 
01659 
01660 //===========================================================================
01661 //
01663 //
01664 
01666 void gnuplot(const MT::Array<double>& X);
01667 
01668 //void write(const MT::Array<double>& X,MT::Array<double>& Y,const char* name);
01669 //void write(const MT::Array<double>& X,const char* name);
01670 
01671 //===========================================================================
01672 //
01674 //
01675 
01676 namespace MT{
01678 class Permutation:public Array<uint>{
01679 private:
01680   Array<uint> storage;
01681 
01682 public:
01684   Permutation(){ init(0); };
01685 
01686 public:
01687 
01688   void init(uint n){ resize(n); for(uint i=0;i<N;i++) elem(i)=i; }
01690   void random(uint n){ init(n); random(); }
01692   void reverse(uint n){ resize(n); for(uint i=0;i<N;i++) elem(N-1-i)=i; }
01693 
01694 public:
01695 
01696   void push(int offset=1){for(uint i=0;i<N;i++) elem(i)=((int)elem(i)+offset)%N;}
01698   void store(){ storage.operator=(*this); }
01700   void restore(){ (Array<uint>&)(*this)=storage; }
01702   void permute(uint i,uint j){ uint x=elem(i); elem(i)=elem(j); elem(j)=x; }
01704   void random(){
01705     int j,r;
01706     for(j=N-1;j>=1;j--){
01707       r=rnd(j+1);
01708       permute(r,j);
01709     }
01710   }
01712   template<class T>
01713   void permute(Array<T>& a){
01714     CHECK(N<=a.N,"array smaller than permutation ("<<a.N<<"<"<<N<<")");
01715     Array<T> b=a;
01716     for(uint i=0;i<N;i++) a(i)=b(p[i]);
01717   }
01719   template<class T>
01720   void invpermute(Array<T>& a){
01721     CHECK(N<=a.N,"array smaller than permutation ("<<a.N<<"<"<<N<<")");
01722     Array<T> b=a;
01723     for(uint i=0;i<N;i++) a(p[i])=b(i);
01724   }
01725 };
01726 }
01727 
01728 
01729 //===========================================================================
01730 //
01732 //
01733 
01734 typedef MT::Array<double> doubleA;
01735 typedef MT::Array<uint> uintA;
01736 typedef MT::Array<int> intA;
01737 typedef MT::Array<char> charA;
01738 typedef MT::Array<byte> byteA;
01739 typedef MT::Array<bool> boolA;
01740 
01741 //using MT::Array;
01742 
01743 
01744 //===========================================================================
01745 //
01746 // implementations
01747 //
01748 
01749 template<class T> char MT::Array<T>::memMoveInit=-1;
01750 template<class T> int MT::Array<T>::sizeT=-1;
01751 
01752 #ifdef MT_IMPLEMENTATION
01753 #  include"array.cpp"
01754 #endif
01755 
01756 
01757 #endif
[]