text:
#!/usr/bin/env python

import os
import numpy as np

from Parsers.VCF import ReferenceGenome, CollectionVCF, ref_alt_variants
from Parsers.CCF import CollectionCCF
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib as mpl

if __name__ == "__main__":
    workdir = "/media/mahajrod/d9e6e5ee-1bf7-4dba-934e-3f898d9611c8/Data/LAN2xx/combined_vcf/"

    sample_set_names_list = ["PmCDA1_3d",
                             "HAP",
                             "PmCDA1_sub1_3d",
                             "PmCDA1_6d",
                             "HAP_sub1",
                             "PmCDA1_sub1_6d",
    ]
    """
    sample_set_names_list = [#"HAP",
                             "HAP_sub1"]
    """

    clustering_dir = "clustering/"
    rainfall_dir = "rainfall"
    distance_threshold = 1000
    os.chdir(workdir)
    samples_list = sorted(os.listdir("."))
    min_cluster_size = 3
    # bad_regions_file = "/home/mahajrod/genetics/desaminases/data/LAN210_v0.10m/masked_regions/LAN210_v0.10m_masked_all.gff"

    annotations_dict = {}
    annotation_synonym_dict = {"three_prime_UTR": "3'_UTR",
                               "five_prime_UTR": "5'_UTR",
                               "snoRNA": "ncRNA",
                               "snRNA": "ncRNA"
    }
    clusters_dict = {}
    clusters_dict_3 = {}
    heatmap_dir = "heatmaps/"
    data = {}
    data_3 = {}
    y_name = "Power(Size/Median)"
    x_names_list = ["Size", "Median", "Homogeneity"]
    """
    for sample_set_name in sample_set_names_list:
        print("Handling %s" % sample_set_name)

        os.chdir(workdir)
        #os.system("mkdir -p %s" % sample_set_name)
        os.chdir(sample_set_name)
        #os.system("mkdir -p %s %s" % (clustering_dir, rainfall_dir))
        #os.system("pwd")
        clusters_dict[sample_set_name] = CollectionCCF(from_file=True, input_file="%s%s_adjusted_size_3+.ccf" %
                                                                                (clustering_dir,
                                                                                sample_set_name))

        data[sample_set_name] = clusters_dict[sample_set_name].get_data_for_stat(additional_data=("Median", "Power")) \
            if "HAP" in sample_set_name \
            else clusters_dict[sample_set_name].get_data_for_stat(additional_data=("Median", "Power", "Homogeneity"))

        data[sample_set_name] = [data[sample_set_name][:, i]
                                 for i in range(0, 4 if "HAP" in sample_set_name else 5)]

    os.chdir(workdir)
    os.system("mkdir -p %s" % heatmap_dir)
    os.chdir(heatmap_dir)

    x_names_dict = {"Length": 0, "Size": 1, "Median": 2, "Power": 3, "Homogeneity": 4}
    figure = plt.figure(1, dpi=150, figsize=(36, 18))
    subplot_dict = {}
    for i in range(0, len(x_names_list)):
        subplot_dict[x_names_list[i]] = []
        for j in range(0, len(sample_set_names_list)):
            if i == 2 and "HAP" in sample_set_names_list[j]:
                continue

            if i == 0 and j == 0:
                subplot_dict[x_names_list[i]].append(plt.subplot(3, 6, i * 6 + j + 1))
            else:
                if j == 0:
                    subplot_dict[x_names_list[i]].append(plt.subplot(3, 6, i * 6 + j + 1,
                                                                 sharey=subplot_dict[x_names_list[0]][0]))
                else:
                    subplot_dict[x_names_list[i]].append(plt.subplot(3, 6, i * 6 + j + 1,
                                                                 sharey=subplot_dict[x_names_list[0]][0],
                                                                 sharex=subplot_dict[x_names_list[i]][0]))
            x_data_index = x_names_dict[x_names_list[i]]
            x_data = data[sample_set_names_list[j]][x_data_index]
            y_data = data[sample_set_names_list[j]][3]

            max_y = max(y_data)
            max_x = max(x_data)
            n_y = int(max_y / 0.05) + 1

            if x_names_list[i] == "Size":
                n_x_bins = np.linspace(3, 30, 27 + 1)
                #print(n_x_bins)
                plt.xlim(xmax=30)
            elif x_names_list[i] == "Median":
                plt.xlim(xmax=1500)
                n_x_bins = np.linspace(0, 800, 16 + 1)
            elif x_names_list[i] == "Homogeneity":
                plt.xlim(xmax=1.0)
                n_x_bins = np.linspace(0.5, 1.0, 10 + 1)
            plt.xlabel(x_names_list[i])
            plt.ylabel("Power(Size/Median dist.)")
            plt.title(sample_set_names_list[j] + " (%i clusters)" % len(y_data))

            n_y_bins = np.linspace(0, 0.05 * n_y, n_y + 1)

            #print(i, j, sample_set_names_list[j])
            #plt.
            plt.hist2d(x_data, y_data, (n_x_bins, n_y_bins), cmin=1)
            plt.ylim(ymax=3)
            plt.colorbar()

    #figure.subplots_adjust(right=0.8)
    #cbar_ax = figure.add_axes([0.85, 0.15, 0.05, 0.7])
    #plt.colorbar(cax=cbar_ax)

    plt.savefig("heatmap_3+.svg")
    plt.close()
    """
    """
    for sample_set_name in sample_set_names_list:
        print("Handling %s" % sample_set_name)

        os.chdir(workdir)
        #os.system("mkdir -p %s" % sample_set_name)
        os.chdir(sample_set_name)
        #os.system("mkdir -p %s %s" % (clustering_dir, rainfall_dir))
        #os.system("pwd")
        clusters_dict[sample_set_name] = CollectionCCF(from_file=True, input_file="%s%s_adjusted_size_3+_power_0.05+.ccf" %
                                                                                (clustering_dir,
                                                                                sample_set_name))

        data[sample_set_name] = clusters_dict[sample_set_name].get_data_for_stat(additional_data=("Median", "Power")) \
            if "HAP" in sample_set_name \
            else clusters_dict[sample_set_name].get_data_for_stat(additional_data=("Median", "Power", "Homogeneity"))

        data[sample_set_name] = [data[sample_set_name][:, i]
                                 for i in range(0, 4 if "HAP" in sample_set_name else 5)]

    os.chdir(workdir)
    os.system("mkdir -p %s" % heatmap_dir)
    os.chdir(heatmap_dir)

    x_names_dict = {"Length": 0, "Size": 1, "Median": 2, "Power": 3, "Homogeneity": 4}
    figure = plt.figure(1, dpi=150, figsize=(36, 18))
    subplot_dict = {}
    for i in range(0, len(x_names_list)):
        subplot_dict[x_names_list[i]] = []
        for j in range(0, len(sample_set_names_list)):
            if i == 2 and "HAP" in sample_set_names_list[j]:
                continue

            if i == 0 and j == 0:
                subplot_dict[x_names_list[i]].append(plt.subplot(3, 6, i * 6 + j + 1))
            else:
                if j == 0:
                    subplot_dict[x_names_list[i]].append(plt.subplot(3, 6, i * 6 + j + 1,
                                                                 sharey=subplot_dict[x_names_list[0]][0]))
                else:
                    subplot_dict[x_names_list[i]].append(plt.subplot(3, 6, i * 6 + j + 1,
                                                                 sharey=subplot_dict[x_names_list[0]][0],
                                                                 sharex=subplot_dict[x_names_list[i]][0]))
            x_data_index = x_names_dict[x_names_list[i]]
            x_data = data[sample_set_names_list[j]][x_data_index]
            y_data = data[sample_set_names_list[j]][3]

            max_y = max(y_data)
            max_x = max(x_data)

            if x_names_list[i] == "Size":
                n_x_bins = np.linspace(3, 30, 27 + 1)
                #print(n_x_bins)
                plt.xlim(xmax=30)
            elif x_names_list[i] == "Median":
                plt.xlim(xmax=300)
                n_x_bins = np.linspace(0, 200, 10 + 1)
            elif x_names_list[i] == "Homogeneity":
                plt.xlim(xmax=1.0)
                n_x_bins = np.linspace(0.5, 1.0, 10 + 1)
            plt.xlabel(x_names_list[i])
            plt.ylabel("Power(Size/Median dist.)")
            plt.title(sample_set_names_list[j] + " (%i clusters)" % len(y_data))
            n_y = int(max_y / 0.05) + 1
            n_y_bins = np.linspace(0, 0.05 * n_y, n_y + 1)

            #print(i, j, sample_set_names_list[j])
            #plt.
            plt.hist2d(x_data, y_data, (n_x_bins, n_y_bins), cmin=1)
            plt.ylim(ymax=3)
            plt.colorbar()

    #figure.subplots_adjust(right=0.8)
    #cbar_ax = figure.add_axes([0.85, 0.15, 0.05, 0.7])
    #plt.colorbar(cax=cbar_ax)

    plt.savefig("heatmap_3+_power_0.05+.svg")
    plt.savefig("heatmap_3+_power_0.05+.eps")
    plt.close()
    """

    for sample_set_name in sample_set_names_list:
        print("Handling %s" % sample_set_name)

        os.chdir(workdir)
        #os.system("mkdir -p %s" % sample_set_name)
        os.chdir(sample_set_name)
        #os.system("mkdir -p %s %s" % (clustering_dir, rainfall_dir))
        #os.system("pwd")
        clusters_dict[sample_set_name] = CollectionCCF(from_file=True, input_file="%s%s_adjusted_size_3+_power_0.1+.ccf" %
                                                                                (clustering_dir,
                                                                                sample_set_name))

        data[sample_set_name] = clusters_dict[sample_set_name].get_data_for_stat(additional_data=("Median", "Power")) \
            if "HAP" in sample_set_name \
            else clusters_dict[sample_set_name].get_data_for_stat(additional_data=("Median", "Power", "Homogeneity"))

        data[sample_set_name] = [data[sample_set_name][:, i]
                                 for i in range(0, 4 if "HAP" in sample_set_name else 5)]

    os.chdir(workdir)
    os.system("mkdir -p %s" % heatmap_dir)
    os.chdir(heatmap_dir)

    x_names_dict = {"Length": 0, "Size": 1, "Median": 2, "Power": 3, "Homogeneity": 4}
    figure = plt.figure(1, dpi=150, figsize=(36, 18))
    subplot_dict = {}
    for i in range(0, len(x_names_list)):
        subplot_dict[x_names_list[i]] = []
        for j in range(0, len(sample_set_names_list)):
            if i == 2 and "HAP" in sample_set_names_list[j]:
                continue

            if i == 0 and j == 0:
                subplot_dict[x_names_list[i]].append(plt.subplot(3, 6, i * 6 + j + 1))
            else:
                if j == 0:
                    subplot_dict[x_names_list[i]].append(plt.subplot(3, 6, i * 6 + j + 1,
                                                                 sharey=subplot_dict[x_names_list[0]][0]))
                else:
                    subplot_dict[x_names_list[i]].append(plt.subplot(3, 6, i * 6 + j + 1,
                                                                 sharey=subplot_dict[x_names_list[0]][0],
                                                                 sharex=subplot_dict[x_names_list[i]][0]))
            x_data_index = x_names_dict[x_names_list[i]]
            x_data = data[sample_set_names_list[j]][x_data_index]
            y_data = data[sample_set_names_list[j]][3]

            max_y = max(y_data)
            max_x = max(x_data)

            if x_names_list[i] == "Size":
                n_x_bins = np.linspace(3, 30, 27 + 1)
                #print(n_x_bins)
                plt.xlim(xmax=30)
            elif x_names_list[i] == "Median":
                plt.xlim(xmax=300)
                n_x_bins = np.linspace(0, 200, 10 + 1)
            elif x_names_list[i] == "Homogeneity":
                plt.xlim(xmax=1.0)
                n_x_bins = np.linspace(0.5, 1.0, 10 + 1)
            plt.xlabel(x_names_list[i])
            plt.ylabel("Power(Size/Median dist.)")
            plt.title(sample_set_names_list[j] + " (%i clusters)" % len(y_data))
            n_y = int(max_y / 0.05) + 1
            n_y_bins = np.linspace(0, 0.05 * n_y, n_y + 1)

            #print(i, j, sample_set_names_list[j])
            #plt.
            plt.hist2d(x_data, y_data, (n_x_bins, n_y_bins), cmin=1)
            plt.ylim(ymax=3)
            plt.colorbar()

    #figure.subplots_adjust(right=0.8)
    #cbar_ax = figure.add_axes([0.85, 0.15, 0.05, 0.7])
    #plt.colorbar(cax=cbar_ax)

    plt.savefig("heatmap_3+_power_0.1+.svg")
    plt.savefig("heatmap_3+_power_0.1+.eps")
    plt.close()

