Line data Source code
1 : /*
2 : * File: OpGradient.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on September 16, 2015, 1:35 AM
6 : *
7 : * ForBES is free software: you can redistribute it and/or modify
8 : * it under the terms of the GNU Lesser General Public License as published by
9 : * the Free Software Foundation, either version 3 of the License, or
10 : * (at your option) any later version.
11 : *
12 : * ForBES is distributed in the hope that it will be useful,
13 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 : * GNU Lesser General Public License for more details.
16 : *
17 : * You should have received a copy of the GNU Lesser General Public License
18 : * along with ForBES. If not, see <http://www.gnu.org/licenses/>.
19 : */
20 :
21 : #include "OpGradient.h"
22 : #include <sstream>
23 :
24 : void call_1d(Matrix & Tx, Matrix& x, const size_t n, double alpha);
25 : void call_1d(Matrix & Tx, Matrix& x, const size_t n, double alpha, double gamma);
26 : void callAdjoint_1d(Matrix& Tstar_x, Matrix& y, const size_t n, double alpha);
27 : void callAdjoint_1d(Matrix& Tstar_x, Matrix& y, const size_t n, double alpha, double gammna);
28 :
29 0 : OpGradient::OpGradient() : m_dimension(_EMPTY_OP_DIM) {
30 0 : }
31 :
32 3 : OpGradient::OpGradient(size_t n) : LinearOperator(), m_dimension(_VECTOR_OP_DIM(n)) {
33 :
34 3 : }
35 :
36 6 : OpGradient::~OpGradient() {
37 6 : }
38 :
39 : /**
40 : * Calls OpGradient when the input is 1D
41 : * @param Tx matrix (vector) to store the result
42 : * @param x input vector x
43 : * @param n size of x
44 : */
45 4 : void call_1d(Matrix & Tx, Matrix& x, const size_t n, double alpha) {
46 190 : for (size_t i = 0; i < n - 1; i++) {
47 186 : Tx[i] = alpha * (x[i + 1] - x[i]);
48 : }
49 4 : }
50 :
51 0 : void call_1d(Matrix & Tx, Matrix& x, const size_t n, double alpha, double gamma) {
52 0 : for (size_t i = 0; i < n - 1; i++) {
53 0 : Tx[i] = gamma * Tx[i] + alpha * (x[i + 1] - x[i]);
54 : }
55 0 : }
56 :
57 4 : void callAdjoint_1d(Matrix& Tstar_x, Matrix& y, const size_t n, double alpha) {
58 4 : Tstar_x[0] = -alpha * y[0];
59 186 : for (size_t i = 1; i < n - 1; i++) {
60 182 : Tstar_x[i] = alpha * (y[i - 1] - y[i]);
61 : }
62 4 : Tstar_x[n - 1] = alpha * y[n - 2];
63 4 : }
64 :
65 0 : void callAdjoint_1d(Matrix& Tstar_x, Matrix& y, const size_t n, double alpha, double gamma) {
66 0 : Tstar_x.set(0, 0, gamma * Tstar_x.get(0, 0) - alpha * y.get(0, 0));
67 0 : for (size_t i = 1; i < n - 1; i++) {
68 0 : Tstar_x[i] = gamma * Tstar_x[i] + alpha * (y[i - 1] - y[i]);
69 : }
70 0 : Tstar_x[n - 1] = gamma * Tstar_x[n - 1] + alpha * y[n - 2];
71 0 : }
72 :
73 4 : int OpGradient::call(Matrix& y, double alpha, Matrix& x, double gamma) {
74 4 : const size_t n = x.getNrows();
75 4 : if (m_dimension.first != 0 && n != m_dimension.first) {
76 0 : std::ostringstream oss;
77 0 : oss << "[call] OpGradient operator with dimension " << m_dimension.first
78 0 : << "; argument is of incompatible dimensions " << n << "x" << x.getNcols();
79 0 : throw std::invalid_argument(oss.str().c_str());
80 : }
81 4 : if (n <= 1) {
82 0 : y = alpha*x;
83 0 : return ForBESUtils::STATUS_NUMERICAL_PROBLEMS;
84 : }
85 4 : call_1d(y, x, n, alpha);
86 4 : return ForBESUtils::STATUS_OK;
87 : }
88 :
89 4 : int OpGradient::callAdjoint(Matrix& y, double alpha, Matrix& x, double gamma) {
90 4 : callAdjoint_1d(y, x, x.getNrows() + 1, alpha);
91 4 : return ForBESUtils::STATUS_OK;
92 : }
93 :
94 26 : std::pair<size_t, size_t> OpGradient::dimensionIn() {
95 26 : return m_dimension;
96 : }
97 :
98 13 : std::pair<size_t, size_t> OpGradient::dimensionOut() {
99 13 : std::pair<size_t, size_t> dims;
100 13 : dims.second = static_cast<size_t> (1);
101 13 : if (m_dimension.first == 0) {
102 0 : dims.first = static_cast<size_t> (0);
103 : } else {
104 13 : dims.first = static_cast<size_t> (dimensionIn().first - 1);
105 : }
106 13 : return dims;
107 : }
108 :
109 1 : bool OpGradient::isSelfAdjoint() {
110 1 : return false;
111 : }
112 :
113 :
114 :
115 :
116 :
117 :
|
