Line data Source code
1 : /*
2 : * File: IndBox.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on July 26, 2015, 5:22 PM
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 :
22 : #include "IndBox.h"
23 : #include <cmath>
24 : #include <assert.h>
25 :
26 7 : IndBox::IndBox(double& uniform_lb, double& uniform_ub) : Function() {
27 7 : m_uniform_lb = &uniform_lb;
28 7 : m_uniform_ub = &uniform_ub;
29 7 : m_lb = NULL;
30 7 : m_ub = NULL;
31 7 : }
32 :
33 6 : IndBox::IndBox(Matrix& lb, Matrix& ub) : Function() {
34 4 : if (!lb.isColumnVector()) {
35 1 : throw std::invalid_argument("LB must be a vector");
36 : }
37 3 : if (!ub.isColumnVector()) {
38 1 : throw std::invalid_argument("UB must be a vector");
39 : }
40 2 : m_lb = &lb;
41 2 : m_ub = &ub;
42 2 : m_uniform_lb = NULL;
43 2 : m_uniform_ub = NULL;
44 2 : }
45 :
46 17 : IndBox::~IndBox() {
47 17 : }
48 :
49 11 : int IndBox::call(Matrix& x, double& f) {
50 11 : if (!x.isColumnVector()) {
51 1 : throw std::invalid_argument("x must be a vector");
52 : }
53 10 : bool isInside = true;
54 10 : volatile size_t i = 0;
55 :
56 10 : if (m_uniform_lb != NULL && !isinf(-(*m_uniform_lb))) { /* there's a uniform LB and this is not -inf*/
57 15 : while (i < x.getNrows() && isInside) {
58 7 : if (m_uniform_lb != NULL) {
59 7 : isInside = isInside && (x[i] >= *m_uniform_lb);
60 : }
61 7 : i++;
62 : }
63 4 : i = 0;
64 : }
65 :
66 10 : if (m_uniform_ub != NULL && !isinf(*m_uniform_ub)) {/* there's a uniform UB and this is not +inf*/
67 14 : while (i < x.getNrows() && isInside) {
68 6 : if (m_uniform_lb != NULL) {
69 6 : isInside = isInside && x[i] <= *m_uniform_ub;
70 : }
71 6 : i++;
72 : }
73 4 : i = 0;
74 : }
75 :
76 10 : if (m_lb != NULL) {
77 132 : while (i < x.getNrows() && isInside) {
78 120 : isInside = isInside && x[i] >= m_lb->get(i);
79 120 : i++;
80 : }
81 6 : i = 0;
82 : }
83 :
84 10 : if (m_ub != NULL) {
85 102 : while (i < x.getNrows() && isInside) {
86 90 : isInside = isInside && x[i] <= m_ub->get(i);
87 90 : i++;
88 : }
89 : }
90 10 : f = isInside ? 0.0 : INFINITY;
91 10 : return ForBESUtils::STATUS_OK;
92 : }
93 :
94 23406 : int IndBox::callProx(Matrix& x, double gamma, Matrix& prox, double& f_at_prox) {
95 23406 : f_at_prox = 0.0;
96 23406 : assert(x.isColumnVector());
97 :
98 : // Step 1
99 : // prox = max(LB, x)
100 117026 : for (size_t i = 0; i < x.getNrows(); i++) {
101 93620 : if (m_lb != NULL) {
102 0 : prox[i] = std::max(m_lb->get(i), x[i]);
103 93620 : } else if (m_uniform_lb != NULL) {
104 93620 : prox[i] = std::max(*m_uniform_lb, x[i]);
105 : }
106 : }
107 :
108 : // Step 2
109 : // prox = min(UB, prox)
110 117026 : for (size_t i = 0; i < x.getNrows(); i++) {
111 93620 : if (m_ub != NULL) {
112 0 : prox[i] = std::min(m_ub->get(i), prox[i]);
113 93620 : } else if (m_uniform_ub != NULL) {
114 93620 : prox[i] = std::min(*m_uniform_ub, prox[i]);
115 : }
116 : }
117 23406 : return ForBESUtils::STATUS_OK;
118 : }
119 :
120 1 : int IndBox::callProx(Matrix& x, double gamma, Matrix& prox) {
121 : double val;
122 1 : return callProx(x, gamma, prox, val);
123 : }
124 :
125 3 : FunctionOntologicalClass IndBox::category() {
126 3 : FunctionOntologicalClass ind_box_category = FunctionOntologyRegistry::indicator();
127 3 : ind_box_category.set_defines_conjugate(true);
128 3 : return ind_box_category;
129 : }
130 :
131 1 : int IndBox::callConj(Matrix& x, double& f_star) {
132 1 : f_star = 0.0;
133 : /*
134 : * Note: either m_lb or m_uniform_lb will be non-NULL.
135 : * Check out the two constructors of this class.
136 : */
137 1 : if (m_uniform_lb == NULL) {
138 0 : for (size_t i = 0; i < x.getNrows(); i++) {
139 0 : double val = x[i];
140 0 : f_star += std::max(val * m_lb->get(i), val * m_ub->get(i));
141 : }
142 : } else {
143 6 : for (size_t i = 0; i < x.getNrows(); i++) {
144 5 : double val = x[i];
145 5 : f_star += std::max(val * (*m_uniform_lb), val * (*m_uniform_ub));
146 : }
147 : }
148 :
149 1 : return ForBESUtils::STATUS_OK;
150 12 : }
|
