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 |
}; |