tmp

example/CMakeLists.txt

@@ -10,13 +10,14 @@ macro(add_example name switch)
     endif()
 endmacro()
 
-add_example(ex1 ON)
-add_example(ex2 ON)
-add_example(ex3 ON)
-add_example(ex4 ON)
-add_example(ex5 ON)
-add_example(ex6 ON)
-add_example(ex7 ON)
-add_example(ex8 ON)
-add_example(ex9 ON)
-add_example(ex10 ON)
+add_example(ex1 OFF)
+add_example(ex2 OFF)
+add_example(ex3 OFF)
+add_example(ex4 OFF)
+add_example(ex5 OFF)
+add_example(ex6 OFF)
+add_example(ex7 OFF)
+add_example(ex8 OFF)
+add_example(ex9 OFF)
+add_example(ex10 OFF)
+add_example(cfg_ex ON)

example/cfg_ex.cpp

@@ -0,0 +1,146 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * 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);
+}

example/data/forward.sh

@@ -0,0 +1,7 @@
+#!/bin/bash
+
+# 系统里有一个隐藏的错误设置 不然不需要在这里手动设置动态库的地址
+export DYLD_LIBRARY_PATH=/usr/local/lib:$DYLD_LIBRARY_PATH
+
+#g++ gaussian2d.cpp -o gaussian2d -lgctl -lnetcdfcxx_legacy -std=c++17 -I/opt/homebrew/include
+./gaussian2d gauss_model.txt

example/data/gauss_model.txt

@@ -0,0 +1,2 @@
+32.0 66.0 24.0 22.0 -0.1
+32.0 10.0 26.0 24.0 0.6

example/data/gauss_model2.txt

@@ -0,0 +1,2 @@
+75.0 25.0 22.0 28.0 0.5
+32.0 66.0 24.0 22.0 -0.1

example/data/gaussian2d.cpp

@@ -0,0 +1,70 @@
+/**
+ * @defgroup   MAIN main
+ *
+ * @brief      测试非线形最优化中的局部极小值问题
+ *
+ * @author     Zhangyi
+ * @date       2020
+ */
+
+#include "gctl/core.h"
+#include "gctl/io.h"
+
+// 声明我们所使用的二维高斯分布的参数
+static std::vector<gctl::gaussian_para2d> p;
+
+int main(int argc, char *argv[])
+{
+	std::ifstream infile;
+	gctl::open_infile(infile, argv[1], "");
+
+	gctl::gaussian_para2d tmp_p;
+	std::string tmp_str;
+	std::stringstream tmp_ss;
+	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;
+			p.push_back(tmp_p);
+		}
+	}
+
+	infile.close();
+
+	// 我们先来确认高斯分布的参数 生成解空间
+	int m = 101, n = 101;
+	gctl::array<double> gauss_dist(m*n, 0.0);
+	for (int t = 0; t < p.size(); t++)
+	{
+		for (int i = 0; i < n; i++)
+		{
+			for (int j = 0; j < m; j++)
+			{
+				gauss_dist[j + i*m] += -1e+3*gctl::gaussian_dist2d(1.0*j, 1.0*i, p[t]);
+			}
+		}
+	}
+
+	// 输出最小值及对应的解
+	double f_min = 1e+30, x_min, y_min;
+	for (int i = 0; i < n; i++)
+	{
+		for (int j = 0; j < m; j++)
+		{
+			if (gauss_dist[j + i*m] < f_min)
+			{
+				f_min = gauss_dist[j + i*m];
+				x_min = j; y_min = i;
+			}
+		}
+	}
+	std::cout << "f_min(" << x_min << ", " << y_min << ") = " << std::setprecision(10) << f_min << std::endl;
+
+	// 保存计算结果
+	gctl::save_netcdf_grid("gauss_dist", gauss_dist, m, n, 0, 1, 0, 1);
+	return 0;
+}

example/data/plot_movie.sh

