Rev Author Line No. Line
178 kaklik 1 namespace mimas {
2  
3 // Convolution of multi-arrays.
4 // Jan Wedekind (jan@wedesoft.de)
5 // Fri Jul 15 21:37:32 UTC 2005
6  
7 template< typename _T >
8 struct correlate_element
9 {
10 typedef correlate_element< _T > type;
11 typedef boost::detail::multi_array::size_type size_type;
12 typedef boost::detail::multi_array::index index_type;
13 void operator()( _T *result, const _T *a, const _T *b,
14 const size_type *, const size_type *, const size_type *,
15 const index_type *, const index_type *,
16 const index_type * ) {
17 *result += *a * *b;
18 };
19 };
20  
21 template< typename _T, int _N >
22 struct correlate_array
23 {
24 typedef correlate_array< _T, _N > type;
25 typedef boost::detail::multi_array::size_type size_type;
26 typedef boost::detail::multi_array::index index_type;
27 void operator()( _T *result, const _T *a, const _T *b,
28 const size_type *sizeR, const size_type *sizeA,
29 const size_type *sizeB,
30 const index_type *stepR, const index_type *stepA,
31 const index_type *stepB );
32 typename boost::mpl::eval_if< boost::mpl::bool_< _N == 1 >,
33 correlate_element< _T >,
34 correlate_array< _T, _N - 1 > >::type f;
35 };
36  
37 template< typename _T, int _N >
38 void correlate_array< _T, _N >::operator()
39 ( _T *result, const _T *a, const _T *b,
40 const size_type *sizeR, const size_type *sizeA, const size_type *sizeB,
41 const index_type *stepR, const index_type *stepA, const index_type *stepB )
42 {
43 size_type
44 curSizeA = *sizeA,
45 curSizeB = *sizeB,
46 curSizeR = *sizeR;
47 index_type
48 curStepA = *stepA,
49 curStepB = *stepB,
50 curStepR = *stepR;
51 int
52 fringe1 = ( ( curSizeA + curSizeB - 1 ) - curSizeR ) / 2,
53 fringe2 = ( ( curSizeA + curSizeB ) - curSizeR ) / 2;
54 assert( fringe1 >= 0 );
55 assert( fringe2 >= 0 );
56 for ( int j=0; j<fringe1; j++ )
57 for ( int k=fringe1-j; k<(signed)curSizeB; k++ ) {
58 int j_plus_k = j - fringe1 + k;
59 assert( j_plus_k >= 0 );
60 assert( j_plus_k < (signed)curSizeR );
61 f( result + j_plus_k * curStepR,
62 a + j * curStepA,
63 b + k * curStepB,
64 sizeR + 1, sizeA + 1, sizeB + 1,
65 stepR + 1, stepA + 1, stepB + 1 );
66 };
67 for ( int j=fringe1; j<(signed)curSizeA - fringe2; j++ )
68 for ( int k=0; k<(signed)curSizeB; k++ ) {
69 int j_plus_k = j - fringe1 + k;
70 assert( j_plus_k >= 0 );
71 assert( j_plus_k < (signed)curSizeR );
72 f( result + j_plus_k * curStepR,
73 a + j * curStepA,
74 b + k * curStepB,
75 sizeR + 1, sizeA + 1, sizeB + 1,
76 stepR + 1, stepA + 1, stepB + 1 );
77 };
78 for ( int j=(signed)curSizeA - fringe2; j<(signed)curSizeA; j++ )
79 for ( int k=0; k<(signed)curSizeB - j + (signed)curSizeA - fringe2 - 1;
80 k++ ){
81 int j_plus_k = j - fringe1 + k;
82 assert( j_plus_k >= 0 );
83 assert( j_plus_k < (signed)curSizeR );
84 f( result + j_plus_k * curStepR,
85 a + j * curStepA,
86 b + k * curStepB,
87 sizeR + 1, sizeA + 1, sizeB + 1,
88 stepR + 1, stepA + 1, stepB + 1 );
89 };
90 };
91  
92 template< class _Array1, class _Array2 >
93 boost::multi_array< typename _Array1::element, _Array1::dimensionality >
94 correlate( const _Array1 &x, const _Array2 &y )
95 {
96 boost::array< std::size_t, _Array1::dimensionality > shape;
97 // for ( int i=0; i<(signed)_Array1::dimensionality; i++ )
98 // shape[ i ] = x.shape()[ i ] + y.shape()[ i ] - 1;
99 std::copy( x.shape(), x.shape() + _Array1::dimensionality, shape.begin() );
100 boost::multi_array< typename _Array1::element, _Array1::dimensionality >
101 retVal( shape );
102  
103 typedef boost::detail::multi_array::size_type size_type;
104 typedef boost::detail::multi_array::index index_type;
105  
106 std::multimap< size_type, int > sizes;
107 for ( int i=0; i<(signed)_Array1::dimensionality; i++ )
108 sizes.insert( std::make_pair( y.shape()[i], i ) );
109  
110 size_type
111 retValShape[ _Array1::dimensionality ],
112 xShape[ _Array1::dimensionality ],
113 yShapeSorted[ _Array1::dimensionality ];
114 index_type
115 retValStrides[ _Array1::dimensionality ],
116 xStrides[ _Array1::dimensionality ],
117 yStrides[ _Array1::dimensionality ];
118  
119 // Sort the loops, such that the inner loops for the filter 'y' have
120 // the highest number of cycles.
121 int j=0;
122 for ( std::multimap< size_type, int >::const_iterator i=sizes.begin();
123 i != sizes.end(); i++, j++ ) {
124 int index = i->second;
125 retValShape [ j ] = retVal.shape()[ index ];
126 xShape [ j ] = x .shape()[ index ];
127 yShapeSorted[ j ] = y .shape()[ index ];
128 retValStrides[ j ] = retVal.strides()[ index ];
129 xStrides [ j ] = x .strides()[ index ];
130 yStrides [ j ] = y .strides()[ index ];
131 };
132  
133 correlate_array< typename _Array1::element, _Array1::dimensionality > f;
134 f( retVal.origin(), x.origin(), y.origin(),
135 &retValShape[0], &xShape[0], &yShapeSorted[0],
136 &retValStrides[0], &xStrides[0], &yStrides[0] );
137 return retVal;
138 };
139  
140 #ifdef HAVE_LIBLAPACK
141 template< typename T >
142 std::vector< boost::multi_array< T, 2 > > separate
143 ( const boost::const_multi_array_ref< T, 2 > &array )
144 {
145 boost::numeric::ublas::matrix< T > A( array.shape()[1], array.shape()[0] );
146 std::copy( array.data(), array.data() + array.num_elements(),
147 A.data().begin() );
148 boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > u, v;
149 boost::numeric::ublas::banded_matrix< T > d( gesdd< T >( A, &u, &v ) );
150 boost::numeric::ublas::vector< T >
151 s( boost::numeric::ublas::matrix_column< boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > >
152 ( u, 0 ) ),
153 t( boost::numeric::ublas::matrix_row< boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > >
154 ( v, 0 ) );
155 std::vector< boost::multi_array< T, 2 > > retval;
156 boost::multi_array< T, 2 > sa( boost::extents[3][1] );
157 std::copy( s.data().begin(), s.data().begin() + s.size(), sa.data() );
158 retval.push_back( sa * d( 0, 0 ) );
159 boost::multi_array< T, 2 > ta( boost::extents[1][3] );
160 std::copy( t.data().begin(), t.data().begin() + t.size(), ta.data() );
161 retval.push_back( ta );
162 return retval;
163 }
164 #endif
165  
166 };