Line data Source code
1 : /*
2 : * File: SumOfNorm2.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on November 6, 2015, 1:23 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 <algorithm>
22 : #include <complex>
23 :
24 : #include "SumOfNorm2.h"
25 : #include "MatrixFactory.h"
26 : #include "Norm2.h"
27 :
28 2 : SumOfNorm2::SumOfNorm2(size_t k) : Norm(), m_partition_length(k) {
29 2 : m_mu = 1.0;
30 2 : m_norm2 = new Norm2();
31 2 : }
32 :
33 5 : SumOfNorm2::SumOfNorm2(double mu, size_t k) : Norm(), m_mu(mu), m_partition_length(k) {
34 5 : m_norm2 = new Norm2();
35 5 : }
36 :
37 21 : SumOfNorm2::~SumOfNorm2() {
38 7 : if (m_norm2 != NULL) {
39 7 : delete m_norm2;
40 7 : m_norm2 = NULL;
41 : }
42 14 : }
43 :
44 7 : int SumOfNorm2::call(Matrix& x, double& f) {
45 7 : f = 0.0;
46 7 : const size_t n = x.getNrows();
47 7 : if (n % m_partition_length != 0) {
48 1 : throw std::invalid_argument("Given vector cannot be partitioned");
49 : }
50 : /* number of chunks */
51 6 : const size_t n_chunks = n / m_partition_length;
52 6 : int status = ForBESUtils::STATUS_OK;
53 69 : for (size_t j = 0; j < n_chunks && status == ForBESUtils::STATUS_OK; j++) {
54 : double f_temp;
55 63 : Matrix x_part = MatrixFactory::ShallowVector(x, m_partition_length, j * m_partition_length);
56 63 : status = m_norm2->call(x_part, f_temp);
57 63 : f += f_temp;
58 63 : }
59 6 : f *= m_mu;
60 6 : return status;
61 : }
62 :
63 2 : int SumOfNorm2::callProx(Matrix& v, double gamma, Matrix& prox) {
64 2 : const size_t n = v.getNrows();
65 2 : if (n % m_partition_length != 0) {
66 1 : throw std::invalid_argument("Given vector cannot be partitioned");
67 : }
68 :
69 1 : const size_t n_chunks = n / m_partition_length;
70 1 : Matrix v_part = MatrixFactory::ShallowVector();
71 :
72 1 : int status = ForBESUtils::STATUS_OK;
73 5 : for (size_t j = 0; j < n_chunks && status == ForBESUtils::STATUS_OK; j++) {
74 : double norm_temp;
75 4 : v_part = MatrixFactory::ShallowVector(v, m_partition_length, j * m_partition_length);
76 4 : Matrix prox_part = MatrixFactory::ShallowVector(prox, m_partition_length, j * m_partition_length);
77 4 : status = m_norm2->call(v_part, norm_temp);
78 4 : double factor = 1 - (gamma * m_mu / norm_temp);
79 16 : for (size_t i = 0; i < m_partition_length; i++) {
80 12 : prox_part[i] = factor * v_part[i];
81 : }
82 4 : }
83 1 : return status;
84 : }
85 :
86 2 : int SumOfNorm2::callProx(Matrix& v, double gamma, Matrix& prox, double& f_at_prox) {
87 2 : const size_t n = v.getNrows();
88 2 : if (n % m_partition_length != 0) {
89 0 : throw std::invalid_argument("Given vector cannot be partitioned");
90 : }
91 :
92 2 : f_at_prox = 0.0;
93 2 : const size_t n_chunks = n / m_partition_length;
94 2 : Matrix v_part = MatrixFactory::ShallowVector();
95 :
96 2 : int status = ForBESUtils::STATUS_OK;
97 55 : for (size_t j = 0; j < n_chunks && status == ForBESUtils::STATUS_OK; j++) {
98 : double norm_temp;
99 53 : v_part = MatrixFactory::ShallowVector(v, m_partition_length, j * m_partition_length);
100 53 : Matrix prox_part = MatrixFactory::ShallowVector(prox, m_partition_length, j * m_partition_length);
101 53 : status = m_norm2->call(v_part, norm_temp);
102 53 : double factor = std::max<double>(0.0, 1.0 - (gamma * m_mu / norm_temp));
103 112 : for (size_t i = 0; i < m_partition_length; i++) {
104 59 : prox_part[i] = factor * v_part[i];
105 : }
106 53 : f_at_prox += std::abs(factor * norm_temp);
107 53 : }
108 2 : f_at_prox *= m_mu;
109 2 : return status;
110 : }
111 :
112 5 : FunctionOntologicalClass SumOfNorm2::category() {
113 5 : FunctionOntologicalClass meta("SumOfNorm2");
114 5 : meta.set_defines_f(true);
115 5 : meta.set_defines_conjugate(true);
116 5 : meta.set_defines_prox(true);
117 5 : meta.add_superclass(FunctionOntologyRegistry::norm());
118 5 : return meta;
119 : }
120 :
121 3 : int SumOfNorm2::dualNorm(Matrix& x, double& norm) {
122 : /* this implementation is very much similar to #call */
123 3 : norm = 0.0;
124 3 : const size_t n = x.getNrows();
125 3 : if (n % m_partition_length != 0) {
126 1 : throw std::invalid_argument("Given vector cannot be partitioned");
127 : }
128 : /* number of chunks */
129 2 : const size_t n_chunks = n / m_partition_length;
130 56 : for (size_t j = 0; j < n_chunks; j++) {
131 : double f_temp;
132 54 : Matrix x_part = MatrixFactory::ShallowVector(x, m_partition_length, j * m_partition_length);
133 54 : m_norm2->call(x_part, f_temp);
134 54 : norm = std::max(norm, f_temp);
135 54 : }
136 2 : norm /= m_mu;
137 2 : return ForBESUtils::STATUS_OK;
138 3 : }
139 :
140 :
|
