From b4b4275576452397c158b5b48f93b38749731784 Mon Sep 17 00:00:00 2001 From: Yi Zhang Date: Tue, 10 Sep 2024 20:04:47 +0800 Subject: [PATCH] initial upload --- .gitignore | 37 +- CMakeLists.txt | 28 + GCTL_OPTIMIZATIONConfig.cmake.in | 25 + README.md | 17 +- config.h.in | 1 + example/CMakeLists.txt | 20 + example/ex1.cpp | 164 ++ example/ex2.cpp | 175 ++ example/ex3.cpp | 69 + example/ex4.cpp | 91 + example/ex5.cpp | 78 + example/ex6.cpp | 90 + example/ex7.cpp | 124 ++ example/ex8.cpp | 233 +++ installer | 54 + lib/CMakeLists.txt | 64 + lib/optimization.h | 43 + lib/optimization/cholesky.cpp | 127 ++ lib/optimization/cholesky.h | 55 + lib/optimization/clcg.cpp | 373 ++++ lib/optimization/clcg.h | 166 ++ lib/optimization/dwa.cpp | 129 ++ lib/optimization/dwa.h | 114 ++ lib/optimization/gctl_optimization_config.h | 1 + lib/optimization/gradnorm.cpp | 328 ++++ lib/optimization/gradnorm.h | 159 ++ lib/optimization/lbfgs.cpp | 1897 +++++++++++++++++++ lib/optimization/lbfgs.h | 559 ++++++ lib/optimization/lcg.cpp | 1002 ++++++++++ lib/optimization/lcg.h | 387 ++++ lib/optimization/lgd.cpp | 505 +++++ lib/optimization/lgd.h | 212 +++ lib/optimization/loss_func.cpp | 100 + lib/optimization/loss_func.h | 61 + lib/optimization/lu.cpp | 152 ++ lib/optimization/lu.h | 56 + lib/optimization/sgd.cpp | 634 +++++++ lib/optimization/sgd.h | 201 ++ lib/optimization/svd.cpp | 184 ++ lib/optimization/svd.h | 71 + 40 files changed, 8751 insertions(+), 35 deletions(-) create mode 100644 CMakeLists.txt create mode 100644 GCTL_OPTIMIZATIONConfig.cmake.in create mode 100644 config.h.in create mode 100644 example/CMakeLists.txt create mode 100644 example/ex1.cpp create mode 100644 example/ex2.cpp create mode 100644 example/ex3.cpp create mode 100644 example/ex4.cpp create mode 100644 example/ex5.cpp create mode 100644 example/ex6.cpp create mode 100644 example/ex7.cpp create mode 100644 example/ex8.cpp create mode 100755 installer create mode 100644 lib/CMakeLists.txt create mode 100644 lib/optimization.h create mode 100644 lib/optimization/cholesky.cpp create mode 100644 lib/optimization/cholesky.h create mode 100644 lib/optimization/clcg.cpp create mode 100644 lib/optimization/clcg.h create mode 100644 lib/optimization/dwa.cpp create mode 100644 lib/optimization/dwa.h create mode 100644 lib/optimization/gctl_optimization_config.h create mode 100644 lib/optimization/gradnorm.cpp create mode 100644 lib/optimization/gradnorm.h create mode 100644 lib/optimization/lbfgs.cpp create mode 100644 lib/optimization/lbfgs.h create mode 100644 lib/optimization/lcg.cpp create mode 100644 lib/optimization/lcg.h create mode 100644 lib/optimization/lgd.cpp create mode 100644 lib/optimization/lgd.h create mode 100644 lib/optimization/loss_func.cpp create mode 100644 lib/optimization/loss_func.h create mode 100644 lib/optimization/lu.cpp create mode 100644 lib/optimization/lu.h create mode 100644 lib/optimization/sgd.cpp create mode 100644 lib/optimization/sgd.h create mode 100644 lib/optimization/svd.cpp create mode 100644 lib/optimization/svd.h diff --git a/.gitignore b/.gitignore index e257658..f0bbcbb 100644 --- a/.gitignore +++ b/.gitignore @@ -1,34 +1,3 @@ -# ---> C++ -# Prerequisites -*.d - -# Compiled Object files -*.slo -*.lo -*.o -*.obj - -# Precompiled Headers -*.gch -*.pch - -# Compiled Dynamic libraries -*.so -*.dylib -*.dll - -# Fortran module files -*.mod -*.smod - -# Compiled Static libraries -*.lai -*.la -*.a -*.lib - -# Executables -*.exe -*.out -*.app - +.DS_Store +build/ +.vscode/ \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..b529f84 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,28 @@ +cmake_minimum_required(VERSION 3.15.2) +# 设置项目名称与语言 +project(GCTL_OPTIMIZATION VERSION 1.0) +# 添加配置配件编写的函数 +include(CMakePackageConfigHelpers) + +message(STATUS "Platform: " ${CMAKE_HOST_SYSTEM_NAME}) +message(STATUS "Install prefix: " ${CMAKE_INSTALL_PREFIX}) +message(STATUS "Processor: " ${CMAKE_HOST_SYSTEM_PROCESSOR}) + +option(GCTL_OPTIMIZATION_TOML "Use the TOML library" ON) + +message(STATUS "[GCTL_OPTIMIZATION] Use the TOML library: " ${GCTL_OPTIMIZATION_TOML}) + +find_package(GCTL REQUIRED) +include_directories(${GCTL_INC_DIR}) + +# 加入一个头文件配置，让cmake对源码进行操作 +configure_file( + "${PROJECT_SOURCE_DIR}/config.h.in" + "${PROJECT_SOURCE_DIR}/lib/optimization/gctl_optimization_config.h" + ) + +# 添加库源文件地址 +add_subdirectory(lib) + +# 去掉注释编译示例 +add_subdirectory(example) \ No newline at end of file diff --git a/GCTL_OPTIMIZATIONConfig.cmake.in b/GCTL_OPTIMIZATIONConfig.cmake.in new file mode 100644 index 0000000..5596957 --- /dev/null +++ b/GCTL_OPTIMIZATIONConfig.cmake.in @@ -0,0 +1,25 @@ +@PACKAGE_INIT@ + +set(@PROJECT_NAME@_VERSION "@PROJECT_VERSION@") +set_and_check(@PROJECT_NAME@_INSTALL_PREFIX "${PACKAGE_PREFIX_DIR}") +set_and_check(@PROJECT_NAME@_INC_DIR "${PACKAGE_PREFIX_DIR}/include") +set_and_check(@PROJECT_NAME@_INCLUDE_DIR "${PACKAGE_PREFIX_DIR}/include") +set_and_check(@PROJECT_NAME@_LIB_DIR "${PACKAGE_PREFIX_DIR}/lib") +set_and_check(@PROJECT_NAME@_LIBRARY_DIR "${PACKAGE_PREFIX_DIR}/lib") + +set(@PROJECT_NAME@_LIB gctl_optimization) +set(@PROJECT_NAME@_LIBRARY gctl_optimization) + +set(@PROJECT_NAME@_TOML @GCTL_OPTIMIZATION_TOML@) + +message(STATUS "[GCTL_OPTIMIZATION] Use the TOML library: " @GCTL_OPTIMIZATION_TOML@) + +if(NOT GCTL_FOUND) + find_package(GCTL REQUIRED) + include_directories(${GCTL_INC_DIR}) +endif() + +# include target information +include("${CMAKE_CURRENT_LIST_DIR}/@PROJECT_NAME@Targets.cmake") + +check_required_components(@PROJECT_NAME@) \ No newline at end of file diff --git a/README.md b/README.md index 8e28d1a..eaecad9 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,17 @@ -# gctl_optimization +## lcg_solver 共轭梯度求解器 +### 求解器参数设置 +#### 1. 从toml文件读取参数 +用户可以从toml文件中读取并设置求解器参数。所有的参数都定义在名为lcg的顶级表格下，可设置的参数及类型如下所示： + +```toml +[lcg] +max_iterations= +epsilon= +abs_diff=0|1 +restart_epsilon= +step= +sigma= +beta= +maxi_m= +``` \ No newline at end of file diff --git a/config.h.in b/config.h.in new file mode 100644 index 0000000..f495026 --- /dev/null +++ b/config.h.in @@ -0,0 +1 @@ +#cmakedefine GCTL_OPTIMIZATION_TOML \ No newline at end of file diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt new file mode 100644 index 0000000..f7e0244 --- /dev/null +++ b/example/CMakeLists.txt @@ -0,0 +1,20 @@ +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3") +set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin/examples) + +macro(add_example name switch) + if(${switch}) + add_executable(${name} ${name}.cpp) + set_target_properties(${name} PROPERTIES CXX_STANDARD 17 CXX_STANDARD_REQUIRED ON) + target_link_libraries(${name} PRIVATE ${GCTL_LIB}) + target_link_libraries(${name} PRIVATE gctl_optimization) + 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) \ No newline at end of file diff --git a/example/ex1.cpp b/example/ex1.cpp new file mode 100644 index 0000000..1fce63c --- /dev/null +++ b/example/ex1.cpp @@ -0,0 +1,164 @@ +/******************************************************** + * ██████╗ ███████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗███████╗ ██║ ██║ + * ██║ ██║╚════██║ ██║ ██║ + * ╚██████╔╝███████║ ██║ ███████╗ + * ╚═════╝ ╚══════╝ ╚═╝ ╚══════╝ + * Generic Scientific Template Library + * + * Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn) + * + * The GSTL 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 . + * + * 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, + * 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" + +#define M 1000 +#define N 800 + +double max_diff(const gctl::array &a, const gctl::array &b) +{ + double max = -1.0; + for (size_t i = 0; i < a.size(); i++) + { + max = std::max(sqrt((a[i] - b[i])*(a[i] - b[i])), max); + } + return max; +} + +class ex1 : public gctl::lcg_solver +{ +public: + ex1(); + virtual ~ex1(); + + // 计算共轭梯度的B项 + void cal_partb(const gctl::array &x, gctl::array &B); + + //定义共轭梯度中Ax的算法 + virtual void LCG_Ax(const gctl::array &x, gctl::array &ax); + virtual void LCG_Mx(const gctl::array &x, gctl::array &mx); + +private: + gctl::matrix kernel; // 普通二维数组做核矩阵 + gctl::array tmp_arr; // 中间结果数组 + gctl::array p; // 预优矩阵 +}; + +ex1::ex1() +{ + kernel.resize(M, N); + kernel.random(-1.0, 1.0, gctl::RdUniform); + + tmp_arr.resize(M); + p.resize(N); + for (size_t i = 0; i < N; i++) + { + p[i] = 1.0; + } + + double diag; + for (size_t i = 0; i < N; i++) + { + diag = 0.0; + for (size_t j = 0; j < M; j++) + { + diag += kernel[j][i]*kernel[j][i]; + } + p[i] = 1.0/diag; + } +} + +ex1::~ex1(){} + +void ex1::cal_partb(const gctl::array &x, gctl::array &B) +{ + LCG_Ax(x, B); + return; +} + +void ex1::LCG_Ax(const gctl::array &x, gctl::array &ax) +{ + matvec(tmp_arr, kernel, x); + matvec(ax, kernel, tmp_arr, gctl::Trans); + return; +} + +void ex1::LCG_Mx(const gctl::array &x, gctl::array &mx) +{ + vecmul(mx, p, x); + return; +} + +int main(int argc, char const *argv[]) +{ + // 生成一组正演解 + gctl::array fm(N); + gctl::random(fm, 1.0, 2.0, gctl::RdUniform); + + ex1 test; + + // 计算共轭梯度B项 + gctl::array B(N); + test.cal_partb(fm, B); + + // 声明一组解 + gctl::array m(N, 0.0); + + test.set_lcg_message(gctl::LCG_SOLUTION); + + std::ofstream ofile("log.txt"); + test.LCG_Minimize(m, B, gctl::LCG_CG, ofile); + ofile << "maximal difference: " << max_diff(fm, m) << std::endl; + + m.assign_all(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); + + test.LCG_Minimize(m, B, gctl::LCG_CGS, ofile); + ofile << "maximal difference: " << max_diff(fm, m) << std::endl; + ofile.close(); + + test.set_lcg_message(gctl::LCG_SOLUTION); + + m.assign_all(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); + + test.LCG_Minimize(m, B, gctl::LCG_BICGSTAB2); + std::clog << "maximal difference: " << max_diff(fm, m) << std::endl; + + gctl::array low(N, 1.0); + gctl::array hig(N, 2.0); + + m.assign_all(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); + + test.LCG_MinimizeConstrained(m, B, low, hig, gctl::LCG_SPG); + std::clog << "maximal difference: " << max_diff(fm, m) << std::endl; + return 0; +} \ No newline at end of file diff --git a/example/ex2.cpp b/example/ex2.cpp new file mode 100644 index 0000000..6b9a637 --- /dev/null +++ b/example/ex2.cpp @@ -0,0 +1,175 @@ +/******************************************************** + * ██████╗ ███████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗███████╗ ██║ ██║ + * ██║ ██║╚════██║ ██║ ██║ + * ╚██████╔╝███████║ ██║ ███████╗ + * ╚═════╝ ╚══════╝ ╚═╝ ╚══════╝ + * Generic Scientific Template Library + * + * Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn) + * + * The GSTL 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 . + * + * 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, + * 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" + +#define M 100 +#define N 90 + +// get random floating points +double random_double(double l, double t) +{ + return (t-l)*rand()*1.0/RAND_MAX + l; +} + +// get random integral numbers +int random_int(int small, int big) +{ + return (rand() % (big - small)) + small; +} + +class ex2 : public gctl::lgd_solver +{ +protected: + gctl::matrix kernel; + gctl::array obs, tmp; + +protected: + double LGD_Evaluate(const gctl::array &x, gctl::array &g); + +public: + ex2(); + virtual ~ex2(){} + + void CalObs(const gctl::array &x); +}; + +ex2::ex2() +{ + kernel.resize(M, N); + tmp.resize(M); + obs.resize(M); + + srand(time(0)); + // 添加一些大数 + int tmp_id, tmp_size; + double tmp_val; + for (int i = 0; i < M; i++) + { + tmp_size = random_int(25, 35); + for (int j = 0; j < tmp_size; j++) + { + tmp_id = random_int(0, N); + tmp_val = random_double(-10, 10); + + kernel[i][tmp_id] = tmp_val; + } + } +} + +double ex2::LGD_Evaluate(const gctl::array &x, gctl::array &g) +{ + for (int i = 0; i < M; i++) + { + tmp[i] = 0.0; + for (int j = 0; j < N; j++) + { + tmp[i] += kernel[i][j] * x[j]; + } + tmp[i] -= obs[i]; + //tmp[i] /= 1e-1; + } + + for (int j = 0; j < N; j++) + { + g[j] = 0.0; + for (int i = 0; i < M; i++) + { + g[j] += kernel[i][j]*tmp[i]; + } + g[j] *= 2.0/M; + } + + double sum = 0.0; + for (int i = 0; i < M; i++) + { + sum += tmp[i]*tmp[i]; + } + return sum/M; +} + +void ex2::CalObs(const gctl::array &x) +{ + // 计算正演值 + for (int i = 0; i < M; i++) + { + obs[i] = 0.0; + for (int j = 0; j < N; j++) + { + obs[i] += kernel[i][j]*x[j]; + } + // 添加噪声 + obs[i] += random_double(-1e-3, 1e-3); + } +} + +int main(int argc, char const *argv[]) +{ + gctl::array m(N, 0.0), mean_m(N, 0.0), stddev_m(N, 0.0), low(N), hig(N); + + // 生成一组正演解 包含一些大值和一些小值 + gctl::array fm(N); + int N2 = (int) N/2; + for (int i = 0; i < N2; i++) + { + //fm[i] = random_double(5, 10); + fm[i] = 10.0; + } + + for (int i = N2; i < N; i++) + { + //fm[i] = random_double(1, 2); + fm[i] = 1.0; + } + + for (int i = 0; i < N2; i++) + { + low[i] = 9.0; // 对解的范围进行约束 + hig[i] = 11.0; + } + + for (int i = N2; i < N; i++) + { + low[i] = 0.0; + hig[i] = 2.0; + } + + ex2 e; + e.CalObs(fm); + + gctl::lgd_para my_para = e.default_lgd_para(); + my_para.flight_times = 20000; + my_para.batch = 100; + e.set_lgd_para(my_para); + e.LGD_Minimize(m, mean_m, stddev_m, low, hig); + + for (int i = 0; i < N; i++) + { + std::cout << fm[i] << " " << m[i] << " " << mean_m[i] << " " << stddev_m[i] << " " << fabs(mean_m[i] - fm[i]) << std::endl; + } + return 0; +} \ No newline at end of file diff --git a/example/ex3.cpp b/example/ex3.cpp new file mode 100644 index 0000000..474eb94 --- /dev/null +++ b/example/ex3.cpp @@ -0,0 +1,69 @@ +#include "../lib/optimization.h" + +class TEST_FUNC : public gctl::lbfgs_solver +{ +public: + TEST_FUNC(); + ~TEST_FUNC(); + + virtual double LBFGS_Evaluate(const gctl::array &x, gctl::array &g); + + void Routine(); + +private: + gctl::array m_x; +}; + +TEST_FUNC::TEST_FUNC() +{ + m_x.resize(3, 0.0); +} + +TEST_FUNC::~TEST_FUNC(){} + +// test functions +// 3 = 3*x1 + x2 + 2*x3*x3 +// 1 = -3*x1 + 5*x2*x2 + 2*x1*x3 +// -12 = 25*x1*x2 + 20*x3 +double TEST_FUNC::LBFGS_Evaluate(const gctl::array &x, gctl::array &g) +{ + double f0,f1,f2,temp; + f0 = 3*x[0] + x[1] + 2*x[2]*x[2] - 3.012; //这里添加一点噪声 + f1 = -3*x[0] + 5*x[1]*x[1] + 2*x[0]*x[2] - 1.04252; + f2 = 25*x[0]*x[1] + 20*x[2] + 12.12479; + temp = sqrt(f0*f0+f1*f1+f2*f2); + + g[0] = 0.5*(6*f0+2*f1*(2*x[2]-3)+50*f2*x[1])/temp; + g[1] = 0.5*(2*f0+20*f1*x[1]+50*f2*x[0])/temp; + g[2] = 0.5*(8*f0*x[2]+4*f1*x[0]+40*f2)/temp; + return temp; +} + +void TEST_FUNC::Routine() +{ + gctl::lbfgs_para self_para = default_lbfgs_para(); + self_para.m = 10; + self_para.past = 5; + self_para.residual = 1e-10; + //self_para.min_step = 1e-30; + //self_para.max_linesearch = 40; + //self_para.linesearch = gctl::LBFGS_LINESEARCH_BACKTRACKING_WOLFE; + + set_lbfgs_para(self_para); + + std::ofstream ofile("log.txt"); + show_lbfgs_para(ofile); + + double fx = LBFGS_Minimize(m_x, ofile); + ofile.close(); + + m_x.show(); + return; +} + +int main(int argc, char const *argv[]) +{ + TEST_FUNC test; + test.Routine(); + return 0; +} \ No newline at end of file diff --git a/example/ex4.cpp b/example/ex4.cpp new file mode 100644 index 0000000..58cdef9 --- /dev/null +++ b/example/ex4.cpp @@ -0,0 +1,91 @@ +/******************************************************** + * ██████╗ ██████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗██║ ██║ ██║ + * ██║ ██║██║ ██║ ██║ + * ╚██████╔╝╚██████╗ ██║ ███████╗ + * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ + * 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 . + * + * 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 "gctl/core.h" +#include "gctl/algorithm.h" +#include "../lib/optimization.h" +#include "iostream" +#include "iomanip" + +using std::cout; +using std::endl; +using std::setw; + +int main(int argc, char const *argv[]) +{ + gctl::matrix A(4, 3); + for (int i = 0; i < A.row_size(); i++) + { + for (int j = 0; j < A.col_size(); j++) + { + A[i][j] = 3*(i+1) + j - 2; + } + } + A[3][1] = 1; + + cout<<"A(" << A.row_size() << ", " << A.col_size() << ") = " < A(5, 5); + for (int i = 0; i < 5; i++) + { + for (int j = 0; j < 5; j++) + { + A[i][j] = 3*(i+1) + j - 2; + } + } + // 注意A要满秩 + A[1][2] = 3.4; + A[4][1] = 2.1; + A[3][4] = 9.7; + A[2][3] = 2.7; + + std::cout<<"A(5, 5) = " < m(5, 0.5), x(5, 0.0); + gctl::array B(5); + for (int i = 0; i < 5; i++) + { + B[i] = 0.0; + for (int j = 0; j < 5; j++) + { + B[i] += A[i][j] * m[j]; + } + } + + gctl::lu glu(A); + glu.decompose(); + glu.solve(B, x); + + for (size_t i = 0; i < 5; i++) + { + std::cout << m[i] << " " << x[i] << std::endl; + } + + return 0; +} \ No newline at end of file diff --git a/example/ex6.cpp b/example/ex6.cpp new file mode 100644 index 0000000..6d92b1e --- /dev/null +++ b/example/ex6.cpp @@ -0,0 +1,90 @@ +/******************************************************** + * ██████╗ ██████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗██║ ██║ ██║ + * ██║ ██║██║ ██║ ██║ + * ╚██████╔╝╚██████╗ ██║ ███████╗ + * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ + * 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 . + * + * 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 "gctl/core.h" +#include "gctl/algorithm.h" +#include "../lib/optimization.h" + +// get random floating points +double random_double(double l, double t) +{ + return (t-l)*rand()*1.0/RAND_MAX + l; +} + +int main(int argc, char const *argv[]) +{ + srand(time(0)); + + gctl::matrix A(5, 5); + for (int i = 0; i < 5; i++) + { + for (int j = i; j < 5; j++) + { + if (i == j) A[i][j] = random_double(1.0, 3.0); + else A[i][j] = random_double(0.1, 1.0); + } + } + + for (int i = 0; i < 5; i++) + { + for (int j = i; j < 5; j++) + { + A[j][i] = A[i][j]; + } + } + + std::cout<<"A(5, 5) = " < m(5, 0.5), x(5, 0.0); + gctl::array B(5); + for (int i = 0; i < 5; i++) + { + B[i] = 0.0; + for (int j = 0; j < 5; j++) + { + B[i] += A[i][j] * m[j]; + } + } + + gctl::cholesky gck(A); + gck.decompose(); + gck.solve(B, x); + + for (size_t i = 0; i < 5; i++) + { + std::cout << m[i] << " " << x[i] << std::endl; + } + + return 0; +} \ No newline at end of file diff --git a/example/ex7.cpp b/example/ex7.cpp new file mode 100644 index 0000000..c91d867 --- /dev/null +++ b/example/ex7.cpp @@ -0,0 +1,124 @@ +/******************************************************** + * ██████╗ ██████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗██║ ██║ ██║ + * ██║ ██║██║ ██║ ██║ + * ╚██████╔╝╚██████╗ ██║ ███████╗ + * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ + * 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 . + * + * 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 "gctl/core.h" +#include "gctl/algorithm.h" +#include "../lib/optimization.h" + +typedef gctl::array> cd_array; + +#define N 1000 + +double max_diff(const cd_array &a, const cd_array &b) +{ + double max = -1; + std::complex t; + for (size_t i = 0; i < a.size(); i++) + { + t = a[i] - b[i]; + max = std::max(std::norm(t), max); + } + return max; +} + +class ex7 : public gctl::clcg_solver +{ +public: + ex7(); + virtual ~ex7(); + virtual void CLCG_Ax(const cd_array &x, cd_array &ax, gctl::matrix_layout_e layout, gctl::conjugate_type_e conj); + + // 计算共轭梯度的B项 + void cal_partb(const cd_array &x, cd_array &B); + +private: + gctl::matrix> kernel; // 普通二维数组做核矩阵 +}; + +ex7::ex7() +{ + gctl::array tmp(round(0.5*(N+1)*N)); + gctl::random(tmp, 1.0, 2.0, gctl::RdUniform); + + size_t c = 0; + kernel.resize(N, N); + for (int i = 0; i < N; i++) + { + for (int j = i; j < N; j++) + { + kernel[i][j] = tmp[c]; + kernel[j][i] = kernel[i][j]; + c++; + } + } +} + +ex7::~ex7(){} + +void ex7::cal_partb(const cd_array &x, cd_array &B) +{ + gctl::matvec(B, kernel, x); + return; +} + +void ex7::CLCG_Ax(const cd_array &x, cd_array &ax, gctl::matrix_layout_e layout, gctl::conjugate_type_e conj) +{ + gctl::matvec(ax, kernel, x, layout, conj); + return; +} + +int main(int argc, char const *argv[]) +{ + // 生成一组正演解 + gctl::array tmp(2*N); + gctl::random(tmp, 1.0, 2.0, gctl::RdUniform); + + cd_array fm(N); + for (size_t i = 0; i < N; i++) + { + fm[i].real(tmp[2*i]); + fm[i].imag(tmp[2*i + 1]); + } + + ex7 test; + + // 计算共轭梯度B项 + cd_array B(N); + test.cal_partb(fm, B); + + // 声明一组解 + cd_array m(N, std::complex(0.0, 0.0)); + + gctl::clcg_para my_para = test.default_clcg_para(); + my_para.abs_diff = 1; + + test.set_clcg_para(my_para); + + test.CLCG_Minimize(m, B, gctl::CLCG_BICG_SYM); + std::clog << "maximal difference: " << max_diff(fm, m) << std::endl; + return 0; +} \ No newline at end of file diff --git a/example/ex8.cpp b/example/ex8.cpp new file mode 100644 index 0000000..56f8efa --- /dev/null +++ b/example/ex8.cpp @@ -0,0 +1,233 @@ +/******************************************************** + * ██████╗ ██████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗██║ ██║ ██║ + * ██║ ██║██║ ██║ ██║ + * ╚██████╔╝╚██████╗ ██║ ███████╗ + * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ + * 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 . + * + * 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 "gctl/core.h" +#include "gctl/algorithm.h" +#include "../lib/optimization.h" + +#define M 1000 +#define N 900 + +// get random floating points +double random_double(double l, double t) +{ + return (t-l)*rand()*1.0/RAND_MAX + l; +} + +// get random integral numbers +int random_int(int small, int big) +{ + return (rand() % (big - small)) + small; +} + +double max_diff(const gctl::_1d_array &a, const gctl::_1d_array &b) +{ + double max = -1.0; + for (size_t i = 0; i < a.size(); i++) + { + max = std::max(fabs(a[i] - b[i]), max); + } + return max; +} + +class ex8 : public gctl::lbfgs_solver, public gctl::grad_norm +{ +public: + ex8(); + virtual ~ex8(); + virtual double LBFGS_Evaluate(const gctl::_1d_array &x, gctl::_1d_array &g); + virtual int LBFGS_Progress(const gctl::_1d_array &x, const gctl::_1d_array &g, const double fx, + const double converge, const double rate, const gctl::lbfgs_para param, int k, int ls, std::ostream &ss); + + void CalTarget(const gctl::_1d_array &x); + +private: + gctl::_1d_array obs1, obs2, obs3, tmp, grad; + gctl::_2d_matrix k1, k2, k3; +}; + +ex8::ex8() +{ + srand(time(0)); + + tmp.resize(M); + grad.resize(N); + + k1.resize(M, N); + obs1.resize(M); + // 添加一些大数 + int tmp_id, tmp_size; + double tmp_val; + for (int i = 0; i < M; i++) + { + tmp_size = random_int(25, 35); + for (int j = 0; j < tmp_size; j++) + { + tmp_id = random_int(0, N); + tmp_val = random_double(-1.0, 1.0); + + k1[i][tmp_id] = tmp_val; + } + } + + k2.resize(M, N); + obs2.resize(M); + // 添加一些大数 + for (int i = 0; i < M; i++) + { + tmp_size = random_int(25, 35); + for (int j = 0; j < tmp_size; j++) + { + tmp_id = random_int(0, N); + tmp_val = random_double(-200.0, 200.0); + + k2[i][tmp_id] = tmp_val; + } + } + + k3.resize(M, N); + obs3.resize(M); + // 添加一些大数 + for (int i = 0; i < M; i++) + { + tmp_size = random_int(25, 35); + for (int j = 0; j < tmp_size; j++) + { + tmp_id = random_int(0, N); + tmp_val = random_double(-0.01, 0.01); + + k3[i][tmp_id] = tmp_val; + } + } +} + +ex8::~ex8(){} + +double ex8::LBFGS_Evaluate(const gctl::_1d_array &x, gctl::_1d_array &g) +{ + gctl::matvec(tmp, k1, x); + tmp -= obs1; + + gctl::matvec(grad, k1, tmp, gctl::Trans); + gctl::scale(grad, 2.0/M); + + AddSingleLoss(gctl::power2(gctl::module(tmp, gctl::L2))/M, grad); + + gctl::matvec(tmp, k2, x); + tmp -= obs2; + + gctl::matvec(grad, k2, tmp, gctl::Trans); + gctl::scale(grad, 2.0/M); + + AddSingleLoss(gctl::power2(gctl::module(tmp, gctl::L2))/M, grad); + + gctl::matvec(tmp, k3, x); + tmp -= obs3; + + gctl::matvec(grad, k3, tmp, gctl::Trans); + gctl::scale(grad, 2.0/M); + + AddSingleLoss(gctl::power2(gctl::module(tmp, gctl::L2))/M, grad); + + return GradNormLoss(g); +} + +int ex8::LBFGS_Progress(const gctl::_1d_array &x, const gctl::_1d_array &g, const double fx, + const double converge, const double rate, const gctl::lbfgs_para param, int k, int ls, std::ostream &ss) +{ + UpdateWeights(); + + return gctl::lbfgs_solver::LBFGS_Progress(x, g, fx, converge, rate, param, k, ls, ss); +} + +void ex8::CalTarget(const gctl::_1d_array &x) +{ + // 计算正演值 + gctl::matvec(obs1, k1, x); + for (int i = 0; i < M; i++) + { + // 添加噪声 + obs1[i] += random_double(-1e-3, 1e-3); + } + + gctl::matvec(obs2, k2, x); + for (int i = 0; i < M; i++) + { + // 添加噪声 + obs2[i] += random_double(-1e-3, 1e-3); + } + + gctl::matvec(obs3, k3, x); + for (int i = 0; i < M; i++) + { + // 添加噪声 + obs3[i] += random_double(-1e-3, 1e-3); + } + return; +} + +int main(int argc, char const *argv[]) +{ + // 生成一组正演解 + gctl::_1d_array fm(N); + random(fm, 1.0, 2.0, gctl::RdUniform); + + ex8 test; + + // 计算拟合目标项 + test.CalTarget(fm); + + // 声明一组解 + gctl::_1d_array m(N, 0.0); + + gctl::lbfgs_para self_para = test.default_lbfgs_para(); + self_para.linesearch = gctl::LBFGS_LINESEARCH_BACKTRACKING_STRONG_WOLFE; + self_para.epsilon = 1e-6; + + test.set_lbfgs_para(self_para); + test.show_lbfgs_para(); + + test.InitGradNorm(3, N); + test.set_control_weight(1.0); + test.set_weight_step(0.00001); + + double fx = test.LBFGS_Minimize(m); + + std::clog << "maximal difference: " << max_diff(fm, m) << std::endl; + + gctl::_1d_array records; + test.get_records(records); + for (size_t i = 0; i < records.size(); i++) + { + if ((i+1)%3 == 0) + { + std::cout << records[i] << "\n"; + } + else std::cout << records[i] << " "; + } + return 0; +} \ No newline at end of file diff --git a/installer b/installer new file mode 100755 index 0000000..4226075 --- /dev/null +++ b/installer @@ -0,0 +1,54 @@ +#!/bin/bash + +if [[ $# == 0 || ${1} == "help" ]]; then + echo "Compiles executables/libraries and maintains installed files. Two tools 'Cmake' and 'stow' are empolyed here. For more information, see https://cmake.org and https://www.gnu.org/software/stow/." + echo "" + echo "School of Earth Sciences, Zhejiang University" + echo "Yi Zhang (yizhang-geo@zju.edu.cn)" + echo "" + echo "Usage: ./installer [option] [Cmake options]" + echo "" + echo "Options:" + echo "(1) configure: Configure Cmake project(s). This option could take extra Cmake options as in