text:
#!/usr/bin/python
#demography.py

'''
Functions to implement demographic operations (birth and death).
'''

#geonomics imports
from geonomics.ops.selection import _calc_prob_death
from geonomics.utils.viz import _check_display

#other imports
import numpy as np
import matplotlib as mpl
_check_display()
import matplotlib.pyplot as plt
from numpy import random as r
from scipy import interpolate
from sklearn.preprocessing import normalize
from collections import OrderedDict as OD
from operator import itemgetter

######################################
# -----------------------------------#
# CLASSES ---------------------------#
# -----------------------------------#
######################################

class _DebugPlotter:
    def __init__(self, land, spp, timestep):
        self.subplot = 241
        self.timestep = timestep
        self.spp = spp
        self.land = land

    def _next_plot(self, varname, var):
        if self.subplot <= 248:
            if not plt.get_fignums():
                fig = plt.figure()
                plt.suptitle('TIMESTEP: %i' % self.timestep)
            else:
                fig = plt.gcf()
            plt.subplot(self.subplot)
            plt.title(varname)
            self.spp._plot(land = self.land, lyr_num = 0, size = 8, hide_land = True)
            plt.imshow(var, interpolation = 'nearest', cmap = 'terrain')
            plt.colorbar()
            self.subplot += 1
        else:
            pass

######################################
# -----------------------------------#
# FUNCTIONS -------------------------#
# -----------------------------------#
######################################

def _calc_n_pairs(pairs, spp):
    #if there are pairs
    if pairs.size > 0:
        #get their coordinates
        pairs_coords = spp._get_coords(individs = pairs.flatten())
        #take just the x and y coordinates from pairs_coords, reshape them to
        #match pairs.shape (i.e.nx2, where n = len(spp)), and then mean across
        #axis 1 (to get the pair's midpoint
        #x and y coordinates)
        p_x = pairs_coords[:,0].reshape(pairs.shape).mean(axis = 1)
        p_y = pairs_coords[:,1].reshape(pairs.shape).mean(axis = 1)

        #use the spp._dens_grids to calculate a raster of the expected number
        #of pairs
        #NOTE: because the window_width of that object is much wider than
        #dispersal_mu in the parameterizations I've typically been running,
        #this produces much wider areas of high values in the N_b, N_d, and
        #hence d rasters than the old method did (which had its window_width
        #set to dispersal_mu, leading of course to huge lags in the old
        #calc_density algorithm) need to think about whether I need to attend
        #to this, and if so, how
        n_pairs = np.clip(spp._dens_grids._calc_density(p_x, p_y), a_min = 0,
                                                                a_max = None)

        #remove NaNs
        n_pairs[np.isnan(n_pairs)] = 0

    #otherwise just return a 0 raster
    else:
        n_pairs = np.zeros(spp._land_dim)

    return n_pairs

#the logistic eqxn
def _calc_logistic_growth(R, N, K):
    dNdt = R*(1-(N/K))*N
    return dNdt

def _calc_logistic_soln(x0, R, t):
        return 1/(1+((1/x0)-1)*np.e**(-1*R*t))

def _calc_dNdt(R, N, K, pop_growth_eq = 'logistic'):
    with np.errstate(divide='ignore', invalid='ignore'):
        #if using the logistic equation for pop growth
        #NOTE: For now this is the only option and the default, but could
        #easily offer other equations later if desired
        if pop_growth_eq == 'logistic':
            #use logistic eqxn, with spp intrinsic growth rate (spp.R)
            #to generate current growth-rate raster
            dNdt = _calc_logistic_growth(R, N, K)
    #coerce extreme negative values, NaNs, and infs to a reasonable value
    #(extreme negative values are produced in areas where the K raster
    #has extremely small values)
    dNdt = np.clip(dNdt, a_min = -1*N.max(), a_max = None)
    dNdt[np.isnan(dNdt)] = -1*N.max()
    dNdt[np.isinf(dNdt)] = -1*N.max()
    return dNdt

def _calc_N_b(b, n_births_distr_lambda, n_pairs):
    #Use n_pairs, spp.b (i.e. probability that a pair gives birth),
    #and spp.n_births_distr_lambda (i.e. expected number of births
    #per pair) to calculate the N_b raster (expected number of births in each
    #cell) according to N_b = b*E[N_b/pair]*n_pairs where E[N_b/pair] =
    #n_births_distr_lambda (b/c births Poission distributed)
    #NOTE: Would be good to in the future to use neg_binom dist to allow
    #births to be modeled as Poisson (using correct parameterization of
    #neg_binom to approximate that) or overdispersed (to be parmeterized
    #with field data)

    #NOTE: Does it matter that I'm not attempting to account for the fact
    #that pairs will disperse out of the cell? Almost certainly a larger effect
    #with increasing dispersal-kernel distance. But nonetheless, would it even
    #make sense (and even BE TRACTABLE) to calculate it any other way? My
    #intuition is NO, but worth considering further later on...

    #NOTE: Currently, birth rate is intrinsic to the species and cannot be
    #made a function of density or spatial fitness. Might be worth adding
    #this at some point in the future though.
    N_b = b * n_births_distr_lambda * n_pairs
    return N_b

def _calc_Nd(N_b, dNdt):
    #Use N_b - dN/dt to calculate N_d raster (expected number of deaths in
    #each cell)
    N_d = N_b - dNdt
    return N_d

def _calc_d(N_d, N, d_min, d_max):
    #Use N_d/N to calculate d raster (probability of death of individuals
    #in each cell) Calculate d as simply N_d/N, then control for some
    #numerical artefacts that an inappropriate
    # NOTE: ignore x/0. and 0./0. warnings
    with np.errstate(divide='ignore', invalid='ignore'):
        d = N_d/N
    #fix infinties and NaNs and negatives if they arise
    #(they occur where N ==0 or where N_d == 0)
    d[np.isnan(d)] = 0
    #constrain to the min and max d values
    d = np.clip(a = d, a_min = d_min, a_max = d_max)

    #NOTE: DEH: 12-08-18: Got rid of the infinity check, because it was
    #coercing all infinities to 0, but in many cases infinities should actually
    #have been 1s, so it created numerous problems that just clipping to d_min
    #and d_max takes care of anyhow; commented out the NaN check below, but not
    #entirely clear that I don't need it...
    #d[np.isnan(d)] = 0
    return d

def _do_mortality(spp, death_probs):
    deaths = np.array([*spp])[np.bool8(r.binomial(n = 1, p = death_probs))]
    if len(deaths) > 0:
        ig = itemgetter(*deaths)
        [spp.pop(ind) for ind in deaths];
    return len(deaths)

