Skip to content

Commit

Permalink
🔧 数值计算模块及最优化算法库使用 std::valarray 代替 std::vector
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhaoxi-scut
zhaoxi-scut committed Dec 30, 2024
1 parent 81f3f9a commit 8a07c2a
Show file tree
Hide file tree
Showing 8 changed files with 113 additions and 276 deletions.
8 changes: 4 additions & 4 deletions extra/compensator/src/gravity_compensator/gravity_compensator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,10 +35,10 @@ GravityCompensator::Impl::Impl() : _yaw_static_com(para::gravity_compensator_par
double a42 = -a24, a44 = a22;

Odes fs(4);
fs[0] = [](double, const std::vector<double> &x) { return x[1]; };
fs[1] = [=](double, const std::vector<double> &x) { return a22 * x[1] + a24 * x[3]; };
fs[2] = [](double, const std::vector<double> &x) { return x[3]; };
fs[3] = [=](double, const std::vector<double> &x) { return a42 * x[1] + a44 * x[3] - gravity_compensator_param.g; };
fs[0] = [](double, const std::valarray<double> &x) { return x[1]; };
fs[1] = [=](double, const std::valarray<double> &x) { return a22 * x[1] + a24 * x[3]; };
fs[2] = [](double, const std::valarray<double> &x) { return x[3]; };
fs[3] = [=](double, const std::valarray<double> &x) { return a42 * x[1] + a44 * x[3] - gravity_compensator_param.g; };
_rk = std::make_unique<RungeKutta2>(fs);
}

Expand Down
150 changes: 0 additions & 150 deletions modules/algorithm/include/rmvl/algorithm/math.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -402,156 +402,6 @@ typename ForwardIterator::value_type calculateModeNum(ForwardIterator first, For
->first;
}

/**
* @brief 使用 `std::vector` 表示的向量加法
*
* @tparam T 数据类型
* @param[in] vec1 向量 1
* @param[in] vec2 向量 2
* @return 和向量
*/
template <typename T>
inline std::vector<T> operator+(const std::vector<T> &vec1, const std::vector<T> &vec2)
{
std::vector<T> retval(vec1.size());
std::transform(vec1.cbegin(), vec1.cend(), vec2.cbegin(), retval.begin(), std::plus<T>());
return retval;
}

/**
* @brief 使用 `std::vector` 表示的向量减法
*
* @tparam T 数据类型
* @param[in] vec1 向量 1
* @param[in] vec2 向量 2
* @return 差向量
*/
template <typename T>
inline std::vector<T> operator-(const std::vector<T> &vec1, const std::vector<T> &vec2)
{
std::vector<T> retval(vec1.size());
std::transform(vec1.cbegin(), vec1.cend(), vec2.cbegin(), retval.begin(), std::minus<T>());
return retval;
}

/**
* @brief 使用 `std::vector` 表示的向量自加
*
* @tparam T 数据类型
* @param[in] vec1 向量 1
* @param[in] vec2 向量 2
* @return 和向量
*/
template <typename T>
inline std::vector<T> &operator+=(std::vector<T> &vec1, const std::vector<T> &vec2)
{
std::transform(vec1.cbegin(), vec1.cend(), vec2.cbegin(), vec1.begin(), std::plus<T>());
return vec1;
}

/**
* @brief 使用 `std::vector` 表示的向量自减
*
* @tparam T 数据类型
* @param[in] vec1 向量 1
* @param[in] vec2 向量 2
* @return 差向量
*/
template <typename T>
inline std::vector<T> &operator-=(std::vector<T> &vec1, const std::vector<T> &vec2)
{
std::transform(vec1.cbegin(), vec1.cend(), vec2.cbegin(), vec1.begin(), std::minus<T>());
return vec1;
}

/**
* @brief 使用 `std::vector` 表示的向量取反
*
* @tparam T 数据类型
* @param[in] vec 向量
* @return 向量取反后的结果
*/
template <typename T>
inline std::vector<T> operator-(const std::vector<T> &vec)
{
std::vector<T> retval(vec.size());
std::transform(vec.cbegin(), vec.cend(), retval.begin(), std::negate<T>());
return retval;
}

/**
* @brief 使用 `std::vector` 表示的向量乘法（数乘）
*
* @tparam T 数据类型
* @param[in] vec 向量
* @param[in] val 数乘因子
* @return 向量乘法（数乘）
*/
template <typename T>
inline std::vector<T> operator*(const std::vector<T> &vec, T val)
{
std::vector<T> retval(vec.size());
std::transform(vec.cbegin(), vec.cend(), retval.begin(), [val](const T &x) { return x * val; });
return retval;
}

