11979115482

Committed 22 Nov 2024 07:45PM UTC coverage: 65.499% (-0.2%) from 65.737%

Build # 11979115482

Build Type

Pull #403

github

Committed by

web-flow

Commit Message

Merge 0d3dd536f into 86cc76cd9

Pull Request Pull Request #403: Improved version tracking and deprecated SVN

Run Details

36 of 54 new or added lines in 2 files covered. (66.67%)

43 existing lines in 5 files now uncovered.

9534 of 14556 relevant lines covered (65.5%)

0.65 hits per line

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

75.77

/EXOSIMS/SurveySimulation/linearJScheduler_det_only.py

from EXOSIMS.SurveySimulation.linearJScheduler import linearJScheduler
import astropy.units as u
import numpy as np
import time
import astropy.constants as const
from EXOSIMS.util._numpy_compat import copy_if_needed


class linearJScheduler_det_only(linearJScheduler):
    """linearJScheduler_det_only - linearJScheduler Detections Only

    This class implements the linear cost function scheduler described
    in Savransky et al. (2010).

    This scheduler inherits from the linearJScheduler module but performs only
    detections.

    Args:
        specs:
            user specified values

    """

    def __init__(self, **specs):

        linearJScheduler.__init__(self, **specs)

    def run_sim(self):
        """Performs the survey simulation"""

        OS = self.OpticalSystem
        TL = self.TargetList
        SU = self.SimulatedUniverse
        Obs = self.Observatory
        TK = self.TimeKeeping

        # TODO: start using this self.currentSep
        # set occulter separation if haveOcculter
        if OS.haveOcculter:
            self.currentSep = Obs.occulterSep

        # choose observing modes selected for detection (default marked with a flag)
        allModes = OS.observingModes
        det_mode = list(filter(lambda mode: mode["detectionMode"], allModes))[0]

        # begin Survey, and loop until mission is finished
        log_begin = "OB%s: survey beginning." % (TK.OBnumber)
        self.logger.info(log_begin)
        self.vprint(log_begin)
        t0 = time.time()
        sInd = None
        ObsNum = 0
        while not TK.mission_is_over(OS, Obs, det_mode):

            # acquire the NEXT TARGET star index and create DRM
            old_sInd = sInd  # used to save sInd if returned sInd is None
            DRM, sInd, det_intTime, waitTime = self.next_target(sInd, det_mode)

            if sInd is not None:
                ObsNum += (
                    1  # we're making an observation so increment observation number
                )

                if OS.haveOcculter:
                    # advance to start of observation
                    # (add slew time for selected target)
                    _ = TK.advanceToAbsTime(TK.currentTimeAbs.copy() + waitTime)

                # beginning of observation, start to populate DRM
                DRM["star_ind"] = sInd
                DRM["star_name"] = TL.Name[sInd]
                DRM["arrival_time"] = TK.currentTimeNorm.to("day").copy()
                DRM["OB_nb"] = TK.OBnumber
                DRM["ObsNum"] = ObsNum
                pInds = np.where(SU.plan2star == sInd)[0]
                DRM["plan_inds"] = pInds.astype(int)
                log_obs = (
                    "  Observation #%s, star ind %s (of %s) with %s planet(s), "
                    + "mission time at Obs start: %s"
                ) % (
                    ObsNum,
                    sInd,
                    TL.nStars,
                    len(pInds),
                    TK.currentTimeNorm.to("day").copy().round(2),
                )
                self.logger.info(log_obs)
                self.vprint(log_obs)

                # PERFORM DETECTION and populate revisit list attribute
                (
                    detected,
                    det_fZ,
                    det_systemParams,
                    det_SNR,
                    FA,
                ) = self.observation_detection(sInd, det_intTime, det_mode)
                # update the occulter wet mass
                if OS.haveOcculter:
                    DRM = self.update_occulter_mass(DRM, sInd, det_intTime, "det")
                # populate the DRM with detection results
                DRM["det_time"] = det_intTime.to("day")
                DRM["det_status"] = detected
                DRM["det_SNR"] = det_SNR
                DRM["det_fZ"] = det_fZ.to("1/arcsec2")
                DRM["det_params"] = det_systemParams

                # populate the DRM with observation modes
                DRM["det_mode"] = dict(det_mode)
                del DRM["det_mode"]["inst"], DRM["det_mode"]["syst"]

                DRM["exoplanetObsTime"] = TK.exoplanetObsTime.copy()

                # append result values to self.DRM
                self.DRM.append(DRM)

                # handle case of inf OBs and missionPortion < 1
                if np.isinf(TK.OBduration) and (TK.missionPortion < 1):
                    self.arbitrary_time_advancement(
                        TK.currentTimeNorm.to("day").copy() - DRM["arrival_time"]
                    )

            else:  # sInd == None
                sInd = old_sInd  # Retain the last observed star
                if (
                    TK.currentTimeNorm.copy() >= TK.OBendTimes[TK.OBnumber]
                ):  # currentTime is at end of OB
                    # Conditional Advance To Start of Next OB
                    if not TK.mission_is_over(
                        OS, Obs, det_mode
                    ):  # as long as the mission is not over
                        TK.advancetToStartOfNextOB()  # Advance To Start of Next OB
                elif waitTime is not None:
                    # CASE 1: Advance specific wait time
                    _ = TK.advanceToAbsTime(TK.currentTimeAbs.copy() + waitTime)
                    self.vprint("waitTime is not None")
                else:
                    startTimes = (
                        TK.currentTimeAbs.copy() + np.zeros(TL.nStars) * u.d
                    )  # Start Times of Observations
                    observableTimes = Obs.calculate_observableTimes(
                        TL,
                        np.arange(TL.nStars),
                        startTimes,
                        self.koMap,
                        self.koTimes,
                        det_mode,
                    )[0]

                    # CASE 2 If There are no observable targets for the rest
                    # of the mission
                    if (
                        observableTimes[
                            (
                                TK.missionFinishAbs.copy().value * u.d
                                > observableTimes.value * u.d
                            )
                            * (
                                observableTimes.value * u.d
                                >= TK.currentTimeAbs.copy().value * u.d
                            )
                        ].shape[0]
                    ) == 0:
                        self.vprint(
                            (
                                "No Observable Targets for Remainder of mission at "
                                "currentTimeNorm = {}"
                            ).format(TK.currentTimeNorm)
                        )
                        # Manually advancing time to mission end
                        TK.currentTimeNorm = TK.missionLife
                        TK.currentTimeAbs = TK.missionFinishAbs
                    else:
                        # CASE 3  nominal wait time if at least 1 target is still
                        # in list and observable
                        # TODO: ADD ADVANCE TO WHEN FZMIN OCURS
                        inds1 = np.arange(TL.nStars)[
                            observableTimes.value * u.d
                            > TK.currentTimeAbs.copy().value * u.d
                        ]
                        # apply intTime filter
                        inds2 = np.intersect1d(self.intTimeFilterInds, inds1)
                        # apply revisit Filter #NOTE this means stars you added
                        # to the revisit list
                        inds3 = self.revisitFilter(
                            inds2, TK.currentTimeNorm.copy() + self.dt_max.to(u.d)
                        )
                        self.vprint(
                            "Filtering %d stars from advanceToAbsTime"
                            % (TL.nStars - len(inds3))
                        )
                        oTnowToEnd = observableTimes[inds3]
                        # there is at least one observableTime between now and the
                        # end of the mission
                        if not oTnowToEnd.value.shape[0] == 0:
                            tAbs = np.min(oTnowToEnd)  # advance to that observable time
                        else:
                            # advance to end of mission
                            tAbs = TK.missionStart + TK.missionLife
                        tmpcurrentTimeNorm = TK.currentTimeNorm.copy()
                        # Advance Time to this time OR start of next OB following
                        # this time
                        _ = TK.advanceToAbsTime(tAbs)
                        self.vprint(
                            (
                                "No Observable Targets a currentTimeNorm= {:.2f} "
                                "Advanced To currentTimeNorm= {:.2f}"
                            ).format(
                                tmpcurrentTimeNorm.to("day"),
                                TK.currentTimeNorm.to("day"),
                            )
                        )
        else:  # TK.mission_is_over()
            dtsim = (time.time() - t0) * u.s
            log_end = (
                "Mission complete: no more time available.\n"
                + "Simulation duration: %s.\n" % dtsim.astype("int")
                + "Results stored in SurveySimulation.DRM (Design Reference Mission)."
            )
            self.logger.info(log_end)
            print(log_end)

    def next_target(self, old_sInd, mode):
        """Finds index of next target star and calculates its integration time.

        This method chooses the next target star index based on which
        stars are available, their integration time, and maximum completeness.
        Returns None if no target could be found.

        Args:
            old_sInd (integer):
                Index of the previous target star
            mode (dict):
                Selected observing mode for detection

        Returns:
            tuple:
                DRM (dict):
                    Design Reference Mission, contains the results of one complete
                    observation (detection and characterization)
                sInd (integer):
                    Index of next target star. Defaults to None.
                intTime (astropy Quantity):
                    Selected star integration time for detection in units of day.
                    Defaults to None.
                waitTime (astropy Quantity):
                    a strategically advantageous amount of time to wait in the case of
                    an occulter for slew times

        """
        OS = self.OpticalSystem
        TL = self.TargetList
        Obs = self.Observatory
        TK = self.TimeKeeping

        # create DRM
        DRM = {}

        # selecting appropriate koMap
        koMap = self.koMaps[mode["syst"]["name"]]

        # allocate settling time + overhead time
        tmpCurrentTimeAbs = (
            TK.currentTimeAbs.copy() + Obs.settlingTime + mode["syst"]["ohTime"]
        )
        tmpCurrentTimeNorm = (
            TK.currentTimeNorm.copy() + Obs.settlingTime + mode["syst"]["ohTime"]
        )

        # look for available targets
        # 1. initialize arrays
        slewTimes = np.zeros(TL.nStars) * u.d
        # fZs = np.zeros(TL.nStars) / u.arcsec**2
        dV = np.zeros(TL.nStars) * u.m / u.s
        intTimes = np.zeros(TL.nStars) * u.d
        obsTimes = np.zeros([2, TL.nStars]) * u.d
        sInds = np.arange(TL.nStars)

        # 2. find spacecraft orbital START positions (if occulter, positions
        # differ for each star) and filter out unavailable targets
        sd = None
        if OS.haveOcculter:
            sd = Obs.star_angularSep(TL, old_sInd, sInds, tmpCurrentTimeAbs)
            obsTimes = Obs.calculate_observableTimes(
                TL, sInds, tmpCurrentTimeAbs, self.koMaps, self.koTimes, mode
            )
            slewTimes = Obs.calculate_slewTimes(
                TL, old_sInd, sInds, sd, obsTimes, tmpCurrentTimeAbs
            )

        # 2.1 filter out totTimes > integration cutoff
        if len(sInds.tolist()) > 0:
            sInds = np.intersect1d(self.intTimeFilterInds, sInds)

        # start times, including slew times
        startTimes = tmpCurrentTimeAbs.copy() + slewTimes
        startTimesNorm = tmpCurrentTimeNorm.copy() + slewTimes

        # 2.5 Filter stars not observable at startTimes
        try:
            koTimeInd = np.where(
                np.round(startTimes[0].value) - self.koTimes.value == 0
            )[0][
                0
            ]  # find indice where koTime is startTime[0]
            sInds = sInds[
                np.where(np.transpose(koMap)[koTimeInd].astype(bool)[sInds])[0]
            ]  # filters inds by koMap #verified against v1.35
        except:  # noqa: E722 If there are no target stars to observe
            sInds = np.asarray([], dtype=int)

        # 3. filter out all previously (more-)visited targets, unless in
        if len(sInds.tolist()) > 0:
            sInds = self.revisitFilter(sInds, tmpCurrentTimeNorm)

        # 4.1 calculate integration times for ALL preselected targets
        (
            maxIntTimeOBendTime,
            maxIntTimeExoplanetObsTime,
            maxIntTimeMissionLife,
        ) = TK.get_ObsDetectionMaxIntTime(Obs, mode)
        maxIntTime = min(
            maxIntTimeOBendTime, maxIntTimeExoplanetObsTime, maxIntTimeMissionLife
        )  # Maximum intTime allowed

        if len(sInds.tolist()) > 0:
            intTimes[sInds] = self.calc_targ_intTime(sInds, startTimes[sInds], mode)
            sInds = sInds[
                np.where(intTimes[sInds] <= maxIntTime)
            ]  # Filters targets exceeding end of OB
            endTimes = startTimes + intTimes

            if maxIntTime.value <= 0:
                sInds = np.asarray([], dtype=int)

        # 5.1 TODO Add filter to filter out stars entering and exiting keepout
        # between startTimes and endTimes

        # 5.2 find spacecraft orbital END positions (for each candidate target),
        # and filter out unavailable targets
        if len(sInds.tolist()) > 0 and Obs.checkKeepoutEnd:
            try:
                # endTimes may exist past koTimes so we have an exception to
                # handle this case
                # koTimeInd[0][0]  # find indice where koTime is endTime[0]
                koTimeInd = np.where(
                    np.round(endTimes[0].value) - self.koTimes.value == 0
                )[0][0]
                # filters inds by koMap #verified against v1.35
                sInds = sInds[
                    np.where(np.transpose(koMap)[koTimeInd].astype(bool)[sInds])[0]
                ]
            except:  # noqa: E722
                sInds = np.asarray([], dtype=int)

        # 6. choose best target from remaining
        if len(sInds.tolist()) > 0:
            # choose sInd of next target
            sInd, waitTime = self.choose_next_target(
                old_sInd, sInds, slewTimes, intTimes[sInds]
            )

            # Should Choose Next Target decide there are no stars it wishes to
            # observe at this time.
            if (sInd is None) and (waitTime is not None):
                self.vprint(
                    (
                        "There are no stars Choose Next Target would like to Observe. "
                        "Waiting {}"
                    ).format(waitTime)
                )
                return DRM, None, None, waitTime
            elif (sInd is None) and (waitTime is not None):
                self.vprint(
                    (
                        "There are no stars Choose Next Target would like to Observe "
                        "and waitTime is None"
                    )
                )
                return DRM, None, None, waitTime
            # store selected star integration time
            intTime = intTimes[sInd]

        # if no observable target, advanceTime to next Observable Target
        else:
            self.vprint(
                "No Observable Targets at currentTimeNorm= "
                + str(TK.currentTimeNorm.copy())
            )
            return DRM, None, None, None

        # update visited list for selected star
        self.starVisits[sInd] += 1
        # store normalized start time for future completeness update
        self.lastObsTimes[sInd] = startTimesNorm[sInd]

        # populate DRM with occulter related values
        if OS.haveOcculter:
            DRM = Obs.log_occulterResults(
                DRM, slewTimes[sInd], sInd, sd[sInd], dV[sInd]
            )
            return DRM, sInd, intTime, slewTimes[sInd]

        return DRM, sInd, intTime, waitTime

    def choose_next_target(self, old_sInd, sInds, slewTimes, intTimes):
        """Choose next target based on truncated depth first search
        of linear cost function.

        Args:
            old_sInd (integer):
                Index of the previous target star
            sInds (integer array):
                Indices of available targets
            slewTimes (astropy quantity array):
                slew times to all stars (must be indexed by sInds)
            intTimes (astropy Quantity array):
                Integration times for detection in units of day

        Returns:
            sInd (integer):
                Index of next target star

        """

        Comp = self.Completeness
        TL = self.TargetList
        TK = self.TimeKeeping
        OS = self.OpticalSystem
        Obs = self.Observatory
        allModes = OS.observingModes

        # cast sInds to array
        sInds = np.array(sInds, ndmin=1, copy=copy_if_needed)

        # current star has to be in the adjmat
        if (old_sInd is not None) and (old_sInd not in sInds):
            sInds = np.append(sInds, old_sInd)

        # calculate dt since previous observation
        dt = TK.currentTimeNorm + slewTimes[sInds] - self.lastObsTimes[sInds]
        # get dynamic completeness values
        comps = Comp.completeness_update(TL, sInds, self.starVisits[sInds], dt)

        # if first target, or if only 1 available target,
        # choose highest available completeness
        nStars = len(sInds)
        if (old_sInd is None) or (nStars == 1):
            sInd = np.random.choice(sInds[comps == max(comps)])
            return sInd, slewTimes[sInd]

        # define adjacency matrix
        A = np.zeros((nStars, nStars))

        # only consider slew distance when there's an occulter
        if OS.haveOcculter:
            r_ts = TL.starprop(sInds, TK.currentTimeAbs)
            u_ts = (
                r_ts.to("AU").value.T / np.linalg.norm(r_ts.to("AU").value, axis=1)
            ).T
            angdists = np.arccos(np.clip(np.dot(u_ts, u_ts.T), -1, 1))
            A[np.ones((nStars), dtype=bool)] = angdists
            A = self.coeffs[0] * (A) / np.pi

        # add factor due to completeness
        A = A + self.coeffs[1] * (1 - comps)

        # add factor due to unvisited ramp
        f_uv = np.zeros(nStars)
        unvisited = self.starVisits[sInds] == 0
        f_uv[unvisited] = float(TK.currentTimeNorm.copy() / TK.missionLife.copy()) ** 2
        A = A - self.coeffs[2] * f_uv

        # add factor due to revisited ramp
        f2_uv = 1 - (np.in1d(sInds, self.starRevisit[:, 0]))
        A = A + self.coeffs[3] * f2_uv

        # kill diagonal
        A = A + np.diag(np.ones(nStars) * np.inf)

        # take two traversal steps
        step1 = np.tile(A[sInds == old_sInd, :], (nStars, 1)).flatten("F")
        step2 = A[np.array(np.ones((nStars, nStars)), dtype=bool)]
        tmp = np.nanargmin(step1 + step2)
        sInd = sInds[int(np.floor(tmp / float(nStars)))]

        waitTime = slewTimes[sInd]
        # Check if exoplanetObsTime would be exceeded
        mode = list(filter(lambda mode: mode["detectionMode"], allModes))[0]
        (
            maxIntTimeOBendTime,
            maxIntTimeExoplanetObsTime,
            maxIntTimeMissionLife,
        ) = TK.get_ObsDetectionMaxIntTime(Obs, mode)
        maxIntTime = min(
            maxIntTimeOBendTime, maxIntTimeExoplanetObsTime, maxIntTimeMissionLife
        )  # Maximum intTime allowed
        intTimes2 = self.calc_targ_intTime(
            np.array([sInd]), TK.currentTimeAbs.copy(), mode
        )
        if (
            intTimes2 > maxIntTime
        ):  # check if max allowed integration time would be exceeded
            self.vprint("max allowed integration time would be exceeded")
            sInd = None
            waitTime = 1.0 * u.d

        return sInd, waitTime

    def revisitFilter(self, sInds, tmpCurrentTimeNorm):
        """Helper method for Overloading Revisit Filtering

        Args:
            sInds - indices of stars still in observation list
            tmpCurrentTimeNorm (MJD) - the simulation time after overhead
            was added in MJD form

        Returns:

            sInds - indices of stars still in observation list
        """
        tovisit = np.zeros(
            self.TargetList.nStars, dtype=bool
        )  # tovisit is a boolean array containing the
        if len(sInds) > 0:  # so long as there is at least 1 star left in sInds
            tovisit[sInds] = (self.starVisits[sInds] == min(self.starVisits[sInds])) & (
                self.starVisits[sInds] < self.nVisitsMax
            )  # Checks that no star has exceeded the number of revisits
            if (
                self.starRevisit.size != 0
            ):  # There is at least one revisit planned in starRevisit
                dt_rev = (
                    self.starRevisit[:, 1] * u.day - tmpCurrentTimeNorm
                )  # absolute temporal spacing between revisit and now.

                ind_rev2 = [
                    int(x)
                    for x in self.starRevisit[dt_rev < 0 * u.d, 0]
                    if (x in sInds)
                ]
                tovisit[ind_rev2] = self.starVisits[ind_rev2] < self.nVisitsMax
            sInds = np.where(tovisit)[0]
        return sInds

    def scheduleRevisit(self, sInd, smin, det, pInds):
        """A Helper Method for scheduling revisits after observation detection
        Args:
            sInd - sInd of the star just detected
            smin - minimum separation of the planet to star of planet just detected
            det -
            pInds - Indices of planets around target star
        Return:
            updates self.starRevisit attribute
        """
        TK = self.TimeKeeping
        TL = self.TargetList
        SU = self.SimulatedUniverse
        # in both cases (detection or false alarm), schedule a revisit
        # based on minimum separation
        Ms = TL.MsTrue[sInd]
        if (
            smin is not None and smin is not np.nan
        ):  # smin is None if no planet was detected
            sp = smin
            if np.any(det):
                pInd_smin = pInds[det][np.argmin(SU.s[pInds[det]])]
                Mp = SU.Mp[pInd_smin]
            else:
                Mp = SU.Mp.mean()
            mu = const.G * (Mp + Ms)
            T = 2.0 * np.pi * np.sqrt(sp**3 / mu)
            t_rev = TK.currentTimeNorm + T / 2.0
        # otherwise, revisit based on average of population semi-major axis and mass
        else:
            sp = SU.s.mean()
            Mp = SU.Mp.mean()
            mu = const.G * (Mp + Ms)
            T = 2.0 * np.pi * np.sqrt(sp**3 / mu)
            t_rev = TK.currentTimeNorm + 0.75 * T

        t_rev = TK.currentTimeNorm.copy() + self.revisit_wait
        # finally, populate the revisit list (NOTE: sInd becomes a float)
        revisit = np.array([sInd, t_rev.to("day").value])
        if self.starRevisit.size == 0:  # If starRevisit has nothing in it
            self.starRevisit = np.array([revisit])  # initialize starRevisit
        else:
            revInd = np.where(self.starRevisit[:, 0] == sInd)[
                0
            ]  # indices of the first column of the starRevisit list containing sInd
            if revInd.size == 0:
                self.starRevisit = np.vstack((self.starRevisit, revisit))
            else:
                self.starRevisit[revInd, 1] = revisit[1]  # over