max_stars_repo_name: mahajrod/MAVR

max_issues_repo_name: mahajrod/MAVR

max_forks_repo_name: mahajrod/MAVR


text:
import cv2
import numpy as np
from matplotlib import pyplot as plt
from datetime import datetime
import operator

sigma = 0.33
spaceBetweenLinesDistance = 275
deltaDistance = 25

eps = 0.01
theta1 = 0
theta2 = 1.57

photoName = 'outputFiles/' + datetime.now().strftime('%Y-%m-%d-%H-%M-%S')

def drawLines(lines, img):
    for line in lines:
        (rho,theta) = line
        if(abs(theta-theta1) < eps):
            drawLine(rho, theta, img)
        elif(abs(theta-theta2) < eps):
            drawLine(rho, theta, img)

def drawLine(rho, theta, img):
            a = np.cos(theta)
            b = np.sin(theta)
            x0 = a*rho
            y0 = b*rho
            x1 = int(x0 + 3000*(-b))
            y1 = int(y0 + 3000*(a))
            x2 = int(x0 - 3000*(-b))
            y2 = int(y0 - 3000*(a))
            cv2.line(img,(x1,y1),(x2,y2),(0,0,255),20)

def inRange(distance):
    if((distance > spaceBetweenLinesDistance - deltaDistance) and (distance < spaceBetweenLinesDistance + deltaDistance)):
        return True
    return False

# chess lines will be either vertical or horizontal - 18 lines in total
# the difference between 2 consecutive lines will be spaceBetweenLinesDistance +- deltaDistance (250,300)

def findChessLines(lines):
    if(len(lines) < 9):
        return []
    for index, item in enumerate(lines):
        chessLines = []
        chessLines.append(item)
        currentIndex = index + 1
        while(currentIndex < len(verticalLines)):
            if(inRange(lines[currentIndex][0] - chessLines[len(chessLines) - 1][0])):
                chessLines.append(lines[currentIndex])
                currentIndex += 1
                if(len(chessLines) == 9):
                    return chessLines
            else:
                break
    return []

img = cv2.imread('inputFiles\\raw.jpg')
imgCopy = img
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

v = np.median(img)

edges = cv2.Canny(img,150,200)

lines = cv2.HoughLines(edges,10,np.pi/90,800)
linesCount = 0

verticalLines = []
horizontalLines = []
for line in lines:
    for rho,theta in line:
        if(abs(theta-theta1) < eps):
            verticalLines.append((rho,theta))
            linesCount += 1
            drawLine(rho, theta, imgCopy)
        elif(abs(theta-theta2) < eps):
            horizontalLines.append((rho,theta))
            linesCount += 1
            drawLine(rho, theta, imgCopy)

print(str(len(verticalLines)))
print(str(len(horizontalLines)))

if(len(verticalLines) == 0 or len(horizontalLines) == 0):
    exit()

verticalLines.sort(key=operator.itemgetter(0))
horizontalLines.sort(key=operator.itemgetter(0))

chessBoardVerticalLines = findChessLines(verticalLines)
chessBoardHorizontalLines = findChessLines(horizontalLines)

print('Chessboard lines count')
print(str(len(chessBoardVerticalLines)))
print(str(len(chessBoardHorizontalLines)))

print('Vertical lines: ')
print(*chessBoardVerticalLines)
print('Horizontal lines: ')
print(*chessBoardHorizontalLines)

drawLines(chessBoardVerticalLines, img)
drawLines(chessBoardHorizontalLines, img)

print('Lines count: ' + str(linesCount))
cv2.imwrite(photoName + '-edges.jpg',edges)
cv2.imwrite(photoName + '-rawLines.jpg',img)
cv2.imwrite(photoName + '-lines.jpg',img)

# plt.subplot(121),plt.imshow(img,cmap = 'gray')
# plt.title('Original Image'), plt.xticks([]), plt.yticks([])
# plt.subplot(122),plt.imshow(edges,cmap = 'gray')
# plt.title('Edge Image'), plt.xticks([]), plt.yticks([])
# plt.show()

ext: py

max_stars_repo_name: DragomirIulian/Gambinator

max_issues_repo_name: DragomirIulian/Gambinator

max_forks_repo_name: DragomirIulian/Gambinator

reason: import numpy


text:
# This is open-source software licensed under a BSD license.
# Please see the file LICENSE.txt for details.
#
from __future__ import print_function, absolute_import, unicode_literals, division

