astrobase/periodbase/kbls.py
def _get_bls_stats(stimes,
                   smags,
                   serrs,
                   thistransdepth,
                   thistransduration,
                   ingressdurationfraction,
                   nphasebins,
                   thistransingressbin,
                   thistransegressbin,
                   thisbestperiod,
                   thisnphasebins,
                   magsarefluxes=False,
                   verbose=False):
    '''
    Actually calculates the stats.

    '''

    try:

        # try getting the minimum light epoch using the phase bin method
        me_epochbin = int((thistransegressbin +
                           thistransingressbin)/2.0)

        me_phases = (
            (stimes - stimes.min())/thisbestperiod -
            npfloor((stimes - stimes.min())/thisbestperiod)
        )
        me_phases_sortind = npargsort(me_phases)
        me_sorted_phases = me_phases[me_phases_sortind]
        me_sorted_times = stimes[me_phases_sortind]

        me_bins = nplinspace(0.0, 1.0, thisnphasebins)
        me_bininds = npdigitize(me_sorted_phases, me_bins)

        me_centertransit_ind = me_bininds == me_epochbin
        me_centertransit_phase = (
            npmedian(me_sorted_phases[me_centertransit_ind])
        )
        me_centertransit_timeloc = npwhere(
            npabs(me_sorted_phases - me_centertransit_phase) ==
            npmin(npabs(me_sorted_phases - me_centertransit_phase))
        )
        me_centertransit_time = me_sorted_times[
            me_centertransit_timeloc
        ]

        if me_centertransit_time.size > 1:
            LOGWARNING('multiple possible times-of-center transits '
                       'found for period %.7f, picking the first '
                       'one from: %s' %
                       (thisbestperiod, repr(me_centertransit_time)))

        thisminepoch = me_centertransit_time[0]

    except Exception as e:

        LOGEXCEPTION(
            'could not determine the center time of transit for '
            'the phased LC, trying SavGol fit instead...'
        )
        # fit a Savitsky-Golay instead and get its minimum
        savfit = savgol_fit_magseries(stimes, smags, serrs,
                                      thisbestperiod,
                                      magsarefluxes=magsarefluxes,
                                      verbose=verbose,
                                      sigclip=None)
        thisminepoch = savfit['fitinfo']['fitepoch']


    if isinstance(thisminepoch, npndarray):
        if verbose:
            LOGWARNING('minimum epoch is actually an array:\n'
                       '%s\n'
                       'instead of a float, '
                       'are there duplicate time values '
                       'in the original input? '
                       'will use the first value in this array.'
                       % repr(thisminepoch))
        thisminepoch = thisminepoch[0]

    # set up trapezoid transit model to fit for this LC
    transitparams = [
        thisbestperiod,
        thisminepoch,
        thistransdepth,
        thistransduration,
        ingressdurationfraction*thistransduration
    ]

    modelfit = traptransit_fit_magseries(
        stimes,
        smags,
        serrs,
        transitparams,
        sigclip=None,
        magsarefluxes=magsarefluxes,
        verbose=verbose
    )

    # if the model fit succeeds, calculate SNR using the trapezoid model fit
    if modelfit and modelfit['fitinfo']['finalparams'] is not None:

        fitparams = modelfit['fitinfo']['finalparams']
        fiterrs = modelfit['fitinfo']['finalparamerrs']
        modelmags, actualmags, modelphase = (
            modelfit['fitinfo']['fitmags'],
            modelfit['magseries']['mags'],
            modelfit['magseries']['phase']
        )
        subtractedmags = actualmags - modelmags
        subtractedrms = npstd(subtractedmags)
        fit_period, fit_epoch, fit_depth, fit_duration, fit_ingress_dur = (
            fitparams
        )

        npts_in_transit = modelfit['fitinfo']['ntransitpoints']
        transit_snr = (
            npsqrt(npts_in_transit) * npabs(fit_depth/subtractedrms)
        )

        if verbose:

            LOGINFO('refit best period: %.6f, '
                    'refit center of transit: %.5f' %
                    (fit_period, fit_epoch))

            LOGINFO('npoints in transit: %s' % npts_in_transit)

            LOGINFO('transit depth (delta): %.5f, '
                    'frac transit length (q): %.3f, '
                    ' SNR: %.3f' %
                    (fit_depth,
                     fit_duration,
                     transit_snr))

        return {'period':fit_period,
                'epoch':fit_epoch,
                'snr':transit_snr,
                'transitdepth':fit_depth,
                'transitduration':fit_duration,
                'nphasebins':nphasebins,
                'transingressbin':thistransingressbin,
                'transegressbin':thistransegressbin,
                'npoints_in_transit':npts_in_transit,
                'blsmodel':modelmags,
                'subtractedmags':subtractedmags,
                'phasedmags':actualmags,
                'phases':modelphase,
                'fitparams':fitparams,
                'fiterrs':fiterrs,
                'fitinfo':modelfit}


    # if the model fit doesn't work, then do the SNR calculation the old way
    else:

        # phase using this epoch
        phased_magseries = phase_magseries_with_errs(stimes,
                                                     smags,
                                                     serrs,
                                                     thisbestperiod,
                                                     thisminepoch,
                                                     wrap=False,
                                                     sort=True)

        tphase = phased_magseries['phase']
        tmags = phased_magseries['mags']

        # use the transit depth and duration to subtract the BLS transit
        # model from the phased mag series. we're centered about 0.0 as the
        # phase of the transit minimum so we need to look at stuff from
        # [0.0, transitphase] and [1.0-transitphase, 1.0]
        transitphase = thistransduration/2.0

        transitindices = ((tphase < transitphase) |
                          (tphase > (1.0 - transitphase)))

        # this is the BLS model
        # constant = median(tmags) outside transit
        # constant = thistransitdepth inside transit
        blsmodel = npfull_like(tmags, npmedian(tmags))

        if magsarefluxes:

            # eebls.f returns +ve transit depth for fluxes
            # so we need to subtract here to get fainter fluxes in transit
            blsmodel[transitindices] = (
                blsmodel[transitindices] - thistransdepth
            )
        else:

            # eebls.f returns -ve transit depth for magnitudes
            # so we need to subtract here to get fainter mags in transits
            blsmodel[transitindices] = (
                blsmodel[transitindices] - thistransdepth
            )

        # see __init__/get_snr_of_dip docstring for description of transit
        # SNR equation, which is what we use for `thissnr`.
        subtractedmags = tmags - blsmodel
        subtractedrms = npstd(subtractedmags)
        npts_in_transit = len(tmags[transitindices])
        thissnr = (
            npsqrt(npts_in_transit) * npabs(thistransdepth/subtractedrms)
        )

        # tell user about stuff if verbose = True
        if verbose:

            LOGINFO('refit best period: %.6f, '
                    'refit center of transit: %.5f' %
                    (thisbestperiod, thisminepoch))

            LOGINFO('transit ingress phase = %.3f to %.3f' % (1.0 -
                                                              transitphase,
                                                              1.0))
            LOGINFO('transit egress phase = %.3f to %.3f' % (0.0,
                                                             transitphase))
            LOGINFO('npoints in transit: %s' % tmags[transitindices].size)

            LOGINFO('transit depth (delta): %.5f, '
                    'frac transit length (q): %.3f, '
                    ' SNR: %.3f' %
                    (thistransdepth,
                     thistransduration,
                     thissnr))

        return {'period':thisbestperiod,
                'epoch':thisminepoch,
                'snr':thissnr,
                'transitdepth':thistransdepth,
                'transitduration':thistransduration,
                'nphasebins':nphasebins,
                'transingressbin':thistransingressbin,
                'transegressbin':thistransegressbin,
                'blsmodel':blsmodel,
                'subtractedmags':subtractedmags,
                'phasedmags':tmags,
                'phases':tphase}
