LinearOperator.cpp 11.5 KB
Newer Older
Tobias Lasser's avatar
Tobias Lasser committed
1
2
3
#include "LinearOperator.h"

#include <stdexcept>
Tobias Lasser's avatar
Tobias Lasser committed
4
#include <typeinfo>
Tobias Lasser's avatar
Tobias Lasser committed
5
6
7
8

namespace elsa
{
    template <typename data_t>
Jens Petit's avatar
Jens Petit committed
9
10
    LinearOperator<data_t>::LinearOperator(const DataDescriptor& domainDescriptor,
                                           const DataDescriptor& rangeDescriptor)
Tobias Lasser's avatar
Tobias Lasser committed
11
        : _domainDescriptor{domainDescriptor.clone()}, _rangeDescriptor{rangeDescriptor.clone()}
Jens Petit's avatar
Jens Petit committed
12
13
    {
    }
Tobias Lasser's avatar
Tobias Lasser committed
14

Tobias Lasser's avatar
Tobias Lasser committed
15
16
    template <typename data_t>
    LinearOperator<data_t>::LinearOperator(const LinearOperator<data_t>& other)
Jens Petit's avatar
Jens Petit committed
17
18
19
20
21
22
        : _domainDescriptor{other._domainDescriptor->clone()},
          _rangeDescriptor{other._rangeDescriptor->clone()},
          _isLeaf{other._isLeaf},
          _isAdjoint{other._isAdjoint},
          _isComposite{other._isComposite},
          _mode{other._mode}
Tobias Lasser's avatar
Tobias Lasser committed
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
    {
        if (_isLeaf)
            _lhs = other._lhs->clone();

        if (_isComposite) {
            _lhs = other._lhs->clone();
            _rhs = other._rhs->clone();
        }
    }

    template <typename data_t>
    LinearOperator<data_t>& LinearOperator<data_t>::operator=(const LinearOperator<data_t>& other)
    {
        if (*this != other) {
            _domainDescriptor = other._domainDescriptor->clone();
            _rangeDescriptor = other._rangeDescriptor->clone();
            _isLeaf = other._isLeaf;
            _isAdjoint = other._isAdjoint;
            _isComposite = other._isComposite;
            _mode = other._mode;

            if (_isLeaf)
                _lhs = other._lhs->clone();

            if (_isComposite) {
                _lhs = other._lhs->clone();
                _rhs = other._rhs->clone();
            }
        }

        return *this;
    }

Tobias Lasser's avatar
Tobias Lasser committed
56
57
58
59
60
61
62
63
64
65
66
67
68
    template <typename data_t>
    const DataDescriptor& LinearOperator<data_t>::getDomainDescriptor() const
    {
        return *_domainDescriptor;
    }

    template <typename data_t>
    const DataDescriptor& LinearOperator<data_t>::getRangeDescriptor() const
    {
        return *_rangeDescriptor;
    }

    template <typename data_t>
69
    DataContainer<data_t> LinearOperator<data_t>::apply(const DataContainer<data_t>& x) const
Tobias Lasser's avatar
Tobias Lasser committed
70
71
72
73
74
    {
        DataContainer<data_t> result(*_rangeDescriptor);
        apply(x, result);
        return result;
    }
Jens Petit's avatar
Jens Petit committed
75

Tobias Lasser's avatar
Tobias Lasser committed
76
    template <typename data_t>
Jens Petit's avatar
Jens Petit committed
77
78
    void LinearOperator<data_t>::apply(const DataContainer<data_t>& x,
                                       DataContainer<data_t>& Ax) const
Tobias Lasser's avatar
Tobias Lasser committed
79
80
81
82
83
    {
        _apply(x, Ax);
    }

