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
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
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
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);
00203 sparse[1](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
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){
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
00395
00396
00397
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
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
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
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){
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
00735 if(binary){
00736 os.write((char*)p,sizeT*N);
00737 }else{
00738 if(nd==1){
00739
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){
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
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;
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
00943
00944
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
00958
00959
00960 UnaryOperator( - )
00961 #undef UnaryOperator
00962
00963
00964
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
01007 BinaryOperator( & )
01008 #ifndef MT_EXPRESSIONS
01009 BinaryOperator( + )
01010 BinaryOperator( - )
01011 #endif
01012
01013 BinaryOperator( / )
01014 BinaryOperator( % )
01015 #undef BinaryOperator
01016 #endif
01017
01018
01019
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
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
01066 UnaryFunction( acos )
01067 UnaryFunction( asin )
01068 UnaryFunction( atan )
01069 UnaryFunction( cos )
01070 UnaryFunction( sin )
01071 UnaryFunction( tan )
01072
01073
01074 UnaryFunction( cosh )
01075 UnaryFunction( sinh )
01076 UnaryFunction( tanh )
01077 UnaryFunction( acosh )
01078 UnaryFunction( asinh )
01079 UnaryFunction( atanh )
01080
01081
01082 UnaryFunction( exp )
01083 UnaryFunction( log )
01084 UnaryFunction( log10 )
01085
01086
01087 UnaryFunction( sqrt )
01088 UnaryFunction( cbrt )
01089
01090
01091 UnaryFunction( ceil )
01092 UnaryFunction( fabs )
01093 UnaryFunction( floor )
01094 #undef UnaryFunction
01095
01096
01097
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
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
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
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
01481
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
01560
01561
01562
01563
01564
01565 if(y.nd==2 && z.nd==2){
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
01576
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){
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
01591
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){
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){
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){
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){
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)
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
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
02015
02016
02017
02018
02019
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