def _do_pop_dynamics(spp, land, with_selection = True, burn = False,
    births_before_deaths = False, asserts = True, debug = False):
    '''Generalized function for implementation population dynamics.
    Will carry out one round of mating and death, according to parameterization
    laid out in params dict (which were grabbed as Species attributes).

       If with_selection == False, only density-dependent death will occur.

       If burn == True, selection will be coerced to False, the Species.K
       attribute will be updated, stationarity metrics will be assessed, and
       the function will return the decision it has reached regarding
       stationarity.

    '''

    #NOTE: Would be good to write in option to ONLY mate or ONLY die when the
    #function is run (so that a spp could be run with mating every x timesteps
    #but mortality each timestep, for example)

    #create the debug-plotter, id debug == True
    if debug:
        timestep = -99999
        dp = _DebugPlotter(land, spp, timestep)

    #find mating pairs
    pairs = spp._find_mating_pairs()

    #calc num_pairs raster (use the calc_pop_density function on the
    #midpoints of the mating pairs)
    n_pairs = _calc_n_pairs(pairs = pairs, spp = spp)
    #run checks on n_pairs
    if asserts:
        assert n_pairs.min() >= 0, 'n_pairs.min() == %0.2f' %(n_pairs.min())
        assert not np.any(np.isnan(n_pairs))
        assert not np.any(np.isinf(n_pairs))
    #add debug plot
    if debug:
        dp._next_plot('n_pairs', n_pairs)

    #if births should happen before (and thus be included in the calculation of)
    #deaths, then mate and disperse babies now
    if births_before_deaths:
        #Feed the land and mating pairs to spp.do_mating, to produce and
        #disperse zygotes
        spp._do_mating(land, pairs, burn)

    #calc N raster, set it as spp.N, then get it
    spp._calc_density(set_N = True)
    N = spp.N
    #run checks on N
    if asserts:
        assert N.min() >= 0
        assert not np.any(np.isnan(N))
        assert not np.any(np.isinf(N))
    #add debug plot
    if debug:
        dp._next_plot('N', N)

    #get K raster
    K = spp.K
    #run checks on K
    if asserts:
        assert K.min() >= 0
        assert not np.any(np.isnan(K))
        assert not np.any(np.isinf(K))
    #add debug plot
    if debug:
        dp._next_plot('K', K)

    #calc dNdt
    dNdt = _calc_dNdt(R = spp.R, N = N, K = K,
                                    pop_growth_eq = 'logistic')
    #run checks on dNdt
    if asserts:
        assert not np.any(np.isnan(dNdt))
        assert not np.any(np.isinf(dNdt)), ('The following cells are '
            'infinite: \n\t%s') % str([i for i, n in enumerate(
                                        dNdt.ravel()) if np.isinf(n)])
    #add debug plot
    if debug:
        dp._next_plot('dNdt', dNdt)

    #calculate N_b (raster of estimated births)
    N_b = _calc_N_b(b = spp.b,
        n_births_distr_lambda = spp.n_births_distr_lambda, n_pairs = n_pairs)
    #run checks on N_b
    if asserts:
        assert N_b.min() >= 0
        assert not np.any(np.isnan(N_b))
        assert not np.any(np.isinf(N_b))
    #add debug plot
    if debug:
        dp._next_plot('N_b', N_b)

    #calc N_d (raster of deaths)
    N_d = _calc_Nd(N_b = N_b, dNdt = dNdt)
    #run checks on N_d
    if asserts:
        assert not np.any(np.isnan(N_d))
        assert not np.any(np.isinf(N_d))
    #add debug plot
    if debug:
        dp._next_plot('N_d', N_d)

    #calc d (raster of probabilities of density-dependent death)
    d = _calc_d(N_d = N_d, N = N, d_min = spp.d_min, d_max = spp.d_max)

    #run checks on d
    if asserts:
        assert d.min() >= 0, 'd.min() is %0.2f, at %s' % (d.min(),
                                                        str(d.argmin()))
        assert d.max() <= 1, 'd.max() is %0.2f' % d.max()
        assert not np.any(np.isnan(d))
        assert not np.any(np.isinf(d))
    #add debug plot
    if debug:
        dp._next_plot('d', d)

    #If births should happen after (and thus not be included in the
    #calculation of) deaths, then instead of having mated and dispersed
    #babies up above, do it now (i.e. now that the d raster has been
    #calculated)
    #NOTE: The way I'm looking at it now, calculating likely births from number
    #of existing pairs, then calculating death raster, then having births take
    #place, then having per-individual death probabilies determined and having
    #deaths take place (such that newborns could also immediately die after
    #birth, based on the death raster calculated before they were born, and
    #even based on individs who weren't their parents, assuming they dispersed
    #across cell lines), which is what I'm doing currently, is the only
    #obvious/reasonable way to do this. Of course there will always be an
    #inherent rub between reality and discrete-time simulations because of
    #the continuous/discrete time conflict. But in this case, is there any
    #strong reason to CHANGE THIS? Or are there any foreseeable and undesirable
    #results/effects?... NEED TO PUT MORE THOUGHT INTO THIS LATER.
    if not births_before_deaths:
        #Feed the land and mating pairs to the mating functions, to produce
        #and disperse zygotes
        spp._do_mating(land, pairs, burn)
    #Get death probabilities
    death_probs = d[spp.cells[:,1], spp.cells[:,0]]
    #If with_selection (i.e. if death probs should account for fitness),
    #then use the d raster and individuals' fitnesses to calculate
    #per-individual probabilities of death
    if with_selection:
        death_probs = _calc_prob_death(spp, death_probs)
    #run checks on death_probs
    if asserts:
        assert np.alltrue(death_probs >= 0)
        assert np.alltrue(death_probs <= 1)

    #kill (and track kills) based on max_age
    num_killed_age = 0
    if spp.max_age is not None:
        death_probs[spp._get_age() > spp.max_age] = 1
        num_killed_age = np.sum(death_probs == 1)

    #Use the per-individual death probabilities to carry out mortality
    num_deaths = _do_mortality(spp, death_probs)
    spp._set_coords_and_cells()
    spp.n_deaths.append(num_deaths)

    #Check if extinct, and return result
    extinct = spp._check_extinct()
    return(extinct)

reason: import numpy,from numpy,from scipy


text:
from database import COLMAPDatabase, blob_to_array, array_to_blob
from read_write_model import read_cameras_binary, read_images_binary
import argparse, os
import numpy as np

def main(args):
    destionation_db = COLMAPDatabase.connect(args.destionation)
    c = destionation_db.cursor()
    images = read_images_binary(os.path.join(args.source,'images.bin'))
    cameras = read_cameras_binary(os.path.join(args.source,'cameras.bin'))
    lookup_image = dict([(images[i][4],images[i][0]) for i in images])
    c.execute('SELECT image_id,name,camera_id FROM images')
    image_list = c.fetchall()
    image_reorder = []
    camera_lookup = {}
    camera_reorder = []
    #image reorder
    for db_img in image_list:
        img_id, qvec, tvec, cam_id, name, xys, point3D_ids = images[lookup_image[db_img[1]]]
        new_image = [img_id, qvec, tvec, cam_id, name, xys, point3D_ids]
        new_image[0] = db_img[0]
        if new_image[3] != db_img[2]:
            camera_lookup[new_image[3]] = db_img[2]
            new_image[3] = db_img[2]
        image_reorder.append(new_image)
    #camera reorder
    for i in cameras:
        cam_id, model, width, height, params = cameras[i]
        new_camera = [cam_id, model, width, height, params]
        if i in camera_lookup:
            new_camera[0] = camera_lookup[i]
        camera_reorder.append(new_camera)
    destionation_db.execute('DELETE FROM images')
    destionation_db.execute('DELETE FROM cameras')
    for i in range(len(image_reorder)):
        img_id, qvec, tvec, cam_id, name, xys, point3D_ids = image_reorder[i]
        destionation_db.add_image( name, cam_id, prior_q=qvec, prior_t=tvec, image_id=img_id)
    for i in range(len(camera_reorder)):
        cam_id, model, width, height, params = camera_reorder[i]
        known_focal = params[0]
        destionation_db.add_camera(model, width, height, params[:-1], prior_focal_length=known_focal, camera_id=cam_id)
    destionation_db.commit()
    destionation_db.close()
    return None

if __name__ == '__main__':
    parser = argparse.ArgumentParser(description='colmap_copy.py - Colmap camera intrisic & Extrinsic copyer')
    parser.add_argument('-s', '--source', default='penguinguy', type=str,
        help='sparse model directory that already do structure from motion to be source file')
    parser.add_argument('-d', '--destionation', default='pen_matched_400img.db', type=str,
        help='database file from feature mapping to set intrisic and extrinsic')
    main(parser.parse_args())

ext: py

max_stars_repo_name: pureexe/colmap-copy-camera

max_issues_repo_name: pureexe/colmap-copy-camera

max_forks_repo_name: pureexe/colmap-copy-camera

reason: import numpy


text:
#!/usr/bin/env python
# encoding:utf-8
"""
Author : Yuanqing Mei
Date : 2021/4/8
Time: 15:42
File: RocMethod.py
HomePage : http://github.com/yuanqingmei
Email : dg1533019@smail.nju.edu.cn

This script find out the cutoff of a metric value by maximizing the AUC value and ROC、BPP、MFM、GM methods.

References：
[1]  Bender, R. Quantitative risk assessment in epidemiological studies investigating threshold effects.
     Biometrical Journal, 41 (1999), 305-319.（计算VARL的SE（标准误）的参考文献P310）
[2]  Zhou, Y., et al. "An in-depth study of the potentially confounding effect of class size in fault prediction."
     ACM Trans. Softw. Eng. Methodol. (2014) 23(1): 1-51. (计算BPP、MFM(F1)值为阈值)
[3]  Shatnawi, R. (2018). Identifying Threshold Values of Change-Prone Modules.
     (计算sum(Sensitivity+Specificity)=sum(TPR+TNR)值为阈值)
"""
import time

