Prepare arguments and call the RRTGM_SW driver to calculate
        radiative fluxes and heating rates
191668
def _compute_heating_rates(self):
        ''''''
        (ncol, nlay, icld, iaer, permuteseed, irng,
         play, plev, tlay, tlev, tsfc,
         h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
         aldif, aldir, asdif, asdir, coszen, adjes, dyofyr, scon, isolvar,
         indsolvar, bndsolvar, solcycfrac,
         inflgsw, iceflgsw, liqflgsw,
         cldfrac, ciwp, clwp, reic, relq, tauc, ssac, asmc, fsfc,
         tauaer, ssaaer, asmaer, ecaer,) = self._prepare_sw_arguments()
        if icld == 0:  # clear-sky only
            cldfmcl = np.zeros((ngptsw,ncol,nlay))
            ciwpmcl = np.zeros((ngptsw,ncol,nlay))
            clwpmcl = np.zeros((ngptsw,ncol,nlay))
            reicmcl = np.zeros((ncol,nlay))
            relqmcl = np.zeros((ncol,nlay))
            taucmcl = np.zeros((ngptsw,ncol,nlay))
            ssacmcl = np.zeros((ngptsw,ncol,nlay))
            asmcmcl = np.zeros((ngptsw,ncol,nlay))
            fsfcmcl = np.zeros((ngptsw,ncol,nlay))
        else:
            #  Call the Monte Carlo Independent Column Approximation (McICA, Pincus et al., JC, 2003)
            (cldfmcl, ciwpmcl, clwpmcl, reicmcl, relqmcl, taucmcl,
            ssacmcl, asmcmcl, fsfcmcl) = _rrtmg_sw.climlab_mcica_subcol_sw(
                            ncol, nlay, icld, permuteseed, irng, play,
                            cldfrac, ciwp, clwp, reic, relq, tauc, ssac, asmc, fsfc)
        #  Call the RRTMG_SW driver to compute radiative fluxes
        (swuflx, swdflx, swhr, swuflxc, swdflxc, swhrc) = \
            _rrtmg_sw.climlab_rrtmg_sw(ncol, nlay, icld, iaer,
                play, plev, tlay, tlev, tsfc,
                h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr,
                asdir, asdif, aldir, aldif,
                coszen, adjes, dyofyr, scon, isolvar,
                inflgsw, iceflgsw, liqflgsw, cldfmcl,
                taucmcl, ssacmcl, asmcmcl, fsfcmcl,
                ciwpmcl, clwpmcl, reicmcl, relqmcl,
                tauaer, ssaaer, asmaer, ecaer,
                bndsolvar, indsolvar, solcycfrac)
        #  Output is all (ncol,nlay+1) or (ncol,nlay)
        self.SW_flux_up = _rrtm_to_climlab(swuflx) + 0.*self.SW_flux_up
        self.SW_flux_down = _rrtm_to_climlab(swdflx) + 0.*self.SW_flux_down
        self.SW_flux_up_clr = _rrtm_to_climlab(swuflxc) + 0.*self.SW_flux_up_clr
        self.SW_flux_down_clr = _rrtm_to_climlab(swdflxc) + 0.*self.SW_flux_down_clr
        #  Compute quantities derived from fluxes, including ASR
        self._compute_SW_flux_diagnostics()
        #  calculate heating rates from flux divergence
        SWheating_Wm2 = np.array(-np.diff(self.SW_flux_net, axis=-1)) + 0.*self.Tatm
        SWheating_clr_Wm2 = np.array(-np.diff(self.SW_flux_net_clr, axis=-1)) + 0.*self.Tatm
        self.heating_rate['Ts'] = np.array(self.SW_flux_net[..., -1, np.newaxis]) + 0.*self.Ts
        self.heating_rate['Tatm'] = SWheating_Wm2
        #  Convert to K / day
        Catm = self.Tatm.domain.heat_capacity
        self.TdotSW = SWheating_Wm2 / Catm * const.seconds_per_day
        self.TdotSW_clr = SWheating_clr_Wm2 / Catm * const.seconds_per_day
1290567
def run(self):
        '''
        

        '''

        try:

            # Load raw data
            log.info("Loading target data...")
            self.load_tpf()
            self.mask_planets()
            self.plot_aperture([self.dvs.top_right() for i in range(4)])
            self.init_kernel()
            M = self.apply_mask(np.arange(len(self.time)))
            self.cdppr_arr = self.get_cdpp_arr()
            self.cdpp_arr = np.array(self.cdppr_arr)
            self.cdppv_arr = np.array(self.cdppr_arr)
            self.cdppr = self.get_cdpp()
            self.cdpp = self.cdppr
            self.cdppv = self.cdppr

            log.info("%s (Raw): CDPP = %s" % (self.name, self.cdpps))
            self.plot_lc(self.dvs.left(), info_right='Raw', color='k')

            # Loop
            for n in range(self.pld_order):
                self.lam_idx += 1
                self.get_outliers()
                if n > 0 and self.optimize_gp:
                    self.update_gp()
                self.cross_validate(self.dvs.right(), info='CV%d' % n)
                self.cdpp_arr = self.get_cdpp_arr()
                self.cdppv_arr *= self.cdpp_arr
                self.cdpp = self.get_cdpp()
                self.cdppv = np.nanmean(self.cdppv_arr)
                log.info("%s (%d/%d): CDPP = %s" %
                         (self.name, n + 1, self.pld_order, self.cdpps))
                self.plot_lc(self.dvs.left(), info_right='LC%d' % (
                    n + 1), info_left='%d outliers' % len(self.outmask))

            # Save
            self.finalize()
            self.plot_final(self.dvs.top_left())
            self.plot_info(self.dvs)
            self.save_model()

        except:

            self.exception_handler(self.debug)
1370692
def rungtd1d(time: Union[datetime, str, np.ndarray],
             altkm: np.ndarray,
             glat: float, glon: float) -> xarray.Dataset:
    
    time = todt64(time)
    # %% get solar parameters for date
    f107Ap = gi.getApF107(time, smoothdays=81)
    f107a = f107Ap['f107s'].item()
    f107 = f107Ap['f107'].item()
    Ap = f107Ap['Ap'].item()
