search

blitzpp / blitz

Blitz++ Multi-Dimensional Array Library for C++
https://github.com/blitzpp/blitz/wiki
Other
405 stars 84 forks source link

vtext to blitz #81

Closed slayoo closed 4 years ago

slayoo commented 5 years ago

A patch proposed on SF by Suresh kumar (@mdsuresh on github?) back in March 2003: https://sourceforge.net/p/blitz/patches/6/

@mdsuresh could you comment more on the content of the patch and on the feasibility of incorporating it into the current version of Blitz?

define BZ_DEBUG

include

include

include<blitz/array.h>

include <blitz/array/indirect.h>

include

ifdef BZ_NAMESPACES

using namespace blitz; using namespace std;

endif

namespace blitz {

inline float abs( float input) { return (float) std::fabs( (float) input); } inline double abs( double input) { return (float) std::fabs( (double) input); }

/* double log(double input) { return (double) std::log( (double) input); }

template T sqrt( T input) { return (T) std::sqrt( (long double) input); } template T log( T input) { return (T) std::log( (long double) input); }

template T pow( T a, int n) { return (T) std::pow( (long double) a, (long double) n); } template T floor( T input) { return (T) std::log( (long double) input); }

template T atan( T input) { return (T) std::atan( (long double) input); } template T max( T a, T b) { return (T) std::max( (long double) a, (long double) b); } */ using std::max; using std::atan; using std::floor; using std::sqrt; using std::log; using std::exp; using std::pow; using std::min; using std::max;

inline int pow(int a, int n) { int q = 1; for(int i=0; i < n; i++) q*=a; return q; }

}

include <blitz/alg/port.h>

include <blitz/alg/matdefault.h>

include <blitz/alg/cholesky.cc>

include <blitz/alg/eigen.cc>

include <blitz/alg/balance.cc>

include <blitz/alg/gaussj.cc>

include <blitz/alg/jacobi.cc>

include <blitz/alg/fourier.cc>

include <blitz/alg/elmhes.cc>

/*

include <blitz/alg/mprove.cc>

*/

include <blitz/alg/walsh.cc>

include <blitz/alg/fwt.cc>

include <blitz/alg/fwt_filter.cc>

include <blitz/alg/svd.cc>

include <blitz/alg/lu.cc>

include <blitz/alg/vander.cc>

struct Sample { int a, b, c;

Sample() { a =0; b =0; c =0; }

};

std::ostream& operator<<(std::ostream& out, const Sample& x) { out << x.a << "," << x.b << "," << x.c << "\n"; return out; }

BEGIN_EXPOSE_BLOCK(Sample,int) 
EXPOSE_MEMBER(a) 
EXPOSE_MEMBER(b) 
EXPOSE_MEMBER(c)
END_EXPOSE_BLOCK

//EXPOSE_STRUCTURE_C3(Sample, int, a,b,c)