def roc_threshold(working_dir="F:\\NJU\\MTmeta\\experiments\\supervised\\trainingData\\",
                  result_dir="F:\\NJU\\MTmeta\\experiments\\supervised\\",
                  training_list="List.txt"):

    import os
    import csv
    import numpy as np
    import pandas as pd
    import statsmodels.api as sm
    # from sklearn import metrics
    from sklearn.metrics import recall_score, precision_score, f1_score, roc_curve, auc, roc_auc_score, confusion_matrix

    # 显示所有列
    pd.set_option('display.max_columns', None)

    # 显示所有行
    pd.set_option('display.max_rows', None)

    # the item of row of dataframe
    pd.set_option('display.width', 5000)

    working_directory = working_dir
    result_directory = result_dir
    os.chdir(working_directory)

    with open(working_dir + training_list) as l:
        lines = l.readlines()

    for line in lines:
        file = line.replace("\n", "")
        print('the file is ', file)

        # 分别处理每一个项目: f1取出要被处理的项目;
        #                  f2:用于存储每一个项目的阈值信息,f2用csv.writer写数据时没有newline参数，会多出一空行;
        #                  deletedList: 用于存储项目中某个度量logistic回归时，系数不显著或系数为零的度量及该项目名
        with open(working_directory + file, 'r', encoding="ISO-8859-1") as f1, \
                open(result_directory + "RocThreshold\\ROC_Thresholds.csv", 'a+', encoding="utf-8", newline='') as f2, \
                open(result_directory + "RocThreshold\\deletedList.csv", 'a+', encoding="utf-8") as deletedList:

            reader = csv.reader(f1)
            writer = csv.writer(f2)
            writer_deletedList = csv.writer(deletedList)
            # receives the first line of a file and convert to dict generator
            fieldnames = next(reader)
            # exclude the non metric fields (12 items) and metric values including undef and undefined (17 items)
            non_metric = ["relName", "className", "bug"]

            # metric_data stores the metric fields (102 items)
            def fun_1(m):
                return m if m not in non_metric else None

            metric_data = filter(fun_1, fieldnames)

            df = pd.read_csv(file)
            # drop all rows that have any NaN values,删除表中含有任何NaN的行,并重新设置行号
            df = df.dropna(axis=0, how='any', inplace=False).reset_index(drop=True)

            if os.path.getsize(result_directory + "RocThreshold\\ROC_Thresholds.csv") == 0:
                writer.writerow(["fileName", "metric", "Corr_metric_bug", "B_0", "B_0_pValue", "B_1", "B_1_pValue",
                                 "cov11", "cov12", "cov22", "BaseProbability_1",
                                 "auc_threshold", "auc_threshold_variance", "auc_max_value", "i_auc_max",
                                 "gm_threshold", "gm_threshold_variance", "gm_max_value", "i_gm_max",
                                 "bpp_threshold", "bpp_threshold_variance", "bpp_max_value", "i_bpp_max",
                                 "mfm_threshold", "mfm_threshold_variance", "f1_max_value", "i_f1_max",
                                 "roc_threshold", "roc_threshold_variance", "roc_max_value", "i_roc_max",
                                 "varl_threshold", "varl_threshold_variance"])

            if os.path.getsize(result_directory + "RocThreshold\\deletedList.csv") == 0:
                writer_deletedList.writerow(["fileName", "metric", "B_0_pValue", "B_0",
                                             "auc_max_value", "i_auc_max", "gm_max_value", "i_gm_max", "bpp_max_value",
                                             "i_bpp_max", "f1_max_value", "i_f1_max", "roc_max_value", "i_roc_max"])

            for metric in metric_data:

                print("the current file is ", file)
                print("the current metric is ", metric)

                # 由于bug中存储的是缺陷个数,转化为二进制存储,若x>2,则可预测bug为3个以上的阈值,其他类推
                df['bugBinary'] = df.bug.apply(lambda x: 1 if x > 0 else 0)

                # 依次用该度量的每一个值作为阈值计算出auc和GM,然后选择auc最大值的那个度量值作为阈值，即断点回归的cutoff
                # 同时计算BPP(Balanced-pf-pd)、MFM(F1)和ROC(Sensitivity+Specificity)=(TPR+TNR)值，
                # 分别定义存入五个值list,最大值和取最大值的下标值
                AUCs = []
                GMs = []
                BPPs = []
                MFMs = []
                ROCs = []

                auc_max_value = 0
                gm_max_value = 0
                bpp_max_value = 0
                f1_max_value = 0
                roc_max_value = 0

                i_auc_max = 0
                i_gm_max = 0
                i_bpp_max = 0
                i_f1_max = 0
                i_roc_max = 0

                # 判断每个度量与bug之间的关系,因为该关系会影响到断点回归时,相关系数大于零,则LATE估计值大于零,反之,则LATE估计值小于零
                Corr_metric_bug = df.loc[:, [metric, 'bug']].corr('spearman')

                # the i value in this loop, is the subscript value in the list of AUCs, GMs etc.
                for i in range(len(df)):

                    t = df.loc[i, metric]
                    if Corr_metric_bug[metric][1] < 0:
                        df['predictBinary'] = df[metric].apply(lambda x: 1 if x <= t else 0)
                    else:
                        df['predictBinary'] = df[metric].apply(lambda x: 1 if x >= t else 0)

                    # confusion_matrix()函数中需要给出label, 0和1，否则该函数算不出TP,因为不知道哪个标签是poistive.
                    c_matrix = confusion_matrix(df["bugBinary"], df['predictBinary'], labels=[0, 1])
                    tn, fp, fn, tp = c_matrix.ravel()

                    if (tn + fp) == 0:
                        tnr_value = 0
                    else:
                        tnr_value = tn / (tn + fp)

                    if (fp + tn) == 0:
                        fpr = 0
                    else:
                        fpr = fp / (fp + tn)
                    # fpr, tpr, thresholds = roc_curve(df['bugBinary'], df['predictBinary'])
                    # AUC = auc(fpr, tpr)
                    auc_value = roc_auc_score(df['bugBinary'], df['predictBinary'])
                    recall_value = recall_score(df['bugBinary'], df['predictBinary'], labels=[0, 1])
                    precision_value = precision_score(df['bugBinary'], df['predictBinary'], labels=[0, 1])
                    f1_value = f1_score(df['bugBinary'], df['predictBinary'], labels=[0, 1])

                    gm_value = (recall_value * tnr_value) ** 0.5
                    pfr = recall_value
                    pdr = fpr  # fp / (fp + tn)
                    bpp_value = 1 - (((0 - pfr) ** 2 + (1 - pdr) ** 2) * 0.5) ** 0.5
                    roc_value = recall_value + tnr_value

                    AUCs.append(auc_value)
                    GMs.append(gm_value)
                    BPPs.append(bpp_value)
                    MFMs.append(f1_value)
                    ROCs.append(roc_value)

                    # 求出上述五个list中最大值，及对应的i值，可能会有几个值相同，且为最大值，则取第一次找到那个值(i)为阈值
                    if auc_value > auc_max_value:
                        auc_max_value = auc_value
                        i_auc_max = i

                    if gm_value > gm_max_value:
                        gm_max_value = gm_value
                        i_gm_max = i

                    if bpp_value > bpp_max_value:
                        bpp_max_value = bpp_value
                        i_bpp_max = i

                    if f1_value > f1_max_value:
                        f1_max_value = f1_value
                        i_f1_max = i

                    if roc_value > roc_max_value:
                        roc_max_value = roc_value
                        i_roc_max = i

                print("auc_max_value is ", auc_max_value)
                print("gm_max_value is ", gm_max_value)
                print("bpp_max_value is ", bpp_max_value)
                print("f1_max_value is ", f1_max_value)
                print("roc_max_value is ", roc_max_value)

                print("i_auc_max is ", i_auc_max)
                print("i_gm_max is ", i_gm_max)
                print("i_bpp_max is ", i_bpp_max)
                print("i_f1_max is ", i_f1_max)
                print("i_roc_max is ", i_roc_max)

                df['intercept'] = 1.0

                # 通过 statsmodels.api 逻辑回归分类; 指定作为训练变量的列，不含目标列`bug`
                logit = sm.Logit(df['bugBinary'], df.loc[:, [metric, 'intercept']])
                # 拟合模型,disp=1 用于显示结果
                result = logit.fit(method='bfgs', disp=0)
                print(result.summary())

                pValueLogit = result.pvalues
                if pValueLogit[0] > 0.05:  # 自变量前的系数
                    writer_deletedList.writerow(
                        [file, metric, pValueLogit[0], B[0], auc_max_value, i_auc_max, gm_max_value,
                         i_gm_max, bpp_max_value, i_bpp_max, f1_max_value, i_f1_max, roc_max_value,
                         i_roc_max])
                    continue  # 若训练数据LOGIT回归系数的P值大于0.05，放弃该数据。

                B = result.params  # logit回归系数
                if B[0] == 0:  # 自变量前的系数
                    writer_deletedList.writerow(
                        [file, metric, pValueLogit[0], B[0], auc_max_value, i_auc_max, gm_max_value,
                         i_gm_max, bpp_max_value, i_bpp_max, f1_max_value, i_f1_max, roc_max_value,
                         i_roc_max])
                    continue  # 若训练数据LOGIT回归系数等于0，放弃该数据。

                # 计算auc阈值及标准差,包括其他四个类型阈值
                auc_threshold = df.loc[i_auc_max, metric]
                gm_threshold = df.loc[i_gm_max, metric]
                bpp_threshold = df.loc[i_bpp_max, metric]
                mfm_threshold = df.loc[i_f1_max, metric]
                roc_threshold = df.loc[i_roc_max, metric]
                # 计算LOGIT回归系数矩阵的协方差矩阵,因为计算aucThreshold的标准差要用到,见参考文献[1],
                # 此处借鉴VARL方法，本质上VARL也是度量值中的一个
                cov = result.cov_params()
                cov11 = cov.iloc[0, 0]
                cov12 = cov.iloc[0, 1]
                cov22 = cov.iloc[1, 1]
                auc_threshold_se = ((cov.iloc[0, 0] + 2 * auc_threshold * cov.iloc[0, 1]
                                     + auc_threshold * auc_threshold * cov.iloc[1, 1]) ** 0.5) / B[0]
                auc_threshold_variance = auc_threshold_se ** 2

                gm_threshold_se = ((cov.iloc[0, 0] + 2 * gm_threshold * cov.iloc[0, 1]
                                    + gm_threshold * gm_threshold * cov.iloc[1, 1]) ** 0.5) / B[0]
                gm_threshold_variance = gm_threshold_se ** 2

                bpp_threshold_se = ((cov.iloc[0, 0] + 2 * bpp_threshold * cov.iloc[0, 1]
                                     + bpp_threshold * bpp_threshold * cov.iloc[1, 1]) ** 0.5) / B[0]
                bpp_threshold_variance = bpp_threshold_se ** 2

                mfm_threshold_se = ((cov.iloc[0, 0] + 2 * mfm_threshold * cov.iloc[0, 1]
                                     + mfm_threshold * mfm_threshold * cov.iloc[1, 1]) ** 0.5) / B[0]
                mfm_threshold_variance = mfm_threshold_se ** 2

                roc_threshold_se = ((cov.iloc[0, 0] + 2 * roc_threshold * cov.iloc[0, 1]
                                     + roc_threshold * roc_threshold * cov.iloc[1, 1]) ** 0.5) / B[0]
                roc_threshold_variance = roc_threshold_se ** 2

                # 求VARL作为阈值，此处未用10折交叉验证的方法 VARL.threshold = (log(Porbability[1]/Porbability[2])-B[1])/B[2]
                valueOfbugBinary = df["bugBinary"].value_counts()  # 0 和 1 的各自的个数
                print("the value of valueOfbugBinary[0] is ", valueOfbugBinary[0])
                print("the value of valueOfbugBinary[1] is ", valueOfbugBinary[1])

                # 用缺陷为大于0的模块数占所有模块之比
                BaseProbability_1 = valueOfbugBinary[1] / (valueOfbugBinary[0] + valueOfbugBinary[1])

                # 计算VARL阈值及标准差
                varl_threshold = (np.log(BaseProbability_1 / (1 - BaseProbability_1)) - B[1]) / B[0]
                varl_threshold_se = ((cov.iloc[0, 0] + 2 * varl_threshold * cov.iloc[0, 1]
                                      + varl_threshold * varl_threshold * cov.iloc[1, 1]) ** 0.5) / B[0]
                varl_threshold_variance = varl_threshold_se ** 2

                # 输出每一度量的结果
                writer.writerow([file, metric, Corr_metric_bug[metric][1], B[0], pValueLogit[0], B[1], pValueLogit[1],
                                 cov11, cov12, cov22, BaseProbability_1,
                                 auc_threshold, auc_threshold_variance, auc_max_value, i_auc_max,
                                 gm_threshold, gm_threshold_variance, gm_max_value, i_gm_max,
                                 bpp_threshold, bpp_threshold_variance, bpp_max_value, i_bpp_max,
                                 mfm_threshold, mfm_threshold_variance, f1_max_value, i_f1_max,
                                 roc_threshold, roc_threshold_variance, roc_max_value, i_roc_max,
                                 varl_threshold, varl_threshold_variance])

        # break