# %% dimensions
    altkm = np.atleast_1d(altkm)
    assert altkm.ndim == 1
    assert isinstance(glon, (int, float))
    assert isinstance(glat, (int, float))
    assert isinstance(time, np.datetime64) or (time.size == 1 and isinstance(
        time[0], np.datetime64)), 'if you have multiple times, for loop over them'

# don't check ap, too complicated
    assert isinstance(MASS, (float, int))
    assert len(TSELECOPS) == 25
# %%
    gtd7.tselec(TSELECOPS)  # like the msis_driver example

    iyd, utsec, stl = datetime2gtd(time, glon)
    altkm = np.atleast_1d(altkm)

    gtd7.meters(1)  # makes output in m^-3 and kg/m^-3
# %%
    if isinstance(Ap, (float, int)):
        Ap = [Ap]*7  # even if SW(9) == 1 due to f2py needs for array

    dens = np.empty((altkm.size, len(species)))
    temp = np.empty((altkm.size, len(ttypes)))
    for i, a in enumerate(altkm):
        dens[i, :], temp[i, :] = gtd7.gtd7(iyd, utsec, a, glat, glon, stl, f107a, f107, Ap, MASS)

    dsf = {k: (('time', 'alt_km', 'lat', 'lon'), v[None, :, None, None]) for (k, v) in zip(species, dens.T)}
    dsf.update({'Tn':  (('time', 'alt_km', 'lat', 'lon'), temp[:, 1][None, :, None, None]),
                'Texo': (('time', 'alt_km', 'lat', 'lon'), temp[:, 0][None, :, None, None])})

    atmos = xarray.Dataset(dsf,
                           coords={'time': time.astype(datetime), 'alt_km': altkm, 'lat': [glat], 'lon': [glon], },
                           attrs={'Ap': Ap, 'f107': f107, 'f107a': f107a,
                                  'species': species})

    return atmos
1160267
def build_hstwcs(crval1, crval2, crpix1, crpix2, naxis1, naxis2, pscale, orientat):
    
    wcsout = wcsutil.HSTWCS()
    wcsout.wcs.crval = np.array([crval1, crval2])
    wcsout.wcs.crpix = np.array([crpix1, crpix2])
    wcsout.naxis1 = naxis1
    wcsout.naxis2 = naxis2
    wcsout.wcs.cd = buildRotMatrix(orientat) * [-1, 1] * pscale / 3600.0
    # Synchronize updates with astropy.wcs objects
    wcsout.wcs.set()
    wcsout.setPscale()
    wcsout.setOrient()
    wcsout.wcs.ctype = ['RA---TAN', 'DEC--TAN']

    return wcsout
1294563
def prepare_stateseries(self, ramflag: bool = True) -> None:
        """"""
        for element in printtools.progressbar(self):
            element.prepare_stateseries(ramflag)
