17651900636

Committed 11 Sep 2025 05:08PM UTC coverage: 35.141% (+4.8%) from 30.363%

Build # 17651900636

Build Type

Pull #190

github

Committed by

weaverba137

Commit Message

move requirements into setup.cfg

Pull Request Pull Request #190: Update test matrix for NumPy 2 compatibility

Run Details

1892 of 5384 relevant lines covered (35.14%)

0.35 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

61.03

/py/desisurvey/optimize.py

"""Optimize future DESI observations.
"""
from __future__ import print_function, division

import pkg_resources

import numpy as np
import scipy.special

import astropy.table
import astropy.coordinates
import astropy.units as u

import desiutil.log

import desisurvey.config


def wrap(angle, offset):
    """Wrap values in the range [0, 360] to [offset, offset+360].
    """
    return np.fmod(angle - offset + 360, 360) + offset


class Optimizer(object):
    """Initialize the hour angle assignments for specified tiles.

    See `DESI-3060
    <https://desi.lbl.gov/DocDB/cgi-bin/private/ShowDocument?docid=3060>`__
    for details.

    Parameters
    ----------
    condition : 'DARK', 'GRAY' or 'BRIGHT'
        Which obsconditions to optimize.
    lst_edges : array
        Array of N+1 LST bin edges.
    lst_hist : array
        Array of N bins giving the available LST distribution in units of
        sidereal hours per bin.
    subset : array or None
        An array of tile ID values to optimize.
        Optimizes all tiles with the relevant conditions if None.
    start : date or None
        Only consider available LST starting from this date.  Use the
        nominal survey start date if None.
    stop : date or None
        Only consider available LST before this end date.  Use the nominal
        survey stop date if None.
    init : 'zero', 'flat' or 'array'
        Method for initializing tile hour angles: 'zero' sets all hour angles
        to zero, 'flat' matches the CDF of available LST to planned LST (without
        accounting for exposure time), 'array' initializes from the initial_ha
        argument.
    initial_ha : array or None
        Only used when init is 'array'. The subset arg must also be provided
        to specify which tile each HA applies to.
    stretch : float
        Amount to stretch exposure times to account for factors other than
        dust and airmass (i.e., seeing, transparency, moon). This does not
        have a big effect on the results so can be approximate.
    smoothing_radius : :class:`astropy.units.Quantity`
        Gaussian sigma for calculating smoothing weights with angular units.
    center : float or None
        Used by the 'flat' initialization method to specify the starting
        DEC for the CDF balancing algorithm. When None, the 'flat' method
        scans over a grid of center values and picks the best one, but this
        is relatively slow.  Ignored unless init is 'flat'.
    seed : int or None
        Random number seed to use for stochastic elements of the optimizer.
        Do not use None if reproducible results are required.
    weights : array
        Array of relative weights to use when selecting which LST bin to
        optimize next.  Candidate bins are ordered by an estimated improvement.
        The length of the weights array determines how many candidates to
        consider, in decreasing order, and the weight values determines their
        relative weight. The next bin to optimize is then selected at random.
    completed : array
        Array of tileid, donefrac_{bright/gray/dark} indicating what fraction
        of particular tiles has already been observed in a particular
        condition.
    """
    def __init__(self, condition, lst_edges, lst_hist, subset=None, start=None, stop=None,
                 init='flat', initial_ha=None, stretch=1.0, smoothing_radius=10,
                 center=None, seed=123, weights=[5, 4, 3, 2, 1],
                 completed=None):

        tiles = desisurvey.tiles.get_tiles()
        if not isinstance(smoothing_radius, u.Quantity):
            smoothing_radius = smoothing_radius * u.deg
        self.log = desiutil.log.get_logger()
        config = desisurvey.config.Configuration()
        self.gen = np.random.RandomState(seed)
        self.cum_weights = np.asarray(weights, float).cumsum()
        self.cum_weights /= self.cum_weights[-1]
        self.stretch = stretch

        if start is None:
            start = config.first_day()
        else:
            start = desisurvey.utils.get_date(start)
        if stop is None:
            stop = config.last_day()
        else:
            stop = desisurvey.utils.get_date(stop)
        if start >= stop:
            raise ValueError('Expected start < stop.')

        self.lst_hist = np.asarray(lst_hist)
        self.lst_hist_sum = self.lst_hist.sum()
        self.nbins = len(lst_hist)
        self.lst_edges = np.asarray(lst_edges)
        if self.lst_edges.shape != (self.nbins + 1,):
            raise ValueError('Inconsistent lengths of lst_hist and lst_edges.')
        self.lst_centers = 0.5 * (self.lst_edges[1:] + self.lst_edges[:-1])
        self.origin = self.lst_edges[0]
        self.binsize = 360. / self.nbins

        # Select the tiles to plan.
        if subset is not None:
            idx = tiles.index(subset)
            # Check that all tiles in the subset are observable in these
            # conditions.
            if not np.all(tiles.allowed_in_conditions(condition)[idx]):
                raise ValueError('Subset contains non-{} tiles.'.format(condition))
            tile_sel = np.zeros(tiles.ntiles, bool)
            tile_sel[idx] = True
        else:
            # Use all tiles in the program by default.
            tile_sel = tiles.allowed_in_conditions(condition)

        # Get nominal exposure time for this program,
        # converted to LST equivalent in degrees.
        texp_nom = u.Quantity([
            getattr(config.programs, program).efftime()
            for program in tiles.tileprogram[tile_sel]])
        moon_up_factor = getattr(config.conditions, condition).moon_up_factor()
        texp_nom *= moon_up_factor
        if completed is not None:
            # completed = astropy.table.Table.read(completed)
            # donefrac = np.zeros(tiles.ntiles, dtype='f4')
            from astropy.io import fits
            completed = fits.getdata(completed)
            idx, mask = tiles.index(completed['TILEID'], return_mask=True)
            idx = idx[mask]
            boost = np.ones(tiles.ntiles, dtype='f4')
            boost[idx] = (completed['NNIGHT_NEEDED'][mask] -
                          completed['NNIGHT'][mask])
            boost = np.clip(boost, 0, np.inf)
            # donefrac[idx] = completed['donefrac'][mask]
            # boost = np.clip(1-donefrac, 0, 1)
            texp_nom *= boost[tile_sel]

        self.dlst_nom = 360 * texp_nom.to(u.day).value / 0.99726956583

        # Save arrays for the tiles to plan.
        self.ra = wrap(tiles.tileRA[tile_sel], self.origin)
        self.dec = tiles.tileDEC[tile_sel]
        self.tid = tiles.tileID[tile_sel]
        self.idx = np.where(tile_sel)[0]
        self.ntiles = len(self.ra)

        self.max_abs_ha = tiles.max_abs_ha[tile_sel]

        # Lookup static dust exposure factors for each tile.
        self.dust_factor = tiles.dust_factor[tile_sel]

        # Initialize smoothing weights.
        self.init_smoothing(smoothing_radius)

        self.log.info(
            '{0}: {1:.1f}h for {2} tiles (texp_nom {3:.1f}).'
            .format(condition, self.lst_hist_sum, self.ntiles,
                    np.median(texp_nom)))

        # Initialize metric histories.
        self.scale_history = []
        self.loss_history = []
        self.RMSE_history = []

        # Initialize improve() counters.
        self.nslow = 0
        self.nimprove = 0
        self.nsmooth = 0

        # Calculate schedule plan with HA=0 asignments to establish
        # the smallest possible total exposure time.
        self.plan_tiles = self.get_plan(np.zeros(self.ntiles))
        self.use_plan(save_history=False)
        self.min_total_time = self.plan_hist.sum()

        self.plan_tiles = self.get_plan(np.zeros(self.ntiles))
        self.use_plan()
        self.plan_hist_ha0 = self.plan_hist.copy()

        # Initialize HA assignments for each tile.
        if init == 'zero':
            self.ha = np.zeros(self.ntiles)
        elif init == 'array':
            if subset is None:
                raise ValueError('Must specify subset when init is "array".')
            if len(initial_ha) != self.ntiles:
                raise ValueError('Array initial_ha has wrong length.')
            self.ha = np.asarray(initial_ha)
        elif init == 'flat':
            if center is None:
                # Try 72 equally spaced centers.
                icenters = np.arange(0, self.nbins, self.nbins // 72)
            else:
                # Find the closest bin edge to the requested value.
                icenters = [np.argmin(np.abs(center - lst_edges))]
            min_score = np.inf
            scores = []
            for icenter in icenters:
                # Rotate the available LST histogram to this new center.
                hist = np.roll(self.lst_hist, -icenter)
                center = self.lst_edges[icenter]
                edges = np.linspace(center, center + 360, self.nbins + 1)
                # Calculate the CDF of available LST.
                lst_cdf = np.zeros_like(edges)
                lst_cdf[1:] = np.cumsum(hist)
                lst_cdf /= lst_cdf[-1]
                lst_cen = (edges[1:]+edges[:-1])/2
                lst_cdf = (lst_cdf[1:]+lst_cdf[:-1])/2
                # Calculate the CDF of planned LST usage relative to the same
                # central LST, assuming HA=0. Instead of spreading each exposure
                # over multiple LST bins, add its entire HA=0 exposure time at
                # LST=RA.
                exptime, _ = self.get_exptime(ha=np.zeros(self.ntiles))
                tile_ra = wrap(tiles.tileRA[tile_sel], center)
                sort_idx = np.argsort(tile_ra)
                tile_cdf = np.cumsum(exptime[sort_idx])
                tile_cdf /= tile_cdf[-1]
                tile_cdf = (np.concatenate([[0], tile_cdf[:-1]])+tile_cdf)/2
                # Use linear interpolation to find an LST for each tile that
                # matches the plan CDF to the available LST CDF.
                new_lst = np.interp(tile_cdf, lst_cdf, lst_cen)
                # Calculate each tile's HA as the difference between its HA=0
                # LST and its new LST after CDF matching.
                ha = np.empty(self.ntiles)
                ha[sort_idx] = np.fmod(new_lst - tile_ra[sort_idx], 360)
                # Clip tiles to their airmass limits.
                ha = np.clip(ha, -self.max_abs_ha, +self.max_abs_ha)
                # Calculate the score for this HA assignment.
                self.plan_tiles = self.get_plan(ha)
                self.use_plan(save_history=False)
                scores.append(self.eval_score(self.plan_hist))
                # Keep track of the best score found so far.
                if scores[-1] < min_score:
                    self.ha = ha.copy()
                    min_score = scores[-1]
                    center_best = center
            self.log.info(
                'Center flat initial HA assignments at LST {:.0f} deg.'
                .format(center_best))
        else:
            raise ValueError('Invalid init option: {0}.'.format(init))

        # Check initial HA assignments against airmass limits.
        ha_clipped = np.clip(self.ha, -self.max_abs_ha, +self.max_abs_ha)
        if not np.all(self.ha == ha_clipped):
            delta = np.abs(self.ha - ha_clipped)
            idx = np.argmax(delta)
            self.log.warn('Clipped {0} HA assignments to airmass limits.'
                          .format(np.count_nonzero(delta)))
            self.log.warn('Max clip is {0:.1f} deg for tile {1}.'
                          .format(delta[idx], self.tid[idx]))
            self.ha = ha_clipped

        # Calculate schedule plan with initial HA asignments.
        self.plan_tiles = self.get_plan(self.ha)
        self.use_plan()
        self.ha_initial = self.ha.copy()
        self.num_adjustments = np.zeros(self.ntiles, int)

    def get_exptime(self, ha, subset=None):
        """Estimate exposure times for the specified tiles.

        Estimates account for airmass and dust extinction only.

        Parameters
        ----------
        ha : array
            Array of hour angle assignments in degrees.
        subset : array or None
            Restrict calculation to a subset of tiles specified by the
            indices in this array, or use all tiles pass to the constructor
            if None.

        Returns
        -------
        tuple
            Tuple (exptime, subset) where exptime is an array of estimated
            exposure times in degrees and subset is the input subset or
            else a slice initialized for all tiles.
        """
        # Calculate for all tiles if no subset is specified.
        if subset is None:
            subset = slice(None)
        # Calculate best-case exposure times (i.e., using only airmass & dust)
        # in degrees for the specified HA assignments.
        tiles = desisurvey.tiles.get_tiles()
        X = tiles.airmass(ha, self.idx[subset])
        exptime = (self.dlst_nom[subset] *
                   desisurvey.etc.airmass_exposure_factor(X) *
                   self.dust_factor[subset]) * self.stretch
        return exptime, subset

    def get_plan(self, ha, subset=None):
        """Calculate an LST usage plan for specified hour angle assignments.

        Parameters
        ----------
        ha : array
            Array of hour angle assignments in degrees.
        subset : array or None
            Restrict calculation to a subset of tiles specified by the
            indices in this array, or use all tiles pass to the constructor
            if None.

        Returns
        -------
        array
            Array of shape (ntiles, nbins) giving the exposure time in hours
            that each tile needs in each LST bin. When a subset is specified,
            ntiles only indexes tiles in the subset.
        """
        exptime, subset = self.get_exptime(ha, subset)
        # Calculate LST windows for each tile's exposure.
        lst_mid = self.ra[subset] + ha
        lst_min = np.fmod(
            lst_mid - 0.5 * exptime - self.origin + 360, 360) + self.origin
        ##assert np.all(lst_min >= self.origin)
        ##assert np.all(lst_min < self.origin + 360)
        lst_max = np.fmod(
            lst_mid + 0.5 * exptime - self.origin + 360, 360) + self.origin
        ##assert np.all(lst_max >= self.origin)
        ##assert np.all(lst_max < self.origin + 360)
        # Calculate each exposure's overlap with each LST bin.
        lo = np.clip(
            self.lst_edges[1:] - lst_min[:, np.newaxis], 0, self.binsize)
        hi = np.clip(
            lst_max[:, np.newaxis] - self.lst_edges[:-1], 0, self.binsize)
        plan = lo + hi
        plan[lst_max >= lst_min] -= self.binsize
        ##assert np.allclose(plan.sum(axis=1), exptime)
        # Convert from degrees to sidereal hours.
        return plan * 24. / 360.

    def eval_score(self, plan_hist):
        """Evaluate the score that improve() tries to minimize.

        Score is calculated as 100 * RMSE + 100 * loss.

        Parameters
        ----------
        plan_hist : array
            Histogram of planned LST usage for all tiles.

        Returns
        -------
        float
            Score value.
        """
        return (
            100 * self.eval_RMSE(plan_hist) +
            100 * self.eval_loss(plan_hist))

    def eval_RMSE(self, plan_hist):
        """Evaluate the mean-squared error metric for the specified plan.

        This is the metric that :meth:`optimize` attempts to improve. It
        measures the similarity of the available and planned LST histogram
        shapes, but not their normalizations.  A separate :meth:`eval_scale`
        metric measures how efficiently the plan uses the available LST.

        The histogram of available LST is rescaled to the same total time (area)
        before calculating residuals relative to the planned LST usage.

        RMSE values are scaled by (10K / ntiles) so the absolute metric value
        is more consistent when ntiles is varied.

        Parameters
        ----------
        plan_hist : array
            Histogram of planned LST usage for all tiles.

        Returns
        -------
        float
            Mean squared error value.
        """
        # Rescale the available LST total time to the plan total time.
        plan_sum = plan_hist.sum()
        scale = plan_sum / self.lst_hist_sum
        residuals = plan_hist - scale * self.lst_hist
        return np.sqrt(residuals.dot(residuals) * self.nbins) / plan_sum

    def eval_loss(self, plan_hist):
        """Evaluate relative loss of current plan relative to HA=0 plan.

        Calculated as (T-T0)/T0 where T is the total exposure time of the
        current plan and T0 is the total exposure time of an HA=0 plan.

        Parameters
        ----------
        plan_hist : array
            Histogram of planned LST usage for all tiles.

        Returns
        -------
        float
            Loss factor.
        """
        return plan_hist.sum() / self.min_total_time - 1.0

    def eval_scale(self, plan_hist):
        """Evaluate the efficiency of the specified plan.

        Calculates the minimum scale factor applied to the available
        LST histogram so that the planned LST usage is always <= the scaled
        available LST histogram.  This value can be intepreted as the fraction
        of the available time required to complete all tiles.

        This metric is only loosely correlated with the RMSE metric, so provides
        a useful independent check that the optimization is producing the
        desired results.

        This metric is not well defined if any bin of the available LST
        histogram is empty, which indicates that some tiles will not be
        observable during the [start:stop] range being optimized.  In this case,
        only bins with some available LST are included in the scale calculation.

        Parameters
        ----------
        plan_hist : array
            Histogram of planned LST usage for all tiles.

        Returns
        -------
        float
            Scale factor.
        """
        nonzero = self.lst_hist > 0
        return (plan_hist[nonzero] / self.lst_hist[nonzero]).max()

    def use_plan(self, save_history=True):
        """Use the current plan and update internal arrays.

        Calculates the `plan_hist` arrays from the per-tile `plan_tiles` array,
        and records the current values of the RMSE and scale metrics.
        """
        self.plan_hist = self.plan_tiles.sum(axis=0)
        if not np.all(np.isfinite(self.plan_hist)):
            raise RuntimeError('Found invalid plan_tiles in use_plan().')
        if save_history:
            self.scale_history.append(self.eval_scale(self.plan_hist))
            self.loss_history.append(self.eval_loss(self.plan_hist))
            self.RMSE_history.append(self.eval_RMSE(self.plan_hist))

    def next_bin(self):
        """Select which LST bin to adjust next.

        The algorithm determines which bin of the planned LST usage histogram
        should be decreased in order to maximize the decrease of the score,
        assuming that the decrease is moved to one of the neighboring bins.

        Since each tile's contribution to the plan can, in general, span several
        LST bins and can change its area (exposure time) when its HA is
        adjusted, the assumptions of this algorithm are not valid in detail
        but it usually does a good job anyway.

        This algorithm has a stochastic component controlled by the `weights`
        parameter passed to our constructor, in order to avoid getting stuck
        in a local minimum.

        Returns
        -------
        tuple
            Tuple (idx, dha_sign) where idx is the LST bin index that should
            be decreased (by moving one of the tiles contributing to it) and
            dha_sign gives the sign +/-1 of the HA adjustment required.
        """
        # Rescale the available LST total time to the current plan total time.
        A = self.lst_hist * self.plan_hist.sum() / self.lst_hist_sum
        # Calculate residuals in each LST bin.
        P = self.plan_hist
        res = P - A
        # Calculate the change from moving tiles between adjacent bins.
        nbins = len(self.lst_hist)
        adjacent = np.roll(np.arange(nbins), -1)
        dres = res[adjacent] - res
        # Cannot move tiles from an empty bin.
        empty = (P == 0)
        dres[empty & (dres < 0)] = 0.
        dres[empty[adjacent] & (dres > 0)] = 0.
        # Select the movements that reduce the RMSE between A and P by
        # the largest amounts.
        order = np.argsort(np.abs(dres))[::-1][:len(self.cum_weights)]
        # Randomly select one of these movements, according to our weights.
        which = np.searchsorted(self.cum_weights, self.gen.uniform())
        idx = order[which]
        if dres[idx] == 0:
            raise RuntimeError('Cannot improve RMSE.')
        elif dres[idx] > 0:
            idx = adjacent[idx]
            dha_sign = -1
        else:
            dha_sign = +1
        return idx, dha_sign

    def improve(self, frac=1.):
        """Perform one iteration of improving the hour angle assignments.

        Each call will adjust the HA of a single tile with a magnitude \|dHA\|
        specified by the `frac` parameter.

        Parameters
        ----------
        frac : float
            Mean fraction of an LST bin to adjust the selected tile's HA by.
            Actual HA adjustments are randomly distributed around this mean
            to smooth out adjustments.
        """
        # Randomly perturb the size of the HA adjustment.  This adds some
        # noise but also makes it possible to get out of dead ends.
        frac = np.minimum(2., self.gen.rayleigh(
            scale=np.sqrt(2 / np.pi) * frac))
        # Calculate the initial score.
        initial_score = self.eval_score(self.plan_hist)
        # Try a fast method first, then fall back to a slower method.
        for method in 'next_bin', 'any_bin':
            if method == 'next_bin':
                # Select which bin to move a tile from and in which direction.
                ibin, dha_sign = self.next_bin()
                # Find tiles using more time in this LST bin than in the
                # adjacent bin they would be moving away from.
                ibin_from = (ibin - dha_sign + self.nbins) % self.nbins
                sel = (self.plan_tiles[:, ibin] > 0) & (
                    self.plan_tiles[:, ibin] >= self.plan_tiles[:, ibin_from])
            else:
                # Try a 20% random subset of all tiles.
                self.nslow += 1
                sel = self.gen.permutation(self.ntiles) < (self.ntiles // 5)
                # Randomly select a direction to shift the next tile.
                dha_sign = +1 if self.gen.uniform() > 0.5 else -1

            # Do not move any tiles that are already at their |HA| limits.
            dha = 360. / self.nbins * frac * dha_sign
            veto = np.abs(self.ha + dha) >= self.max_abs_ha
            sel[veto] = False
            # Are there any tiles available to adjust?
            nsel = np.count_nonzero(sel)
            if nsel == 0:
                self.log.debug('No tiles available for {0} method.'
                               .format(method))
                continue
            # How many times have these tiles already been adjusted?
            nadj = self.num_adjustments[sel]
            if np.min(nadj) < np.max(nadj):
                # Do not adjust a tile that already has received the max
                # number of adjustments.  This has the effect of smoothing
                # the spatial distribution of HA adjustments.
                sel = sel & (self.num_adjustments < np.max(nadj))
                nsel = np.count_nonzero(sel)
            subset = np.where(sel)[0]
            # Calculate how the plan changes by moving each selected tile.
            scenario = self.get_plan(self.ha[subset] + dha, subset)
            new_score = np.zeros(nsel)
            for i, itile in enumerate(subset):
                # Calculate the (downsampled) plan when this tile is moved.
                new_plan_hist = (
                    self.plan_hist - self.plan_tiles[itile] + scenario[i])
                new_score[i]= self.eval_score(new_plan_hist)
            i = np.argmin(new_score)
            if new_score[i] > initial_score:
                # All candidate adjustments give a worse score.
                continue
            # Accept the tile that gives the smallest score.
            itile = subset[i]
            self.num_adjustments[itile] += 1
            # Update the plan.
            self.ha[itile] = self.ha[itile] + dha
            assert np.abs(self.ha[itile]) < self.max_abs_ha[itile]
            self.plan_tiles[itile] = scenario[i]
            # No need to try additional methods.
            break
        self.use_plan()
        self.nimprove += 1

    def init_smoothing(self, radius):
        """Calculate and save smoothing weights.

        Weights for each pair of tiles [i,j] are calculated as::

            wgt[i,j] = exp(-0.5 * (sep[i,j]/radius) ** 2)

        where sep[i,j] is the separation angle between the tile centers.

        Parameters
        ----------
        radius : astropy.units.Quantity
            Gaussian sigma for calculating weights with angular units.
        """
        separations = desisurvey.utils.separation_matrix(
            self.ra, self.dec, self.ra, self.dec)
        ratio = separations / radius.to(u.deg).value
        self.smoothing_weights = np.exp(-0.5 * ratio ** 2)
        # Set self weight to zero.
        self_weights = np.diag(self.smoothing_weights)
        assert np.allclose(self_weights, 1.)
        self.smoothing_weights -= np.diag(self_weights)
        self.smoothing_sums = self.smoothing_weights.sum(axis=1)

    def smooth(self, alpha=0.1):
        """Smooth the current HA assignments.

        Each HA is replaced with a smoothed value::

            (1-alpha) * HA + alpha * HA[avg]

        where HA[avg] is the weighted average of all other tile HA assignments.
        """
        avg_ha = self.smoothing_weights.dot(self.ha) / self.smoothing_sums
        self.ha = (1 - alpha) * self.ha + alpha * avg_ha
        self.plan_tiles = self.get_plan(self.ha)
        self.use_plan()

    def plot(self, save=None, relative=True):
        """Plot the current optimzation status.

        Requires that matplotlib is installed.

        Parameters
        ----------
        save : str or None
            Filename where the generated plot will be saved.
        """
        import matplotlib.pyplot as plt
        fig, ax = plt.subplots(2, 2, figsize=(12, 8))
        ax = ax.flatten()
        scale = self.scale_history[-1]
        ax[0].hist(self.lst_centers, bins=self.lst_edges,
                 weights=self.lst_hist * scale, histtype='stepfilled',
                 fc=(1,0,0,0.25), ec='r')
        ax[0].hist(self.lst_centers, bins=self.lst_edges,
                 weights=self.plan_hist, histtype='stepfilled',
                 fc=(0.7,0.7,1), ec='b')
        RMSE_scale = self.plan_hist.sum() / self.lst_hist_sum
        ax[0].hist(self.lst_centers, bins=self.lst_edges,
                 weights=self.lst_hist * RMSE_scale, histtype='step',
                 fc=(1,0,0,0.25), ec='r', ls='--')

        nbins = len(self.lst_centers)
        idx, dha = self.next_bin()
        ax[0].axvline(self.lst_centers[idx], c='b', ls='-')
        ax[0].axvline(
            self.lst_centers[(idx + nbins + dha) % nbins], c='b', ls='--')

        ax[0].set_xlim(self.lst_edges[0], self.lst_edges[-1])
        ax[0].set_ylim(0, None)
        ax[0].set_xlabel('Local Sidereal Time [deg]')
        ax[0].set_ylabel('Hours / LST Bin')

        ax[1].plot(self.RMSE_history, 'r.', ms=10, label='RMSE')
        ax[1].legend(loc='lower left', numpoints=1)
        rhs = ax[1].twinx()
        rhs.plot(self.scale_history, 'kx', ms=5, label='Scale')
        #rhs.plot(self.loss_history, 'g+', ms=5, label='Loss')
        rhs.legend(loc='upper right', numpoints=1)
        ax[1].set_xlim(-0.1, len(self.RMSE_history) - 0.9)
        ax[1].set_xlabel('Iterations')

        ax[2].hist(self.ha, bins=50, histtype='stepfilled')
        ax[2].set_xlabel('Tile Design Hour Angle [deg]')

        if relative:
            c = self.ha - self.ha_initial
            clabel = 'Tile Hour Angle Adjustments [deg]'
        else:
            c = self.ha
            clabel = 'Tile Hour Angle [deg]'
        s = ax[3].scatter(self.ra, self.dec, c=c, s=12, lw=0, cmap='jet')
        ax[3].set_xlim(self.lst_edges[0] - 5, self.lst_edges[-1] + 5)
        ax[3].set_xticks([])
        ax[3].set_xlabel(clabel)
        ax[3].set_ylim(-20, 80)
        ax[3].set_yticks([])
        cbar = plt.colorbar(s, ax=ax[3], orientation='vertical',
                            fraction=0.05, pad=0.01, format='%.1f')
        cbar.ax.tick_params(labelsize=9)
        plt.tight_layout()
        if save:
            plt.savefig(save)
        else:
            plt.show()


if __name__ == '__main__':
    """This should eventually be made into a first-class script entry point.
    """
    opt = Optimizer(scheduler, 'GRAY')
    for i in range(10):
        opt.improve(frac=0.25, verbose=True)

desihub / desisurvey / 17651900636

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous