How would i write a function similiar to np.logspace but with a given first interval and a flexible base?
I am trying to write a function that returns an array that has in the first section linearly increasing elements and in the second section ever increasing distances between the elements. As inputs, I would like to give the starting value, the final value and the length of the array. This would be solveable with np.logspace(). However, I would like the transition from the first section to the second section to be smooth, therefore I had the idea to fix the interval between the first and second element in the second section to the distance between two ajacent elements in the first section. My comprehension is that this is not doable with np.logspace().
Maybe for some context: The goal is to run optimisations in a loop with a running variable, for which the first part of the running variable is more interesting than the last part. Hence I am looking for an elegant way to devide the sections smoothly.
I read and tried a lot but could not find a ready to use solution. My approach was hence the following, where x_n is the final value, n is the length of the array and x_0 describes the threshold between the two sections.
The underlying function for the exponential section is the following:
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))
This mostly works as intended. However, if I now increase or decrease the threshold x_0 to a large or small value relative to x_n (e.g. 15/0.005 for n = 100) this approach does work anymore. In my application, this is mostly and issue for the case that x_0 is small compared to x_n, so I tried to square i in the function, which however does not achieve the desired result.
Is there an easy fix to my problem or are there other solutions to achieve the desired result?
Thanks in advance!
There are several changes in order:
- You should frame the first segment as a half-open interval, including 0 but excluding the threshold.
- The second segment should include both the threshold and the endpoint.
- You should enforce the threshold, endpoint, and a constraint of first-order differential continuity.
Because of the exponential, this is very sensitive to initial conditions and a poor initial guess will not solve. This works for the order of magnitude of your example inputs.
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]
To do better, you need to solve with a Lambert W:
import numpy as np
import scipy.special
from matplotlib import pyplot as plt
def exp_function(a: float, b: float, c: float, i: np.ndarray) -> np.ndarray:
return a*b**i + c
def make_series(
n: int, # length of array
x_0: float, # threshold
x_n: float, # final value
) -> 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)
# Exponential region with decreasing resolution
# Closed interval: [x_0, n]
i_n = n - n_half - 1
'''
d0 = alnb # 1: gradient continuity
x0 = a + c # 2: threshold continuity
xn = ab^in + c # 3: endpoint
x0-a = xn-ab^in # 2,3 for c
a(b^in - 1) = xn - x0
b^in = 1 + (xn - x0)/a
b = (1 + (xn - x0)/a)^(1/in)
d0 = aln( (1 + (xn - x0)/a)^(1/in) ) # 1 for b
exp(d0in/a) = 1 + (xn - x0)/a
'''
p = first_step*i_n
q = x_n - x_0
a = -p*q/(
p + q*scipy.special.lambertw(
-p/q*np.exp(-p/q),
k=-1, tol=1e-16,
).real
)
b = (1 + q/a)**(1/i_n)
c = x_0 - a
print(f'a={a:.3e}, b={b:.3f}, c={c:.3e}')
print(f'Threshold: {x_0} ~ {a + c:.3e}')
print(f'Endpoint: {x_n} ~ {a*b**i_n + c:.3f}')
exp_section = exp_function(a, b, c, np.arange(1 + i_n))
return np.concatenate((lin_section, exp_section))
def demo() -> None:
n = 100
series = make_series(x_0=5e-3, x_n=15, n=n)
fig, ax = plt.subplots()
ax.semilogy(np.arange(n), series)
plt.show()
if __name__ == '__main__':
demo()
a=4.728e-04, b=1.236, c=4.527e-03
Threshold: 0.005 ~ 5.000e-03
Endpoint: 15 ~ 15.000