1	from EXOSIMS.SurveySimulation.linearJScheduler import linearJScheduler	1✔
2	import astropy.units as u	1✔
3	import numpy as np	1✔
4	import time	1✔
5	import astropy.constants as const	1✔
6	from EXOSIMS.util._numpy_compat import copy_if_needed	1✔
7
8
9	class linearJScheduler_det_only(linearJScheduler):	1✔
10	"""linearJScheduler_det_only - linearJScheduler Detections Only
11
12	This class implements the linear cost function scheduler described
13	in Savransky et al. (2010).
14
15	This scheduler inherits from the linearJScheduler module but performs only
16	detections.
17
18	Args:
19	specs:
20	user specified values
21
22	"""
23
24	def __init__(self, **specs):	1✔
25
26	linearJScheduler.__init__(self, **specs)	1✔
27
28	def run_sim(self):	1✔
29	"""Performs the survey simulation"""
30
31	OS = self.OpticalSystem	1✔
32	TL = self.TargetList	1✔
33	SU = self.SimulatedUniverse	1✔
34	Obs = self.Observatory	1✔
35	TK = self.TimeKeeping	1✔
36
37	# TODO: start using this self.currentSep
38	# set occulter separation if haveOcculter
39	if OS.haveOcculter:	1✔
40	self.currentSep = Obs.occulterSep	×
41
42	# choose observing modes selected for detection (default marked with a flag)
43	allModes = OS.observingModes	1✔
44	det_mode = list(filter(lambda mode: mode["detectionMode"], allModes))[0]	1✔
45
46	# begin Survey, and loop until mission is finished
47	log_begin = "OB%s: survey beginning." % (TK.OBnumber)	1✔
48	self.logger.info(log_begin)	1✔
49	self.vprint(log_begin)	1✔
50	t0 = time.time()	1✔
51	sInd = None	1✔
52	ObsNum = 0	1✔
53	while not TK.mission_is_over(OS, Obs, det_mode):	1✔
54
55	# acquire the NEXT TARGET star index and create DRM
56	old_sInd = sInd # used to save sInd if returned sInd is None	1✔
57	DRM, sInd, det_intTime, waitTime = self.next_target(sInd, det_mode)	1✔
58
59	if sInd is not None:	1✔
60	ObsNum += (	1✔
61	1 # we're making an observation so increment observation number
62	)
63
64	if OS.haveOcculter:	1✔
65	# advance to start of observation
66	# (add slew time for selected target)
67	_ = TK.advanceToAbsTime(TK.currentTimeAbs.copy() + waitTime)	×
68
69	# beginning of observation, start to populate DRM
70	DRM["star_ind"] = sInd	1✔
71	DRM["star_name"] = TL.Name[sInd]	1✔
72	DRM["arrival_time"] = TK.currentTimeNorm.to("day").copy()	1✔
73	DRM["OB_nb"] = TK.OBnumber	1✔
74	DRM["ObsNum"] = ObsNum	1✔
75	pInds = np.where(SU.plan2star == sInd)[0]	1✔
76	DRM["plan_inds"] = pInds.astype(int)	1✔
77	log_obs = (	1✔
78	" Observation #%s, star ind %s (of %s) with %s planet(s), "
79	+ "mission time at Obs start: %s"
80	) % (
81	ObsNum,
82	sInd,
83	TL.nStars,
84	len(pInds),
85	TK.currentTimeNorm.to("day").copy().round(2),
86	)
87	self.logger.info(log_obs)	1✔
88	self.vprint(log_obs)	1✔
89
90	# PERFORM DETECTION and populate revisit list attribute
91	(	1✔
92	detected,
93	det_fZ,
94	det_systemParams,
95	det_SNR,
96	FA,
97	) = self.observation_detection(sInd, det_intTime, det_mode)
98	# update the occulter wet mass
99	if OS.haveOcculter:	1✔
100	DRM = self.update_occulter_mass(DRM, sInd, det_intTime, "det")	×
101	# populate the DRM with detection results
102	DRM["det_time"] = det_intTime.to("day")	1✔
103	DRM["det_status"] = detected	1✔
104	DRM["det_SNR"] = det_SNR	1✔
105	DRM["det_fZ"] = det_fZ.to("1/arcsec2")	1✔
106	DRM["det_params"] = det_systemParams	1✔
107
108	# populate the DRM with observation modes
109	DRM["det_mode"] = dict(det_mode)	1✔
110	del DRM["det_mode"]["inst"], DRM["det_mode"]["syst"]	1✔
111
112	DRM["exoplanetObsTime"] = TK.exoplanetObsTime.copy()	1✔
113
114	# append result values to self.DRM
115	self.DRM.append(DRM)	1✔
116
117	# handle case of inf OBs and missionPortion < 1
118	if np.isinf(TK.OBduration) and (TK.missionPortion < 1):	1✔
119	self.arbitrary_time_advancement(	1✔
120	TK.currentTimeNorm.to("day").copy() - DRM["arrival_time"]
121	)
122
123	else: # sInd == None
124	sInd = old_sInd # Retain the last observed star	×
125	if (	×
126	TK.currentTimeNorm.copy() >= TK.OBendTimes[TK.OBnumber]
127	): # currentTime is at end of OB
128	# Conditional Advance To Start of Next OB
129	if not TK.mission_is_over(	×
130	OS, Obs, det_mode
131	): # as long as the mission is not over
132	TK.advancetToStartOfNextOB() # Advance To Start of Next OB	×
133	elif waitTime is not None:	×
134	# CASE 1: Advance specific wait time
135	_ = TK.advanceToAbsTime(TK.currentTimeAbs.copy() + waitTime)	×
136	self.vprint("waitTime is not None")	×
137	else:
138	startTimes = (	×
139	TK.currentTimeAbs.copy() + np.zeros(TL.nStars) * u.d
140	) # Start Times of Observations
141	observableTimes = Obs.calculate_observableTimes(	×
142	TL,
143	np.arange(TL.nStars),
144	startTimes,
145	self.koMap,
146	self.koTimes,
147	det_mode,
148	)[0]
149
150	# CASE 2 If There are no observable targets for the rest
151	# of the mission
152	if (	×
153	observableTimes[
154	(
155	TK.missionFinishAbs.copy().value * u.d
156	> observableTimes.value * u.d
157	)
158	* (
159	observableTimes.value * u.d
160	>= TK.currentTimeAbs.copy().value * u.d
161	)
162	].shape[0]
163	) == 0:
164	self.vprint(	×
165	(
166	"No Observable Targets for Remainder of mission at "
167	"currentTimeNorm = {}"
168	).format(TK.currentTimeNorm)
169	)
170	# Manually advancing time to mission end
171	TK.currentTimeNorm = TK.missionLife	×
172	TK.currentTimeAbs = TK.missionFinishAbs	×
173	else:
174	# CASE 3 nominal wait time if at least 1 target is still
175	# in list and observable
176	# TODO: ADD ADVANCE TO WHEN FZMIN OCURS
177	inds1 = np.arange(TL.nStars)[	×
178	observableTimes.value * u.d
179	> TK.currentTimeAbs.copy().value * u.d
180	]
181	# apply intTime filter
182	inds2 = np.intersect1d(self.intTimeFilterInds, inds1)	×
183	# apply revisit Filter #NOTE this means stars you added
184	# to the revisit list
185	inds3 = self.revisitFilter(	×
186	inds2, TK.currentTimeNorm.copy() + self.dt_max.to(u.d)
187	)
188	self.vprint(	×
189	"Filtering %d stars from advanceToAbsTime"
190	% (TL.nStars - len(inds3))
191	)
192	oTnowToEnd = observableTimes[inds3]	×
193	# there is at least one observableTime between now and the
194	# end of the mission
195	if not oTnowToEnd.value.shape[0] == 0:	×
196	tAbs = np.min(oTnowToEnd) # advance to that observable time	×
197	else:
198	# advance to end of mission
199	tAbs = TK.missionStart + TK.missionLife	×
200	tmpcurrentTimeNorm = TK.currentTimeNorm.copy()	×
201	# Advance Time to this time OR start of next OB following
202	# this time
203	_ = TK.advanceToAbsTime(tAbs)	×
204	self.vprint(	×
205	(
206	"No Observable Targets a currentTimeNorm= {:.2f} "
207	"Advanced To currentTimeNorm= {:.2f}"
208	).format(
209	tmpcurrentTimeNorm.to("day"),
210	TK.currentTimeNorm.to("day"),
211	)
212	)
213	else: # TK.mission_is_over()
214	dtsim = (time.time() - t0) * u.s	1✔
215	log_end = (	1✔
216	"Mission complete: no more time available.\n"
217	+ "Simulation duration: %s.\n" % dtsim.astype("int")
218	+ "Results stored in SurveySimulation.DRM (Design Reference Mission)."
219	)
220	self.logger.info(log_end)	1✔
221	print(log_end)	1✔
222
223	def next_target(self, old_sInd, mode):	1✔
224	"""Finds index of next target star and calculates its integration time.
225
226	This method chooses the next target star index based on which
227	stars are available, their integration time, and maximum completeness.
228	Returns None if no target could be found.
229
230	Args:
231	old_sInd (integer):
232	Index of the previous target star
233	mode (dict):
234	Selected observing mode for detection
235
236	Returns:
237	tuple:
238	DRM (dict):
239	Design Reference Mission, contains the results of one complete
240	observation (detection and characterization)
241	sInd (integer):
242	Index of next target star. Defaults to None.
243	intTime (astropy Quantity):
244	Selected star integration time for detection in units of day.
245	Defaults to None.
246	waitTime (astropy Quantity):
247	a strategically advantageous amount of time to wait in the case of
248	an occulter for slew times
249
250	"""
251	OS = self.OpticalSystem	1✔
252	TL = self.TargetList	1✔
253	Obs = self.Observatory	1✔
254	TK = self.TimeKeeping	1✔
255
256	# create DRM
257	DRM = {}	1✔
258
259	# selecting appropriate koMap
260	koMap = self.koMaps[mode["syst"]["name"]]	1✔
261
262	# allocate settling time + overhead time
263	tmpCurrentTimeAbs = (	1✔
264	TK.currentTimeAbs.copy() + Obs.settlingTime + mode["syst"]["ohTime"]
265	)
266	tmpCurrentTimeNorm = (	1✔
267	TK.currentTimeNorm.copy() + Obs.settlingTime + mode["syst"]["ohTime"]
268	)
269
270	# look for available targets
271	# 1. initialize arrays
272	slewTimes = np.zeros(TL.nStars) * u.d	1✔
273	# fZs = np.zeros(TL.nStars) / u.arcsec**2
274	dV = np.zeros(TL.nStars) * u.m / u.s	1✔
275	intTimes = np.zeros(TL.nStars) * u.d	1✔
276	obsTimes = np.zeros([2, TL.nStars]) * u.d	1✔
277	sInds = np.arange(TL.nStars)	1✔
278
279	# 2. find spacecraft orbital START positions (if occulter, positions
280	# differ for each star) and filter out unavailable targets
281	sd = None	1✔
282	if OS.haveOcculter:	1✔
283	sd = Obs.star_angularSep(TL, old_sInd, sInds, tmpCurrentTimeAbs)	×
284	obsTimes = Obs.calculate_observableTimes(	×
285	TL, sInds, tmpCurrentTimeAbs, self.koMaps, self.koTimes, mode
286	)
287	slewTimes = Obs.calculate_slewTimes(	×
288	TL, old_sInd, sInds, sd, obsTimes, tmpCurrentTimeAbs
289	)
290
291	# 2.1 filter out totTimes > integration cutoff
292	if len(sInds.tolist()) > 0:	1✔
293	sInds = np.intersect1d(self.intTimeFilterInds, sInds)	1✔
294
295	# start times, including slew times
296	startTimes = tmpCurrentTimeAbs.copy() + slewTimes	1✔
297	startTimesNorm = tmpCurrentTimeNorm.copy() + slewTimes	1✔
298
299	# 2.5 Filter stars not observable at startTimes
300	try:	1✔
301	koTimeInd = np.where(	1✔
302	np.round(startTimes[0].value) - self.koTimes.value == 0
303	)[0][
304	0
305	] # find indice where koTime is startTime[0]
306	sInds = sInds[	1✔
307	np.where(np.transpose(koMap)[koTimeInd].astype(bool)[sInds])[0]
308	] # filters inds by koMap #verified against v1.35
309	except: # noqa: E722 If there are no target stars to observe	×
310	sInds = np.asarray([], dtype=int)	×
311
312	# 3. filter out all previously (more-)visited targets, unless in
313	if len(sInds.tolist()) > 0:	1✔
314	sInds = self.revisitFilter(sInds, tmpCurrentTimeNorm)	1✔
315
316	# 4.1 calculate integration times for ALL preselected targets
317	(	1✔
318	maxIntTimeOBendTime,
319	maxIntTimeExoplanetObsTime,
320	maxIntTimeMissionLife,
321	) = TK.get_ObsDetectionMaxIntTime(Obs, mode)
322	maxIntTime = min(	1✔
323	maxIntTimeOBendTime, maxIntTimeExoplanetObsTime, maxIntTimeMissionLife
324	) # Maximum intTime allowed
325
326	if len(sInds.tolist()) > 0:	1✔
327	intTimes[sInds] = self.calc_targ_intTime(sInds, startTimes[sInds], mode)	1✔
328	sInds = sInds[	1✔
329	np.where(intTimes[sInds] <= maxIntTime)
330	] # Filters targets exceeding end of OB
331	endTimes = startTimes + intTimes	1✔
332
333	if maxIntTime.value <= 0:	1✔
334	sInds = np.asarray([], dtype=int)	×
335
336	# 5.1 TODO Add filter to filter out stars entering and exiting keepout
337	# between startTimes and endTimes
338
339	# 5.2 find spacecraft orbital END positions (for each candidate target),
340	# and filter out unavailable targets
341	if len(sInds.tolist()) > 0 and Obs.checkKeepoutEnd:	1✔
342	try:	1✔
343	# endTimes may exist past koTimes so we have an exception to
344	# handle this case
345	# koTimeInd[0][0] # find indice where koTime is endTime[0]
346	koTimeInd = np.where(	1✔
347	np.round(endTimes[0].value) - self.koTimes.value == 0
348	)[0][0]
349	# filters inds by koMap #verified against v1.35
350	sInds = sInds[	1✔
351	np.where(np.transpose(koMap)[koTimeInd].astype(bool)[sInds])[0]
352	]
353	except: # noqa: E722	×
354	sInds = np.asarray([], dtype=int)	×
355
356	# 6. choose best target from remaining
357	if len(sInds.tolist()) > 0:	1✔
358	# choose sInd of next target
359	sInd, waitTime = self.choose_next_target(	1✔
360	old_sInd, sInds, slewTimes, intTimes[sInds]
361	)
362
363	# Should Choose Next Target decide there are no stars it wishes to
364	# observe at this time.
365	if (sInd is None) and (waitTime is not None):	1✔
366	self.vprint(	×
367	(
368	"There are no stars Choose Next Target would like to Observe. "
369	"Waiting {}"
370	).format(waitTime)
371	)
372	return DRM, None, None, waitTime	×
373	elif (sInd is None) and (waitTime is not None):	1✔
374	self.vprint(	×
375	(
376	"There are no stars Choose Next Target would like to Observe "
377	"and waitTime is None"
378	)
379	)
380	return DRM, None, None, waitTime	×
381	# store selected star integration time
382	intTime = intTimes[sInd]	1✔
383
384	# if no observable target, advanceTime to next Observable Target
385	else:
386	self.vprint(	×
387	"No Observable Targets at currentTimeNorm= "
388	+ str(TK.currentTimeNorm.copy())
389	)
390	return DRM, None, None, None	×
391
392	# update visited list for selected star
393	self.starVisits[sInd] += 1	1✔
394	# store normalized start time for future completeness update
395	self.lastObsTimes[sInd] = startTimesNorm[sInd]	1✔
396
397	# populate DRM with occulter related values
398	if OS.haveOcculter:	1✔
399	DRM = Obs.log_occulterResults(	×
400	DRM, slewTimes[sInd], sInd, sd[sInd], dV[sInd]
401	)
402	return DRM, sInd, intTime, slewTimes[sInd]	×
403
404	return DRM, sInd, intTime, waitTime	1✔
405
406	def choose_next_target(self, old_sInd, sInds, slewTimes, intTimes):	1✔
407	"""Choose next target based on truncated depth first search
408	of linear cost function.
409
410	Args:
411	old_sInd (integer):
412	Index of the previous target star
413	sInds (integer array):
414	Indices of available targets
415	slewTimes (astropy quantity array):
416	slew times to all stars (must be indexed by sInds)
417	intTimes (astropy Quantity array):
418	Integration times for detection in units of day
419
420	Returns:
421	sInd (integer):
422	Index of next target star
423
424	"""
425
426	Comp = self.Completeness	1✔
427	TL = self.TargetList	1✔
428	TK = self.TimeKeeping	1✔
429	OS = self.OpticalSystem	1✔
430	Obs = self.Observatory	1✔
431	allModes = OS.observingModes	1✔
432
433	# cast sInds to array
434	sInds = np.array(sInds, ndmin=1, copy=copy_if_needed)	1✔
435
436	# current star has to be in the adjmat
437	if (old_sInd is not None) and (old_sInd not in sInds):	1✔
438	sInds = np.append(sInds, old_sInd)	1✔
439
440	# calculate dt since previous observation
441	dt = TK.currentTimeNorm + slewTimes[sInds] - self.lastObsTimes[sInds]	1✔
442	# get dynamic completeness values
443	comps = Comp.completeness_update(TL, sInds, self.starVisits[sInds], dt)	1✔
444
445	# if first target, or if only 1 available target,
446	# choose highest available completeness
447	nStars = len(sInds)	1✔
448	if (old_sInd is None) or (nStars == 1):	1✔
449	sInd = np.random.choice(sInds[comps == max(comps)])	1✔
450	return sInd, slewTimes[sInd]	1✔
451
452	# define adjacency matrix
453	A = np.zeros((nStars, nStars))	1✔
454
455	# only consider slew distance when there's an occulter
456	if OS.haveOcculter:	1✔
457	r_ts = TL.starprop(sInds, TK.currentTimeAbs)	1✔
458	u_ts = (	1✔
459	r_ts.to("AU").value.T / np.linalg.norm(r_ts.to("AU").value, axis=1)
460	).T
461	angdists = np.arccos(np.clip(np.dot(u_ts, u_ts.T), -1, 1))	1✔
462	A[np.ones((nStars), dtype=bool)] = angdists	1✔
463	A = self.coeffs[0] * (A) / np.pi	1✔
464
465	# add factor due to completeness
466	A = A + self.coeffs[1] * (1 - comps)	1✔
467
468	# add factor due to unvisited ramp
469	f_uv = np.zeros(nStars)	1✔
470	unvisited = self.starVisits[sInds] == 0	1✔
471	f_uv[unvisited] = float(TK.currentTimeNorm.copy() / TK.missionLife.copy()) ** 2	1✔
472	A = A - self.coeffs[2] * f_uv	1✔
473
474	# add factor due to revisited ramp
475	f2_uv = 1 - (np.in1d(sInds, self.starRevisit[:, 0]))	1✔
476	A = A + self.coeffs[3] * f2_uv	1✔
477
478	# kill diagonal
479	A = A + np.diag(np.ones(nStars) * np.inf)	1✔
480
481	# take two traversal steps
482	step1 = np.tile(A[sInds == old_sInd, :], (nStars, 1)).flatten("F")	1✔
483	step2 = A[np.array(np.ones((nStars, nStars)), dtype=bool)]	1✔
484	tmp = np.nanargmin(step1 + step2)	1✔
485	sInd = sInds[int(np.floor(tmp / float(nStars)))]	1✔
486
487	waitTime = slewTimes[sInd]	1✔
488	# Check if exoplanetObsTime would be exceeded
489	mode = list(filter(lambda mode: mode["detectionMode"], allModes))[0]	1✔
490	(	1✔
491	maxIntTimeOBendTime,
492	maxIntTimeExoplanetObsTime,
493	maxIntTimeMissionLife,
494	) = TK.get_ObsDetectionMaxIntTime(Obs, mode)
495	maxIntTime = min(	1✔
496	maxIntTimeOBendTime, maxIntTimeExoplanetObsTime, maxIntTimeMissionLife
497	) # Maximum intTime allowed
498	intTimes2 = self.calc_targ_intTime(	1✔
499	np.array([sInd]), TK.currentTimeAbs.copy(), mode
500	)
501	if (	1✔
502	intTimes2 > maxIntTime
503	): # check if max allowed integration time would be exceeded
504	self.vprint("max allowed integration time would be exceeded")	×
505	sInd = None	×
506	waitTime = 1.0 * u.d	×
507
508	return sInd, waitTime	1✔
509
510	def revisitFilter(self, sInds, tmpCurrentTimeNorm):	1✔
511	"""Helper method for Overloading Revisit Filtering
512
513	Args:
514	sInds - indices of stars still in observation list
515	tmpCurrentTimeNorm (MJD) - the simulation time after overhead
516	was added in MJD form
517
518	Returns:
519
520	sInds - indices of stars still in observation list
521	"""
522	tovisit = np.zeros(	1✔
523	self.TargetList.nStars, dtype=bool
524	) # tovisit is a boolean array containing the
525	if len(sInds) > 0: # so long as there is at least 1 star left in sInds	1✔
526	tovisit[sInds] = (self.starVisits[sInds] == min(self.starVisits[sInds])) & (	1✔
527	self.starVisits[sInds] < self.nVisitsMax
528	) # Checks that no star has exceeded the number of revisits
529	if (	1✔
530	self.starRevisit.size != 0
531	): # There is at least one revisit planned in starRevisit
532	dt_rev = (	1✔
533	self.starRevisit[:, 1] * u.day - tmpCurrentTimeNorm
534	) # absolute temporal spacing between revisit and now.
535
536	ind_rev2 = [	1✔
537	int(x)
538	for x in self.starRevisit[dt_rev < 0 * u.d, 0]
539	if (x in sInds)
540	]
541	tovisit[ind_rev2] = self.starVisits[ind_rev2] < self.nVisitsMax	1✔
542	sInds = np.where(tovisit)[0]	1✔
543	return sInds	1✔
544
545	def scheduleRevisit(self, sInd, smin, det, pInds):	1✔
546	"""A Helper Method for scheduling revisits after observation detection
547	Args:
548	sInd - sInd of the star just detected
549	smin - minimum separation of the planet to star of planet just detected
550	det -
551	pInds - Indices of planets around target star
552	Return:
553	updates self.starRevisit attribute
554	"""
555	TK = self.TimeKeeping	1✔
556	TL = self.TargetList	1✔
557	SU = self.SimulatedUniverse	1✔
558	# in both cases (detection or false alarm), schedule a revisit
559	# based on minimum separation
560	Ms = TL.MsTrue[sInd]	1✔
561	if (	1✔
562	smin is not None and smin is not np.nan
563	): # smin is None if no planet was detected
UNCOV 564	sp = smin	×
UNCOV 565	if np.any(det):	×
UNCOV 566	pInd_smin = pInds[det][np.argmin(SU.s[pInds[det]])]	×
UNCOV 567	Mp = SU.Mp[pInd_smin]	×
568	else:
569	Mp = SU.Mp.mean()	×
UNCOV 570	mu = const.G * (Mp + Ms)	×
UNCOV 571	T = 2.0 * np.pi * np.sqrt(sp**3 / mu)	×
UNCOV 572	t_rev = TK.currentTimeNorm + T / 2.0	×
573	# otherwise, revisit based on average of population semi-major axis and mass
574	else:
575	sp = SU.s.mean()	1✔
576	Mp = SU.Mp.mean()	1✔
577	mu = const.G * (Mp + Ms)	1✔
578	T = 2.0 * np.pi * np.sqrt(sp**3 / mu)	1✔
579	t_rev = TK.currentTimeNorm + 0.75 * T	1✔
580
581	t_rev = TK.currentTimeNorm.copy() + self.revisit_wait	1✔
582	# finally, populate the revisit list (NOTE: sInd becomes a float)
583	revisit = np.array([sInd, t_rev.to("day").value])	1✔
584	if self.starRevisit.size == 0: # If starRevisit has nothing in it	1✔
585	self.starRevisit = np.array([revisit]) # initialize starRevisit	1✔
586	else:
587	revInd = np.where(self.starRevisit[:, 0] == sInd)[	1✔
588	0
589	] # indices of the first column of the starRevisit list containing sInd
590	if revInd.size == 0:	1✔
591	self.starRevisit = np.vstack((self.starRevisit, revisit))	1✔
592	else:
593	self.starRevisit[revInd, 1] = revisit[1] # over	×

dsavransky / EXOSIMS / 11979115482

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