/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * 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 "cmn_grad.h"

gctl::common_gradient::common_gradient()
{
    set_lcg_message(LCG_THROW);
}

gctl::common_gradient::common_gradient(size_t Ln, size_t Mn)
{
    init(Ln, Mn);
}

gctl::common_gradient::~common_gradient(){}

void gctl::common_gradient::LCG_Ax(const array<double> &x, array<double> &ax)
{
    matvec(t_, G_, x, NoTrans);
    vecmul(t_, t_, w_);
    matvec(ax, G_, t_, Trans);
    return;
}

void gctl::common_gradient::set_solver(const lcg_para &para)
{
    set_lcg_para(para);
    return;
}

void gctl::common_gradient::set_weights(const _1d_array &w)
{
    if (w.size()!= Ln_) throw std::runtime_error("[gctl::common_gradient] Invalid array size.");

    for (size_t i = 0; i < Ln_; i++)
    {
        w_[i] = 1.0/w[i];
    }
    return;
}

void gctl::common_gradient::init(size_t Ln, size_t Mn)
{
    Ln_ = Ln; 
    Mn_ = Mn;
    g_.resize(Mn_);
    B_.resize(Mn_);
    G_.resize(Ln_, Mn_);
    t_.resize(Ln_);
    gm_.resize(Ln_);
    w_.resize(Ln_, 1.0);
    x_.resize(Ln_, 1.0);
    filled_.resize(Ln_, false);
    return;
}

void gctl::common_gradient::fill_model_gradient(size_t id, const _1d_array &g)
{
    if (id >= Ln_) throw std::runtime_error("[gctl::common_gradient] Invalid index.");
    if (g.size() != Mn_) throw std::runtime_error("[gctl::common_gradient] Invalid array size.");

    G_.fill_row(id, g);
    filled_[id] = true;
    gm_[id] = g.module();
    return;
}

const gctl::_1d_array &gctl::common_gradient::get_common_gradient(bool normalized)
{
    for (size_t i = 0; i < Ln_; i++)
    {
        if (!filled_[i]) throw std::runtime_error("[gctl::common_gradient] Unfilled model gradient.");
    }
    filled_.assign(false);

    G_.normalize(RowMajor);
    matvec(B_, G_, x_, Trans);

    g_.random_float(-1.0, 1.0, RdUniform);
    g_.normalize();
    LCG_Minimize(g_, B_, LCG_CG);

    if (normalized) g_.normalize();
    else
    {
        g_.normalize();
        matvec(t_, G_, g_, NoTrans);

        double tmod = t_.module(L1);
        double gmod = 0.0;
        for (size_t i = 0; i < Ln_; i++)
        {
            gmod += gm_[i]*abs(t_[i])/tmod;
        }
        g_.scale(gmod);
    }
    return g_;
}