import numpy as np
from astropy.coordinates.matrix_utilities import rotation_matrix

def calculate_sky_offset(xoff, yoff, sky_pa):
    """
    Convert pixel offset to sky offset in arcseconds

    Arguments
    ---------
    xoff, yoff : float
        desired pixel offsets
    sky_pa : float
        current instrument PA
    """
    px_scale = 0.081  # arcseconds per pixel
    flipEW = True  # is E to the right in the image?
    EofN = True  # does increasing rotator PA move us to E?
    paOff = 69.3  # pa when rotator = 0

    if EofN:
        theta = sky_pa - paOff
    else:
        theta = -sky_pa - paOff
    # only take 3d of 3d rotation matrix
    rmat = rotation_matrix(theta)[:2, :2]

    # +ve shifts should move stars right and up
    # offset for skyPA = 0 are:
    ra_shift_arcsecs = -xoff*px_scale if flipEW else xoff*px_scale
    dec_shift_arcsecs = -yoff*px_scale
    pix_shift = np.array([ra_shift_arcsecs, dec_shift_arcsecs])
    return rmat.dot(pix_shift)

max_stars_repo_name: StuartLittlefair/scam-drivers

max_issues_repo_name: StuartLittlefair/scam-drivers

max_forks_repo_name: StuartLittlefair/scam-drivers

reason: import numpy,from astropy


text:
from scipy.stats import f
import numpy
from scipy.special import erf
from scipy.optimize import fsolve
from PyAstronomy.pyaC import pyaErrors as PE

def chi2test():
    pass

def bayesFactors(model1, model2):
    raise(PE.PyAUnclassifiedError(
        "The function 'bayesFactors' has moved to funcFit.utils."))

def ftest(chi1, chi2, dof1, dof2, compute_ratio_of_variance=False):
    """
    Performs an F-test.

    Parameters
    ----------
    chi1 : float
        The chi-squared statistic or variance for model 1.
    chi2 : float
        The chi-squared statistic or variance for model 2.
        Should be better (lower) than `chi1`.
    dof1 : integer
        The number of degrees of freedom of model 1
    dof2 : integer
        The number of degrees of freedom of model 2.
        Should be lower than `dof1`.
    compute_ratio_of_variance : boolean, optional
        Determines how the F-statistics is computed.
        If False (default), the "ANOVA F-test in regression analysis
        of nested nonlinear models" is applied.
        If True, "Fisher's classical test for the equality
        of two variances" is used.

    Returns
    -------
    Test results : dictionary
        - "F statistic" - Value of F statistic.
        - "p-value" - Probability to obtain an F-value as large as
          or larger than observed assuming H_0 to be true
        - "Gaussian sigma level" - Same probability converted to
          Gaussian standard deviations ("sigmas")

    Notes
    -----

    The "method of least squares" is widely used in parameter estimation.
    Much of its appeal lies in its simplicity: Minimizing the sum of squared
    differences between the model and the data.

    The backbone of least squares is the classical multiple regression
    analysis using linear models to relate several independent variables to
    a dependent variable. However, many physical models are nonlinear;
    this is the realm of nonlinear least squares regression methods.

    The aim of model fitting is to explain as much of the variation in the
    dependent variable as possible from information contained in the
    independent variables. The contributions of the independent variables
    to the model are measured by partitions of the total sum of squares of
    the difference of data and model ("analysis of variance", ANOVA).

    **ANOVA F-test in regression analysis for nested nonlinear models**

    The sum of squares for any hypothesis can be determined from the
    difference between the residual sums of squares (:math:`RSS`) of two models: The
    so-called "full" model, in the context of which the null hypothesis is to be
    tested, and the "reduced" model, which is derived from
    the full model by imposing additional constraints specified by the null hypothesis;
    setting one parameter to zero is a common example.
    The reduced model is a special case of the full model and, hence, its
    residual sum of squares must not be lower than the residual
    sum of squares for the full model.

    Accounting for :math:`k` independent constraints on the full model and :math:`N` data
    points, we compute

    .. math:: Q = RSS(reduced) - RSS(full)

    which has :math:`k` degrees of freedom. Here, :math:`k` corresponds to the
    difference in the number of parameters---or equivalently, the difference
    in the degrees of freedoms---between the full and reduced model.
    The :math:`F`-test of the null hypothesis can
    then be computed as

    .. math:: \\hat{F} = (Q/k) / s^2

    where :math:`s^2` is an unbiased estimate of :math:`\\sigma^2`, e.g., derived from the
    full model. In the case of error-weighted data, the :math:`F` statistic reads

    .. math:: \\hat{F} = \\frac{(\chi^2_{reduced}-\chi^2_{full})/(\\nu_{reduced}-\\nu_{full})}{\chi^2_{full}/\\nu_{full}}

    where :math:`\\nu` denotes the number of degrees of freedom, which may be calculated as
    :math:`N-p-1` given a model with :math:`p` parameters and an additional constant.
    The expectation value of :math:`F` is :math:`1`, thus, if :math:`\\hat{F}=1`
    there is no significant difference
    between the RSS of both models. If, however, :math:`F` deviates from :math:`1`, we can compute
    the probability for :math:`F` to equal or exceed the obtained value :math:`\\hat{F}` by

    .. math:: Prob(F\\geq \\hat{F}) = CDF(F( \\nu_{reduced}-\\nu_{full}, \\nu_{full} ))

    For details on the principles of the formalism,
    see Sect. 4.5.4 in Rawlings' "Applied Regression Analysis".

    In case of nonlinear models, the distribution of the least-square
    parameter estimates only approximately follows the normal distribution.
    Then, the :math:`F` statistic (in this case called "Wald statistic") also
    holds approximately.
    The validity of the method depends on how well the model is represented
    by a linear approximation in the parameters ("parameter effects curvature").
    For details, we refer to Sect. 15.3 in Rawlings' "Applied Regression Analysis"
    and references therein.

    **Fisher's classical test for the equality of two variances**

    If the two vying models cannot be treated as *nested*, in the sense that
    the full model encompasses the reduced model, we can directly compare
    the model variances with an F-test.

    Let :math:`s_1` and :math:`s_2` denote unbiased estimates of the variances
    of two independent, normally distributed populations of size :math:`N_1` and :math:`N_2`,
    respectively. Under the null hypothesis
    both variances are equal:

    .. math:: H_0: s_1^2=s_2^2, \hspace{20pt} H_1: s_1^2>s_2^2

    Then the quantity

    .. math:: \\hat{F} = \\frac{s_1^2}{s_2^2}

    is :math:`F`-distributed with :math:`N_1-1` and :math:`N_2-1` degrees of freedom.
    In the case of weighted errors we have

    .. math:: \\hat{F} = \\frac{\\chi_1^2/\\nu_1}{\\chi_2^2/\\nu_2}

    and :math:`\\hat{F}` is distributed according to an
    :math:`F`-distribution with :math:`\\nu_1` and :math:`\\nu_2`
    degrees of freedom.

    **References**

    John O. Rawlings, Sastry G. Pantula, David A. Dickey
    *Applied Regression Analysis: A Research Tool*, Second Edition
    Springer, New York, NY, 1998

    """
    # Some consistency checks
    if not dof2 < dof1:
        raise(PE.PyAValError("The number of degrees of freedom in the new fit, dof2, should be less" +
                             " than the number of degrees of freedom in the original fit, dof1."))
    if not chi2 < chi1:
        raise(PE.PyAValError(
            "The new fit (2) should have lower chi-squared than the original fit (1)."))

    # Compute the test statistic and the probability
    if compute_ratio_of_variance:
        F = (chi1 / dof1) / (chi2 / dof2)
        prob = f.cdf(F, dof1, dof2)
    else:
        F = ((chi1 - chi2) / (dof1 - dof2)) / (chi2 / dof2)
        prob = f.cdf(F, dof1 - dof2, dof2)

    # Convert probability into Gaussian 'sigmas'
    def err(n, prob):
        return erf(n / numpy.sqrt(2)) - prob
    sigmalvl = float(fsolve(err, 1.0, args=(prob,)))

    return {"F statistic": F, "p-value": 1. - prob, "Gaussian sigma level": sigmalvl}

max_stars_repo_name: mirofedurco/PyAstronomy

max_issues_repo_name: mirofedurco/PyAstronomy

max_forks_repo_name: mirofedurco/PyAstronomy

reason: import numpy,from scipy


text:
# This is to poly fit the lanelines by sliding window method
import numpy as np
import cv2
import csv
from CLine import Line

