New SolarProximityCnst constraint, extend obs bug correction

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"
 

nenupy/schedule/__init__.py

@@ -14,6 +14,7 @@
     LocalTimeCnst,
     TimeRangeCnst,
     NightTimeCnst,
+    SolarProximityCnst,
     Constraints
 )
 from .obsblocks import Block, ObsBlock, ReservedBlock

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 <SolarProximityCnst(min_separation_deg="
+                f"{self.min_separation_deg}, max_separation_deg={self.max_separation_deg})> "
+                "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,17 +1219,27 @@ 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
         for indices_window in test_indices:
 
             # 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
 

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}.'
                 )

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 <ObsBlock> #{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))