@@ -0,0 +1,27 @@
+#!/usr/bin/env bash
+
+# 1. Create files needed in the loop
+cat << 'EOF' > pre.sh
+gmt begin pre
+	gmt set FONT_ANNOT_PRIMARY=15p,Times-Roman,black
+	gmt set FONT_LABEL=15p,Times-Roman,black
+	gmt set MAP_GRID_CROSS_SIZE_PRIMARY=5p
+	gmt set MAP_FRAME_PEN=thinnest,black
+	gmt set MAP_TICK_LENGTH_PRIMARY=4p/2p
+	gmt grd2cpt gauss_dist.nc -Clapaz -R0/100/0/100 -Z -D
+	gmt grdimage gauss_dist.nc -R0/100/0/100 -Bxag+l"x (m)" -Byag+l"y (m)" -JX15c/15c -X4.5c -Y1.5c
+gmt end
+EOF
+
+# 2. Set up the main frame script
+cat << 'EOF' > main.sh
+gmt begin
+	# Plot smooth blue curve and dark red dots at all steps so far
+	gmt convert trace.txt -qi0:${MOVIE_FRAME} > data.txt
+	gmt plot data.txt -W0.05p,white -R0/100/0/100 -JX15c/15c -X4.5c -Y1.5c
+	gmt plot data.txt -Sc0.05i -Gred
+gmt end
+EOF
+
+# 3. Run the movie
+gmt movie main.sh -Sbpre.sh -Cxga -Ttrace.txt -Vi -D24 -Zs -Nlgd_anim -Fmp4

example/ex1.cpp

@@ -1,31 +1,32 @@
 /********************************************************
- *  ██████╗ ███████╗████████╗██╗     
+ *  ██████╗  ██████╗████████╗██╗
  * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗███████╗   ██║   ██║     
- * ██║   ██║╚════██║   ██║   ██║     
- * ╚██████╔╝███████║   ██║   ███████╗
- *  ╚═════╝ ╚══════╝   ╚═╝   ╚══════╝                       
- * Generic Scientific Template Library
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
  *
  * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
  *
- * The GSTL is distributed under a dual licensing scheme. You can redistribute 
+ * 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 (LGPL) along with 
- * this program. If not, see <http://www.gnu.org/licenses/>.
+ * 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 GSTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GSTL's original author. As a rule, 
+ * 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"
+#include "../lib/optimization/lcg.h"
+#include "gctl/graphic/gnuplot.h"
 
 #define M 1000
 #define N 800
@@ -125,12 +126,21 @@ int main(int argc, char const *argv[])
 	test.LCG_Minimize(m, B, gctl::LCG_CG, ofile);
 	ofile << "maximal difference: " << max_diff(fm, m) << std::endl;
 
-	m.assign_all(0.0);
+	test.save_convergence("convergence");
+	gctl::gnuplot gt;
+	gt.to_buffer();
+	gt.send("set terminal png size 800,600");
+	gt.send("set output \"convergence.png\"");
+	gt.send("plot \"convergence.txt\" using 1:2 with lines");
+	gt.send("set output");
+	gt.send_buffer();
+
+	m.assign(0.0);
 	
 	test.LCG_Minimize(m, B, gctl::LCG_PCG, ofile);
 	ofile << "maximal difference: " << max_diff(fm, m) << std::endl;
 
-	m.assign_all(0.0);
+	m.assign(0.0);
 
 	test.LCG_Minimize(m, B, gctl::LCG_CGS, ofile);
 	ofile << "maximal difference: " << max_diff(fm, m) << std::endl;
@@ -138,12 +148,12 @@ int main(int argc, char const *argv[])
 
 	test.set_lcg_message(gctl::LCG_SOLUTION);
 
-	m.assign_all(0.0);
+	m.assign(0.0);
 
 	test.LCG_Minimize(m, B, gctl::LCG_BICGSTAB);
 	std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
 
-	m.assign_all(0.0);
+	m.assign(0.0);
 
 	test.LCG_Minimize(m, B, gctl::LCG_BICGSTAB2);
 	std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
@@ -151,12 +161,12 @@ int main(int argc, char const *argv[])
 	gctl::array<double> low(N, 1.0);
 	gctl::array<double> hig(N, 2.0);
 
-	m.assign_all(0.0);
+	m.assign(0.0);
 
 	test.LCG_MinimizeConstrained(m, B, low, hig, gctl::LCG_PG);
 	std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
 
-	m.assign_all(0.0);
+	m.assign(0.0);
 
 	test.LCG_MinimizeConstrained(m, B, low, hig, gctl::LCG_SPG);
 	std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;

example/ex10.cpp

@@ -25,11 +25,11 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#include "gctl/core.h"
-#include "gctl/algorithm.h"
+#include <iostream>
+#include <iomanip>
 #include "../lib/optimization.h"
-#include "iostream"
-#include "iomanip"
+#include "gctl/core.h"
+#include "gctl/algorithms.h"
 
 using std::cout;
 using std::endl;

example/ex2.cpp

@@ -1,24 +1,24 @@
 /********************************************************
- *  ██████╗ ███████╗████████╗██╗     
+ *  ██████╗  ██████╗████████╗██╗
  * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗███████╗   ██║   ██║     
- * ██║   ██║╚════██║   ██║   ██║     
- * ╚██████╔╝███████║   ██║   ███████╗
- *  ╚═════╝ ╚══════╝   ╚═╝   ╚══════╝                       
- * Generic Scientific Template Library
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
  *
  * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
  *
- * The GSTL is distributed under a dual licensing scheme. You can redistribute 
+ * 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 (LGPL) along with 
- * this program. If not, see <http://www.gnu.org/licenses/>.
+ * 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 GSTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GSTL's original author. As a rule, 
+ * 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. 

example/ex4.cpp

@@ -25,11 +25,11 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#include "gctl/core.h"
-#include "gctl/algorithm.h"
+#include <iostream>
+#include <iomanip>
 #include "../lib/optimization.h"
-#include "iostream"
-#include "iomanip"
+#include "gctl/core.h"
+#include "gctl/algorithms.h"
 
 using std::cout;
 using std::endl;

example/ex5.cpp

@@ -25,9 +25,9 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#include "gctl/core.h"
-#include "gctl/algorithm.h"
 #include "../lib/optimization.h"
+#include "gctl/core.h"
+#include "gctl/algorithms.h"
 
 int main(int argc, char const *argv[])
 {

example/ex6.cpp

@@ -25,9 +25,9 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#include "gctl/core.h"
-#include "gctl/algorithm.h"
 #include "../lib/optimization.h"
+#include "gctl/core.h"
+#include "gctl/algorithms.h"
 
 // get random floating points
 double random_double(double l, double t)

example/ex7.cpp

@@ -25,9 +25,9 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#include "gctl/core.h"
-#include "gctl/algorithm.h"
 #include "../lib/optimization.h"
+#include "gctl/core.h"
+#include "gctl/algorithms.h"
 
 typedef gctl::array<std::complex<double>> cd_array;
 
@@ -121,15 +121,15 @@ int main(int argc, char const *argv[])
     test.CLCG_Minimize(m, B, gctl::CLCG_BICG_SYM);
     std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
 
-    m.assign_all(std::complex<double>(0.0, 0.0));
+    m.assign(std::complex<double>(0.0, 0.0));
     test.CLCG_Minimize(m, B, gctl::CLCG_BICG);
     std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
 
-    m.assign_all(std::complex<double>(0.0, 0.0));
+    m.assign(std::complex<double>(0.0, 0.0));
     test.CLCG_Minimize(m, B, gctl::CLCG_CGS);
     std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
 
-    m.assign_all(std::complex<double>(0.0, 0.0));
+    m.assign(std::complex<double>(0.0, 0.0));
     test.CLCG_Minimize(m, B, gctl::CLCG_BICGSTAB);
     std::clog << "maximal difference: " << max_diff(fm, m) << std::endl;
 	return 0;

example/ex8.cpp

@@ -25,9 +25,9 @@
  * Also add information on how to contact you by electronic and paper mail.
  ******************************************************/
 