/**
* @brief 使用 `std::vector` 表示的向量乘法（数乘）
*
* @tparam T 数据类型
* @param[in] val 数乘因子
* @param[in] vec 向量
* @return 向量乘法（数乘）
*/
template <typename T>
inline std::vector<T> operator*(T val, const std::vector<T> &vec) { return vec * val; }

/**
* @brief 使用 `std::vector` 表示的向量乘法（数乘）
*
* @tparam T 数据类型
* @param[in] vec 向量
* @param[in] val 数乘因子
* @return 向量乘法（数乘）
*/
template <typename T>
inline std::vector<T> &operator*=(std::vector<T> &vec, T val)
{
std::transform(vec.cbegin(), vec.cend(), vec.begin(), [val](const T &x) { return x * val; });
return vec;
}

/**
* @brief 使用 `std::vector` 表示的向量除法（数乘）
*
* @tparam T 数据类型
* @param[in] vec 向量
* @param[in] val 除数
* @return 向量除法（数乘）
*/
template <typename T>
inline std::vector<T> operator/(const std::vector<T> &vec, T val)
{
std::vector<T> retval(vec.size());
std::transform(vec.cbegin(), vec.cend(), retval.begin(), [val](const T &x) { return x / val; });
return retval;
}

/**
* @brief 使用 `std::vector` 表示的向量除法（数乘）
*
* @tparam T 数据类型
* @param[in] val 除数
* @param[in] vec 向量
* @return 向量除法（数乘）
*/
template <typename T>
inline std::vector<T> &operator/=(std::vector<T> &vec, T val)
{
std::transform(vec.cbegin(), vec.cend(), vec.begin(), [val](const T &x) { return x / val; });
return vec;
}

// ------------------------【数学模型算法】------------------------

//! 熵权 TOPSIS 算法
Expand Down
41 changes: 21 additions & 20 deletions modules/algorithm/include/rmvl/algorithm/numcal.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include <bitset>
#include <cstdint>
#include <functional>
#include <valarray>
#include <vector>

#if __cplusplus >= 202302L
Expand Down Expand Up @@ -182,9 +183,9 @@ class RMVL_EXPORTS_W NonlinearSolver
////////////// 常微分方程（组）数值解 //////////////

//! 常微分方程
using Ode = std::function<double(double, const std::vector<double> &)>;
using Ode = std::function<double(double, const std::valarray<double> &)>;
//! 常微分方程组
using Odes = std::vector<std::function<double(double, const std::vector<double> &)>>;
using Odes = std::vector<std::function<double(double, const std::valarray<double> &)>>;

/**
* @brief Butcher 表形式的常微分方程（组）数值求解器
Expand All @@ -204,23 +205,23 @@ class RMVL_EXPORTS_W RungeKutta
* @param[in] lam Butcher 表 \f$\pmb\lambda\f$ 向量
* @param[in] r Butcher 表 \f$R\f$ 矩阵
*/
RMVL_W RungeKutta(const Odes &fs, const std::vector<double> &p, const std::vector<double> &lam, const std::vector<std::vector<double>> &r);
RMVL_W RungeKutta(const Odes &fs, const std::valarray<double> &p, const std::valarray<double> &lam, const std::valarray<std::valarray<double>> &r);

/**
* @brief 设置常微分方程（组）的初值
*
* @param[in] t0 初始位置的自变量 \f$t_0\f$
* @param[in] x0 初始位置的因变量 \f$\pmb x(t_0)\f$
*/
RMVL_W inline void init(double t0, const std::vector<double> &x0) { _t0 = t0, _x0 = x0; }
RMVL_W inline void init(double t0, const std::valarray<double> &x0) { _t0 = t0, _x0 = x0; }

/**
* @brief 设置常微分方程（组）的初值
*
* @param[in] t0 初始位置的自变量 \f$t_0\f$
* @param[in] x0 初始位置的因变量 \f$\pmb x(t_0)\f$
*/
inline void init(double t0, std::vector<double> &&x0) { _t0 = t0, _x0 = std::move(x0); }
inline void init(double t0, std::valarray<double> &&x0) { _t0 = t0, _x0 = std::move(x0); }

/**
* @brief 计算常微分方程（组）的数值解
Expand All @@ -230,7 +231,7 @@ class RMVL_EXPORTS_W RungeKutta
*
* @return 从初始位置开始迭代 \f$n\f$ 次后共 \f$n+1\f$ 个数值解，自变量可通过 \f$t_0+ih\f$ 计算得到
*/
RMVL_W std::vector<std::vector<double>> solve(double h, std::size_t n);
RMVL_W std::vector<std::valarray<double>> solve(double h, std::size_t n);