if __name__ == '__main__':

    s_time = time.time()
    roc_threshold()
    e_time = time.time()
    execution_time = e_time - s_time

    print("The __name__ is ", __name__, ". This is end of RocMethod.py!\n",
          "The execution time of Bender.py script is ", execution_time)

max_stars_repo_name: meiyuanqing/MetaThreshold

max_issues_repo_name: meiyuanqing/MetaThreshold

max_forks_repo_name: meiyuanqing/MetaThreshold

reason: import numpy,import statsmodels


text:
from typing import Tuple, Union
import math
import numpy as np
import cv2
import pytesseract
import re

import app.lib.pyimgscan.tools as tools

#pytesseract.pytesseract.tesseract_cmd = '/usr/share/tesseract-ocr'

def rotate(image=np.ndarray, angle=float, background=Union[int, Tuple[int, int, int]]):
    old_width, old_height = image.shape[:2]
    angle_radian = math.radians(angle)
    width = abs(np.sin(angle_radian) * old_height) + abs(np.cos(angle_radian) * old_width)
    height = abs(np.sin(angle_radian) * old_width) + abs(np.cos(angle_radian) * old_height)

    image_center = tuple(np.array(image.shape[1::-1]) / 2)
    rot_mat = cv2.getRotationMatrix2D(image_center, angle, 1.0)
    rot_mat[1, 2] += (width - old_width) / 2
    rot_mat[0, 2] += (height - old_height) / 2
    return cv2.warpAffine(image, rot_mat, (int(round(height)), int(round(width))), borderValue=background)

def findAngle(image, documentType, center = None, scale = 1.0):
    image = cv2.resize(image, None, fx=3, fy=3, interpolation=cv2.INTER_CUBIC)
    documentconfig = ''

    # Set the config of Tesseract depending on the type of document
    if documentType == 'document':
        documentconfig = ''
    else:
        documentconfig = '--dpi 300 --psm 0 -c min_characters_to_try=1'

    angle = 360-int(re.search('(?<=Rotate: )\d+', pytesseract.image_to_osd(image, config=documentconfig)).group(0))

    return angle

#USED
def deskew_img(image, documentType):
    grayscale = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

    # Blur the image
    #grayscale = cv2.GaussianBlur(grayscale, (9, 9), 0)

    # Using deskew (old solution for identifying the angle)
    # angle = determine_skew(grayscale)

    # Tresseract solution for finding the angle of the image
    angle = findAngle(image, documentType)
    rotated = rotate(grayscale, angle, (0, 0, 0))
    deskewed_img = rotated.astype(np.uint8)
    padded_img = np.pad(deskewed_img, pad_width=30, mode='constant', constant_values=0)
    avg_brightness = np.mean(grayscale)

    return padded_img, avg_brightness

#USED
def pyscanimg(image, avg_brightness):
    """
    Main Proccess of the Program
    """
    img_adj, scale, _, img_edge = tools.preprocess(image, avg_brightness)
    img_hull = tools.gethull(img_edge)
    corners = tools.getcorners(img_hull)
    corners = corners.reshape(4, 2) * scale
    img_corrected = tools.perspective_transform(img_adj, corners)
    return img_corrected

#USED
def remove_shadow(img):
    rgb_planes = cv2.split(img)

    result_planes = []
    result_norm_planes = []
    for plane in rgb_planes:
        dilated_img = cv2.dilate(plane, np.ones((5,5), np.uint8))
        bg_img = cv2.medianBlur(dilated_img, 21)
        diff_img = 255 - cv2.absdiff(plane, bg_img)
        norm_img = cv2.normalize(diff_img, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8UC1)
        result_planes.append(diff_img)
        result_norm_planes.append(norm_img)

    return cv2.merge(result_planes)

#USED
def sharpen(img_grey):
    kernel = np.array([[0, -1, 0], \
                      [-1, 5,-1], \
                      [0, -1, 0]])
    return cv2.filter2D(img_grey, -1, kernel)

def align(img, documentType="", color="bw"):
    color = color.lower()
    documentType = documentType.lower()
    image_rotated, avg_brightness = deskew_img(img, documentType)
    image_corrected = pyscanimg(image_rotated, avg_brightness)
    image_filtered = remove_shadow(image_corrected)
    image_sharpened = sharpen(image_filtered)
    if color == "bw":
        (_, blackAndWhiteImage) = cv2.threshold(image_sharpened, 240, 255, cv2.THRESH_BINARY+cv2.THRESH_OTSU)
        return blackAndWhiteImage
    elif color == "grey" or color == "gray":
        return image_sharpened
    else:
        return None

max_stars_repo_name: JosiahVillanueva/kara

max_issues_repo_name: JosiahVillanueva/kara

max_forks_repo_name: JosiahVillanueva/kara


text:
'''
Author: Siyun WANG
'''
import numpy as np

class findPeaks():
    def __init__(self, signal, maxDepth, minWidth, percentile=50):
        '''
        Recursively finds all the peaks by finding their starting point and the width using
        signal: array-like object, input signal whose local maximum to be found.
        maxDepth: positive integer, maximum depth of the binary search tree. The larger the more local maximum will be found, but may also include more noise.
        minWidth: positive integer, minimum width of the found "peak". This also constrains the depth of the binary search tree.
        percentile: integer between 0 and 100, optional, default to 50.
                    The threshold greater than which values are considered to be potential peaks and will be scanned by the binary searching process
        '''
        self.signal = signal
        self.maxDepth = maxDepth
        self.minWidth = minWidth
        self.percentile = percentile
        self.rslt = []

    def scan_peaks_dicho(self, startInd, width, depth):
        val = np.percentile(self.signal[startInd:startInd+width],self.percentile)
        isPeak = (self.signal[startInd:startInd+width]>=val)

        subStartInds = []
        subWidths = []
        counter = 0
        # get contiguous sub-sequences of values greater than the threshold defined by percentile
        for j in range(len(isPeak)):
            if isPeak[j]:
                if counter == 0:
                    subStartInds.append(j)
                counter +=1
            else:
                if counter > 0:
                    subWidths.append(counter)
                counter = 0
        if len(subWidths) == len(subStartInds)-1:
            subWidths.append(counter)
        subStartInds,subWidths = np.array(subStartInds)+startInd,np.array(subWidths)
        # here expands the "binary search tree"
        for i in range(len(subStartInds)):
            if (subWidths[i] <= self.minWidth) | (depth+1 > self.maxDepth):
                self.rslt.append((subStartInds[i],subWidths[i]))
                continue
            self.scan_peaks_dicho(subStartInds[i], subWidths[i], depth+1)

