diff --git a/requirements.txt b/requirements.txt index 1a430a6c..d4af55be 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,7 +2,7 @@ astropy hypothesis matplotlib more-itertools -numpy>=1.23.* +numpy>=1.23.*, <1.24 openpyxl pandas pytest diff --git a/scheduler/core/calculations/targetinfo.py b/scheduler/core/calculations/targetinfo.py index 952cfe23..965a5312 100644 --- a/scheduler/core/calculations/targetinfo.py +++ b/scheduler/core/calculations/targetinfo.py @@ -4,6 +4,7 @@ from dataclasses import dataclass from typing import Dict, Tuple +import numpy as np import numpy.typing as npt from astropy.coordinates import Angle, SkyCoord from astropy.time import TimeDelta @@ -52,6 +53,9 @@ class TargetInfo: rem_visibility_time: TimeDelta rem_visibility_frac: float + def mean_airmass(self, interval: npt.NDArray[int]): + return np.mean(self.airmass[interval]) + # Type aliases for TargetInfo information. TargetInfoNightIndexMap = Dict[NightIndex, TargetInfo] diff --git a/scheduler/core/components/optimizer/__init__.py b/scheduler/core/components/optimizer/__init__.py index d7b92107..918633c5 100644 --- a/scheduler/core/components/optimizer/__init__.py +++ b/scheduler/core/components/optimizer/__init__.py @@ -7,6 +7,7 @@ from scheduler.core.plans import Plans import numpy.typing as npt +from lucupy.minimodel import Program # Convenient type alias for Interval Interval = npt.NDArray[int] @@ -24,6 +25,7 @@ def __init__(self, selection: Selection, algorithm=None): # TODO: Assumes that all sites schedule the same amount of nights # if num_nights_optimize is None: self.period = len(list(self.night_events.values())[0].time_grid) + self.selection: Selection = selection # else: # self.period = num_nights_optimize @@ -32,3 +34,15 @@ def schedule(self) -> List[Plans]: nights = [Plans(self.night_events, night) for night in range(self.period)] self.algorithm.schedule(nights) return nights + + def _update_score(self, program: Program) -> None: + """Update the scores of the incomplete groups in the scheduled program""" + program_calculations = self.selection.score_program(program) + + for unique_group_id in program_calculations.top_level_groups: + group_data = program_calculations.group_data_map[unique_group_id] + group, group_info = group_data + schedulable_group = self.selection.schedulable_groups[unique_group_id] + # update scores in schedulable_groups if the group is not completely observed + if schedulable_group.group.exec_time() >= schedulable_group.group.total_used(): + schedulable_group.group_info.scores[:] = group_info.scores[:] \ No newline at end of file diff --git a/scheduler/core/components/optimizer/base.py b/scheduler/core/components/optimizer/base.py index 0af2fef2..3283f114 100644 --- a/scheduler/core/components/optimizer/base.py +++ b/scheduler/core/components/optimizer/base.py @@ -2,9 +2,9 @@ # For license information see LICENSE or https://opensource.org/licenses/BSD-3-Clause from abc import ABC, abstractmethod -from datetime import datetime -from typing import Mapping, List, Optional, Tuple - +from dataclasses import dataclass +from datetime import timedelta +from typing import Mapping, List, Optional, Union from lucupy.minimodel.program import ProgramID from scheduler.core.calculations.groupinfo import GroupData @@ -14,6 +14,20 @@ from . import Interval +@dataclass(frozen=True) +class MaxGroup: + """ + Store information about the selected group (max score) + """ + group_data: GroupData + max_score: float + interval: Interval + n_min: int + n_slots_remaining: int + n_std: int + exec_sci_nir: timedelta + + class BaseOptimizer(ABC): """ Base class for all Optimizer components. @@ -40,5 +54,6 @@ def setup(self, program_info: Mapping[ProgramID, ProgramInfo]): ... @abstractmethod - def add(self, group: GroupData, plans: Plans, interval: Optional[Interval] = None): + def add(self, night: int, max_group_info: Union[GroupData, MaxGroup]): ... + diff --git a/scheduler/core/components/optimizer/greedymax.py b/scheduler/core/components/optimizer/greedymax.py index bb9ef035..c5fe0f76 100644 --- a/scheduler/core/components/optimizer/greedymax.py +++ b/scheduler/core/components/optimizer/greedymax.py @@ -5,13 +5,13 @@ from dataclasses import dataclass from datetime import datetime, timedelta -from typing import Dict, FrozenSet, List, Optional, Tuple +from typing import Dict, FrozenSet, List, Optional, Tuple, Union from scheduler.core.calculations.selection import Selection from scheduler.core.calculations import GroupData from scheduler.core.plans import Plan, Plans from scheduler.core.components.optimizer.timeline import Timelines -from .base import BaseOptimizer +from .base import BaseOptimizer, MaxGroup from . import Interval from lucupy.minimodel import Program, Group, Observation, Sequence @@ -33,20 +33,6 @@ class ObsPlanData: visit_score: float -@dataclass(frozen=True) -class MaxGroup: - """ - Store information about the selected group (max score) - """ - group_data: GroupData - max_score: float - interval: Interval - n_min: int - n_slots_remaining: int - n_std: int - exec_sci_nir: timedelta - - class GreedyMaxOptimizer(BaseOptimizer): """ GreedyMax is an optimizer that schedules the visits for the rest of the night in a greedy fashion. @@ -82,11 +68,6 @@ def setup(self, selection: Selection) -> GreedyMaxOptimizer: self.obs_in_plan[site] = {} return self - # @staticmethod - # def progid(obsid) -> str: - # """Return program ID string from observation ID string""" - # return obsid[0:obsid.rfind('-')] - @staticmethod def non_zero_intervals(scores: npt.NDArray[float]) -> npt.NDArray[int]: """ @@ -135,7 +116,6 @@ def num_nir_standards(exec_sci, wavelengths=None, mode='spectroscopy') -> int: Calculated the number of NIR standards from the length of the NIR science and the mode """ n_std = 0 - time_per_standard = timedelta(0) # TODO: need mode or other info to distinguish imaging from spectroscopy if mode == 'imaging': @@ -167,11 +147,11 @@ def _exec_time_remaining(self, group: Group, verbose=False) -> Tuple[timedelta, Resource('IGRINS')] nsci = nprt = 0 - exec_remain = exec_remain_min = timedelta(0) - exec_sci = exec_sci_min = exec_sci_nir = timedelta(0) + + exec_sci_min = exec_sci_nir = timedelta(0) exec_prt = timedelta(0) time_per_standard = timedelta(0) - time_remain = time_remain_min = sci_times = timedelta(0) + sci_times = timedelta(0) n_std = 0 part_times = [] sci_times_min = [] @@ -272,7 +252,7 @@ def _min_slots_remaining(self, group: Group) -> Tuple[int, int, int, timedelta]: return n_min, n_slots_remaining, n_std, exec_sci_nir - def _find_max_group(self, plans: Plans) -> MaxGroup: + def _find_max_group(self, plans: Plans) -> Optional[MaxGroup]: """ Find the group with the max score in an open interval Returns None if there is no such group. @@ -293,8 +273,6 @@ def _find_max_group(self, plans: Plans) -> MaxGroup: n_min_list = [] n_std_list = [] exec_nir_list = [] - n_std = None - exec_time_nir = timedelta(0) # Make a list of scores in the remaining groups for group_data in self.group_data_list: @@ -349,7 +327,7 @@ def _find_max_group(self, plans: Plans) -> MaxGroup: max_slots_remaining = None max_n_std = None max_exec_nir = timedelta(0) - iscore_sort = None + if len(max_scores) > 0: # sort scores from high to low iscore_sort = np.flip(np.argsort(max_scores)) @@ -428,40 +406,29 @@ def _integrate_score(night_idx: int, start = idx end = start + max_group_info.n_slots_remaining - 1 - # print(f"Initial start end: {start} {end} {n_time} {end - start + 1}") - # Shift to window boundary if within minimum block time of edge. # If near both boundaries, choose boundary with higher score. score_start = scores[start] # score at start score_end = scores[end-1] # score at end delta_start = start - max_group_info.interval[0] # difference between start of window and block delta_end = max_group_info.interval[-1] - end # difference between end of window and block - # n_min, n_time_remaining = self._min_slots_remaining(max_group_info.group_data.group) - # print(f"delta_start: {delta_start}, delta_end: {delta_end}") - # print(f"score_start: {score_start}, score_end: {score_end}") + if delta_start < max_group_info.n_min and delta_end < max_group_info.n_min: if score_start > score_end and score_start > 0.0: - # print('a') start = max_group_info.interval[0] end = start + max_group_info.n_slots_remaining - 1 elif score_end > 0.0: - # print('b') start = max_group_info.interval[-1] - max_group_info.n_slots_remaining + 1 end = max_group_info.interval[-1] elif delta_start < max_group_info.n_min and score_start > 0.0: - # print('c') start = max_group_info.interval[0] end = start + max_group_info.n_slots_remaining - 1 elif delta_end < max_group_info.n_min and score_end > 0: - # print('d') start = max_group_info.interval[-1] - max_group_info.n_slots_remaining + 1 end = max_group_info.interval[-1] - # print(f"Shifted start end: {start} {end} {end - start + 1}") - # Make final list of indices for the highest scoring shifted sub-interval best_interval = np.arange(start=start, stop=end+1) - # print(f"len(best_interval): {len(best_interval)}") return best_interval @@ -475,8 +442,6 @@ def _find_group_position(self, night_idx: int, max_group_info: MaxGroup) -> Inte # n_time_remaining = int(np.ceil((time_remaining / self.time_slot_length))) # number of time slots # n_min, n_slots_remaining = self._min_slots_remaining(max_group_info.group_data.group) - # print('_find_group_position') - # print(max_group_info.n_slots_remaining, max_group_info.interval) if max_group_info.n_slots_remaining < len(max_group_info.interval): # Determine position based on max integrated score # If we don't end up here, then the group will have to be split later @@ -496,17 +461,10 @@ def nir_slots(self, science_obs, n_slots_filled, len_interval) -> Tuple[int, int slot_start_nir = None slot_end_nir = None slot_start = 0 - slot_end = 0 obs_id_nir = None for obs in science_obs: obs_id = obs.id - # print(f"Science observation: {obs_id.id}") - # index in observation list, for timeline - # iobs = self.obs_group_ids.index(obs.to_unique_group_id) - cumul_seq = self.cumulative_seq_exec_times(obs.sequence) - # print(f'len(cumm_seq) = {len(cumm_seq)}') - exec_time = timedelta(0) atom_start = self.first_nonzero_time(cumul_seq) atom_end = atom_start @@ -517,11 +475,9 @@ def nir_slots(self, science_obs, n_slots_filled, len_interval) -> Tuple[int, int while n_slots_filled + visit_length <= len_interval and atom_end <= len(cumul_seq) - 2: atom_end += 1 visit_length = n_slots_acq + \ - Plan.time2slots(self.time_slot_length, cumul_seq[atom_end]) + Plan.time2slots(self.time_slot_length, cumul_seq[atom_end]) slot_end = slot_start + visit_length - 1 - # print(f'slot_start = {slot_start} slot_end = {slot_end}') - # NIR science time for to determine the number of tellurics if any(inst in obs.required_resources() for inst in nir_inst): if slot_start_nir is None: @@ -554,17 +510,10 @@ def place_standards(self, night, interval, science_obs, partner_obs, n_std, verb standards = [] placement = [] - exec_sci_nir = timedelta(0) - - # obs_science = max_group_info.group_data.group.science_observations() - # obs_partner = max_group_info.group_data.group.partner_observations() - if verbose: print(f'Running place_standards') xdiff_min = xdiff_before_min = xdiff_after_min = 99999. - place_before = False - best_placement = False std_before = None std_after = None # If only one standard needed, try before or after, use best airmass match @@ -679,9 +628,6 @@ def _charge_time(observation: Observation, atom_start: int = 0, atom_end: int = if atom_end < 0: atom_end += seq_length - # print(f"Internal time charging") - # print(observation.id.id, atom_start, atom_end) - # Update observation status if atom_end == seq_length: observation.status = ObservationStatus.OBSERVED @@ -690,8 +636,6 @@ def _charge_time(observation: Observation, atom_start: int = 0, atom_end: int = for n_atom in range(atom_start, atom_end + 1): # "Charge" the expected program and partner times for the atoms: - # print(observation.id.id, n_atom, observation.sequence[n_atom].prog_time, - # observation.sequence[n_atom].part_time) observation.sequence[n_atom].program_used = observation.sequence[n_atom].prog_time observation.sequence[n_atom].partner_used = observation.sequence[n_atom].part_time @@ -716,7 +660,6 @@ def plot_airmass(self, obsid, interval=None, night=0) -> None: :return: None """ programid = obsid.program_id() - # print(obsid.id, programid.id) airmass = self.selection.program_info[programid].target_info[obsid][night].airmass @@ -761,7 +704,6 @@ def plot_timelines(self, night: int = 0) -> None: # fig, ax1 = plt.subplots(nrows=1, ncols=1, figsize=(10, 6)) obs_order = self.timelines[night][site].get_observation_order() for idx, istart, iend in obs_order: - obs_id = 'Unscheduled' if idx != -1: unique_group_id = self.obs_group_ids[idx] obs_id = ObservationID(unique_group_id.id) @@ -780,8 +722,6 @@ def plot_timelines(self, night: int = 0) -> None: ax2.plot(x[istart:iend + 1], np.log10(scores[istart:iend + 1]), linewidth=4, color=colour, label=obs_id.id) - # print(idx, obs_id, istart, iend) - # ax1.plot(self.timelines[night][site].time_slots) # ax1.axhline(0.0, xmax=1.0, color='black') ax1.axhline(2.0, xmax=1.0, color='black') @@ -859,8 +799,7 @@ def _add_visit(self, night, obs, max_group_info, best_interval, n_slots_filled) iobs = self.obs_group_ids.index(obs.to_unique_group_id) cumul_seq = self.cumulative_seq_exec_times(obs.sequence) - # print(f"len(cumm_seq) = {len(cumul_seq)}") - exec_time = timedelta(0) + atom_start = self.first_nonzero_time(cumul_seq) atom_end = atom_start @@ -876,8 +815,7 @@ def _add_visit(self, night, obs, max_group_info, best_interval, n_slots_filled) start_time_slot, start = timeline.add(iobs, visit_length, best_interval) # Get visit score and store information for the output plans - visit_score = sum(max_group_info.group_data.group_info.scores[night][start_time_slot:start_time_slot - + visit_length]) + visit_score = sum(max_group_info.group_data.group_info.scores[night][start_time_slot:start_time_slot+visit_length]) self.obs_in_plan[site][start_time_slot] = ObsPlanData( obs=obs, obs_start=start, @@ -893,7 +831,7 @@ def _add_visit(self, night, obs, max_group_info, best_interval, n_slots_filled) return n_slots_filled - def add(self, night: int, max_group_info: MaxGroup) -> bool: + def add(self, night: int, max_group_info: Union[GroupData, MaxGroup]) -> bool: """ Add a group to a Plan - find the best location within the interval (maximize the score) and select standards """ @@ -908,7 +846,6 @@ def add(self, night: int, max_group_info: MaxGroup) -> bool: timeline = self.timelines[night][site] program = self.selection.program_info[max_group_info.group_data.group.program_id].program # visit = [] # list of observations in visit - best_interval = None result = False print(f"Interval start end: {max_group_info.interval[0]} {max_group_info.interval[-1]}") @@ -916,8 +853,6 @@ def add(self, night: int, max_group_info: MaxGroup) -> bool: if not timeline.is_full: # Find the best location in timeline for the group best_interval = self._find_group_position(night, max_group_info) - # print(f"Interval start end: {max_group_info.interval[0]} {max_group_info.interval[-1]}") - # print(f"Best interval start end: {best_interval[0]} {best_interval[-1]}") if self.show_plots: self._plot_interval(max_group_info.group_data.group_info.scores[night], max_group_info.interval, @@ -927,18 +862,16 @@ def add(self, night: int, max_group_info: MaxGroup) -> bool: # When/If we eventually support splitting NIR observations, then we need to calculate the # NIR science time in best_interval and the number of basecal (e.g. telluric standards) needed. # This may require the calibration service for selecting standards. - # For now, we assume that NIR observations are not split and we use the telluric standards provided. + # For now, we assume that NIR observations are not split, and we use the telluric standards provided. # Pick standard(s) if needed n_slots_cal = 0 before_std = None after_std = None prog_obs = max_group_info.group_data.group.program_observations() - # print(f'sci_obs: {[o.id.id for o in prog_obs]}') part_obs = max_group_info.group_data.group.partner_observations() - # print(f'part_obs: {[o.id.id for o in part_obs]}') + if max_group_info.n_std > 0: - # print(f'\t n_std: {max_group_info.n_std} exec_sci_nir: {max_group_info.exec_sci_nir}') if max_group_info.exec_sci_nir > timedelta(0): standards, place_before = self.place_standards(night, best_interval, prog_obs, part_obs, max_group_info.n_std) @@ -958,7 +891,6 @@ def add(self, night: int, max_group_info: MaxGroup) -> bool: print(f"Adding before_std: {obs.to_unique_group_id} {obs.id.id}") n_slots_filled = self._add_visit(night, obs, max_group_info, best_interval, n_slots_filled) - # if n_slots_remaining > len(best_interval): # this would avoid splitting # split at atoms # Reserve space for the cals, otherwise the science observes will fill the interval n_slots_filled = n_slots_cal @@ -971,7 +903,7 @@ def add(self, night: int, max_group_info: MaxGroup) -> bool: print(f"Adding after_std: {obs.to_unique_group_id} {obs.id.id}") n_slots_filled = self._add_visit(night, obs, max_group_info, best_interval, n_slots_filled) - # TOTO: Shift to remove any gaps in the plan? + # TODO: Shift to remove any gaps in the plan? # Re-score program (pseudo time accounting) self._update_score(program, night=night) @@ -997,8 +929,6 @@ def output_plans(self, plans: Plans) -> None: # print('output_plans') # print(f'idx = {idx} len(obs_group_ids) = {len(self.obs_group_ids)}') # print(self.obs_group_ids) - unique_group_id = self.obs_group_ids[idx] - obs_id = ObservationID(unique_group_id.id) # Add visit to final plan obs_in_plan = self.obs_in_plan[timeline.site][start_time_slot] diff --git a/scheduler/core/plans/__init__.py b/scheduler/core/plans/__init__.py index 652ccfeb..a3fc46fd 100644 --- a/scheduler/core/plans/__init__.py +++ b/scheduler/core/plans/__init__.py @@ -4,10 +4,13 @@ from dataclasses import dataclass from datetime import datetime, timedelta from math import ceil -from typing import List, Mapping, Optional +from typing import List, Mapping, Optional, Tuple from lucupy.minimodel import Observation, ObservationID, Site, Conditions, Band, Resource +import numpy as np +import numpy.typing as npt +from scheduler.core.calculations import TargetInfoNightIndexMap from scheduler.core.calculations.nightevents import NightEvents @@ -22,6 +25,7 @@ class Visit: score: float instrument: Optional[Resource] + @dataclass(frozen=True) class NightStats: timeloss: str @@ -50,6 +54,139 @@ def __post_init__(self): self.night_stats: Optional[NightStats] = None self.alt_degs: List[List[float]] = [] + def nir_slots(self, + science_obs: List[Observation], + n_slots_filled: int, + len_interval: int) -> Tuple[int, int, ObservationID]: + """ + Return the starting and ending timeline slots (indices) for the NIR science observations + """ + # TODO: This should probably be moved to a more general location + nir_inst = [Resource('Flamingos2'), Resource('GNIRS'), Resource('NIRI'), Resource('NIFS'), + Resource('IGRINS')] + + # science, split at atom + slot_start_nir = None + slot_end_nir = None + slot_start = 0 + obs_id_nir = None + + for obs in science_obs: + obs_id = obs.id + + # TODO: currently in lucupy + cumul_seq = [] + atom_start = 0 + atom_end = atom_start + + n_slots_acq = Plan.time2slots(self.time_slot_length, obs.acq_overhead) + visit_length = n_slots_acq + Plan.time2slots(self.time_slot_length, + cumul_seq[atom_end]) + + # TODO: can this be done w/o a loop? convert cumm_seq to slots, and find the value that fits + while n_slots_filled + visit_length <= len_interval and atom_end <= len(cumul_seq) - 2: + atom_end += 1 + visit_length = n_slots_acq + \ + Plan.time2slots(self.time_slot_length, cumul_seq[atom_end]) + + slot_end = slot_start + visit_length - 1 + # NIR science time for to determine the number of tellurics + if any(inst in obs.required_resources() for inst in nir_inst): + if slot_start_nir is None: + slot_start_nir = slot_start + n_slots_acq # start of the science sequence, after acq + slot_end_nir = slot_end + obs_id_nir = obs_id + + n_slots_filled += visit_length + + slot_start = slot_end + 1 # for the next iteration + + return slot_start_nir, slot_end_nir, obs_id_nir + + def _place_standards(self, + interval: npt.NDArray[int], + science_obs: List[Observation], + partner_obs: List[Observation], + target_info: TargetInfoNightIndexMap, + night: int, + n_std) -> Tuple[List, List]: + """Pick the standards that best match the NIR science + observations by airmass + """ + + standards = [] + placement = [] + + # print(f'Running place_standards') + + xdiff_min = xdiff_before_min = xdiff_after_min = 99999. + std_before = None + std_after = None + # If only one standard needed, try before or after, use best airmass match + # TODO: Any preference to taking the standard before or after? + # TODO: Check scores to confirm that the observations are scheduleable (?) + for partcal_obs in partner_obs: + # Need the length of the calibration sequence only + n_slots_cal = Plan.time2slots(self.time_slot_length, partcal_obs.exec_time()) + n_slots_acq = Plan.time2slots(self.time_slot_length, partcal_obs.acq_overhead) + + # Try std first + # Mean std airmass + slot_start = n_slots_acq + slot_end = n_slots_cal - 1 + + xmean_cal = target_info[night][partcal_obs.id].mean_airmass(interval[slot_start:slot_end + 1]) + + # Mean NIR science airmass + idx_start_nir, idx_end_nir, obs_id_nir = self.nir_slots(science_obs, n_slots_cal, len(interval)) + slot_start_nir = slot_end + idx_start_nir + slot_end_nir = slot_end + idx_end_nir + + xmean_nir = target_info[night][obs_id_nir].mean_airmass(interval[slot_start_nir:slot_end_nir + 1]) + xdiff_before = np.abs(xmean_nir - xmean_cal) + + # Try std last + # Mean std airmass + len_int = len(interval) + slot_start = len_int - 1 - n_slots_cal + n_slots_acq + slot_end = slot_start + n_slots_cal - n_slots_acq - 1 + + xmean_cal = target_info[night][partcal_obs.id][night].mean_airmass(interval[slot_start:slot_end + 1]) + + # Mean NIR science airmass + slot_start_nir = idx_start_nir + slot_end_nir = idx_end_nir + + xmean_nir = target_info[night][obs_id_nir].mean_airmass(interval[slot_start_nir:slot_end_nir + 1]) + xdiff_after = np.abs(xmean_nir - xmean_cal) + + if n_std == 1: + if xdiff_before <= xdiff_after: + xdiff = xdiff_before + place_before = True # before + else: + xdiff = xdiff_after + place_before = False # after + + if xdiff < xdiff_min: + xdiff_min = xdiff + placement = [place_before] + standards = [partcal_obs] + else: + if xdiff_before < xdiff_before_min: + xdiff_before_min = xdiff_before + std_before = partcal_obs + + if xdiff_after < xdiff_after_min: + xdiff_after_min = xdiff_after + std_after = partcal_obs + # TODO: This should be added directly to the plan + if n_std > 1: + placement = [True, False] + standards = [std_before, std_after] + + return standards, placement + @staticmethod def time2slots(time_slot_length: timedelta, time: timedelta) -> int: # return ceil((time.total_seconds() / self.time_slot_length.total_seconds()) / 60) @@ -73,12 +210,12 @@ def add(self, self.visits.append(visit) self._time_slots_left -= time_slots - def __contains__(self, obs: Observation) -> bool: - return any(visit.obs_id == obs.id for visit in self.visits) - def time_left(self) -> int: return self._time_slots_left + def __contains__(self, obs: Observation) -> bool: + return any(visit.obs_id == obs.id for visit in self.visits) + class Plans: """