int main() {

Array<double, 1> vector(8); vector = 1.0,-1.0, 1.0, -1.0, 1.0 , 1.0, 1.0, -1.0;
Array<double, 1> result = blitz::walsh(vector); // char c; Array<Sample,1> arr_test(8); arr_test = Sample(); Array<int , 1> arr_a = arr_test->b(); arr_a = 10;

cout << "Sample arr=" << arr_test << "\n";

std::cout << vector << "result" << iwalsh(result) << "\n";

// //Test<Sample, 1> q; // //Array<int, 1> arr_a = q->a(); //
//
//
// //std::vector A(3); // // A = {0, 4 , 5}; //
// //result[ indexSet(A) ] = vector; Array<complex , 2> fft_input(8,2); //
fft_input = 4,4, -3,-3, 2,2 , 0,0 , -1,-1, -2,-2, 3,3, 1,1 ;

Array<complex<double> , 1> fft_input_uno(8);

//
fft_input_uno = 4, -3, 2, 0, -1, -2, 3, 1 ; //
cout << "1 dim fft = " << blitz::fft(fft_input_uno) << "\n"; //
//
Array<double,2> a(6,6), b(6,6);

a =  delim<'['>(),  10 ,    5,     9,     4,     1,     0,
                   2 ,     0,     7,     9,     2,     7,
                   6 ,    8,     2,     1,     2,      4,
                   5 ,     4,     4,     4,     6,     9,
                   9 ,     6,     9,     8,     3,     5,
                   8 ,     8,     9,     0,     2,     4, delim<']'>();

Array<double,2> k(6,8);

k =  delim<'['>(),  10 ,    5,     9,     4,     1,    delim<';'>(),
                   2 ,     0,     7,     9,     2,     delim<';'>(),
                   6 ,    8,     2,     1,     2,      delim<';'>(),
                   5 ,     4,     4,     4,     6,     delim<';'>(),
                   9 ,     6,     9,     8,     3,     delim<';'>(),
                   8 ,     8,     9,     0,     2,     delim<';'>(),
                   delim<']'>();

cout << "k=" << k;

 b = a.copy();

//
// cout << a; //
//
///*
Array<complex , 2> fft_;

fft_.reference( blitz::fft(fft_input));

cout << fft << "\n"; // cout << "walsh" << result; //
// cout << "walsh inverse"; //
// cout << blitz::iwalsh(result); //
//
cout << "inv" << blitz::ifft(fft) << "\n";

cout << "end"; //
// //cout << result[ indexSet( A)] << "\n"; //
//
//
// // cin >> c; 

   a = blitz::inv(a);
  cout << "Matrix_Inverse=" ;
  cout << a;

//
Array<double, 2> wavelet_input(8,2), wavelet_output; wavelet_input = 4,4, -3,-3, 2,2 , 0,0 , -1,-1, -2,-2, 3,3, 1,1 ;

Array<double,1>  filter;
filter.reference( blitz::makeOnFilter("Daubechies",4, double())) ;
cout << "out";
cout << filter << "\n";

cout << blitz::downlo(wavelet_input,filter) << "\n";
cout << blitz::downhi(wavelet_input,blitz::mirrorfilter(filter)) << "\n";

wavelet_output.reference(blitz::fwt_po(wavelet_input,1,filter));

cout << wavelet_output << "\n" ;

Array<double,2> f(8,8);

Array<double,2> i_f(8,8);

f=  3,     8,     8,     8,     7,     9,     1,     3,
 8,     6,     0,    10,     6,     8,     3,     5,
 5,     1,     9,     5,     4,     9,     9,     1,
 2,     9,     7,     0,     2,     6,     6,     3,
 2,     7,     5,     0,     5,     6,    10,     3,
 1,     4,     7,    10,     7,     0,     2,     7,
 2,     0,     4,     4,     9,     6,    10,     5,
 5,     3,     9,     8,     3,     5,    10,     8;

i_f=  blitz::fwt2_po( f, 0, filter);

cout   << "fwt2_po =" <<  i_f << "\n";
cout   << "iwt2_po =" <<  blitz::iwt2_po( i_f, 0, filter);

//
Array<double,2> back = blitz::iwt_po(wavelet_output,1, filter); //
cout << back << "\n"; //
Array<double,2> S , D; Array<double,1> V;

lhs(ref(S),ref(V),ref(D)) = blitz::svd(a);

 cout << "S=" << S << "\n"
      << "V=" << V << "\n"
      << "D=" << D ;

 //cout << "pinv=" << blitz::pinv(a,0.0001);

lhs(ref(S),ref(V),ref(D)) = blitz::svd(b);

 cout << "S=" << S << "\n"
      << "V=" << V << "\n"
      << "D=" << D ;

//
cout << "pinv=" << blitz::pinv(a, 0.00001); //
Array<double,2> lu_decomp(6,6), ch_decomp(6,6), ch_mat(6,6); Array<double,1> ch_p(6); lu_decomp = a; cout << "lu = " << b; lu_decomp = blitz::luinv(lu_decomp);

cout   << "lu_decomp=" << lu_decomp;

ch_decomp  = sum(  b(blitz::tensor::i, blitz::tensor::k) 
                 * b(blitz::tensor::j, blitz::tensor::k),  tensor::k);
cout << ch_decomp << "\n";

blitz::lhs( blitz::ref(ch_mat) , ch_p) = blitz::chol( ch_decomp);

cout << "ch_decomp=" << ch_mat << "p=" << ch_p;

Array<double,2> eigen_1(6,6), eigen_2;
Array<double,1> eigen_e(6);

lhs(ref(eigen_2), ref(eigen_e) ) = blitz::tridiag(ch_decomp);

cout << "eigen_e=" << eigen_e << "\n" 
     << "eigen_v=" << eigen_2 << "\n";

lhs(ref(eigen_e), ref(eigen_2) )  = jacobi(ch_decomp);

cout << "eigen_e=" << eigen_e << "\n" 
     << "eigen_v=" << eigen_2 << "\n";
cout << blitz::balance(b);

cout << blitz::hess( ch_decomp );

Array<double,1> x(6), y(6);
x = 83 ,  50,    71,    43,    30,    19;
y = 19 ,  68,    30,    54,    15,    70;

cout << "coeff = " << blitz::vander(x,y);

//cout << "eigen =" << eigen_1;

//
// Array<double,2> fwt_mat_test(4,4), ifwt_mat_test; //
// fwt_mat_test = 8 tensor::i + 4tensor::j; // cout << fwt_mat_test << "\n"; //
// ifwt_mat_test.reference(blitz::fwt2_po(fwt_mat_test, 1, filter)); //
// cout << ifwt_mat_test; //
// fwt_mat_test.reference(blitz::iwt2_po(ifwt_mat_test, 1, filter)); // cout << fwt_mat_test; //
// Array<double,3> slice_test(5,5,2);

 slice_test = 0;
 slice_test(3, Range::all() ) = 1.0; 

 //slice_test(4, Range::all() ) 

 cout << slice_test;

 char c;
 cin >> c;
 return 0;

}