class CLocateLines:
    def __init__(self, nWindows = 9, nMargin = 100, nMinPixel = 50, shape = [720, 1280],
    xm_per_pix = 3.7/700, ym_per_pix = 30/720, data_file = 'data.txt'):
        self.__nWindows = nWindows      # Choose the number of sliding windows
        self.__nMargin = nMargin        # Set the width of the windows +/- margin
        self.__nMinPixel = nMinPixel    # Set minimum number of pixels found to recenter window
        self.__bImg = None
        self.__left_lane = Line(xm_per_pix, ym_per_pix)
        self.__right_lane = Line(xm_per_pix, ym_per_pix)
        # file = open(data_file, mode = 'w')
        # self.__csvwriter = csv.writer(file, delimiter='\t', quotechar='|', quoting=csv.QUOTE_MINIMAL)
        self.process_method = 'None'

    def update(self, binary_img):
        self.__set_binary_img(binary_img)

        if self.__left_lane.detected is False or self.__right_lane.detected is False:
            lanes_xy = self.__sliding_window()
            self.process_method = 'FirstSliding'
            self.__left_lane.detected = True
            self.__right_lane.detected = True
        else:
            lanes_xy = self.__poly_fit_prev(self.__left_lane.current_fit, self.__right_lane.current_fit)
            self.process_method = 'PrevFit'
            if lanes_xy is None or lanes_xy.__len__() == 0:
                lanes_xy = self.__sliding_window()
                self.process_method = 'ExceptionSliding'

        l_current_fit, l_radius, l_distance = self.__left_lane.try_fitting(lanes_xy[0], lanes_xy[1])
        r_current_fit, r_radius, r_distance = self.__right_lane.try_fitting(lanes_xy[2], lanes_xy[3])

        lane_width = abs(l_distance - r_distance)

        if l_radius == 0.0 or r_radius == 0.0:
            radius_ratio = 1.0
        elif l_radius < r_radius:
            radius_ratio = l_radius / r_radius
        else:
            radius_ratio = r_radius / l_radius

        if radius_ratio < 0.1 or (lane_width < 3.0 or lane_width > 4.5):
            # result is too worse, do sliding window again
            lanes_xy = self.__sliding_window()
            l_current_fit, l_radius, l_distance = self.__left_lane.try_fitting(lanes_xy[0], lanes_xy[1])
            r_current_fit, r_radius, r_distance = self.__right_lane.try_fitting(lanes_xy[2], lanes_xy[3])
            self.process_method = 'WorseResultSliding'
        else:
            # TODO use EMWA to find out a better fit
            if self.__left_lane.current_fit is not None:
                l_current_fit = self.__left_lane.current_fit * (1.0 - radius_ratio) + l_current_fit * radius_ratio
            if self.__right_lane.current_fit is not None:
                r_current_fit = self.__right_lane.current_fit * (1.0 - radius_ratio) + r_current_fit * radius_ratio

        # write debug data
        # self.__csvwriter.writerow([l_radius, r_radius, radius_ratio, lane_width, self.process_method])

        self.__left_lane.update(lanes_xy[0], lanes_xy[1], l_current_fit, l_radius, l_distance)
        self.__right_lane.update(lanes_xy[2], lanes_xy[3], r_current_fit, r_radius, r_distance)

    def visualize(self):
        ## Visualization ##
        # Create an image to draw on and an image to show the selection window
        out_img = np.dstack((self.__bImg, self.__bImg, self.__bImg))*255
        window_img = np.zeros_like(out_img)

        # Generate x and y values for plotting
        ploty = np.linspace(0, self.__bImg.shape[0] - 1, self.__bImg.shape[0])
        ### TO-DO: Calc both polynomials using ploty, left_fit and right_fit ###
        left_fitx = self.__left_lane.current_fit[0]*ploty**2 + self.__left_lane.current_fit[1]*ploty + self.__left_lane.current_fit[2]
        right_fitx = self.__right_lane.current_fit[0]*ploty**2 + self.__right_lane.current_fit[1]*ploty + self.__right_lane.current_fit[2]

        # Generate a polygon to illustrate the lane line area
        # And recast the x and y points into usable format for cv2.fillPoly()
        left_line_pts = np.array([np.transpose(np.vstack([left_fitx, ploty]))])
        right_line_pts = np.array([np.flipud(np.transpose(np.vstack([right_fitx, ploty])))])
        lane_line_poly = np.hstack((left_line_pts, right_line_pts))

        # fillPoly: shape[1, 1440, 2] -> xy coordinates with 720 * 2 size

        # Draw the lane onto the warped blank image
        cv2.fillPoly(window_img, np.int_(lane_line_poly), (0,255, 0))

        return window_img

    def get_mean_curvature(self):
        left_radius = self.__left_lane.radius_of_curvature
        right_radius = self.__right_lane.radius_of_curvature
        return (left_radius + right_radius) // 2

    def get_radius(self):
        return [self.__left_lane.radius_of_curvature, self.__right_lane.radius_of_curvature]

    def get_line_pos(self):
        lane_width = abs(self.__left_lane.line_base_pos - self.__right_lane.line_base_pos)
        return [self.__left_lane.line_base_pos, self.__right_lane.line_base_pos, lane_width]

    def get_deviation(self):
        # get deviated dist from center
        return -1.0 * (self.__left_lane.line_base_pos + self.__right_lane.line_base_pos) / 2.0

    def __set_binary_img(self, binary_img):
        self.__bImg = binary_img

    def __sliding_window(self):
        """
        fit polynomail with sliding windows
        """
        # Take a histogram of the bottom half of the image
        histogram = np.sum(self.__bImg[self.__bImg.shape[0] // 2:,:], axis = 0)
        # Find the peak of the left and right halves of the histogram
        # These will be the starting point for the left and right lines
        midpoint = np.int(histogram.shape[0]//2)
        leftx_base = np.argmax(histogram[:midpoint])
        rightx_base = np.argmax(histogram[midpoint:]) + midpoint

        # Set height of windows - based on nwindows above and image shape
        window_height = np.int(self.__bImg.shape[0] // self.__nWindows)
        # Identify the x and y positions of all nonzero pixels in the image
        nonzero = self.__bImg.nonzero()
        nonzeroy = np.array(nonzero[0])
        nonzerox = np.array(nonzero[1])
        # Current positions to be updated later for each window in nwindows
        leftx_current = leftx_base
        rightx_current = rightx_base

        # Create empty lists to receive left and right lane pixel indices
        left_lane_inds = []
        right_lane_inds = []

        # Step through the windows one by one
        for window in range(self.__nWindows):
            # Identify window boundaries in x and y (and right and left)
            win_y_low = self.__bImg.shape[0] - (window+1)*window_height
            win_y_high = self.__bImg.shape[0] - window*window_height
            ### TO-DO: Find the four below boundaries of the window ###
            win_xleft_low = leftx_current - self.__nMargin  # Update this
            win_xleft_high = leftx_current + self.__nMargin  # Update this
            win_xright_low = rightx_current - self.__nMargin  # Update this
            win_xright_high = rightx_current + self.__nMargin

            ### TO-DO: Identify the nonzero pixels in x and y within the window ###
            good_left_inds = ((nonzeroy >= win_y_low) & (nonzeroy < win_y_high) & (nonzerox >= win_xleft_low) & (nonzerox < win_xleft_high)).nonzero()[0]
            good_right_inds = ((nonzeroy >= win_y_low) & (nonzeroy < win_y_high) & (nonzerox >= win_xright_low) &  (nonzerox < win_xright_high)).nonzero()[0]

            # Append these indices to the lists
            left_lane_inds.append(good_left_inds)
            right_lane_inds.append(good_right_inds)

            if len(good_left_inds) > self.__nMinPixel:
                leftx_current = np.int(np.mean(nonzerox[good_left_inds]))
            if len(good_right_inds) > self.__nMinPixel:
                rightx_current = np.int(np.mean(nonzerox[good_right_inds]))

        # Concatenate the arrays of indices (previously was a list of lists of pixels)
        try:
            left_lane_inds = np.concatenate(left_lane_inds)
            right_lane_inds = np.concatenate(right_lane_inds)
        except ValueError:
            # Avoids an error if the above is not implemented fully
            pass

        # Extract left and right line pixel positions
        leftx = nonzerox[left_lane_inds]
        lefty = nonzeroy[left_lane_inds]
        rightx = nonzerox[right_lane_inds]
        righty = nonzeroy[right_lane_inds]

        return [leftx, lefty, rightx, righty]

    def __poly_fit_prev(self, left_fit, right_fit):
        """
        fit polynomail with previous coefficients
        """
        # Grab activated pixels
        nonzero = self.__bImg.nonzero()
        nonzeroy = np.array(nonzero[0])
        nonzerox = np.array(nonzero[1])

        ### TO-DO: Set the area of search based on activated x-values ###
        ### within the +/- margin of our polynomial function ###
        ### Hint: consider the window areas for the similarly named variables ###
        ### in the previous quiz, but change the windows to our new search area ###
        leftXMin = left_fit[0]*nonzeroy**2 + left_fit[1]*nonzeroy + left_fit[2] - self.__nMargin
        leftXMax = left_fit[0]*nonzeroy**2 + left_fit[1]*nonzeroy + left_fit[2] + self.__nMargin
        rightXMin = right_fit[0]*nonzeroy**2 + right_fit[1]*nonzeroy + right_fit[2] - self.__nMargin
        rightXMax = right_fit[0]*nonzeroy**2 + right_fit[1]*nonzeroy + right_fit[2] + self.__nMargin

        left_lane_inds = ((nonzerox >= leftXMin) & (nonzerox < leftXMax)).nonzero()[0]
        right_lane_inds = ((nonzerox >= rightXMin) & (nonzerox < rightXMax)).nonzero()[0]

        # Again, extract left and right line pixel positions
        leftx = nonzerox[left_lane_inds]
        lefty = nonzeroy[left_lane_inds]
        rightx = nonzerox[right_lane_inds]
        righty = nonzeroy[right_lane_inds]

        if leftx.size == 0 or lefty.size == 0 or rightx.size == 0 or righty.size == 0:
            return None

        # Fit new polynomials
        return [leftx, lefty, rightx, righty]

max_stars_repo_name: chifai/CarND-Advanced-Lane-Lines

max_issues_repo_name: chifai/CarND-Advanced-Lane-Lines

max_forks_repo_name: chifai/CarND-Advanced-Lane-Lines


text:
import math
import tensorflow
import keras
import numpy as np
import os
import sys
from keras.models import load_model
from keras import backend, optimizers
from keras.utils import to_categorical

###Note some code from professor template and from class####
def generate_adversary(BlackboxModel, TrainData, TargetPreds):
    TrainData = TrainData / 4.5
    cost = backend.binary_crossentropy(TargetPreds, BlackboxModel.output)
    grad = backend.gradients(cost,BlackboxModel.input)
    calculategrads = backend.function(BlackboxModel.input, grad)
    gradients = calculategrads(TrainData)[0]
    adversaries = TrainData + epsilon * np.sign(gradients)
    return(adversaries)
learning_rate = 0.001
epsilon = 0.0625
split = 200
opt =  optimizers.Adam(learning_rate = learning_rate)

TestData = sys.argv[1]
TestLabels = sys.argv[2]
targetmodelfile = sys.argv[3]

BlackBox = sys.argv[4]
target_model = load_model(targetmodelfile)

TestData = np.load(TestData)
TestData = TestData.astype(np.float32) / (255.0)
TestLabels = np.load(TestLabels)

TestData = TestData[(TestLabels == 0) | (TestLabels == 1)]
TestLabels = TestLabels[(TestLabels == 0) | (TestLabels == 1)]
TestData = TestData - np.mean(TestData , axis = 0)

testdata = TestData
testlabels = to_categorical(TestLabels, num_classes = 2)

black_box_model = load_model(BlackBox)

target_model.compile(optimizer = opt, loss = 'binary_crossentropy',metrics = ['accuracy'])

TargetPreds = target_model.predict_classes(testdata)

adversaries = generate_adversary(black_box_model, testdata, to_categorical(TargetPreds, num_classes = 2))

TargetModelTestAccuracy = target_model.evaluate(testdata, testlabels)[1]
AdversarialTestAccuracy = target_model.evaluate(adversaries , testlabels)[1]

print("Target Model Accuracy: %s" % TargetModelTestAccuracy)
print("Adversary Test Accuracy: %s" % AdversarialTestAccuracy)
print("Accuracy Drop:", (TargetModelTestAccuracy - AdversarialTestAccuracy))


text:
# Copyright (c) Facebook, Inc. and its affiliates.

"""
Example command:
python fairmotion/tasks/clustering/generate_features.py \
    --type kinetic \
    --folder $PATH_WITH_BVH_FILES \
    --output-folder $PATH_TO_STORE_FEATURES \
    --cpu 10
"""

import argparse
import numpy as np
import os
from multiprocessing import Pool
from fairmotion.tasks.clustering.features import kinetic, manual
from fairmotion.tasks.clustering.features import utils as feat_utils
from fairmotion.data import bvh

thresholds = None

def extract_manual_features(motion):
    features = []
    f = manual.ManualFeatures(motion, feat_utils.Skeleton.PFNN)
    for _ in range(1, motion.num_frames(), 30):
        pose_features = []
        pose_features.append(
            f.f_nmove("neck", "rhip", "lhip", "rwrist", 1.8 * f.hl)
        )
        pose_features.append(
            f.f_nmove("neck", "lhip", "rhip", "lwrist", 1.8 * f.hl)
        )
        pose_features.append(
            f.f_nplane("chest", "neck", "neck", "rwrist", 0.2 * f.hl)
        )
        pose_features.append(
            f.f_nplane("chest", "neck", "neck", "lwrist", 0.2 * f.hl)
        )
        pose_features.append(
            f.f_move("belly", "chest", "chest", "rwrist", 1.8 * f.hl)
        )
        pose_features.append(
            f.f_move("belly", "chest", "chest", "lwrist", 1.8 * f.hl)
        )
        pose_features.append(
            f.f_angle("relbow", "rshoulder", "relbow", "rwrist", [0, 110])
        )
        pose_features.append(
            f.f_angle("lelbow", "lshoulder", "lelbow", "lwrist", [0, 110])
        )
        pose_features.append(
            f.f_nplane(
                "lshoulder", "rshoulder", "lwrist", "rwrist", 2.5 * f.sw
            )
        )
        pose_features.append(
            f.f_move("lwrist", "rwrist", "rwrist", "lwrist", 1.4 * f.hl)
        )
        pose_features.append(
            f.f_move("rwrist", "root", "lwrist", "root", 1.4 * f.hl)
        )
        pose_features.append(
            f.f_move("lwrist", "root", "rwrist", "root", 1.4 * f.hl)
        )
        pose_features.append(f.f_fast("rwrist", 2.5 * f.hl))
        pose_features.append(f.f_fast("lwrist", 2.5 * f.hl))
        pose_features.append(
            f.f_plane("root", "lhip", "ltoes", "rankle", 0.38 * f.hl)
        )
        pose_features.append(
            f.f_plane("root", "rhip", "rtoes", "lankle", 0.38 * f.hl)
        )
        pose_features.append(
            f.f_nplane("zero", "y_unit", "y_min", "rankle", 1.2 * f.hl)
        )
        pose_features.append(
            f.f_nplane("zero", "y_unit", "y_min", "lankle", 1.2 * f.hl)
        )
        pose_features.append(
            f.f_nplane("lhip", "rhip", "lankle", "rankle", 2.1 * f.hw)
        )
        pose_features.append(
            f.f_angle("rknee", "rhip", "rknee", "rankle", [0, 110])
        )
        pose_features.append(
            f.f_angle("lknee", "lhip", "lknee", "lankle", [0, 110])
        )
        pose_features.append(f.f_fast("rankle", 2.5 * f.hl))
        pose_features.append(f.f_fast("lankle", 2.5 * f.hl))
        pose_features.append(
            f.f_angle("neck", "root", "rshoulder", "relbow", [25, 180])
        )
        pose_features.append(
            f.f_angle("neck", "root", "lshoulder", "lelbow", [25, 180])
        )
        pose_features.append(
            f.f_angle("neck", "root", "rhip", "rknee", [50, 180])
        )
        pose_features.append(
            f.f_angle("neck", "root", "lhip", "lknee", [50, 180])
        )
        pose_features.append(
            f.f_plane("rankle", "neck", "lankle", "root", 0.5 * f.hl)
        )
        pose_features.append(
            f.f_angle("neck", "root", "zero", "y_unit", [70, 110])
        )
        pose_features.append(
            f.f_nplane("zero", "minus_y_unit", "y_min", "rwrist", -1.2 * f.hl)
        )
        pose_features.append(
            f.f_nplane("zero", "minus_y_unit", "y_min", "lwrist", -1.2 * f.hl)
        )
        pose_features.append(f.f_fast("root", 2.3 * f.hl))
        features.append(pose_features)
        f.next_frame()
    return features

def extract_kinetic_features(motion, thresholds, up_vec):
    features = kinetic.KineticFeatures(
        motion, 1 / motion.fps, thresholds, up_vec,
    )
    kinetic_feature_vector = []
    for i in range(motion.skel.num_joints()):
        positions_mean, positions_stddev = features.local_position_stats(i)
        feature_vector = np.hstack(
            [
                features.average_kinetic_energy_horizontal(i),
                features.average_kinetic_energy_vertical(i),
                features.average_energy_expenditure(i),
            ]
        )
        kinetic_feature_vector.extend(feature_vector)
    return kinetic_feature_vector

def extract_features(filepath, feature_type, thresholds=None, up_vec="z"):
    motion = bvh.load(filepath)
    if feature_type == "manual":
        return extract_manual_features(motion)
    else:
        return extract_kinetic_features(motion, thresholds, up_vec)

def wrapper_extract_features(inputs):
    global thresholds
    filepath = inputs[0]
    feature_type = inputs[1]
    up_vec = inputs[2]
    output_folder = inputs[3]

    features = extract_features(filepath, feature_type, thresholds, up_vec)
    filename = filepath.split("/")[-1]
    with open(
        os.path.join(output_folder, "features.tsv"), "a",
    ) as all_features:
        if feature_type == "manual":
            all_features.write(
                filename
                + ":"
                + "\t".join(
                    [
                        str(int(datapoint))
                        for pose_features in features
                        for datapoint in pose_features
                    ]
                )
                + "\n"
            )
        else:
            all_features.write(
                filename + ":" + "\t".join([str(f) for f in features]) + "\n"
            )
    print(f"Finished processing {filepath}")

def main(args):
    global thresholds
    thresholds = None
    os.makedirs(args.output_folder, exist_ok=True)
    pool = Pool(args.cpu)
    for root, _, files in os.walk(args.folder, topdown=False):
        pool.map(
            wrapper_extract_features,
            [
                (os.path.join(root, filename), args.type, args.up_vec, args.output_folder)
                for filename in files
            ],
        )

if __name__ == "__main__":
    parser = argparse.ArgumentParser(
        description="Extract features from BVH files"
    )
    parser.add_argument(
        "--type", nargs="?", choices=["manual", "kinetic"], default="manual",
    )
    parser.add_argument(
        "--folder",
        type=str,
        help="Folder with files to extract features from",
        required=True,
    )
    parser.add_argument("--output-folder", type=str, required=True)
    parser.add_argument("--up-vec", nargs="?", choices=["y", "z"], default="z")
    parser.add_argument("--cpu", type=int, default=40)
    args = parser.parse_args()
    main(args)

max_stars_repo_name: bmahlbrand/fairmotion

max_issues_repo_name: bmahlbrand/fairmotion

max_forks_repo_name: bmahlbrand/fairmotion


text:
"""
* affine triangulation
* display point clouds
* DEM projection

Copyright (C) 2018, Enric Meinhardt-Llopis <enric.meinhardt@cmla.ens-cachan.fr>
Copyright (C) 2018, Carlo de Franchis <carlo.de-franchis@ens-cachan.fr>
Copyright (C) 2018, Gabriele Facciolo <facciolo@cmla.ens-cachan.fr>
"""

import numpy as np

def triangulation_affine(PA, PB, x1, y1, x2, y2):
    """
    Triangulate a (list of) match(es) between two images of affine cameras.

    Arguments:
        PA, PB : affine (projection) camera matrices of the two images
        x1, y1 : pixel coordinates in the domain of the first image
        x2, y2 : pixel coordinates in the domain of the second image

    Return value: a 4-tuple (lon, lat, h, e)
        lon, lat, h, e : coordinates of the 3D point(s), reprojection error
    """

    # build projection and localization matrices as 4x4 homogeneous maps
    # (x,y,h,1) <-> (lon,lat,h,1)
    PA = np.vstack([ PA[0:2], [0,0,1,0], [0,0,0,1]])  # pick only first two rows
    PB = np.vstack([ PB[0:2], [0,0,1,0], [0,0,0,1]])  # pick only first two rows
    LA = np.linalg.inv(PA)  # inverse of a 4x4 matrix

    # affine epipolar correspondence
    E = PB @ LA

    # Now, the linear system E * (x1, y1, h, 1) = (x2, y2, h, 1)
    # has two different equations and one unknown h.  We solve it by
    # least squares (a simple projection, in this case).

    # give names to the 8 non-trivial coefficients of E
    a, b, p, r = E[0]
    c, d, q, s = E[1]

    # coefficients of the affine triangulation
    f = [-p*a - q*c, -p*b - q*d, p, q, -p*r - q*s ] / (p*p + q*q)

    # apply the triangulation (first use of the input points)
    h = f[0]*x1 + f[1]*y1 + f[2]*x2 + f[3]*y2 + f[4]

    # finish the computation and return the 4 required numbers (or vectors)
    lon = LA[0,0] * x1 + LA[0,1] * y1 + LA[0,2] * h + LA[0,3]
    lat = LA[1,0] * x1 + LA[1,1] * y1 + LA[1,2] * h + LA[1,3]
    ex = E[0,0] * x1 + E[0,1] * y1 + E[0,2] * h + E[0,3] - x2
    ey = E[1,0] * x1 + E[1,1] * y1 + E[1,2] * h + E[1,3] - y2
    e = ex * ex + ey * ey
    return lon, lat, h, e

def triangulate_disparities(dmap, rpc1, rpc2, S1, S2, PA, PB,):
    """
    Triangulate a disparity map

    Arguments:
        dmap : a disparity map between two rectified images
        rpc1, rpc2 : calibration data of each image
        S1, S2 : rectifying affine maps (from the domain of the fullsize images)
        PA, PB : the affine approximations of rpc1 and rpc2 (not always used)

    Return:
        xyz : a matrix of size Nx3 (where N is the number of finite disparites
              in dmap) this matrix contains the coordinates of the 3d points in
              "lon,lat,h" or "e,n,h"
    """
    from utils import utm_from_lonlat

    # 1. unroll all the valid (finite) disparities of dmap into a vector
    m = np.isfinite(dmap.flatten())
    x = np.argwhere(np.isfinite(dmap))[:,1]  # attention to order of the indices
    y = np.argwhere(np.isfinite(dmap))[:,0]
    d = dmap.flatten()[m]

    # 2. for each disparity

    # 2.1. produce a pair in the original domain by composing with S1 and S2
    p = np.linalg.inv(S1) @ np.vstack( (x+0, y, np.ones(len(d))) )
    q = np.linalg.inv(S2) @ np.vstack( (x+d, y, np.ones(len(d))) )

    # 2.2. triangulate the pair of image points to find a 3D point (in UTM)
    lon, lat, h, e = triangulation_affine(PA,PB, p[0,:],p[1,:], q[0,:],q[1,:])
    #east, north, _, _ = np.vectorize(utm.from_latlon)(lat, lon)
    east, north = utm_from_lonlat(lon, lat)

    # 2.3. add this point to the output list
    xyz = np.vstack((east, north, h)).T

    return xyz

# code for projection into UTM coordinates
from numba import jit

# utility function for computing the average at each cell
@jit
def reduceavg(w, h,  ix, iy, z):
    D_sum = np.zeros((h,w))  #-np.ones((h,w))*np.inf
    D_cnt = np.zeros((h,w))
    for t in range(len(ix)):
        if iy[t]<0 or ix[t]<0 or iy[t]>=h or ix[t]>=w:
            continue
        ty = iy[t]
        tx = ix[t]
        D_sum[ty,tx] += z[t]
        D_cnt[ty,tx] += 1

    D_sum /= D_cnt*(D_cnt>0) + 1*(D_cnt==0) # needed for computing average

    return D_sum

# projects them into a grid defined by the bunding box (emin, emax, nmin, nmax)
def project_cloud_into_utm_grid(xyz, emin, emax, nmin, nmax, resolution=1):
    """
    Project a point cloud into an utm grid to produce a DEM
    The algorithm averages all the points that fall into each square of the grid

    Arguments:
        xyz : a Nx3 matrix representing a point cloud in (lon,lat,h) coordinates
        emin,emax,nmin,nmax : a bounding box in UTM coordinates
        resolution : the target resolution in meters (by default, 1 meter)

    Return:
        dem : a 2D array of heights in meters
    """

    # width and height of the image domain
    w = int(np.ceil((emax - emin)/resolution))
    h = int(np.ceil((nmax - nmin)/resolution))

    # extract and quantize columns
    x = xyz[:,0]
    y = xyz[:,1]
    z = xyz[:,2]

    ix = np.asarray((x - emin)/resolution, dtype="int")
    iy = np.asarray((nmax - y)/resolution, dtype="int")

    # per-cell average
    dem = reduceavg (w,h,  ix, iy, z )

    return dem

reason: import numpy,from numba


text:
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
"""
Interleaved RB analysis class.
"""
from typing import List, Union

import numpy as np

import qiskit_experiments.curve_analysis as curve
from qiskit_experiments.framework import AnalysisResultData, FitVal
from .rb_analysis import RBAnalysis

class InterleavedRBAnalysis(RBAnalysis):
    r"""A class to analyze interleaved randomized benchmarking experiment.

    # section: overview
        This analysis takes only two series for standard and interleaved RB curve fitting.
        From the fit :math:`\alpha` and :math:`\alpha_c` value this analysis estimates
        the error per Clifford (EPC) of the interleaved gate.

        The EPC estimate is obtained using the equation

        .. math::

            r_{\mathcal{C}}^{\text{est}} =
                \frac{\left(d-1\right)\left(1-\alpha_{\overline{\mathcal{C}}}/\alpha\right)}{d}

        The systematic error bounds are given by

        .. math::

            E = \min\left\{
                \begin{array}{c}
                    \frac{\left(d-1\right)\left[\left|\alpha-\alpha_{\overline{\mathcal{C}}}\right|
                    +\left(1-\alpha\right)\right]}{d} \\
                    \frac{2\left(d^{2}-1\right)\left(1-\alpha\right)}
                    {\alpha d^{2}}+\frac{4\sqrt{1-\alpha}\sqrt{d^{2}-1}}{\alpha}
                \end{array}
            \right.

        See Ref. [1] for more details.

    # section: fit_model
        The fit is based on the following decay functions:

        Fit model for standard RB

        .. math::

            F(x) = a \alpha^{x} + b

        Fit model for interleaved RB

        .. math::

            F(x) = a (\alpha_c \alpha)^{x_2} + b

    # section: fit_parameters
        defpar a:
            desc: Height of decay curve.
            init_guess: Determined by the average :math:`a` of the standard and interleaved RB.
            bounds: [0, 1]
        defpar b:
            desc: Base line.
            init_guess: Determined by the average :math:`b` of the standard and interleaved RB.
                Usually equivalent to :math:`(1/2)^n` where :math:`n` is number of qubit.
            bounds: [0, 1]
        defpar \alpha:
            desc: Depolarizing parameter.
            init_guess: Determined by the slope of :math:`(y - b)^{-x}` of the first and the
                second data point of the standard RB.
            bounds: [0, 1]
        defpar \alpha_c:
            desc: Ratio of the depolarizing parameter of interleaved RB to standard RB curve.
            init_guess: Estimate :math:`\alpha' = \alpha_c \alpha` from the
                interleaved RB curve, then divide this by the initial guess of :math:`\alpha`.
            bounds: [0, 1]

    # section: reference
        .. ref_arxiv:: 1 1203.4550

    """

    __series__ = [
        curve.SeriesDef(
            name="Standard",
            fit_func=lambda x, a, alpha, alpha_c, b: curve.fit_function.exponential_decay(
                x, amp=a, lamb=-1.0, base=alpha, baseline=b
            ),
            filter_kwargs={"interleaved": False},
            plot_color="red",
            plot_symbol=".",
            plot_fit_uncertainty=True,
            model_description=r"a \alpha^{x} + b",
        ),
        curve.SeriesDef(
            name="Interleaved",
            fit_func=lambda x, a, alpha, alpha_c, b: curve.fit_function.exponential_decay(
                x, amp=a, lamb=-1.0, base=alpha * alpha_c, baseline=b
            ),
            filter_kwargs={"interleaved": True},
            plot_color="orange",
            plot_symbol="^",
            plot_fit_uncertainty=True,
            model_description=r"a (\alpha_c\alpha)^{x} + b",
        ),
    ]

    @classmethod
    def _default_options(cls):
        """Default analysis options."""
        default_options = super()._default_options()
        default_options.result_parameters = ["alpha", "alpha_c"]
        return default_options

    def _generate_fit_guesses(
        self, user_opt: curve.FitOptions
    ) -> Union[curve.FitOptions, List[curve.FitOptions]]:
        """Compute the initial guesses.

        Args:
            user_opt: Fit options filled with user provided guess and bounds.

        Returns:
            List of fit options that are passed to the fitter function.
        """
        user_opt.bounds.set_if_empty(
            a=(0, 1),
            alpha=(0, 1),
            alpha_c=(0, 1),
            b=(0, 1),
        )

        # for standard RB curve
        std_curve = self._data(series_name="Standard")
        opt_std = user_opt.copy()
        opt_std = self._initial_guess(opt_std, std_curve.x, std_curve.y, self._num_qubits)

        # for interleaved RB curve
        int_curve = self._data(series_name="Interleaved")
        opt_int = user_opt.copy()
        if opt_int.p0["alpha_c"] is not None:
            opt_int.p0["alpha"] = opt_std.p0["alpha"] * opt_int.p0["alpha_c"]
        opt_int = self._initial_guess(opt_int, int_curve.x, int_curve.y, self._num_qubits)

        user_opt.p0.set_if_empty(
            a=np.mean([opt_std.p0["a"], opt_int.p0["a"]]),
            alpha=opt_std.p0["alpha"],
            alpha_c=min(opt_int.p0["alpha"] / opt_std.p0["alpha"], 1),
            b=np.mean([opt_std.p0["b"], opt_int.p0["b"]]),
        )

        return user_opt

    def _extra_database_entry(self, fit_data: curve.FitData) -> List[AnalysisResultData]:
        """Calculate EPC."""
        nrb = 2 ** self._num_qubits
        scale = (nrb - 1) / nrb

        alpha = fit_data.fitval("alpha")
        alpha_c = fit_data.fitval("alpha_c")

        # Calculate epc_est (=r_c^est) - Eq. (4):
        epc = FitVal(value=scale * (1 - alpha_c.value), stderr=scale * alpha_c.stderr)

        # Calculate the systematic error bounds - Eq. (5):
        systematic_err_1 = scale * (abs(alpha.value - alpha_c.value) + (1 - alpha.value))
        systematic_err_2 = (
            2 * (nrb * nrb - 1) * (1 - alpha.value) / (alpha.value * nrb * nrb)
            + 4 * (np.sqrt(1 - alpha.value)) * (np.sqrt(nrb * nrb - 1)) / alpha.value
        )
        systematic_err = min(systematic_err_1, systematic_err_2)
        systematic_err_l = epc.value - systematic_err
        systematic_err_r = epc.value + systematic_err

        extra_data = AnalysisResultData(
            name="EPC",
            value=epc,
            chisq=fit_data.reduced_chisq,
            quality=self._evaluate_quality(fit_data),
            extra={
                "EPC_systematic_err": systematic_err,
                "EPC_systematic_bounds": [max(systematic_err_l, 0), systematic_err_r],
            },
        )

        return [extra_data]

max_stars_repo_name: jakelishman/qiskit-experiments

max_issues_repo_name: jakelishman/qiskit-experiments

max_forks_repo_name: jakelishman/qiskit-experiments


text:
import os
import pickle
import struct
import sys
import tarfile

import numpy as np
import requests
from scipy import io
import argparse

import tensorflow as tf
import cv2
import scipy.io
import numpy as np
from PIL import Image

def download_file(url, save_path):
    with open(save_path, 'wb') as f:
        response = requests.get(url, stream=True)
        total_length = response.headers.get('content-length')

        print('Downloading file of size {} MB.'.format(str(long(total_length) / (1024 * 1024))))

        if total_length is None:
            f.write(response.content)
        else:
            dl = 0
            total_length = int(total_length)
            for data in response.iter_content(chunk_size=4096):
                dl += len(data)
                f.write(data)
                done = int(50 * dl / total_length)
                sys.stdout.write('\r[{}{}]'.format('=' * done, ' ' * (50 - done)))
                sys.stdout.flush()
    sys.stdout.write('\n')

def download_svhn_dataset(save_dir):
    train_path = os.path.join(save_dir, 'train_32x32.mat')
    test_path = os.path.join(save_dir, 'test_32x32.mat')

    # print 'Downloading Svhn Train...'
    # download_file('http://ufldl.stanford.edu/housenumbers/train_32x32.mat', train_path)
    #
    # print 'Downloading Svhn Test...'
    # download_file('http://ufldl.stanford.edu/housenumbers/test_32x32.mat', test_path)

    # np.set_printoptions(threshold=np.nan)

    train_mat = scipy.io.loadmat(train_path)
    x_train = np.array(train_mat["X"]).transpose()
    y_train = np.array(train_mat["y"])

    test_mat = scipy.io.loadmat(test_path)
    x_test = np.array(test_mat["X"]).transpose()
    y_test = np.array(test_mat["y"])

    print("""
            SVHN dataset summary

            x_train shape   :   {}
            y_train shape   :   {}
            x_test shape    :   {}
            x_test shape    :   {}
            """.format(x_train.shape, y_train.shape, x_test.shape, y_test.shape))

    file_num = 0

    # for i in range(train_y[0:]:
    #     img = Image.fromarray(train_x[:, :, :, i], 'RGB')
    #     img.save("svhn" + str(i) + ".jpg")
    #     print "saving: " + str(i)
    #
    #     break

def create_mnist_dataset(save_dir):
    (x_train, y_train), (x_test, y_test) = tf.keras.datasets.mnist.load_data()

    print("""
        MNIST dataset summary

        x_train shape   :   {}
        y_train shape   :   {}
        x_test shape    :   {}
        x_test shape    :   {}
        """.format(x_train.shape, y_train.shape, x_test.shape, y_test.shape))

    file_num = 0

    for x, y in zip(x_train, y_train):
        current_file_path = os.path.join(os.path.join(save_dir, 'train'), str(y))

        if not os.path.exists(current_file_path):
            os.makedirs(current_file_path)

        cv2.imwrite('{}.png'.format(os.path.join(current_file_path, str(file_num))), x)
        file_num = file_num + 1

    print('Saved training dataset. Started saving test dataset.')

    for x, y in zip(x_test, y_test):
        current_file_path = os.path.join(os.path.join(save_dir, 'test'), str(y))

        if not os.path.exists(current_file_path):
            os.makedirs(current_file_path)

        cv2.imwrite('{}.png'.format(os.path.join(current_file_path, str(file_num))), x)
        file_num = file_num + 1

    print('Saved test dataset.')

def download_cifar10_dataset(save_dir):
    (x_train, y_train), (x_test, y_test) = tf.keras.datasets.cifar10.load_data()

    print("""
            CIFAR10 dataset summary

            x_train shape   :   {}
            y_train shape   :   {}
            x_test shape    :   {}
            x_test shape    :   {}
            """.format(x_train.shape, y_train.shape, x_test.shape, y_test.shape))

    file_num = 0

    for x, y in zip(x_train, y_train):
        current_file_path = os.path.join(os.path.join(save_dir, 'train'), str(y[0]))

        if not os.path.exists(current_file_path):
            os.makedirs(current_file_path)

        cv2.imwrite('{}.png'.format(os.path.join(current_file_path, str(file_num))), x)
        file_num = file_num + 1

    print('Saved training dataset. Started saving test dataset.')

    for x, y in zip(x_test, y_test):
        current_file_path = os.path.join(os.path.join(save_dir, 'test'), str(y[0]))

        if not os.path.exists(current_file_path):
            os.makedirs(current_file_path)

        cv2.imwrite('{}.png'.format(os.path.join(current_file_path, str(file_num))), x)
        file_num = file_num + 1

    print('Saved test dataset.')
#
#
# def save_dataset_as_images(x_train, y_train, x_test, y_test, save_dir):
#     train_dir = os.path.join(save_dir, 'training')
#     test_dir = os.path.join(save_dir, 'test')
#
#     print 'Labels train', np.unique(y_train)
#     print 'Labels test', np.unique(y_test)
#
#     x_train = x_train.astype(float)
#     y_train = y_train.astype(float)
#     x_test = x_test.astype(float)
#     y_test = y_test.astype(float)
#
#
# def delete_irrelavent_files(save_dir):
#     pass

def create_dataset(dataset_name, save_dir):
    if dataset_name == 'cifar10':
        download_cifar10_dataset(save_dir)
    elif dataset_name == 'svhn':
        download_svhn_dataset(save_dir)
    elif dataset_name == 'mnist':
        create_mnist_dataset(save_dir)
    else:
        print('Cannot download the dataset \'{}\', please dowwnload it manually!')

    print("Creating dataset completed!!")

if __name__ == '__main__':
    parser = argparse.ArgumentParser()

    dataset_options = ['cifar10', 'svhn', 'mnist']

    parser.add_argument('-dataset', choices=dataset_options, help='Which data set')
    parser.add_argument('-save_folder', help='Where to save images')

    args = parser.parse_args()

    save_folder = os.path.join(args.save_folder, args.dataset)

    os.system('rm -rf "{}/train"'.format(save_folder))
    os.system('rm -rf "{}/test"'.format(save_folder))

    if not os.path.exists(save_folder):
        os.makedirs(save_folder)

    create_dataset(args.dataset, save_folder)

# import tensorflow as tf
# import cv2
#
# (x_train, y_train), (x_test, y_test) = tf.cnn.datasets.mnist.load_data()
#
# print("""
#     x_train shape   :   {}
#     y_train shape   :   {}
#     x_test shape    :   {}
#     x_test shape    :   {}
#     """.format(x_train.shape, y_train.shape, x_test.shape, y_test.shape))
#
# file_num = 0
# for x,y in zip(x_train,y_train):
#     print(x.shape)
#     print(y)
#     cv2.imwrite('{}_{}.png'.format(y,file_num), x)
#     file_num = file_num + 1
#
#     break

# import cv2
#
# cv2.imwrite('out.png', data_point)

reason: import numpy,import scipy,from scipy


text:
import pandas as pd
import numpy as np
import itertools
import argparse

class AMS:

"""
Author: Tania Lopez-Cantu
E-mail: tlopez@andrew.cmu.edu

#############################
Purpose

This class contains methods to clean a timeseries, and an implementation of
the fixed maxima and sliding
maxima algorithms provided in

Papalexiou, S. M., Dialynas, Y. G., & Grimaldi, S. (2016).
Hershfield factor revisited: Correcting annual maximum precipitation.
Journal of Hydrology, 542, 884–895. https://doi.org/10.1016/j.jhydrol.2016.09.058

These two functions extract the annual maximum precipitation (AMS)
from a precipitation time series. The two approaches arise from the need to account for the
fact that precipitation is systematicall recorded. For example, at a meteorological station,
someone will check the tipping bucket pluviographs at some fixed local time each day which is
then recorded as a "daily rainfall time series" at a particular location. This case results
in "fixed" records, but "fixed" records have been shown inappropriate for
estimating rainfall maxima. Because rainfall is a continuous variable, discretizing it can result
in biases when estimating extreme rainfall, so it is advised to estimate
annual maximum series using the sliding maxima approach.

"""

    def __init__(self, *args, **kwargs):

        pd.read_csv.__init__(self, *args, **kwargs)

    def reformat(self, path):
        ts = pd.read_csv(path, index_col=0, parse_dates=['date'])
        self.reformat = ts.drop(['qflags', 'mflags'], axis=1)
        return self.reformat

    def sliding_max(self, k, data):
        """
        Function to extract AMS from a time series using the sliding maxima approach.

        Parameters
        ----------
        Input:
            k: int, duration over which AMS will be computed. For example, if the data is at an hourly resolution, and
                    want to compute the 3-hour AMS, then k=3.
            data: DataFrame, rainfall time series in a pandas two column dataframe format. One column should be the "date"
                    of the record (i.e. "1960-05-24 00:00:00" if hourly records), and the second must be the rainfall values.
        """
        tp = data.values
        period = 24*365
        agg_values = []
        start_j = 1
        end_j = k*int(np.floor(period/k))
        for j in range(start_j, end_j + 1):
            start_i = j - 1
            end_i = j + k + 1
            agg_values.append(np.nansum(tp[start_i:end_i]))
        self.sliding_max = max(agg_values)
        return self.sliding_max

    def fixed_max(self, k, data):
        """
        Function to extract AMS from a time series using the fixed maxima approach.

        Parameters
        ----------
        Input:
            k: int, duration over which AMS will be computed. For example, if the data is at an hourly resolution, and
                    want to compute the 3-hour AMS, then k=3.
            data: DataFrame, rainfall time series in a pandas two column dataframe format. One column should be the "date"
                    of the record (i.e. "1960-05-24 00:00:00" if hourly records), and the second must be the rainfall values.
        """
        tp = data.values
        period = 24*365
        agg_values = []
        start_j = 1
        end_j = int(np.floor(period/k))
        for j in range(start_j, end_j + 1):
            start_i = (j - 1)*k
            end_i = (j*k)
            agg_values.append(np.nansum(tp[start_i:end_i]))
        self.fixed_max = maxx(agg_values)
        return self.fixed_max

max_stars_repo_name: tplopez/historical-IDF

max_issues_repo_name: tplopez/historical-IDF

max_forks_repo_name: tplopez/historical-IDF