def _get_bls_stats(stimes,
                   smags,
                   serrs,
                   thistransdepth,
                   thistransduration,
                   ingressdurationfraction,
                   nphasebins,
                   thistransingressbin,
                   thistransegressbin,
                   thisbestperiod,
                   thisnphasebins,
                   magsarefluxes=False,
                   verbose=False):
    '''
    Actually calculates the stats.

    '''

    try:

        # try getting the minimum light epoch using the phase bin method
        me_epochbin = int((thistransegressbin +
                           thistransingressbin)/2.0)

        me_phases = (
            (stimes - stimes.min())/thisbestperiod -
            npfloor((stimes - stimes.min())/thisbestperiod)
        )
        me_phases_sortind = npargsort(me_phases)
        me_sorted_phases = me_phases[me_phases_sortind]
        me_sorted_times = stimes[me_phases_sortind]

        me_bins = nplinspace(0.0, 1.0, thisnphasebins)
        me_bininds = npdigitize(me_sorted_phases, me_bins)

        me_centertransit_ind = me_bininds == me_epochbin
        me_centertransit_phase = (
            npmedian(me_sorted_phases[me_centertransit_ind])
        )
        me_centertransit_timeloc = npwhere(
            npabs(me_sorted_phases - me_centertransit_phase) ==
            npmin(npabs(me_sorted_phases - me_centertransit_phase))
        )
        me_centertransit_time = me_sorted_times[
            me_centertransit_timeloc
        ]

        if me_centertransit_time.size > 1:
            LOGWARNING('multiple possible times-of-center transits '
                       'found for period %.7f, picking the first '
                       'one from: %s' %
                       (thisbestperiod, repr(me_centertransit_time)))

        thisminepoch = me_centertransit_time[0]

    except Exception as e:

        LOGEXCEPTION(
            'could not determine the center time of transit for '
            'the phased LC, trying SavGol fit instead...'
        )
        # fit a Savitsky-Golay instead and get its minimum
        savfit = savgol_fit_magseries(stimes, smags, serrs,
                                      thisbestperiod,
                                      magsarefluxes=magsarefluxes,
                                      verbose=verbose,
                                      sigclip=None)
        thisminepoch = savfit['fitinfo']['fitepoch']


    if isinstance(thisminepoch, npndarray):
        if verbose:
            LOGWARNING('minimum epoch is actually an array:\n'
                       '%s\n'
                       'instead of a float, '
                       'are there duplicate time values '
                       'in the original input? '
                       'will use the first value in this array.'
                       % repr(thisminepoch))
        thisminepoch = thisminepoch[0]

    # set up trapezoid transit model to fit for this LC
    transitparams = [
        thisbestperiod,
        thisminepoch,
        thistransdepth,
        thistransduration,
        ingressdurationfraction*thistransduration
    ]

    modelfit = traptransit_fit_magseries(
        stimes,
        smags,
        serrs,
        transitparams,
        sigclip=None,
        magsarefluxes=magsarefluxes,
        verbose=verbose
    )

    # if the model fit succeeds, calculate SNR using the trapezoid model fit
    if modelfit and modelfit['fitinfo']['finalparams'] is not None:

        fitparams = modelfit['fitinfo']['finalparams']
        fiterrs = modelfit['fitinfo']['finalparamerrs']
        modelmags, actualmags, modelphase = (
            modelfit['fitinfo']['fitmags'],
            modelfit['magseries']['mags'],
            modelfit['magseries']['phase']
        )
        subtractedmags = actualmags - modelmags
        subtractedrms = npstd(subtractedmags)
        fit_period, fit_epoch, fit_depth, fit_duration, fit_ingress_dur = (
            fitparams
        )

        npts_in_transit = modelfit['fitinfo']['ntransitpoints']
        transit_snr = (
            npsqrt(npts_in_transit) * npabs(fit_depth/subtractedrms)
        )

        if verbose:

            LOGINFO('refit best period: %.6f, '
                    'refit center of transit: %.5f' %
                    (fit_period, fit_epoch))

            LOGINFO('npoints in transit: %s' % npts_in_transit)

            LOGINFO('transit depth (delta): %.5f, '
                    'frac transit length (q): %.3f, '
                    ' SNR: %.3f' %
                    (fit_depth,
                     fit_duration,
                     transit_snr))

        return {'period':fit_period,
                'epoch':fit_epoch,
                'snr':transit_snr,
                'transitdepth':fit_depth,
                'transitduration':fit_duration,
                'nphasebins':nphasebins,
                'transingressbin':thistransingressbin,
                'transegressbin':thistransegressbin,
                'npoints_in_transit':npts_in_transit,
                'blsmodel':modelmags,
                'subtractedmags':subtractedmags,
                'phasedmags':actualmags,
                'phases':modelphase,
                'fitparams':fitparams,
                'fiterrs':fiterrs,
                'fitinfo':modelfit}


    # if the model fit doesn't work, then do the SNR calculation the old way
    else:

        # phase using this epoch
        phased_magseries = phase_magseries_with_errs(stimes,
                                                     smags,
                                                     serrs,
                                                     thisbestperiod,
                                                     thisminepoch,
                                                     wrap=False,
                                                     sort=True)

        tphase = phased_magseries['phase']
        tmags = phased_magseries['mags']

        # use the transit depth and duration to subtract the BLS transit
        # model from the phased mag series. we're centered about 0.0 as the
        # phase of the transit minimum so we need to look at stuff from
        # [0.0, transitphase] and [1.0-transitphase, 1.0]
        transitphase = thistransduration/2.0

        transitindices = ((tphase < transitphase) |
                          (tphase > (1.0 - transitphase)))

        # this is the BLS model
        # constant = median(tmags) outside transit
        # constant = thistransitdepth inside transit
        blsmodel = npfull_like(tmags, npmedian(tmags))

        if magsarefluxes:

            # eebls.f returns +ve transit depth for fluxes
            # so we need to subtract here to get fainter fluxes in transit
            blsmodel[transitindices] = (
                blsmodel[transitindices] - thistransdepth
            )
        else:

            # eebls.f returns -ve transit depth for magnitudes
            # so we need to subtract here to get fainter mags in transits
            blsmodel[transitindices] = (
                blsmodel[transitindices] - thistransdepth
            )

        # see __init__/get_snr_of_dip docstring for description of transit
        # SNR equation, which is what we use for `thissnr`.
        subtractedmags = tmags - blsmodel
        subtractedrms = npstd(subtractedmags)
        npts_in_transit = len(tmags[transitindices])
        thissnr = (
            npsqrt(npts_in_transit) * npabs(thistransdepth/subtractedrms)
        )

        # tell user about stuff if verbose = True
        if verbose:

            LOGINFO('refit best period: %.6f, '
                    'refit center of transit: %.5f' %
                    (thisbestperiod, thisminepoch))

            LOGINFO('transit ingress phase = %.3f to %.3f' % (1.0 -
                                                              transitphase,
                                                              1.0))
            LOGINFO('transit egress phase = %.3f to %.3f' % (0.0,
                                                             transitphase))
            LOGINFO('npoints in transit: %s' % tmags[transitindices].size)

            LOGINFO('transit depth (delta): %.5f, '
                    'frac transit length (q): %.3f, '
                    ' SNR: %.3f' %
                    (thistransdepth,
                     thistransduration,
                     thissnr))

        return {'period':thisbestperiod,
                'epoch':thisminepoch,
                'snr':thissnr,
                'transitdepth':thistransdepth,
                'transitduration':thistransduration,
                'nphasebins':nphasebins,
                'transingressbin':thistransingressbin,
                'transegressbin':thistransegressbin,
                'blsmodel':blsmodel,
                'subtractedmags':subtractedmags,
                'phasedmags':tmags,
                'phases':tphase}