max_stars_repo_name: ChatCatKatzeNeko/FindPeaks

max_issues_repo_name: ChatCatKatzeNeko/FindPeaks

max_forks_repo_name: ChatCatKatzeNeko/FindPeaks


text:
"""
Chi-squared and related functions
"""
#***************************************************************************************************
# Copyright 2015, 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights
# in this software.
# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
# in compliance with the License.  You may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory.
#***************************************************************************************************

import numpy as _np

from pygsti.tools.legacytools import deprecate as _deprecated_fn

def chi2(model, dataset, circuits=None,
         min_prob_clip_for_weighting=1e-4, prob_clip_interval=(-10000, 10000),
         op_label_aliases=None, mdc_store=None, comm=None, mem_limit=None):
    """
    Computes the total (aggregate) chi^2 for a set of circuits.

    The chi^2 test statistic obtained by summing up the
    contributions of a given set of circuits or all
    the circuits available in a dataset.  For the gradient or
    Hessian, see the :function:`chi2_jacobian` and
    :function:`chi2_hessian` functions.

    Parameters
    ----------
    model : Model
        The model used to specify the probabilities and SPAM labels

    dataset : DataSet
        The data used to specify frequencies and counts

    circuits : list of Circuits or tuples, optional
        List of circuits whose terms will be included in chi^2 sum.
        Default value (None) means "all strings in dataset".

    min_prob_clip_for_weighting : float, optional
        defines the clipping interval for the statistical weight.

    prob_clip_interval : tuple, optional
        A `(min, max)` tuple that specifies the minium (possibly negative) and maximum values
        allowed for probabilities generated by the model.  If the model gives probabilities
        outside this range they are clipped to `min` or `max`.  These values can be quite
        generous, as the optimizers are quite tolerant of badly behaved probabilities.

    op_label_aliases : dictionary, optional
        Dictionary whose keys are operation label "aliases" and whose values are tuples
        corresponding to what that operation label should be expanded into before querying
        the dataset. Defaults to the empty dictionary (no aliases defined)
        e.g. op_label_aliases['Gx^3'] = ('Gx','Gx','Gx')

    mdc_store : ModelDatasetCircuitsStore, optional
        An object that bundles cached quantities along with a given model, dataset, and circuit
        list.  If given, `model` and `dataset` and `circuits` should be set to None.

    comm : mpi4py.MPI.Comm, optional
        When not None, an MPI communicator for distributing the computation
        across multiple processors.

    mem_limit : int, optional
        A rough memory limit in bytes which restricts the amount of intermediate
        values that are computed and stored.

    Returns
    -------
    chi2 : float
        chi^2 value, equal to the sum of chi^2 terms from all specified circuits
    """
    from ..objectivefns import objectivefns as _objfns
    return _objfns._objfn(_objfns.Chi2Function, model, dataset, circuits,
                          {'min_prob_clip_for_weighting': min_prob_clip_for_weighting},
                          {'prob_clip_interval': prob_clip_interval},
                          op_label_aliases, comm, mem_limit, ('fn',), (), mdc_store).fn()  # gathers internally

def chi2_per_circuit(model, dataset, circuits=None,
                     min_prob_clip_for_weighting=1e-4, prob_clip_interval=(-10000, 10000),
                     op_label_aliases=None, mdc_store=None, comm=None, mem_limit=None):
    """
    Computes the per-circuit chi^2 contributions for a set of cirucits.

    This function returns the same value as :func:`chi2` except the
    contributions from different circuits are not summed but
    returned as an array (the contributions of all the outcomes of a
    given cirucit *are* summed together).

    Parameters
    ----------
    model : Model
        The model used to specify the probabilities and SPAM labels

    dataset : DataSet
        The data used to specify frequencies and counts

    circuits : list of Circuits or tuples, optional
        List of circuits whose terms will be included in chi^2 sum.
        Default value (None) means "all strings in dataset".

    min_prob_clip_for_weighting : float, optional
        defines the clipping interval for the statistical weight.

    prob_clip_interval : tuple, optional
        A `(min, max)` tuple that specifies the minium (possibly negative) and maximum values
        allowed for probabilities generated by the model.  If the model gives probabilities
        outside this range they are clipped to `min` or `max`.  These values can be quite
        generous, as the optimizers are quite tolerant of badly behaved probabilities.

    op_label_aliases : dictionary, optional
        Dictionary whose keys are operation label "aliases" and whose values are tuples
        corresponding to what that operation label should be expanded into before querying
        the dataset. Defaults to the empty dictionary (no aliases defined)
        e.g. op_label_aliases['Gx^3'] = ('Gx','Gx','Gx')

    mdc_store : ModelDatasetCircuitsStore, optional
        An object that bundles cached quantities along with a given model, dataset, and circuit
        list.  If given, `model` and `dataset` and `circuits` should be set to None.

    comm : mpi4py.MPI.Comm, optional
        When not None, an MPI communicator for distributing the computation
        across multiple processors.

    mem_limit : int, optional
        A rough memory limit in bytes which restricts the amount of intermediate
        values that are computed and stored.

    Returns
    -------
    chi2 : numpy.ndarray
        Array of length either `len(circuits)` or `len(dataset.keys())`.
        Values are the chi2 contributions of the corresponding circuit
        aggregated over outcomes.
    """
    from ..objectivefns import objectivefns as _objfns
    obj = _objfns._objfn(_objfns.Chi2Function, model, dataset, circuits,
                         {'min_prob_clip_for_weighting': min_prob_clip_for_weighting},
                         {'prob_clip_interval': prob_clip_interval},
                         op_label_aliases, comm, mem_limit, ('percircuit',), (), mdc_store)
    return obj.layout.allgather_local_array('c', obj.percircuit())

def chi2_jacobian(model, dataset, circuits=None,
                  min_prob_clip_for_weighting=1e-4, prob_clip_interval=(-10000, 10000),
                  op_label_aliases=None, mdc_store=None, comm=None, mem_limit=None):
    """
    Compute the gradient of the chi^2 function computed by :function:`chi2`.

    The returned value holds the derivatives of the chi^2 function with
    respect to `model`'s parameters.

    Parameters
    ----------
    model : Model
        The model used to specify the probabilities and SPAM labels

    dataset : DataSet
        The data used to specify frequencies and counts

    circuits : list of Circuits or tuples, optional
        List of circuits whose terms will be included in chi^2 sum.
        Default value (None) means "all strings in dataset".

    min_prob_clip_for_weighting : float, optional
        defines the clipping interval for the statistical weight.

    prob_clip_interval : tuple, optional
        A `(min, max)` tuple that specifies the minium (possibly negative) and maximum values
        allowed for probabilities generated by the model.  If the model gives probabilities
        outside this range they are clipped to `min` or `max`.  These values can be quite
        generous, as the optimizers are quite tolerant of badly behaved probabilities.

    op_label_aliases : dictionary, optional
        Dictionary whose keys are operation label "aliases" and whose values are tuples
        corresponding to what that operation label should be expanded into before querying
        the dataset. Defaults to the empty dictionary (no aliases defined)
        e.g. op_label_aliases['Gx^3'] = ('Gx','Gx','Gx')

    mdc_store : ModelDatasetCircuitsStore, optional
        An object that bundles cached quantities along with a given model, dataset, and circuit
        list.  If given, `model` and `dataset` and `circuits` should be set to None.

    comm : mpi4py.MPI.Comm, optional
        When not None, an MPI communicator for distributing the computation
        across multiple processors.

    mem_limit : int, optional
        A rough memory limit in bytes which restricts the amount of intermediate
        values that are computed and stored.

    Returns
    -------
    numpy array
        The gradient vector of length `model.num_params`, the number of model parameters.
    """
    from ..objectivefns import objectivefns as _objfns
    obj = _objfns._objfn(_objfns.Chi2Function, model, dataset, circuits,
                         {'min_prob_clip_for_weighting': min_prob_clip_for_weighting},
                         {'prob_clip_interval': prob_clip_interval},
                         op_label_aliases, comm, mem_limit, ('jacobian',), (), mdc_store)
    return obj.layout.allgather_local_array('ep', obj.jacobian())

