00001
00002
00003
00004
00007 #ifndef MT_array_h
00008 #define MT_array_h
00009
00010
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
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
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:
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){
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
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
00329
00330
00331
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
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
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
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
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){
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
00657 if(binary){
00658 os.write((char*)p,sizeT*N);
00659 }else{
00660 if(nd==1){
00661
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){
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){
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
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
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;
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
00828
00829
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
00843
00844
00845 UnaryOperator( - )
00846 #undef UnaryOperator
00847
00848
00849
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
00892 BinaryOperator( & )
00893 #ifndef MT_EXPRESSIONS
00894 BinaryOperator( + )
00895 BinaryOperator( - )
00896 #endif
00897
00898 BinaryOperator( / )
00899 BinaryOperator( % )
00900 #undef BinaryOperator
00901 #endif
00902
00903
00904
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
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
00951 UnaryFunction( acos )
00952 UnaryFunction( asin )
00953 UnaryFunction( atan )
00954 UnaryFunction( cos )
00955 UnaryFunction( sin )
00956 UnaryFunction( tan )
00957
00958
00959 UnaryFunction( cosh )
00960 UnaryFunction( sinh )
00961 UnaryFunction( tanh )
00962 UnaryFunction( acosh )
00963 UnaryFunction( asinh )
00964 UnaryFunction( atanh )
00965
00966
00967 UnaryFunction( exp )
00968 UnaryFunction( log )
00969 UnaryFunction( log10 )
00970
00971
00972 UnaryFunction( sqrt )
00973 UnaryFunction( cbrt )
00974
00975
00976 UnaryFunction( ceil )
00977 UnaryFunction( fabs )
00978 UnaryFunction( floor )
00979 #undef UnaryFunction
00980
00981
00982
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
01203
01204
01205
01206
01207
01208 if(y.nd==2 && z.nd==2){
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
01216
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){
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
01231
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){
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){
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){
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){
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)
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
01669
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
01742
01743
01744
01745
01746
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