#if __cpp_lib_generator >= 202207L
/**
Expand All @@ -240,22 +241,22 @@ class RMVL_EXPORTS_W RungeKutta
* @param[in] n 迭代次数
* @return 从初始位置开始迭代计算的生成器，初值不会被 `co_yield`，共生成 \f$n\f$ 个数值解
*/
std::generator<std::vector<double>> generate(double h, std::size_t n);
std::generator<std::valarray<double>> generate(double h, std::size_t n);
#endif

private:
// 加权系数，`k[i][j]`: `i` 表示第 `i` 个加权系数组，`j` 表示来自第 `j` 条方程
std::vector<std::vector<double>> _ks;
std::vector<std::valarray<double>> _ks;

protected:
//! 一阶常微分方程组的函数对象 \f$\dot{\pmb x}=\pmb F(t, \pmb x)\f$
Odes _fs;
double _t0; //!< 初值的自变量 \f$t\f$
std::vector<double> _x0; //!< 初值的因变量 \f$\pmb x(t)\f$
double _t0; //!< 初值的自变量 \f$t\f$
std::valarray<double> _x0; //!< 初值的因变量 \f$\pmb x(t)\f$

std::vector<double> _p; //!< Butcher 表 \f$\pmb p\f$ 向量
std::vector<double> _lambda; //!< Butcher 表 \f$\pmb\lambda\f$ 向量
std::vector<std::vector<double>> _r; //!< Butcher 表 \f$R\f$ 矩阵
std::valarray<double> _p; //!< Butcher 表 \f$\pmb p\f$ 向量
std::valarray<double> _lambda; //!< Butcher 表 \f$\pmb\lambda\f$ 向量
std::valarray<std::valarray<double>> _r; //!< Butcher 表 \f$R\f$ 矩阵
};

//! 2 阶 2 级 Runge-Kutta 求解器
Expand Down Expand Up @@ -304,9 +305,9 @@ using Func1d = std::function<double(double)>;
//! 一元函数组
using Func1ds = std::vector<std::function<double(double)>>;
//! 多元函数
using FuncNd = std::function<double(const std::vector<double> &)>;
using FuncNd = std::function<double(const std::valarray<double> &)>;
//! 多元函数组
using FuncNds = std::vector<std::function<double(const std::vector<double> &)>>;
using FuncNds = std::vector<std::function<double(const std::valarray<double> &)>>;

