Line data    Source code

       1             : /*
       2             :  * File:   IndSOC.cpp
       3             :  * Author: Lorenzo Stella
       4             :  *
       5             :  * Created on September 21, 2015, 11:08 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 "IndSOC.h"
      22             : 
      23           2 : IndSOC::IndSOC(int n) : Function() {
      24           2 :     this->m_n = n;
      25           2 : }
      26             : 
      27           4 : IndSOC::~IndSOC() {
      28           4 : }
      29             : 
      30           3 : int IndSOC::call(Matrix& x, double& f) {
      31           3 :     if (!x.isColumnVector()) {
      32           0 :         throw std::invalid_argument("x must be a vector");
      33             :     }
      34           3 :     volatile bool isInside = true;
      35           3 :     volatile size_t i = 0;
      36           3 :     size_t n = x.getNrows();
      37           3 :     double squaredNorm = 0;
      38           3 :     double t = x[n - 1];
      39             : 
      40          18 :     while (i < n - 1 && isInside) {
      41             :         double xi;
      42          12 :         xi = x[i];
      43          12 :         squaredNorm += xi*xi;
      44          12 :         i++;
      45             :     }
      46           3 :     f = (t >= sqrt(squaredNorm)) ? 0.0 : INFINITY;
      47           3 :     return ForBESUtils::STATUS_OK;
      48             : }
      49             : 
      50           0 : int IndSOC::callProx(Matrix& x, double gamma, Matrix& prox) {
      51             :     double f_at_prox;
      52           0 :     return callProx(x, gamma, prox, f_at_prox);
      53             : }
      54             : 
      55           3 : int IndSOC::callProx(Matrix& x, double gamma, Matrix& prox, double& f_at_prox) {
      56           3 :     f_at_prox = 0.0;
      57           3 :     if (!x.isColumnVector()) {
      58           0 :         throw std::invalid_argument("x must be a vector");
      59             :     }
      60           3 :     size_t i = 0;
      61           3 :     size_t n = x.getNrows();
      62           3 :     double norm = 0, xi, scal;
      63           3 :     double t = x[n-1];
      64             : 
      65          18 :     while (i < n - 1) {
      66          12 :         xi = x[i];
      67          12 :         norm += xi*xi;
      68          12 :         i++;
      69             :     }
      70             : 
      71           3 :     norm = sqrt(norm);
      72             : 
      73           3 :     if (t > norm) {
      74           1 :         prox = Matrix(x); // prox = x
      75           2 :     } else if (t < -norm) { // prox = zero vector
      76           6 :         for (size_t j = 0; j < n; j++) {
      77           5 :             prox[j] = 0.0;
      78             :         }
      79             :     } else {
      80             :         /* perform actual projection here */
      81           1 :         scal = (1 + t / norm) / 2;
      82           5 :         for (size_t j = 0; j < n - 1; j++) {
      83           4 :             xi = x[j];
      84           4 :             prox[j] = scal * xi;
      85             :         }
      86           1 :         prox[n - 1] = (norm + t) / 2.0;
      87             :     }
      88           3 :     return ForBESUtils::STATUS_OK;
      89             : }
      90             : 
      91           2 : FunctionOntologicalClass IndSOC::category() {
      92           2 :     return FunctionOntologyRegistry::indicator();
      93           3 : }