//#include <blitz/alg/mat_eigen.c>


* output from the example code:

$ ./func.exe Sample arr=8 [ 0,10,0 0,10,0 0,10,0 0,10,0 0,10,0 0,10,0 0,10,0

      0,10,0

] 8 [ 1 -1 1 -1 1 1 1 -1 ]result8 [ 1 -1 1 -1 1 1 1 -1 ] 1 dim fft = 8 [ (4,0) (5,2.41421) (-2,6) (5,0.414214) (12,0) (5,-0.414214) (-2,-6 ) (5,-2.41421) ] k=6 x 8 [ 10 5 9 4 1 0 0 0 2 0 7 9 2 0 0 0 6 8 2 1 2 0 0 0 5 4 4 4 6 0 0 0 9 6 9 8 3 0 0 0 8 8 9 0 2 0 0 0 ] 8 x 2 [ (4,0) (4,0) (5,2.41421) (5,2.41421) (-2,6) (-2,6) (5,0.414214) (5,0.414214) (12,0) (12,0) (5,-0.414214) (5,-0.414214) (-2,-6) (-2,-6) (5,-2.41421) (5,-2.41421) ]

inv8 x 2 [ (4,-0) (4,-0) (-3,-9.09344e-17) (-3,-9.09344e-17) (2,7.85098e-17) (2,7.85098e-17) (5.55112e-17,2.47954e-16) (5.55112e-17,2.47954e-16) (-1,0) (-1,0) (-2,-1.69458e-16) (-2,-1.69458e-16) (3,-7.85098e-17) (3,-7.85098e-17) (1,1.24382e-17) (1,1.24382e-17) ]

endMatrix_Inverse=6 x 6 [ 0.711543 0.343495 0.395073 0.159595 -0.871385 -0.266046 -0.519156 -0.285243 -0.143858 -0.166351 0.681635 0.16528 -0.231029 -0.111065 -0.263986 -0.0578397 0.302381 0.210513 -0.194117 -0.0904672 -0.0577573 -0.0856884 0.373898 -0.0584988 -0.66392 -0.647277 -0.624537 0.0991184 1.08865 0.173437 0.467002 0.457032 0.403807 0.0940924 -0.845184 -0.108841 ] out4 [ 0.482963 0.836516 0.224144 -0.12941 ] 4 x 2 [ -0.12941 -0.12941 1.0006 1.0006 -1.61297 -1.61297 3.57021 3.57021 ]

4 x 2 [ -4.18258 -4.18258 -1.82783 -1.82783 0.12941 0.12941 -2.60428 -2.60428 ]