-#include "gctl/core.h"
-#include "gctl/algorithm.h"
 #include "../lib/optimization.h"
+#include "gctl/core.h"
+#include "gctl/algorithms.h"
 
 #define M 1000
 #define N 900

example/ex9.cpp

@@ -1,24 +1,24 @@
 /********************************************************
- *  ██████╗ ███████╗████████╗██╗     
+ *  ██████╗  ██████╗████████╗██╗
  * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗███████╗   ██║   ██║     
- * ██║   ██║╚════██║   ██║   ██║     
- * ╚██████╔╝███████║   ██║   ███████╗
- *  ╚═════╝ ╚══════╝   ╚═╝   ╚══════╝                       
- * Generic Scientific Template Library
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
  *
  * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
  *
- * The GSTL is distributed under a dual licensing scheme. You can redistribute 
+ * 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 (LGPL) along with 
- * this program. If not, see <http://www.gnu.org/licenses/>.
+ * 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 GSTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GSTL's original author. As a rule, 
+ * 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. 

lib/optimization.h

@@ -28,16 +28,17 @@
 #ifndef _GCTL_OPTIMIZATION_H
 #define _GCTL_OPTIMIZATION_H
 
+#include "optimization/lcg.h"
 #include "optimization/loss_func.h"
 #include "optimization/lu.h"
 #include "optimization/cholesky.h"
 #include "optimization/svd.h"
-#include "optimization/lcg.h"
 #include "optimization/clcg.h"
 #include "optimization/lgd.h"
 #include "optimization/lbfgs.h"
 #include "optimization/sgd.h"
 #include "optimization/gradnorm.h"
 #include "optimization/dwa.h"
+#include "optimization/cmn_grad.h"
 
 #endif // _GCTL_OPTIMIZATION_H

lib/optimization/cholesky.h

@@ -50,6 +50,6 @@ namespace gctl
         void operator=(const gctl::cholesky&) = delete;
         matrix<double> &decomposedMatrix;
     };
-}
+};
 
 #endif // _CHOLESKY_H

lib/optimization/clcg.h

@@ -28,9 +28,10 @@
 #ifndef _GCTL_CLCG_H
 #define _GCTL_CLCG_H
 
+#include "gctl/utility.h"
 #include "gctl/core.h"
 #include "gctl/maths.h"
-#include "gctl/algorithm.h"
+#include "gctl/algorithms.h"
 
 #include "gctl_optimization_config.h"
 #ifdef GCTL_OPTIMIZATION_TOML
@@ -184,6 +185,6 @@ namespace gctl
         void clbicgstab(array<std::complex<double> > &m, const array<std::complex<double> > &B, std::ostream &ss = std::clog);
         void cltfqmr(array<std::complex<double> > &m, const array<std::complex<double> > &B, std::ostream &ss = std::clog);
     };
-}
+};
 
 #endif // _GCTL_CLCG_H

lib/optimization/cmn_grad.cpp

@@ -0,0 +1,167 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * 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)
+{
+    // Ax product of the linear system
+    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)
+{
+    zero_iter_ = true;
+    Ln_ = Ln; 
+    Mn_ = Mn;
+    g_.resize(Mn_);
+    B_.resize(Mn_);
+    G_.resize(Ln_, Mn_);
+    t_.resize(Ln_);
+    gm_.resize(Ln_);
+    fx_.resize(Ln_);
+    fx0_.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, double fx, 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);
+    gm_[id] = g.module();
+    fx_[id] = fx;
+    if (zero_iter_) fx0_[id] = fx;
+    filled_[id] = true;
+    return;
+}
+
+const gctl::_1d_array &gctl::common_gradient::get_common_gradient(bool normalized, bool fixed_w)
+{
+    zero_iter_ = false; // set to false after the first evaluation.
+
+    for (size_t i = 0; i < Ln_; i++)
+    {
+        if (!filled_[i]) throw std::runtime_error("[gctl::common_gradient] Unfilled model gradient.");
+    }
+    filled_.assign(false);
+
+    if (!fixed_w)
+    {
+        // 计算权重
+        double a;
+        for (size_t i = 0; i < Ln_; i++)
+        {
+            // the bigger gradient module is, the bigger weight is
+            // the faster convergence is, the smaller weight is
+            a = gm_[i]*fx_[i]/fx0_[i];
+            w_[i] = 1.0/a;
+        }
+    }
+
+    rcd_wgts_.push_back(w_);
+    rcd_fxs_.push_back(fx_);
+
+    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();
+    return g_;
+}
+
+void gctl::common_gradient::save_records(std::string file)
+{
+    std::ofstream fout;
+    open_outfile(fout, file, ".csv");
+
+    fout << "num";
+    for (size_t j = 0; j < Ln_; j++)
+    {
+        fout << ",l" << j;
+    }
+
+    for (size_t j = 0; j < Ln_; j++)
+    {
+        fout << ",w" << j;
+    }
+
+    for (size_t i = 0; i < rcd_wgts_.size(); i++)
+    {
+        fout << std::endl << i;
+        for (size_t j = 0; j < Ln_; j++)
+        {
+            fout << "," << rcd_fxs_[i][j];
+        }
+
+        for (size_t j = 0; j < Ln_; j++)
+        {
+            fout << "," << rcd_wgts_[i][j];
+        }
+    }
+    fout.close();
+    return;
+}

lib/optimization/cmn_grad.h

@@ -0,0 +1,114 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * 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.
+ ******************************************************/
+
+#ifndef _GCTL_COMMON_GRADIENT_H
+#define _GCTL_COMMON_GRADIENT_H
+
+#include "lcg.h"
+
+namespace gctl
+{
+    class common_gradient : public lcg_solver
+    {
+    public:
+        common_gradient(); ///< 构造函数
+
+        /**
+         * @brief Construct a new common_gradient object
+         * 
+         * @param Ln Number of loss functions
+         * @param Mn Number of model parameters
+         */
+        common_gradient(size_t Ln, size_t Mn);
+
+        virtual ~common_gradient(); ///< 析构函数
+
+        virtual void LCG_Ax(const array<double> &x, array<double> &ax); ///< 计算Ax
+
+        /**
+         * @brief Configure the solver's setups
+         * 
+         * @param para LCG solver parameters
+         */
+        void set_solver(const lcg_para &para);
+
+        /**
+         * @brief Set the weights for the loss functions.
+         * 
+         * The number of weights equal to the number of the loss functions.
+         * The bigger weights is the calculated gradient is more dependent
+         * on the corresponding gradients.
+         */
+        void set_weights(const _1d_array &w);
+
+        /**
+         * @brief Initialize the common_gradient object
+         * 
+         * @param Ln Number of loss functions
+         * @param Mn Number of model parameters
+         */
+        void init(size_t Ln, size_t Mn);
+
+        /**
+         * @brief Fill the model gradient
+         *
+         * @param id Loss function index
+         * @param fx Objective value
+         * @param g Model gradient
+         */
+        void fill_model_gradient(size_t id, double fx, const _1d_array &g);
+
+        /**
+         * @brief Get the conflict free gradient
+         * 
+         * @param normalized Normalize the output gradient
+         * @param fixed_w Fixed weights
+         * @return Calculated model gradient
+         */
+        const _1d_array &get_common_gradient(bool normalized = true, bool fixed_w = true);
+
+        /**
+         * @brief Save the recorded weights.
+         * 
+         * @param file Output file name
+         */
+        void save_records(std::string file);
+
+    private:
+        bool zero_iter_;
+        size_t Ln_, Mn_; // Ln_: loss_func number，Mn_: model number
+        _2d_matrix G_; // kernel martix
+        _1d_array B_, g_, t_, x_; // variables of the linear system
+        _1d_array gm_, w_; // gradient module and functions' weight
+        _1d_array fx_, fx0_; // functions' value, initial functions' value
+        array<bool> filled_; // new gradient filled for the current round of evaluation
+        std::vector<array<double> > rcd_wgts_; // weights records
+        std::vector<array<double> > rcd_fxs_; // fx records
+    };
+};
+
+#endif // _GCTL_COMMON_GRADIENT_H

