我如何编写一个类似于 np.logspace 的函数，但具有给定的第一个间隔和灵活的基数？

Question

我正在尝试编写一个函数，该函数返回一个数组，该数组在第一部分中线性增加元素，在第二部分中元素之间的距离不断增加。作为输入，我想给出起始值、最终值和数组的长度。这可以通过 np.logspace() 来解决。但是，我希望从第一部分到第二部分的过渡能够平滑，因此我有想法将第二部分中的第一个和第二个元素之间的间隔固定为第一部分中两个相邻元素之间的距离。我的理解是，这对于 np.logspace() 是不可能的。

也许对于某些上下文：目标是在带有运行变量的循环中运行优化，其中运行变量的第一部分比最后部分更有趣。因此，我正在寻找一种优雅的方式来顺利地划分各个部分。

我阅读并尝试了很多，但找不到现成的解决方案。因此，我的方法如下，其中 x_n 是最终值，n 是数组的长度，x_0 描述两个部分之间的阈值。

指数部分的基础函数如下：

x_i = x_0 + b*(y^i-1)

def list_creation(x_n, n, x_0, y_initial_guess = 1.1):
    
    n_half = int((n + 1) // 2)  # Adjust for odd/even cases

    # Linear difference, region with detailed resolution
    first_section = np.linspace(0,x_0, n_half)
    first_step = first_section[1] - first_section[0]

    # Exponention decreasing resolution section
    x_0 =  x_0 + first_step         # This is necessary to avoid doubled values

    # Function to solve for 
    def equation_for_y(y):
        return x_0 + first_step * (y**(n_half-1) - 1) / (y - 1) - x_n

    y_solution = fsolve(equation_for_y, y_initial_guess)[0]

    # Calculating scaling factor b
    b = first_step / (y_solution - 1)

    # Calculating array given the calculated variables
    second_section = np.array([x_0 + b * (y_solution**i - 1) for i in range(n_half)])

    return np.concatenate((first_section, second_section))

这基本上可以按预期工作。但是，如果我现在将阈值 x_0 增加或减少到相对于 x_n 较大或较小的值（例如，n = 100 时为 15/0.005），则这种方法不再有效。在我的应用程序中，这主要是 x_0 与 x_n 相比较小的情况下的问题，因此我尝试在函数中对 i 进行平方，但这并没有达到预期的结果。

有一个简单的方法可以解决我的问题，还是有其他解决方案可以达到预期的结果？

提前致谢！

Answer 1

顺序有几处变化：

• 您应该将第一段构建为半开区间，包括 0 但不包括阈值。
• 第二段应包括阈值和端点。
• 您应该强制执行阈值、端点、和一阶微分连续性的约束。

由于指数，这对初始条件非常敏感，糟糕的初始猜测将无法解决。这适用于示例输入的数量级。

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import fsolve, least_squares


def exp_function(a: float, b: float, c: float, i: np.ndarray) -> np.ndarray:
    return a*b**i + c


def b_to_abc(
    b: float,
    x_0: float,
    x_n: float,
    n_half: int,
    n: int,
    first_step: float,
) -> tuple[float, float, float]:
    # gradient constraint
    # a = first_step/np.log(b_est)

    # endpoint constraint
    # x0 - a == xn - a*b**(n - nhalf - 1)
    a = (x_n - x_0)/(b**(n - n_half - 1) - 1)

    # threshold constraint
    c = x_0 - a

    return a, b, c


def exp_equations(
    b: float,
    x_0: float,
    x_n: float,
    n_half: int,
    n: int,
    first_step: float,
) -> float:
    a, b, c = b_to_abc(b=b, x_0=x_0, x_n=x_n, n_half=n_half, n=n, first_step=first_step)
    # dx/di = a*ln(b) * b^i: gradient starts at same value of linear section
    differential_error = np.exp(first_step/a) - b
    return differential_error


def print_params(
    b: float, args: tuple[int | float, ...]) -> None:
    x_0, x_n, n_half, n, first_step = args
    abc = b_to_abc(b, *args)
    a, b, c = abc

    print('abc =', abc)
    print(f'threshold: {x_0} ~ {exp_function(*abc, 0)}')
    print(f'endpoint: {x_n} ~ {exp_function(*abc, n - n_half - 1)}')
    print(f'differential: {first_step} ~ {a*np.log(b)}')
    print('differential error =', exp_equations(b, *args))


def make_series(
    n: int,      # length of array
    x_0: float,  # threshold
    x_n: float,  # final value
    use_fsolve: bool = False,
) -> np.ndarray:
    n_half = n//2
    first_step = x_0/n_half

    # Linear region with detailed resolution
    # Half-open interval: [0, x_0)
    lin_section = np.linspace(start=0, stop=x_0*(n_half - 1)/n_half, num=n_half)

    # Analytic solution would require a lambert W. Do the easier thing and call a solver.
    b_est = 1.2
    args = x_0, x_n, n_half, n, first_step
    print('Estimate:')
    print_params(b_est, args)
    print()

    if use_fsolve:
        (b,), infodict, success, message = fsolve(
            func=exp_equations, x0=b_est, args=args, maxfev=10000, full_output=True,
        )
        assert success == 1, message
    else:
        result = least_squares(
            fun=exp_equations, x0=b_est, args=args, max_nfev=10000,
        )
        assert result.success, result.message
        b, = result.x

    print('Fit:')
    print_params(b, args)
    print()

    # Exponential region with decreasing resolution
    # Closed interval: [x_0, n]
    abc = b_to_abc(b, *args)
    exp_section = exp_function(*abc, np.arange(n - n_half))

    return np.concatenate((lin_section, exp_section))


def demo() -> None:
    n = 100
    series = make_series(x_0=5e-3, x_n=15, n=n, use_fsolve=False)
    print(series)
    fig, ax = plt.subplots()
    ax.semilogy(np.arange(n), series)
    plt.show()


if __name__ == '__main__':
    demo()

Estimate:
abc = (0.0019775281955880866, 1.2, 0.0030224718044119135)
threshold: 0.005 ~ 0.005
endpoint: 15 ~ 15.0
differential: 0.0001 ~ 0.00036054601922355986
differential error = -0.14813142362076936

Fit:
abc = (0.00047275979928614655, 1.2355595042735252, 0.004527240200713854)
threshold: 0.005 ~ 0.005
endpoint: 15 ~ 14.999999999999998
differential: 0.0001 ~ 9.99999999999992e-05
differential error = 1.9984014443252818e-15

[0.00000000e+00 1.00000000e-04 2.00000000e-04 3.00000000e-04
 4.00000000e-04 5.00000000e-04 6.00000000e-04 7.00000000e-04
 8.00000000e-04 9.00000000e-04 1.00000000e-03 1.10000000e-03
 1.20000000e-03 1.30000000e-03 1.40000000e-03 1.50000000e-03
 1.60000000e-03 1.70000000e-03 1.80000000e-03 1.90000000e-03
 2.00000000e-03 2.10000000e-03 2.20000000e-03 2.30000000e-03
 2.40000000e-03 2.50000000e-03 2.60000000e-03 2.70000000e-03
 2.80000000e-03 2.90000000e-03 3.00000000e-03 3.10000000e-03
 3.20000000e-03 3.30000000e-03 3.40000000e-03 3.50000000e-03
 3.60000000e-03 3.70000000e-03 3.80000000e-03 3.90000000e-03
 4.00000000e-03 4.10000000e-03 4.20000000e-03 4.30000000e-03
 4.40000000e-03 4.50000000e-03 4.60000000e-03 4.70000000e-03
 4.80000000e-03 4.90000000e-03 5.00000000e-03 5.11136306e-03
 5.24895876e-03 5.41896642e-03 5.62902101e-03 5.88855595e-03
 6.20922681e-03 6.60543474e-03 7.09497322e-03 7.69982714e-03
 8.44716014e-03 9.37053454e-03 1.05114186e-02 1.19210486e-02
 1.36627305e-02 1.58146821e-02 1.84735463e-02 2.17587312e-02
 2.58177727e-02 3.08329600e-02 3.70295223e-02 4.46857437e-02
 5.41454609e-02 6.58335044e-02 8.02747776e-02 9.81178299e-02
 1.20163983e-01 1.47403317e-01 1.81059134e-01 2.22642900e-01
 2.74022116e-01 3.37504196e-01 4.15940083e-01 5.12852288e-01
 6.32593084e-01 7.80539963e-01 9.63337135e-01 1.18919392e+00
 1.46825341e+00 1.81304803e+00 2.23906229e+00 2.76542825e+00
 3.41578473e+00 4.21933885e+00 5.21217778e+00 6.43888936e+00
 7.95456452e+00 9.82727136e+00 1.21411121e+01 1.50000000e+01]

为了做得更好，您需要使用 Lambert W 来求解。我稍后会看看是否可以用这种方法生成一个程序。