astrobase/fakelcs/recovery.py
def get_recovered_variables_for_magbin(simbasedir,
                                       magbinmedian,
                                       stetson_stdev_min=2.0,
                                       inveta_stdev_min=2.0,
                                       iqr_stdev_min=2.0,
                                       statsonly=True):
    '''This runs variability selection for the given magbinmedian.

    To generate a full recovery matrix over all magnitude bins, run this
    function for each magbin over the specified stetson_stdev_min and
    inveta_stdev_min grid.

    Parameters
    ----------

    simbasedir : str
        The input directory of fake LCs.

    magbinmedian : float
        The magbin to run the variable recovery for. This is an item from the
        dict from `simbasedir/fakelcs-info.pkl: `fakelcinfo['magrms'][magcol]`
        list for each magcol and designates which magbin to get the recovery
        stats for.

    stetson_stdev_min : float
        The minimum sigma above the trend in the Stetson J variability index
        distribution for this magbin to use to consider objects as variable.

    inveta_stdev_min : float
        The minimum sigma above the trend in the 1/eta variability index
        distribution for this magbin to use to consider objects as variable.

    iqr_stdev_min : float
        The minimum sigma above the trend in the IQR variability index
        distribution for this magbin to use to consider objects as variable.

    statsonly : bool
        If this is True, only the final stats will be returned. If False, the
        full arrays used to generate the stats will also be returned.

    Returns
    -------

    dict
        The returned dict contains statistics for this magbin and if requested,
        the full arrays used to calculate the statistics.

    '''


    # get the info from the simbasedir
    with open(os.path.join(simbasedir, 'fakelcs-info.pkl'),'rb') as infd:
        siminfo = pickle.load(infd)

    objectids = siminfo['objectid']
    varflags = siminfo['isvariable']
    sdssr = siminfo['sdssr']

    # get the column defs for the fakelcs
    timecols = siminfo['timecols']
    magcols = siminfo['magcols']
    errcols = siminfo['errcols']

    # register the fakelc pklc as a custom lcproc format
    # now we should be able to use all lcproc functions correctly
    fakelc_formatkey = 'fake-%s' % siminfo['lcformat']
    lcproc.register_lcformat(
        fakelc_formatkey,
        '*-fakelc.pkl',
        timecols,
        magcols,
        errcols,
        'astrobase.lcproc',
        '_read_pklc',
        magsarefluxes=siminfo['magsarefluxes']
    )

    # make the output directory if it doesn't exit
    outdir = os.path.join(simbasedir, 'recvar-threshold-pkls')
    if not os.path.exists(outdir):
        os.mkdir(outdir)

    # run the variability search
    varfeaturedir = os.path.join(simbasedir, 'varfeatures')
    varthreshinfof = os.path.join(
        outdir,
        'varthresh-magbinmed%.2f-stet%.2f-inveta%.2f.pkl' % (magbinmedian,
                                                             stetson_stdev_min,
                                                             inveta_stdev_min)
    )
    varthresh = varthreshold.variability_threshold(
        varfeaturedir,
        varthreshinfof,
        lcformat=fakelc_formatkey,
        min_stetj_stdev=stetson_stdev_min,
        min_inveta_stdev=inveta_stdev_min,
        min_iqr_stdev=iqr_stdev_min,
        verbose=False
    )

    # get the magbins from the varthresh info
    magbins = varthresh['magbins']

    # get the magbininds
    magbininds = np.digitize(sdssr, magbins)

    # bin the objects according to these magbins
    binned_objectids = []
    binned_actualvars = []
    binned_actualnotvars = []

    # go through all the mag bins and bin up the objectids, actual variables,
    # and actual not-variables
    for mbinind, _magi in zip(np.unique(magbininds),
                              range(len(magbins)-1)):

        thisbinind = np.where(magbininds == mbinind)

        thisbin_objectids = objectids[thisbinind]
        thisbin_varflags = varflags[thisbinind]

        thisbin_actualvars = thisbin_objectids[thisbin_varflags]
        thisbin_actualnotvars = thisbin_objectids[~thisbin_varflags]

        binned_objectids.append(thisbin_objectids)
        binned_actualvars.append(thisbin_actualvars)
        binned_actualnotvars.append(thisbin_actualnotvars)


    # this is the output dict
    recdict = {
        'simbasedir':simbasedir,
        'timecols':timecols,
        'magcols':magcols,
        'errcols':errcols,
        'magsarefluxes':siminfo['magsarefluxes'],
        'stetj_min_stdev':stetson_stdev_min,
        'inveta_min_stdev':inveta_stdev_min,
        'iqr_min_stdev':iqr_stdev_min,
        'magbinmedian':magbinmedian,
    }


    # now, for each magcol, find the magbin corresponding to magbinmedian, and
    # get its stats
    for magcol in magcols:

        # this is the index of the matching magnitude bin for the magbinmedian
        # provided
        magbinind = np.where(
            np.array(varthresh[magcol]['binned_sdssr_median']) == magbinmedian
        )

        magbinind = np.asscalar(magbinind[0])

        # get the objectids, actual vars and actual notvars in this magbin
        thisbin_objectids = binned_objectids[magbinind]
        thisbin_actualvars = binned_actualvars[magbinind]
        thisbin_actualnotvars = binned_actualnotvars[magbinind]

        # stetson recovered variables in this magbin
        stet_recoveredvars = varthresh[magcol][
            'binned_objectids_thresh_stetsonj'
        ][magbinind]

        # calculate TP, FP, TN, FN
        stet_recoverednotvars = np.setdiff1d(thisbin_objectids,
                                             stet_recoveredvars)

        stet_truepositives = np.intersect1d(stet_recoveredvars,
                                            thisbin_actualvars)
        stet_falsepositives = np.intersect1d(stet_recoveredvars,
                                             thisbin_actualnotvars)
        stet_truenegatives = np.intersect1d(stet_recoverednotvars,
                                            thisbin_actualnotvars)
        stet_falsenegatives = np.intersect1d(stet_recoverednotvars,
                                             thisbin_actualvars)

        # calculate stetson recall, precision, Matthews correl coeff
        stet_recall = recall(stet_truepositives.size,
                             stet_falsenegatives.size)

        stet_precision = precision(stet_truepositives.size,
                                   stet_falsepositives.size)

        stet_mcc = matthews_correl_coeff(stet_truepositives.size,
                                         stet_truenegatives.size,
                                         stet_falsepositives.size,
                                         stet_falsenegatives.size)


        # inveta recovered variables in this magbin
        inveta_recoveredvars = varthresh[magcol][
            'binned_objectids_thresh_inveta'
        ][magbinind]
        inveta_recoverednotvars = np.setdiff1d(thisbin_objectids,
                                               inveta_recoveredvars)

        inveta_truepositives = np.intersect1d(inveta_recoveredvars,
                                              thisbin_actualvars)
        inveta_falsepositives = np.intersect1d(inveta_recoveredvars,
                                               thisbin_actualnotvars)
        inveta_truenegatives = np.intersect1d(inveta_recoverednotvars,
                                              thisbin_actualnotvars)
        inveta_falsenegatives = np.intersect1d(inveta_recoverednotvars,
                                               thisbin_actualvars)

        # calculate inveta recall, precision, Matthews correl coeff
        inveta_recall = recall(inveta_truepositives.size,
                               inveta_falsenegatives.size)

        inveta_precision = precision(inveta_truepositives.size,
                                     inveta_falsepositives.size)

        inveta_mcc = matthews_correl_coeff(inveta_truepositives.size,
                                           inveta_truenegatives.size,
                                           inveta_falsepositives.size,
                                           inveta_falsenegatives.size)


        # iqr recovered variables in this magbin
        iqr_recoveredvars = varthresh[magcol][
            'binned_objectids_thresh_iqr'
        ][magbinind]
        iqr_recoverednotvars = np.setdiff1d(thisbin_objectids,
                                            iqr_recoveredvars)

        iqr_truepositives = np.intersect1d(iqr_recoveredvars,
                                           thisbin_actualvars)
        iqr_falsepositives = np.intersect1d(iqr_recoveredvars,
                                            thisbin_actualnotvars)
        iqr_truenegatives = np.intersect1d(iqr_recoverednotvars,
                                           thisbin_actualnotvars)
        iqr_falsenegatives = np.intersect1d(iqr_recoverednotvars,
                                            thisbin_actualvars)

        # calculate iqr recall, precision, Matthews correl coeff
        iqr_recall = recall(iqr_truepositives.size,
                            iqr_falsenegatives.size)

        iqr_precision = precision(iqr_truepositives.size,
                                  iqr_falsepositives.size)

        iqr_mcc = matthews_correl_coeff(iqr_truepositives.size,
                                        iqr_truenegatives.size,
                                        iqr_falsepositives.size,
                                        iqr_falsenegatives.size)


        # calculate the items missed by one method but found by the other
        # methods
        stet_missed_inveta_found = np.setdiff1d(inveta_truepositives,
                                                stet_truepositives)
        stet_missed_iqr_found = np.setdiff1d(iqr_truepositives,
                                             stet_truepositives)

        inveta_missed_stet_found = np.setdiff1d(stet_truepositives,
                                                inveta_truepositives)
        inveta_missed_iqr_found = np.setdiff1d(iqr_truepositives,
                                               inveta_truepositives)

        iqr_missed_stet_found = np.setdiff1d(stet_truepositives,
                                             iqr_truepositives)
        iqr_missed_inveta_found = np.setdiff1d(inveta_truepositives,
                                               iqr_truepositives)

        if not statsonly:

            recdict[magcol] = {
                # stetson J alone
                'stet_recoveredvars':stet_recoveredvars,
                'stet_truepositives':stet_truepositives,
                'stet_falsepositives':stet_falsepositives,
                'stet_truenegatives':stet_truenegatives,
                'stet_falsenegatives':stet_falsenegatives,
                'stet_precision':stet_precision,
                'stet_recall':stet_recall,
                'stet_mcc':stet_mcc,
                # inveta alone
                'inveta_recoveredvars':inveta_recoveredvars,
                'inveta_truepositives':inveta_truepositives,
                'inveta_falsepositives':inveta_falsepositives,
                'inveta_truenegatives':inveta_truenegatives,
                'inveta_falsenegatives':inveta_falsenegatives,
                'inveta_precision':inveta_precision,
                'inveta_recall':inveta_recall,
                'inveta_mcc':inveta_mcc,
                # iqr alone
                'iqr_recoveredvars':iqr_recoveredvars,
                'iqr_truepositives':iqr_truepositives,
                'iqr_falsepositives':iqr_falsepositives,
                'iqr_truenegatives':iqr_truenegatives,
                'iqr_falsenegatives':iqr_falsenegatives,
                'iqr_precision':iqr_precision,
                'iqr_recall':iqr_recall,
                'iqr_mcc':iqr_mcc,
                # true positive variables missed by one method but picked up by
                # the others
                'stet_missed_inveta_found':stet_missed_inveta_found,
                'stet_missed_iqr_found':stet_missed_iqr_found,
                'inveta_missed_stet_found':inveta_missed_stet_found,
                'inveta_missed_iqr_found':inveta_missed_iqr_found,
                'iqr_missed_stet_found':iqr_missed_stet_found,
                'iqr_missed_inveta_found':iqr_missed_inveta_found,
                # bin info
                'actual_variables':thisbin_actualvars,
                'actual_nonvariables':thisbin_actualnotvars,
                'all_objectids':thisbin_objectids,
                'magbinind':magbinind,

            }

        # if statsonly is set, then we only return the numbers but not the
        # arrays themselves
        else:

            recdict[magcol] = {
                # stetson J alone
                'stet_recoveredvars':stet_recoveredvars.size,
                'stet_truepositives':stet_truepositives.size,
                'stet_falsepositives':stet_falsepositives.size,
                'stet_truenegatives':stet_truenegatives.size,
                'stet_falsenegatives':stet_falsenegatives.size,
                'stet_precision':stet_precision,
                'stet_recall':stet_recall,
                'stet_mcc':stet_mcc,
                # inveta alone
                'inveta_recoveredvars':inveta_recoveredvars.size,
                'inveta_truepositives':inveta_truepositives.size,
                'inveta_falsepositives':inveta_falsepositives.size,
                'inveta_truenegatives':inveta_truenegatives.size,
                'inveta_falsenegatives':inveta_falsenegatives.size,
                'inveta_precision':inveta_precision,
                'inveta_recall':inveta_recall,
                'inveta_mcc':inveta_mcc,
                # iqr alone
                'iqr_recoveredvars':iqr_recoveredvars.size,
                'iqr_truepositives':iqr_truepositives.size,
                'iqr_falsepositives':iqr_falsepositives.size,
                'iqr_truenegatives':iqr_truenegatives.size,
                'iqr_falsenegatives':iqr_falsenegatives.size,
                'iqr_precision':iqr_precision,
                'iqr_recall':iqr_recall,
                'iqr_mcc':iqr_mcc,
                # true positive variables missed by one method but picked up by
                # the others
                'stet_missed_inveta_found':stet_missed_inveta_found.size,
                'stet_missed_iqr_found':stet_missed_iqr_found.size,
                'inveta_missed_stet_found':inveta_missed_stet_found.size,
                'inveta_missed_iqr_found':inveta_missed_iqr_found.size,
                'iqr_missed_stet_found':iqr_missed_stet_found.size,
                'iqr_missed_inveta_found':iqr_missed_inveta_found.size,
                # bin info
                'actual_variables':thisbin_actualvars.size,
                'actual_nonvariables':thisbin_actualnotvars.size,
                'all_objectids':thisbin_objectids.size,
                'magbinind':magbinind,
            }


    #
    # done with per magcol
    #

    return recdict