//! 梯度/导数计算模式
enum class DiffMode : uint8_t
Expand Down Expand Up @@ -362,7 +363,7 @@ RMVL_EXPORTS_W double derivative(Func1d func, double x, DiffMode mode = DiffMode
* @param[in] dx 计算偏导数时，坐标的微小增量，默认为 `1e-3`
* @return 函数在指定点的梯度向量
*/
RMVL_EXPORTS_W std::vector<double> grad(FuncNd func, const std::vector<double> &x, DiffMode mode = DiffMode::Central, double dx = 1e-3);
RMVL_EXPORTS_W std::valarray<double> grad(FuncNd func, const std::valarray<double> &x, DiffMode mode = DiffMode::Central, double dx = 1e-3);

/**
* @brief 采用进退法确定搜索区间
Expand Down Expand Up @@ -392,7 +393,7 @@ RMVL_EXPORTS_W std::pair<double, double> fminbnd(Func1d func, double x1, double
* @param[in] options 优化选项，可供设置的有 `fmin_mode`、`max_iter`、`tol` 和 `dx`
* @return `[x, fval]` 最小值点和最小值
*/
RMVL_EXPORTS_W std::pair<std::vector<double>, double> fminunc(FuncNd func, const std::vector<double> &x0, const OptimalOptions &options = {});
RMVL_EXPORTS_W std::pair<std::valarray<double>, double> fminunc(FuncNd func, const std::valarray<double> &x0, const OptimalOptions &options = {});

/**
* @brief 有约束多维函数的最小值搜索
Expand All @@ -404,7 +405,7 @@ RMVL_EXPORTS_W std::pair<std::vector<double>, double> fminunc(FuncNd func, const
* @param[in] options options 优化选项，可供设置的有 `exterior`、`fmin_mode`、`max_iter`、`tol` 和 `dx`
* @return `[x, fval]` 最小值点和最小值
*/
RMVL_EXPORTS_W std::pair<std::vector<double>, double> fmincon(FuncNd func, const std::vector<double> &x0, FuncNds c, FuncNds ceq, const OptimalOptions &options = {});
RMVL_EXPORTS_W std::pair<std::valarray<double>, double> fmincon(FuncNd func, const std::valarray<double> &x0, FuncNds c, FuncNds ceq, const OptimalOptions &options = {});

/**
* @brief 非线性最小二乘求解，实现与 \cite Agarwal23 类似的算法
Expand All @@ -415,7 +416,7 @@ RMVL_EXPORTS_W std::pair<std::vector<double>, double> fmincon(FuncNd func, const
* @param[in] options 优化选项，可供设置的有 `lsq_mode`、`max_iter`、`tol` 和 `dx`
* @return 最小二乘解
*/
RMVL_EXPORTS_W std::vector<double> lsqnonlin(const FuncNds &funcs, const std::vector<double> &x0, const OptimalOptions &options = {});
RMVL_EXPORTS_W std::valarray<double> lsqnonlin(const FuncNds &funcs, const std::valarray<double> &x0, const OptimalOptions &options = {});

//! Robust 核函数
enum class RobustMode : uint8_t
Expand All @@ -436,7 +437,7 @@ enum class RobustMode : uint8_t
* @param[in] options 优化选项，可供设置的有 `lsq_mode`、`max_iter`、`tol` 和 `dx`
* @return 最小二乘解
*/
RMVL_EXPORTS_W std::vector<double> lsqnonlinRKF(const FuncNds &funcs, const std::vector<double> &x0, RobustMode rb, const OptimalOptions &options = {});
RMVL_EXPORTS_W std::valarray<double> lsqnonlinRKF(const FuncNds &funcs, const std::valarray<double> &x0, RobustMode rb, const OptimalOptions &options = {});

//! @} algorithm_optimal

Expand Down
12 changes: 6 additions & 6 deletions modules/algorithm/perf/perf_optimal.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ namespace rm_test

/////////////////////// Quadratic Function ///////////////////////

static double quadraticFunc(const std::vector<double> &x)
static double quadraticFunc(const std::valarray<double> &x)
{
const auto &x1 = x[0], &x2 = x[1];
return 60 - 10 * x1 - 4 * x2 + x1 * x1 + x2 * x2 - x1 * x2;
Expand Down Expand Up @@ -57,9 +57,9 @@ static void cg_quadratic_cv(benchmark::State &state)
BENCHMARK(cg_quadratic_rmvl)->Name("fminunc (conj_grad, quadratic) - by rmvl ")->Iterations(50);
BENCHMARK(cg_quadratic_cv)->Name("fminunc (conj_grad, quadratic) - by opencv")->Iterations(50);

static inline double cle1(const std::vector<double> &x) { return -x[0] - x[1] + 10; }
static inline double cle2(const std::vector<double> &x) { return 2 * x[0] + x[1] - 30; }
static inline double cle3(const std::vector<double> &x) { return -x[0] + x[1] - 5; }
static inline double cle1(const std::valarray<double> &x) { return -x[0] - x[1] + 10; }
static inline double cle2(const std::valarray<double> &x) { return 2 * x[0] + x[1] - 30; }
static inline double cle3(const std::valarray<double> &x) { return -x[0] + x[1] - 5; }

static void com_quadratic_rmvl(benchmark::State &state)
{
Expand All @@ -71,7 +71,7 @@ BENCHMARK(com_quadratic_rmvl)->Name("fmincon (conj_grad, quadratic) - by rmvl "

/////////////////////// Rosenbrock Function ///////////////////////

static double rosenbrockFunc(const std::vector<double> &x)
static double rosenbrockFunc(const std::valarray<double> &x)
{
const auto &x1 = x[0], &x2 = x[1];
return 100 * (x2 - x1 * x1) * (x2 - x1 * x1) + (1 - x1) * (1 - x1);
Expand Down Expand Up @@ -127,7 +127,7 @@ void lsqnonlin_rmvl(benchmark::State &state)
{
rm::FuncNds lsq_sine(20);
for (std::size_t i = 0; i < lsq_sine.size(); ++i)
lsq_sine[i] = [=](const std::vector<double> &x) { return x[0] * std::sin(x[1] * i + x[2]) + x[3] - real_f(i); };
lsq_sine[i] = [=](const std::valarray<double> &x) { return x[0] * std::sin(x[1] * i + x[2]) + x[3] - real_f(i); };

auto x = rm::lsqnonlin(lsq_sine, {1, 0.02, 0, 1.09});
}
Expand Down
Loading

0 comments on commit 8a07c2a

Please sign in to comment.