"""
@author: BACON Philippe, APC/IN2P3/CNRS
"""
from itertools import product
import logging
import warnings
import numpy as np
from numba import njit
from ecpi.common.instru.model_geom import InstruECLAIRs
from ecpi.common.num.integration_tools import solidangle_rectangle_ana_numba
from ecpi.common.num.array_square_cell import ArraySquareCell
EPS_LIMIT = 2.3e-8
# warnings.simplefilter('ignore', category=NumbaPendingDeprecationWarning)
[docs]def plot(arr_pos):
"""..warning: for diagnostic purposes
"""
import matplotlib.pyplot as plt
plt.figure()
for x,y in arr_pos:
plt.scatter(x,y)
plt.show()
[docs]class SolidAngleComputation(object):
"""Solid angle computation for the ECLAIRs instrument.
"""
def __init__(self, instru_ecl : InstruECLAIRs, use_numba=False):
"""**constructor**
..note :: numba mode only concerns the calculator method
for the moment
:param instru_ecl: description of the mask of ECLAIRs
:type instru_ecl: InstruECLAIRs object
:param use_numba: whether or not to optimize the detector-mask loops
:type use_numba: bool. Default is False.
"""
assert isinstance(instru_ecl, InstruECLAIRs)
self.l_submasks = [instru_ecl._l_mask[m_idx] for m_idx in range(4)]
self.pixpitch = instru_ecl.pitch_pix_size
self.detpix = instru_ecl.p_nb_pix
self.detmas = instru_ecl.dist_mask_detect
self.massize = instru_ecl.mask_total_size
self.masnerv = instru_ecl.half_cross * 2
self.maspsize = instru_ecl.mask_pix_size
self.maspix = instru_ecl.mask_nb_pix
self.detsize = self.pixpitch * self.detpix
self.detpsize = instru_ecl.det_pix_size
self.instru_ecl = instru_ecl
# declare logger
self.logger = logging.getLogger(__name__)
self.logger.info(f'creating an instance of {__class__}')
# decide which function to call when numba is (in)active
if use_numba:
self.logger.info('running optimized version')
self.loops_computation_func = _run_over_detector_mask_pixels_numba
else:
self.logger.info('running default version: this may take a while...')
self.loops_computation_func = _run_over_detector_mask_pixels
[docs] def set_mask_configuration(self, sub_mask_arr, fully_opened_mask):
"""Sets the mask configuration: whether or not to compute
solid angle with the real mask geometry or with open submasks.
:param sub_mask_arr: array specifying submask geometry
:type sub_mask_arr: 2D array(int). Filled with 0s and 1s.
:param fully_opened_mask: whether or not to consider the opened mask.
The ECLAIRs maks cross remains !
:type fully_opened_mask: bool
"""
if fully_opened_mask:
self.logger.info(f'Using mask configuration: fully opened submasks')
sub_mask_holes_row, sub_mask_holes_col = np.array(np.where(sub_mask_arr != 2))
else:
self.logger.info(f'Using mask configuration: real mask')
sub_mask_holes_row, sub_mask_holes_col = np.array(np.where(sub_mask_arr == 1))
return sub_mask_holes_row, sub_mask_holes_col
[docs] def run_calculator(self, fully_opened_mask=False, debug=False):
"""Compute solid angle map without the Earth in the FOV.
documentation:
* simulation_correction_cxb_earth_15_05_2020_v2_alg.pdf (version 2) p.10->11
* simulation_correction_pipeline_cxb_earth_15_05_2020_v2_pb.pdf (version 2) p.4
:param fully_opened_mask: whether or not to consider the opened mask.
The ECLAIRs mask cross remains ! Default is False.
:type fully_opened_mask: bool
:return: solid angle map
:rtype: 2D array(float) with shape 80x80
..note: the computation includes the cos theta effect. It follows the sum of the
solid angle is 2199.0922 instead of 2452.0531 (without cos theta case).
"""
solid_angle_with_costheta = np.zeros((self.detpix, self.detpix))
self.logger.info('starting computation of solid angle with Earth')
for sub_mask in self.l_submasks:
sub_mask_arr = sub_mask.ar_val
sub_mask_holes_row, sub_mask_holes_col = self.set_mask_configuration(
sub_mask_arr, fully_opened_mask
)
sub_mask_pos = np.array([sub_mask.idx2pos(np.array([k, l])) \
for (k, l) in zip(sub_mask_holes_row, sub_mask_holes_col)])
self.logger.info(f'processing submask {sub_mask.name}...')
l_detpix = list(product(range(self.detpix), repeat=2))
sa_array = self.loops_computation_func(l_detpix, sub_mask_pos, self.detpix,
self.pixpitch, self.detmas, self.maspsize, self.detpsize)
solid_angle_with_costheta += sa_array
if debug:
print(sub_mask_pos)
plot(sub_mask_pos)
return np.flipud(solid_angle_with_costheta)
def _run_over_detector_mask_pixels(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize
): # pragma: no cover
"""Loop over detector-mask pixels so as to effectively compute the solid angle.
Pure python implementation.
"""
ox_unit = np.array([1, 0, 0]) # define (Ox) axis normal to detector plane.
sa_array = np.zeros((detpix, detpix))
for (i,j) in l_detpix:
# determine physical center of detector pixel d(i,j)
y_det_c = (i+0.5) * pixpitch - 0.5*(detpix*pixpitch)
z_det_c = (j+0.5) * pixpitch - 0.5*(detpix*pixpitch)
for (y_msk_c, z_msk_c) in l_submaskpos: # loop over mask holes
# compute solid angle subtended by a pixel m(k,l) of mask as seen by
# the center of detector pixel d(i,j)
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
xp = detmas #+ 0.03
# compute solid angle without Earth.
sa_mask_pix = solidangle_rectangle_ana_numba(maspsize, maspsize, yp, zp, xp)
# compute direction from detector pixel to mask pixel.
v_unit = np.array([xp, yp, zp]) / np.linalg.norm([xp, yp, zp])
# compute angles.
cos_theta = np.abs(np.dot(v_unit, ox_unit))
sa_array[i, j] += sa_mask_pix * cos_theta
return sa_array
@njit
def _run_over_detector_mask_pixels_numba(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize
): # pragma: no cover
"""Loop over detector-mask pixels so as to effectively compute the solid angle.
Python/C implementation (numba).
..warning: seems upgrading to numba==0.48.XXX could cure use of linalg.norm and cos. TBC
"""
sa_array = np.zeros((detpix, detpix))
for (i,j) in l_detpix:
y_det_c = (i+0.5) * pixpitch - 0.5*(detpix*pixpitch)
z_det_c = (j+0.5) * pixpitch - 0.5*(detpix*pixpitch)
for k in range(len(l_submaskpos)):
y_msk_c = l_submaskpos[k][0]
z_msk_c = l_submaskpos[k][1]
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
xp = detmas
sa_mask_pix = solidangle_rectangle_ana_numba(maspsize, maspsize, yp, zp, xp)
v_unit = np.array([xp, yp, zp]) / (xp**2 + yp**2 + zp**2)**0.5
cos_theta = v_unit[0]
sa_array[i, j] += sa_mask_pix * cos_theta
return sa_array
[docs]class SolidAngleComputationWithEarth(SolidAngleComputation):
"""Solid angle computation for the ECLAIRs instrument
in the presence of Earth in the FOV.
"""
def __init__(self, instru_ecl : InstruECLAIRs, use_numba=False):
"""**constructor**
..note :: numba mode only concerns the calculator method
for the moment
:param instru_ecl: description of the mask of ECLAIRs
:type instru_ecl: InstruECLAIRs object
:param use_numba: whether or not to optimize the detector-mask loops
:type use_numba: bool. Default is False.
"""
super().__init__(instru_ecl)
# decide which function to call when numba is (in)active
if use_numba:
print('running optimized version')
self.logger.info('running optimized version')
self.loops_computation_func = _run_over_detector_mask_pixels_with_earth_numba
else:
print('running default version: this may take a while...')
self.logger.info('running default version: this may take a while...')
self.loops_computation_func = _run_over_detector_mask_pixels_with_earth
self.use_numba = use_numba
[docs] def run_calculator(self, limb_angle, u_earth, fully_opened_mask=False, debug=False):
"""Compute solid angle map with the Earth in the FOV using
the mean CXB shape calculation.
..note: when angle(vec direction, vec Earth) very close to limb angle
may lead to numerical errors in comparsion. introduced a bias in
solid angle computation. Check has been added that raises a warning.
:param limb_angle: limb angle (cf. doc) in degree
:type limb_angle: float
:param u_earth: unit vector pointing towards Earth center
:type u_earth: [float, float, float]
:param fully_opened_mask: whether or not to consider the opened mask.
The ECLAIRs maks cross remains !
:type fully_opened_mask: bool
:return: solid angle map
:rtype: 2D array(float) with shape 80x80
"""
print("DEBUG: run SolidAngleComputationWithEarth")
limb = np.deg2rad(limb_angle)
solid_angle_with_costheta = np.zeros((self.detpix, self.detpix))
self.logger.info('starting computation of solid angle with Earth')
for sub_mask in self.l_submasks:
sub_mask_arr = sub_mask.ar_val
sub_mask_holes_row, sub_mask_holes_col = self.set_mask_configuration(
sub_mask_arr, fully_opened_mask
)
sub_mask_pos = np.array([sub_mask.idx2pos(np.array([k, l])) \
for (k, l) in zip(sub_mask_holes_row, sub_mask_holes_col)])
logger = self.logger
self.logger.info(f'processing submask {sub_mask.name}...')
l_detpix = list(product(range(self.detpix), repeat=2))
if self.use_numba:
sa_array = self.loops_computation_func(
l_detpix, sub_mask_pos, self.detpix, self.pixpitch,
self.detmas, self.maspsize, self.detpsize, limb, u_earth
)
else:
sa_array = self.loops_computation_func(
l_detpix, sub_mask_pos, self.detpix, self.pixpitch,
self.detmas, self.maspsize, self.detpsize, limb, u_earth, logger
)
solid_angle_with_costheta += sa_array
if debug:
print(sub_mask_pos)
plot(sub_mask_pos)
return np.flipud(solid_angle_with_costheta)
def _run_over_detector_mask_pixels_with_earth(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize,
limb, u_earth, logger
): # pragma: no cover
"""Loop over detector-mask pixels so as to effectively compute the solid angle
in presence of Earth in FOV.
Pure python implementation.
"""
ox_unit = np.array([1, 0, 0]) # define (Ox) axis normal to detector plane.
sa_array = np.zeros((detpix, detpix))
for (i,j) in l_detpix:
# determine physical center of detector pixel d(i,j)
y_det_c = (i+0.5) * pixpitch - 0.5*(detpix*pixpitch)
z_det_c = (j+0.5) * pixpitch - 0.5*(detpix*pixpitch)
for (y_msk_c, z_msk_c) in l_submaskpos: # loop over mask holes
# compute solid angle subtended by a pixel m(k,l) of mask as seen by
# the center of detector pixel d(i,j)
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
xp = detmas
# compute solid angle without Earth.
sa_mask_pix = solidangle_rectangle_ana_numba(maspsize, maspsize, yp, zp, xp)
sa_array[i, j] += sa_mask_pix
# compute direction from detector pixel to mask pixel.
v_unit = np.array([xp, yp, zp]) / np.linalg.norm([xp, yp, zp])
# compute angles.
cos_alpha = np.dot(v_unit, u_earth)
cos_theta = np.abs(np.dot(v_unit, ox_unit))
alpha = np.arccos(cos_alpha)
if np.abs(limb - alpha) < EPS_LIMIT:
logger.warning(f'Detected pixel with alpha close to limb at pixel ({i}, {j})')
logger.warning(f'limb angle={limb} rad | alpha angle={alpha} rad')
if not alpha < limb:
sa_array[i, j] += sa_mask_pix * cos_theta
return sa_array
@njit
def _run_over_detector_mask_pixels_with_earth_numba(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize,
limb, u_earth
): # pragma: no cover
"""Loop over detector-mask pixels so as to effectively compute the solid angle
in presence of Earth in FOV.
Python/C implementation (numba).
..warning: seems upgrading to numba==0.48.XXX could cure use of linalg.norm and cos. TBC
"""
#TODO PB: cannot communicate with logger (numba).
sa_array = np.zeros((detpix, detpix))
for (i,j) in l_detpix:
y_det_c = (i+0.5) * pixpitch - 0.5*(detpix*pixpitch)
z_det_c = (j+0.5) * pixpitch - 0.5*(detpix*pixpitch)
for k in range(len(l_submaskpos)):
y_msk_c = l_submaskpos[k][0]
z_msk_c = l_submaskpos[k][1]
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
xp = detmas
sa_mask_pix = solidangle_rectangle_ana_numba(maspsize, maspsize, yp, zp, xp)
v_unit = np.array([xp, yp, zp]) / (xp**2 + yp**2 + zp**2)**0.5
cos_alpha = v_unit[0]*u_earth[0] + v_unit[1]*u_earth[1] + v_unit[2]*u_earth[2] #np.dot(v_unit, u_earth)
cos_theta = v_unit[0]
alpha = np.arccos(cos_alpha)
if not alpha < limb:
sa_array[i, j] += sa_mask_pix * cos_theta
return sa_array
[docs]class SolidAngleWithEarthOrNot(SolidAngleComputation):
"""Solid angle computation for the ECLAIRs instrument
in the presence of Earth in the FOV.
"""
def __init__(self, instru_ecl : InstruECLAIRs, fully_opened_mask=False):
"""**constructor**
..note :: numba mode only concerns the calculator method
for the moment
:param instru_ecl: description of the mask of ECLAIRs
:type instru_ecl: InstruECLAIRs object
:param use_numba: whether or not to optimize the detector-mask loops
:type use_numba: bool. Default is False.
"""
super().__init__(instru_ecl)
# liste de couple d'indice des pixels detecteur
self.l_detpix = list(product(range(self.detpix), repeat=2))
# list des solid angle by hole mask
self.l_sabhm = []
self.l_sub_mask_pos = []
self.do_solid_angle_by_hole(fully_opened_mask=fully_opened_mask)
[docs] def do_solid_angle_by_hole(self, fully_opened_mask=False):
self.logger.info('starting computation of solid angle for each pixel det-mask')
for sub_mask in self.l_submasks:
sub_mask_arr = sub_mask.ar_val
sub_mask_holes_row, sub_mask_holes_col = self.set_mask_configuration(
sub_mask_arr, fully_opened_mask
)
sub_mask_pos = np.array([sub_mask.idx2pos(np.array([k, l])) \
for (k, l) in zip(sub_mask_holes_row, sub_mask_holes_col)])
self.l_sub_mask_pos.append(sub_mask_pos)
sabhm = solid_angle_for_mask_hole_numba(self.l_detpix,
sub_mask_pos,
self.detpix,
self.pixpitch,
self.detmas,
self.maspsize,
self.detpsize)
self.l_sabhm.append(sabhm)
[docs] def run_calculator(self, limb_angle, u_earth, fully_opened_mask=False, debug=False):
"""Compute solid angle map with the Earth in the FOV using
the mean CXB shape calculation.
..note: when angle(vec direction, vec Earth) very close to limb angle
may lead to numerical errors in comparsion. introduced a bias in
solid angle computation. Check has been added that raises a warning.
:param limb_angle: limb angle (cf. doc) in degree
:type limb_angle: float
:param u_earth: unit vector pointing towards Earth center
:type u_earth: [float, float, float]
:param fully_opened_mask: whether or not to consider the opened mask.
The ECLAIRs maks cross remains !
:type fully_opened_mask: bool
:return: solid angle map
:rtype: 2D array(float) with shape 80x80
"""
limb = np.deg2rad(limb_angle)
self.logger.info('starting computation of solid angle with Earth')
sa_array = solid_angle_with_earth_numba(self.l_detpix,
self.l_sub_mask_pos,
self.detpix,
self.pixpitch,
self.detmas,
self.maspsize,
self.detpsize,
limb, u_earth,
self.l_sabhm)
return np.flipud(sa_array)
@njit
def solid_angle_for_mask_hole_numba(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize
): # pragma: no cover
"""Loop over detector-mask pixels so as to effectively compute the solid angle
in presence of Earth in FOV.
Python/C implementation (numba).
..warning: seems upgrading to numba==0.48.XXX could cure use of linalg.norm and cos. TBC
"""
#TODO PB: cannot communicate with logger (numba).
sa_array = np.zeros((detpix, detpix, len(l_submaskpos)),dtype=np.float32)
for (i,j) in l_detpix:
y_det_c = (i+0.5) * pixpitch - 0.5*(detpix*pixpitch)
z_det_c = (j+0.5) * pixpitch - 0.5*(detpix*pixpitch)
for k in range(len(l_submaskpos)):
y_msk_c = l_submaskpos[k][0]
z_msk_c = l_submaskpos[k][1]
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
xp = detmas
sa_mask_pix = solidangle_rectangle_ana_numba(maspsize, maspsize, yp, zp, xp)
v_unit = np.array([xp, yp, zp]) / np.sqrt(xp**2 + yp**2 + zp**2)
cos_theta = v_unit[0]
sa_array[i, j, k] = sa_mask_pix * cos_theta
return sa_array
@njit(cache=True)
def solid_angle_with_earth_numba(
l_detpix, l_submaskpos, detpix,
pixpitch, detmas, maspsize, detpsize,
limb, u_earth, l_sabhm
): # pragma: no cover
"""Loop over detector and all submask.
..note: by Colley JM, optimization of _run_over_detector_mask_pixels_with_earth_numba,
* used precomputing of solid angle of each pixel
det-mask, it's a constant of the system
* include all submask
* earth presence: test on cos angle not angle directly (remove arccos)
* => ratio time ~ 10
"""
sa_array = np.zeros((detpix, detpix))
cos_limb = np.cos(limb)
xp = detmas
xp2 = xp**2
for sabhm, submaskpos in zip(l_sabhm, l_submaskpos):
for (i,j) in l_detpix:
y_det_c = (i+0.5) * pixpitch - 0.5*(detpix*pixpitch)
z_det_c = (j+0.5) * pixpitch - 0.5*(detpix*pixpitch)
for k in range(len(submaskpos)):
y_msk_c = submaskpos[k][0]
z_msk_c = submaskpos[k][1]
yp = y_msk_c - y_det_c
zp = z_msk_c - z_det_c
norm_p = np.sqrt(xp2 + yp**2 + zp**2)
# scalar product : cos(angle)*norm_p
cos_alpha_m = xp*u_earth[0] + yp*u_earth[1] + zp*u_earth[2]
if cos_alpha_m < cos_limb*norm_p:
sa_array[i, j] += sabhm[i,j,k]
return sa_array