1233466
def main(arguments=None):
    
    # setup the command-line util settings
    su = tools(
        arguments=arguments,
        docString=__doc__,
        logLevel="WARNING",
        options_first=False,
        projectName="qubits"
    )
    arguments, settings, log, dbConn = su.setup()

    # unpack remaining cl arguments using `exec` to setup the variable names
    # automatically
    for arg, val in arguments.iteritems():
        if arg[0] == "-":
            varname = arg.replace("-", "") + "Flag"
        else:
            varname = arg.replace("<", "").replace(">", "")
        if varname == "import":
            varname = "iimport"
        if isinstance(val, str) or isinstance(val, unicode):
            exec(varname + " = '%s'" % (val,))
        else:
            exec(varname + " = %s" % (val,))
        if arg == "--dbConn":
            dbConn = val
        log.debug('%s = %s' % (varname, val,))

    ## START LOGGING ##
    startTime = times.get_now_sql_datetime()
    log.info(
        '--- STARTING TO RUN THE cl_utils.py AT %s' %
        (startTime,))

    if init:
        from . import workspace
        ws = workspace(
            log=log,
            pathToWorkspace=pathToWorkspace
        )
        ws.setup()
        return

    # IMPORT THE SIMULATION SETTINGS
    (allSettings,
     programSettings,
     limitingMags,
     sampleNumber,
     peakMagnitudeDistributions,
     explosionDaysFromSettings,
     extendLightCurveTail,
     relativeSNRates,
     lowerRedshiftLimit,
     upperRedshiftLimit,
     redshiftResolution,
     restFrameFilter,
     kCorrectionTemporalResolution,
     kCorPolyOrder,
     kCorMinimumDataPoints,
     extinctionType,
     extinctionConstant,
     hostExtinctionDistributions,
     galacticExtinctionDistribution,
     surveyCadenceSettings,
     snLightCurves,
     surveyArea,
     CCSNRateFraction,
     transientToCCSNRateFraction,
     extraSurveyConstraints,
     lightCurvePolyOrder,
     logLevel) = cu.read_in_survey_parameters(
        log,
        pathToSettingsFile=pathToSettingsFile
    )

    logFilePath = pathToOutputDirectory + "/qubits.log"
    del log
    log = _set_up_command_line_tool(
        level=str(logLevel),
        logFilePath=logFilePath
    )

    # dbConn, log = cu.settings(
    #     pathToSettingsFile=pathToSettingsFile,
    #     dbConn=False,
    #     log=True
    # )

    ## START LOGGING ##
    startTime = dcu.get_now_sql_datetime()
    log.info('--- STARTING TO RUN THE qubits AT %s' % (startTime,))

    resultsDict = {}

    pathToOutputPlotDirectory = pathToOutputDirectory + "/plots/"
    dcu.dryx_mkdir(
        log,
        directoryPath=pathToOutputPlotDirectory
    )

    pathToResultsFolder = pathToOutputDirectory + "/results/"
    dcu.dryx_mkdir(
        log,
        directoryPath=pathToResultsFolder
    )

    if not programSettings['Extract Lightcurves from Spectra'] and not programSettings['Generate KCorrection Database'] and not programSettings['Run the Simulation'] and not programSettings['Compile and Plot Results']:
        print "All stages of the simulatation have been switched off. Please switch on at least one stage of the simulation under the 'Programming Settings' in the settings file `%(pathToSettingsFile)s`" % locals()

    # GENERATE THE DATA FOR SIMULATIONS
    if programSettings['Extract Lightcurves from Spectra']:
        log.info('generating the Lightcurves')
        dg.generate_model_lightcurves(
            log=log,
            pathToSpectralDatabase=pathToSpectralDatabase,
            pathToOutputDirectory=pathToOutputDirectory,
            pathToOutputPlotDirectory=pathToOutputPlotDirectory,
            explosionDaysFromSettings=explosionDaysFromSettings,
            extendLightCurveTail=extendLightCurveTail,
            polyOrder=lightCurvePolyOrder
        )
        print "The lightcurve file can be found here: %(pathToOutputDirectory)stransient_light_curves.yaml" % locals()
        print "The lightcurve plots can be found in %(pathToOutputPlotDirectory)s" % locals()

    if programSettings['Generate KCorrection Database']:
        log.info('generating the kcorrection data')
        dg.generate_kcorrection_listing_database(
            log,
            pathToOutputDirectory=pathToOutputDirectory,
            pathToSpectralDatabase=pathToSpectralDatabase,
            restFrameFilter=restFrameFilter,
            temporalResolution=kCorrectionTemporalResolution,
            redshiftResolution=redshiftResolution,
            redshiftLower=lowerRedshiftLimit,
            redshiftUpper=upperRedshiftLimit + redshiftResolution)
        log.info('generating the kcorrection polynomials')
        dg.generate_kcorrection_polynomial_database(
            log,
            pathToOutputDirectory=pathToOutputDirectory,
            restFrameFilter=restFrameFilter,
            kCorPolyOrder=kCorPolyOrder,  # ORDER OF THE POLYNOMIAL TO FIT
            kCorMinimumDataPoints=kCorMinimumDataPoints,
            redshiftResolution=redshiftResolution,
            redshiftLower=lowerRedshiftLimit,
            redshiftUpper=upperRedshiftLimit + redshiftResolution,
            plot=programSettings['Generate KCorrection Plots'])

        print "The k-correction database has been generated here: %(pathToOutputDirectory)sk_corrections" % locals()
        if programSettings['Generate KCorrection Plots']:
            print "The k-correction polynomial plots can also be found in %(pathToOutputDirectory)sk_corrections" % locals()

    if programSettings['Run the Simulation']:
        # CREATE THE OBSERVABLE UNIVERSE!
        log.info('generating the redshift array')
        redshiftArray = u.random_redshift_array(
            log,
            sampleNumber,
            lowerRedshiftLimit,
            upperRedshiftLimit,
            redshiftResolution=redshiftResolution,
            pathToOutputPlotDirectory=pathToOutputPlotDirectory,
            plot=programSettings['Plot Simulation Helper Plots'])
        resultsDict['Redshifts'] = redshiftArray.tolist()

        log.info('generating the SN type array')
        snTypesArray = u.random_sn_types_array(
            log,
            sampleNumber,
            relativeSNRates,
            pathToOutputPlotDirectory=pathToOutputPlotDirectory,
            plot=programSettings['Plot Simulation Helper Plots'])
        resultsDict['SN Types'] = snTypesArray.tolist()

        log.info('generating peak magnitudes for the SNe')
        peakMagnitudesArray = u.random_peak_magnitudes(
            log,
            peakMagnitudeDistributions,
            snTypesArray,
            plot=programSettings['Plot Simulation Helper Plots'])

        log.info('generating the SN host extictions array')
        hostExtinctionArray = u.random_host_extinction(
            log,
            sampleNumber,
            extinctionType,
            extinctionConstant,
            hostExtinctionDistributions,
            plot=programSettings['Plot Simulation Helper Plots'])

        log.info('generating the SN galactic extictions array')
        galacticExtinctionArray = u.random_galactic_extinction(
            log,
            sampleNumber,
            extinctionType,
            extinctionConstant,
            galacticExtinctionDistribution,
            plot=programSettings['Plot Simulation Helper Plots'])

        log.info('generating the raw lightcurves for the SNe')
        rawLightCurveDict = u.generate_numpy_polynomial_lightcurves(
            log,
            snLightCurves=snLightCurves,
            pathToOutputDirectory=pathToOutputDirectory,
            pathToOutputPlotDirectory=pathToOutputPlotDirectory,
            plot=programSettings['Plot Simulation Helper Plots'])

        log.info('generating the k-correction array for the SNe')
        kCorrectionArray = u.build_kcorrection_array(
            log,
            redshiftArray,
            snTypesArray,
            snLightCurves,
            pathToOutputDirectory=pathToOutputDirectory,
            plot=programSettings['Plot Simulation Helper Plots'])

        log.info('generating the observed lightcurves for the SNe')
        observedFrameLightCurveInfo, peakAppMagList = u.convert_lightcurves_to_observered_frame(
            log,
            snLightCurves=snLightCurves,
            rawLightCurveDict=rawLightCurveDict,
            redshiftArray=redshiftArray,
            snTypesArray=snTypesArray,
            peakMagnitudesArray=peakMagnitudesArray,
            kCorrectionArray=kCorrectionArray,
            hostExtinctionArray=hostExtinctionArray,
            galacticExtinctionArray=galacticExtinctionArray,
            restFrameFilter=restFrameFilter,
            pathToOutputDirectory=pathToOutputDirectory,
            pathToOutputPlotDirectory=pathToOutputPlotDirectory,
            polyOrder=lightCurvePolyOrder,
            plot=programSettings['Plot Simulation Helper Plots'])

        log.info('generating the survey observation cadence')
        cadenceDictionary = ss.survey_cadence_arrays(
            log,
            surveyCadenceSettings,
            pathToOutputDirectory=pathToOutputDirectory,
            pathToOutputPlotDirectory=pathToOutputPlotDirectory,
            plot=programSettings['Plot Simulation Helper Plots'])

        log.info('determining if the SNe are discoverable by the survey')
        discoverableList = ss.determine_if_sne_are_discoverable(
            log,
            redshiftArray=redshiftArray,
            limitingMags=limitingMags,
            observedFrameLightCurveInfo=observedFrameLightCurveInfo,
            pathToOutputDirectory=pathToOutputDirectory,
            pathToOutputPlotDirectory=pathToOutputPlotDirectory,
            plot=programSettings['Plot Simulation Helper Plots'])

        log.info(
            'determining the day (if and) when each SN is first discoverable by the survey')
        ripeDayList = ss.determine_when_sne_are_ripe_for_discovery(
            log,
            redshiftArray=redshiftArray,
            limitingMags=limitingMags,
            discoverableList=discoverableList,
            observedFrameLightCurveInfo=observedFrameLightCurveInfo,
            plot=programSettings['Plot Simulation Helper Plots'])

        # log.info('determining the day when each SN is disappears fainter than the survey limiting mags')
        # disappearDayList = determine_when_discovered_sne_disappear(
        #     log,
        #     redshiftArray=redshiftArray,
        #     limitingMags=limitingMags,
        #     ripeDayList=ripeDayList,
        #     observedFrameLightCurveInfo=observedFrameLightCurveInfo,
        #     plot=programSettings['Plot Simulation Helper Plots'])

        log.info('determining if and when each SN is discovered by the survey')
        lightCurveDiscoveryDayList, surveyDiscoveryDayList, snCampaignLengthList = ss.determine_if_sne_are_discovered(
            log,
            limitingMags=limitingMags,
            ripeDayList=ripeDayList,
            cadenceDictionary=cadenceDictionary,
            observedFrameLightCurveInfo=observedFrameLightCurveInfo,
            extraSurveyConstraints=extraSurveyConstraints,
            plot=programSettings['Plot Simulation Helper Plots'])

        resultsDict[
            'Discoveries Relative to Peak Magnitudes'] = lightCurveDiscoveryDayList
        resultsDict[
            'Discoveries Relative to Survey Year'] = surveyDiscoveryDayList
        resultsDict['Campaign Length'] = snCampaignLengthList
        resultsDict['Cadence Dictionary'] = cadenceDictionary
        resultsDict['Peak Apparent Magnitudes'] = peakAppMagList

        now = datetime.now()
        now = now.strftime("%Y%m%dt%H%M%S")
        fileName = pathToOutputDirectory + \
            "/simulation_results_%s.yaml" % (now,)
        stream = file(fileName, 'w')
        yamlContent = dict(allSettings.items() + resultsDict.items())
        yaml.dump(yamlContent, stream, default_flow_style=False)
        stream.close()

        print "The simulation output file can be found here: %(fileName)s. Remember to update your settings file 'Simulation Results File Used for Plots' parameter with this filename before compiling the results." % locals()
        if programSettings['Plot Simulation Helper Plots']:
            print "The simulation helper-plots found in %(pathToOutputPlotDirectory)s" % locals()

    # COMPILE AND PLOT THE RESULTS
    if programSettings['Compile and Plot Results']:
        pathToYamlFile = pathToOutputDirectory + "/" + \
            programSettings['Simulation Results File Used for Plots']
        result_log = r.log_the_survey_settings(log, pathToYamlFile)
        snSurveyDiscoveryTimes, lightCurveDiscoveryTimes, snTypes, redshifts, cadenceDictionary, peakAppMagList, snCampaignLengthList = r.import_results(
            log, pathToYamlFile)
        snRatePlotLink, totalRate, tooFaintRate, shortCampaignRate = r.determine_sn_rate(
            log,
            lightCurveDiscoveryTimes,
            snSurveyDiscoveryTimes,
            redshifts,
            surveyCadenceSettings=surveyCadenceSettings,
            lowerRedshiftLimit=lowerRedshiftLimit,
            upperRedshiftLimit=upperRedshiftLimit,
            redshiftResolution=redshiftResolution,
            surveyArea=surveyArea,
            CCSNRateFraction=CCSNRateFraction,
            transientToCCSNRateFraction=transientToCCSNRateFraction,
            peakAppMagList=peakAppMagList,
            snCampaignLengthList=snCampaignLengthList,
            extraSurveyConstraints=extraSurveyConstraints,
            pathToOutputPlotFolder=pathToOutputPlotDirectory)
        result_log += """
## Results ##

This simulated survey discovered a total of **%s** transients per year. An extra **%s** transients were detected but deemed too faint to constrain a positive transient identification and a further **%s** transients where detected but an observational campaign of more than **%s** days could not be completed to ensure identification. See below for the various output plots.

        """ % (totalRate, tooFaintRate, shortCampaignRate, extraSurveyConstraints["Observable for at least ? number of days"])
        cadenceWheelLink = r.plot_cadence_wheel(
            log,
            cadenceDictionary,
            pathToOutputPlotFolder=pathToOutputPlotDirectory)
        result_log += """%s""" % (cadenceWheelLink,)
        discoveryMapLink = r.plot_sn_discovery_map(
            log,
            snSurveyDiscoveryTimes,
            peakAppMagList,
            snCampaignLengthList,
            redshifts,
            extraSurveyConstraints,
            pathToOutputPlotFolder=pathToOutputPlotDirectory)
        result_log += """%s""" % (discoveryMapLink,)
        ratioMapLink = r.plot_sn_discovery_ratio_map(
            log,
            snSurveyDiscoveryTimes,
            redshifts,
            peakAppMagList,
            snCampaignLengthList,
            extraSurveyConstraints,
            pathToOutputPlotFolder=pathToOutputPlotDirectory)
        result_log += """%s""" % (ratioMapLink,)
        result_log += """%s""" % (snRatePlotLink,)

        now = datetime.now()
        now = now.strftime("%Y%m%dt%H%M%S")
        mdLogPath = pathToResultsFolder + \
            "simulation_result_log_%s.md" % (now,)
        mdLog = open(mdLogPath, 'w')
        mdLog.write(result_log)
        mdLog.close()

        dmd.convert_to_html(
            log=log,
            pathToMMDFile=mdLogPath,
            css="amblin"
        )

        print "Results can be found here: %(pathToResultsFolder)s" % locals()
        html = mdLogPath.replace(".md", ".html")
        print "Open this file in your browser: %(html)s" % locals()

    if "dbConn" in locals() and dbConn:
        dbConn.commit()
        dbConn.close()
    ## FINISH LOGGING ##
    endTime = times.get_now_sql_datetime()
    runningTime = times.calculate_time_difference(startTime, endTime)
    log.info('-- FINISHED ATTEMPT TO RUN THE cl_utils.py AT %s (RUNTIME: %s) --' %
             (endTime, runningTime, ))

    return
