LinearOperator.cpp 11.9 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
#include "DescriptorUtils.h"

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

Tobias Lasser's avatar
Tobias Lasser committed
17
18
    template <typename data_t>
    LinearOperator<data_t>::LinearOperator(const LinearOperator<data_t>& other)
19
20
        : Cloneable<LinearOperator<data_t>>(),
          _domainDescriptor{other._domainDescriptor->clone()},
Jens Petit's avatar
Jens Petit committed
21
22
23
24
25
          _rangeDescriptor{other._rangeDescriptor->clone()},
          _isLeaf{other._isLeaf},
          _isAdjoint{other._isAdjoint},
          _isComposite{other._isComposite},
          _mode{other._mode}
Tobias Lasser's avatar
Tobias Lasser committed
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
56
57
58
    {
        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
59
60
61
62
63
64
65
66
67
68
69
70
71
    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>
72
    DataContainer<data_t> LinearOperator<data_t>::apply(const DataContainer<data_t>& x) const
Tobias Lasser's avatar
Tobias Lasser committed
73
    {
74
        DataContainer<data_t> result(*_rangeDescriptor, x.getDataHandlerType());
Tobias Lasser's avatar
Tobias Lasser committed
75
76
77
        apply(x, result);
        return result;
    }
Jens Petit's avatar
Jens Petit committed
78

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

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

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

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

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

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

127
            if (_mode == CompositeMode::MULT) {
Tobias Lasser's avatar
Tobias Lasser committed
128
                // sanity check the arguments for the intended evaluation tree leaf operation
Jens Petit's avatar
Jens Petit committed
129
130
                if (_rhs->getDomainDescriptor().getNumberOfCoefficients() != x.getSize()
                    || _lhs->getRangeDescriptor().getNumberOfCoefficients() != Ax.getSize())
131
                    throw InvalidArgumentError(
Jens Petit's avatar
Jens Petit committed
132
                        "LinearOperator::apply: incorrect input/output sizes for mult leaf");
Tobias Lasser's avatar
Tobias Lasser committed
133

134
                DataContainer<data_t> temp(_rhs->getRangeDescriptor(), x.getDataHandlerType());
Tobias Lasser's avatar
Tobias Lasser committed
135
136
137
138
139
140
                _rhs->apply(x, temp);
                _lhs->apply(temp, Ax);
                return;
            }
        }

141
        throw LogicError("LinearOperator: apply called on ill-formed object");
Tobias Lasser's avatar
Tobias Lasser committed
142
143
144
    }

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

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

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

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

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

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

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

200
            if (_mode == CompositeMode::MULT) {
Tobias Lasser's avatar
Tobias Lasser committed
201
                // sanity check the arguments for the intended evaluation tree leaf operation
Jens Petit's avatar
Jens Petit committed
202
203
                if (_lhs->getRangeDescriptor().getNumberOfCoefficients() != y.getSize()
                    || _rhs->getDomainDescriptor().getNumberOfCoefficients() != Aty.getSize())
204
                    throw InvalidArgumentError(
Jens Petit's avatar
Jens Petit committed
205
                        "LinearOperator::applyAdjoint: incorrect input/output sizes for mult leaf");
Tobias Lasser's avatar
Tobias Lasser committed
206

207
                DataContainer<data_t> temp(_lhs->getDomainDescriptor(), y.getDataHandlerType());
Tobias Lasser's avatar
Tobias Lasser committed
208
209
210
211
212
213
                _lhs->applyAdjoint(y, temp);
                _rhs->applyAdjoint(temp, Aty);
                return;
            }
        }

214
        throw LogicError("LinearOperator: applyAdjoint called on ill-formed object");
Tobias Lasser's avatar
Tobias Lasser committed
215
216
217
218
219
    }

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

        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
    {
232
233
234
        if (typeid(other) != typeid(*this))
            return false;

Jens Petit's avatar
Jens Petit committed
235
236
        if (*_domainDescriptor != *other._domainDescriptor
            || *_rangeDescriptor != *other._rangeDescriptor)
Tobias Lasser's avatar
Tobias Lasser committed
237
238
            return false;

239
240
241
        if (_isLeaf ^ other._isLeaf || _isComposite ^ other._isComposite)
            return false;

Tobias Lasser's avatar
Tobias Lasser committed
242
        if (_isLeaf)
243
            return (_isAdjoint == other._isAdjoint) && (*_lhs == *other._lhs);
Tobias Lasser's avatar
Tobias Lasser committed
244
245

        if (_isComposite)
246
            return _mode == other._mode && (*_lhs == *other._lhs) && (*_rhs == *other._rhs);
Tobias Lasser's avatar
Tobias Lasser committed
247
248
249
250
251

        return true;
    }

    template <typename data_t>
Tobias Lasser's avatar
Tobias Lasser committed
252
    LinearOperator<data_t>::LinearOperator(const LinearOperator<data_t>& op, bool isAdjoint)
Jens Petit's avatar
Jens Petit committed
253
254
255
256
257
258
259
        : _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
260
261
    {
    }
Tobias Lasser's avatar
Tobias Lasser committed
262
263
264

    template <typename data_t>
    LinearOperator<data_t>::LinearOperator(const LinearOperator<data_t>& lhs,
265
                                           const LinearOperator<data_t>& rhs, CompositeMode mode)
266
267
        : _domainDescriptor{mode == CompositeMode::MULT
                                ? rhs.getDomainDescriptor().clone()
268
269
270
271
                                : bestCommon(*lhs._domainDescriptor, *rhs._domainDescriptor)},
          _rangeDescriptor{mode == CompositeMode::MULT
                               ? lhs.getRangeDescriptor().clone()
                               : bestCommon(*lhs._rangeDescriptor, *rhs._rangeDescriptor)},
Jens Petit's avatar
Jens Petit committed
272
273
274
275
          _lhs{lhs.clone()},
          _rhs{rhs.clone()},
          _isComposite{true},
          _mode{mode}
Tobias Lasser's avatar
Tobias Lasser committed
276
277
278
    {
        // sanity check the descriptors
        switch (_mode) {
279
            case CompositeMode::ADD:
280
                /// feasibility checked by bestCommon()
Tobias Lasser's avatar
Tobias Lasser committed
281
282
                break;

283
            case CompositeMode::MULT:
Tobias Lasser's avatar
Tobias Lasser committed
284
                // for multiplication, domain of _lhs should match range of _rhs
285
286
                if (_lhs->getDomainDescriptor().getNumberOfCoefficients()
                    != _rhs->getRangeDescriptor().getNumberOfCoefficients())
287
                    throw InvalidArgumentError(
Jens Petit's avatar
Jens Petit committed
288
                        "LinearOperator: composite mult domain/range mismatch");
Tobias Lasser's avatar
Tobias Lasser committed
289
290
291
                break;

            default:
292
                throw LogicError("LinearOperator: unknown composition mode");
Tobias Lasser's avatar
Tobias Lasser committed
293
294
295
296
297
298
299
300
301
302
303
        }
    }

    // ------------------------------------------
    // 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