def get_recovered_variables_for_magbin(simbasedir,
                                       magbinmedian,
                                       stetson_stdev_min=2.0,
                                       inveta_stdev_min=2.0,
                                       iqr_stdev_min=2.0,
                                       statsonly=True):
    '''This runs variability selection for the given magbinmedian.

    To generate a full recovery matrix over all magnitude bins, run this
    function for each magbin over the specified stetson_stdev_min and
    inveta_stdev_min grid.

    Parameters
    ----------

    simbasedir : str
        The input directory of fake LCs.

    magbinmedian : float
        The magbin to run the variable recovery for. This is an item from the
        dict from `simbasedir/fakelcs-info.pkl: `fakelcinfo['magrms'][magcol]`
        list for each magcol and designates which magbin to get the recovery
        stats for.

    stetson_stdev_min : float
        The minimum sigma above the trend in the Stetson J variability index
        distribution for this magbin to use to consider objects as variable.

    inveta_stdev_min : float
        The minimum sigma above the trend in the 1/eta variability index
        distribution for this magbin to use to consider objects as variable.

    iqr_stdev_min : float
        The minimum sigma above the trend in the IQR variability index
        distribution for this magbin to use to consider objects as variable.

    statsonly : bool
        If this is True, only the final stats will be returned. If False, the
        full arrays used to generate the stats will also be returned.

    Returns
    -------

    dict
        The returned dict contains statistics for this magbin and if requested,
        the full arrays used to calculate the statistics.

    '''


    # get the info from the simbasedir
    with open(os.path.join(simbasedir, 'fakelcs-info.pkl'),'rb') as infd:
        siminfo = pickle.load(infd)

    objectids = siminfo['objectid']
    varflags = siminfo['isvariable']
    sdssr = siminfo['sdssr']

    # get the column defs for the fakelcs
    timecols = siminfo['timecols']
    magcols = siminfo['magcols']
    errcols = siminfo['errcols']

    # register the fakelc pklc as a custom lcproc format
    # now we should be able to use all lcproc functions correctly
    fakelc_formatkey = 'fake-%s' % siminfo['lcformat']
    lcproc.register_lcformat(
        fakelc_formatkey,
        '*-fakelc.pkl',
        timecols,
        magcols,
        errcols,
        'astrobase.lcproc',
        '_read_pklc',
        magsarefluxes=siminfo['magsarefluxes']
    )

    # make the output directory if it doesn't exit
    outdir = os.path.join(simbasedir, 'recvar-threshold-pkls')
    if not os.path.exists(outdir):
        os.mkdir(outdir)

    # run the variability search
    varfeaturedir = os.path.join(simbasedir, 'varfeatures')
    varthreshinfof = os.path.join(
        outdir,
        'varthresh-magbinmed%.2f-stet%.2f-inveta%.2f.pkl' % (magbinmedian,
                                                             stetson_stdev_min,
                                                             inveta_stdev_min)
    )
    varthresh = varthreshold.variability_threshold(
        varfeaturedir,
        varthreshinfof,
        lcformat=fakelc_formatkey,
        min_stetj_stdev=stetson_stdev_min,
        min_inveta_stdev=inveta_stdev_min,
        min_iqr_stdev=iqr_stdev_min,
        verbose=False
    )

    # get the magbins from the varthresh info
    magbins = varthresh['magbins']

    # get the magbininds
    magbininds = np.digitize(sdssr, magbins)

    # bin the objects according to these magbins
    binned_objectids = []
    binned_actualvars = []
    binned_actualnotvars = []

    # go through all the mag bins and bin up the objectids, actual variables,
    # and actual not-variables
    for mbinind, _magi in zip(np.unique(magbininds),
                              range(len(magbins)-1)):

        thisbinind = np.where(magbininds == mbinind)

        thisbin_objectids = objectids[thisbinind]
        thisbin_varflags = varflags[thisbinind]

        thisbin_actualvars = thisbin_objectids[thisbin_varflags]
        thisbin_actualnotvars = thisbin_objectids[~thisbin_varflags]

        binned_objectids.append(thisbin_objectids)
        binned_actualvars.append(thisbin_actualvars)
        binned_actualnotvars.append(thisbin_actualnotvars)


    # this is the output dict
    recdict = {
        'simbasedir':simbasedir,
        'timecols':timecols,
        'magcols':magcols,
        'errcols':errcols,
        'magsarefluxes':siminfo['magsarefluxes'],
        'stetj_min_stdev':stetson_stdev_min,
        'inveta_min_stdev':inveta_stdev_min,
        'iqr_min_stdev':iqr_stdev_min,
        'magbinmedian':magbinmedian,
    }


    # now, for each magcol, find the magbin corresponding to magbinmedian, and
    # get its stats
    for magcol in magcols:

        # this is the index of the matching magnitude bin for the magbinmedian
        # provided
        magbinind = np.where(
            np.array(varthresh[magcol]['binned_sdssr_median']) == magbinmedian
        )

        magbinind = np.asscalar(magbinind[0])

        # get the objectids, actual vars and actual notvars in this magbin
        thisbin_objectids = binned_objectids[magbinind]
        thisbin_actualvars = binned_actualvars[magbinind]
        thisbin_actualnotvars = binned_actualnotvars[magbinind]

        # stetson recovered variables in this magbin
        stet_recoveredvars = varthresh[magcol][
            'binned_objectids_thresh_stetsonj'
        ][magbinind]

        # calculate TP, FP, TN, FN
        stet_recoverednotvars = np.setdiff1d(thisbin_objectids,
                                             stet_recoveredvars)

        stet_truepositives = np.intersect1d(stet_recoveredvars,
                                            thisbin_actualvars)
        stet_falsepositives = np.intersect1d(stet_recoveredvars,
                                             thisbin_actualnotvars)
        stet_truenegatives = np.intersect1d(stet_recoverednotvars,
                                            thisbin_actualnotvars)
        stet_falsenegatives = np.intersect1d(stet_recoverednotvars,
                                             thisbin_actualvars)

        # calculate stetson recall, precision, Matthews correl coeff
        stet_recall = recall(stet_truepositives.size,
                             stet_falsenegatives.size)

        stet_precision = precision(stet_truepositives.size,
                                   stet_falsepositives.size)

        stet_mcc = matthews_correl_coeff(stet_truepositives.size,
                                         stet_truenegatives.size,
                                         stet_falsepositives.size,
                                         stet_falsenegatives.size)


        # inveta recovered variables in this magbin
        inveta_recoveredvars = varthresh[magcol][
            'binned_objectids_thresh_inveta'
        ][magbinind]
        inveta_recoverednotvars = np.setdiff1d(thisbin_objectids,
                                               inveta_recoveredvars)

        inveta_truepositives = np.intersect1d(inveta_recoveredvars,
                                              thisbin_actualvars)
        inveta_falsepositives = np.intersect1d(inveta_recoveredvars,
                                               thisbin_actualnotvars)
        inveta_truenegatives = np.intersect1d(inveta_recoverednotvars,
                                              thisbin_actualnotvars)
        inveta_falsenegatives = np.intersect1d(inveta_recoverednotvars,
                                               thisbin_actualvars)

        # calculate inveta recall, precision, Matthews correl coeff
        inveta_recall = recall(inveta_truepositives.size,
                               inveta_falsenegatives.size)

        inveta_precision = precision(inveta_truepositives.size,
                                     inveta_falsepositives.size)

        inveta_mcc = matthews_correl_coeff(inveta_truepositives.size,
                                           inveta_truenegatives.size,
                                           inveta_falsepositives.size,
                                           inveta_falsenegatives.size)


        # iqr recovered variables in this magbin
        iqr_recoveredvars = varthresh[magcol][
            'binned_objectids_thresh_iqr'
        ][magbinind]
        iqr_recoverednotvars = np.setdiff1d(thisbin_objectids,
                                            iqr_recoveredvars)

        iqr_truepositives = np.intersect1d(iqr_recoveredvars,
                                           thisbin_actualvars)
        iqr_falsepositives = np.intersect1d(iqr_recoveredvars,
                                            thisbin_actualnotvars)
        iqr_truenegatives = np.intersect1d(iqr_recoverednotvars,
                                           thisbin_actualnotvars)
        iqr_falsenegatives = np.intersect1d(iqr_recoverednotvars,
                                            thisbin_actualvars)

        # calculate iqr recall, precision, Matthews correl coeff
        iqr_recall = recall(iqr_truepositives.size,
                            iqr_falsenegatives.size)

        iqr_precision = precision(iqr_truepositives.size,
                                  iqr_falsepositives.size)

        iqr_mcc = matthews_correl_coeff(iqr_truepositives.size,
                                        iqr_truenegatives.size,
                                        iqr_falsepositives.size,
                                        iqr_falsenegatives.size)


        # calculate the items missed by one method but found by the other
        # methods
        stet_missed_inveta_found = np.setdiff1d(inveta_truepositives,
                                                stet_truepositives)
        stet_missed_iqr_found = np.setdiff1d(iqr_truepositives,
                                             stet_truepositives)

        inveta_missed_stet_found = np.setdiff1d(stet_truepositives,
                                                inveta_truepositives)
        inveta_missed_iqr_found = np.setdiff1d(iqr_truepositives,
                                               inveta_truepositives)

        iqr_missed_stet_found = np.setdiff1d(stet_truepositives,
                                             iqr_truepositives)
        iqr_missed_inveta_found = np.setdiff1d(inveta_truepositives,
                                               iqr_truepositives)

        if not statsonly:

            recdict[magcol] = {
                # stetson J alone
                'stet_recoveredvars':stet_recoveredvars,
                'stet_truepositives':stet_truepositives,
                'stet_falsepositives':stet_falsepositives,
                'stet_truenegatives':stet_truenegatives,
                'stet_falsenegatives':stet_falsenegatives,
                'stet_precision':stet_precision,
                'stet_recall':stet_recall,
                'stet_mcc':stet_mcc,
                # inveta alone
                'inveta_recoveredvars':inveta_recoveredvars,
                'inveta_truepositives':inveta_truepositives,
                'inveta_falsepositives':inveta_falsepositives,
                'inveta_truenegatives':inveta_truenegatives,
                'inveta_falsenegatives':inveta_falsenegatives,
                'inveta_precision':inveta_precision,
                'inveta_recall':inveta_recall,
                'inveta_mcc':inveta_mcc,
                # iqr alone
                'iqr_recoveredvars':iqr_recoveredvars,
                'iqr_truepositives':iqr_truepositives,
                'iqr_falsepositives':iqr_falsepositives,
                'iqr_truenegatives':iqr_truenegatives,
                'iqr_falsenegatives':iqr_falsenegatives,
                'iqr_precision':iqr_precision,
                'iqr_recall':iqr_recall,
                'iqr_mcc':iqr_mcc,
                # true positive variables missed by one method but picked up by
                # the others
                'stet_missed_inveta_found':stet_missed_inveta_found,
                'stet_missed_iqr_found':stet_missed_iqr_found,
                'inveta_missed_stet_found':inveta_missed_stet_found,
                'inveta_missed_iqr_found':inveta_missed_iqr_found,
                'iqr_missed_stet_found':iqr_missed_stet_found,
                'iqr_missed_inveta_found':iqr_missed_inveta_found,
                # bin info
                'actual_variables':thisbin_actualvars,
                'actual_nonvariables':thisbin_actualnotvars,
                'all_objectids':thisbin_objectids,
                'magbinind':magbinind,

            }

        # if statsonly is set, then we only return the numbers but not the
        # arrays themselves
        else:

            recdict[magcol] = {
                # stetson J alone
                'stet_recoveredvars':stet_recoveredvars.size,
                'stet_truepositives':stet_truepositives.size,
                'stet_falsepositives':stet_falsepositives.size,
                'stet_truenegatives':stet_truenegatives.size,
                'stet_falsenegatives':stet_falsenegatives.size,
                'stet_precision':stet_precision,
                'stet_recall':stet_recall,
                'stet_mcc':stet_mcc,
                # inveta alone
                'inveta_recoveredvars':inveta_recoveredvars.size,
                'inveta_truepositives':inveta_truepositives.size,
                'inveta_falsepositives':inveta_falsepositives.size,
                'inveta_truenegatives':inveta_truenegatives.size,
                'inveta_falsenegatives':inveta_falsenegatives.size,
                'inveta_precision':inveta_precision,
                'inveta_recall':inveta_recall,
                'inveta_mcc':inveta_mcc,
                # iqr alone
                'iqr_recoveredvars':iqr_recoveredvars.size,
                'iqr_truepositives':iqr_truepositives.size,
                'iqr_falsepositives':iqr_falsepositives.size,
                'iqr_truenegatives':iqr_truenegatives.size,
                'iqr_falsenegatives':iqr_falsenegatives.size,
                'iqr_precision':iqr_precision,
                'iqr_recall':iqr_recall,
                'iqr_mcc':iqr_mcc,
                # true positive variables missed by one method but picked up by
                # the others
                'stet_missed_inveta_found':stet_missed_inveta_found.size,
                'stet_missed_iqr_found':stet_missed_iqr_found.size,
                'inveta_missed_stet_found':inveta_missed_stet_found.size,
                'inveta_missed_iqr_found':inveta_missed_iqr_found.size,
                'iqr_missed_stet_found':iqr_missed_stet_found.size,
                'iqr_missed_inveta_found':iqr_missed_inveta_found.size,
                # bin info
                'actual_variables':thisbin_actualvars.size,
                'actual_nonvariables':thisbin_actualnotvars.size,
                'all_objectids':thisbin_objectids.size,
                'magbinind':magbinind,
            }


    #
    # done with per magcol
    #

    return recdict