777910
def _prepare_general_arguments(RRTMGobject):
    ''''''
    tlay = _climlab_to_rrtm(RRTMGobject.Tatm)
    tlev = _climlab_to_rrtm(interface_temperature(**RRTMGobject.state))
    play = _climlab_to_rrtm(RRTMGobject.lev * np.ones_like(tlay))
    plev = _climlab_to_rrtm(RRTMGobject.lev_bounds * np.ones_like(tlev))
    ncol, nlay = tlay.shape
    tsfc = _climlab_to_rrtm_sfc(RRTMGobject.Ts, RRTMGobject.Ts)
    # GASES -- put them in proper dimensions and units
    vapor_mixing_ratio = mmr_to_vmr(RRTMGobject.specific_humidity, gas='H2O')
    h2ovmr   = _climlab_to_rrtm(vapor_mixing_ratio * np.ones_like(RRTMGobject.Tatm))
    o3vmr    = _climlab_to_rrtm(RRTMGobject.absorber_vmr['O3'] * np.ones_like(RRTMGobject.Tatm))
    co2vmr   = _climlab_to_rrtm(RRTMGobject.absorber_vmr['CO2'] * np.ones_like(RRTMGobject.Tatm))
    ch4vmr   = _climlab_to_rrtm(RRTMGobject.absorber_vmr['CH4'] * np.ones_like(RRTMGobject.Tatm))
    n2ovmr   = _climlab_to_rrtm(RRTMGobject.absorber_vmr['N2O'] * np.ones_like(RRTMGobject.Tatm))
    o2vmr    = _climlab_to_rrtm(RRTMGobject.absorber_vmr['O2'] * np.ones_like(RRTMGobject.Tatm))
    cfc11vmr = _climlab_to_rrtm(RRTMGobject.absorber_vmr['CFC11'] * np.ones_like(RRTMGobject.Tatm))
    cfc12vmr = _climlab_to_rrtm(RRTMGobject.absorber_vmr['CFC12'] * np.ones_like(RRTMGobject.Tatm))
    cfc22vmr = _climlab_to_rrtm(RRTMGobject.absorber_vmr['CFC22'] * np.ones_like(RRTMGobject.Tatm))
    ccl4vmr  = _climlab_to_rrtm(RRTMGobject.absorber_vmr['CCL4'] * np.ones_like(RRTMGobject.Tatm))
    #  Cloud parameters
    cldfrac = _climlab_to_rrtm(RRTMGobject.cldfrac * np.ones_like(RRTMGobject.Tatm))
    ciwp = _climlab_to_rrtm(RRTMGobject.ciwp * np.ones_like(RRTMGobject.Tatm))
    clwp = _climlab_to_rrtm(RRTMGobject.clwp * np.ones_like(RRTMGobject.Tatm))
    relq = _climlab_to_rrtm(RRTMGobject.r_liq * np.ones_like(RRTMGobject.Tatm))
    reic = _climlab_to_rrtm(RRTMGobject.r_ice * np.ones_like(RRTMGobject.Tatm))

    return (ncol, nlay, play, plev, tlay, tlev, tsfc,
            h2ovmr, o3vmr, co2vmr, ch4vmr, n2ovmr, o2vmr, cfc11vmr,
            cfc12vmr, cfc12vmr, cfc22vmr, ccl4vmr,
            cldfrac, ciwp, clwp, relq, reic)
