Улучшение времени выполнения алгоритма аналитической трассировки лучей. ⇐ Python

[Anonymous]

Справочная информация
Я написал класс Python, предназначенный для вычисления времени, необходимого лучу света для распространения между двумя точками (

Код: Выделить всё

init_point

и term_point) в сложной среде (смоделированной с помощью параметров medium_model = [A, B, C]).
Чтобы вычислить этот свет время прохождения, нам нужны две вещи:

• Вышеупомянутые начальная и конечная точки.
• Угол (

  Код: Выделить всё

  launch_angle

  ), который луч света образует по оси X в начальной точке.

У меня есть выражение для координаты z луч в терминах координат x и y (которые мы рассматриваем как одну координату «x» (

Код: Выделить всё

x = sqrt(x**2 + y**2)

), launch_angle, а также начальную и конечную точки, поэтому мы можем определить последнюю, просто оптимизируя разницу между известной z-координатой терминала и тем, что мы получаем для различных значений launch_angle >.

Код

Код: Выделить всё

from dataclasses import dataclass
import warnings
import numpy as np

# Suppress warnings (sometimes there is overflow in _calculate_z_coord, etc.)
warnings.filterwarnings("ignore", category=RuntimeWarning)

@dataclass
class RayTracer:
medium_model: np.ndarray
SPEED_OF_LIGHT = 299792458  # Speed of light in m/s

def _calculate_z_coord(self, x: np.ndarray, launch_angle: np.ndarray, x0: np.ndarray, z0: np.ndarray) -> np.ndarray:
"""Calculate the z-coordinate based on launch angle and other parameters."""
A, B, C = self.medium_model

exp_Cz0 = np.exp(C * z0)
cos_launch_angle = np.cos(launch_angle)
beta = (A - B * exp_Cz0) * cos_launch_angle

sqrt_A2_beta2 = np.sqrt(A**2 - beta**2)
K = C * sqrt_A2_beta2 / beta

term1 = A**2 - beta**2
term2 = A * B * exp_Cz0
sqrt_term = np.sqrt(term1 + 2 * term2 + B**2 * exp_Cz0**2)

# Precompute for efficiency
log_arg = term1 + term2 + sqrt_A2_beta2 * sqrt_term
t = (sqrt_A2_beta2 * C * x0 - beta * C * z0 + beta * np.log(log_arg)) / (sqrt_A2_beta2 * C)

exp_Kx = np.exp(K * x)
exp_Kt = np.exp(K * t)
log_term_num = 2 * term1 * np.exp(K * (t + x))
log_term_den = beta**2 * B**2 * exp_Kx**2 - 2 * A * B * exp_Kt + exp_Kt**2
log_term = log_term_num / log_term_den

return (1 / C) * np.log(log_term)

def _find_launch_angle(self, init_points: np.ndarray, term_points: np.ndarray, num_steps: int = 1000) -> np.ndarray:
"""Find the optimal launch angle."""
x0 = np.hypot(init_points[:, 0], init_points[:, 1])
x1 = np.hypot(term_points[:, 0], term_points[:, 1])

# Coarse search with a fine search
launch_angles = np.linspace(-np.pi, np.pi, num_steps)

# Precompute z_coords for all launch angles
z_coords = self._calculate_z_coord(x1[:, np.newaxis], launch_angles, x0[:, np.newaxis], init_points[:, 2, np.newaxis])

term_z = term_points[:, 2][:, np.newaxis]
objective_values = (z_coords - term_z)**2

# Find the best angles
best_indices = np.nanargmin(objective_values, axis=1)
best_angles = launch_angles[best_indices]

return best_angles

def transit_time(self, init_points: np.ndarray, term_points: np.ndarray) -> np.ndarray:
"""Calculate the transit time."""
A, B, C = self.medium_model

# Vectorized launch angle search
launch_angles = self._find_launch_angle(init_points, term_points)

exp_Cz_init = np.exp(C * init_points[:, 2])
beta = np.abs((A - B * exp_Cz_init) * np.cos(launch_angles))
sqrt_A2_beta2 = np.sqrt(A**2 - beta**2)
K = C * sqrt_A2_beta2 / beta

def time_expression(z: np.ndarray, beta: np.ndarray, K: np.ndarray) ->  np.ndarray:
exp_Cz = np.exp(C * z)
t = np.sqrt((A + B * exp_Cz - beta) / (A + B * exp_Cz + beta))
alpha = np.sqrt((A - beta) / (A + beta))
log_expr = np.log(np.abs((t - alpha) / (t + alpha)))
return A * np.sqrt((C**2 + K**2) / (C**2 * K**2)) * log_expr

time_diff = time_expression(term_points[:, 2], beta, K) - time_expression(init_points[:, 2], beta, K)

return time_diff / self.SPEED_OF_LIGHT

# Example usage
if __name__ == "__main__":
ray_tracer = RayTracer(np.array([1.78, 0.454, 0.0132]))

# Generate random init_points and term_points
num_points = 100000
init_points = np.random.uniform(low=-50, high=50, size=(num_points, 3))
term_points = np.random.uniform(low=-50, high=50, size=(num_points, 3))

# I'm including a simple but rough measure of the execution time using `time` here, for your reference.
import time

start_time = time.time()
transit_times = ray_tracer.transit_time(init_points, term_points)
end_time = time.time()

# Print results
print("Transit times:", transit_times)
print(f"Elapsed time: {end_time - start_time:.6f} seconds")

Код векторизован. Инициализация RayTracer выполняется следующим образом:

Код: Выделить всё

my_ray_tracer = RayTracer(np.array([A, B, C]))

Функция Transition_time() векторизована и принимает массив начальных и конечных точек, например:

Код: Выделить всё

init_points = np.random.uniform(low=-50, high=50, size=(num_points, 3))
term_points = np.random.uniform(low=-50, high=50, size=(num_points, 3))

transit_times = ray_tracer.transit_time(init_points, term_points)

Моя «проблема»
Мне в некоторых приложениях хотелось бы иметь возможность быстро вычислять время распространения тысячи, десятки тысяч или даже сотни тысяч лучей. Хотя этот код довольно быстр, для некоторых моих приложений он еще недостаточно быстр. Я применил все знакомые мне приемы (например, векторизацию) и значительно сократил среднее время выполнения, но все же недостаточно значительно, чтобы это было полезно для того, что мне нужно.
Мои требования следующие:

• Меня не особо волнует эффективность использования памяти, только время выполнения transient_time().
• Я добавляю это в базу кода, зависимостями которой являются только NumPy и SciPy. Я не могу использовать какие-либо дополнительные внешние библиотеки.

Как я могу еще больше улучшить время выполнения этого кода?

Подробнее здесь: https://stackoverflow.com/questions/790 ... -algorithm