    template <typename data_t>
Jens Petit's avatar
Jens Petit committed
84
85
    void LinearOperator<data_t>::_apply(const DataContainer<data_t>& x,
                                        DataContainer<data_t>& Ax) const
Tobias Lasser's avatar
Tobias Lasser committed
86
    {
Tobias Lasser's avatar
Tobias Lasser committed
87
88
89
        if (_isLeaf) {
            if (_isAdjoint) {
                // sanity check the arguments for the intended evaluation tree leaf operation
Jens Petit's avatar
Jens Petit committed
90
91
92
93
                if (_lhs->getRangeDescriptor().getNumberOfCoefficients() != x.getSize()
                    || _lhs->getDomainDescriptor().getNumberOfCoefficients() != Ax.getSize())
                    throw std::invalid_argument(
                        "LinearOperator::apply: incorrect input/output sizes for adjoint leaf");
Tobias Lasser's avatar
Tobias Lasser committed
94
95
96

                _lhs->applyAdjoint(x, Ax);
                return;
Jens Petit's avatar
Jens Petit committed
97
            } else {
Tobias Lasser's avatar
Tobias Lasser committed
98
                // sanity check the arguments for the intended evaluation tree leaf operation
Jens Petit's avatar
Jens Petit committed
99
100
101
102
                if (_lhs->getDomainDescriptor().getNumberOfCoefficients() != x.getSize()
                    || _lhs->getRangeDescriptor().getNumberOfCoefficients() != Ax.getSize())
                    throw std::invalid_argument(
                        "LinearOperator::apply: incorrect input/output sizes for leaf");
Tobias Lasser's avatar
Tobias Lasser committed
103
104
105
106

                _lhs->apply(x, Ax);
                return;
            }
Tobias Lasser's avatar
Tobias Lasser committed
107
108
109
110
111
        }

        if (_isComposite) {
            if (_mode == compositeMode::add) {
                // sanity check the arguments for the intended evaluation tree leaf operation
Jens Petit's avatar
Jens Petit committed
112
113
114
115
116
117
                if (_rhs->getDomainDescriptor().getNumberOfCoefficients() != x.getSize()
                    || _rhs->getRangeDescriptor().getNumberOfCoefficients() != Ax.getSize()
                    || _lhs->getDomainDescriptor().getNumberOfCoefficients() != x.getSize()
                    || _lhs->getRangeDescriptor().getNumberOfCoefficients() != Ax.getSize())
                    throw std::invalid_argument(
                        "LinearOperator::apply: incorrect input/output sizes for add leaf");
Tobias Lasser's avatar
Tobias Lasser committed
118
119
120
121
122
123
124
125

                _rhs->apply(x, Ax);
                Ax += _lhs->apply(x);
                return;
            }

            if (_mode == compositeMode::mult) {
                // sanity check the arguments for the intended evaluation tree leaf operation
Jens Petit's avatar
Jens Petit committed
126
127
128
129
                if (_rhs->getDomainDescriptor().getNumberOfCoefficients() != x.getSize()
                    || _lhs->getRangeDescriptor().getNumberOfCoefficients() != Ax.getSize())
                    throw std::invalid_argument(
                        "LinearOperator::apply: incorrect input/output sizes for mult leaf");
Tobias Lasser's avatar
Tobias Lasser committed
130
131
132
133
134
135
136
137
138
139
140
141

                DataContainer<data_t> temp(_rhs->getRangeDescriptor());
                _rhs->apply(x, temp);
                _lhs->apply(temp, Ax);
                return;
            }
        }

        throw std::logic_error("LinearOperator: apply called on ill-formed object");
    }

    template <typename data_t>
142
    DataContainer<data_t> LinearOperator<data_t>::applyAdjoint(const DataContainer<data_t>& y) const
Tobias Lasser's avatar
Tobias Lasser committed
143
144
145
146
147
    {
        DataContainer<data_t> result(*_domainDescriptor);
        applyAdjoint(y, result);
        return result;
    }
Jens Petit's avatar
Jens Petit committed
148

Tobias Lasser's avatar
Tobias Lasser committed
149
    template <typename data_t>
Jens Petit's avatar
Jens Petit committed
150
151
    void LinearOperator<data_t>::applyAdjoint(const DataContainer<data_t>& y,
                                              DataContainer<data_t>& Aty) const
Tobias Lasser's avatar
Tobias Lasser committed
152
153
154
155
156
    {
        _applyAdjoint(y, Aty);
    }