702598
def prepare(self):
        
        global print_prepare
        T = self.tm.trace
        if T: T.write(self.trace_message(u'Task.prepare()', self.node))

        # Now that it's the appropriate time, give the TaskMaster a
        # chance to raise any exceptions it encountered while preparing
        # this task.
        self.exception_raise()

        if self.tm.message:
            self.display(self.tm.message)
            self.tm.message = None

        # Let the targets take care of any necessary preparations.
        # This includes verifying that all of the necessary sources
        # and dependencies exist, removing the target file(s), etc.
        #
        # As of April 2008, the get_executor().prepare() method makes
        # sure that all of the aggregate sources necessary to build this
        # Task's target(s) exist in one up-front check.  The individual
        # target t.prepare() methods check that each target's explicit
        # or implicit dependencies exists, and also initialize the
        # .sconsign info.
        executor = self.targets[0].get_executor()
        if executor is None:
            return
        executor.prepare()
        for t in executor.get_action_targets():
            if print_prepare:
                print("Preparing target %s..."%t)
                for s in t.side_effects:
                    print("...with side-effect %s..."%s)
            t.prepare()
            for s in t.side_effects:
                if print_prepare:
                    print("...Preparing side-effect %s..."%s)
                s.prepare()
288020
def run_analysis(self, argv):
        """"""
        args = self._parser.parse_args(argv)

        if not HAVE_ST:
            raise RuntimeError(
                "Trying to run fermipy analysis, but don't have ST")

        gta = GTAnalysis(args.config, logging={'verbosity': 3},
                         fileio={'workdir_regex': '\.xml$|\.npy$'})

        gta.setup(overwrite=False)

        baseline_roi_fit(gta, make_plots=args.make_plots,
                         minmax_npred=[1e3, np.inf])

        localize_sources(gta, nstep=5, dtheta_max=0.5, update=True,
                         prefix='base', make_plots=args.make_plots)

        gta.find_sources(sqrt_ts_threshold=5.0, search_skydir=gta.roi.skydir,
                         search_minmax_radius=[1.0, np.nan])
        gta.optimize()
        gta.print_roi()
        gta.print_params()

        gta.free_sources(skydir=gta.roi.skydir, distance=1.0, pars='norm')
        gta.fit(covar=True)
        gta.print_roi()
        gta.print_params()

        gta.write_roi(args.roi_baseline, make_plots=args.make_plots)
1145554
def prepare_simseries(self, ramflag: bool = True) -> None:
        """"""
        for node in printtools.progressbar(self):
            node.prepare_simseries(ramflag)
