Skip to content

Commit

Permalink
Merge pull request #606 from VChristiaens/master
Browse files Browse the repository at this point in the history 
Updated hciplot requirement to solve broken dependencies
  • Loading branch information
@VChristiaens
VChristiaens authored Jun 28, 2023
2 parents 781d7bf + 20b1938 commit 6cbea8d
Show file tree
Hide file tree
Showing 2 changed files with 44 additions and 28 deletions.
2 changes: 1 addition & 1 deletion requirements.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ psutil
pyprind
munch
nbsphinx
hciplot>=0.2.1
hciplot>=0.2.4
typing
dataclass_builder
urllib3 < 2.0
70 changes: 43 additions & 27 deletions vip_hci/invprob/fmmf.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,7 +90,8 @@ class FMMFParams:
model: str = "KLIP"
var: Enum = VarEstim.FR
param: dict = field(
default_factory=lambda: {"ncomp": 20, "tolerance": 5e-3, "delta_rot": 0.5}
default_factory=lambda: {"ncomp": 20,
"tolerance": 5e-3, "delta_rot": 0.5}
)
crop: int = 5
imlib: Enum = Imlib.VIPFFT
Expand Down Expand Up @@ -214,7 +215,8 @@ def fmmf(*all_args, **all_kwargs: dict):
if algo_params.nproc is None:
algo_params.nproc = cpu_count() // 2

add_params = {"ann_center": iterable(range(algo_params.min_r, algo_params.max_r))}
add_params = {"ann_center": iterable(
range(algo_params.min_r, algo_params.max_r))}

func_params = setup_parameters(
params_obj=algo_params,
Expand All @@ -230,8 +232,10 @@ def fmmf(*all_args, **all_kwargs: dict):
*func_params,
)

flux_matrix = np.zeros((algo_params.cube.shape[1], algo_params.cube.shape[2]))
snr_matrix = np.zeros((algo_params.cube.shape[1], algo_params.cube.shape[2]))
flux_matrix = np.zeros(
(algo_params.cube.shape[1], algo_params.cube.shape[2]))
snr_matrix = np.zeros(
(algo_params.cube.shape[1], algo_params.cube.shape[2]))