def chi2_hessian(model, dataset, circuits=None,
                 min_prob_clip_for_weighting=1e-4, prob_clip_interval=(-10000, 10000),
                 op_label_aliases=None, mdc_store=None, comm=None, mem_limit=None):
    """
    Compute the Hessian matrix of the :func:`chi2` function.

    Parameters
    ----------
    model : Model
        The model used to specify the probabilities and SPAM labels

    dataset : DataSet
        The data used to specify frequencies and counts

    circuits : list of Circuits or tuples, optional
        List of circuits whose terms will be included in chi^2 sum.
        Default value (None) means "all strings in dataset".

    min_prob_clip_for_weighting : float, optional
        defines the clipping interval for the statistical weight.

    prob_clip_interval : tuple, optional
        A `(min, max)` tuple that specifies the minium (possibly negative) and maximum values
        allowed for probabilities generated by the model.  If the model gives probabilities
        outside this range they are clipped to `min` or `max`.  These values can be quite
        generous, as the optimizers are quite tolerant of badly behaved probabilities.

    op_label_aliases : dictionary, optional
        Dictionary whose keys are operation label "aliases" and whose values are tuples
        corresponding to what that operation label should be expanded into before querying
        the dataset. Defaults to the empty dictionary (no aliases defined)
        e.g. op_label_aliases['Gx^3'] = ('Gx','Gx','Gx')

    mdc_store : ModelDatasetCircuitsStore, optional
        An object that bundles cached quantities along with a given model, dataset, and circuit
        list.  If given, `model` and `dataset` and `circuits` should be set to None.

    comm : mpi4py.MPI.Comm, optional
        When not None, an MPI communicator for distributing the computation
        across multiple processors.

    mem_limit : int, optional
        A rough memory limit in bytes which restricts the amount of intermediate
        values that are computed and stored.

    Returns
    -------
    numpy array
        The Hessian matrix of shape (nModelParams, nModelParams), where
        nModelParams = `model.num_params`.
    """
    from ..objectivefns import objectivefns as _objfns
    obj = _objfns._objfn(_objfns.Chi2Function, model, dataset, circuits,
                         {'min_prob_clip_for_weighting': min_prob_clip_for_weighting},
                         {'prob_clip_interval': prob_clip_interval},
                         op_label_aliases, comm, mem_limit, ('hessian',), (), mdc_store)
    return obj.layout.allgather_local_array('epp', obj.hessian())

def chi2_approximate_hessian(model, dataset, circuits=None,
                             min_prob_clip_for_weighting=1e-4, prob_clip_interval=(-10000, 10000),
                             op_label_aliases=None, mdc_store=None, comm=None, mem_limit=None):
    """
    Compute and approximate Hessian matrix of the :func:`chi2` function.

    This approximation neglects terms proportional to the Hessian of the
    probabilities w.r.t. the model parameters (which can take a long time
    to compute).  See `logl_approximate_hessian` for details on the analogous approximation
    for the log-likelihood Hessian.

    Parameters
    ----------
    model : Model
        The model used to specify the probabilities and SPAM labels

    dataset : DataSet
        The data used to specify frequencies and counts

    circuits : list of Circuits or tuples, optional
        List of circuits whose terms will be included in chi^2 sum.
        Default value (None) means "all strings in dataset".

    min_prob_clip_for_weighting : float, optional
        defines the clipping interval for the statistical weight.

    prob_clip_interval : tuple, optional
        A `(min, max)` tuple that specifies the minium (possibly negative) and maximum values
        allowed for probabilities generated by the model.  If the model gives probabilities
        outside this range they are clipped to `min` or `max`.  These values can be quite
        generous, as the optimizers are quite tolerant of badly behaved probabilities.

    op_label_aliases : dictionary, optional
        Dictionary whose keys are operation label "aliases" and whose values are tuples
        corresponding to what that operation label should be expanded into before querying
        the dataset. Defaults to the empty dictionary (no aliases defined)
        e.g. op_label_aliases['Gx^3'] = ('Gx','Gx','Gx')

    mdc_store : ModelDatasetCircuitsStore, optional
        An object that bundles cached quantities along with a given model, dataset, and circuit
        list.  If given, `model` and `dataset` and `circuits` should be set to None.

    comm : mpi4py.MPI.Comm, optional
        When not None, an MPI communicator for distributing the computation
        across multiple processors.

    mem_limit : int, optional
        A rough memory limit in bytes which restricts the amount of intermediate
        values that are computed and stored.

    Returns
    -------
    numpy array
        The Hessian matrix of shape (nModelParams, nModelParams), where
        nModelParams = `model.num_params`.
    """
    from ..objectivefns import objectivefns as _objfns
    obj = _objfns._objfn(_objfns.Chi2Function, model, dataset, circuits,
                         {'min_prob_clip_for_weighting': min_prob_clip_for_weighting},
                         {'prob_clip_interval': prob_clip_interval},
                         op_label_aliases, comm, mem_limit, ('approximate_hessian',), (), mdc_store)
    return obj.layout.allgather_local_array('epp', obj.approximate_hessian())

def chialpha(alpha, model, dataset, circuits=None,
             pfratio_stitchpt=1e-2, pfratio_derivpt=1e-2, prob_clip_interval=(-10000, 10000),
             radius=None, op_label_aliases=None,
             mdc_store=None, comm=None, mem_limit=None):
    """
    Compute the chi-alpha objective function.

    Parameters
    ----------
    alpha : float
        The alpha parameter, which lies in the interval (0,1].

    model : Model
        The model used to specify the probabilities and SPAM labels

    dataset : DataSet
        The data used to specify frequencies and counts

    circuits : list of Circuits or tuples, optional
        List of circuits whose terms will be included in chi-alpha sum.
        Default value (None) means "all strings in dataset".

    pfratio_stitchpt : float, optional
        The x-value (x = probility/frequency ratio) below which the chi-alpha
        function is replaced with it second-order Taylor expansion.

    pfratio_derivpt : float, optional
        The x-value at which the Taylor expansion derivatives are evaluated at.

    prob_clip_interval : tuple, optional
        A `(min, max)` tuple that specifies the minium (possibly negative) and maximum values
        allowed for probabilities generated by `model`.  If the `model` gives probabilities
        outside this range they are clipped to `min` or `max`.  These values can be quite
        generous, as the optimizers are quite tolerant of badly behaved probabilities.

    radius : float, optional
        If `radius` is not None then a "harsh" method of regularizing the zero-frequency
        terms (where the local function = `N*p`) is used.  If `radius` is None, then
        `fmin` is used to handle the zero-frequency terms.

    op_label_aliases : dictionary, optional
        Dictionary whose keys are operation label "aliases" and whose values are tuples
        corresponding to what that operation label should be expanded into before querying
        the dataset. Defaults to the empty dictionary (no aliases defined)
        e.g. op_label_aliases['Gx^3'] = ('Gx','Gx','Gx')

    mdc_store : ModelDatasetCircuitsStore, optional
        An object that bundles cached quantities along with a given model, dataset, and circuit
        list.  If given, `model` and `dataset` and `circuits` should be set to None.

    comm : mpi4py.MPI.Comm, optional
        When not None, an MPI communicator for distributing the computation
        across multiple processors.

    mem_limit : int, optional
        A rough memory limit in bytes which restricts the amount of intermediate
        values that are computed and stored.

    Returns
    -------
    float
    """
    from ..objectivefns import objectivefns as _objfns
    return _objfns._objfn(_objfns.ChiAlphaFunction, model, dataset, circuits,
                          {'pfratio_stitchpt': pfratio_stitchpt,
                           'pfratio_derivpt': pfratio_derivpt,
                           'radius': radius},
                          {'prob_clip_interval': prob_clip_interval},
                          op_label_aliases, comm, mem_limit, ('fn',), (), mdc_store, alpha=alpha).fn()