8 x 2 [ -0.0490381 -0.0490381 2.04904 2.04904 1.18301 1.18301 3.28109 3.28109 -4.18258 -4.18258 -1.82783 -1.82783 0.12941 0.12941 -2.60428 -2.60428 ]

fwt2_po =8 x 8 [ 42.25 -1.07804 -2.4482 -5.64211 -1.92196 5.23085 2.00657 3.02123 2.22756 0.237139 7.57616 0.777283 2.7778 -5.57115 3.57684 -5.26554 -4.70723 0.79337 -1.54611 1.75725 4.3403 -1.09271 -3.66506 -0.575962 0.626342 3.78752 -0.744391 -2.07139 1.43301 -3.59808 0.996715 -3.69287 -2.28109 -1.65521 -1.39383 4.98205 -2.73205 3.34808 -0.353766 -2.87828 3.84687 -2.07356 -2.59479 -1.06907 3.69078 3.35256 -1.30561 3.83582 -1.80801 5.87588 -3.33253 -0.753285 -1.03678 0.61843 0.128285 4.35705 3.67524 3.23205 -4.00897 0.358253 0.68093 -0.470994 3.38397 1.38157 ]

iwt2_po =8 x 8 [ 3 8 8 8 7 9 1 3 8 6 3.72835e-11 10 6 8 3 5 5 1 9 5 4 9 9 1 2 9 7 3.1354e-11 2 6 6 3 2 7 5 3.45368e-11 5 6 10 3 1 4 7 10 7 3.25819e-11 2 7 2 3.18927e-11 4 4 9 6 10 5 5 3 9 8 3 5 10 8 ] 8 x 2 [ 4 4 -3 -3 2 2 1.00053e-12 1.00053e-12 -1 -1 -2 -2 3 3 1 1 ]

S=6 x 6 [ -0.49333 0.449302 -0.329626 0.528266 -0.10841 -0.39406 0.357912 -0.548519 -0.207353 0.397923 0.383244 -0.472031 0.190638 -0.0549546 0.674312 0.531672 -0.47196 -0.0227877 0.158921 -0.20925 -0.569756 0.336852 -0.374675 0.593705 0.606027 0.655458 -0.0151774 0.201146 0.365145 0.170541 -0.44683 -0.144306 0.262304 0.354948 0.587293 0.489709 ]

V=6 [ 2.562 0.380109 0.26599 0.0312606 0.117209 0.0868117 ] D=6 x 6 [ -0.477265 0.408348 -0.148533 0.425327 -0.533194 -0.343816 -0.352686 -0.406166 0.19654 0.334389 -0.240429 0.708796 -0.337554 -0.485705 -0.489115 0.299328 0.489507 -0.285833 -0.0565538 0.619431 0.0559518 0.378519 0.595527 0.334735 0.71363 -0.0650522 -0.38143 0.531429 -0.221452 0.0977497 0.146362 -0.210821 0.742601 0.438168 0.120782 -0.419639 ] S=6 x 6 [ -0.425327 0.343816 -0.533194 0.477265 -0.148533 -0.408348 -0.334389 -0.708796 -0.240429 0.352686 0.19654 0.406166 -0.299328 0.285832 0.489507 0.337554 -0.489115 0.485705 -0.378519 -0.334735 0.595527 0.0565538 0.0559518 -0.619431 -0.531428 -0.0977496 -0.221452 -0.71363 -0.38143 0.0650522 -0.438168 0.419639 0.120783 -0.146362 0.742601 0.210821 ]