def get_recovered_variables_for_magbin ( simbasedir , magbinmedian , stetson_stdev_min = 2.0 , inveta_stdev_min = 2.0 , iqr_stdev_min = 2.0 , statsonly = True ) : # get the info from the simbasedir with open ( os . path . join ( simbasedir , 'fakelcs-info.pkl' ) , 'rb' ) as infd : siminfo = pickle . load ( infd ) objectids = siminfo [ 'objectid' ] varflags = siminfo [ 'isvariable' ] sdssr = siminfo [ 'sdssr' ] # get the column defs for the fakelcs timecols = siminfo [ 'timecols' ] magcols = siminfo [ 'magcols' ] errcols = siminfo [ 'errcols' ] # register the fakelc pklc as a custom lcproc format # now we should be able to use all lcproc functions correctly fakelc_formatkey = 'fake-%s' % siminfo [ 'lcformat' ] lcproc . register_lcformat ( fakelc_formatkey , '*-fakelc.pkl' , timecols , magcols , errcols , 'astrobase.lcproc' , '_read_pklc' , magsarefluxes = siminfo [ 'magsarefluxes' ] ) # make the output directory if it doesn't exit outdir = os . path . join ( simbasedir , 'recvar-threshold-pkls' ) if not os . path . exists ( outdir ) : os . mkdir ( outdir ) # run the variability search varfeaturedir = os . path . join ( simbasedir , 'varfeatures' ) varthreshinfof = os . path . join ( outdir , 'varthresh-magbinmed%.2f-stet%.2f-inveta%.2f.pkl' % ( magbinmedian , stetson_stdev_min , inveta_stdev_min ) ) varthresh = varthreshold . variability_threshold ( varfeaturedir , varthreshinfof , lcformat = fakelc_formatkey , min_stetj_stdev = stetson_stdev_min , min_inveta_stdev = inveta_stdev_min , min_iqr_stdev = iqr_stdev_min , verbose = False ) # get the magbins from the varthresh info magbins = varthresh [ 'magbins' ] # get the magbininds magbininds = np . digitize ( sdssr , magbins ) # bin the objects according to these magbins binned_objectids = [ ] binned_actualvars = [ ] binned_actualnotvars = [ ] # go through all the mag bins and bin up the objectids, actual variables, # and actual not-variables for mbinind , _magi in zip ( np . unique ( magbininds ) , range ( len ( magbins ) - 1 ) ) : thisbinind = np . where ( magbininds == mbinind ) thisbin_objectids = objectids [ thisbinind ] thisbin_varflags = varflags [ thisbinind ] thisbin_actualvars = thisbin_objectids [ thisbin_varflags ] thisbin_actualnotvars = thisbin_objectids [ ~ thisbin_varflags ] binned_objectids . append ( thisbin_objectids ) binned_actualvars . append ( thisbin_actualvars ) binned_actualnotvars . append ( thisbin_actualnotvars ) # this is the output dict recdict = { 'simbasedir' : simbasedir , 'timecols' : timecols , 'magcols' : magcols , 'errcols' : errcols , 'magsarefluxes' : siminfo [ 'magsarefluxes' ] , 'stetj_min_stdev' : stetson_stdev_min , 'inveta_min_stdev' : inveta_stdev_min , 'iqr_min_stdev' : iqr_stdev_min , 'magbinmedian' : magbinmedian , } # now, for each magcol, find the magbin corresponding to magbinmedian, and # get its stats for magcol in magcols : # this is the index of the matching magnitude bin for the magbinmedian # provided magbinind = np . where ( np . array ( varthresh [ magcol ] [ 'binned_sdssr_median' ] ) == magbinmedian ) magbinind = np . asscalar ( magbinind [ 0 ] ) # get the objectids, actual vars and actual notvars in this magbin thisbin_objectids = binned_objectids [ magbinind ] thisbin_actualvars = binned_actualvars [ magbinind ] thisbin_actualnotvars = binned_actualnotvars [ magbinind ] # stetson recovered variables in this magbin stet_recoveredvars = varthresh [ magcol ] [ 'binned_objectids_thresh_stetsonj' ] [ magbinind ] # calculate TP, FP, TN, FN stet_recoverednotvars = np . setdiff1d ( thisbin_objectids , stet_recoveredvars ) stet_truepositives = np . intersect1d ( stet_recoveredvars , thisbin_actualvars ) stet_falsepositives = np . intersect1d ( stet_recoveredvars , thisbin_actualnotvars ) stet_truenegatives = np . intersect1d ( stet_recoverednotvars , thisbin_actualnotvars ) stet_falsenegatives = np . intersect1d ( stet_recoverednotvars , thisbin_actualvars ) # calculate stetson recall, precision, Matthews correl coeff stet_recall = recall ( stet_truepositives . size , stet_falsenegatives . size ) stet_precision = precision ( stet_truepositives . size , stet_falsepositives . size ) stet_mcc = matthews_correl_coeff ( stet_truepositives . size , stet_truenegatives . size , stet_falsepositives . size , stet_falsenegatives . size ) # inveta recovered variables in this magbin inveta_recoveredvars = varthresh [ magcol ] [ 'binned_objectids_thresh_inveta' ] [ magbinind ] inveta_recoverednotvars = np . setdiff1d ( thisbin_objectids , inveta_recoveredvars ) inveta_truepositives = np . intersect1d ( inveta_recoveredvars , thisbin_actualvars ) inveta_falsepositives = np . intersect1d ( inveta_recoveredvars , thisbin_actualnotvars ) inveta_truenegatives = np . intersect1d ( inveta_recoverednotvars , thisbin_actualnotvars ) inveta_falsenegatives = np . intersect1d ( inveta_recoverednotvars , thisbin_actualvars ) # calculate inveta recall, precision, Matthews correl coeff inveta_recall = recall ( inveta_truepositives . size , inveta_falsenegatives . size ) inveta_precision = precision ( inveta_truepositives . size , inveta_falsepositives . size ) inveta_mcc = matthews_correl_coeff ( inveta_truepositives . size , inveta_truenegatives . size , inveta_falsepositives . size , inveta_falsenegatives . size ) # iqr recovered variables in this magbin iqr_recoveredvars = varthresh [ magcol ] [ 'binned_objectids_thresh_iqr' ] [ magbinind ] iqr_recoverednotvars = np . setdiff1d ( thisbin_objectids , iqr_recoveredvars ) iqr_truepositives = np . intersect1d ( iqr_recoveredvars , thisbin_actualvars ) iqr_falsepositives = np . intersect1d ( iqr_recoveredvars , thisbin_actualnotvars ) iqr_truenegatives = np . intersect1d ( iqr_recoverednotvars , thisbin_actualnotvars ) iqr_falsenegatives = np . intersect1d ( iqr_recoverednotvars , thisbin_actualvars ) # calculate iqr recall, precision, Matthews correl coeff iqr_recall = recall ( iqr_truepositives . size , iqr_falsenegatives . size ) iqr_precision = precision ( iqr_truepositives . size , iqr_falsepositives . size ) iqr_mcc = matthews_correl_coeff ( iqr_truepositives . size , iqr_truenegatives . size , iqr_falsepositives . size , iqr_falsenegatives . size ) # calculate the items missed by one method but found by the other # methods stet_missed_inveta_found = np . setdiff1d ( inveta_truepositives , stet_truepositives ) stet_missed_iqr_found = np . setdiff1d ( iqr_truepositives , stet_truepositives ) inveta_missed_stet_found = np . setdiff1d ( stet_truepositives , inveta_truepositives ) inveta_missed_iqr_found = np . setdiff1d ( iqr_truepositives , inveta_truepositives ) iqr_missed_stet_found = np . setdiff1d ( stet_truepositives , iqr_truepositives ) iqr_missed_inveta_found = np . setdiff1d ( inveta_truepositives , iqr_truepositives ) if not statsonly : recdict [ magcol ] = { # stetson J alone 'stet_recoveredvars' : stet_recoveredvars , 'stet_truepositives' : stet_truepositives , 'stet_falsepositives' : stet_falsepositives , 'stet_truenegatives' : stet_truenegatives , 'stet_falsenegatives' : stet_falsenegatives , 'stet_precision' : stet_precision , 'stet_recall' : stet_recall , 'stet_mcc' : stet_mcc , # inveta alone 'inveta_recoveredvars' : inveta_recoveredvars , 'inveta_truepositives' : inveta_truepositives , 'inveta_falsepositives' : inveta_falsepositives , 'inveta_truenegatives' : inveta_truenegatives , 'inveta_falsenegatives' : inveta_falsenegatives , 'inveta_precision' : inveta_precision , 'inveta_recall' : inveta_recall , 'inveta_mcc' : inveta_mcc , # iqr alone 'iqr_recoveredvars' : iqr_recoveredvars , 'iqr_truepositives' : iqr_truepositives , 'iqr_falsepositives' : iqr_falsepositives , 'iqr_truenegatives' : iqr_truenegatives , 'iqr_falsenegatives' : iqr_falsenegatives , 'iqr_precision' : iqr_precision , 'iqr_recall' : iqr_recall , 'iqr_mcc' : iqr_mcc , # true positive variables missed by one method but picked up by # the others 'stet_missed_inveta_found' : stet_missed_inveta_found , 'stet_missed_iqr_found' : stet_missed_iqr_found , 'inveta_missed_stet_found' : inveta_missed_stet_found , 'inveta_missed_iqr_found' : inveta_missed_iqr_found , 'iqr_missed_stet_found' : iqr_missed_stet_found , 'iqr_missed_inveta_found' : iqr_missed_inveta_found , # bin info 'actual_variables' : thisbin_actualvars , 'actual_nonvariables' : thisbin_actualnotvars , 'all_objectids' : thisbin_objectids , 'magbinind' : magbinind , } # if statsonly is set, then we only return the numbers but not the # arrays themselves else : recdict [ magcol ] = { # stetson J alone 'stet_recoveredvars' : stet_recoveredvars . size , 'stet_truepositives' : stet_truepositives . size , 'stet_falsepositives' : stet_falsepositives . size , 'stet_truenegatives' : stet_truenegatives . size , 'stet_falsenegatives' : stet_falsenegatives . size , 'stet_precision' : stet_precision , 'stet_recall' : stet_recall , 'stet_mcc' : stet_mcc , # inveta alone 'inveta_recoveredvars' : inveta_recoveredvars . size , 'inveta_truepositives' : inveta_truepositives . size , 'inveta_falsepositives' : inveta_falsepositives . size , 'inveta_truenegatives' : inveta_truenegatives . size , 'inveta_falsenegatives' : inveta_falsenegatives . size , 'inveta_precision' : inveta_precision , 'inveta_recall' : inveta_recall , 'inveta_mcc' : inveta_mcc , # iqr alone 'iqr_recoveredvars' : iqr_recoveredvars . size , 'iqr_truepositives' : iqr_truepositives . size , 'iqr_falsepositives' : iqr_falsepositives . size , 'iqr_truenegatives' : iqr_truenegatives . size , 'iqr_falsenegatives' : iqr_falsenegatives . size , 'iqr_precision' : iqr_precision , 'iqr_recall' : iqr_recall , 'iqr_mcc' : iqr_mcc , # true positive variables missed by one method but picked up by # the others 'stet_missed_inveta_found' : stet_missed_inveta_found . size , 'stet_missed_iqr_found' : stet_missed_iqr_found . size , 'inveta_missed_stet_found' : inveta_missed_stet_found . size , 'inveta_missed_iqr_found' : inveta_missed_iqr_found . size , 'iqr_missed_stet_found' : iqr_missed_stet_found . size , 'iqr_missed_inveta_found' : iqr_missed_inveta_found . size , # bin info 'actual_variables' : thisbin_actualvars . size , 'actual_nonvariables' : thisbin_actualnotvars . size , 'all_objectids' : thisbin_objectids . size , 'magbinind' : magbinind , } # # done with per magcol # return recdict
This runs variability selection for the given magbinmedian.

    To generate a full recovery matrix over all magnitude bins, run this
    function for each magbin over the specified stetson_stdev_min and
    inveta_stdev_min grid.

    Parameters
    ----------

    simbasedir : str
        The input directory of fake LCs.

    magbinmedian : float
        The magbin to run the variable recovery for. This is an item from the
        dict from `simbasedir/fakelcs-info.pkl: `fakelcinfo['magrms'][magcol]`
        list for each magcol and designates which magbin to get the recovery
        stats for.

    stetson_stdev_min : float
        The minimum sigma above the trend in the Stetson J variability index
        distribution for this magbin to use to consider objects as variable.

    inveta_stdev_min : float
        The minimum sigma above the trend in the 1/eta variability index
        distribution for this magbin to use to consider objects as variable.

    iqr_stdev_min : float
        The minimum sigma above the trend in the IQR variability index
        distribution for this magbin to use to consider objects as variable.

    statsonly : bool
        If this is True, only the final stats will be returned. If False, the
        full arrays used to generate the stats will also be returned.

    Returns
    -------

    dict
        The returned dict contains statistics for this magbin and if requested,
        the full arrays used to calculate the statistics.
This runs variability selection for the given magbinmedian .