lib/optimization/dwa.cpp

@@ -99,7 +99,7 @@ double gctl::dwa::DWALoss(array<double> &g)
 
     fx_c_ = 0;
     multi_fx_ = 0.0;
-    grad_.assign_all(0.0);
+    grad_.assign(0.0);
     return fx;
 }
 
@@ -115,15 +115,27 @@ void gctl::dwa::set_normal_sum(double k)
     return;
 }
 
-void gctl::dwa::get_records(array<double> &logs)
+void gctl::dwa::save_records(std::string file)
 {
-    logs.resize(fx_n_*rcd_wgts_.size());
+    std::ofstream fout;
+    open_outfile(fout, file, ".csv");
+
+    fout << "Num";
+    for (size_t j = 0; j < fx_n_; j++)
+    {
+        fout << ",l" << j;
+    }
+    fout << std::endl;
+
     for (size_t i = 0; i < rcd_wgts_.size(); i++)
     {
+        fout << i;
         for (size_t j = 0; j < fx_n_; j++)
         {
-            logs[i*fx_n_ + j] = rcd_wgts_[i][j];
+            fout << "," << rcd_wgts_[i][j];
         }
+        fout << std::endl;
     }
+    fout.close();
     return;
 }

lib/optimization/dwa.h

@@ -29,6 +29,7 @@
 #define _GCTL_DWA_H
 
 #include "gctl/core.h"