def chialpha_per_circuit(alpha, model, dataset, circuits=None,
                         pfratio_stitchpt=1e-2, pfratio_derivpt=1e-2, prob_clip_interval=(-10000, 10000),
                         radius=None, op_label_aliases=None,
                         mdc_store=None, comm=None, mem_limit=None):
    """
    Compute the per-circuit chi-alpha objective function.

    Parameters
    ----------
    alpha : float
        The alpha parameter, which lies in the interval (0,1].

    model : Model
        The model used to specify the probabilities and SPAM labels

    dataset : DataSet
        The data used to specify frequencies and counts

    circuits : list of Circuits or tuples, optional
        List of circuits whose terms will be included in chi-alpha sum.
        Default value (None) means "all strings in dataset".

    pfratio_stitchpt : float, optional
        The x-value (x = probility/frequency ratio) below which the chi-alpha
        function is replaced with it second-order Taylor expansion.

    pfratio_derivpt : float, optional
        The x-value at which the Taylor expansion derivatives are evaluated at.

    prob_clip_interval : tuple, optional
        A `(min, max)` tuple that specifies the minium (possibly negative) and maximum values
        allowed for probabilities generated by `model`.  If the `model` gives probabilities
        outside this range they are clipped to `min` or `max`.  These values can be quite
        generous, as the optimizers are quite tolerant of badly behaved probabilities.

    radius : float, optional
        If `radius` is not None then a "harsh" method of regularizing the zero-frequency
        terms (where the local function = `N*p`) is used.  If `radius` is None, then
        `fmin` is used to handle the zero-frequency terms.

    op_label_aliases : dictionary, optional
        Dictionary whose keys are operation label "aliases" and whose values are tuples
        corresponding to what that operation label should be expanded into before querying
        the dataset. Defaults to the empty dictionary (no aliases defined)
        e.g. op_label_aliases['Gx^3'] = ('Gx','Gx','Gx')

    mdc_store : ModelDatasetCircuitsStore, optional
        An object that bundles cached quantities along with a given model, dataset, and circuit
        list.  If given, `model` and `dataset` and `circuits` should be set to None.

    comm : mpi4py.MPI.Comm, optional
        When not None, an MPI communicator for distributing the computation
        across multiple processors.

    mem_limit : int, optional
        A rough memory limit in bytes which restricts the amount of intermediate
        values that are computed and stored.

    Returns
    -------
    numpy.ndarray
        Array of length either `len(circuits)` or `len(dataset.keys())`.
        Values are the chi-alpha contributions of the corresponding circuit
        aggregated over outcomes.
    """
    from ..objectivefns import objectivefns as _objfns
    obj = _objfns._objfn(_objfns.ChiAlphaFunction, model, dataset, circuits,
                         {'pfratio_stitchpt': pfratio_stitchpt,
                          'pfratio_derivpt': pfratio_derivpt,
                          'radius': radius},
                         {'prob_clip_interval': prob_clip_interval},
                         op_label_aliases, comm, mem_limit, ('percircuit',), (), mdc_store, alpha=alpha)
    return obj.layout.allgather_local_array('c', obj.percircuit())

@_deprecated_fn('This function will be removed soon.  Use chi2fn(...) with `p` and `1-p`.')
def chi2fn_2outcome(n, p, f, min_prob_clip_for_weighting=1e-4):
    """
    Computes chi^2 for a 2-outcome measurement.

    The chi-squared function for a 2-outcome measurement using
    a clipped probability for the statistical weighting.

    Parameters
    ----------
    n : float or numpy array
        Number of samples.

    p : float or numpy array
        Probability of 1st outcome (typically computed).

    f : float or numpy array
        Frequency of 1st outcome (typically observed).

    min_prob_clip_for_weighting : float, optional
        Defines clipping interval (see return value).

    Returns
    -------
    float or numpy array
        n(p-f)^2 / (cp(1-cp)),
        where cp is the value of p clipped to the interval
        (min_prob_clip_for_weighting, 1-min_prob_clip_for_weighting)
    """
    cp = _np.clip(p, min_prob_clip_for_weighting, 1 - min_prob_clip_for_weighting)
    return n * (p - f)**2 / (cp * (1 - cp))

@_deprecated_fn('This function will be removed soon.')
def chi2fn_2outcome_wfreqs(n, p, f):
    """
    Computes chi^2 for a 2-outcome measurement using frequency-weighting.

    The chi-squared function for a 2-outcome measurement using
    the observed frequency in the statistical weight.

    Parameters
    ----------
    n : float or numpy array
        Number of samples.

    p : float or numpy array
        Probability of 1st outcome (typically computed).

    f : float or numpy array
        Frequency of 1st outcome (typically observed).

    Returns
    -------
    float or numpy array
        n(p-f)^2 / (f*(1-f*)),
        where f* = (f*n+1)/n+2 is the frequency value used in the
        statistical weighting (prevents divide by zero errors)
    """
    f1 = (f * n + 1) / (n + 2)
    return n * (p - f)**2 / (f1 * (1 - f1))

@_deprecated_fn('Use RawChi2Function object instead')
def chi2fn(n, p, f, min_prob_clip_for_weighting=1e-4):
    """
    Computes the chi^2 term corresponding to a single outcome.

    The chi-squared term for a single outcome of a multi-outcome
    measurement using a clipped probability for the statistical
    weighting.

    Parameters
    ----------
    n : float or numpy array
        Number of samples.

    p : float or numpy array
        Probability of 1st outcome (typically computed).

    f : float or numpy array
        Frequency of 1st outcome (typically observed).

    min_prob_clip_for_weighting : float, optional
        Defines clipping interval (see return value).

    Returns
    -------
    float or numpy array
        n(p-f)^2 / cp ,
        where cp is the value of p clipped to the interval
        (min_prob_clip_for_weighting, 1-min_prob_clip_for_weighting)
    """
    from ..objectivefns import objectivefns as _objfns
    rawfn = _objfns.RawChi2Function({'min_prob_clip_for_weighting': min_prob_clip_for_weighting})
    return rawfn.terms(p, n * f, n, f)

@_deprecated_fn('Use RawFreqWeightedChi2Function object instead')
def chi2fn_wfreqs(n, p, f, min_freq_clip_for_weighting=1e-4):
    """
    Computes the frequency-weighed chi^2 term corresponding to a single outcome.

    The chi-squared term for a single outcome of a multi-outcome
    measurement using the observed frequency in the statistical weight.

    Parameters
    ----------
    n : float or numpy array
        Number of samples.

    p : float or numpy array
        Probability of 1st outcome (typically computed).

    f : float or numpy array
        Frequency of 1st outcome (typically observed).

    min_freq_clip_for_weighting : float, optional
        The minimum frequency weighting used in the weighting,
        i.e. the largest weighting factor is `1 / fmin_freq_clip_for_weighting`.

    Returns
    -------
    float or numpy array
    """
    from ..objectivefns import objectivefns as _objfns
    rawfn = _objfns.RawFreqWeightedChi2Function({'min_freq_clip_for_weighting': min_freq_clip_for_weighting})
    return rawfn.terms(p, n * f, n, f)

max_stars_repo_name: colibri-coruscans/pyGSTi

max_issues_repo_name: colibri-coruscans/pyGSTi

max_forks_repo_name: colibri-coruscans/pyGSTi


text:
"""
008. Generated AGN halo pairs

This script counts single AGN pairs. This is calculated with the AGN catalogue that
accounts for close pairs, and also for the one that doesn't.

1. Defining input parameters
2. Parameters for creating the Major Merger catalogue
3. Find pairs for chosen pixels
    3.1. Iteration: finding pairs for every radius bin

Script written by: Soumya Shreeram
Project supervised by: Johan Comparat
Date: 22nd June 2021
"""

# scipy modules
from scipy.spatial import cKDTree
from scipy.interpolate import interp1d

# astropy modules
import astropy.units as u
import astropy.io.fits as fits

from astropy.table import Table, Column
from astropy.coordinates import SkyCoord
from astropy.cosmology import FlatLambdaCDM, z_at_value

import numpy as np

# system imports
import os
import sys
import importlib as ib
import glob

# plotting imports
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns

# Load the imported file(s) that contains all the functions used in this notebooks
sys.path.append('../imported_files/')
import Exploring_DM_Halos as edh
import Agn_incidence_from_Major_Mergers as aimm
import Comparison_simulation_with_literature_data as cswl
import plotting_cswl05 as pt
import All_sky as sky
import Scaling_relations as sr

"""
1. Defining input parameters
"""

# pixel number from the simulation file
ll, ul = int(sys.argv[1]), int(sys.argv[2])
pixel_no_cont_arr = sky.allPixelNames()
pixel_no_arr = pixel_no_cont_arr[ll:ul]

# Define cosmology used in the notebook
cosmo = FlatLambdaCDM(H0=67.77*u.km/u.s/u.Mpc, Om0=0.307115)
h = 0.6777

# look back into redshifts until...
redshift_limit = 2

data_dir = os.path.join('/data24s', 'comparat', 'simulation', 'UNIT', 'ROCKSTAR_HALOS',\
         'fixedAmp_InvPhase_001')

"""
2. Parameters for creating the Major Merger catalogue

These criteria were also applied when building the model.
So for fair compairson the validation is also applied with the same criteria.

"""

# get shell volume and projected radius bins [Mpc]
r_p, shell_volume = aimm.shellVolume()

# max mass ratio to classify as a major merger [dimensionless]
mass_max = 3

# defining the redshift bin for a merger in terms of dv = c*z [km/s]
dz_cut =  0.001

# keywords can be: 'mm and dv', 'dv' or 'all'
keyword = 'mm and dv'
major_mergers_only, delta_v_cut = cswl.decideBools(keyword = keyword)

cp_agn_cat = False

# fraction of close pair agns added to the cat_AGN_all
frac_cp_agn = float(sys.argv[3])
if frac_cp_agn == 0.15:
    frac_str = '%.2f'%frac_cp_agn
else:
    frac_str = '%.1f'%frac_cp_agn

tmm_min , tmm_max = float(sys.argv[4]), float(sys.argv[5])
xoff_min, xoff_max = float(sys.argv[6]), float(sys.argv[7])
param_names = '_%.1f_tmm_%.1f_%.2f_xoff_%.2f'%(tmm_min, tmm_max ,xoff_min, xoff_max)

