Source code for ctapipe.reco.impact

#!/usr/bin/env python3

import math

import numpy as np
import as ma
from astropy import units as u
from astropy.coordinates import SkyCoord, AltAz
from iminuit import Minuit
from scipy.optimize import minimize, least_squares
from scipy.stats import norm

from ctapipe.core import Component
from ctapipe.coordinates import (
from ctapipe.image import neg_log_likelihood, mean_poisson_likelihood_gaussian
from ctapipe.instrument import get_atmosphere_profile_functions
from ctapipe.containers import (
from ctapipe.utils.template_network_interpolator import (

__all__ = ["ImPACTReconstructor", "energy_prior", "xmax_prior", "guess_shower_depth"]

[docs]def guess_shower_depth(energy): """ Simple estimation of depth of shower max based on the expected gamma-ray elongation rate. Parameters ---------- energy: float Energy of the shower in TeV Returns ------- float: Expected depth of shower maximum """ x_max_exp = 300 + 93 * np.log10(energy) return x_max_exp
[docs]def energy_prior(energy, index=-1): return -2 * np.log(energy**index)
[docs]def xmax_prior(energy, xmax, width=100): x_max_exp = guess_shower_depth(energy) diff = xmax - x_max_exp return -2 * np.log(norm.pdf(diff / width))
[docs]class ImPACTReconstructor(Component): """This class is an implementation if the impact_reco Monte Carlo Template based image fitting method from parsons14. This method uses a comparision of the predicted image from a library of image templates to perform a maximum likelihood fit for the shower axis, energy and height of maximum. Because this application is computationally intensive the usual advice to use astropy units for all quantities is ignored (as these slow down some computations), instead units within the class are fixed: - Angular units in radians - Distance units in metres - Energy units in TeV References ---------- .. [parsons14] Parsons & Hinton, Astroparticle Physics 56 (2014), pp. 26-34 """ # For likelihood calculation we need the with of the # pedestal distribution for each pixel # currently this is not availible from the calibration, # so for now lets hard code it in a dict ped_table = { "LSTCam": 2.8, "NectarCam": 2.3, "FlashCam": 2.3, "CHEC": 0.5, "DUMMY": 0, } spe = 0.5 # Also hard code single p.e. distribution width def __init__( self, root_dir=".", minimiser="minuit", prior="", template_scale=1.0, xmax_offset=0, use_time_gradient=False, ): """ Create a new instance of ImPACTReconstructor """ # First we create a dictionary of image template interpolators # for each telescope type self.root_dir = root_dir self.priors = prior self.minimiser_name = minimiser self.file_names = { "CHEC": ["GCT_05deg_ada.template.gz", "GCT_05deg_time.template.gz"], "LSTCam": ["LST_05deg.template.gz", "LST_05deg_time.template.gz"], "NectarCam": ["MST_05deg.template.gz", "MST_05deg_time.template.gz"], "FlashCam": ["MST_xm_full.fits"], } # We also need a conversion function from height above ground to # depth of maximum To do this we need the conversion table from CORSIKA ( self.thickness_profile, self.altitude_profile, ) = get_atmosphere_profile_functions("paranal", with_units=False) # Next we need the position, area and amplitude from each pixel in the event # making this a class member makes passing them around much easier self.pixel_x, self.pixel_y = None, None self.image, self.time = None, None self.tel_types, self.tel_id = None, None # We also need telescope positions self.tel_pos_x, self.tel_pos_y = None, None # And the peak of the images self.peak_x, self.peak_y, self.peak_amp = None, None, None self.hillas_parameters, self.ped = None, None self.prediction = dict() self.time_prediction = dict() self.array_direction = None self.array_return = False self.nominal_frame = None # For now these factors are required to fix problems in templates self.template_scale = template_scale self.xmax_offset = xmax_offset self.use_time_gradient = use_time_gradient
[docs] def initialise_templates(self, tel_type): """Check if templates for a given telescope type has been initialised and if not do it and add to the dictionary Parameters ---------- tel_type: dictionary Dictionary of telescope types in event Returns ------- boolean: Confirm initialisation """ for t in tel_type: if tel_type[t] in self.prediction.keys() or tel_type[t] == "DUMMY": continue self.prediction[tel_type[t]] = TemplateNetworkInterpolator( self.root_dir + "/" + self.file_names[tel_type[t]][0] ) if self.use_time_gradient: self.time_prediction[tel_type[t]] = TimeGradientInterpolator( self.root_dir + "/" + self.file_names[tel_type[t]][1] ) return True
[docs] def get_hillas_mean(self): """This is a simple function to find the peak position of each image in an event which will be used later in the Xmax calculation. Peak is found by taking the average position of the n hottest pixels in the image. """ peak_x = np.zeros([len(self.pixel_x)]) # Create blank arrays for peaks # rather than a dict (faster) peak_y = np.zeros(peak_x.shape) peak_amp = np.zeros(peak_x.shape) # Loop over all tels to take weighted average of pixel # positions This loop could maybe be replaced by an array # operation by a numpy wizard # Maybe a vectorize? tel_num = 0 for hillas in self.hillas_parameters: peak_x[tel_num] = # Fill up array peak_y[tel_num] = peak_amp[tel_num] = hillas.intensity tel_num += 1 self.peak_x = peak_x # * unit # Add to class member self.peak_y = peak_y # * unit self.peak_amp = peak_amp
# This function would be useful elsewhere so probably be implemented in a # more general form
[docs] def get_shower_max(self, source_x, source_y, core_x, core_y, zen): """Function to calculate the depth of shower maximum geometrically under the assumption that the shower maximum lies at the brightest point of the camera image. Parameters ---------- source_x: float Event source position in nominal frame source_y: float Event source position in nominal frame core_x: float Event core position in telescope tilted frame core_y: float Event core position in telescope tilted frame zen: float Zenith angle of event Returns ------- float: Depth of maximum of air shower """ # Calculate displacement of image centroid from source position (in # rad) disp = np.sqrt((self.peak_x - source_x) ** 2 + (self.peak_y - source_y) ** 2) # Calculate impact parameter of the shower impact = np.sqrt( (self.tel_pos_x - core_x) ** 2 + (self.tel_pos_y - core_y) ** 2 ) # Distance above telescope is ratio of these two (small angle) height = impact / disp weight = np.power(self.peak_amp, 0.0) # weight average by sqrt amplitude # sqrt may not be the best option... # Take weighted mean of estimates mean_height = np.sum(height * weight) / np.sum(weight) # This value is height above telescope in the tilted system, # we should convert to height above ground mean_height *= np.cos(zen) # Add on the height of the detector above sea level mean_height += 2150 if mean_height > 100000 or np.isnan(mean_height): mean_height = 100000 # Lookup this height in the depth tables, the convert Hmax to Xmax x_max = self.thickness_profile(mean_height) # Convert to slant depth x_max /= np.cos(zen) return x_max + self.xmax_offset
[docs] @staticmethod def rotate_translate(pixel_pos_x, pixel_pos_y, x_trans, y_trans, phi): """ Function to perform rotation and translation of pixel lists Parameters ---------- pixel_pos_x: ndarray Array of pixel x positions pixel_pos_y: ndarray Array of pixel x positions x_trans: float Translation of position in x coordinates y_trans: float Translation of position in y coordinates phi: float Rotation angle of pixels Returns ------- ndarray,ndarray: Transformed pixel x and y coordinates """ cosine_angle = np.cos(phi[..., np.newaxis]) sin_angle = np.sin(phi[..., np.newaxis]) pixel_pos_trans_x = (x_trans - pixel_pos_x) * cosine_angle - ( y_trans - pixel_pos_y ) * sin_angle pixel_pos_trans_y = (pixel_pos_x - x_trans) * sin_angle + ( pixel_pos_y - y_trans ) * cosine_angle return pixel_pos_trans_x, pixel_pos_trans_y
[docs] def image_prediction(self, tel_type, energy, impact, x_max, pix_x, pix_y): """Creates predicted image for the specified pixels, interpolated from the template library. Parameters ---------- tel_type: string Telescope type specifier energy: float Event energy (TeV) impact: float Impact diance of shower (metres) x_max: float Depth of shower maximum (num bins from expectation) pix_x: ndarray X coordinate of pixels pix_y: ndarray Y coordinate of pixels Returns ------- ndarray: predicted amplitude for all pixels """ return self.prediction[tel_type](energy, impact, x_max, pix_x, pix_y)
[docs] def predict_time(self, tel_type, energy, impact, x_max): """Creates predicted image for the specified pixels, interpolated from the template library. Parameters ---------- tel_type: string Telescope type specifier energy: float Event energy (TeV) impact: float Impact diance of shower (metres) x_max: float Depth of shower maximum (num bins from expectation) Returns ------- ndarray: predicted amplitude for all pixels """ return self.time_prediction[tel_type](energy, impact, x_max)
[docs] def get_likelihood( self, source_x, source_y, core_x, core_y, energy, x_max_scale, goodness_of_fit=False, ): """Get the likelihood that the image predicted at the given test position matches the camera image. Parameters ---------- source_x: float Source position of shower in the nominal system (in deg) source_y: float Source position of shower in the nominal system (in deg) core_x: float Core position of shower in tilted telescope system (in m) core_y: float Core position of shower in tilted telescope system (in m) energy: float Shower energy (in TeV) x_max_scale: float Scaling factor applied to geometrically calculated Xmax goodness_of_fit: boolean Determines whether expected likelihood should be subtracted from result Returns ------- float: Likelihood the model represents the camera image at this position """ # First we add units back onto everything. Currently not # handled very well, maybe in future we could just put # everything in the correct units when loading in the class # and ignore them from then on zenith = self.array_direction.zen.to_value(u.rad) # Geometrically calculate the depth of maximum given this test position x_max = self.get_shower_max(source_x, source_y, core_x, core_y, zenith) x_max *= x_max_scale # Calculate expected Xmax given this energy x_max_exp = guess_shower_depth(energy) # / np.cos(20*u.deg) # Convert to binning of Xmax x_max_bin = x_max - x_max_exp # Check for range if x_max_bin > 200: x_max_bin = 200 if x_max_bin < -100: x_max_bin = -100 # Calculate impact distance for all telescopes impact = np.sqrt( (self.tel_pos_x - core_x) ** 2 + (self.tel_pos_y - core_y) ** 2 ) # And the expected rotation angle phi = np.arctan2((self.tel_pos_x - core_x), (self.tel_pos_y - core_y)) * u.rad # Rotate and translate all pixels such that they match the # template orientation pix_y_rot, pix_x_rot = self.rotate_translate( self.pixel_x, self.pixel_y, source_x, source_y, phi ) # In the interpolator class we can gain speed advantages by using masked arrays # so we need to make sure here everything is masked prediction = ma.zeros(self.image.shape) prediction.mask = ma.getmask(self.image) time_gradients = np.zeros((self.image.shape[0], 2)) # Loop over all telescope types and get prediction for tel_type in np.unique(self.tel_types).tolist(): type_mask = self.tel_types == tel_type prediction[type_mask] = self.image_prediction( tel_type, energy * np.ones_like(impact[type_mask]), impact[type_mask], x_max_bin * np.ones_like(impact[type_mask]), -np.rad2deg(pix_x_rot[type_mask]), np.rad2deg(pix_y_rot[type_mask]), ) if self.use_time_gradient: time_gradients[type_mask] = self.predict_time( tel_type, energy * np.ones_like(impact[type_mask]), impact[type_mask], x_max_bin * np.ones_like(impact[type_mask]), ) if self.use_time_gradient: time_mask = np.logical_and(np.invert(ma.getmask(self.image)), self.time > 0) weight = np.sqrt(self.image) * time_mask rv = norm() sx = pix_x_rot * weight sxx = pix_x_rot * pix_x_rot * weight sy = self.time * weight sxy = self.time * pix_x_rot * weight d = weight.sum(axis=1) * sxx.sum(axis=1) - sx.sum(axis=1) * sx.sum(axis=1) time_fit = ( weight.sum(axis=1) * sxy.sum(axis=1) - sx.sum(axis=1) * sy.sum(axis=1) ) / d time_fit /= -1 * (180 / math.pi) chi2 = -2 * np.log( rv.pdf((time_fit - time_gradients.T[0]) / time_gradients.T[1]) ) # Likelihood function will break if we find a NaN or a 0 prediction[np.isnan(prediction)] = 1e-8 prediction[prediction < 1e-8] = 1e-8 prediction *= self.template_scale # Get likelihood that the prediction matched the camera image like = neg_log_likelihood(self.image, prediction, self.spe, self.ped) like[np.isnan(like)] = 1e9 like *= np.invert(ma.getmask(self.image)) like = ma.MaskedArray(like, mask=ma.getmask(self.image)) array_like = like if goodness_of_fit: return np.sum( like - mean_poisson_likelihood_gaussian(prediction, self.spe, self.ped) ) prior_pen = 0 # Add prior penalities if we have them array_like += 1e-8 if "energy" in self.priors: prior_pen += energy_prior(energy, index=-1) if "xmax" in self.priors: prior_pen += xmax_prior(energy, x_max) array_like += prior_pen / float(len(array_like)) if self.array_return: array_like = array_like.ravel() return array_like[np.invert(ma.getmask(array_like))] final_sum = array_like.sum() if self.use_time_gradient: final_sum += chi2.sum() # * np.sum(ma.getmask(self.image)) return final_sum
[docs] def get_likelihood_min(self, x): """Wrapper class around likelihood function for use with scipy minimisers Parameters ---------- x: ndarray Array of minimisation parameters Returns ------- float: Likelihood value of test position """ val = self.get_likelihood(x[0], x[1], x[2], x[3], x[4], x[5]) return val
[docs] def get_likelihood_nlopt(self, x, grad): """Wrapper class around likelihood function for use with scipy minimisers Parameters ---------- x: ndarray Array of minimisation parameters Returns ------- float: Likelihood value of test position """ val = self.get_likelihood(x[0], x[1], x[2], x[3], x[4], x[5]) return val
[docs] def set_event_properties( self, image, time, pixel_x, pixel_y, type_tel, tel_x, tel_y, array_direction, hillas, ): """The setter class is used to set the event properties within this class before minimisation can take place. This simply copies a bunch of useful properties to class members, so that we can use them later without passing all this information around. Parameters ---------- image: dict Amplitude of pixels in camera images time: dict Time information per each pixel in camera images pixel_x: dict X position of pixels in nominal system pixel_y: dict Y position of pixels in nominal system type_tel: dict Type of telescope tel_x: dict X position of telescope in TiltedGroundFrame tel_y: dict Y position of telescope in TiltedGroundFrame array_direction: SkyCoord[AltAz] Array pointing direction in the AltAz Frame hillas: dict dictionary with telescope IDs as key and HillasParametersContainer instances as values Returns ------- None """ # First store these parameters in the class so we can use them # in minimisation For most values this is simply copying self.image = image self.tel_pos_x = np.zeros(len(tel_x)) self.tel_pos_y = np.zeros(len(tel_x)) self.ped = np.zeros(len(tel_x)) self.tel_types, self.tel_id = list(), list() max_pix_x = 0 px, py, pa, pt = list(), list(), list(), list() self.hillas_parameters = list() # So here we must loop over the telescopes for x, i in zip(tel_x, range(len(tel_x))): px.append(pixel_x[x].to(u.rad).value) if len(px[i]) > max_pix_x: max_pix_x = len(px[i]) py.append(pixel_y[x].to(u.rad).value) pa.append(image[x]) pt.append(time[x]) self.tel_pos_x[i] = tel_x[x].to(u.m).value self.tel_pos_y[i] = tel_y[x].to(u.m).value self.ped[i] = self.ped_table[type_tel[x]] self.tel_types.append(type_tel[x]) self.tel_id.append(x) self.hillas_parameters.append(hillas[x]) # Most interesting stuff is now copied to the class, but to remove our requirement # for loops we must copy the pixel positions to an array with the length of the # largest image # First allocate everything shape = (len(tel_x), max_pix_x) self.pixel_x, self.pixel_y = ma.zeros(shape), ma.zeros(shape) self.image, self.time, self.ped = ( ma.zeros(shape), ma.zeros(shape), ma.zeros(shape), ) self.tel_types = np.array(self.tel_types) # Copy everything into our masked arrays for i in range(len(tel_x)): array_len = len(px[i]) self.pixel_x[i][:array_len] = px[i] self.pixel_y[i][:array_len] = py[i] self.image[i][:array_len] = pa[i] self.time[i][:array_len] = pt[i] self.ped[i][:array_len] = self.ped_table[self.tel_types[i]] # Set the image mask mask = self.image == 0.0 self.pixel_x[mask], self.pixel_y[mask] = ma.masked, ma.masked self.image[mask] = ma.masked self.time[mask] = ma.masked self.array_direction = array_direction self.nominal_frame = NominalFrame(origin=self.array_direction) # Finally run some functions to get ready for the event self.get_hillas_mean() self.initialise_templates(type_tel)
[docs] def reset_interpolator(self): """ This function is needed in order to reset some variables in the interpolator at each new event. Without this reset, a new event starts with information from the previous event. """ list(self.prediction.values())[0].reset()
[docs] def predict(self, shower_seed, energy_seed): """Predict method for the ImPACT reconstructor. Used to calculate the reconstructed ImPACT shower geometry and energy. Parameters ---------- shower_seed: ReconstructedShowerContainer Seed shower geometry to be used in the fit energy_seed: ReconstructedEnergyContainer Seed energy to be used in fit Returns ------- ReconstructedShowerContainer, ReconstructedEnergyContainer: """ self.reset_interpolator() horizon_seed = SkyCoord(, alt=shower_seed.alt, frame=AltAz()) nominal_seed = horizon_seed.transform_to(self.nominal_frame) source_x = nominal_seed.fov_lon.to_value(u.rad) source_y = nominal_seed.fov_lat.to_value(u.rad) ground = GroundFrame(x=shower_seed.core_x, y=shower_seed.core_y, z=0 * u.m) tilted = ground.transform_to( TiltedGroundFrame(pointing_direction=self.array_direction) ) tilt_x = tilt_y = zenith = self.array_direction.zen seeds = spread_line_seed( self.hillas_parameters, self.tel_pos_x, self.tel_pos_y, source_x, source_y, tilt_x, tilt_y,, shift_frac=[1], )[0] # Perform maximum likelihood fit fit_params, errors, like = self.minimise( params=seeds[0], step=seeds[1], limits=seeds[2], minimiser_name=self.minimiser_name, ) # Create a container class for reconstructed shower # Convert the best fits direction and core to Horizon and ground systems and # copy to the shower container nominal = SkyCoord( fov_lon=fit_params[0] * u.rad, fov_lat=fit_params[1] * u.rad, frame=self.nominal_frame, ) horizon = nominal.transform_to(AltAz()) tilted = TiltedGroundFrame( x=fit_params[2] * u.m, y=fit_params[3] * u.m, z=0 * u.m, pointing_direction=self.array_direction, ) ground = project_to_ground(tilted) h_max = ( fit_params[5] * np.cos(zenith) * self.get_shower_max( fit_params[0], fit_params[1], fit_params[2], fit_params[3],, ) ) shower_result = ReconstructedGeometryContainer( alt=horizon.alt,, core_x=ground.x, core_y=ground.y, is_valid=True, alt_uncert=np.nan, az_uncert=np.nan, core_uncert=np.nan, # Copy reconstructed hmax h_max=h_max, h_max_uncert=errors[5] * h_max, goodness_of_fit=like, prefix=self.__class__.__name__, ) # Create a container class for reconstructed energy energy_result = ReconstructedEnergyContainer( prefix=self.__class__.__name__, energy=fit_params[4] * u.TeV, energy_uncert=errors[4] * u.TeV, is_valid=True, ) return shower_result, energy_result
[docs] def minimise(self, params, step, limits, minimiser_name="minuit", max_calls=0): """ Parameters ---------- params: ndarray Seed parameters for fit step: ndarray Initial step size in the fit limits: ndarray Fit bounds minimiser_name: str Name of minimisation method max_calls: int Maximum number of calls to minimiser Returns ------- tuple: best fit parameters and errors """ limits = np.asarray(limits) if minimiser_name == "minuit": self.min = Minuit( self.get_likelihood, print_level=1, source_x=params[0], error_source_x=step[0], limit_source_x=limits[0], fix_source_x=False, source_y=params[1], error_source_y=step[1], limit_source_y=limits[1], fix_source_y=False, core_x=params[2], error_core_x=step[2], limit_core_x=limits[2], fix_core_x=False, core_y=params[3], error_core_y=step[3], limit_core_y=limits[3], fix_core_y=False, energy=params[4], error_energy=step[4], limit_energy=limits[4], fix_energy=False, x_max_scale=params[5], error_x_max_scale=step[5], limit_x_max_scale=limits[5], fix_x_max_scale=False, goodness_of_fit=False, fix_goodness_of_fit=True, errordef=1, ) self.min.tol *= 1000 self.min.set_strategy(1) self.min.migrad() fit_params = self.min.values errors = self.min.errors return ( ( fit_params["source_x"], fit_params["source_y"], fit_params["core_x"], fit_params["core_y"], fit_params["energy"], fit_params["x_max_scale"], ), ( errors["source_x"], errors["source_y"], errors["core_x"], errors["core_x"], errors["energy"], errors["x_max_scale"], ), self.min.fval, ) elif "nlopt" in minimiser_name: import nlopt opt = nlopt.opt(nlopt.LN_BOBYQA, 6) opt.set_min_objective(self.get_likelihood_nlopt) opt.set_initial_step(step) opt.set_lower_bounds(np.asarray(limits).T[0]) opt.set_upper_bounds(np.asarray(limits).T[1]) opt.set_xtol_rel(1e-3) if max_calls: opt.set_maxeval(max_calls) x = opt.optimize(np.asarray(params)) return x, (0, 0, 0, 0, 0, 0), self.get_likelihood_min(x) elif minimiser_name in ("lm", "trf", "dogleg"): self.array_return = True min = least_squares( self.get_likelihood_min, params, method=minimiser_name, x_scale=step, xtol=1e-10, ftol=1e-10, ) return min.x, (0, 0, 0, 0, 0, 0), self.get_likelihood_min(min.x) else: min = minimize( self.get_likelihood_min, np.array(params), method=minimiser_name, bounds=limits, options={"disp": False}, tol=1e-5, ) return np.array(min.x), (0, 0, 0, 0, 0, 0), self.get_likelihood_min(min.x)
def spread_line_seed( hillas, tel_x, tel_y, source_x, source_y, tilt_x, tilt_y, energy, shift_frac=[2, 1.5, 1, 0.5, 0, -0.5, -1, -1.5], ): """ Parameters ---------- hillas: list Hillas parameters in event tel_x: list telescope X positions in tilted system tel_y: list telescope Y positions in tilted system source_x: float Source X position in nominal system (radians) source_y:float Source Y position in nominal system (radians) tilt_x: float Core X position in tilited system (radians) tilt_y: float Core Y position in tilited system (radians) energy: float Energy in TeV shift_frac: list Fractional values to shist source and core positions Returns ------- list of seed positions to try """ centre_x, centre_y, amp = list(), list(), list() for tel_hillas in hillas: centre_x.append( centre_y.append( amp.append(tel_hillas.intensity) centre_x = np.average(centre_x, weights=amp) centre_y = np.average(centre_y, weights=amp) centre_tel_x = np.average(tel_x, weights=amp) centre_tel_y = np.average(tel_y, weights=amp) diff_x = source_x - centre_x diff_y = source_y - centre_y diff_tel_x = tilt_x - centre_tel_x diff_tel_y = tilt_y - centre_tel_y seed_list = list() for shift in shift_frac: seed_list.append( create_seed( centre_x + (diff_x * shift), centre_y + (diff_y * shift), centre_tel_x + (diff_tel_x * shift), centre_tel_y + (diff_tel_y * shift), energy, ) ) return seed_list def create_seed(source_x, source_y, tilt_x, tilt_y, energy): """ Function for creating seed, step and limits for a given position Parameters ---------- source_x: float Source X position in nominal system (radians) source_y:float Source Y position in nominal system (radians) tilt_x: float Core X position in tilited system (radians) tilt_y: float Core Y position in tilited system (radians) energy: float Energy in TeV Returns ------- tuple of seed, steps size and fit limits """ lower_en_limit = energy * 0.5 en_seed = energy # If our energy estimate falls outside of the range of our templates set it to # the edge if lower_en_limit < 0.01: lower_en_limit = 0.01 en_seed = 0.01 # Take the seed from Hillas-based reconstruction seed = (source_x, source_y, tilt_x, tilt_y, en_seed, 1) # Take a reasonable first guess at step size step = [0.04 / 57.3, 0.04 / 57.3, 5, 5, en_seed * 0.1, 0.05] # And some sensible limits of the fit range limits = [ [source_x - 0.1, source_x + 0.1], [source_y - 0.1, source_y + 0.1], [tilt_x - 100, tilt_x + 100], [tilt_y - 100, tilt_y + 100], [lower_en_limit, en_seed * 2], [0.5, 2], ] return seed, step, limits