+#include "gctl/utility.h"
 
 namespace gctl
 {
@@ -103,12 +104,12 @@ namespace gctl
         void set_normal_sum(double k);
 
         /**
-         * @brief Get the recorded weights. Size of the log equals the function size times iteration times.
+         * @brief Save the recorded weights.
          * 
-         * @param logs Output log
+         * @param file Output file name
          */
-        void get_records(array<double> &logs);
+        void save_records(std::string file);
     };
-}
+};
 
 #endif // _GCTL_DWA_H

lib/optimization/gradnorm.cpp

@@ -213,7 +213,7 @@ double gctl::grad_norm::GradNormLoss(array<double> &g)
 
     fx_c_ = 0;
     multi_fx_ = 0.0;
-    grad_.assign_all(0.0);
+    grad_.assign(0.0);
     return fx;
 }
 

lib/optimization/gradnorm.h

@@ -154,6 +154,6 @@ namespace gctl
          */
         void save_records(std::string file);
     };
-}
+};
 
 #endif // _GCTL_GRADNORM_H

lib/optimization/lbfgs.h

@@ -28,9 +28,10 @@
 #ifndef _GCTL_LBFGS_H
 #define _GCTL_LBFGS_H
 
+#include "gctl/utility.h"
 #include "gctl/core.h"
 #include "gctl/maths.h"