"""
3. Find pairs for chosen pixels
"""

count_pairs_arr2D = np.zeros((0, len(r_p)))

# iterate over the pixels of the simulation
for pixel_no in pixel_no_arr:
    all_mm_dv_idx_all = []
    count_pairs_arr = []
    # open halo file
    hd_agn, hd_halo, _ = edh.getHeaders(pixel_no, np.array([ 'agn', 'halo']))

    # extracting the conditions for choosing halos
    _, _, conditions_halo = edh.getGalaxyData(hd_halo, '', redshift_limit)
    hd_halo_z = hd_halo[conditions_halo]
    print("Halos: %d"%(len(hd_halo_z) ))

    if cp_agn_cat:
        # load the close pairs catAGN
        catAGN_dir_name = 'CP_10_sigma_1.0_frac_' + frac_str + '_%.1f_tmm_*_*_xoff_*'%tmm_min
        catAGN_dir = glob.glob(os.path.join(data_dir, catAGN_dir_name ))[0]
        catAGN_dir_file = os.path.join(data_dir, catAGN_dir, pixel_no+'.fit')
        hd_cp_agn = Table.read(catAGN_dir_file, format='fits')

        hd_agn_z = hd_cp_agn[hd_cp_agn['redshift_R']< redshift_limit]

    else:
        hd_agn_z = hd_agn[hd_agn['redshift_R']<redshift_limit]

    """
    3.1 Iteration: finding pairs for every radius bin
    """

    for r in range(len(r_p)):
        print('\n ---- pairs within radius %.3f Mpc ---'%r_p[r])
        pos_spherical_agn = aimm.getSphericalCoord(hd_agn_z)
        pos_spherical_halo = aimm.getSphericalCoord(hd_halo_z)

        # create tree
        tree_agn = cKDTree(np.transpose(np.abs(pos_spherical_agn)), leafsize=1000.0)
        tree_halo = cKDTree(np.transpose(np.abs(pos_spherical_halo)), leafsize=1000.0)

        # list of lists of all pair indicies per DM halo
        all_idx = tree_agn.query_ball_tree(tree_halo, r=r_p[r], p=2)
        count_pairs = cswl.countSelectedPairs(all_idx, string = 'All pairs: ')

        if delta_v_cut:
            # (1) choose  pairs that satisfy the delta v criteria
            all_dv_idx, count_pairs = cswl.deltaVelSelection(hd_halo_z, all_idx, dz_cut=dz_cut)

        if major_mergers_only and delta_v_cut:
            # (2) choose major pairs (major mergers) with the delta v criteria
            all_mm_dv_idx, count_pairs = cswl.majorMergerSelection(hd_halo_z, all_dv_idx)

        if major_mergers_only and not delta_v_cut:
            # (3) choose only major pairs (major mergers)
            all_mm_idx, count_pairs = cswl.majorMergerSelection(hd_halo_z, all_idx, keyword=keyword)

        pairs_selected = cswl.tuplePairArr(np.array(all_mm_dv_idx))
        count_pairs_arr.append(len(pairs_selected))
    count_pairs_arr2D = np.append(count_pairs_arr2D, [count_pairs_arr], axis=0)

local_data_dir = '../Data/pairs_z%.1f/cat_AGN_halo_pairs/'%redshift_limit
if cp_agn_cat:
    np.save(local_data_dir+'np_%s-%s_fracAGN%.2f'%(pixel_no_arr[0], pixel_no_arr[-1], frac_cp_agn)+param_names, count_pairs_arr2D, allow_pickle=True)
else:
    np.save(local_data_dir+'cat_without_CP/np_pixels_%s-%s.npy'%(pixel_no_arr[0], pixel_no_arr[-1]), count_pairs_arr2D, allow_pickle=True)

reason: import numpy,from scipy,import astropy,from astropy


text:
#!/usr/bin/env python3
#
# This program
# 1- pulls data from canvas
# 2- unpacks the files
# 3- checks for errors
#

# echo a date in a nice format
import datetime
import sys
from canvasFunctions import *
import numpy as np

logFile = datetime.datetime.now().strftime("%Y%m%d-%H%M")
sys.stderr = open('/data/informatica/canvas/log/' + logFile + '.error', 'w')
sys.stdout = open('/data/informatica/canvas/log/' + logFile + '.log', 'w')

# 1 - Sync the data from Canvas
sync_data_from_canvas()

#  2 - Get Table Lists
canvas_table = get_canvas_table_list()
master_table = get_master_table_list()

#  3 - unpack each file
unpack_canvas_data(master_table)

clean_canvas_files(master_table)

#  4 - error checking
# check_for_issues(canvas_table, master_table)

max_stars_repo_name: byu-oit/canvas-data-download

max_issues_repo_name: byu-oit/canvas-data-download

max_forks_repo_name: byu-oit/canvas-data-download


text:
import os
from PIL import Image
import numpy as np
import pandas as pd
import cv2
from tqdm import tqdm

input_path = 'Christof/assets/train_rgb_1024/'
data = pd.read_csv('Christof/assets/train.csv')
target_size = 256
target_path = 'Christof/assets/train_rgb_1024_9crop/ur/'

ids = data['Id']

for id in tqdm(ids):
    image = cv2.imread(input_path + id + '.png', cv2.IMREAD_UNCHANGED)
    cv2.imwrite(target_path + id + '.png',image[:512, 512:,:])

input_path = 'Christof/assets/test_rgb_1024/'
data = pd.read_csv('Christof/assets/sample_submission.csv')
target_path = 'Christof/assets/test_rgb_1024_9crop/ur/'

ids = data['Id']

for id in tqdm(ids):
    image = cv2.imread(input_path + id +'.png', cv2.IMREAD_UNCHANGED)
    cv2.imwrite(target_path + id + '.png',image[:512, 512:,:])

input_path = 'Christof/assets/ext_tomomi_rgb_1024/'
data_ext1 = pd.read_csv('Christof/assets/train_ext1.csv')

ids = [fn[:-5] for fn in data_ext1['Id']]

target_path = 'Christof/assets/ext_tomomi_rgb_1024_9crop/ur/'

for id in tqdm(ids):
    image = cv2.imread(input_path + id +'.png', cv2.IMREAD_UNCHANGED)
    cv2.imwrite(target_path + id + '.png',image[:512, 512:,:])

max_stars_repo_name: guitarmind/HPA-competition-solutions

max_issues_repo_name: guitarmind/HPA-competition-solutions

max_forks_repo_name: guitarmind/HPA-competition-solutions


text:
# Copyright (c) 2017 The Khronos Group Inc.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

from __future__ import division, print_function, absolute_import

import typing

import numpy as np

from nnef_tools.core import utils
from nnef_tools.core.base_graph import *

class CaffeTensor(BaseTensor['CaffeGraph', 'CaffeOperation']):

    def __init__(self,
                 graph,  # type: CaffeGraph
                 name=None,  # type: typing.Optional[str]
                 shape=None,  # type: typing.Optional[typing.List[int]]
                 dtype=None,  # type: typing.Optional[str]
                 data=None,  # type: typing.Union[None, np.ndarray, typing.Any]
                 ):
        # type: (...)->None
        super(CaffeTensor, self).__init__(graph=graph, name=name, shape=shape, dtype=dtype, data=data)

class CaffeOperation(BaseOperation['CaffeGraph', 'CaffeTensor']):

    def __init__(self, graph, name=None, inputs=None, outputs=None, attribs=None, label=None):
        super(CaffeOperation, self).__init__(graph, name=name, inputs=inputs, outputs=outputs, attribs=attribs)
        self.label = label

    def _str_dict(self):
        _dict = super(CaffeOperation, self)._str_dict()
        _dict.update([('label', self.label)])
        return _dict

class CaffeGraph(BaseGraph['CaffeTensor', 'CaffeOperation']):

    def __init__(self, name=None):
        # type: (typing.Optional[str])->None
        super(CaffeGraph, self).__init__(name=name)

    def generate_missing_names(self):
        name_generator = utils.NameGenerator(used_names=set(t.name for t in self.tensors if t.name))

        for t in self.tensors:
            if not t.name:
                if t.is_variable:
                    t.name = name_generator.get_new_name('variable')
                elif t.is_constant:
                    t.name = name_generator.get_new_name('constant')
                else:
                    assert t.producer
                    t.name = name_generator.get_new_name(t.producer.name)

        label_generator = utils.NameGenerator(used_names=set(op.label for op in self.operations if op.label))
        for op in self.operations:
            if not op.label:
                op.label = label_generator.get_new_name(op.outputs[0].name)

__all__ = [
    'CaffeTensor',
    'CaffeOperation',
    'CaffeGraph'
]

max_stars_repo_name: rgiduthuri/NNEF-Tools

max_issues_repo_name: rgiduthuri/NNEF-Tools

max_forks_repo_name: rgiduthuri/NNEF-Tools