367946
def run_analysis(self, argv):
        """"""
        args = self._parser.parse_args(argv)

        if not HAVE_ST:
            raise RuntimeError(
                "Trying to run fermipy analysis, but don't have ST")

        workdir = os.path.dirname(args.config)
        _config_file = self._clone_config_and_srcmaps(args.config, args.seed)

        gta = GTAnalysis(_config_file, logging={'verbosity': 3},
                         fileio={'workdir_regex': '\.xml$|\.npy$'})
        gta.load_roi(args.roi_baseline)

        simfile = os.path.join(workdir, 'sim_%s_%s.yaml' %
                               (args.sim, args.sim_profile))

        mcube_file = "%s_%s_%06i" % (args.sim, args.sim_profile, args.seed)
        sim_config = utils.load_yaml(simfile)

        injected_source = sim_config.get('injected_source', None)
        if injected_source is not None:
            src_dict =  injected_source['source_model']
            src_dict['ra'] = gta.config['selection']['ra']
            src_dict['dec'] = gta.config['selection']['dec']
            injected_name = injected_source['name']
            gta.add_source(injected_name, src_dict)
            gta.write_model_map(mcube_file)
            mc_spec_dict = dict(true_counts=gta.model_counts_spectrum(injected_name),
                                energies=gta.energies,
                                model=src_dict)
            mcspec_file = os.path.join(workdir,
                                       "mcspec_%s_%06i.yaml" % (mcube_file, args.seed))
            utils.write_yaml(mc_spec_dict, mcspec_file)
        else:
            injected_name = None

        gta.write_roi('sim_baseline_%06i' % args.seed)

        test_sources = []
        for profile in args.profiles:
            profile_path = os.path.join(workdir, 'profile_%s.yaml' % profile)
            test_source = load_yaml(profile_path)
            test_sources.append(test_source)
            first = args.seed
            last = first + args.nsims
            for seed in range(first, last):
                self._run_simulation(gta, args.roi_baseline,
                                     injected_name, test_sources, first, seed,
                                     non_null_src=args.non_null_src)
869557
def run_analysis(self, argv):
        """"""
        args = self._parser.parse_args(argv)

        if not HAVE_ST:
            raise RuntimeError(
                "Trying to run fermipy analysis, but don't have ST")

        if args.load_baseline:
            gta = GTAnalysis.create(args.roi_baseline,
                                    args.config)
        else:
            gta = GTAnalysis(args.config,
                             logging={'verbosity': 3},
                             fileio={'workdir_regex': '\.xml$|\.npy$'})
            gta.setup()
            if is_not_null(args.input_pars):
                gta.load_parameters_from_yaml(args.input_pars)
            gta.write_roi(args.roi_baseline,
                          save_model_map=True,
                          save_weight_map=True,
                          make_plots=args.make_plots)

        src_list = get_src_names(gta)
        plotter = plotting.AnalysisPlotter(gta.config['plotting'],
                                           fileio=gta.config['fileio'],
                                           logging=gta.config['logging'])

        if is_null(args.fit_strategy):
            return

        fit_strategy = load_yaml(args.fit_strategy)
        npred_current = None
        npred_prev = None
        
        plots_only = False

        for fit_stage in fit_strategy:
            mask = fit_stage.get('mask', None)
            npred_threshold = fit_stage.get('npred_threshold', 1.0e4)
            frac_threshold = fit_stage.get('frac_threshold', 0.5)
            npred_frac = fit_stage.get('npred_frac', 0.9999)

            if plots_only:
                gta.load_roi("%s.npy" % fit_stage['key'])
                npred_current =  set_wts_get_npred_wt(gta, mask)
                skip_list_region = get_unchanged(src_list,
                                                 npred_current,
                                                 npred_prev,
                                                 frac_threshold=frac_threshold)
            else:
                npred_current =  set_wts_get_npred_wt(gta, mask)
                skip_list_region = get_unchanged(src_list,
                                                 npred_current,
                                                 npred_prev,
                                                 frac_threshold=frac_threshold)     
                gta.optimize(npred_frac=npred_frac, 
                             npred_threshold=npred_threshold,
                             skip=skip_list_region)
            
            snapshot(gta, plotter, fit_stage['key'], make_plots=args.make_plots)
            npred_prev = npred_current
            npred_current = build_srcdict(gta, 'npred_wt')
1125755
def main():
    
    options = handle_options()

    # read in observed and synthetic data
    elecs, d_obs = readin_volt(options.d_obs)
    elecs, d_est = readin_volt(options.d_est)
    elecs, d_estTC = readin_volt(options.d_estTC)
    # calculate corrected data
    volt_corr = calc_correction(d_obs,
                                d_est,
                                d_estTC,
                                )
    # save data
    save_volt(elecs,
              volt_corr,
              options.output,
              )
912278
def _run_rtg_eval(vrn_file, rm_file, rm_interval_file, base_dir, data, validate_method):
    
    out_dir = os.path.join(base_dir, "rtg")
    if not utils.file_exists(os.path.join(out_dir, "done")):
        if os.path.exists(out_dir):
            shutil.rmtree(out_dir)
        vrn_file, rm_file, interval_bed = _prepare_inputs(vrn_file, rm_file, rm_interval_file, base_dir, data)

        rtg_ref = tz.get_in(["reference", "rtg"], data)
        if isinstance(rtg_ref, dict) and "base" in rtg_ref:
            rtg_ref = os.path.dirname(rtg_ref["base"])
        assert rtg_ref and os.path.exists(rtg_ref), ("Did not find rtg indexed reference file for validation:\n%s\n"
                                                     "Run bcbio_nextgen.py upgrade --data --aligners rtg" % rtg_ref)
        # handle CWL where we have a reference to a single file in the RTG directory
        if os.path.isfile(rtg_ref):
            rtg_ref = os.path.dirname(rtg_ref)

        # get core and memory usage from standard configuration
        threads = min(dd.get_num_cores(data), 6)
        resources = config_utils.get_resources("rtg", data["config"])
        memory = config_utils.adjust_opts(resources.get("jvm_opts", ["-Xms500m", "-Xmx1500m"]),
                                          {"algorithm": {"memory_adjust": {"magnitude": threads,
                                                                           "direction": "increase"}}})
        jvm_stack = [x for x in memory if x.startswith("-Xms")]
        jvm_mem = [x for x in memory if x.startswith("-Xmx")]
        jvm_stack = jvm_stack[0] if len(jvm_stack) > 0 else "-Xms500m"
        jvm_mem = jvm_mem[0].replace("-Xmx", "") if len(jvm_mem) > 0 else "3g"
        cmd = ["rtg", "vcfeval", "--threads", str(threads),
               "-b", rm_file, "--bed-regions", interval_bed,
               "-c", vrn_file, "-t", rtg_ref, "-o", out_dir]
        if validate_method == "rtg-squash-ploidy":
            cmd += ["--squash-ploidy"]
        rm_samples = vcfutils.get_samples(rm_file)
        if len(rm_samples) > 1 and dd.get_sample_name(data) in rm_samples:
            cmd += ["--sample=%s" % dd.get_sample_name(data)]
        cmd += ["--vcf-score-field='%s'" % (_pick_best_quality_score(vrn_file))]
        mem_export = "%s export RTG_JAVA_OPTS='%s' && export RTG_MEM=%s" % (utils.local_path_export(),
                                                                            jvm_stack, jvm_mem)
        cmd = mem_export + " && " + " ".join(cmd)
        do.run(cmd, "Validate calls using rtg vcfeval", data)
    out = {"fp": os.path.join(out_dir, "fp.vcf.gz"),
           "fn": os.path.join(out_dir, "fn.vcf.gz")}
    tp_calls = os.path.join(out_dir, "tp.vcf.gz")
    tp_baseline = os.path.join(out_dir, "tp-baseline.vcf.gz")
    if os.path.exists(tp_baseline):
        out["tp"] = tp_baseline
        out["tp-calls"] = tp_calls
    else:
        out["tp"] = tp_calls
    return out