def _get_bls_stats ( stimes , smags , serrs , thistransdepth , thistransduration , ingressdurationfraction , nphasebins , thistransingressbin , thistransegressbin , thisbestperiod , thisnphasebins , magsarefluxes = False , verbose = False ) : try : # try getting the minimum light epoch using the phase bin method me_epochbin = int ( ( thistransegressbin + thistransingressbin ) / 2.0 ) me_phases = ( ( stimes - stimes . min ( ) ) / thisbestperiod - npfloor ( ( stimes - stimes . min ( ) ) / thisbestperiod ) ) me_phases_sortind = npargsort ( me_phases ) me_sorted_phases = me_phases [ me_phases_sortind ] me_sorted_times = stimes [ me_phases_sortind ] me_bins = nplinspace ( 0.0 , 1.0 , thisnphasebins ) me_bininds = npdigitize ( me_sorted_phases , me_bins ) me_centertransit_ind = me_bininds == me_epochbin me_centertransit_phase = ( npmedian ( me_sorted_phases [ me_centertransit_ind ] ) ) me_centertransit_timeloc = npwhere ( npabs ( me_sorted_phases - me_centertransit_phase ) == npmin ( npabs ( me_sorted_phases - me_centertransit_phase ) ) ) me_centertransit_time = me_sorted_times [ me_centertransit_timeloc ] if me_centertransit_time . size > 1 : LOGWARNING ( 'multiple possible times-of-center transits ' 'found for period %.7f, picking the first ' 'one from: %s' % ( thisbestperiod , repr ( me_centertransit_time ) ) ) thisminepoch = me_centertransit_time [ 0 ] except Exception as e : LOGEXCEPTION ( 'could not determine the center time of transit for ' 'the phased LC, trying SavGol fit instead...' ) # fit a Savitsky-Golay instead and get its minimum savfit = savgol_fit_magseries ( stimes , smags , serrs , thisbestperiod , magsarefluxes = magsarefluxes , verbose = verbose , sigclip = None ) thisminepoch = savfit [ 'fitinfo' ] [ 'fitepoch' ] if isinstance ( thisminepoch , npndarray ) : if verbose : LOGWARNING ( 'minimum epoch is actually an array:\n' '%s\n' 'instead of a float, ' 'are there duplicate time values ' 'in the original input? ' 'will use the first value in this array.' % repr ( thisminepoch ) ) thisminepoch = thisminepoch [ 0 ] # set up trapezoid transit model to fit for this LC transitparams = [ thisbestperiod , thisminepoch , thistransdepth , thistransduration , ingressdurationfraction * thistransduration ] modelfit = traptransit_fit_magseries ( stimes , smags , serrs , transitparams , sigclip = None , magsarefluxes = magsarefluxes , verbose = verbose ) # if the model fit succeeds, calculate SNR using the trapezoid model fit if modelfit and modelfit [ 'fitinfo' ] [ 'finalparams' ] is not None : fitparams = modelfit [ 'fitinfo' ] [ 'finalparams' ] fiterrs = modelfit [ 'fitinfo' ] [ 'finalparamerrs' ] modelmags , actualmags , modelphase = ( modelfit [ 'fitinfo' ] [ 'fitmags' ] , modelfit [ 'magseries' ] [ 'mags' ] , modelfit [ 'magseries' ] [ 'phase' ] ) subtractedmags = actualmags - modelmags subtractedrms = npstd ( subtractedmags ) fit_period , fit_epoch , fit_depth , fit_duration , fit_ingress_dur = ( fitparams ) npts_in_transit = modelfit [ 'fitinfo' ] [ 'ntransitpoints' ] transit_snr = ( npsqrt ( npts_in_transit ) * npabs ( fit_depth / subtractedrms ) ) if verbose : LOGINFO ( 'refit best period: %.6f, ' 'refit center of transit: %.5f' % ( fit_period , fit_epoch ) ) LOGINFO ( 'npoints in transit: %s' % npts_in_transit ) LOGINFO ( 'transit depth (delta): %.5f, ' 'frac transit length (q): %.3f, ' ' SNR: %.3f' % ( fit_depth , fit_duration , transit_snr ) ) return { 'period' : fit_period , 'epoch' : fit_epoch , 'snr' : transit_snr , 'transitdepth' : fit_depth , 'transitduration' : fit_duration , 'nphasebins' : nphasebins , 'transingressbin' : thistransingressbin , 'transegressbin' : thistransegressbin , 'npoints_in_transit' : npts_in_transit , 'blsmodel' : modelmags , 'subtractedmags' : subtractedmags , 'phasedmags' : actualmags , 'phases' : modelphase , 'fitparams' : fitparams , 'fiterrs' : fiterrs , 'fitinfo' : modelfit } # if the model fit doesn't work, then do the SNR calculation the old way else : # phase using this epoch phased_magseries = phase_magseries_with_errs ( stimes , smags , serrs , thisbestperiod , thisminepoch , wrap = False , sort = True ) tphase = phased_magseries [ 'phase' ] tmags = phased_magseries [ 'mags' ] # use the transit depth and duration to subtract the BLS transit # model from the phased mag series. we're centered about 0.0 as the # phase of the transit minimum so we need to look at stuff from # [0.0, transitphase] and [1.0-transitphase, 1.0] transitphase = thistransduration / 2.0 transitindices = ( ( tphase < transitphase ) | ( tphase > ( 1.0 - transitphase ) ) ) # this is the BLS model # constant = median(tmags) outside transit # constant = thistransitdepth inside transit blsmodel = npfull_like ( tmags , npmedian ( tmags ) ) if magsarefluxes : # eebls.f returns +ve transit depth for fluxes # so we need to subtract here to get fainter fluxes in transit blsmodel [ transitindices ] = ( blsmodel [ transitindices ] - thistransdepth ) else : # eebls.f returns -ve transit depth for magnitudes # so we need to subtract here to get fainter mags in transits blsmodel [ transitindices ] = ( blsmodel [ transitindices ] - thistransdepth ) # see __init__/get_snr_of_dip docstring for description of transit # SNR equation, which is what we use for `thissnr`. subtractedmags = tmags - blsmodel subtractedrms = npstd ( subtractedmags ) npts_in_transit = len ( tmags [ transitindices ] ) thissnr = ( npsqrt ( npts_in_transit ) * npabs ( thistransdepth / subtractedrms ) ) # tell user about stuff if verbose = True if verbose : LOGINFO ( 'refit best period: %.6f, ' 'refit center of transit: %.5f' % ( thisbestperiod , thisminepoch ) ) LOGINFO ( 'transit ingress phase = %.3f to %.3f' % ( 1.0 - transitphase , 1.0 ) ) LOGINFO ( 'transit egress phase = %.3f to %.3f' % ( 0.0 , transitphase ) ) LOGINFO ( 'npoints in transit: %s' % tmags [ transitindices ] . size ) LOGINFO ( 'transit depth (delta): %.5f, ' 'frac transit length (q): %.3f, ' ' SNR: %.3f' % ( thistransdepth , thistransduration , thissnr ) ) return { 'period' : thisbestperiod , 'epoch' : thisminepoch , 'snr' : thissnr , 'transitdepth' : thistransdepth , 'transitduration' : thistransduration , 'nphasebins' : nphasebins , 'transingressbin' : thistransingressbin , 'transegressbin' : thistransegressbin , 'blsmodel' : blsmodel , 'subtractedmags' : subtractedmags , 'phasedmags' : tmags , 'phases' : tphase }
Actually calculates the stats.
Actually calculates the stats .