From d62724076283449e5b63b3b01a96e94ab9971c94 Mon Sep 17 00:00:00 2001 From: Alan Loh <24960512+AlanLoh@users.noreply.github.com> Date: Wed, 11 Dec 2024 11:53:10 +0100 Subject: [PATCH] New SolarProximityCnst constraint, extend obs bug correction --- nenupy/__init__.py | 2 +- nenupy/schedule/__init__.py | 1 + nenupy/schedule/constraints.py | 98 +++++++++++++++++++++++++++++++--- nenupy/schedule/geneticalgo.py | 2 +- nenupy/schedule/schedule.py | 94 +++++++++++++++++++++++--------- 5 files changed, 161 insertions(+), 36 deletions(-) diff --git a/nenupy/__init__.py b/nenupy/__init__.py index c44afa3..8ff993a 100644 --- a/nenupy/__init__.py +++ b/nenupy/__init__.py @@ -5,7 +5,7 @@ __copyright__ = "Copyright 2023, nenupy" __credits__ = ["Alan Loh"] __license__ = "MIT" -__version__ = "2.7.13" +__version__ = "2.7.14" __maintainer__ = "Alan Loh" __email__ = "alan.loh@obspm.fr" diff --git a/nenupy/schedule/__init__.py b/nenupy/schedule/__init__.py index 8b06a27..1a4d423 100644 --- a/nenupy/schedule/__init__.py +++ b/nenupy/schedule/__init__.py @@ -14,6 +14,7 @@ LocalTimeCnst, TimeRangeCnst, NightTimeCnst, + SolarProximityCnst, Constraints ) from .obsblocks import Block, ObsBlock, ReservedBlock diff --git a/nenupy/schedule/constraints.py b/nenupy/schedule/constraints.py index 55aac48..ba66eda 100644 --- a/nenupy/schedule/constraints.py +++ b/nenupy/schedule/constraints.py @@ -176,7 +176,7 @@ from copy import copy from typing import List -from nenupy.schedule.targets import _Target +from nenupy.schedule.targets import _Target, SSTarget import logging log = logging.getLogger(__name__) @@ -1024,6 +1024,74 @@ def _evaluate(self, sun_elevation): # ============================================================= # # ============================================================= # +# ============================================================= # +# -------------------- SolarProximityCnst --------------------- # +# ============================================================= # +class SolarProximityCnst(Constraint): + """ """ + + def __init__(self, + closer: bool = False, + min_separation_deg: float = None, + max_separation_deg: float = None, + weight: float = 1 + ): + super().__init__(weight=weight) + self.closer = closer + self.min_separation_deg = min_separation_deg + self.max_separation_deg = max_separation_deg + + # --------------------------------------------------------- # + # ------------------------ Methods ------------------------ # + def get_score(self, indices: np.ndarray) -> np.ndarray: + """ """ + return np.mean( + np.where( + self.score[indices] > 0., + 1, + 0 + ) + ) + + # --------------------------------------------------------- # + # ----------------------- Internal ------------------------ # + def _evaluate(self, target: _Target, sun_position: SSTarget): + """ + """ + separation = target.target.separation(sun_position).deg + + # Set to NaN intervals where the separation is outside requested values + if not (self.min_separation_deg is None): + separation[separation < self.min_separation_deg] = np.nan + if not (self.max_separation_deg is None): + separation[separation > self.max_separation_deg] = np.nan + + if all(np.isnan(separation)): + log.warning( + "Constraint " + "cannot be satisfied over the given time range." + ) + self.score = separation[1:] * np.nan + + else: + sep_min = np.nanmin(separation) + sep_max = np.nanmax(separation) + # farther_normalized_score = (separation - sep_min) / (sep_max - sep_min) + # if self.closer: + # self.score = 1 - farther_normalized_score + # else: + # self.score = farther_normalized_score + if self.closer: + # Add .1 to avoid 0 + self.score = (sep_max + 0.1 - separation) / np.nanmax(sep_max + 0.1 - separation) + else: + self.score = separation / sep_max + + return self.score +# ============================================================= # +# ============================================================= # + # ============================================================= # # ------------------------ Constraints ------------------------ # @@ -1119,7 +1187,7 @@ def weights(self): # --------------------------------------------------------- # # ------------------------ Methods ------------------------ # - def evaluate(self, target, time, nslots=1, sun_elevation=None): + def evaluate(self, target, time, nslots=1, sun_elevation=None, sun_position=None): """ """ cnts = np.zeros((self.size, time.size - 1)) @@ -1131,6 +1199,8 @@ def evaluate(self, target, time, nslots=1, sun_elevation=None): cnts[i, :] = cnt(tuple(time), nslots) elif isinstance(cnt, NightTimeCnst): cnts[i, :] = cnt(sun_elevation) + elif isinstance(cnt, SolarProximityCnst): + cnts[i, :] = cnt(target, sun_position) else: unEvaluatedCnst += 1 if unEvaluatedCnst == self.size: @@ -1149,7 +1219,7 @@ def evaluate(self, target, time, nslots=1, sun_elevation=None): return self.score - def evaluate_various_nslots(self, target: _Target, time: Time, test_indices: List[np.ndarray], sun_elevation=None) -> List[np.ndarray]: + def evaluate_various_nslots(self, target: _Target, time: Time, test_indices: List[np.ndarray], sun_elevation=None, sun_position=None) -> List[np.ndarray]: scores = [] # Loop over each set of indices @@ -1157,9 +1227,19 @@ def evaluate_various_nslots(self, target: _Target, time: Time, test_indices: Lis # Copy target and only keep the time/coordinates around the boundaries of test_indices to optimize the computations min_idx, max_idx = np.min(indices_window), np.max(indices_window) + 1 - sliced_target = None if target is None else target[min_idx : max_idx] # Perform a copy without damaging the original object - time = time[min_idx : max_idx] - sun_elevation = None if sun_elevation is None else sun_elevation[min_idx : max_idx] + if target is None: + sliced_target = None + else: + sliced_target = target[min_idx : max_idx] # Perform a copy without damaging the original object + sliced_time = time[min_idx : max_idx] + if sun_elevation is None: + sun_elev = None + else: + sun_elev = sun_elevation[min_idx : max_idx] + if sun_position is None: + sun_pos = None + else: + sun_pos = sun_position[min_idx : max_idx] nslots = indices_window.size - 1 @@ -1171,9 +1251,11 @@ def evaluate_various_nslots(self, target: _Target, time: Time, test_indices: Lis if isinstance(constraint, TargetConstraint) and (sliced_target is not None): constraint_scores[constraint_i, :] = constraint(sliced_target, nslots) elif isinstance(constraint, ScheduleConstraint): - constraint_scores[constraint_i, :] = constraint(tuple(time), nslots) + constraint_scores[constraint_i, :] = constraint(tuple(sliced_time), nslots) elif isinstance(constraint, NightTimeCnst): - constraint_scores[constraint_i, :] = constraint(sun_elevation) + constraint_scores[constraint_i, :] = constraint(sun_elev[:-1]) # already taken at mid dt + elif isinstance(constraint, SolarProximityCnst) and (sliced_target is not None): + constraint_scores[constraint_i, :] = constraint(sliced_target, sun_pos[:-1]) else: n_unevaluated_constraints += 1 diff --git a/nenupy/schedule/geneticalgo.py b/nenupy/schedule/geneticalgo.py index a5ca0ed..439eba8 100644 --- a/nenupy/schedule/geneticalgo.py +++ b/nenupy/schedule/geneticalgo.py @@ -246,7 +246,7 @@ def evolve(self, # Show status if self.generation%modGen == 0: - log.debug( + log.info( f'Generation {self.generation}, ' f'best score: {self._bestScore}.' ) diff --git a/nenupy/schedule/schedule.py b/nenupy/schedule/schedule.py index d92c733..aa6847c 100644 --- a/nenupy/schedule/schedule.py +++ b/nenupy/schedule/schedule.py @@ -25,6 +25,7 @@ from astropy.table import Table import numpy as np from copy import deepcopy, copy +import itertools from nenupy.schedule.obsblocks import ( Block, @@ -160,7 +161,7 @@ def statusColor(self): # --------------------------------------------------------- # # ------------------------ Methods ------------------------ # - def evaluate_score(self, time, sun_elevation): + def evaluate_score(self, time, sun_elevation, sun_position): """ """ # Evaluate the target positions over time and compute the @@ -175,7 +176,8 @@ def evaluate_score(self, time, sun_elevation): target=self.target, time=time, nslots=self.nSlots, - sun_elevation=sun_elevation + sun_elevation=sun_elevation, + sun_position=sun_position ) log.debug( @@ -591,8 +593,8 @@ def _compute_slots_time_range(time_min: Time, time_max: Time, dt: TimeDelta): ) period = time_max - time_min - time_steps = int( np.ceil(period/dt) ) - dt_shifts = np.arange(time_steps)*dt + time_steps = int(np.floor(period / dt)) + dt_shifts = np.arange(time_steps) * dt slot_starts = time_min + dt_shifts slot_stops = slot_starts + dt @@ -676,19 +678,22 @@ def __init__( self.observation_blocks = ScheduleBlocks(dt=self.dt) self.reserved_blocks = None - # Store the Sun's elevation - sun = SSTarget.fromName('Sun') - sun.computePosition( - Time( - np.append( - self._startsJD, - self._startsJD[-1] + self.dt.jd - ), - format='jd' - ) - ) - elevation = sun.elevation.deg - self.sun_elevation = (elevation[1:] + elevation[:-1])/2 + # Store the Sun's elevation and position + sun = SSTarget.fromName("Sun") + # sun.computePosition( + # Time( + # np.append( + # self._startsJD, + # self._startsJD[-1] + self.dt.jd + # ), + # format="jd" + # ) + # ) + # elevation = sun.elevation.deg + # self.sun_elevation = (elevation[1:] + elevation[:-1]) / 2 + sun.computePosition(Time(self._startsJD + self.dt.jd / 2, format="jd")) + self.sun_elevation = sun.elevation.deg + self.sun_position = sun._fk5 self._obs_extension_done = False @@ -1009,7 +1014,8 @@ def book(self, optimize=False, **kwargs): """ # Ré-initialize the booking, in case the user calls this method several times - self._reset_observation_block_bookings() + if kwargs.get("reset_booking", True): + self._reset_observation_block_bookings() # Pre-compute the constraints scores self._compute_block_constraint_scores(**kwargs) @@ -1028,6 +1034,17 @@ def book(self, optimize=False, **kwargs): f'Fitting {sum(self._toSchedule)} observation blocks...' ) + # Remove blocks if they cannot be fully inserted + if kwargs.get("very_strict", False): + for i, blk in enumerate(self.observation_blocks): + if not self._toSchedule[i]: + continue + non_null_constraints = self._cnstScores[i] > 0 + consecutive_dt = np.array([ sum( 1 for _ in group ) for key, group in itertools.groupby(non_null_constraints) if key ]) + if not np.any(consecutive_dt >= blk.nSlots): + log.warning(f"Removed #{blk.blockIdx} '{blk.name}': no consecutive {blk.nSlots} slots found, 'very_strict' config!") + self._toSchedule[i] = False + if optimize: # All the observation blocks are booked at # once with the genetic algorithm. @@ -1094,6 +1111,9 @@ def book(self, optimize=False, **kwargs): elif kwargs.get('sort_by_availability', False): sort_idx = np.argsort(np.sum(self._cnstScores > 0, axis=1)) block_indices = block_indices[sort_idx] + elif kwargs.get("sort_by_nslots", False): + sort_idx = np.argsort([blk.nSlots for blk in self.observation_blocks])[::-1] + block_indices = block_indices[sort_idx] # Block are booked iteratively # for i, blk in enumerate(self.observation_blocks): @@ -1149,9 +1169,9 @@ def book(self, optimize=False, **kwargs): ) - def extend_scheduled_observations(self, relative_score_threshold: float = 0.8) -> None: + def extend_scheduled_observations(self, relative_score_threshold: float = 1., max_duration: TimeDelta = None) -> None: """ """ - if self._obs_extension_done: + if self._obs_extension_done and (max_duration is None): log.warning("Scheduled observations have already been extended.") return @@ -1173,8 +1193,15 @@ def extend_scheduled_observations(self, relative_score_threshold: float = 0.8) - for block in self.observation_blocks: if not block.isBooked: continue - extdt = block.max_extended_duration + + # If max_duration argument is set, use it instead of the ObsBlock value + if max_duration is None: + extdt = block.max_extended_duration + else: + extdt = max_duration + extend_mask.append(extdt is not None) + try: n_slots = int(np.ceil(extdt.sec/self.dt.sec)) except AttributeError: # extdt is None, and doesnt have .sec @@ -1184,8 +1211,13 @@ def extend_scheduled_observations(self, relative_score_threshold: float = 0.8) - # slot_extension.append(n_slots if extdt is not None else 0) # indices.append(block.blockIdx) + if block.nSlots == n_slots: + # Block is already at maximum size + slot_extension_indices.append([]) + continue + # Compute the constraint score for each possibility (i.e. from self.nSlots to n_slots) - extended_windows_size = np.arange(block.nSlots + 1, n_slots + 1) + 1 # Add 1 to evaluate the constraints + extended_windows_size = np.arange(block.nSlots + 1, n_slots + 1)# + 1 # Add 1 to evaluate the constraints extended_windows_slots = [np.arange(window_size) - int(window_size/2) + block.startIdx + int(block.nSlots/2) for window_size in extended_windows_size] # Deal with edges (don't check indices outside of the Schedule range) #extended_windows_slots = [window_indices[(window_indices >= 0) * (window_indices < self.size - 1)] for window_indices in extended_windows_slots] # self.size - 1 to compensate for the constraints evaluation @@ -1195,7 +1227,8 @@ def extend_scheduled_observations(self, relative_score_threshold: float = 0.8) - target=block.target, time=times, test_indices=extended_windows_slots, - sun_elevation=self.sun_elevation + sun_elevation=copy(self.sun_elevation), + sun_position=copy(self.sun_position) ) # Compute the mean score @@ -1209,12 +1242,17 @@ def treat(data): # Filter out the constraint score to discard those that diminish the block score # Add a relative threshold in order to filter... - positive_gradient_mask = constraint_score >= constraint_score[0] * relative_score_threshold + # Do not put >= as a score of 0 would be extended + positive_gradient_mask = constraint_score > constraint_score[0] * relative_score_threshold slot_extension_indices.append([slot_indices for slot_indices, better_option in zip(extended_windows_slots, positive_gradient_mask) if better_option]) # Iteratively try to increase the observation blocks iteration_i = 0 - max_iteration = max([len(indices) for indices in slot_extension_indices]) + try: + max_iteration = max([len(indices) for indices in slot_extension_indices]) + except ValueError: + # Nothing in slot_extension_indices + max_iteration = 0 while iteration_i < max_iteration: booked_block_i = 0 @@ -1648,7 +1686,11 @@ def _compute_block_constraint_scores(self, **kwargs): if blk.constraints.score is None: # If != 1, the constraint score has already been computed # Append the last time slot to get the last slot top - blk.evaluate_score(time=times, sun_elevation=self.sun_elevation) + blk.evaluate_score( + time=times, + sun_elevation=self.sun_elevation, + sun_position=self.sun_position + ) # Set other attributes relying on the schedule blk.nSlots = int(np.ceil(blk.duration/self.dt))