Line data Source code
1 : /*
2 : * File: QuadOverAffine.cpp
3 : * Author: Pantelis Sopasakis
4 : *
5 : * Created on July 24, 2015, 4:55 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 : #include "QuadOverAffine.h"
22 : #include "LDLFactorization.h"
23 :
24 : void checkConstructorArguments(const Matrix& Q, const Matrix& q, const Matrix& A, const Matrix& b);
25 :
26 0 : QuadOverAffine::QuadOverAffine() : Function() {
27 0 : m_A = NULL;
28 0 : m_F = NULL;
29 0 : m_Fsolver = NULL;
30 0 : m_Q = NULL;
31 0 : m_q = NULL;
32 0 : m_b = NULL;
33 0 : m_sigma = NULL;
34 0 : }
35 :
36 11 : QuadOverAffine::~QuadOverAffine() {
37 4 : if (m_Fsolver != NULL) {
38 4 : delete m_Fsolver;
39 : }
40 4 : if (m_F != NULL) {
41 4 : delete m_F;
42 : }
43 4 : if (m_sigma != NULL) {
44 4 : delete m_sigma;
45 : }
46 7 : }
47 :
48 4 : void checkConstructorArguments(const Matrix& Q, const Matrix& q, const Matrix& A, const Matrix& b) {
49 4 : if (Q.getNrows() != Q.getNcols()) {
50 0 : throw std::invalid_argument("Matrix Q is not square");
51 : }
52 4 : if (!q.isColumnVector()) {
53 0 : throw std::invalid_argument("q is not a column vector");
54 : }
55 4 : if (!b.isColumnVector()) {
56 0 : throw std::invalid_argument("b is not a column vector");
57 : }
58 4 : if (Q.getNcols() != q.getNrows()) {
59 0 : throw std::invalid_argument("Q and q have incompatible dimensions");
60 : }
61 4 : if (Q.getNcols() != A.getNcols()) {
62 0 : throw std::invalid_argument("Q and A have incompatible dimensions");
63 : }
64 4 : if (A.getNrows() != b.getNrows()) { // A and b must have the same number of rows
65 0 : throw std::invalid_argument("A and b have incompatible dimensions");
66 : }
67 4 : }
68 :
69 4 : QuadOverAffine::QuadOverAffine(Matrix& Q, Matrix& q, Matrix& A, Matrix& b) {
70 4 : checkConstructorArguments(Q, q, A, b);
71 :
72 4 : m_F = NULL;
73 4 : m_Fsolver = NULL;
74 4 : m_sigma = NULL;
75 :
76 4 : this->m_Q = &Q;
77 4 : this->m_q = &q;
78 4 : this->m_A = &A;
79 4 : this->m_b = &b;
80 4 : size_t n = Q.getNrows();
81 4 : size_t s = A.getNrows();
82 4 : size_t nF = n + s;
83 4 : if (Q.getType() == Matrix::MATRIX_DENSE) {
84 4 : m_F = new Matrix(nF, nF, Matrix::MATRIX_DENSE);
85 : /*
86 : * F = [Q * ; * *]
87 : */
88 36 : for (size_t i = 0; i < n; i++) {
89 288 : for (size_t j = 0; j < n; j++) {
90 256 : m_F->set(i, j, Q.get(i, j));
91 : }
92 : }
93 : /*
94 : * F = [Q A' ; A *]
95 : */
96 20 : for (size_t i = 0; i < s; i++) {
97 144 : for (size_t j = 0; j < n; j++) {
98 128 : m_F->set(i + n, j, A.get(i, j));
99 128 : m_F->set(j, i + n, A.get(i, j));
100 : }
101 : }
102 : }
103 4 : if (m_F != NULL) {
104 4 : m_Fsolver = new LDLFactorization(*m_F);
105 4 : int status = m_Fsolver -> factorize();
106 4 : if (ForBESUtils::STATUS_OK != status) {
107 0 : throw std::invalid_argument("LDL factorization failed for matrix F = [Q A'; A 0] (dense) - invalid arguments Q and A");
108 : }
109 : }
110 4 : m_sigma = new Matrix(nF, 1, Matrix::MATRIX_DENSE);
111 20 : for (size_t i = 0; i < s; i++) {
112 16 : m_sigma->set(i + n, 0, b.get(i, 0));
113 : }
114 4 : }
115 :
116 0 : int QuadOverAffine::callConj(Matrix& y, double& f_star) {
117 0 : Matrix grad(y.getNrows(), y.getNcols());
118 0 : return callConj(y, f_star, grad);
119 : }
120 :
121 3 : int QuadOverAffine::callConj(Matrix& y, double& f_star, Matrix& grad) {
122 : /* Update sigma(x) */
123 27 : for (size_t i = 0; i < m_Q->getNrows(); i++) {
124 24 : m_sigma->getData()[i] = y[i] - m_q->get(i);
125 : }
126 : /* Solve F*grad = sigma */
127 3 : int status = m_Fsolver->solve(*m_sigma, grad);
128 : /* Take the first n elements of grad */
129 3 : grad.reshape(m_Q->getNrows(), 1);
130 : /* f_star = grad' * Q * grad / 2.0 */
131 3 : f_star = m_Q->quad(grad);
132 27 : for (size_t i = 0; i < grad.getNrows(); i++) { /* Dot product */
133 24 : f_star += grad[i] * (m_q->get(i) - y[i]);
134 : }
135 3 : f_star = -f_star;
136 3 : return status;
137 : }
138 :
139 3 : FunctionOntologicalClass QuadOverAffine::category() {
140 3 : FunctionOntologicalClass meta("QuadraticOverAffine");
141 3 : meta.set_defines_conjugate(true);
142 3 : meta.set_defines_conjugate_grad(true);
143 3 : meta.add_superclass(FunctionOntologyRegistry::conj_quadratic());
144 3 : return meta;
145 6 : }
146 :
147 :
|