1187554
def run_analysis(self, argv):
        
        args = self._parser.parse_args(argv)
        components = Component.build_from_yamlfile(args.comp)
        NAME_FACTORY.update_base_dict(args.data)
        model_dict = make_library(**args.__dict__)
        model_manager = model_dict['ModelManager']
        models = load_yaml(args.models)
        data = args.data
        hpx_order = args.hpx_order
        for modelkey in models:
            model_manager.make_srcmap_manifest(modelkey, components, data)
            model_manager.make_fermipy_config_yaml(modelkey, components, data,
                                                   hpx_order=hpx_order,
                                                   irf_ver=NAME_FACTORY.irf_ver())
1339560
def set_plot_tcrhoc(self,linestyle=['-'],burn_limit=0.997,marker=['o'],markevery=500,end_model=[-1],deg_line=True): 
		'''
			
		'''
		if len(linestyle)==0:
			linestyle=200*['-']

                for i in range(len(self.runs_H5_surf)):
                        m1p65_last=se(self.runs_H5_out[i])
			t1_model=-1
			if end_model[i] != -1:
				t1_model=end_model[i]

			model=m1p65_last.se.cycles
			model_list=[]
			for k in range(0,len(model),5):
				model_list.append(model[k])
                        rho1=m1p65_last.get(model_list,'rho')   #[:(t1_model-t0_model)]
                        T1=m1p65_last.get(model_list,'temperature')#[:(t1_model-t0_model)]
			T_unit=m1p65_last.get('temperature_unit')
			mass=m1p65_last.get(model_list,'mass')
			mini=m1p65_last.get('mini')
			zini=m1p65_last.get('zini')
			#info=m1p65_last.burnstage_upgrade()
			#burn_info=[]
			#burn_cycle=[]
			#for k in range(len(info[0])):
			#    if 'start' in info[3][k]:
			#	burn_info.append( info[3][k])
			#	burn_cycle.append(info[0][k])
			#print burn_info
			#print burn_cycle
			'''
			#H
			#Maybe use get_elemental_abunds(), too slow if individually taken
			h=m1p65_last.get(model_list,'H-1')
			h_ini=mini*h[0][-1]
			h_limit=h_ini*burn_limit
			print h_limit
			#He
			he=m1p65_last.get(model_list,'He-4')
			he+=m1p65_last.get(model_list,'He-3')
			he_ini=mini*he[0][-1]
			he_limit=he_ini*burn_limit
			print he_limit
			if mini>5:
				c=m1p65_last.get(model_list,'C-12')
				c_ini=mini*c[0][-1]
				c_limit=c_ini*burn_limit
			if mini>10:
				ne=m1p65_last.get(model_list,'Ne-20')
				ne_ini=	mini*ne[0][-1]
				ne_limit=ne_ini*burn_limit
				o=m1p65_last.get(model_list,'O-16')
				o_ini=mini*o[0][-1]
				o_limit=o_ini*burn_limit
				si=m1p65_last.get(model_list,'Si-28')
				si_ini=mini*si[0][-1]
				si_limit=si_ini*burn_limit
			print '#############################3'
			print h_ini,he_ini
			rho=[]
			T=[]
			he_depl_idx=-1
			h_depl_idx=-1
			c_depl_idx=-1
			ne_depl_idx=-1
			o_depl_idx=-1
			si_depl_idx=-1
			'''
			rho=[]
			T=[]
			for k in range(len(model_list)):
				rho_center=rho1[k][0]
				T_center=T1[k][0]
				rho.append(rho_center)
				T.append(T_center)
				'''
				mass11=np.array([0]+list(mass[k]))
				delta_mass=mass11[1:]-mass11[:-1]
				#for the low-mass + massive AGB
				if (sum(np.array(he[k])*np.array(delta_mass))<he_limit) and  (he_depl_idx==-1):
					he_depl_idx=k			
                                if (sum(np.array(h[k])*np.array(delta_mass))<h_limit) and  (h_depl_idx==-1):
                                        h_depl_idx=k
				#for the SAGB + massive
				if mini>5:
					if (sum(np.array(c[k])*np.array(delta_mass))<c_limit) and  (c_depl_idx==-1):
                                        	c_depl_idx=k
				#for the massive stars
				if mini>10:
					if (sum(np.array(ne[k])*np.array(delta_mass))<ne_limit) and  (ne_depl_idx==-1):
                                        	ne_depl_idx=k
					if (sum(np.array(o[k])*np.array(delta_mass))<o_limit) and  (o_depl_idx==-1):
                                        	o_depl_idx=k
					if (sum(np.array(si[k])*np.array(delta_mass))<si_limit) and  (si_depl_idx==-1):
                                        	si_depl_idx=k
				'''

			T=np.log10(np.array(T)*T_unit)
			#T_degeneracy=np.log10(np.array(rho)**(2./3.))
			rho=np.log10(np.array(rho))
                        figure(1)
                        #plot(logTeff,logL,marker=symbs[i],label=self.run_label[i],linestyle=linestyle[i],markevery=markevery)
			#pl.plot(rho,T,marker=symbs[i],label=self.run_label[i],linestyle=linestyle[i],markevery=markevery)
			label=str(mini)+'$M_{\odot}$, Z='+str(zini)
			plt.plot(rho,T,label=label,linestyle=linestyle[i],marker=marker[i],markevery=markevery)
			'''
			#plt.plot(rho,T_degeneracy,color='b',linestyle='-',label='$P_e = P_{e,deg}$')
			print burn_info
			print burn_cycle
			burn_cycle=np.array(burn_cycle)/100.
			print 'changed cycles:', burn_cycle
			print len(rho)
			#For different burning stages:
			if 'H_start' in burn_info:
				plt.plot(rho[burn_cycle[0]],T[burn_cycle[0]],marker='o',color='b')
			if 'He_start' in burn_info:
				plt.plot(rho[burn_cycle[1]],T[burn_cycle[1]],marker='o',color='r')
			if 'C_start' in burn_info:
				plt.plot(rho[burn_cycle[2]],T[burn_cycle[2]],marker='o',color='g')
			if 'Ne_start' in burn_info:
				plt.plot(rho[burn_cycle[3]],T[burn_cycle[3]],marker='D',color='b')
			if 'O_start' in burn_info:
				plt.plot(rho[burn_cycle[4]],T[burn_cycle[4]],marker='D',color='r')
			if 'Si_start' in burn_info:
				plt.plot(rho[burn_cycle[5]],T[burn_cycle[5]],marker='D',color='g')
			'''
		
			#ax = plt.gca()
			#ax.invert_xaxis()
			plt.rcParams.update({'font.size': 16})
			plt.rc('xtick', labelsize=16)
			plt.rc('ytick', labelsize=16)
			legend(loc=4)
		        plt.xlabel('log $\\rho_{\\rm c}$',fontsize=18)
        		plt.ylabel('log $T_{\\rm c}$',fontsize=18)
			'''
			#plt.gca().invert_xaxis()	
			figure(0)
			#pl.plot(rho,T,marker=symbs[i],label=self.run_label[i],linestyle=linestyle[i],markevery=markevery)
                        plt.plot(rho,T,label=self.extra_label[i],linestyle=linestyle[i])
			plt.plot(rho,T_degeneracy,color='b',linestyle='-',label='$P_e = P_{e,deg}$')
			#For different burning stages:
			if not h_depl_idx ==-1:
				plt.plot(rho[h_depl_idx],T[h_depl_idx],marker='o',color='b')
			if not he_depl_idx ==-1:
				plt.plot(rho[he_depl_idx],T[he_depl_idx],marker='o',color='r')
			if not c_depl_idx ==-1:
				plt.plot(rho[c_depl_idx],T[c_depl_idx],marker='o',color='g')
			if not ne_depl_idx ==-1:
				plt.plot(rho[ne_depl_idx],T[ne_depl_idx],marker='D',color='b')
			if not o_depl_idx ==-1:
				plt.plot(rho[o_depl_idx],T[o_depl_idx],marker='D',color='r')
			if not si_depl_idx ==-1:
				plt.plot(rho[si_depl_idx],T[si_depl_idx],marker='D',color='g')
			

			#ax = plt.gca()
			#ax.invert_xaxis()
			plt.rcParams.update({'font.size': 16})
                	plt.rc('xtick', labelsize=16)
                	plt.rc('ytick', labelsize=16)
                	legend(loc=4)
		        plt.xlabel('log $\\rho_{\\rm c}$',fontsize=18)
        		plt.ylabel('log $T_{\\rm c}$',fontsize=18)
			#plt.gca().invert_xaxis()
			'''
			i+=1	
		if deg_line==True:
			rho=np.arange(0,9,0.01)
			T_degeneracy=2./3. *rho +np.log10(1.207e5 * 1.8/(2.**(5./3.)))
			#T_degeneracy=np.log10( (10**np.array(rho))**(2./3.))
			plt.plot(rho,T_degeneracy,color='b',linestyle='-',label='$P_e = P_{e,deg}$')
		plt.legend(loc=2)