-#include "gctl/algorithm.h"
+#include "gctl/algorithms.h"
 
 #include "gctl_optimization_config.h"
 #ifdef GCTL_OPTIMIZATION_TOML
@@ -554,6 +555,6 @@ namespace gctl
          */
         double LBFGS_MinimizePreconditioned(array<double> &m, std::ostream &ss = std::clog, bool err_throw = false);
     };
-}
+};
 
 #endif // _GCTL_LBFGS_H

lib/optimization/lcg.cpp

@@ -93,6 +93,19 @@ void gctl::lcg_solver::set_lcg_para(const lcg_para &in_param)
     return;
 }
 
+void gctl::lcg_solver::save_convergence(std::string filename)
+{
+    std::ofstream ofs;
+    open_outfile(ofs, filename, ".txt");
+
+    for (size_t i = 0; i < rcd.size(); ++i)
+    {
+        ofs << i << " " << rcd[i] << std::endl;
+    }
+    ofs.close();
+    return;
+}
+
 gctl::lcg_para gctl::lcg_solver::default_lcg_para()
 {
     lcg_para dp = lcg_defparam;
@@ -345,6 +358,7 @@ void gctl::lcg_solver::lcg(array<double> &m, const array<double> &B, std::ostrea
     gk.resize(n_size);
     dk.resize(n_size);
     Adk.resize(n_size);
+    if (!rcd.empty()) rcd.clear();
 
     LCG_Ax(m, Adk);
 
@@ -361,6 +375,7 @@ void gctl::lcg_solver::lcg(array<double> &m, const array<double> &B, std::ostrea
     {
         if (lcg_param_.abs_diff) residual = sqrt(gk_mod)/n_size;
         else residual = gk_mod/g0_mod;
+        rcd.push_back(residual);
 
         if (LCG_Progress(m, residual, lcg_param_, t, ss))
         {
@@ -415,6 +430,7 @@ void gctl::lcg_solver::lpcg(array<double> &m, const array<double> &B, std::ostre
     // locate memory
     rk.resize(n_size); zk.resize(n_size);
     dk.resize(n_size); Adk.resize(n_size);
+    if (!rcd.empty()) rcd.clear();
 
     LCG_Ax(m, Adk);
 
@@ -436,6 +452,7 @@ void gctl::lcg_solver::lpcg(array<double> &m, const array<double> &B, std::ostre
     {
         if (lcg_param_.abs_diff) residual = sqrt(rk_mod)/n_size;
         else residual = rk_mod/r0_mod;
+        rcd.push_back(residual);
 
         if (LCG_Progress(m, residual, lcg_param_, t, ss))
         {
@@ -499,6 +516,7 @@ void gctl::lcg_solver::lcgs(array<double> &m, const array<double> &B, std::ostre
     uk.resize(n_size);
     qk.resize(n_size);
     wk.resize(n_size);
+    if (!rcd.empty()) rcd.clear();
 
     LCG_Ax(m, Apx);
 
@@ -520,6 +538,7 @@ void gctl::lcg_solver::lcgs(array<double> &m, const array<double> &B, std::ostre
     {
         if (lcg_param_.abs_diff) residual = sqrt(rk_mod)/n_size;
         else residual = rk_mod/r0_mod;
+        rcd.push_back(residual);
 
         if (LCG_Progress(m, residual, lcg_param_, t, ss))
         {
@@ -584,6 +603,7 @@ void gctl::lcg_solver::lbicgstab(array<double> &m, const array<double> &B, std::
 	pk.resize(n_size); Adk.resize(n_size);
 	sk.resize(n_size); Apx.resize(n_size);
     yk.resize(n_size);
+    if (!rcd.empty()) rcd.clear();
 
 	LCG_Ax(m, Adk);
 
@@ -602,6 +622,7 @@ void gctl::lcg_solver::lbicgstab(array<double> &m, const array<double> &B, std::
 	{
 		if (lcg_param_.abs_diff) residual = sqrt(rk_mod)/n_size;
 		else residual = rk_mod/r0_mod;
+        rcd.push_back(residual);
 
 		if (LCG_Progress(m, residual, lcg_param_, t, ss))
         {
@@ -669,6 +690,7 @@ void gctl::lcg_solver::lbicgstab2(array<double> &m, const array<double> &B, std:
 	pk.resize(n_size); Adk.resize(n_size);
 	sk.resize(n_size); Apx.resize(n_size);
     yk.resize(n_size);
+    if (!rcd.empty()) rcd.clear();
 
 	LCG_Ax(m, Adk);
 
@@ -688,6 +710,7 @@ void gctl::lcg_solver::lbicgstab2(array<double> &m, const array<double> &B, std:
 	{
 		if (lcg_param_.abs_diff) residual = sqrt(rk_mod)/n_size;
 		else residual = rk_mod/r0_mod;
+        rcd.push_back(residual);
 
 		if (LCG_Progress(m, residual, lcg_param_, t, ss))
         {
@@ -804,6 +827,8 @@ void gctl::lcg_solver::lpg(array<double> &m, const array<double> &B,
     gk_new.resize(n_size);
     sk.resize(n_size);
     yk.resize(n_size);
+    if (!rcd.empty()) rcd.clear();
+
     double alpha_k = lcg_param_.step;
 
     vectop(m, hig);
@@ -823,6 +848,7 @@ void gctl::lcg_solver::lpg(array<double> &m, const array<double> &B,
     {
         if (lcg_param_.abs_diff) residual = sqrt(gk_mod)/n_size;
         else residual = gk_mod/g0_mod;
+        rcd.push_back(residual);
 
         if (LCG_Progress(m, residual, lcg_param_, t, ss))
         {
@@ -889,6 +915,8 @@ void gctl::lcg_solver::lspg(array<double> &m, const array<double> &B,
     yk.resize(n_size);
     dk.resize(n_size);
     qk_m.resize(lcg_param_.maxi_m);
+    if (!rcd.empty()) rcd.clear();
+
     double lambda_k = lcg_param_.step;
     double qk = 0;
 
@@ -923,6 +951,7 @@ void gctl::lcg_solver::lspg(array<double> &m, const array<double> &B,
     {
         if (lcg_param_.abs_diff) residual = sqrt(gk_mod)/n_size;
         else residual = gk_mod/g0_mod;
+        rcd.push_back(residual);
 
         if (LCG_Progress(m, residual, lcg_param_, t, ss))
         {

lib/optimization/lcg.h

@@ -28,9 +28,7 @@
 #ifndef _GCTL_LCG_H
 #define _GCTL_LCG_H
 
-#include "gctl/core.h"
-#include "gctl/maths.h"
-#include "gctl/algorithm.h"
+#include <iostream>
 
 #include "gctl_optimization_config.h"
 #ifdef GCTL_OPTIMIZATION_TOML
@@ -41,6 +39,10 @@
 #include "windows.h"
 #endif // _WINDOWS || __WIN32__
 
+#include "gctl/utility/stream_t.h"
+#include "gctl/utility/stream.h"
+#include "gctl/maths/linear_algebra.h"
+
 namespace gctl
 {
 	/**
@@ -199,6 +201,7 @@ namespace gctl
     	array<double> Apx, uk, qk, qk_m, wk;
 		array<double> m_new, gk_new;
     	array<double> sk, yk;
+		std::vector<double> rcd;
 
 		/**
 		 * @brief Display info of a given return code. This is a private function
@@ -265,6 +268,13 @@ namespace gctl
 		 */
 		void set_lcg_para(const lcg_para &param);
 		
+		/**
+		 * @brief Save the convergence history to a file.
+		 * 
+		 * @param filename Filename of the output file.
+		 */
+		void save_convergence(std::string filename);
+
 		/**
 		 * @brief Return a lcg_para object with default values
 		 * 
@@ -383,6 +393,6 @@ namespace gctl
 		 */
 		void lspg(array<double> &m, const array<double> &B, const array<double> &low, const array<double> &hig, std::ostream &ss = std::clog);
 	};
-}
+};
 
 #endif // _GCTL_LCG_H

lib/optimization/lgd.cpp

@@ -500,7 +500,7 @@ gctl::lgd_return_code gctl::lgd_solver::lgd(array<double> &best_m, array<double>
 
             if (lgd_save_trace_)
             {
-                lgd_trace_.append_array(best_m);
+                lgd_trace_.concat(best_m);
                 lgd_trace_times_++;
             }
         }
@@ -671,7 +671,7 @@ gctl::lgd_return_code gctl::lgd_solver::lgd_momentum(array<double> &best_m, arra
 
             if (lgd_save_trace_)
             {
-                lgd_trace_.append_array(best_m);
+                lgd_trace_.concat(best_m);
                 lgd_trace_times_++;
             }
         }

lib/optimization/lgd.h

@@ -28,11 +28,6 @@
 #ifndef _GCTL_LGD_H
 #define _GCTL_LGD_H
 
-#include "gctl/core.h"
-#include "gctl/io.h"
-#include "gctl/maths.h"
-#include "gctl/algorithm.h"
-
 #include "gctl_optimization_config.h"
 #ifdef GCTL_OPTIMIZATION_TOML
 #include "toml.hpp"
@@ -46,6 +41,11 @@
 #include "omp.h"
 #endif // GSTL_OPENMP
 
+#include "gctl/core.h"
+#include "gctl/io.h"
+#include "gctl/maths.h"
+#include "gctl/algorithms.h"
+
 namespace gctl
 {
 	/**
@@ -287,6 +287,6 @@ namespace gctl
         bool lgd_silent_, lgd_has_range_, lgd_has_alpha_, lgd_save_trace_;
         array<double> lgd_low_, lgd_hig_, lgd_alpha_, lgd_trace_;
     };
-}
+};
 
 #endif // _GCTL_LGD_H

lib/optimization/loss_func.h

@@ -120,6 +120,6 @@ namespace gctl
         array<double> tars_, diff_;
         array<double> us_;
     };
-}
+};
 
 #endif // _GCTL_LOSS_FUNC_H

lib/optimization/lu.h

@@ -51,6 +51,6 @@ namespace gctl
         matrix<double> &decomposedMatrix; // Output matrix after decomposition
         array<int> rowPermutation; // Permutation of rows during pivoting
     };
-}
+};
 
 #endif // _GCTL_LU_H

lib/optimization/sgd.h

@@ -28,8 +28,9 @@
 #ifndef _GCTL_SGD_H
 #define _GCTL_SGD_H
 
+#include "gctl/utility.h"
 #include "gctl/core.h"
-#include "gctl/algorithm.h"
+#include "gctl/algorithms.h"
 
 #include "gctl_optimization_config.h"
 #ifdef GCTL_OPTIMIZATION_TOML
@@ -196,6 +197,6 @@ namespace gctl
 
         void SGD_Minimize(array<double> &m, sgd_solver_type solver_id = ADAM, std::ostream &ss = std::clog, bool verbose = true, bool err_throw = false);
     };
-}
+};
 
 #endif // _GCTL_SGD_H

lib/optimization/svd.cpp

@@ -129,8 +129,8 @@ void gctl::svd::decompose(const matrix<double> &src_mat)
 
         for(int iter=0; diff >= epsilon && iter < maxi_iteration; iter++)
         {
-            next_left_vector.assign_all(0.0);
-            next_right_vector.assign_all(0.0);
+            next_left_vector.assign(0.0);
+            next_right_vector.assign(0.0);
 
             for(int i=0;i<M;i++)
                 for(int j=0;j<N;j++)

lib/optimization/svd.h

@@ -29,7 +29,7 @@
 #define _GCTL_SVD_H
 
 #include "gctl/core.h"
-#include "gctl/algorithm.h"
+#include "gctl/algorithms.h"
 
 namespace gctl
 {
@@ -74,6 +74,6 @@ namespace gctl
         int maxi_iteration, K;
         double epsilon;
     };
-}
+};
 
 #endif // _GCTL_SVD_H

tool/levy_planner/levy_planner.cpp

@@ -1,24 +1,24 @@
 /********************************************************
- *  ██████╗ ███████╗████████╗██╗     
+ *  ██████╗  ██████╗████████╗██╗
  * ██╔════╝ ██╔════╝╚══██╔══╝██║
- * ██║  ███╗███████╗   ██║   ██║     
- * ██║   ██║╚════██║   ██║   ██║     
- * ╚██████╔╝███████║   ██║   ███████╗
- *  ╚═════╝ ╚══════╝   ╚═╝   ╚══════╝                       
- * Generic Scientific Template Library
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
  *
  * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
  *
- * The GSTL is distributed under a dual licensing scheme. You can redistribute 
+ * 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 (LGPL) along with 
- * this program. If not, see <http://www.gnu.org/licenses/>.
+ * 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 GSTL, please consider the option to obtain a commercial license for a 
- * fee. These licenses are offered by the GSTL's original author. As a rule, 
+ * 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.