gctl_optimization/lib/optimization/loss_func.cpp
2024-10-08 11:25:52 +08:00

150 lines
4.6 KiB
C++

/********************************************************
* ██████╗ ██████╗████████╗██╗
* ██╔════╝ ██╔════╝╚══██╔══╝██║
* ██║ ███╗██║ ██║ ██║
* ██║ ██║██║ ██║ ██║
* ╚██████╔╝╚██████╗ ██║ ███████╗
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
* Geophysical Computational Tools & Library (GCTL)
*
* Copyright (c) 2023 Yi Zhang (yizhang-geo@zju.edu.cn)
*
* GCTL is distributed under a dual licensing scheme. You can redistribute
* it and/or modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation, either version 2
* of the License, or (at your option) any later version. You should have
* received a copy of the GNU Lesser General Public License along with this
* program. If not, see <http://www.gnu.org/licenses/>.
*
* If the terms and conditions of the LGPL v.2. would prevent you from using
* the GCTL, please consider the option to obtain a commercial license for a
* fee. These licenses are offered by the GCTL's original author. As a rule,
* licenses are provided "as-is", unlimited in time for a one time fee. Please
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
* to include some description of your company and the realm of its activities.
* Also add information on how to contact you by electronic and paper mail.
******************************************************/
#include "loss_func.h"
gctl::loss_func::loss_func()
{
init_ = false;
eps_ = 1e-16;
loss_ = 0.0;
tnum_ = 0;
ntype_ = L2;
}
gctl::loss_func::loss_func(const array<double> &tar, norm_type_e n_type, double p, double eps)
{
init(tar, n_type, p, eps);
}
gctl::loss_func::~loss_func(){}
void gctl::loss_func::init(const array<double> &tar, norm_type_e n_type, double p, double eps)
{
if (p < 1) throw std::runtime_error("[gctl::loss_func] Invalid power number.");
if (eps <= 0) throw std::runtime_error("[gctl::loss_func] Invalid epsilon value.");
init_ = true;
tnum_ = tar.size();
diff_.resize(tnum_);
us_.resize(tnum_, 1.0);
tars_ = tar;
ntype_ = n_type;
eps_ = eps;
p_ = p;
return;
}
void gctl::loss_func::set_uncertainty(double uncer)
{
if (!init_) throw std::runtime_error("[gctl::loss_func] Not initialized.");
us_.resize(tnum_, uncer);
return;
}
void gctl::loss_func::set_uncertainty(const array<double> &uncer)
{
if (!init_) throw std::runtime_error("[gctl::loss_func] Not initialized.");
if (uncer.size() != tnum_) throw std::runtime_error("[gctl::loss_func] Invalid array size.");
us_ = uncer;
return;
}
double gctl::loss_func::evaluate(const array<double> &x, array<double> &g)
{
if (!init_) throw std::runtime_error("[gctl::loss_func] Not initialized.");
if (x.size() != tnum_) throw std::runtime_error("[gctl::loss_func] Invalid array size.");
for (size_t i = 0; i < tnum_; i++)
{
diff_[i] = (x[i] - tars_[i])/us_[i];
}
double loss = 0.0;
g.resize(tnum_);
if (ntype_ == L1)
{
for (size_t i = 0; i < tnum_; i++)
{
loss += fabs(diff_[i]);
if (diff_[i] >= 0.0) g[i] = 1.0/(us_[i]*tnum_);
else g[i] = -1.0/(us_[i]*tnum_);
}
}
else if (ntype_ == L2)
{
for (size_t i = 0; i < tnum_; i++)
{
loss += diff_[i]*diff_[i];
g[i] = 2.0*diff_[i]/(us_[i]*tnum_);
}
}
else if (ntype_ == Lp)
{
for (size_t i = 0; i < tnum_; i++)
{
loss += pow(diff_[i]*diff_[i] + eps_*eps_, 0.5*p_);
g[i] = p_*pow(diff_[i]*diff_[i] + eps_*eps_, 0.5*p_ - 1)*diff_[i]/(us_[i]*tnum_);
}
}
else throw std::runtime_error("[gctl::loss_func] Invalid measurement type.");
return loss/tnum_;
}
double gctl::loss_func::evaluate(double inp, int id)
{
double val = (inp - tars_[id])/us_[id];
if (ntype_ == L1) val = fabs(val);
if (ntype_ == L2) val = val*val;
loss_ += val;
return val/tnum_;
}
double gctl::loss_func::get_loss()
{
double l = loss_;
loss_ = 0.0;
return l;
}
double gctl::loss_func::gradient(double inp, int id)
{
double val;
if (ntype_ == L1 && val >= 0) val = 1.0/(us_[id]*tnum_);
if (ntype_ == L1 && val < 0) val = -1.0/(us_[id]*tnum_);
if (ntype_ == L2) val = 2.0*(inp - tars_[id])/(us_[id]*us_[id]*tnum_);
return val;
}