596118
def build_args():
    """"""

    parser = argparse.ArgumentParser(description=DESCRIPTION, formatter_class=argparse.RawDescriptionHelpFormatter)
    parser.add_argument(u'sensor_graph', type=str, help=u"The sensor graph file to load and run.")
    parser.add_argument(u'--stop', u'-s', action=u"append", default=[], type=str, help=u"A stop condition for when the simulation should end.")
    parser.add_argument(u'--realtime', u'-r', action=u"store_true", help=u"Do not accelerate the simulation, pin the ticks to wall clock time")
    parser.add_argument(u'--watch', u'-w', action=u"append", default=[], help=u"A stream to watch and print whenever writes are made.")
    parser.add_argument(u'--trace', u'-t', help=u"Trace all writes to output streams to a file")
    parser.add_argument(u'--disable-optimizer', action="store_true", help=u"disable the sensor graph optimizer completely")
    parser.add_argument(u"--mock-rpc", u"-m", action=u"append", type=str, default=[], help=u"mock an rpc, format should be <slot id>:<rpc_id> = value.  For example -m \"slot 1:0x500a = 10\"")
    parser.add_argument(u"--port", u"-p", help=u"The port to use to connect to a device if we are semihosting")
    parser.add_argument(u"--semihost-device", u"-d", type=lambda x: int(x, 0), help=u"The device id of the device we should semihost this sensor graph on.")
    parser.add_argument(u"-c", u"--connected", action="store_true", help=u"Simulate with a user connected to the device (to enable realtime outputs)")
    parser.add_argument(u"-i", u"--stimulus", action=u"append", default=[], help="Push a value to an input stream at the specified time (or before starting).  The syntax is [time: ][system ]input X = Y where X and Y are integers")
    return parser