for res_temp in res_full:
indices = get_annulus_segments(algo_params.cube[0], res_temp[2], 1)
Expand Down Expand Up @@ -282,7 +286,8 @@ def _snr_contrast_esti(
# Computation of the reference PSF, and the matrices
# required for the computation of the PSF forward models

pa_threshold = np.rad2deg(2 * np.arctan(delta_rot * fwhm / (2 * (ann_center))))
pa_threshold = np.rad2deg(
2 * np.arctan(delta_rot * fwhm / (2 * (ann_center))))
mid_range = np.abs(np.amax(angle_list) - np.amin(angle_list)) / 2
if pa_threshold >= mid_range - mid_range * 0.1:
pa_threshold = float(mid_range - mid_range * 0.1)
Expand Down Expand Up @@ -395,7 +400,8 @@ def _snr_contrast_esti(
)

psf_map[b, indices[0][0], indices[0][1]] = psf_map_temp
psf_map[b, indices[0][0], indices[0][1]] -= np.mean(psf_map_temp)
psf_map[b, indices[0][0], indices[0]
[1]] -= np.mean(psf_map_temp)

psf_map_der = cube_derotate(
psf_map, angle_list, imlib=imlib, interpolation=interpolation
Expand All @@ -414,7 +420,8 @@ def _snr_contrast_esti(

cube_res_fc = np.zeros_like(model_matrix)

matrix_res_fc = np.zeros((values_fc.shape[0], indices[0][0].shape[0]))
matrix_res_fc = np.zeros(
(values_fc.shape[0], indices[0][0].shape[0]))

for e in range(values_fc.shape[0]):
recon_fc = np.dot(coef_list[e], values_fc[ind_ref_list[e]])
Expand Down Expand Up @@ -455,14 +462,16 @@ def _snr_contrast_esti(
psfm = frame_crop(
psfm_temp[j],
crop,
cenxy=[int(psfm_temp.shape[-1] / 2), int(psfm_temp.shape[-1] / 2)],
cenxy=[int(psfm_temp.shape[-1] / 2),
int(psfm_temp.shape[-1] / 2)],
verbose=False,
)

num.append(
np.multiply(
frame_crop(
mcube[j], crop, cenxy=[poscentx, poscenty], verbose=False
mcube[j], crop, cenxy=[
poscentx, poscenty], verbose=False
),
psfm,
).sum()
Expand All @@ -489,7 +498,8 @@ def _var_esti(mcube, angle_list, var, crop, ann_center):
if var == "FR":
var_f = np.zeros(n)

indices = get_annulus_segments(mcube[0], ann_center - int(crop / 2), crop, 1)
indices = get_annulus_segments(
mcube[0], ann_center - int(crop / 2), crop, 1)

poscentx = indices[0][1]
poscenty = indices[0][0]
Expand All @@ -504,7 +514,8 @@ def _var_esti(mcube, angle_list, var, crop, ann_center):

var_f = np.zeros((len(indicesy), n))

indices = get_annulus_segments(mcube[0], ann_center - int(crop / 2), crop, 1)
indices = get_annulus_segments(
mcube[0], ann_center - int(crop / 2), crop, 1)

for a in range(len(indicesy)):
indc = disk((indicesy[a], indicesx[a]), 3)
Expand Down Expand Up @@ -532,23 +543,28 @@ def _var_esti(mcube, angle_list, var, crop, ann_center):

for a in range(0, len(indicesy)):
radist = np.sqrt(
(indicesx[a] - int(x / 2)) ** 2 + (indicesy[a] - int(y / 2)) ** 2
(indicesx[a] - int(x / 2)) ** 2 +
(indicesy[a] - int(y / 2)) ** 2
)

if (indicesy[a] - int(y / 2)) >= 0:
ang_s = np.arccos((indicesx[a] - int(x / 2)) / radist) / np.pi * 180
ang_s = np.arccos(
(indicesx[a] - int(x / 2)) / radist) / np.pi * 180
else:
ang_s = (
360 - np.arccos((indicesx[a] - int(x / 2)) / radist) / np.pi * 180
360 -
np.arccos((indicesx[a] - int(x / 2)) / radist) / np.pi * 180
)

for b in range(n):
twopi = 2 * np.pi
sigposy = int(
y / 2 + np.sin((ang_s - angle_list[b]) / 360 * twopi) * radist
y / 2 +
np.sin((ang_s - angle_list[b]) / 360 * twopi) * radist
)
sigposx = int(
x / 2 + np.cos((ang_s - angle_list[b]) / 360 * twopi) * radist
x / 2 +
np.cos((ang_s - angle_list[b]) / 360 * twopi) * radist
)

y0 = int(sigposy - int(crop / 2))
Expand Down Expand Up @@ -673,9 +689,8 @@ def _perturb(
return model_sci[None, :] - klipped_oversub - klipped_selfsub


def KLIP_patch(
frame, matrix, numbasis, angle_list, fwhm, pa_threshold, ann_center, nframes=None
):
def KLIP_patch(frame, matrix, numbasis, angle_list, fwhm, pa_threshold,
ann_center, nframes=None):
"""
Function allowing the computation of the reference PSF via KLIP for a
given sub-region of the original ADI sequence. Code inspired by the
Expand Down Expand Up @@ -720,8 +735,8 @@ def KLIP_patch(

# Computation of the eigenvectors/values of the covariance matrix
evals, evecs = la.eigh(covar_psfs)
evals = np.copy(evals[int(tot_basis - max_basis) : int(tot_basis)])
evecs = np.copy(evecs[:, int(tot_basis - max_basis) : int(tot_basis)])
evals = np.copy(evals[int(tot_basis - max_basis): int(tot_basis)])
evecs = np.copy(evecs[:, int(tot_basis - max_basis): int(tot_basis)])
evals = np.copy(evals[::-1])
evecs = np.copy(evecs[:, ::-1])

Expand All @@ -735,7 +750,7 @@ def KLIP_patch(
sci_rows = np.reshape(sci_mean_sub, (1, N_pix))

inner_products = np.dot(sci_rows, KL_basis.T)
inner_products[0, int(max_basis) : :] = 0
inner_products[0, int(max_basis)::] = 0

# Projection of the science image on the selected prinicpal component
# to generate the speckle field model
Expand All @@ -757,9 +772,8 @@ def KLIP_patch(
)


def LOCI_FM(
cube, psf, ann_center, angle_list, asize, fwhm, Tol, delta_rot, pa_threshold
):
def LOCI_FM(cube, psf, ann_center, angle_list, asize, fwhm, Tol, delta_rot,
pa_threshold):
"""
Computation of the optimal factors weigthing the linear combination of
reference frames used to obtain the modeled speckle field for each frame
Expand Down Expand Up @@ -810,7 +824,8 @@ def LOCI_FM(
return cube_res, ind_ref_list, coef_list


def _leastsq_patch_fm(ayxyx, angle_list, fwhm, cube, dist_threshold, tol, psf=None):
def _leastsq_patch_fm(ayxyx, angle_list, fwhm, cube, dist_threshold, tol,
psf=None):
"""
Function allowing th estimation of the optimal factors for the modeled
speckle field estimation via the LOCI framework. The code has been
Expand Down Expand Up @@ -839,7 +854,8 @@ def _leastsq_patch_fm(ayxyx, angle_list, fwhm, cube, dist_threshold, tol, psf=No
n_frames = cube.shape[0]

for i in range(n_frames):
ind_fr_i = _find_indices_adi(angle_list, i, pa_threshold, truncate=False)
ind_fr_i = _find_indices_adi(
angle_list, i, pa_threshold, truncate=False)
if len(ind_fr_i) > 0:
A = values_opt[ind_fr_i]
b = values_opt[i]
Expand Down

0 comments on commit 6cbea8d

Please sign in to comment.