/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2022  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 "../lib/optimization.h"

using namespace gctl;

class cfg : public lgd_solver
{
private:
    std::vector<gaussian_para2d> gs1_, gs2_;
    array<double> g1_, g2_;

    common_gradient cmg_;

    array<double> x_;
    array<double> xm_;
    array<double> xs_;

public:
    cfg(/* args */){}
    ~cfg(){}

    virtual double LGD_Evaluate(const array<double> &x, array<double> &g);
    void read_gaussians(std::string file, std::vector<gaussian_para2d> &gs);
    void routine(int argc, char *argv[]);
};

double cfg::LGD_Evaluate(const array<double> &x, array<double> &g)
{
    double f, fx = 0.0;
	g1_.assign(0.0);
	g2_.assign(0.0);

    f = 0.0;
	for (int i = 0; i < gs1_.size(); i++)
	{
		f += -1e+3*gctl::gaussian_dist2d(x[0], x[1], gs1_[i]);
		g1_[0] += -1e+3*gctl::gaussian_dist2d(x[0], x[1], gs1_[i], gctl::Dx);
		g1_[1] += -1e+3*gctl::gaussian_dist2d(x[0], x[1], gs1_[i], gctl::Dy);
	}

    cmg_.fill_model_gradient(0, f, g1_);
    fx += f;

    f = 0.0;
    for (int i = 0; i < gs2_.size(); i++)
	{
		f += -1e+3*gctl::gaussian_dist2d(x[0], x[1], gs2_[i]);
		g2_[0] += -1e+3*gctl::gaussian_dist2d(x[0], x[1], gs2_[i], gctl::Dx);
		g2_[1] += -1e+3*gctl::gaussian_dist2d(x[0], x[1], gs2_[i], gctl::Dy);
	}

    cmg_.fill_model_gradient(1, f, g2_);
    fx += f;
	fx += 0.666027; // 目标值

    g = cmg_.get_common_gradient();
    return fx;
}

void cfg::read_gaussians(std::string file, std::vector<gaussian_para2d> &gs)
{
    gaussian_para2d tmp_p;
	std::string tmp_str;
	std::stringstream tmp_ss;

    std::ifstream infile;
	gctl::open_infile(infile, file, ".txt");

	while(getline(infile, tmp_str))
	{
		if (tmp_str[0] == '#') continue;
		else
		{
			gctl::str2ss(tmp_str, tmp_ss);
			tmp_ss >> tmp_p.mu_x >> tmp_p.mu_y >> tmp_p.sigma_x >> tmp_p.sigma_y >> tmp_p.rho;
			gs.push_back(tmp_p);
		}
	}
	infile.close();
    return;
}

void cfg::routine(int argc, char *argv[])
{
    cmg_.init(2, 2);
    g1_.resize(2);
    g2_.resize(2);

    read_gaussians("example/data/gauss_model", gs1_);
    read_gaussians("example/data/gauss_model2", gs2_);

    x_.resize(2, 20.0);
    xm_.resize(2, 0.0);
    xs_.resize(2, 0.0);
    array<double> low(2, 0.0);
    array<double> high(2, 100.0);

    lgd_para para = default_lgd_para();
    para.alpha = 0.01;
    para.flight_times = 1000;
    set_lgd_para(para);

    set_lgd_record_trace();
    LGD_Minimize(x_, xm_, xs_, low, high);

    save_lgd_trace("trace.txt");
    std::cout << "x = (" << x_[0] << ", " << x_[1] << ")\n";
    return;
}

int main(int argc, char *argv[]) try
{
    cfg c;
    c.routine(argc, argv);
    return 0;
}
catch (const std::exception& e)
{
    GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
}