    template <typename data_t>
Jens Petit's avatar
Jens Petit committed
157
158
    void LinearOperator<data_t>::_applyAdjoint(const DataContainer<data_t>& y,
                                               DataContainer<data_t>& Aty) const
Tobias Lasser's avatar
Tobias Lasser committed
159
    {
Tobias Lasser's avatar
Tobias Lasser committed
160
161
162
        if (_isLeaf) {
            if (_isAdjoint) {
                // sanity check the arguments for the intended evaluation tree leaf operation
Jens Petit's avatar
Jens Petit committed
163
164
165
166
                if (_lhs->getDomainDescriptor().getNumberOfCoefficients() != y.getSize()
                    || _lhs->getRangeDescriptor().getNumberOfCoefficients() != Aty.getSize())
                    throw std::invalid_argument("LinearOperator::applyAdjoint: incorrect "
                                                "input/output sizes for adjoint leaf");
Tobias Lasser's avatar
Tobias Lasser committed
167
168
169

                _lhs->apply(y, Aty);
                return;
Jens Petit's avatar
Jens Petit committed
170
            } else {
Tobias Lasser's avatar
Tobias Lasser committed
171
                // sanity check the arguments for the intended evaluation tree leaf operation
Jens Petit's avatar
Jens Petit committed
172
173
174
175
                if (_lhs->getRangeDescriptor().getNumberOfCoefficients() != y.getSize()
                    || _lhs->getDomainDescriptor().getNumberOfCoefficients() != Aty.getSize())
                    throw std::invalid_argument(
                        "LinearOperator::applyAdjoint: incorrect input/output sizes for leaf");
Tobias Lasser's avatar
Tobias Lasser committed
176
177
178
179

                _lhs->applyAdjoint(y, Aty);
                return;
            }
Tobias Lasser's avatar
Tobias Lasser committed
180
181
182
183
184
        }

        if (_isComposite) {
            if (_mode == compositeMode::add) {
                // sanity check the arguments for the intended evaluation tree leaf operation
Jens Petit's avatar
Jens Petit committed
185
186
187
188
189
190
                if (_rhs->getRangeDescriptor().getNumberOfCoefficients() != y.getSize()
                    || _rhs->getDomainDescriptor().getNumberOfCoefficients() != Aty.getSize()
                    || _lhs->getRangeDescriptor().getNumberOfCoefficients() != y.getSize()
                    || _lhs->getDomainDescriptor().getNumberOfCoefficients() != Aty.getSize())
                    throw std::invalid_argument(
                        "LinearOperator::applyAdjoint: incorrect input/output sizes for add leaf");
Tobias Lasser's avatar
Tobias Lasser committed
191
192
193
194
195
196
197
198

                _rhs->applyAdjoint(y, Aty);
                Aty += _lhs->applyAdjoint(y);
                return;
            }

            if (_mode == compositeMode::mult) {
                // sanity check the arguments for the intended evaluation tree leaf operation
Jens Petit's avatar
Jens Petit committed
199
200
                if (_lhs->getRangeDescriptor().getNumberOfCoefficients() != y.getSize()
                    || _rhs->getDomainDescriptor().getNumberOfCoefficients() != Aty.getSize())
Tobias Lasser's avatar
Tobias Lasser committed
201
                    throw std::invalid_argument(
Jens Petit's avatar
Jens Petit committed
202
                        "LinearOperator::applyAdjoint: incorrect input/output sizes for mult leaf");
Tobias Lasser's avatar
Tobias Lasser committed
203
204
205
206
207
208
209
210
211
212
213
214
215
216

                DataContainer<data_t> temp(_lhs->getDomainDescriptor());
                _lhs->applyAdjoint(y, temp);
                _rhs->applyAdjoint(temp, Aty);
                return;
            }
        }

        throw std::logic_error("LinearOperator: applyAdjoint called on ill-formed object");
    }