V=6 [ 31.9892 11.5192 8.53176 0.39032 3.75954 2.63082 ] D=6 x 6 [ -0.528266 0.39406 -0.10841 0.49333 -0.329626 -0.449302 -0.397923 0.472031 0.383244 -0.357912 -0.207352 0.548519 -0.531672 0.0227878 -0.47196 -0.190637 0.674312 0.0549546 -0.336852 -0.593705 -0.374675 -0.158922 -0.569756 0.20925 -0.201146 -0.170541 0.365145 -0.606027 -0.0151774 -0.655458 -0.354948 -0.489709 0.587292 0.44683 0.262304 0.144306 ] pinv=6 x 6 [ 10 5 9 4 1 -1.1198e-06 2 -6.04665e-07 7 9 2 7 6 8 2 0.999999 2 4 5 4 4 4 6 9 9 6 9 8 3 5 8 8 9 -2.38878e-06 2 4 ] lu = 6 x 6 [ 10 5 9 4 1 0 2 0 7 9 2 7 6 8 2 1 2 4 5 4 4 4 6 9 9 6 9 8 3 5 8 8 9 0 2 4 ] 6 x 6 [ 0.711543 0.343495 0.395073 0.159595 -0.871385 -0.266046 -0.519156 -0.285243 -0.143858 -0.166351 0.681635 0.16528 -0.231029 -0.111065 -0.263986 -0.0578397 0.302381 0.210513 -0.194117 -0.0904672 -0.0577573 -0.0856884 0.373898 -0.0584988 -0.66392 -0.647277 -0.624537 0.0991184 1.08865 0.173437 0.467002 0.457032 0.403807 0.0940924 -0.845184 -0.108841 ]

6 [ 0.871385 0.681635 0.302381 0.373898 1.08865 0.845184 ] lu_decomp=6 x 6 [ 10 5 9 4 1 6.24121e-16 2 -1.35902e-15 7 9 2 7 6 8 2 1 2 4 5 4 4 4 6 9 9 6 9 8 3 5 8 8 9 -8.88178e-15 2 4 ] 6 x 6 [ 223 121 124 128 236 203 121 187 67 149 194 111 124 67 125 122 154 150 128 149 122 190 200 156 236 194 154 200 296 227 203 111 150 156 227 229 ]

ch_decomp=6 x 6 [ 223 121 124 128 236 203 8.10276 187 67 149 194 111 8.30365 -0.0256463 125 122 154 150 8.57151 7.22126 6.81358 190 200 156 15.8037 5.98657 3.06213 0.104652 296 227 13.5939 0.0773459 4.95862 1.21183 -3.58022 229 ] p=6 [ 14.9332 11.0157 7.48657 4.23766 1.00773 2.30779 ]eigen_e=6 [ 6.92124 0.152351 14.1342 72.7909 132.692 1023.31 ] eigen_v=6 x 6 [ -0.408347 -0.477265 0.148534 -0.533194 -0.343816 0.425327 0.406166 -0.352686 -0.196541 -0.240429 0.708796 0.334389 0.485706 -0.337554 0.489115 0.489507 -0.285832 0.299328 -0.619432 -0.0565536 -0.0559508 0.595527 0.334734 0.378519 0.065053 0.71363 0.38143 -0.221452 0.0977497 0.531428 0.21082 0.146362 -0.742602 0.120783 -0.419639 0.438168 ]

eigen_e=6 [ 1023.31 132.692 72.7909 6.92124 0.152352 14.1341 ] eigen_v=6 x 6 [ 0.425327 -0.343816 -0.533194 0.408348 -0.477265 -0.148533 0.334389 0.708796 -0.240429 -0.406166 -0.352686 0.19654 0.299328 -0.285832 0.489507 -0.485705 -0.337554 -0.489115 0.378519 0.334735 0.595527 0.619431 -0.0565537 0.0559517 0.531428 0.0977496 -0.221452 -0.0650523 0.71363 -0.38143 0.438168 -0.419639 0.120783 -0.210821 0.146362 0.742601 ]

6 x 6 [ 10 5 9 4 2 0 2 0 7 9 4 7 6 8 2 1 4 4 5 4 4 4 12 9 4.5 3 4.5 4 3 2.5 8 8 9 0 4 4 ] 6 x 6 [ 223 607.229 203.819 214.056 101.817 124 236 780.114 327.177 206.598 170.177 154 0.525424 154.62 114.342 31.444 37.8277 38.4746 0.542373 0.254391 -34.2875 97.7105 -52.049 24.182 0.512712 0.655794 -0.78063 -28.237 23.9455 -18.3119 0.860169 -0.172792 0.850909 -0.154698 -3.38177 10.8877 ] coeff = 6 [ -15.3477 -479.797 64.9681 700.48 -420.626 169.323 ]5 x 5 x 2 [ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 ]

slayoo commented 4 years ago

closing as in all likelihood there is too little information and too long time since the patch was proposed for it to be incorporated