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 //
00045 // Array class
00046 //
00047 
00048 namespace MT{
00058 template<class T>
00059 class Array{
00060 public:
00061   T *p;     
00062   T *pstop; 
00063   uint N;   
00064   uint nd;  
00065   uint d0;  
00066   uint d1;  
00067   uint d2;  
00068   uint M;
00069   bool reference;
00070   T **pp;   
00071   MT::Array<uint> *sparse; 
00072   IFEX( Ex *ex; )
00073   bool temp;
00074 
00075 public:
00076 
00081   bool memMove;
00082 
00088   bool flexiMem;
00090 private:
00091   static char memMoveInit;
00092   static int sizeT;
00093   void init(){
00094     reference=false;
00095     memMove=false;
00096     if(sizeT==-1) sizeT=sizeof(T);
00097     if(memMoveInit==-1){
00098       memMoveInit=0;
00099       if(typeid(T)==typeid(bool) ||
00100   typeid(T)==typeid(char) ||
00101   typeid(T)==typeid(unsigned char) ||
00102   typeid(T)==typeid(int) ||
00103   typeid(T)==typeid(unsigned int) ||
00104   typeid(T)==typeid(short) ||
00105   typeid(T)==typeid(unsigned short) ||
00106   typeid(T)==typeid(long) ||
00107   typeid(T)==typeid(unsigned long) ||
00108   typeid(T)==typeid(float) ||
00109   typeid(T)==typeid(double)) memMoveInit=1;
00110     }
00111     memMove=(memMoveInit==1);
00112     flexiMem=true;
00113     p=pstop=0;
00114     M=N=nd=d0=d1=d2=0;
00115     pp=0;
00116     sparse=0;
00117     IFEX( ex=0; );
00118     temp=false;
00119   }
00120 
00121 public:
00122 
00124   Array(){ init(); }
00125 
00127   Array(const Array<T>& a){ init(); operator=(a); }
00128 
00130   Array(uint i){ init(); resize(i); }
00131 
00133   Array(uint i,uint j){ init(); resize(i,j); }
00134 
00136   Array(uint i,uint j,uint k){ init(); resize(i,j,k); }
00137 
00139   Array(const Array<T>& a,uint dim){ init(); referToSubDim(a,dim); }
00140 
00142   Array(const Array<T>& a,uint i,uint j){ init(); referToSubDim(a,i,j); }
00143 
00145   Array(T* p,uint size){ init(); referTo(p,size); }
00146 
00148   Array(const char* str){ init(); set(str); }
00149 
00150   IFEX( Array(Ex *_ex){ ex=_ex; } )
00151 
00152   ~Array(){
00153     IFEX( if(ex){ delete ex; return; } );
00154     freeMEM();
00155   }
00156   
00157 
00158 public:
00159 
00161   void clear(){ freeMEM(); }
00162 
00164   void resize(uint D0){ nd=1; d0=D0; resizeMEM(d0,false); }
00166   void resizeCopy(uint D0){ nd=1; d0=D0; resizeMEM(d0,true); }
00168   void resizeZero(uint D0){ nd=1; d0=D0; resizeMEM(d0,false); setZero(); }
00170   void reshape(uint D0){ CHECK(N==D0,"reshape must preserve total memory size"); nd=1; d0=D0; d1=d2=0; }
00171 
00172   // 2D resize
00173   void resize(uint D0,uint D1){ nd=2; d0=D0; d1=D1; resizeMEM(d0*d1,false); }
00174   void resizeCopy(uint D0,uint D1){ nd=2; d0=D0; d1=D1; resizeMEM(d0*d1,true); }
00175   void reshape(uint D0,uint D1){ CHECK(N==D0*D1,"reshape must preserve total memory size"); nd=2; d0=D0; d1=D1; d2=0; }
00176 
00177   // 3D resize
00178   void resize(uint D0,uint D1,uint D2){ nd=3; d0=D0; d1=D1; d2=D2; resizeMEM(d0*d1*d2,false); }
00179   void resizeCopy(uint D0,uint D1,uint D2){ nd=3; d0=D0; d1=D1; d2=D2; resizeMEM(d0*d1*d2,true); }
00180   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; }
00181 
00182 #ifndef MT_MSVC
00183 
00184   void resize(const Array<T>& a){ nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2; resizeMEM(a.N,false); }
00185   void resizeCopy(const Array<T>& a){ nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2; resizeMEM(a.N,true); }
00186 #endif
00187 
00189   template<class S>
00190   void resize(const Array<S>& a){ nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2; resizeMEM(a.N,false); }
00191   template<class S>
00192   void resizeCopy(const Array<S>& a){ nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2; resizeMEM(a.N,true); }
00193 
00194   void makeSparse(){
00195     CHECK(!sparse,"only once yet");
00196     uint n=0;
00197     if(nd==1){
00198       uint i;
00199       sparse=new MT::Array<uint> [2];
00200       sparse[1].resize(d0); sparse[1]=-1;
00201       for(i=0;i<d0;i++) if(p[i]){
00202   sparse[0].append(i); //list of entries (maps n->i)
00203   sparse[1](i)=n;      //index list to entries (maps i->n)
00204   permute(i,n);
00205   n++;
00206       }
00207       N=n; resizeMEM(n,true);
00208       return;
00209     }
00210     if(nd==2){
00211       uint i,j;
00212       MT::Array<uint> pair(2);
00213       sparse=new MT::Array<uint> [1+d1+d0];
00214       for(i=0;i<d0;i++) for(j=0;j<d1;j++) if(p[i*d1+j]){
00215   pair(0)=i; pair(1)=j; sparse[0].append(pair);   sparse[0].reshape(n+1,2);
00216   permute(i*d1+j,n);
00217   //register entry in columns an row indices
00218   pair(0)=i; pair(1)=n; sparse[1+j]   .append(pair); sparse[1+j]   .reshape(sparse[1+j]   .N/2,2);
00219   pair(0)=j; pair(1)=n; sparse[1+d1+i].append(pair); sparse[1+d1+j].reshape(sparse[1+d1+j].N/2,2);
00220   n++;
00221       }
00222       N=n; resizeMEM(n,true);
00223       return;
00224     }
00225   }
00226 
00228   uint memsize(){ return M*sizeof(T); }
00229 
00230 
00231 
00232 public:
00233 
00234   void resizeMEM(uint n,bool copy){
00235     if(n==N) return;
00236     CHECK(!reference,"real resize of subarray is not allowed! (only a resize without changing memory size)");
00237     uint i;
00238     T *pold=p;
00239     if(!flexiMem){
00240       if(M!=n){
00241   if(n) p=new T [n];
00242   if(!p) HALT("MT::Array failed memory allocation of "<<n*sizeT<<"bytes");
00243   if(copy && !memMove) for(i=N<n?N:n;i--;) p[i]=pold[i];
00244   if(copy && memMove) memmove(p,pold,sizeT*(N<n?N:n));
00245   if(M) delete[] pold;
00246   N=M=n;
00247   pstop=p+N;
00248       }
00249     }else{
00250       if(n>0 && M==0){ //first time
00251   p=new T [n];
00252   if(!p) HALT("MT::Array failed memory allocation of "<<n*sizeT<<"bytes");
00253   M=n;
00254       }else if(n>M || 10+2*n<M/2){
00255   uint oversize=10+2*n;
00256   p=new T [oversize];
00257   if(!p) HALT("MT::Array failed memory allocation of "<<n*sizeT<<"bytes");
00258   if(copy && !memMove) for(i=N<n?N:n;i--;) p[i]=pold[i];
00259   if(copy && memMove) memmove(p,pold,sizeT*(N<n?N:n));
00260   if(M) delete[] pold;
00261   M=oversize;
00262       }
00263       N=n;
00264       pstop=p+N;
00265     }
00266     if(pp){ delete[] pp; pp=0; }
00267   }
00269   void freeMEM(){
00270     if(M) delete[] p;
00271     if(pp) delete[] pp;
00272     if(sparse) delete[] sparse;
00273     p=pstop=0;
00274     M=N=nd=d0=d1=d2=0;
00275     sparse=0;
00276   }
00277 
00278 
00279 public:
00280 
00282   void append(const T& x){
00283     if(nd==2 && d1==1)
00284       resizeCopy(d0+1,d1);
00285     else
00286       resizeCopy(N+1);
00287     p[N-1]=x;
00288   }
00290   void append(const Array<T>& x){
00291     uint oldN=N,i;
00292     if(nd==2 && d1==x.N) 
00293       resizeCopy(d0+1,d1);
00294     else
00295       resizeCopy(N+x.N);
00296     if(!memMove){
00297       for(i=0;i<x.N;i++) p[oldN+i]=x.p[i];
00298     }else memmove(p+oldN,x.p,sizeT*x.N);
00299   }
00301   void append(const T *p,uint n){
00302     uint oldN=N,i;
00303     if(nd==2 && d1==n) 
00304       resizeCopy(d0+1,d1);
00305     else
00306       resizeCopy(N+n);
00307     if(!memMove){
00308       for(i=0;i<n;i++) p[n+i]=p[i];
00309     }else memmove(p+oldN,p,sizeT*n);
00310   }
00311 
00312   void replicate(uint copies){
00313     if(copies<2) return;
00314     uint i,oldN=N;
00315     resizeCopy(copies*N);
00316     if(!memMove){
00317       HALT("not implemented yet");
00318     }else{
00319       for(i=0;i<copies;i++) memmove(p+i*oldN,p,sizeT*oldN);
00320     }
00321   }
00322 
00324   void insert(uint i,const T& x){
00325     CHECK(memMove,"only with memMove");
00326     uint Nold=N;
00327     resizeCopy(Nold+1);
00328     if(i<Nold) memmove(p+i+1,p+i,sizeT*(Nold-i));
00329     p[i]=x;
00330   }
00331 
00333   void remove(uint i,uint n=1){
00334     CHECK(memMove,"only with memMove");
00335     if(N>i+n) memmove(p+i,p+i+n,sizeT*(N-i-n));
00336     resizeCopy(N-n);
00337   }
00338 
00340   void removePerm(uint i){
00341     p[i]=p[N-1];
00342     resizeCopy(N-1);
00343   }
00344 
00346   void replace(uint i,uint n,const Array<T>& x){
00347     CHECK(memMove,"only with memMove");
00348     uint Nold=N;
00349     if(n==x.N){
00350       memmove(p+i,x.p,sizeT*(x.N));
00351     }else if(n>x.N){
00352       memmove(p+i+x.N,p+i+n,sizeT*(Nold-i-n));
00353       if(i+n<Nold) memmove(p+i,x.p,sizeT*(x.N));
00354       resizeCopy(Nold-n+x.N);
00355     }else{
00356       resizeCopy(Nold+x.N-n);
00357       if(i+n<Nold) memmove(p+i+x.N,p+i+n,sizeT*(Nold-i-n));
00358       memmove(p+i,x.p,sizeT*(x.N));
00359     }
00360   }
00361 
00363   void delRow(uint i){
00364     CHECK(memMove,"only with memMove");
00365     CHECK(nd==2,"only for matricies");
00366     CHECK(i<d0,"range check error");
00367     uint n=d1;
00368     if(i+1<d0) memmove(p+i*n,p+(i+1)*n,sizeT*(d0-i-1)*n);
00369     resizeCopy(d0-1,n);
00370   }
00371 
00373   void delColumn(uint i){
00374     CHECK(memMove,"only with memMove");
00375     CHECK(nd==2,"only for matricies");
00376     CHECK(i<d1,"range check error");
00377     uint n=d1;
00378     for(uint j=0;j<d0;j++){
00379       memmove(p+j*(n-1)  ,p+j*n  ,sizeT*i);
00380       memmove(p+j*(n-1)+i,p+j*n+(i+1),sizeT*(n-i-1));
00381     }
00382     resizeCopy(d0,n-1);
00383   }
00384   
00385 
00386 
00387 public:
00388 
00389   T& elem(uint i) const{
00390     CHECK(i<N,"range error ("<<i<<">="<<N<<")");
00391     return p[i];
00392   }
00393 
00394   /* scalar reference (legal iff N==1) */
00395   /*operator T&() const{ 
00396     CHECK(N==1,"scalar reference ("<<N<<"!=1)");
00397     return *e; }*/
00398 
00400   T& operator()(uint i) const{ 
00401     CHECK(nd==1 && i<d0 && i<N,
00402     "1D range error ("<<nd<<"=1,"<<i<<"<"<<d0<<")");
00403     return p[i];
00404   }
00405 
00407   T& operator()(uint i,uint j) const{
00408     CHECK(nd==2 && i<d0 && j<d1 && !sparse,
00409     "2D range error ("<<nd<<"=2,"<<i<<"<"<<d0<<","<<j<<"<"<<d1<<")");
00410     return p[i*d1+j];
00411   }
00412 
00414   T& operator()(uint i,uint j,uint k) const{
00415     CHECK(nd==3 && i<d0 && j<d1 && k<d2 && !sparse,
00416     "3D range error ("<<nd<<"=3,"<<i<<"<"<<d0<<","<<j<<"<"<<d1<<","<<k<<"<"<<d2<<")");
00417     return p[(i*d1+j)*d2+k];
00418   }
00419 
00421   Array<T> operator[](uint i) const{ return Array(*this,i); }
00422 
00424   Array<T>& operator()(){ return (*this); }
00425 
00428   uint maxIndex(){ uint i,m=0; for(i=0;i<N;i++) if(p[i]>p[m]) m=i; return m; }
00429 
00432   void maxIndex(uint& i,uint& j){ CHECK(nd==2,"needs 2D array"); j=maxIndex(); i=j/d1; j=j%d1; }
00433 
00436   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; }
00437 
00439   uint minIndex(){ uint i,m=0; for(i=0;i<N;i++) if(p[i]<p[m]) m=i; return m; }
00440 
00441   bool contains(T& x){ uint i; for(i=0;i<N;i++) if(p[i]==x) return true; return false; }
00442 
00446   Array<T> sub(uint i,int I) const{
00447     CHECK(nd==1,"1D range error ");
00448     Array<T> x;
00449     if(I==-1) I=d0-1;
00450     CHECK(i<=I,"lower limit higher than upper!");
00451     x.resize(I-i+1);
00452     uint k;
00453     for(k=i;k<=I;k++) x(k-i)=operator()(k);
00454     return x;
00455   }
00456 
00460   Array<T> sub(uint i,int I,uint j,int J) const{
00461     CHECK(nd==2,"2D range error ");
00462     Array<T> x;
00463     if(I==-1) I=d0-1;
00464     if(J==-1) J=d1-1;
00465     CHECK(i<(uint)I && j<(uint)J,"lower limit must be higher than upper!");
00466     x.resize(I-i+1,J-j+1);
00467     uint k,l;
00468     for(k=i;k<=(uint)I;k++) for(l=j;l<=(uint)J;l++) x(k-i,l-j)=operator()(k,l);
00469     return x;
00470   }
00471 
00476   Array<T> sub(uint i,int I,Array<uint> cols) const{
00477     CHECK(nd==2,"2D range error ");
00478     Array<T> x;
00479     if(I==-1) I=d0-1;
00480     CHECK(i<=(uint)I,"lower limit higher than upper!");
00481     x.resize(I-i+1,cols.N);
00482     uint k,l;
00483     for(k=i;k<=(uint)I;k++) for(l=0;l<cols.N;l++) x(k-i,l)=operator()(k,cols(l));
00484     return x;
00485   }
00486   
00487 public:
00488 
00490   T** getCarray(){
00491     CHECK(nd==2,"only 2D array gives C-array of type T**");
00492     if(pp) return pp;
00493     pp=new T* [d0];
00494     for(uint i=0;i<d0;i++) pp[i]=p+i*d1;
00495     return pp;
00496   }
00497 
00500   T** getPointers(Array<T*>& array2d) const{
00501     CHECK(nd==2,"only 2D array gives C-array of type T**");
00502     array2d.resize(d0);
00503     for(uint i=0;i<d0;i++) array2d(i)=p+i*d1;
00504     return array2d.p;
00505   }
00506 
00508   T*** getPointers(Array<T**>& array3d,Array<T*>& array2d) const{
00509     CHECK(nd==3,"only 3D array gives C-array of type T*** ");
00510     array2d.resize(d0,d1);
00511     for(uint i=0;i<d0;i++){
00512       for(uint j=0;j<d1;j++) array2d(i,j)=&operator()(i,j,0);
00513       array3d(i)=&array2d(i,0);
00514     }
00515     return array3d.p;
00516   }
00517 
00518 
00519 
00520 public:
00521 
00523   Array<T>& operator=(const T v){
00524     uint i;
00525     //if(memMove && typeid(T)==typeid(T)) memset(p,*((int*)&v),N); else
00526     for(i=0;i<N;i++) p[i]=v;
00527     return *this;
00528   }
00529 
00531   Array<T>& operator=(const Array<T>& a){
00532     IFEX( if(a.ex){ assignEx(a); return *this; } );
00533     //if(a.temp){ takeOver(*((Array<T>*)&a)); return *this; }
00534     resize(a);
00535     uint i;
00536     if(memMove) memmove(p,a.p,sizeT*N);
00537     else for(i=0;i<N;i++) p[i]=a.p[i];
00538     return *this;
00539   }
00540 
00542   void set(const char* str){
00543     std::istringstream s(str);
00544     readRaw(s);
00545   }
00546 
00547 #ifdef MT_EXPRESSIONS
00548   void checkEx() const{ if(ex) assignEx(); }
00549   void assignEx(const Array<T>& a) const;
00550   void assignEx() const;
00551   //Array<T>& operator=(const Ex& ex);
00552 #else
00553   void checkEx() const{}
00554 #endif
00555 
00557   template<class S>
00558   Array<T>& copy(const Array<S>& a){
00559     resize(a);
00560     uint i;
00561     if(memMove && typeid(T)==typeid(S)) memmove(p,a.p,sizeT*N);
00562     else for(i=0;i<N;i++) p[i]=(T)a.p[i];
00563     return *this;
00564   }
00565 
00568   void setZero(byte zero=0){
00569     CHECK(memMove,"can set array's memory to zero only if memMove option is true");
00570     memset(p,zero,sizeT*N);
00571   }
00572 
00575   void setId(int n=-1){
00576     CHECK(n!=-1 || (nd==2 && d0==d1),"need squared matrix to set to identity");
00577     if(n!=-1) resize(n,n);
00578     setZero();
00579     for(uint i=0;i<d0;i++) operator()(i,i)=1.;
00580   }
00581 
00583   template<class S>
00584   void copy(S *buffer,uint D0){
00585     resize(D0);
00586     uint i;
00587     if(memMove && typeid(T)==typeid(S))
00588       memmove(p,buffer,sizeT*d0);
00589     else for(i=0;i<d0;i++) operator()(i)=(T)buffer[i];
00590   }
00591 
00593   template<class S>
00594   void copy(S **buffer,uint D0,uint D1){
00595     resize(D0,D1);
00596     uint i,j;
00597     for(i=0;i<d0;i++){
00598       if(memMove && typeid(T)==typeid(S))
00599   memmove(p+i*d1,buffer[i],sizeT*d1);
00600       else for(j=0;j<d1;j++) operator()(i,j)=(T)buffer[i][j];
00601     }
00602   }
00603 
00606   template<class S>
00607   void copyInto(S *buffer){
00608     CHECK(nd==1,"can only copy 1D Array into 1D C-array");
00609     uint i;
00610     if(memMove && typeid(T)==typeid(S)) memmove(buffer,p,sizeT*d0);
00611     else for(i=0;i<d0;i++) buffer[i]=(S)operator()(i);
00612   }
00613 
00616   template<class S>
00617   void copyInto2D(S **buffer){
00618     CHECK(nd==2,"can only copy 2D Array into 2D C-array");
00619     uint i,j;
00620     for(i=0;i<d0;i++){
00621       if(memMove && typeid(T)==typeid(S)) memmove(buffer[i],p+i*d1,sizeT*d1);
00622       else for(j=0;j<d1;j++) buffer[i][j]=(S)operator()(i,j);
00623     }
00624   }
00625 
00627   void referTo(const T *buffer,uint n){
00628     freeMEM();
00629     reference=true;
00630     nd=1; d0=n; d1=d2=0; N=n;
00631     p=(T*)buffer;
00632     pstop=p+N;
00633   }
00634 
00636   void referTo(const Array<T>& a){
00637     freeMEM();
00638     reference=true; memMove=a.memMove; flexiMem=a.flexiMem;
00639     N=a.N; nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2;
00640     p=a.p;
00641     pstop=a.pstop;
00642   }
00643 
00645   void referToSubRange(const Array<T>& a,uint i,int I){
00646     CHECK(a.nd==1,"not implemented yet");
00647     freeMEM();
00648     reference=true; memMove=a.memMove; flexiMem=a.flexiMem;
00649     if(I==-1) I=a.N-1;
00650     N=I+1-i; nd=1; d0=N; d1=0; d2=0;
00651     p=a.p+i;
00652     pstop=p+N;
00653   }
00654 
00656   void referToSubDim(const Array<T>& a,uint dim){
00657     CHECK(a.nd>1,"can't create subarray of array less than 2 dimensions");
00658     freeMEM();
00659     reference=true; memMove=a.memMove; flexiMem=a.flexiMem;
00660     if(a.nd==2){
00661       nd=1; d0=a.d1; d1=d2=0; N=d0;
00662       p=&a(dim,0);
00663       pstop=p+N;
00664     }
00665     if(a.nd==3){
00666       nd=2; d0=a.d1; d1=a.d2; d2=0; N=d0*d1;
00667       p=&a(dim,0,0);
00668       pstop=p+N;
00669     }
00670   }
00671 
00673   void referToSubDim(const Array<T>& a,uint i,uint j){
00674     CHECK(a.nd>2,"can't create subsubarray of array less than 3 dimensions");
00675     freeMEM();
00676     reference=true; memMove=a.memMove; flexiMem=a.flexiMem;
00677     if(a.nd==3){
00678       nd=1; d0=a.d2; d1=0; d2=0; N=d0;
00679       p=&a(i,j,0);
00680       pstop=p+N;
00681     }
00682   }
00683 
00684   void redirect(const Array<T>& a,uint i){ CHECK(reference && a.nd==2,"can only do that hack with references"); p=a.p+i*a.d1; }
00685 
00687   void takeOver(Array<T>& a){
00688     freeMEM();
00689     memMove=a.memMove; flexiMem=a.flexiMem;
00690     N=a.N; nd=a.nd; d0=a.d0; d1=a.d1; d2=a.d2;
00691     p=a.p;
00692     pstop=a.pstop;
00693     M=a.M;
00694     a.reference=true;
00695     a.M=0;
00696   }
00697 
00698 
00699 public:
00700   void First(T*& i) const{ i=p; }
00701   void Last(T*& i) const{ i=(p+N)-1; }
00702   void Next(T*& i) const{ i++; }
00703   bool CheckNext(T*& i) const{ return i!=pstop; }
00704   void Prev(T*& i) const{ i--; }
00705   bool CheckPrev(T*& i) const{ return i!=p-1; }
00706 
00707 
00708 public:
00709 
00710   void permute(uint i,uint j){ T x=p[i]; p[i]=p[j]; p[j]=x; }
00711 
00713   void randomPermute(){
00714     int j,r;
00715     for(j=N-1;j>=1;j--){
00716       r=rnd(j+1);
00717       permute(r,j);
00718     }
00719   }
00720 
00721 
00722 public:
00723 
00725   void write(std::ostream& os,char *ARRAYSEP=" ",bool binary=false) const{
00726     CHECK(!binary || memMove,"binary write works only for memMoveable data");
00727     uint i,j,k;
00728     if(!IOraw){//write tag
00729       if(nd==0){ os <<"<0>"; return; }
00730       if(nd==1)  os <<"<1:" <<d0;
00731       if(nd==2)  os <<"<2:" <<d0 <<',' <<d1;
00732       if(nd==3)  os <<"<3:" <<d0 <<',' <<d1 <<',' <<d2;
00733     }
00734     //write data
00735     if(binary){
00736       os.write((char*)p,sizeT*N);
00737     }else{
00738       if(nd==1){
00739   //os <<' ';
00740   for(i=0;i<N;i++) os <<ARRAYSEP <<operator()(i);
00741       }
00742       if(nd==2) for(j=0;j<d0;j++){
00743   if(j || !IOraw) os <<"\n";
00744   for(i=0;i<d1;i++) os <<ARRAYSEP <<operator()(j,i);
00745       }
00746       if(nd==3) for(k=0;k<d0;k++){
00747         if(k || !IOraw) os <<"\n";
00748         for(j=0;j<d1;j++){
00749           for(i=0;i<d2;i++) os <<ARRAYSEP <<operator()(k,j,i);
00750           os <<"\n";
00751         }
00752       }
00753     }
00754     if(!IOraw){//write end tag
00755       os <<">";
00756       if(nd>1) os <<std::endl;
00757     }
00758   }
00759 
00761   void read(std::istream& is,bool binary=false){
00762     CHECK(!binary || memMove,"binary write works only for memMoveable data");
00763     uint d,i,j,k;
00764     char c;
00765     MT::skip(is);
00766     is.get(c);
00767     switch(c){
00768     case '<':
00769       is >>d;
00770       if(is.fail()) HALT ("could not read array tag");
00771       if(d==0){ is >>">"; return; }
00772       if(d==1){
00773         is >>":" >>i;
00774         if(is.fail()) HALT ("could not read array's dimensions");
00775   resize(i);
00776       }
00777       if(d==2){
00778         is >>":" >>i >>"," >>j;
00779   if(is.fail()) HALT ("could not read array's dimensions");
00780         resize(i,j);
00781       }
00782       if(d==3){
00783   is >>":" >>i >>"," >>j >>"," >>k;
00784   if(is.fail()) HALT ("could not read array's dimensions");
00785         resize(i,j,k);
00786       }
00787       //read data
00788       if(binary){
00789         is.read((char*)p,sizeT*N);
00790       }else{
00791   for(i=0;i<N;i++){
00792     if(is.fail()) HALT("could not read "<<i<<"-th element of an array");
00793     is >>p[i];
00794         }
00795       }
00796       is >>">";
00797       if(is.fail()) HALT ("could not read array end tag");
00798       break;
00799     case '[':
00800       resize(0);
00801       uint d=0;
00802       T x;
00803       for(;;){
00804   MT::skip(is);
00805   is.get(c);
00806   if(c==']') break;
00807   if(c==';'){ if(!d) d=N; else CHECK(!(N%d),"mis-structured array in row "<<N/d); } else is.putback(c);
00808   is >>x;
00809   if(is.fail()) HALT("error when reading "<<N<<"th element in readRaw");
00810   append(x);
00811       }
00812       if(d){
00813   CHECK(!(N%d),"mis-structured array"); 
00814   reshape(N/d,d);
00815       }
00816       break;
00817     }
00818   }
00819 
00821   const Array<T>& ioraw() const{ IOraw=true; return *this; }
00822 
00824   const Array<T>& ionoraw() const{ IOraw=false; return *this; }
00825 
00827   void readRaw(std::istream& is,uint d1=0,uint d2=0){
00828     resize(0);
00829     T x;
00830     for(;;){
00831       MT::skip(is);
00832       if(!is.good()) break;
00833       is >>x;
00834       if(is.fail()) HALT("error when reading "<<N<<"th element in readRaw");
00835       if(d1 && d2 && N>=d1*d2) break;
00836       append(x);
00837     }
00838     if(d1 && N>d1){
00839       CHECK(!(N%d1),"read "<<N<<"elements - not dividable by "<<d1);
00840       resizeCopy(N/d1,d1);
00841     }
00842   }
00843   void readRaw(const char* filename,uint d1=0,uint d2=0){
00844     std::ifstream is; MT::open(is,filename); readRaw(is,d1,d2);
00845   }
00846   void writeRaw(const char* filename){
00847     std::ofstream os; MT::open(os,filename); os <<(*this).ioraw();
00848   }
00850 
00851 #ifndef MT_EXPRESSIONS
00852 
00853   friend inline Array<T> operator~(const Array<T>& y){ Array<T> x; transpose(x,y); x.temp=true; return x; }
00854   friend inline Array<T>& operator!(Array<T>& y){ return y; }
00855   friend inline const Array<T>& operator!(const Array<T>& y){ return y; }
00857   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; }
00859   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; }
00861   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; }
00863   friend inline Array<T> operator*(const Array<T>& y,T a){ Array<T> x; scalarMultiplication(x,y,a); x.temp=true; return x; }
00865   friend inline Array<T> operator*(T a,const Array<T>& y){ Array<T> x; scalarMultiplication(x,y,a); x.temp=true; return x; }
00866 #endif
00867 };
00868 }
00869 
00870 
00871 //===========================================================================
00872 //
00874 //
00875 
00876 typedef MT::Array<double> arr;
00877 typedef MT::Array<double> doubleA;
00878 typedef MT::Array<uint> uintA;
00879 typedef MT::Array<int> intA;
00880 typedef MT::Array<char> charA;
00881 typedef MT::Array<byte> byteA;
00882 typedef MT::Array<bool> boolA;
00883 
00884 
00885 //===========================================================================
00886 //
00888 //
00889 
00891 template<class T> inline std::istream& operator>>(std::istream& is,MT::Array<T>& x){ x.read(is);return is; }
00892 template<class T> inline MT::Array<T>& operator<<(MT::Array<T>& x,const char* str){ std::istringstream ss(str); ss >>x; return x; }
00893 
00895 template<class T> inline std::ostream& operator<<(std::ostream& os,const MT::Array<T>& x){
00896   IFEX( if(x.ex) x.assignEx(); )
00897   x.write(os); return os;
00898 }
00899 
00901 template<class T,class S>
00902 inline bool samedim(const MT::Array<T>& a,const MT::Array<S>& b){
00903   return (b.nd==a.nd && b.d0==a.d0 && b.d1==a.d1 && b.d2==a.d2);
00904 }
00905 
00907 template<class T>
00908 inline bool operator==(const MT::Array<T>& v,const MT::Array<T>& w){
00909   if(!samedim(v,w)) return false;
00910   const T *iv,*iw;
00911   for(iv=v.p,iw=w.p; iv!=v.pstop; iv++,iw++)
00912     if (*iv != *iw) return false;
00913   return true;
00914 }
00915 
00917 template<class T>
00918 inline bool operator!=(const MT::Array<T>& v,const MT::Array<T>& w){
00919   return !(v==w);
00920 }
00921 
00923 template<class T>
00924 inline bool operator<(const MT::Array<T>& v,const MT::Array<T>& w){
00925   if(v.N==w.N){
00926     for(uint i=0;i<v.N;i++){
00927       if(v.p[i]>w.p[i]) return false;
00928       if(v.p[i]<w.p[i]) return true;
00929     }
00930     return false; //they are equal
00931   }
00932   return v.N<w.N;
00933 }
00934 
00935 
00936 //===========================================================================
00937 //
00939 //
00940 
00941 #ifndef MT_doxy // exclude these macros when generating the documentation
00942                 // (doxygen can't handle them...)
00943 
00944 //---------- unary operators
00945 
00946 #define UnaryOperator( op )                                               \
00947 template<class T>                                                         \
00948 inline MT::Array<T> operator op (const MT::Array<T>& y){                  \
00949   MT::Array<T> x;                                                         \
00950   x.resize(y);                  \
00951   T *ix=x.p;                                                              \
00952   const T *iy=y.p;                                                        \
00953   for(; iy!=y.pstop; iy++,ix++) *ix= op *iy;            \
00954   x.temp=true; return x;                                                  \
00955 }
00956 
00957 //UnaryOperator( ! )
00958 //UnaryOperator( ~ )
00959 //UnaryOperator( + )
00960 UnaryOperator( - )
00961 #undef UnaryOperator
00962 
00963 
00964 //---------- binary operators
00965 
00966 #ifndef MT_EXPRESSIONS
00967 #define BinaryOperator( op )                                              \
00968 template<class T>                                                         \
00969 inline MT::Array<T> operator op(const MT::Array<T>& v,const MT::Array<T>& w){ \
00970   v.checkEx(); w.checkEx(); \
00971   CHECK(v.N==w.N,               \
00972   "binary operator on different array dimensions ("<<v.N<<", "<<w.N<<")");\
00973   MT::Array<T> z;                                                         \
00974   z.resize(v);                  \
00975   T *iz=z.p;                  \
00976   const T *iw=w.p,*iv=v.p;              \
00977   for(; iv!=v.pstop; iv++,iw++,iz++) *iz = *iv op *iw;                    \
00978   z.temp=true; return z;                                                               \
00979 }                                                                         \
00980                                                                           \
00981 template<class T>                                                         \
00982 inline MT::Array<T> operator op (const MT::Array<T>& v,T w){              \
00983   v.checkEx(); \
00984   MT::Array<T> z;                                                         \
00985   z.resize(v);                  \
00986   T *iz=z.p;                  \
00987   const T *iv=v.p;                \
00988   for(; iv!=v.pstop; iv++,iz++) *iz = *iv op w;                           \
00989   z.temp=true; return z;                                                               \
00990 }                                                                         \
00991                                                                           \
00992 template<class T>                                                         \
00993 inline MT::Array<T> operator op ( T v, const MT::Array<T>& w ){           \
00994   w.checkEx(); \
00995   MT::Array<T> z;                                                         \
00996   z.resize(w);                  \
00997   T *iz=z.p;                  \
00998   const T *iw=w.p;                \
00999   for(; iw!=w.pstop; iw++,iz++) *iz = v op *iw;                           \
01000   z.temp=true; return z;                                                               \
01001 }
01002 
01003 BinaryOperator( || )
01004 BinaryOperator( && )
01005 BinaryOperator( | )
01006   //BinaryOperator( ^ ) I reserve this for the outer product!
01007 BinaryOperator( & )
01008 #ifndef MT_EXPRESSIONS
01009 BinaryOperator( + )
01010 BinaryOperator( - )
01011 #endif
01012   //BinaryOperator( * ) I reserve this operator for the inner product!
01013 BinaryOperator( / )
01014 BinaryOperator( % )
01015 #undef BinaryOperator
01016 #endif
01017 
01018 
01019 //---------- compound assignment operators
01020 
01021 #define CompoundAssignmentOperator( op )                                  \
01022 template<class T,class S>                                                 \
01023 inline MT::Array<T>& operator op (MT::Array<T>& x,const MT::Array<S>& y){ \
01024   y.checkEx(); \
01025   CHECK(x.N==y.N,               \
01026   "binary operator on different array dimensions ("<<x.N<<", "<<y.N<<")");\
01027   T *ix=x.p;                  \
01028   const S *iy=y.p;                \
01029   for(; ix!=x.pstop; ix++,iy++) *ix op *iy;                               \
01030   x.temp=true; return x;                                                               \
01031 }                                                                         \
01032                                                                           \
01033 template<class T>                                                         \
01034 inline MT::Array<T>& operator op ( MT::Array<T>& x, T y ){                \
01035   T *ix=x.p;                  \
01036   for(; ix!=x.pstop; ix++) *ix op y;                                      \
01037   x.temp=true; return x;                                                               \
01038 }                                                                         \
01039 
01040 CompoundAssignmentOperator( |= )
01041 CompoundAssignmentOperator( ^= )
01042 CompoundAssignmentOperator( &= )
01043 CompoundAssignmentOperator( += )
01044 CompoundAssignmentOperator( -= )
01045 CompoundAssignmentOperator( *= )
01046 CompoundAssignmentOperator( /= )
01047 CompoundAssignmentOperator( %= )
01048 #undef CompoundAssignmentOperator
01049 
01050 
01051 //---------- unary functions
01052 
01053 #define UnaryFunction( func )                                             \
01054 template<class T>                                                         \
01055 inline MT::Array<T> func (const MT::Array<T>& v){                         \
01056   v.checkEx(); \
01057   MT::Array<T> z;                                                         \
01058   z.resize(v);                                                            \
01059   T *iz=z.p;                  \
01060   const T *iv=v.p;                \
01061   for(; iv!=v.pstop; iv++, iz++) *iz = ::func( *iv );                     \
01062   z.temp=true; return z;                                                               \
01063 }                   \
01064 
01065 // trigonometric functions
01066 UnaryFunction( acos )
01067 UnaryFunction( asin )
01068 UnaryFunction( atan )
01069 UnaryFunction( cos )
01070 UnaryFunction( sin )
01071 UnaryFunction( tan )
01072 
01073 // hyperbolic functions
01074 UnaryFunction( cosh )
01075 UnaryFunction( sinh )
01076 UnaryFunction( tanh )
01077 UnaryFunction( acosh )
01078 UnaryFunction( asinh )
01079 UnaryFunction( atanh )
01080 
01081 // exponential and logarithmic functions
01082 UnaryFunction( exp )
01083 UnaryFunction( log )
01084 UnaryFunction( log10 )
01085 
01086 // power functions
01087 UnaryFunction( sqrt )
01088 UnaryFunction( cbrt )
01089 
01090 // nearest integer and absolute value
01091 UnaryFunction( ceil )
01092 UnaryFunction( fabs )
01093 UnaryFunction( floor )
01094 #undef UnaryFunction
01095 
01096 
01097 //---------- binary functions
01098 
01099 #define BinaryFunction( func )                                            \
01100 template<class T>                                                         \
01101 inline MT::Array<T> func(const MT::Array<T>& v,const MT::Array<T>& w){    \
01102   v.checkEx(); w.checkEx(); \
01103   CHECK(v.N==w.N,               \
01104   "binary operator on different array dimensions ("<<v.N<<", "<<w.N<<")");\
01105   MT::Array<T> z;                                                         \
01106   z.resize(v);                                                            \
01107   for(uint i=z.N;i--; ) z.p[i]= func(v.p[i],w.p[i]);        \
01108   z.temp=true; return z;                                                               \
01109 }                                                                         \
01110                                                                           \
01111 template<class T>                                                         \
01112 inline MT::Array<T> func(const MT::Array<T>& v,T w){                      \
01113   v.checkEx(); \
01114   MT::Array<T> z;                                                         \
01115   z.resize(v);                                                            \
01116   for(uint i=z.N;i--; ) z.p[i]= func(v.p[i],w);       \
01117   z.temp=true; return z;                                                               \
01118 }                                                                         \
01119                                                                           \
01120 template<class T>                                                         \
01121 inline MT::Array<T> func(T v,const MT::Array<T>& w){                      \
01122   w.checkEx(); \
01123   MT::Array<T> z;                                                         \
01124   z.resize(w);                                                            \
01125   for(uint i=z.N;i--; ) z.p[i]= func(v,w.p[i]);       \
01126   z.temp=true; return z;                                                               \
01127 }
01128 
01129 BinaryFunction( atan2 )
01130 BinaryFunction( pow )
01131 BinaryFunction( fmod )
01132 #undef BinaryFunction
01133 
01134 #endif //(doxygen exclusion)
01135 
01136 
01137 
01138 //===========================================================================
01139 //
01141 //
01142 
01143 #ifdef MT_MSVC
01144 #  ifdef plus
01145 #    undef plus
01146 #  endif
01147 #endif
01148 
01150 template<class T>
01151 void plus(MT::Array<T>& x,const MT::Array<T>& y,const MT::Array<T>& z){
01152   CHECK(y.N==z.N,"must have same size for adding!");
01153   uint i,n=y.N;
01154   x.resize(y);
01155   for(i=0;i<n;i++) x.p[i]=y.p[i]+z.p[i];
01156 }
01157 
01159 template<class T>
01160 void scalarPlus(MT::Array<T>& x,const MT::Array<T>& y,T a){
01161   uint i,n=y.N;
01162   x.resize(y);
01163   for(i=0;i<n;i++) x.p[i]=y.p[i]+a;
01164 }
01165 
01167 template<class T>
01168 void plus(MT::Array<T>& x,T a,const MT::Array<T>& y,T b,const MT::Array<T>& z){
01169   CHECK(y.N==z.N,"must have same size for adding!");
01170   uint i,n=y.N;
01171   x.resize(y);
01172   for(i=0;i<n;i++) x.p[i]=a*y.p[i]+b*z.p[i];
01173 }
01174 
01176 template<class T>
01177 void minus(MT::Array<T>& x,const MT::Array<T>& y,const MT::Array<T>& z){
01178   CHECK(y.N==z.N,"must have same size for subtracting!");
01179   uint i,n=y.N;
01180   x.resize(y);
01181   for(i=0;i<n;i++) x.p[i]=y.p[i]-z.p[i];
01182 }
01183 
01185 template<class T>
01186 void mult(MT::Array<T>& x,const MT::Array<T>& y,const MT::Array<T>& z){
01187   CHECK(y.N==z.N,"must have same size for element-wise multiplication!");
01188   uint i,n=y.N;
01189   x.resize(y);
01190   for(i=0;i<n;i++) x.p[i]=y.p[i]*z.p[i];
01191 }
01192 
01194 template<class T>
01195 void div(MT::Array<T>& x,const MT::Array<T>& y,const MT::Array<T>& z){
01196   CHECK(y.N==z.N,"must have same size for element-wise multiplication!");
01197   uint i,n=y.N;
01198   x.resize(y);
01199   for(i=0;i<n;i++) x.p[i]=y.p[i]/z.p[i];
01200 }
01201 
01202 
01203 //===========================================================================
01204 //
01206 //
01207 
01209 inline MT::Array<double> Identity(uint dim){
01210   MT::Array<double> z;
01211   z.setId(dim);
01212   return z;
01213 }
01214 
01216 inline void makeSymmetric(MT::Array<double>& A){
01217   CHECK(A.nd==2 && A.d0==A.d1,"not symmetric");
01218   uint n=A.d0,i,j;
01219   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));
01220 }
01221 
01223 inline void transpose(MT::Array<double>& A){
01224   CHECK(A.nd==2 && A.d0==A.d1,"not symmetric");
01225   uint n=A.d0,i,j;
01226   double z;
01227   for(i=1;i<n;i++) for(j=0;j<i;j++){ z=A(j,i); A(j,i)=A(i,j); A(i,j)=z; }
01228 }
01229 
01231 template<class T>
01232 void transpose(MT::Array<T>& x,const MT::Array<T>& y){
01233   CHECK(y.nd<=2,"can only transpose up to 2D arrays");
01234   if(y.nd==2){
01235     uint i,j,d0=y.d1,d1=y.d0;
01236     x.resize(d0,d1);
01237     for(i=0;i<d0;i++) for(j=0;j<d1;j++) x.p[i*d1+j]=y.p[j*d0+i];
01238     return;
01239   }
01240   if(y.nd==1){
01241     x=y;
01242     x.resizeCopy(1,y.N);
01243     return;
01244   }
01245   HALT("transpose not implemented for this dims");
01246 }
01247 
01249 template<class T>
01250 inline void diag(MT::Array<T>& x,const MT::Array<T>& y){
01251   CHECK(y.nd==2 && y.d0==y.d1,"can only give diagonal of symmetric 2D matrix");
01252   x.resize(y.d0);
01253   uint i;
01254   for(i=0;i<x.d0;i++) x(i)=y(i,i);
01255 }
01256 
01258 template<class T>
01259 inline void setDiagonal(MT::Array<T>& x,const MT::Array<T>& v){
01260   CHECK(v.nd==1,"can only give diagonal of 1D array");
01261   x.resize(v.d0,v.d0);
01262   x.setZero();
01263   uint i;
01264   for(i=0;i<v.d0;i++) x(i,i)=v(i);
01265 }
01266 
01268 template<class T>
01269 inline void setDiagonal(MT::Array<T>& x,uint d,double v){
01270   x.resize(d,d);
01271   x.setZero();
01272   uint i;
01273   for(i=0;i<d;i++) x(i,i)=v;
01274 }
01275 
01277 template<class T>
01278 void inverse2d(MT::Array<T>& invA,const MT::Array<T>& A){
01279   invA.resize(2,2);
01280   invA(0,0)=A(1,1); invA(1,1)=A(0,0); invA(0,1)=-A(0,1); invA(1,0)=-A(1,0);
01281   invA/=A(0,0)*A(1,1)-A(0,1)*A(1,0);
01282 }
01283 
01285 template<class T>
01286 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){
01287   CHECK(A.nd==2 && B.nd==2 && C.nd==2 && D.nd==2,"");
01288   CHECK(A.d0==B.d0 && A.d1==C.d1 && B.d1==D.d1 && C.d0==D.d0,"");
01289   uint i,j,a=A.d0,b=A.d1;
01290   X.resize(A.d0+C.d0,A.d1+B.d1);
01291   for(i=0;i<A.d0;i++) for(j=0;j<A.d1;j++) X(i  ,j  )=A(i,j);
01292   for(i=0;i<B.d0;i++) for(j=0;j<B.d1;j++) X(i  ,j+b)=B(i,j);
01293   for(i=0;i<C.d0;i++) for(j=0;j<C.d1;j++) X(i+a,j  )=C(i,j);
01294   for(i=0;i<D.d0;i++) for(j=0;j<D.d1;j++) X(i+a,j+b)=D(i,j);
01295 }
01296 
01298 template<class T>
01299 void blockMatrix(MT::Array<T>& X,const MT::Array<T>& A,const MT::Array<T>& B,const MT::Array<T>& C,const double& D){
01300   CHECK(A.nd==2 && B.nd==1 && C.nd==1,"");
01301   CHECK(A.d0==B.d0 && A.d1==C.d0,"");
01302   uint i,j,a=A.d0,b=A.d1;
01303   X.resize(A.d0+1,A.d1+1);
01304   for(i=0;i<a;i++) for(j=0;j<b;j++) X(i,j)=A(i,j);
01305   for(i=0;i<a;i++) X(i,b)=B(i);
01306   for(j=0;j<b;j++) X(a,j)=C(j);
01307   X(a,b)=D;
01308 }
01309 
01310 
01311 //===========================================================================
01312 //
01314 //
01315 
01317 template<class T> 
01318 inline T sqrDistance(const MT::Array<T>& v, const MT::Array<T>& w){
01319   CHECK(v.N==w.N, 
01320     "sqrDistance on different array dimensions ("<<v.N<<", "<<w.N<<")");
01321   T d,t(0);
01322   for(uint i=v.N;i--;){ d=v.p[i]-w.p[i]; t+=d*d; }
01323   return t;
01324 }
01325 
01326 template<class T>
01327 inline T maxDiff(const MT::Array<T>& v,const MT::Array<T>& w){
01328   CHECK(v.N==w.N, 
01329     "maxDiff on different array dimensions ("<<v.N<<", "<<w.N<<")");
01330   T d,t(0);
01331   for(uint i=v.N;i--;){ d=fabs(v.p[i]-w.p[i]); if(d>t) t=d; }
01332   return t;
01333 }
01334 
01336 template<class T> 
01337 inline T sqrDistance(const MT::Array<T>& v, const MT::Array<T>& w,const MT::Array<bool>& mask){
01338   CHECK(v.N==w.N, 
01339     "sqrDistance on different array dimensions ("<<v.N<<", "<<w.N<<")");
01340   T d,t(0);
01341   for(uint i=v.N;i--;) if(mask(i)){ d=v.p[i]-w.p[i]; t+=d*d; }
01342   return t;
01343 }
01344 
01346 template<class T>
01347 inline T sqrDistance(const MT::Array<T>& g,const MT::Array<T>& v,const MT::Array<T>& w){
01348   MT::Array<T> d;
01349   minus(d,v,w);
01350   return scalarProduct(g,d,d);
01351 }
01352 
01354 template<class T> 
01355 inline T euclideanDistance(const MT::Array<T>& v, const MT::Array<T>& w){
01356   return ::sqrt(sqrDistance(v,w));
01357 }
01358 
01359 
01360 
01361 //===========================================================================
01362 //
01364 //
01365 
01367 template<class T>
01368 inline T sum(const MT::Array<T>& v){
01369   T t(0);
01370   for(uint i=v.N; i--; t+=v.p[i]);
01371   return t;
01372 }
01373 
01375 template<class T>
01376 inline T sumOfAbs(const MT::Array<T>& v){
01377   T t(0);
01378   for(uint i=v.N; i--; t+=fabs(v.p[i]));
01379   return t;
01380 }
01381 
01383 template<class T>
01384 inline T sumOfSqr(const MT::Array<T>& v){
01385   T t(0);
01386   for(uint i=v.N; i--; t+=v.p[i]*v.p[i]);
01387   return t;
01388 }
01389 
01391 template<class T>
01392 inline T norm(const MT::Array<T>& v){ return sqrt(sumOfSqr(v)); }
01393 
01395 template<class T>
01396 inline T mean(const MT::Array<T>& v){ return sum(v)/v.N; }
01397 
01398 template<class T>
01399 inline T var(const MT::Array<T>& v){ double m=mean(v); return (sumOfSqr(v)-m*m/v.N)/v.N; }
01400 
01402 template<class T>
01403 inline T trace(const MT::Array<T>& v){
01404   CHECK(v.nd==2 && v.d0==v.d1,"only for squared matrix");
01405   T t(0);
01406   for(uint i=0;i<v.d0;i++) t+=v(i,i);
01407   return t;
01408 }
01409 
01410 
01411 //===========================================================================
01412 //
01414 //
01415 
01416 namespace MT{
01418 extern bool useLapack;
01419 }
01420 
01421 #ifdef MT_LAPACK
01422 extern "C"{
01423 #include"cblas.h"
01424 #include"f2c.h"
01425 #include"clapack.h"
01426 }
01427 
01428 //matrix-matrix multiplication
01429 inline void blasMM(void *pX,void *pA,void *pB){
01430   doubleA &X=*(doubleA*)pX,&A=*(doubleA*)pA,&B=*(doubleA*)pB;
01431   CHECK(A.d1==B.d0,"matrix multiplication: wrong dimensions");
01432   X.resize(A.d0,B.d1);
01433   cblas_dgemm(CblasRowMajor,
01434         CblasNoTrans,CblasNoTrans,
01435         A.d0,B.d1,A.d1,
01436         1.,A.p,A.d1,
01437         B.p,B.d1,1.,
01438         X.p,X.d1);
01439 #ifdef MT_Linux
01440   //bad joke? the 10th element is always nan...
01441   if(X.N>9){
01442     X(0,9)=0.; for(uint k=0;k<A.d1;k++) X(0,9) += A(0,k) * B(k,9);
01443   }
01444 #endif
01445 
01446 #if 0 //test
01447   doubleA Y(A.d0,B.d1);
01448   uint i,j,k;
01449   Y.setZero();
01450   for(i=0;i<Y.d0;i++) for(j=0;j<Y.d1;j++) for(k=0;k<A.d1;k++)
01451     Y(i,j) += A(i,k) * B(k,j);
01452   std::cout  <<"blasMM error = " <<sqrDistance(X,Y) <<std::endl;
01453 #endif
01454 }
01455 
01456 inline uint lapackSVD(const doubleA& A,
01457           doubleA& U,
01458           doubleA& d,
01459           doubleA& Vt){
01460   static doubleA Atmp,work;
01461   Atmp=A;
01462   //transpose(Atmp,A);
01463   integer M=A.d0,N=A.d1,D=M<N?M:N;
01464   U.resize(M,D);
01465   d.resize(D);
01466   Vt.resize(D,N);
01467   work.resize(10*(M+N));
01468   integer info,wn=work.N;
01469   dgesvd_("S", "S", &N, &M, Atmp.p, &N,
01470     d.p, Vt.p, &N, U.p, &D,
01471     work.p,&wn, &info);
01472   if(info){
01473     std::cerr <<"LAPACK SVD error info = " <<info <<std::endl;
01474   }
01475 
01476 #if 0 //test!
01477   MT::useLapack=false;
01478   doubleA dD,I;
01479   setDiagonal(dD,d);
01480   //cout <<U <<dD <<Vt;
01481   //Atmp = V * D * U;
01482   Atmp = U * dD * Vt;
01483   std::cout <<"SVD is correct:  " <<sqrDistance(Atmp,A) <<' ' <<endl;
01484   setIdentity(I,D);
01485   std::cout <<"U is orthogonal: " <<sqrDistance(~U * U,I) <<' ' <<endl;
01486   std::cout <<"Vt is orthogonal: " <<sqrDistance(Vt * ~Vt,I) <<endl;
01487   MT::useLapack=true;
01488 #endif
01489   return D;
01490 }
01491 #endif
01492 
01493 
01494 //===========================================================================
01495 //
01497 //
01498 
01499 uint own_svd_implementation(const doubleA& A,
01500           doubleA& U,
01501           doubleA& w,
01502           doubleA& V,
01503           bool sort);
01504 
01509 inline uint svd(const doubleA& A,doubleA& U,doubleA& d,doubleA& V,bool sort=true){
01510 #ifdef MT_LAPACK
01511   if(MT::useLapack) return lapackSVD(A,U,d,V);
01512   V=~V;
01513 #endif
01514   uint r=own_svd_implementation(A,U,d,V,sort);
01515   return r;
01516 }
01517 inline void svd(const doubleA& A,doubleA& U,doubleA& V){
01518   doubleA d,D;
01519   ::svd(A,U,d,V);
01520   D.resize(d.N,d.N); D=0.;
01521   for(uint i=0;i<d.N;i++) D(i,i)=::sqrt(d(i));
01522   U=U*D;
01523   V=V*D;
01524 }
01525 
01527 uint inverse(doubleA& inverse,const doubleA& A);
01528 
01529 inline doubleA inverse(const doubleA& A){ doubleA B; inverse(B,A); return B; }
01530 
01532 double determinant(const doubleA& A);
01533 
01536 double cofactor(const doubleA& A,uint i,uint j);
01537 
01538 
01539 //===========================================================================
01540 //
01542 //
01543 
01545 template<class T>
01546 void scalarMultiplication(MT::Array<T>& x,const MT::Array<T>& y,const T& a){
01547   x.resize(y);
01548   T *ix=x.p;
01549   const T *iy=y.p;
01550   for(; ix!=x.pstop; ix++,iy++) *ix = *iy * a;
01551 }
01552 
01556 template<class T> 
01557 void innerProduct(MT::Array<T>& x,const MT::Array<T>& y, const MT::Array<T>& z){
01558   /*
01559   if(y.nd==2 && z.nd==2 && y.N==z.N && y.d1==1 && z.d1==1){  //elem-wise
01560     HALT("make element-wise multiplication explicite!");
01561     mult(x,y,z);
01562     return;
01563   }
01564   */
01565   if(y.nd==2 && z.nd==2){ //plain matrix multiplication
01566 #ifdef MT_LAPACK
01567     if(MT::useLapack && typeid(T)==typeid(double)){ blasMM(&x,(void*)&y,(void*)&z); return; }
01568 #endif
01569     CHECK(y.d1==z.d0,"wrong dimensions for inner product");
01570     uint i,j,d0=y.d0,d1=z.d1,dk=y.d1;
01571     T *a,*astop,*b,*c;
01572     x.resize(d0,d1); x.setZero();
01573     c=x.p;
01574     for(i=0;i<d0;i++) for(j=0;j<d1;j++){
01575       //for(s=0.,k=0;k<dk;k++) s+=y.p[i*dk+k]*z.p[k*d1+j];
01576       //this is faster:
01577       a=y.p+i*dk; astop=a+dk; b=z.p+j;
01578       for(;a<astop; a++, b+=d1) (*c)+=(*a) * (*b);
01579       c++;
01580     }
01581     return;
01582   }
01583   if(y.nd==2 && z.nd==1){ //matrix x vector -> vector
01584     CHECK(y.d1==z.d0,"wrong dimensions for inner product");
01585     uint i,k,d0=y.d0,dk=y.d1;
01586     T *a,*b;
01587     T s;
01588     x.resize(d0);
01589     for(i=0;i<d0;i++){
01590       //for(s=0.,k=0;k<dk;k++) s+=y.p[i*dk+k]*z.p[k];
01591       //this is faster:
01592       a=y.p+i*dk; b=z.p;
01593       for(s=0,k=0;k<dk;k++){ s+=(*a) * (*b); a++; b++; }
01594       x.p[i]=s;
01595     }
01596     return;
01597   }
01598   if(y.nd==1 && z.nd==2 && z.d0==1){ //vector x vector^T -> matrix (outer product)
01599     CHECK(y.d0==z.d1,"wrong dimensions for inner product");
01600     uint i,j,d0=y.d0,d1=z.d1;
01601     x.resize(d0,d1);
01602     for(i=0;i<d0;i++) for(j=0;j<d1;j++) x(i,j)=y(i)*z(0,j);
01603     return;
01604   }
01605   if(y.nd==1 && z.nd==2){ //vector^T x matrix -> vector^T
01606     CHECK(y.d0==z.d0,"wrong dimensions for inner product");
01607     uint i,k,d0=z.d1,dk=y.d0;
01608     x.resize(d0);
01609     T s;
01610     for(i=0;i<d0;i++){
01611       for(s=0,k=0;k<dk;k++) s+=y.p[k]*z.p[k*d0+i];
01612       x.p[i]=s;
01613     }
01614     return;
01615   }
01616   if(y.nd==1 && z.nd==1 && z.N==1){ //vector multiplied with scalar (disguised as 1D vector)
01617     uint k,dk=y.N;
01618     x.resize(y.N);
01619     for(k=0;k<dk;k++) x.p[k]=y.p[k]*z.p[0];
01620     return;
01621   }
01622   if(y.nd==1 && z.nd==1){ //should be scalar product, but be careful
01623     HALT("what do you want? scalar product or element wise multiplication?");
01624     CHECK(y.d0==z.d0,"wrong dimensions for inner product");
01625     uint k,dk=y.d0;
01626     x.resize(1);
01627     T s;
01628     for(s=0,k=0;k<dk;k++) s+=y.p[k]*z.p[k];
01629     x.p[0]=s;
01630     return;
01631   }
01632   HALT("inner product - not yet implemented for these dimensions");
01633 }
01634 
01637 template<class T> 
01638 void outerProduct(MT::Array<T>& x,const MT::Array<T>& y, const MT::Array<T>& z){
01639   if(y.nd==1 && z.nd==1){
01640     uint i,j,d0=y.d0,d1=z.d0;
01641     x.resize(d0,d1);
01642     for(i=0;i<d0;i++) for(j=0;j<d1;j++) x.p[i*d1+j]=y.p[i]*z.p[j];
01643     return;
01644   }
01645   HALT("outer product - not yet implemented for these dimensions");
01646 }
01647 
01649 template<class T> 
01650 inline T scalarProduct(const MT::Array<T>& v, const MT::Array<T>& w){
01651   CHECK(v.N==w.N, 
01652     "scalar product on different array dimensions ("<<v.N<<", "<<w.N<<")");
01653   T t(0);
01654   for(uint i=v.N; i--; t+=v.p[i]*w.p[i]);
01655   return t;
01656 }
01657 
01659 template<class T> 
01660 inline T scalarProduct(const MT::Array<T>& g,const MT::Array<T>& v, const MT::Array<T>& w){
01661   CHECK(v.N==w.N && g.nd==2 && g.d0==v.N && g.d1==w.N,
01662     "scalar product on different array dimensions ("<<v.N<<", "<<w.N<<")");
01663   T t(0);
01664   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];
01665   return t;
01666 }
01667 
01668 
01669 //===========================================================================
01670 //
01672 //
01673 
01676 template<class T>
01677 inline T entropy(const MT::Array<T>& v){
01678   T t(0);
01679   for(uint i=v.N; i--;) if(v.p[i]) t-=v.p[i]*::log(v.p[i]);
01680   return t/MT_LN2;
01681 }
01682 
01684 template<class T>
01685 inline T normalizeDist(MT::Array<T>& v){
01686   double Z=sum(v);
01687   if(Z) v/=Z; else v=1./(double)v.N;
01688   return Z;
01689 }
01690 
01691 template<class T>
01692 inline void eliminate(MT::Array<T>& x,const MT::Array<T>& y,uint d){
01693   CHECK(y.nd==2,"only implemented for 2D yet");
01694   uint i,j;
01695   if(d==1){
01696     x.resize(y.d0); x=0.;
01697     for(i=0;i<y.d0;i++) for(j=0;j<y.d1;j++) x(i)+=y(i,j);
01698   }
01699   if(d==0){
01700     x.resize(y.d1); x=0.;
01701     for(i=0;i<y.d0;i++) for(j=0;j<y.d1;j++) x(j)+=y(i,j);
01702   }
01703 }
01704 
01705 template<class T>
01706 inline void eliminate(MT::Array<T>& x,const MT::Array<T>& y,uint d,uint e){
01707   CHECK(y.nd==3,"only implemented for 3D yet");
01708   uint i,j,k;
01709   if(d==1 && e==2){
01710     x.resize(y.d0); x=0.;
01711     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);
01712   }
01713   if(d==0 && e==2){
01714     x.resize(y.d1); x=0.;
01715     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);
01716   }
01717   if(d==0 && e==1){
01718     x.resize(y.d2); x=0.;
01719     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);
01720   }
01721 }
01722 
01727 void SUS(const doubleA& p,uint n,uintA& s);
01728 
01730 uint SUS(const doubleA& p);
01731 
01732 
01733 //===========================================================================
01734 //
01736 //
01737 
01739 template<class T>
01740 inline T product(const MT::Array<T>& v){
01741   T t(1);
01742   for(uint i=v.N; i--; t *= v.p[i]);
01743   return t;
01744 }
01745 
01746 
01747 //===========================================================================
01748 //
01750 //
01751 
01753 template<class T>
01754 inline T minA(const MT::Array<T>& v){
01755   CHECK(v.N>0,"");
01756   T t(v.p[0]);
01757   for(uint i=1; i<v.N; ++i) if(v.p[i]<t) t=v.p[i];
01758   return t;
01759 }
01760 
01762 template<class T>
01763 inline T minA(const MT::Array<T>& v, uint & ind, uint start=0, uint end=0){
01764   CHECK(v.N>0,"");
01765   CHECK(v.N>start,"");  
01766   CHECK(v.N>=end,"");
01767   CHECK(end>=start,"");
01768   T t(v(start));
01769   ind=start;
01770   if (end==0) end=v.N;
01771   for(uint i=start; i<end; ++i) if(v.p[i]<t){
01772     t  =v.p[i];
01773     ind=i;
01774   }
01775   return t;
01776 }
01777 
01779 template<class T>
01780 inline T maxA(const MT::Array<T>& v){
01781   CHECK(v.N>0,"");
01782   T t(v.p[0]);
01783   for(uint i=1; i<v.N; ++i) if(v.p[i]>t) t=v.p[i];
01784   return t;
01785 }
01786 
01788 template<class T>
01789 inline T maxA(const MT::Array<T>& v, uint & ind, uint start=0, uint end=0){
01790   CHECK(v.N>0,"");
01791   CHECK(v.N>start,"");
01792   CHECK(v.N>=end,"");
01793   CHECK(end>=start,"");
01794   T t(v(start));
01795   ind=start;
01796   if(end==0){
01797     end=v.N;
01798   }
01799   for(uint i=start; i<end; ++i) if(v.p[i]>t){
01800     t  =v.p[i];
01801     ind=long(i);
01802   }
01803   return t;
01804 }
01805 
01807 template<class T>
01808 inline void minmaxA(const MT::Array<T>& v, T& minVal, T& maxVal){
01809   if(v.N>0){
01810     minVal=maxVal=v.p[0];
01811     for(uint i=1; i<v.N; ++i){
01812       if(v.p[i]<minVal) minVal=v.p[i];
01813       else if(v.p[i]>maxVal) maxVal=v.p[i];
01814     }
01815   }
01816 }
01817 
01818 template<class T>
01819 inline T absMax(const MT::Array<T>& x){
01820   CHECK(x.N>0,"");
01821   T t(x.p[0]);
01822   for(uint i=1; i<x.N; ++i) if(fabs(x.p[i])>t) t=fabs(x.p[i]);
01823   return t;
01824 }
01825 
01826 
01827 //===========================================================================
01828 //
01830 //
01831 
01833 template<class T>
01834 void rndInt(MT::Array<T>& a,int low=0.,int high=1.){
01835   for(uint i=0;i<a.N;i++) a.p[i]=low+(int)rnd.num(1+high-low);
01836 }
01837 
01839 template<class T>
01840 void rndUni(MT::Array<T>& a,double low=0.,double high=1.){
01841   for(uint i=0;i<a.N;i++) a.p[i]=rnd.uni(low,high);
01842 }
01843 
01845 template<class T>
01846 void rndGauss(MT::Array<T>& a,double stdDev,bool add=false){
01847   if(a.nd==2 && a.d0!=a.d1)  //assumes that it is a batch of data
01848     stdDev/=::sqrt((double)a.d1);
01849   else
01850     stdDev/=::sqrt((double)a.N);
01851   if(!add) for(uint i=0;i<a.N;i++) a.p[i]=rnd.gauss(stdDev);
01852   else     for(uint i=0;i<a.N;i++) a.p[i]+=rnd.gauss(stdDev);
01853 }
01854 
01856 template<class T>
01857 MT::Array<T>& rndGauss(double stdDev,uint dim){
01858   static MT::Array<T> z;
01859   stdDev/=::sqrt(dim);
01860   z.resize(dim);
01861   rndGauss(z,stdDev);
01862   return z;
01863 }
01864 
01868 template<class T>
01869 uint softMax(const MT::Array<T>& a,doubleA& soft,double beta){
01870   double norm=0.,r;
01871   uint i; int sel=-1;
01872   soft.resize(a);
01873   for(i=0;i<a.N;i++){
01874     soft(i)=exp(beta*a(i));
01875     norm+=soft(i);
01876   }
01877   r=rnd.uni();
01878   for(i=0;i<a.N;i++){
01879     soft(i)/=norm;
01880     r-=soft(i);
01881     if(sel==-1 && r<0.) sel=i;
01882   }
01883   return sel;
01884 }
01885 
01886 
01887 //===========================================================================
01888 //
01890 //
01891 
01892 template<class T>
01893 void grid(MT::Array<T>& a,uint dim,T lo,T hi,uint steps){
01894   uint i,j;
01895   if(dim==1){
01896     a.resize(steps+1,1);
01897     for(i=0;i<a.d0;i++) a(i,0)=lo+(hi-lo)*i/steps;
01898     return;
01899   }
01900   if(dim==2){
01901     a.resize(steps+1,steps+1,2);
01902     for(i=0;i<a.d0;i++) for(j=0;j<a.d1;j++){
01903       a(i,j,0)=lo+(hi-lo)*j/steps;
01904       a(i,j,1)=lo+(hi-lo)*i/steps;
01905     }
01906     a.reshape(a.d0*a.d1,2);
01907     return;
01908   }
01909   HALT("not implemented yet");
01910 }
01911 
01912 
01913 //===========================================================================
01914 //
01916 //
01917 
01918 #ifdef MT_EXPRESSIONS
01919 
01920 inline doubleA operator*(double a,const doubleA& b){ return doubleA(new MT::Ex((void*)&b,0,a,false)); }
01921 inline doubleA operator*(const doubleA& a,double b){ return doubleA(new MT::Ex((void*)&a,0,b,false)); }
01922 inline doubleA operator/(const doubleA& a,double b){ return doubleA(new MT::Ex((void*)&a,0,1./b,false)); }
01923 inline doubleA operator+(double a,const doubleA& b){ return doubleA(new MT::Ex((void*)&b,a,1,false)); }
01924 inline doubleA operator+(const doubleA& a,double b){ return doubleA(new MT::Ex((void*)&a,b,1,false)); }
01925 inline doubleA operator-(double a,const doubleA& b){ return doubleA(new MT::Ex((void*)&b,-a,1,false)); }
01926 inline doubleA operator-(const doubleA& a,double b){ return doubleA(new MT::Ex((void*)&a,-b,1,false)); }
01927 inline doubleA operator~(const doubleA& a){ return doubleA(new MT::Ex((void*)&a,0,1,true)); }
01928 
01929 inline doubleA operator*(const doubleA& y,const doubleA& z){ return doubleA(new MT::Ex(MT::PROD,(void*)&y,(void*)&z)); }
01930 inline doubleA operator^(const doubleA& y,const doubleA& z){ return doubleA(new MT::Ex(MT::OUT ,(void*)&y,(void*)&z)); }
01931 inline doubleA operator%(const doubleA& y,const doubleA& z){ return doubleA(new MT::Ex(MT::MUL ,(void*)&y,(void*)&z)); }
01932 
01933 inline doubleA operator/(const doubleA& y,const doubleA& z){ return doubleA(new MT::Ex(MT::Div,(void*)&y,(void*)&z)); }
01934 inline doubleA operator+(const doubleA& y,const doubleA& z){ return doubleA(new MT::Ex(MT::PLUS,(void*)&y,(void*)&z)); }
01935 inline doubleA operator-(const doubleA& y,const doubleA& z){ return doubleA(new MT::Ex(MT::MINUS,(void*)&y,(void*)&z)); }
01936 
01937 void assign(doubleA& x,const doubleA& a);
01938 void assign(doubleA& x);
01939 template<class T> inline void MT::Array<T>::assignEx(const MT::Array<T>& a) const{ ::assign(*((doubleA*)this),*((doubleA*)&a)); }
01940 template<class T> inline void MT::Array<T>::assignEx() const{ ::assign(*((doubleA*)this)); }
01941 
01942 #endif //MT_EXPRESSIONS
01943 
01944 
01945 //===========================================================================
01946 //
01948 //
01949 
01951 void gnuplot(const doubleA& X);
01952 
01953 void write(const doubleA& X,doubleA& Y,const char* name);
01954 //void write(const doubleA& X,const char* name);
01955 
01956 //===========================================================================
01957 //
01959 //
01960 
01961 namespace MT{
01963 class Permutation:public Array<uint>{
01964 private:
01965   Array<uint> storage;
01966 
01967 public:
01969   Permutation(){ init(0); };
01970 
01971 public:
01972 
01973   void init(uint n){ resize(n); for(uint i=0;i<N;i++) elem(i)=i; }
01975   void random(uint n){ init(n); random(); }
01977   void reverse(uint n){ resize(n); for(uint i=0;i<N;i++) elem(N-1-i)=i; }
01978 
01979 public:
01980 
01981   void push(int offset=1){for(uint i=0;i<N;i++) elem(i)=((int)elem(i)+offset)%N;}
01983   void store(){ storage.operator=(*this); }
01985   void restore(){ (Array<uint>&)(*this)=storage; }
01987   void permute(uint i,uint j){ uint x=elem(i); elem(i)=elem(j); elem(j)=x; }
01989   void random(){
01990     int j,r;
01991     for(j=N-1;j>=1;j--){
01992       r=rnd(j+1);
01993       permute(r,j);
01994     }
01995   }
01997   template<class T>
01998   void permute(Array<T>& a){
01999     CHECK(N<=a.N,"array smaller than permutation ("<<a.N<<"<"<<N<<")");
02000     Array<T> b=a;
02001     for(uint i=0;i<N;i++) a(i)=b(p[i]);
02002   }
02004   template<class T>
02005   void invpermute(Array<T>& a){
02006     CHECK(N<=a.N,"array smaller than permutation ("<<a.N<<"<"<<N<<")");
02007     Array<T> b=a;
02008     for(uint i=0;i<N;i++) a(p[i])=b(i);
02009   }
02010 };
02011 }
02012 
02013 
02014 //using MT::Array;
02015 
02016 
02017 //===========================================================================
02018 //
02019 // implementations
02020 //
02021 
02022 template<class T> char MT::Array<T>::memMoveInit=-1;
02023 template<class T> int MT::Array<T>::sizeT=-1;
02024 
02025 #ifdef MT_IMPLEMENTATION
02026 #  include"array.cpp"
02027 #endif
02028 
02029 
02030 #endif
[]