    template <typename data_t>
    LinearOperator<data_t>* LinearOperator<data_t>::cloneImpl() const
    {
Tobias Lasser's avatar
Tobias Lasser committed
217
218
        if (_isLeaf)
            return new LinearOperator<data_t>(*_lhs, _isAdjoint);
Tobias Lasser's avatar
Tobias Lasser committed
219
220
221
222
223
224
225
226
227
228

        if (_isComposite)
            return new LinearOperator<data_t>(*_lhs, *_rhs, _mode);

        return new LinearOperator<data_t>(*_domainDescriptor, *_rangeDescriptor);
    }

    template <typename data_t>
    bool LinearOperator<data_t>::isEqual(const LinearOperator<data_t>& other) const
    {
Jens Petit's avatar
Jens Petit committed
229
230
        if (*_domainDescriptor != *other._domainDescriptor
            || *_rangeDescriptor != *other._rangeDescriptor)
Tobias Lasser's avatar
Tobias Lasser committed
231
232
            return false;

Tobias Lasser's avatar
Tobias Lasser committed
233
234
        if (_isLeaf)
            return other._isLeaf && (_isAdjoint == other._isAdjoint) && (*_lhs == *other._lhs);
Tobias Lasser's avatar
Tobias Lasser committed
235
236

        if (_isComposite)
Jens Petit's avatar
Jens Petit committed
237
238
            return other._isComposite && _mode == other._mode && (*_lhs == *other._lhs)
                   && (*_rhs == *other._rhs);
Tobias Lasser's avatar
Tobias Lasser committed
239
240
241
242
243

        return true;
    }

    template <typename data_t>
Tobias Lasser's avatar
Tobias Lasser committed
244
    LinearOperator<data_t>::LinearOperator(const LinearOperator<data_t>& op, bool isAdjoint)
Jens Petit's avatar
Jens Petit committed
245
246
247
248
249
250
251
        : _domainDescriptor{(isAdjoint) ? op.getRangeDescriptor().clone()
                                        : op.getDomainDescriptor().clone()},
          _rangeDescriptor{(isAdjoint) ? op.getDomainDescriptor().clone()
                                       : op.getRangeDescriptor().clone()},
          _lhs{op.clone()},
          _isLeaf{true},
          _isAdjoint{isAdjoint}
Tobias Lasser's avatar
Tobias Lasser committed
252
253
    {
    }
Tobias Lasser's avatar
Tobias Lasser committed
254
255
256

    template <typename data_t>
    LinearOperator<data_t>::LinearOperator(const LinearOperator<data_t>& lhs,
Jens Petit's avatar
Jens Petit committed
257
258
259
260
261
262
263
                                           const LinearOperator<data_t>& rhs, compositeMode mode)
        : _domainDescriptor{rhs.getDomainDescriptor().clone()},
          _rangeDescriptor{lhs.getRangeDescriptor().clone()},
          _lhs{lhs.clone()},
          _rhs{rhs.clone()},
          _isComposite{true},
          _mode{mode}
Tobias Lasser's avatar
Tobias Lasser committed
264
265
266
267
268
    {
        // sanity check the descriptors
        switch (_mode) {
            case compositeMode::add:
                // for addition, both domains and ranges should match
Jens Petit's avatar
Jens Petit committed
269
270
271
272
                if (_lhs->getDomainDescriptor() != _rhs->getDomainDescriptor()
                    || _lhs->getRangeDescriptor() != _rhs->getRangeDescriptor())
                    throw std::invalid_argument(
                        "LinearOperator: composite add domain/range mismatch");
Tobias Lasser's avatar
Tobias Lasser committed
273
274
275
276
277
                break;

            case compositeMode::mult:
                // for multiplication, domain of _lhs should match range of _rhs
                if (_lhs->getDomainDescriptor() != _rhs->getRangeDescriptor())
Jens Petit's avatar
Jens Petit committed
278
279
                    throw std::invalid_argument(
                        "LinearOperator: composite mult domain/range mismatch");
Tobias Lasser's avatar
Tobias Lasser committed
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
                break;

            default:
                throw std::logic_error("LinearOperator: unknown composition mode");
        }
    }

    // ------------------------------------------
    // explicit template instantiation
    template class LinearOperator<float>;
    template class LinearOperator<std::complex<float>>;
    template class LinearOperator<double>;
    template class LinearOperator<std::complex<double>>;

} // namespace elsa