text:
import numpy as np
import pytest

from grappa.kernel_estimation import _geometry_kernel_estimation

@pytest.mark.parametrize('ac, ny, expected_kernel', [
    (  # first basic example with a kernel using only the 2 direct neighbours
        [
            [1, 2, 1],
            [2, 4, 2],
            [3, 6, 3],
            [7, 2, -5],
        ],
        1,
        [[1]*2],
    ),
    (  # second example, a bit more complex, still in single geometry
        [
            [1, 2, 1],
            [2, 4, 2],
            [3, 6, 3],
            [7, 2, -5],
            [0.1, 1.1, 1],
            [3.1, 5.1, 2],
            [3.1, 12.1, 9],
            [7.1, 3.1, -4],
        ],
        3,
        [[0, 0, 1, 1, 0, 0]],
    ),
    (  # third example to test the complex part of kernel estimation
        [
            [1, 1 + 1j, 1],
            [2, 2 + 2j, 2],
            [3, 3 + 3j, 3],
            [7, 7 - 5j, -5],
        ],
        1,
        [[1, 1j]],
    ),
])
def test_geometry_kernel_estimation_singlecoil_simple(ac, ny, expected_kernel):
    ac = np.array(ac)[None, :]
    expected_kernel = np.array(expected_kernel)
    grappa_kernel = _geometry_kernel_estimation(
        ac,
        i_geom=0,
        ny=ny,
        n_geometries=1,
        ncoils=1,
        lamda=1e-8,
    )
    np.testing.assert_array_almost_equal(expected_kernel, grappa_kernel, decimal=5)

@pytest.mark.parametrize('ac, ny, i_geom, expected_kernel', [
    (  # first basic example with a kernel using only the 2 direct sampled neighbours
        [
            [1, 2, 1, 1],
            [2, 4, 2, 2],
            [3, 6, 3, 3],
            [7, 2, -5, -5],
        ],
        1,
        0,
        [[1]*2],
    ),
    (  # first basic example with a kernel using only the 2 direct sampled neighbours
        [
            [1, 2, 1, 1],
            [2, 4, 2, 2],
            [3, 6, 3, 3],
            [7, 2, -5, -5],
        ],
        1,
        1,
        [[0, 1]],
    ),
    (  # second example, a bit more complex
        [
            [1, 2, 1, 1],
            [2, 4, 2, 2],
            [3, 6, 3, 3],
            [7, 2, -5, -5],
            [0.1, 1.1, 1, 1],
            [3.1, 5.1, 2, 2],
            [3.1, 12.1, 9, 9],
            [7.1, 3.1, -4, -4],
        ],
        3,
        0,
        [[0, 0, 1, 1, 0, 0]],
    )
])
def test_geometry_kernel_estimation_singlecoil_double(ac, ny, i_geom, expected_kernel):
    ac = np.array(ac)[None, :]
    expected_kernel = np.array(expected_kernel)
    grappa_kernel = _geometry_kernel_estimation(
        ac,
        i_geom=i_geom,
        ny=ny,
        n_geometries=2,
        ncoils=1,
        lamda=1e-8,
    )
    np.testing.assert_array_almost_equal(expected_kernel, grappa_kernel, decimal=5)

@pytest.mark.parametrize('ac, ny, expected_kernel', [
    (  # first basic example with a kernel using only the 2 direct neighbours
       # from the first coil (either positively or negatively)
        [
            [
                [1, 2, 1],
                [2, 4, 2],
                [3, 6, 3],
                [7, 2, -5],
            ],
            [
                [0.001, -2, 0.0001],
                [0.004, -4, 0.002],
                [-0.001, -6, 0.0002],
                [0.002, -2, 0.001],
            ],
        ],
        1,
        [[1, 1, 0, 0], [-1, -1, 0, 0]],
    ),
    (  # second example, a bit more complex, still in single geometry
        [
            [
                [1, 2, 1],
                [2, 4, 2],
                [3, 6, 3],
                [7, 2, -5],
                [0.1, 1.1, 1],
                [3.1, 5.1, 2],
                [3.1, 12.1, 9],
                [7.1, 3.1, -4],
                [8, 3.5, -4.5],
                [18, 12, -6],
                [-3, 2, 5],
                [-4.2, 9, 13.2],
            ],
            [
                [0.01, -2, 0.0012],
                [0.02, -4, 0.0012],
                [0.001, -6, 0.0011],
                [-0.01, -2, 0.0013],
                [0.003, -1.1, 0.0014],
                [0.005, -5.1, 0.0013],
                [-0.02, -12.1, 0.0019],
                [0.007, -3.1, 0.0013],
                [0.002, -3.5, 0.001],
                [0.004, -12, 0.001],
                [0.00036, -2, 0.003],
                [-0.002, -9, -0.003],
            ],
        ],
        3,
        [
            [0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
            [0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0],
        ],
    )
])
def test_geometry_kernel_estimation_multicoil_simple(ac, ny, expected_kernel):
    ac = np.array(ac)
    expected_kernel = np.array(expected_kernel)
    grappa_kernel = _geometry_kernel_estimation(
        ac,
        i_geom=0,
        ny=ny,
        n_geometries=1,
        ncoils=2,
        lamda=1e-6,
    )
    np.testing.assert_array_almost_equal(expected_kernel, grappa_kernel, decimal=3)

@pytest.mark.parametrize('ac, ny, i_geom, expected_kernel', [
    (  # first basic example with a kernel using only the 2 direct neighbours
       # from the first coil (either positively or negatively)
        [
            [
                [1, 2, 1, 1],
                [2, 4, 2, 2],
                [3, 6, 3, 3],
                [7, 2, -5, -5],
            ],
            [
                [0.001, -2, 0, 0],
                [0.004, -4, 0, 0],
                [-0.001, -6, 0, 0],
                [0.002, -2, 0, 0],
            ],
        ],
        1,
        0,
        [[1, 1, 0, 0], [-1, -1, 0, 0]],
    ),
])
def test_geometry_kernel_estimation_multicoil_double(ac, ny, i_geom, expected_kernel):
    ac = np.array(ac)
    expected_kernel = np.array(expected_kernel)
    grappa_kernel = _geometry_kernel_estimation(
        ac,
        i_geom=i_geom,
        ny=ny,
        n_geometries=2,
        ncoils=2,
    )
    np.testing.assert_array_almost_equal(expected_kernel, grappa_kernel, decimal=3)

max_stars_repo_name: zaccharieramzi/grappa

max_issues_repo_name: zaccharieramzi/grappa

max_forks_repo_name: zaccharieramzi/grappa


text:
# -*- coding: utf-8 -*-
"""
Generators for random intersection graphs.
"""
#    Copyright (C) 2011 by
#    Aric Hagberg <hagberg@lanl.gov>
#    Dan Schult <dschult@colgate.edu>
#    Pieter Swart <swart@lanl.gov>
#    All rights reserved.
#    BSD license.
import random
import networkx as nx
__author__ = "\n".join(['Aric Hagberg (hagberg@lanl.gov)'])

__all__ = ['uniform_random_intersection_graph',
           'k_random_intersection_graph',
           'general_random_intersection_graph',
           ]

def uniform_random_intersection_graph(n, m, p, seed=None):
    """Return a uniform random intersection graph.

    Parameters
    ----------
    n : int
        The number of nodes in the first bipartite set (nodes)
    m : int
        The number of nodes in the second bipartite set (attributes)
    p : float
        Probability of connecting nodes between bipartite sets
    seed : int, optional
        Seed for random number generator (default=None).

    See Also
    --------
    gnp_random_graph

    References
    ----------
    .. [1] K.B. Singer-Cohen, Random Intersection Graphs, 1995,
       PhD thesis, Johns Hopkins University
    .. [2] Fill, J. A., Scheinerman, E. R., and Singer-Cohen, K. B.,
       Random intersection graphs when m = !(n):
       An equivalence theorem relating the evolution of the g(n, m, p)
       and g(n, p) models. Random Struct. Algorithms 16, 2 (2000), 156–176.
    """
    G=nx.bipartite_random_graph(n, m, p, seed=seed)
    return nx.projected_graph(G, range(n))

def k_random_intersection_graph(n,m,k):
    """Return a intersection graph with randomly chosen attribute sets for
    each node that are of equal size (k).

    Parameters
    ----------
    n : int
        The number of nodes in the first bipartite set (nodes)
    m : int
        The number of nodes in the second bipartite set (attributes)
    k : float
        Size of attribute set to assign to each node.
    seed : int, optional
        Seed for random number generator (default=None).

    See Also
    --------
    gnp_random_graph, uniform_random_intersection_graph

    References
    ----------
    .. [1] Godehardt, E., and Jaworski, J.
       Two models of random intersection graphs and their applications.
       Electronic Notes in Discrete Mathematics 10 (2001), 129--132.
    """
    G = nx.empty_graph(n + m)
    mset = range(n,n+m)
    for v in range(n):
        targets = random.sample(mset, k)
        G.add_edges_from(zip([v]*len(targets), targets))
    return nx.projected_graph(G, range(n))

def general_random_intersection_graph(n,m,p):
    """Return a random intersection graph with independent probabilities
    for connections between node and attribute sets.

    Parameters
    ----------
    n : int
        The number of nodes in the first bipartite set (nodes)
    m : int
        The number of nodes in the second bipartite set (attributes)
    p : list of floats of length m
        Probabilities for connecting nodes to each attribute
    seed : int, optional
        Seed for random number generator (default=None).

    See Also
    --------
    gnp_random_graph, uniform_random_intersection_graph

    References
    ----------
    .. [1] Nikoletseas, S. E., Raptopoulos, C., and Spirakis, P. G.
       The existence and efficient construction of large independent sets
       in general random intersection graphs. In ICALP (2004), J. D´ıaz,
       J. Karhum¨aki, A. Lepist¨o, and D. Sannella, Eds., vol. 3142
       of Lecture Notes in Computer Science, Springer, pp. 1029–1040.
    """
    if len(p)!=m:
        raise ValueError("Probability list p must have m elements.")
    G = nx.empty_graph(n + m)
    mset = range(n,n+m)
    for u in range(n):
        for v,q in zip(mset,p):
            if random.random()<q:
                G.add_edge(u,v)
    return nx.projected_graph(G, range(n))

ext: py

max_stars_repo_name: kushgrover/apt-vs-dift

max_issues_repo_name: kushgrover/apt-vs-dift

max_forks_repo_name: kushgrover/apt-vs-dift

reason: import networkx


text:
#!/usr/bin/env python2

import os
import argparse

import numpy as np
from matplotlib import rc, cm, colors
import matplotlib.pyplot as plt
from colorama import init, Fore
from analyze_rrt import kSecPerOutIter

init(autoreset=True)

rc('font', **{'serif': ['Cardo'], 'size': 20})
rc('text', usetex=True)

def readRRTDataIter(verts_fn, edges_fn):
    verts = np.loadtxt(verts_fn)
    assert verts.shape[1] == 3
    verts_converted = np.zeros(verts.shape)
    verts_converted[:, 0] = verts[:, 1]
    verts_converted[:, 1] = -verts[:, 0]
    connected_vert_ids = np.loadtxt(edges_fn)
    assert connected_vert_ids.shape[1] == 2
    edges_start = np.array([verts_converted[int(idx[0]), 0:2]
                            for idx in connected_vert_ids])
    edges_end = np.array([verts_converted[int(idx[1]), 0:2]
                          for idx in connected_vert_ids])

    return verts_converted, edges_start, edges_end

if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument('--res_dir', required=True)
    parser.add_argument('--out_dir', default=None, type=str)
    parser.add_argument('--plot_iter', type=int, default=10)
    parser.add_argument('--plot_iter_step', type=int, default=2)
    parser.add_argument("--points3d", required=True)
    parser.add_argument('--max_z', default=5, type=float)
    parser.add_argument('--no_legend', action='store_false', dest='legend')
    parser.add_argument('--colorbar', action='store_true', dest='colorbar')
    parser.set_defaults(legend=True, colorbar=False)
    args = parser.parse_args()

    out_dir = args.out_dir if args.out_dir else args.res_dir
    nm = os.path.basename(args.res_dir)

    rrt_edges_nms = sorted([v for v in os.listdir(
        args.res_dir) if v.startswith('rrt_edges_iter_')])
    rrt_edges_fns = [os.path.join(args.res_dir, v) for v in rrt_edges_nms]
    rrt_verts_nms = sorted([v for v in os.listdir(
        args.res_dir) if v.startswith('rrt_verts_iter_')])
    rrt_verts_fns = [os.path.join(args.res_dir, v) for v in rrt_verts_nms]

    n_iter = len(rrt_edges_fns)
    assert len(rrt_verts_fns) == n_iter

    print(Fore.YELLOW + "Found RRT data of {} iterations.".format(n_iter))

    print("Read RRT data...")
    verts = []
    edges_start = []
    edges_end = []
    for iter in range(args.plot_iter):
        verts_i, edge_start_i, edge_end_i = readRRTDataIter(
            rrt_verts_fns[iter], rrt_edges_fns[iter])
        print("- iter {}: {} vert.,  {} edges".format(
            iter, verts_i.shape[0], edge_start_i.shape[0]))
        verts.append(verts_i)
        edges_start.append(edge_start_i)
        edges_end.append(edge_end_i)

    print("Load and preprocess landmarks...")
    raw_pts3d = np.loadtxt(args.points3d)
    pts3d = []
    for idx in range(raw_pts3d.shape[0]):
        pt_i = raw_pts3d[idx]
        if pt_i[2] > args.max_z:
            continue
        pts3d.append([pt_i[1], -pt_i[0], pt_i[2]])
    pts3d = np.array(pts3d)

    print("Plot...")
    fig = plt.figure(figsize=(8, 5))
    ax = fig.add_subplot(111)
    ax.scatter(pts3d[:, 0], pts3d[:, 1], c='gray', s=2.0, alpha=0.5)
    ax.scatter([], [], c='gray', s=10.0, label='Landmarks')
    ax.set_xticks([])
    ax.set_yticks([])
    ax.axis('off')
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()

    cmap = cm.get_cmap('jet')
    for iter in np.arange(0, args.plot_iter, args.plot_iter_step):
        print('- iter {}: {} verts'.format(iter, verts[iter].shape))
        c = cmap(iter * 1.0 / (args.plot_iter - 1))
        ax.scatter(verts[iter][:, 0], verts[iter]
                   [:, 1], c=c, s=0.1, zorder=-iter)

    label = 'RRT* Vertices'
    dummy_d = 5 * args.plot_iter
    dummy_offset = 1
    dummy_x = np.arange(xlim[0] - dummy_offset -
                        dummy_d, xlim[0] - dummy_offset + 1e-5, 5)
    dummy_y = np.arange(ylim[0] - dummy_offset -
                        dummy_d, ylim[0] - dummy_offset + 1e-5, 5)
    dummy_c = np.array([v*5 for v in range(args.plot_iter + 1)])
    sc_cb = ax.scatter(dummy_x, dummy_y, c=dummy_c, vmin=0, vmax= 5 *(args.plot_iter - 1),
                       s=10.0, cmap=cmap, label=label)
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)

    if args.legend:
        ax.legend(loc='upper center', ncol=2, fontsize='medium')
    if args.colorbar:
        cb_ax = fig.add_axes([0.13, 0.15, 0.02, 0.70], anchor='C')
        fig.colorbar(sc_cb, cax=cb_ax)

    save_nm = os.path.join(out_dir, 'rrt_evolution_{}.png'.format(nm))
    print("Saving to {}...".format(save_nm))
    fig.savefig(save_nm, bbox_inches='tight')

reason: import numpy


text:
from flask import Flask, render_template, request
from PIL import Image
import numpy as np
import re
import os
import base64
import requests
import subprocess
from PIL import Image
import PIL.ImageOps

# Declare a flask app
app = Flask(__name__)

FILE_NAME = 'output.png'

def parseImage(imgData):
    # parse canvas bytes for the image
    imgstr = re.search(b'base64,(.*)', imgData).group(1)

    # write the canvas image
    with open(FILE_NAME,'wb') as output:
        output.write(base64.b64decode(imgstr))

    # open the output canvas image, invert the colors and save it
    image = Image.open(FILE_NAME)
    inverted_image = PIL.ImageOps.invert(image.convert('RGB'))
    inverted_image.save(FILE_NAME)

@app.route('/')
def index():
    # Main page
    return render_template('index.html')

@app.route('/predict', methods=['GET', 'POST'])
def predict():
    # get data from drawing canvas and save as image
    parseImage(request.get_data())

    # call haskell script and pass in the model structure file, trained model file, and the filepath to the canvas image
    res = subprocess.check_output(['stack', 'build', '--exec', 'image-exe run sigmoid.cfg trainedSigmoid.cfg ' + FILE_NAME])

    # return the prediction of the canvas image
    return (res)

if __name__ == '__main__':
    app.run(host='0.0.0.0', port=9091, debug=True)


text:
#!/usr/bin/env python

from __future__ import division

import matplotlib as mpl
mpl.use('Agg')

import matplotlib.pyplot as plt
plt.style.use('seaborn-paper')
#plt.style.use('fivethirtyeight')
import numpy as np

mpl.rcParams.update({'text.usetex': True,
                     'lines.linewidth': 2.5,
                     'font.size': 16,
                     'xtick.labelsize': 'small',
                     'ytick.labelsize': 'small',
                     'axes.grid': True,
                     'axes.labelsize': 'medium',
                     'grid.alpha': 0.73,
                     'lines.markersize': 12,
                     'legend.borderpad': 0.2,
                     'legend.fancybox': True,
                     'legend.fontsize': 13,
                     'legend.framealpha': 0.7,
                     'legend.handletextpad': 0.1,
                     'legend.labelspacing': 0.2,
                     'legend.loc': 'best',
                     'savefig.dpi': 100,
                     'pdf.compression': 9})

x = np.arange(0.0, 1.0, 0.001)

A    = 1
phi0 = 0
f    = 9
omega = 2*np.pi*f
tau  = 0.5
L    = 1
#s    = A * np.sin(2*np.pi*f*t + phi0)
#dy   = np.exp(-t/tau)
#s   *= dy

fig = plt.figure(11122, figsize=(8,4))
ax = fig.add_subplot(111)

n = 4
t = 0
s = A * np.sin(n*np.pi * x / L - omega*t)
ax.plot(x,  s, alpha=0.85,
            rasterized = True,
            label = "t = 0")

t = 0.01
s = A * np.sin(n*np.pi * x / L - omega*t)
ax.plot(x,  s, alpha=0.7,
            rasterized = True,
            label = "t = " + str(t))

ax.annotate(
     r'$y(x,t) = y(x + \Delta x, t + \frac{\Delta x}{v})$',
     fontsize='x-small',
     xy=(-0.031, -0.97), xycoords='data',
     xytext=(0.5, -0.6), textcoords='offset points')

ax.set_xlabel('Displacement [L]')
ax.set_ylabel('Amplitude [mm]')
ax.set_title(r'$y = sin(\frac{2 \pi}{\lambda} (x - v t))$', fontsize=19)
ax.grid(True)
#ax.yaxis.set_ticks(np.arange(-1, 1.1, 0.5))
ax.legend()
plt.savefig("detail-trav.pdf", bbox_inches='tight')


text:
# dependencies
import numpy as np
import skimage.io

def contrast(input_img, intensity=5, display=False, output_img=''):

    '''
    Adjust the contrast of an image, given the intensity specificication.
    Works by making dark pixels much darker, and lighter pixels only slightly darker.

    Input
    ============
    input_img -- string of path for an image file.
    intensity -- integer in range [0,10] to specificy intensity of contrast enhancement applied to the image,
                 defaults to 5.
    display -- bool. If True, output will be displayed. Defaults to False.
    output_img -- string path for output image with adjusted contrast.

    Returns
    ============
    an image file with adjusted contrast in specified output path.

    Example
    ============
    contrast("../img/image.png", 5, display=False, "../img/contrast_image.png")
    '''

    # load input image
    # lift pixel value restrictions between 0 and 255 temporarily
    try:
        img = skimage.io.imread(input_img).astype(int)

    except AttributeError:
        print("Please use a string to specify the file path for an input image.")
        raise

    except FileNotFoundError:
        print("Cannot find image file.")
        raise

    except OSError:
        print("Please provide a valid image file.")
        raise

    # validate intensity level
    if (intensity < 0 or intensity > 10):
        raise ValueError("Intensity level must be between 0 and 10.")

    # map user specified intensity to a desired level of contrast (C)
    # note that intensity is between the range 0 and 10
    # and that contrast level (C) is between the range 0 and 128
    # attribution: https://stackoverflow.com/questions/345187/math-mapping-numbers
    C = intensity / 10 * 128

    # obtain contrast correction factor (F) from the desired contrast level (C)
    # attribution: https://bit.ly/2EObG6B
    F = (259 * (C + 255)) / (255 * (259 - C))

    # perform contrast enhancement to image
    enhanced_img = F * (img - 128) + 128

    # restrict values to being an integer between 0 and 255 for each pixel
    def restrict(x):
        x = round(x, 0)
        if x > 255:
            return 255
        elif x < 0:
            return 0
        else:
            return x

    np_restrict = np.vectorize(restrict)
    enhanced_img = np_restrict(enhanced_img).astype(int)

    # display or save enhanced image
    if (display == True):
        skimage.io.imshow(enhanced_img)

    if (output_img != ''):
        try:
            skimage.io.imsave(output_img, enhanced_img, check_contrast=False)

        except FileNotFoundError:
            print("The file path for the output image does not exist.")
            raise

        except ValueError:
            print("Please specify a valid output image type.")
            raise

max_stars_repo_name: mkeyim/picfixPy

max_issues_repo_name: mkeyim/picfixPy

max_forks_repo_name: mkeyim/picfixPy


text:
#  Copyright 2022 Feedzai
#
#  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.

import pandas as pd
import numpy as np
import copy
import altair as alt

def plot_feat_barplot(feat_data: pd.DataFrame,
                      top_x_feats: int = 15,
                      plot_features: dict = None
                      ):
    """Plots local feature explanations

    Parameters
    ----------
    feat_data: pd.DataFrame
        Feature explanations

    top_x_feats: int
        The number of feature to display.

    plot_features: dict
        Dict containing mapping between model features and display features
    """
    feat_data = copy.deepcopy(feat_data)
    if plot_features:
        plot_features['Pruned Events'] = 'Pruned Events'
        feat_data['Feature'] = feat_data['Feature'].apply(lambda x: plot_features[x])

    feat_data['sort_col'] = feat_data['Shapley Value'].apply(lambda x: abs(x))

    if top_x_feats is not None and feat_data.shape[0] > top_x_feats:
        sorted_df = feat_data.sort_values('sort_col', ascending=False)
        cutoff_contribution = abs(sorted_df.iloc[4]['Shapley Value'])
        feat_data = feat_data[np.logical_or(feat_data['Explanation'] >= cutoff_contribution, feat_data['Explanation'] <= -cutoff_contribution)]

    a = alt.Chart(feat_data).mark_bar(size=15, thickness=1).encode(
        y=alt.Y("Feature", axis=alt.Axis(title="Feature", labelFontSize=15,
                                         titleFontSize=15, titleX=-61),
                sort=alt.SortField(field='sort_col', order='descending')),
        x=alt.X('Shapley Value', axis=alt.Axis(grid=True, title="Shapley Value",
                                            labelFontSize=15, titleFontSize=15),
                scale=alt.Scale(domain=[-0.1, 0.4])),
    )

    line = alt.Chart(pd.DataFrame({'x': [0]})).mark_rule(
        color='#798184').encode(x='x')

    feature_plot = (a + line).properties(
        width=190,
        height=225
    )
    return feature_plot

def plot_global_feat(feat_data: pd.DataFrame,
                     top_x_feats=12,
                     threshold=None,
                     plot_features=None,
                     **kwargs
                     ):
    """ Plots global feature plots

    Parameters
    ----------
    feat_data: pd.DataFrame
        Feature explanations to plot

    top_x_feats: int
        The number of feature to display.

    threshold: int
        The minimum absolute importance that a feature needs to have to be displayed

    plot_features: dict
        Dict containing mapping between model features and display features
    """
    avg_df = feat_data.groupby('Feature').mean()['Shapley Value']
    if threshold is None and len(avg_df) >= top_x_feats:
        sorted_series = avg_df.abs().sort_values(ascending=False)
        threshold = sorted_series.iloc[top_x_feats-1]
    if threshold:
        avg_df = avg_df[np.logical_or(avg_df <= -threshold, avg_df >= threshold)]
    feat_data = feat_data[feat_data['Feature'].isin(avg_df.index)][['Shapley Value', 'Feature']]

    if threshold:
        avg_df = avg_df.append(pd.Series([0], index=['(...)']))
        #some prunning happened
        feat_data = feat_data.append({'Feature': '(...)', 'Shapley Value': -0.6, }, ignore_index=True)
    feat_data['type'] = 'Shapley Value'

    for index, value in avg_df.items():
        if index == '(...)':
            feat_data = feat_data.append(
                {'Feature': index, 'Shapley Value': None, 'type': 'Mean'},
                ignore_index=True)
        else:
            feat_data = feat_data.append(
                {'Feature': index, 'Shapley Value': value, 'type': 'Mean'},
                ignore_index=True)

    sort_features = list(avg_df.sort_values(ascending=False).index)
    if plot_features:
        plot_features = copy.deepcopy(plot_features)
        plot_features['Pruned Events'] = 'Pruned Events'
        plot_features['(...)'] = '(...)'
        feat_data['Feature'] = feat_data['Feature'].apply(lambda x: plot_features[x])
        sort_features = [plot_features[x] for x in sort_features]

    global_feats_plot = alt.Chart(feat_data).mark_point(stroke='white',
                                                         strokeWidth=.6).encode(
        x=alt.X('Shapley Value', axis=alt.Axis(title='Shapley Value', grid=True),
                scale=alt.Scale(domain=[-0.2, 0.6])),
        y=alt.Y('Feature:O',
                sort=sort_features,
                axis=alt.Axis(labelFontSize=13, titleX=-51)),
        color=alt.Color('type',
                        scale=alt.Scale(domain=['Shapley Value', 'Mean'],
                                        range=["#618FE0", '#d76d58']),
                        legend=alt.Legend(title=None, fillColor="white",
                                          symbolStrokeWidth=0, symbolSize=50,
                                          orient="bottom-right")),
        opacity=alt.condition(alt.datum.type == 'Mean', alt.value(1.0),
                              alt.value(0.1)),
        size=alt.condition(alt.datum.type == 'Mean', alt.value(70),
                           alt.value(30)),
    ).properties(
        width=288,
        height=280
    )
    return global_feats_plot

max_stars_repo_name: feedzai/timeshap

max_issues_repo_name: feedzai/timeshap

max_forks_repo_name: feedzai/timeshap


text:
import numpy as np
from lights.model.associations import AssociationFunctionFeatures
from sklearn.preprocessing import StandardScaler
import numba as nb
from llvmlite import binding
binding.set_option('SVML', '-vector-library=SVML')
from lights.simulation import features_normal_cov_toeplitz

class EstepFunctions:
    """A class to define functions relative to the E-step of the prox_QNMCEM

    Parameters
    ----------
    T_u : `np.ndarray`, shape=(J,)
        The J unique training censored times of the event of interest

    n_long_features :  `int`,
            Number of longitudinal features

    fixed_effect_time_order : `int`
        Order of the higher time monomial considered for the representations
        of the time-varying features corresponding to the fixed effect. The
        dimension of the corresponding design matrix is then equal to
        fixed_effect_time_order + 1

    asso_functions_list : `list` or `str`='all'
        List of association functions wanted or string 'all' to select all
        defined association functions. The available functions are :
            - 'lp' : linear predictor
            - 're' : random effects
            - 'tps' : time dependent slope
            - 'ce' : cumulative effects

    theta : `dict`
        Vector that concatenates all parameters to be inferred in the lights
        model
    """

    def __init__(self, T_u, n_long_features, fixed_effect_time_order,
                 asso_functions_list, theta):
        self.theta = theta
        self.n_long_features = n_long_features
        self.J = len(T_u)
        self.nb_total_asso_features = len(theta["gamma_0"])
        self.fixed_effect_time_order = fixed_effect_time_order
        self.q_l, r_l = fixed_effect_time_order + 1, 2
        alpha, L = self.fixed_effect_time_order, self.n_long_features
        self.F_f, self.F_r = AssociationFunctionFeatures(asso_functions_list,
                                                T_u, alpha, L).get_asso_feat()
        self.asso_funcs = None

    def compute_AssociationFunctions(self, S, simu, cov_corr_rdn_long=None, S_k=None):
        """
        Compute the value of association functions

        Parameters
        ----------
        S : `np.ndarray`, shape=(N_MC, r)
            Set of constructed Monte Carlo samples, with N_MC = 2 * N

        simu : `bool`, defaut=True
            If `True` we comute the asso feat with simulated data.

        cov_corr_rdn_long : `float`
            Correlation coefficient of the toeplitz correlation matrix of
            random longitudinal features.

        S_k : `list`
            Set of nonactive group for 2 classes (will be useful in case of
            simulated data)

        """
        N_MC = S.shape[0]
        K = 2
        L = self.n_long_features
        J, nb_total_asso_features = self.J, self.nb_total_asso_features
        q_l, r_l = self.q_l, 2
        nb_total_asso_param = nb_total_asso_features // L
        nb_asso_param = self.F_f.shape[1]
        nb_noise_param = nb_total_asso_param - nb_asso_param
        self.asso_funcs = np.zeros((J, N_MC, K, nb_total_asso_features))
        beta = np.hstack((self.theta["beta_0"], self.theta["beta_1"])).T
        if simu:
            # Correlation coefficient of the toeplitz correlation matrix
            rho = cov_corr_rdn_long
            for k in range(K):
                for l in range(L):
                    start_idx = nb_total_asso_param * l
                    stop_idx = nb_total_asso_param * l + nb_asso_param
                    if l not in S_k[k]:
                        beta_tmp = beta[k, r_l * l: r_l * (l + 1)]
                        S_tmp = S[:, q_l * l: q_l * (l + 1)].T
                        tmp = (self.F_f.dot(beta_tmp).T +
                               self.F_r.dot(S_tmp).T).T.swapaxes(1, -1)
                        self.asso_funcs[:, :, k, start_idx: stop_idx] = tmp
                        self.asso_funcs[:, :, k, nb_asso_param + nb_total_asso_param * l
                        : nb_total_asso_param * (l + 1)] = features_normal_cov_toeplitz(
                            J * N_MC,nb_noise_param, rho)[0].reshape(J, N_MC, -1)

                    else:
                        self.asso_funcs[:, :, k, nb_total_asso_param * l : nb_total_asso_param * (l + 1)] = \
                        features_normal_cov_toeplitz(J * N_MC, nb_total_asso_param,
                                rho, .1)[0].reshape(J, N_MC, -1)
                # normalize features
                shape = self.asso_funcs[:, :, k].shape
                tmp_reshaped = self.asso_funcs[:, :, k].copy().reshape(-1, shape[-1])
                self.asso_funcs[:, :, k] = StandardScaler().fit_transform(
                    tmp_reshaped).copy().reshape(shape)
        else:
            for k in range(K):
                for l in range(L):
                    start_idx = nb_total_asso_param * l
                    stop_idx = nb_total_asso_param * l + nb_asso_param
                    beta_tmp = beta[k, r_l * l: r_l * (l + 1)]
                    S_tmp = S[:, q_l * l: q_l * (l + 1)].T
                    tmp = (self.F_f.dot(beta_tmp).T +
                            self.F_r.dot(S_tmp).T).T.swapaxes(1, -1)
                    self.asso_funcs[:, :, k, start_idx: stop_idx] = tmp
            pass

    def construct_MC_samples(self, N_MC):
        """Constructs the set of samples used for Monte Carlo approximation

        Parameters
        ----------
        N_MC : `int`
            Number of Monte Carlo samples

        Returns
        -------
        S : `np.ndarray`, shape=(N_MC, r)
            Set of constructed Monte Carlo samples, with N_MC = 2 * N
        """
        D = self.theta["long_cov"]
        C = np.linalg.cholesky(D)
        r = D.shape[0]
        Omega = np.random.multivariate_normal(np.zeros(r), np.eye(r), N_MC)
        b = Omega.dot(C.T)
        S = np.vstack((b, -b))

        return S

    def g1(self, S):
        """Computes g1

        Parameters
        ----------
        S : `np.ndarray`, shape=(N_MC, r)
            Set of constructed Monte Carlo samples

        Returns
        -------
        g1 : `np.ndarray`, shape=(N_MC, r)
            The values of g1 function
        """

        g1 = S
        return g1

    def g2(self, S):
        """Computes g2

        Parameters
        ----------
        S : `np.ndarray`, shape=(N_MC, r)
            Set of constructed Monte Carlo samples

        Returns
        -------
        g2 : `np.ndarray`, shape=(N_MC, r, r)
            The values of g2 function
        """
        g2 = np.array([s.reshape(-1, 1).dot(s.reshape(-1, 1).T) for s in S])
        return g2

    def g3(self):
        """Computes g3
        Returns
        -------
        g3 : `np.ndarray`, shape=(N_MC, J, dim, K)
            The values of g3 function
        """
        g3 = self.asso_funcs.swapaxes(0, 1).swapaxes(2, 3)
        return g3

    def g4(self, gamma_0, gamma_1):
        """Computes g4

        Parameters
        ----------
        gamma_0 : `np.ndarray`, shape=(L * nb_asso_param,)
            Association parameters for low-risk group

        gamma_1 : `np.ndarray`, shape=(L * nb_asso_param,)
            Association parameters for high-risk group

        Returns
        -------
        g4 : `np.ndarray`, shape=(N_MC, J, K)
            The values of g4 function
        """
        gamma = np.hstack((gamma_0, gamma_1)).T
        g4 = np.exp((self.asso_funcs * gamma).sum(axis=-1).swapaxes(0, 1))
        return g4

    def g5(self, gamma_0, gamma_1):
        """Computes g5
        Parameters
        ----------
        gamma_0 : `np.ndarray`, shape=(L * nb_asso_param,)
            Association parameters for low-risk group

        gamma_1 : `np.ndarray`, shape=(L * nb_asso_param,)
            Association parameters for high-risk group
        Returns
        -------
        g5 : `np.ndarray`, shape=(N_MC, J, dim, K)
            The values of g5 function
        """
        g3 = self.g3()
        g4 = self.g4(gamma_0, gamma_1)
        g5 = g4[:, :, np.newaxis, :] * g3
        return g5

    @staticmethod
    def Lambda_g(g, f):
        """Approximated integral (see (15) in the lights paper)

        Parameters
        ----------
        g : `np.ndarray`, shape=(N_MC, ...)
            Values of g function for all subjects, all groups and all Monte
            Carlo samples. Each element could be real or matrices depending on
            Im(\tilde{g}_i)

        f: `np.ndarray`, shape=(n_samples, K, N_MC)
            Values of the density of the observed data given the latent ones and
            the current estimate of the parameters, computed for all subjects,
            all groups and all Monte Carlo samples

        Returns
        -------
        Lambda_g : `np.array`, shape=(n_samples, K, shape(g))
            The approximated integral computed for all subjects, all groups and
            all Monte Carlo samples. Each element could be real or matrices
            depending on Im(\tilde{g}_i)
        """
        if len(g.shape) == 1:
            Lambda_g = Lambda_g_1_nb(g, f)
        elif len(g.shape) == 2:
            Lambda_g = Lambda_g_2_nb(g, f)
        elif len(g.shape) == 3:
            Lambda_g = Lambda_g_3_nb(g, f)
        elif len(g.shape) == 4:
            Lambda_g = Lambda_g_4_nb(g, f)
        else:
            raise ValueError("The shape of g function is not supported")
        return Lambda_g

    def Eg(self, g, Lambda_1, pi_xi, f):
        """Computes approximated expectations of different functions g taking
        random effects as input, conditional on the observed data and the
        current estimate of the parameters. See (13) in the lights paper

        Parameters
        ----------
        g : `np.array`
            The value of g function for all samples

        Lambda_1: `np.ndarray`, shape=(n_samples, K)
            Approximated integral (see (17) in the lights paper) with
            \tilde(g)=1

        pi_xi: `np.ndarray`, shape=(n_samples,)
            The posterior probability of the sample for being on the
            high-risk group given all observed data

        f: `np.ndarray`, shape=(n_samples, K, N_MC)
            The value of the f(Y, T, delta| S, G ; theta)

        Returns
        -------
        Eg : `np.ndarray`, shape=(n_samples, g.shape)
            The approximated expectations for g
        """
        Lambda_g = self.Lambda_g(g, f)
        Eg = (Lambda_g[:, 0].T * (1 - pi_xi) + Lambda_g[:, 1].T * pi_xi) / (
                Lambda_1[:, 0] * (1 - pi_xi) + Lambda_1[:, 1] * pi_xi)
        return Eg.T

@nb.njit(parallel=True, fastmath=True)
def Lambda_g_1_nb(g, f):
    n_samples, K, N_MC = f.shape
    res = np.zeros((n_samples, K))
    for i in nb.prange(n_samples):
        for k in nb.prange(K):
            for j in nb.prange(N_MC):
                res[i, k] += (g[j] * f[i, k, j]) / N_MC
    return res

@nb.njit(parallel=True, fastmath=True)
def Lambda_g_2_nb(g, f):
    n_samples, K, N_MC,  = f.shape
    res = np.zeros((n_samples, K,) + g.shape[1:])
    for i in nb.prange(n_samples):
        for k in nb.prange(K):
            for j in nb.prange(N_MC):
                for m in nb.prange(g.shape[1]):
                    res[i, k, m] += (g[j, m] * f[i, k, j]) / N_MC
    return res

@nb.njit(parallel=True, fastmath=True)
def Lambda_g_3_nb(g, f):
    n_samples, K, N_MC,  = f.shape
    res = np.zeros((n_samples, K,) + g.shape[1:])
    for i in nb.prange(n_samples):
        for k in nb.prange(K):
            for j in nb.prange(N_MC):
                for m in nb.prange(g.shape[1]):
                    for n in nb.prange(g.shape[2]):
                        res[i, k, m, n] += \
                                        (g[j, m, n] * f[i, k, j]) / N_MC
    return res

@nb.njit(parallel=True, fastmath=True)
def Lambda_g_4_nb(g, f):
    n_samples, K, N_MC,  = f.shape
    res = np.zeros((n_samples, K,) + g.shape[1:])
    for i in nb.prange(n_samples):
        for k in nb.prange(K):
            for j in nb.prange(N_MC):
                for m in nb.prange(g.shape[1]):
                    for n in nb.prange(g.shape[2]):
                        for t in nb.prange(g.shape[3]):
                            res[i, k, m, n, t] += \
                                (g[j, m, n, t] * f[i, k, j]) / N_MC
    return res

max_stars_repo_name: Califrais/lights

max_issues_repo_name: Califrais/lights

max_forks_repo_name: Califrais/lights

reason: import numpy,import numba


text:
"""
    Jacob Kaplan
    functions.py

    Implementations of functions called in align.py
"""

import os
import random
import cv2 as cv
import numpy as np

def findImages(directory):
    """
    Take in directory of images, open directory, read in each image,
    add to list of images, return list of images
    """
    images = list()
    for fil in os.listdir(directory):
        if fil.lower().endswith('.jpg'):
            try:
                img = cv.imread(directory + "/" + fil, 1)
            except cv.error:
                print("{} malformed!".format(fil))
                sys.exit()
            images.append((img, fil))
    return images

"""
1. Extract and match keypoints =================================================
"""
def matchKeyPoints(img1, img2):
    """
    Take in two images
    Get a set of keypoints and descriptors for each image
    Matches both sets of keypoints with brute-force matcher
    Perform ratio test to get best matches
    Extract the keypoints that survived ratio test and return them
    """
    kps1, dsc1 = getSIFT(img1)
    kps2, dsc2 = getSIFT(img2)
    allMatches = matchAll(dsc1, dsc2)
    ratioMatches = ratioTest(kps1, kps2, allMatches)
    kps1, kps2 = extractKps(kps1, kps2, ratioMatches)
    return kps1, kps2

def getSIFT(img):
    """
    Take in image
    Return SIFT keyPoints and descriptors
    """
    sift = cv.xfeatures2d.SIFT_create()
    keyPoints, descriptors = sift.detectAndCompute(img, None)
    return keyPoints, descriptors

def matchAll(dsc1, dsc2):
    """
    Take in two sets of descriptors for keypoints
    Use brute-force matcher to get matches
    Return matches
    """
    bfMatcher = cv.BFMatcher()
    matches = bfMatcher.knnMatch(dsc1,dsc2,k=2)
    return matches

def ratioTest(kps1, kps2, matches):
    """
    Take in two sets of keypoints and the set of matches between them
    Perform ratio test to get best matches (60% threshold)
    Return surviving matches
    """
    ratioMatches = list()
    for match in matches:
        m1, m2 = match
        if m1.distance < 0.60*m2.distance:
            ratioMatches.append(m1)
    return ratioMatches

def extractKps(kps1, kps2, ratioMatches):
    """
    Take in two sets of keypoints and the matches that survived ratio test
    Extract all the keypoints in the set of ratio matches
    Return extracted keypoints
    """
    newKps1 = list()
    newKps2 = list()
    for match in ratioMatches:
        newKps1.append(kps1[match.queryIdx].pt)
        newKps2.append(kps2[match.trainIdx].pt)
    return newKps1, newKps2

def outputMatchInfo(kps, name1, name2):
    """
    Take in keypoints and image names and output info about them
    """
    print("\nComparing {} and {}...".format(name1, name2))
    print("Applying ratio test...")
    print("{} matches found...".format(len(kps)))
    print("Drawing best matches...")

"""
2. Fundamental matrix ==========================================================
"""
def checkMatches(kps):
    """
    Take in keypoints
    If they are above 100, there are enough to continue, return True
    Else, False
    """
    if len(kps) >= 100:
        print("There are enough matches to continue...")
        return True
    else:
        print("There are not enough matches to continue.")
        return False

def getFundamental(kps1, kps2):
    """
    Take in two set of keypoints
    Convert to np arrays and get fundamental matrix using RANSAC algorithm
    Use mask of fundamental matrix to get keypoints that survived the estimation
    Return the inlying keypoints
    """
    kps1 = np.int32(kps1)
    kps2 = np.int32(kps2)
    funMat, mask = cv.findFundamentalMat(kps1, kps2, cv.RANSAC)
    ind = np.where(mask.ravel() == 1)[0]
    inl1 = kps1[ind]
    inl2 = kps2[ind]
    return inl1, inl2

def outputFundamentalInfo(kps, inl):
    """
    Take in keypoints before and after fundamental matrix estimation
    and output info about how many survived
    """
    print("Building fundamental matrix...")
    print("Of the {} best matches...".format(len(kps)))
    print("{} survived the fundamental matrix estimation...".format(len(inl)))
    print("Drawing inliers...")

"""
3. Check if there are enough inliers to continue ===============================
"""
def checkInliers(kps, inl):
    """
    Take in keypoints before and after fundamental matrix estimation
    If 80% or more survived, continue with program, return True
    Else, False
    """
    percentRemain = round(100* len(inl) / len(kps))
    print("{}% matches remain...".format(percentRemain))
    if percentRemain >= 80:
        print("That means the images are likely of the same scene...")
        return True
    else:
        print("The images are probably not of the same scene.")
        return False

"""
4. Estimate homography matrix ==================================================
"""
def getHomography(inl1, inl2):
    """
    Take in inliers of fundamental matrix estimation
    Use them to estimate homography matrix using RANSAC algorithm
    Use mask to get keypoints that survived homography estimation
    Return inliers and homography matrix
    """
    hmMat, mask = cv.findHomography(inl1, inl2, cv.RANSAC, 5.0)
    ind = np.where(mask.ravel() == 1)[0]
    hm1 = inl1[ind]
    hm2 = inl2[ind]
    return hm1, hm2, hmMat

def outputHomographyInfo(inl, hm):
    """
    Take in keypoints before and after homography matrix estimation
    Output info about how many remain
    """
    print("Building homography matrix ...")
    print("Of the {} original inliers...".format((len(inl))))
    print("{} survived the homography matrix estimation...".format(len(hm)))
    print("Drawing inliers...")

"""
5. Check alignment =============================================================
"""
def checkAlignment(inl, hm):
    """
    Take in keypoints before and after homography matrix estimation
    Get percent difference and if it's less than 30, return True
    Else, False
    """
    percentDiff = np.abs(round(100*(len(inl) - len(hm)) / len(inl)))
    print("Homography and fundamental matrix had a {}% difference in results..."\
            .format(percentDiff))
    if percentDiff <= 30:
        print("The images can be aligned...")
        return True
    else:
        print("The images cannot be aligned.")
        return False

"""
6. Build mosaic ================================================================
"""
def buildMosaic(img1, img2, homography):
    """
    Python Gods forgive me for this behemoth...

    Take in two images and their homography matrix
    Using the corners of the images, get the dimensions of the shape of the
      mosaic shape based on transform around the homography matrix
    Use the dimensions of the mosaic shape to get the coordinates of where to
      cut the original images (stitchX, stitchY)
    Use homography inverse to map image 2 to the new mosaic shape
    Get the overlap of the mapped image 2 and the original image 1
    Average the color of the overlapping section
    Paste image 1 onto the mapped image 2
    Return mosaic
    """

    print("Building mosaic...")
    # Original dimensions
    m1, n1 = img1.shape[:2]
    m2, n2 = img2.shape[:2]
    corners1 = np.asarray([[0,0], [0,m1], [n1,m1], [n1,0]]).reshape(-1,1,2)
    corners2 = np.asarray([[0,0], [0,m2], [n2,m2], [n2,0]]).reshape(-1,1,2)

    # Mosaic shape and dimensions
    transform = cv.perspectiveTransform(np.float32(corners2), homography)
    mosaicShape = np.vstack((corners1, transform))
    mCornersBot = np.min(mosaicShape, axis=0)
    mCornersTop = np.max(mosaicShape, axis=0)
    x1, y1 = mCornersTop[0,:2]
    x2, y2 = mCornersBot[0,:2]
    stitchX, stitchY = int(round(-1*x2)), int(round(-1*y2))

    # Map image 2 onto mosaic using invere homography
    affine = np.array([[1, 0, stitchX], [0, 1, stitchY], [0,0,1]])
    x1, x2 = int(round(x1)), int(round(x2))
    y1, y2 = int(round(y1)), int(round(y2))
    mosaic = cv.warpPerspective(img2, affine.dot(homography), (x1-x2, y1-y2))

    # Get overlapping section of mapped image 2 and original image 1
    mosaicZeros = np.zeros_like((mosaic))
    ms1 = mosaic.copy()
    ms2 = mosaicZeros.copy()
    ms1[ms1 > 0] = 1
    ms2[stitchY:stitchY+m1, stitchX:stitchX+n1] = 1
    overlap = ms1 + ms2
    mosaicTmp = mosaic.copy()
    mosaicZeros[stitchY:stitchY + m1, stitchX:stitchX + n1] = img1
    mosaic[stitchY:stitchY + m1, stitchX:stitchX + n1] = img1
    # Average the colors
    mosaic[np.where(overlap == 2)] = 0.5*mosaicZeros[np.where(overlap == 2)] \
                                    + 0.5*mosaicTmp[np.where(overlap == 2)]
    print("Mosaic built!")
    return mosaic

"""
   Drawing lines between keypoints =============================================
"""
def drawMatches(img1, img2, kps1, kps2):
    """
    Take in two images and their corresponding sets of keypoints
    Paste the two images onto a larger canvas
    Iterate through both sets of keypoints and extract coordinates
    Draw circle of random color around each keypoint onto new image
    Draw line of random color through each keypoint pair onto new image
    Return new image
    """
    m1, n1 = img1.shape[:2]
    m2, n2 = img2.shape[:2]
    draw = np.zeros((max(m1, m2), n1 + n2, 3), dtype="uint8")
    draw[0:m1, 0:n1] = img1
    draw[0:m2, n2:] = img2

    for i in range(len(kps1)):
        b = random.randint(0,255)
        g = random.randint(0,255)
        r = random.randint(0,255)

        pt1 = (int(kps1[i][0]), int(kps1[i][1]))
        pt2 = (int(kps2[i][0]) + n1, int(kps2[i][1]))

        cv.circle(draw, pt1, 1, (b,g,r), 1)
        cv.circle(draw, pt2, 1, (b,g,r), 1)
        cv.line(draw, pt1, pt2, (b, g, r), 2)

    return draw


text:
import pandas as pd
import numpy as np
import pickle
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.colors import Normalize
from matplotlib.collections import PatchCollection
from mpl_toolkits.basemap import Basemap
from shapely.geometry import Point, Polygon, MultiPoint, MultiPolygon
from shapely.prepared import prep
from pysal.esda.mapclassify import Natural_Breaks as nb
from descartes import PolygonPatch
from itertools import chain
from scipy.interpolate import Rbf
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from osgeo import gdal
import shapefile

# unpickle the dataset:
#df[0] = ID, [1] = name, [2] = elevation, [3] = latitude, [4] = longitude,
#[5] = Monthly averages: [5][station][year][month],
#[6] = Average Maximum Temperature: [6][station][month],
#[7] = Standard Deviation Maximum Temperature: [7][station][month],
#[8] = monthly maximum temperature anomalies: [8][station][year][month],
#[9] = Maximum t-score: [9][station][year][month],
#[10] Monthly minima (see [5]),
#[11] = Average Minimum Temperature (see[6]),
#[12] = Standard Deviation Minimum Temperature (see[7]),
#[13] monthly minimum temperature anomalies (see[8]),
#[14] Minimum t-score (see[9])

df = pd.read_pickle('biglist.pkl')
#there are fifty states and they are abbreviated like this:
states = ['MA', 'VT', 'NH', 'CT', 'RI','CA','ME', 'OR','WA','UT','NM','AZ','ID','CO','TX','WY','MT','ND','SD','NV','NE','KS','OK','LA','AR','MO','IA','MN','MS','TN','IL','KY','IN','WI','MI','OH','WV','AL','FL','GA','SC','NC','VA','PA','DC','MD','DE','NJ', 'NY','IA']

# generate the basemap and load in the shapefiles
m= Basemap(llcrnrlon=-125., llcrnrlat = 24., urcrnrlon = -66., urcrnrlat = 50., resolution=None, projection='merc')
sf = shapefile.Reader("states")

#set figure
fig = plt.figure()
ax = fig.add_subplot(111)

#generate shapefile vertices lists
verticesx = []
verticesy = []
codes = []

#Import per basemap tutorial: http://basemaptutorial.readthedocs.org/en/latest/clip.html
for shape_rec in sf.shapeRecords():
    if shape_rec.record[4] in states:
        pts = shape_rec.shape.points
        prt = list(shape_rec.shape.parts) + [len(pts)]
        for i in range(len(prt) - 1):
            for j in range(prt[i], prt[i+1]):
                verticesx.append(pts[j][0])
                verticesy.append(pts[j][1])
            codes += [Path.MOVETO]
            codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
            codes += [Path.CLOSEPOLY]
        vertice = m(verticesx, verticesy)
        vertices = [[j[i] for j in vertice] for i in range(len(vertice[0]))]
        clip = Path(vertices, codes)
        clip = PathPatch(clip, transform = ax.transData)

#Extract Variable of interest (in this case the maximum temperature t-score - for later: loop and automate, save all graphics)
Zs = [ df[9][x][10][1] for x in range(0, len(df[9]))]

#Because RBF interpolation does not like nan, extract indices of non-nans
inds = pd.notnull(Zs).nonzero()

#collect other parameters of interest
Longs = df[4]
Lats = df[3]

#pare down lists to only non-nan values
zs = [Zs[l] for l in inds[0]]
Longs = [Longs[l] for l in inds[0]]
Lats = [Lats[l] for l in inds[0]]

#prepare longitude and latitude for RBF and mapping to meters coordinates, grid.
lons = np.array([l for l in Longs])
lats = np.array([l for l in Lats])
x , y = m(lons, lats)
longti = np.linspace(min(x)-10000, max(x)+10000, 1000)
latti = np.linspace(min(y)- 30000, max(y)+30000, 543)
XI, YI = np.meshgrid(longti, latti)

#interpolate
rbf = Rbf(x,y,zs,function='thin_plate')
ZI = rbf(XI, YI)

#set levels if desired
levelsz = np.linspace(-5,5,11)
#levelsstd = np.linspace(0,5,6)
#levelsresults= np.linspace(260, 302, 21)
#levelsavg = np.linspace(260,302,21)
#levelsanom=np.linspace(-12,10,23)

#Generate contours from interpolation using specified levels
cs = m.contourf(XI,YI,ZI, levels = levelsz, alpha = 0.5, zorder = 10)

m.colorbar(location='bottom', label='t-score')
plt.title('January 1961 Maximum Temperature Anomaly')
#fill in the state shapes
for contour in cs.collections:
    contour.set_clip_path(clip)

#plot the individual datapoints and draw map boundaries if desired
#m.scatter(x, y, latlon = False,  marker = 'o', c='k',  zorder= 10)
m.readshapefile('states', 'states', drawbounds = True, color = 'w', zorder=11)
#m.savefig('map.tif', dpi=dpi)

plt.show()

max_stars_repo_name: jklobas/weatherproject

max_issues_repo_name: jklobas/weatherproject

max_forks_repo_name: jklobas/weatherproject

reason: import numpy,from scipy


text:
# coding: utf-8
# Copyright (c) MoGroup at UMD.
# Distributed under the terms of the MIT License.

from __future__ import division, unicode_literals
import numpy as np
from monty.json import MSONable
from scipy import stats
from scipy.optimize import curve_fit
from pymatgen.io.vasp.outputs import Vasprun
from pymatgen.util.coord import pbc_diff
from pymatgen.core.structure import Structure
from pymatgen.core.periodic_table import Specie
import csv

__author__ = "Xingfeng He"
__email__ = "albert2@umd.edu"
__version__ = "0.1"
__date__ = "6/6/2017"
__credit__ = "Pymatgen Development Team"

class DiffusivityAnalyzer(MSONable):
    def __init__(self, structure, displacements, specie, temperature,
                 time_step, step_skip, time_intervals_number=1000,
                 spec_dict=None):
        """
        Calculate MSD from pre-processed data, and implemented linear fitting to obtain diffusivity.

        :param structure (Structure): initial structure
        :param displacements (np.array): numpy array, shape is [n_ions, n_steps, axis]
        :param specie (str): species string, can be Li or Li+, make sure structure has oxidation
                state accordingly.
        :param temperature (float): temperature of MD
        :param time_step (float): time step in MD
        :param step_skip (int): Sampling frequency of the displacements (
                time_step is multiplied by this number to get the real time
                between measurements)
        :param time_intervals_number (int): number of time intervals. Default is 1000
                means there are ~1000 time intervals.
        :param spec_dict (dict): spec dict of linear fitting. Default is
                {'lower_bound': 4.5, 'upper_bound': 0.5, 'minimum_msd_diff': 4.5}
                lower_bound is in unit of Angstrom square
                upper_bound is in unit of total time. 0.5 means upper fitting bound is 0.5*t_total
                minimum_msd_diff is in unit of Angstrom square. msd[upper_bound] - msd[lower_bound] should larger
                than minimum_msd_diff to do linear fitting.
        """
        spec_dict = spec_dict if spec_dict is not None else {'lower_bound': 4.5, 'upper_bound': 0.5,
                                                             'minimum_msd_diff': 4.5}
        if not {'lower_bound', 'upper_bound', 'minimum_msd_diff'} <= set(spec_dict.keys()):
            raise Exception("spec_dict does not have enough parameters.")
        time_step_displacements = time_step * step_skip
        # prepare
        indices = []
        framework_indices = []
        for i, site in enumerate(structure):
            if site.specie.symbol == specie:
                indices.append(i)
            else:
                framework_indices.append(i)
        if len(indices) == 0:
            raise Exception("There is no specie {} in the structure".format(specie))
        if len(framework_indices) == 0:
            dc = displacements
        else:
            framework_disp = displacements[framework_indices]
            drift = np.average(framework_disp, axis=0)[None, :, :]
            dc = displacements - drift
        df = structure.lattice.get_fractional_coords(dc)
        displacements_final_diffusion_ions = dc[indices]
        displacements_frac_final_diffusion_ions = df[indices]
        n_ions, n_steps, dim = displacements_final_diffusion_ions.shape

        # time intervals, dt
        dt_indices = np.arange(1, n_steps, max(int((n_steps - 1) / time_intervals_number), 1))
        dt = dt_indices * time_step_displacements

        # calculate msd
        # define functions, algorithm from
        # http://stackoverflow.com/questions/34222272/computing-mean-square-displacement-using-python-and-fft
        def autocorrelation_fft(x):
            N = x.shape[0]
            F = np.fft.fft(x, n=2 * N)
            PSD = F * F.conjugate()
            res = np.fft.ifft(PSD)
            res = (res[:N]).real
            n = N * np.ones(N) - np.arange(N)
            return res / n

        def one_ion_msd_fft(r, dt_indices):
            """
            r (np.array, shape is typically [n_step,3], n_step is number of steps, 3 is 3 dimentions)
            """
            # ------------ S1
            n_step, dim = r.shape
            r_square = np.square(r)
            r_square = np.append(r_square, np.zeros((1, dim)), axis=0)  # (n_step+1, 3)
            S1_component = np.zeros((dim, n_step))  # (dim, n_step)
            r_square_sum = 2 * np.sum(r_square, axis=0)  # (3)
            for i in range(n_step):
                r_square_sum = r_square_sum - r_square[i - 1, :] - r_square[n_step - i, :]
                S1_component[:, i] = r_square_sum / (n_step - i)
            S1 = np.sum(S1_component, axis=0)

            # ------------ S2
            S2_component = np.array([autocorrelation_fft(r[:, i]) for i in range(r.shape[1])])  # (dim, N)
            S2 = np.sum(S2_component, axis=0)

            # ------------ return
            return (S1 - 2 * S2)[dt_indices], (S1_component - 2 * S2_component)[:, dt_indices]

        n_dt = len(dt_indices)
        msd_by_ions = np.empty([0, n_dt])  # shape of n_ions * n_dt
        msd_component_by_ions = np.empty([3, 0, n_dt])  # shape of 3 * n_ions * n_dt
        for i in range(n_ions):
            msd_i, msd_component_i = one_ion_msd_fft(displacements_final_diffusion_ions[i, :, :], dt_indices)
            msd_by_ions = np.append(msd_by_ions,
                                    msd_i.reshape(1, n_dt),
                                    axis=0)
            msd_component_by_ions = np.append(msd_component_by_ions,
                                              msd_component_i.reshape(3, 1, n_dt),
                                              axis=1)
        msd = np.average(msd_by_ions, axis=0)
        msd_component = np.average(msd_component_by_ions, axis=1)

        # further things, 1. determine lower_index, upper_index 2. linear fitting, 3. error bar

        # one headache, how about error in different axis
        lower_bound_index = len(msd[msd < spec_dict['lower_bound']])
        upper_bound_index = int(len(msd) * spec_dict['upper_bound']) - 1
        """
        if lower_bound_index >= upper_bound_index - 2:
            raise Exception("Maximum MSD is {:.2f}. ".format(max(msd)) + \
                            "MSD array has shape of {}. ".format(msd.shape) + \
                            "Lower bound index is {}, upper bound index is {}. ".format(lower_bound_index,
                                                                                upper_bound_index) + \
                            "There is no enough data to fit. " + \
                            "Please consider extending your MD simulation or increasing the temperature.")

        if msd[upper_bound_index] - msd[lower_bound_index] < spec_dict['minimum_msd_diff']:
            raise Exception(
                "Maximum MSD is {:.2f}. ".format(max(msd)) + \
                "MSD at lower bound is {:.2f}, MSD at upper bound is {:.2f}. The MSD fitting range is too small. " \
                .format(msd[lower_bound_index], msd[upper_bound_index]) + \
                "Please consider extending your MD simulation or increasing the temperature.")
        """
        if lower_bound_index >= upper_bound_index - 2 or \
                                msd[upper_bound_index] - msd[lower_bound_index] < spec_dict['minimum_msd_diff']:
            slope = -1
            intercept = -1
            slope_components = np.zeros(dim)
        else:
            slope, intercept, _, _, _ = stats.linregress(dt[lower_bound_index:upper_bound_index + 1],
                                                         msd[lower_bound_index:upper_bound_index + 1])
            slope_components = np.zeros(dim)
            for i in range(dim):
                s, _, _, _, _ = stats.linregress(dt[lower_bound_index:upper_bound_index + 1],
                                                 msd_component[i, :][lower_bound_index:upper_bound_index + 1])
                slope_components[i] = s

        self.structure = structure
        self.indices = indices
        self.framework_indices = framework_indices
        self.drift = drift
        self.drift_maximum = np.max(np.abs(drift), axis=1)[0] # the maximum drift vector of the framework ions, shape is (3,)
        self.disp = displacements
        self.displacements_final_diffusion_ions = displacements_final_diffusion_ions
        self.specie = specie
        self.temperature = temperature
        self.time_step = time_step
        self.step_skip = step_skip
        self.time_step_displacements = time_step_displacements
        self.time_intervals_number = time_intervals_number
        self.spec_dict = spec_dict
        if len(framework_indices) == 0:
            self.max_framework_displacement = 0.0
        else:
            self.max_ion_displacements = np.max(np.sum(
                dc ** 2, axis=-1) ** 0.5, axis=1)
            self.max_framework_displacement = \
                np.max(self.max_ion_displacements[framework_indices])

        self.dt = dt
        self.lower_bound = spec_dict['lower_bound']
        self.upper_bound = spec_dict['upper_bound']
        self.lower_bound_index = lower_bound_index
        self.upper_bound_index = upper_bound_index

        self.msd = msd
        self.msd_by_ions = msd_by_ions
        self.msd_component = msd_component
        self.diffusivity = slope / (20 * dim)
        self.diffusivity_components = slope_components / 20

    def get_summary_dict(self, oxidized_specie=None):
        """
        A summary of information
        :param oxidized_specie (str): specie string with oxidation state. If provided or specie in initial
            function is oxidized, it will calculate conductivity based on nernst-einstein relationship.
        :return: dict of diffusion information
            keys: D, D_components, specie, step_skip, temperature, msd, msd_component, dt, time_intervals_number
                  spec_dict
        """
        d = {"diffusivity": self.diffusivity,
             "diffusivity_components": self.diffusivity_components,
             "specie": self.specie,
             "step_skip": self.step_skip,
             "temperature": self.temperature,
             "msd": self.msd,
             "msd_component": self.msd_component,
             "dt": self.dt,
             "time_intervals_number": self.time_intervals_number,
             "spec_dict": self.spec_dict,
             "drift_maximum": self.drift_maximum
              }
        oxi = False
        if oxidized_specie:
            df_sp = Specie.from_string(oxidized_specie)
            oxi = True
        else:
            try:
                df_sp = Specie.from_string(self.specie)
                oxi = True
            except:
                pass
        if oxi:
            factor = get_conversion_factor(self.structure, df_sp, self.temperature)
            d['conductivity'] = factor * self.diffusivity
            d['conductivity_components'] = factor * self.diffusivity_components
            d['conversion_factor'] = factor
            d['oxidation_state'] = df_sp.oxi_state
        return d

    @classmethod
    def from_structures(cls, structures, specie, temperature,
                        time_step, step_skip, time_intervals_number=1000,
                        spec_dict=None):
        """
        Convenient constructor that takes in a list of Structure objects to
        perform diffusion analysis.

        :param structures ([Structure]): list of Structure objects:
        :param specie (str): species string, like Li, Li+
        :param temperature (float): temperature of MD
        :param time_step (float): time step in MD
        :param step_skip (int): Sampling frequency of the displacements (
                time_step is multiplied by this number to get the real time
                between measurements)
        :param time_intervals_number (int): number of time intervals. Default is 1000
                means there are ~1000 time intervals.
        :param spec_dict (dict): spec dict of linear fitting. Default is
                {'lower_bound': 4.5, 'upper_bound': 0.5, 'minimum_msd_diff': 4.5}
                lower_bound is in unit of Angstrom square
                upper_bound is in unit of total time. 0.5 means upper fitting bound is 0.5*t_total
                minimum_msd_diff is in unit of Angstrom square. msd[upper_bound] - msd[lower_bound] should larger
                than minimum_msd_diff to do linear fitting.
        """
        p = []
        for i, s in enumerate(structures):
            if i == 0:
                structure = s
            p.append(np.array(s.frac_coords)[:, None])

        p.insert(0, p[0])
        p = np.concatenate(p, axis=1)
        dp = p[:, 1:] - p[:, :-1]
        dp = dp - np.round(dp)
        f_disp = np.cumsum(dp, axis=1)
        disp = structure.lattice.get_cartesian_coords(f_disp)

        return cls(structure, disp, specie, temperature,
                   time_step, step_skip=step_skip,
                   time_intervals_number=time_intervals_number,
                   spec_dict=spec_dict)

    @classmethod
    def from_vaspruns(cls, vaspruns, specie,
                      time_intervals_number=1000,
                      spec_dict=None):
        """
        Convenient constructor that takes in a list of Vasprun objects to
        perform diffusion analysis.

        :param vaspruns ([Vasprun]): List of Vaspruns (ordered):
        :param specie (str): species string, like Li, Li+
        :param time_intervals_number (int): number of time intervals. Default is 1000
                means there are ~1000 time intervals.
        :param spec_dict (dict): spec dict of linear fitting. Default is
                {'lower_bound': 4.5, 'upper_bound': 0.5, 'minimum_msd_diff': 4.5}
                lower_bound is in unit of Angstrom square
                upper_bound is in unit of total time. 0.5 means upper fitting bound is 0.5*t_total
                minimum_msd_diff is in unit of Angstrom square. msd[upper_bound] - msd[lower_bound] should larger
                than minimum_msd_diff to do linear fitting.
        """

        def get_structures(vaspruns):
            for i, vr in enumerate(vaspruns):
                if i == 0:
                    step_skip = vr.ionic_step_skip or 1
                    final_structure = vr.initial_structure
                    temperature = vr.parameters['TEEND']
                    time_step = vr.parameters['POTIM']
                    yield step_skip, temperature, time_step
                # check that the runs are continuous
                fdist = pbc_diff(vr.initial_structure.frac_coords,
                                 final_structure.frac_coords)
                if np.any(fdist > 0.001):
                    raise ValueError('initial and final structures do not '
                                     'match.')
                final_structure = vr.final_structure

                assert (vr.ionic_step_skip or 1) == step_skip
                for s in vr.ionic_steps:
                    yield s['structure']

        s = get_structures(vaspruns)
        step_skip, temperature, time_step = next(s)

        return cls.from_structures(structures=s, specie=specie,
                                   temperature=temperature, time_step=time_step, step_skip=step_skip,
                                   time_intervals_number=time_intervals_number, spec_dict=spec_dict)

    @classmethod
    def from_files(cls, filepaths, specie, step_skip=10, ncores=None,
                   time_intervals_number=1000,
                   spec_dict=None):
        """
        Convenient constructor that takes in a list of vasprun.xml paths to
        perform diffusion analysis.

        :param filepaths ([str]): List of paths to vasprun.xml files of runs, ordered.
        :param specie (str): species string, like Li, Li+
        :param step_skip (int): Sampling frequency of the displacements (
                time_step is multiplied by this number to get the real time
                between measurements)
        :param ncores (int): Numbers of cores to use for multiprocessing. Can
                speed up vasprun parsing considerably. Defaults to None,
                which means serial. It should be noted that if you want to
                use multiprocessing, the number of ionic steps in all vasprun
                .xml files should be a multiple of the ionic_step_skip.
                Otherwise, inconsistent results may arise. Serial mode has no
                such restrictions.
        :param time_intervals_number (int): number of time intervals. Default is 1000
                means there are ~1000 time intervals.
        :param spec_dict (dict): spec dict of linear fitting. Default is
                {'lower_bound': 4.5, 'upper_bound': 0.5, 'minimum_msd_diff': 4.5}
                lower_bound is in unit of Angstrom square
                upper_bound is in unit of total time. 0.5 means upper fitting bound is 0.5*t_total
                minimum_msd_diff is in unit of Angstrom square. msd[upper_bound] - msd[lower_bound] should larger
                than minimum_msd_diff to do linear fitting.
        """
        if ncores is not None and len(filepaths) > 1:
            import multiprocessing
            p = multiprocessing.Pool(ncores)
            vaspruns = p.imap(_get_vasprun,
                              [(fp, step_skip) for fp in filepaths])
            analyzer = cls.from_vaspruns(vaspruns, specie=specie,
                                         time_intervals_number=time_intervals_number,
                                         spec_dict=spec_dict)
            p.close()
            p.join()
            return analyzer
        else:
            def vr(filepaths):
                offset = 0
                for p in filepaths:
                    v = Vasprun(p, ionic_step_offset=offset,
                                ionic_step_skip=step_skip)
                    yield v
                    # Recompute offset.
                    offset = (-(v.nionic_steps - offset)) % step_skip

            return cls.from_vaspruns(vr(filepaths), specie=specie,
                                     time_intervals_number=time_intervals_number,
                                     spec_dict=spec_dict)

    def as_dict(self):
        return {
            "@module": self.__class__.__module__,
            "@class": self.__class__.__name__,
            "structure": self.structure.as_dict(),
            "displacements": self.disp.tolist(),
            "specie": self.specie,
            "temperature": self.temperature,
            "time_step": self.time_step,
            "step_skip": self.step_skip,
            "time_intervals_number": self.time_intervals_number,
            "spec_dict": self.spec_dict}

    @classmethod
    def from_dict(cls, d):
        structure = Structure.from_dict(d["structure"])
        return cls(structure, np.array(d["displacements"]), specie=d["specie"],
                   temperature=d["temperature"], time_step=d["time_step"],
                   step_skip=d["step_skip"], time_intervals_number=d["time_intervals_number"],
                   spec_dict=d['spec_dict'])

class ErrorAnalysisFromDiffusivityAnalyzer(object):
    def __init__(self, diffusivity_analyzer, site_distance=3.0):
        """
        Estimate the relative standard deviation (RSD) of D from the equation:
        RSD = 3.43/sqrt(N_jump) + 0.04

        :param diffusivity_analyzer (DiffusivityAnalyzer object):
        :param site_distance (float): the site distance between diffusion ions (averagely)
        """
        n_jump = len(diffusivity_analyzer.indices) * \
                 np.max(diffusivity_analyzer.msd) / (site_distance * site_distance)
        n_jump_component = len(diffusivity_analyzer.indices) * \
                           np.max(diffusivity_analyzer.msd_component, axis=1) / (site_distance * site_distance)
        RSD_D = 3.43 / np.sqrt(n_jump) + 0.04
        RSD_D_component = [None, None, None]
        for i in range(3):
            RSD_D_component[i] = 3.43 / np.sqrt(n_jump_component[i]) + 0.04
        self.diffusivity_analyzer = diffusivity_analyzer
        self.n_jump = n_jump
        self.n_jump_component = n_jump_component
        self.RSD_D = RSD_D
        self.RSD_D_component = np.array(RSD_D_component)

    def get_summary_dict(self, oxidized_specie=None):
        """
        A summary of information
        :param oxidized_specie (str): specie string with oxidation state. If provided or specie in initial
            function is oxidized, it will calculate conductivity based on nernst-einstein relationship.
        :return: dict of diffusion information
        """
        d = self.diffusivity_analyzer.get_summary_dict(oxidized_specie=oxidized_specie)
        d['n_jump'] = self.n_jump
        d['n_jump_component'] = self.n_jump_component
        d['diffusivity_relative_standard_deviation'] = self.RSD_D
        d['diffusivity_standard_deviation'] = self.RSD_D * d['diffusivity']
        d['diffusivity_component_relative_standard_deviation'] = self.RSD_D_component
        d['diffusivity_component_relative_standard_deviation'] = self.RSD_D_component * d['diffusivity_components']
        return d

def _get_vasprun(args):
    """
    Internal method to support multiprocessing.
    """
    return Vasprun(args[0], ionic_step_skip=args[1],
                   parse_dos=False, parse_eigen=False)

class ArreheniusAnalyzer(object):
    def __init__(self, temperatures, diffusivities, diffusivity_errors=None):
        """
        Fitting arrehenius relationship from temperatures, diffusivities and diffusivity_error
        :param temperatures (List): list of temperatures
        :param diffusivities (List): list of diffusivities at different temperatures
        :param diffusivity_errors (List): optional, list of diffusivity error at different temperatures
        """
        slope_to_eV = -8.617e-5 * 1000 * np.log(10)

        def linear(x, k, b):
            return k * x + b

        x = np.array([1000.0 / i for i in temperatures])
        y = np.array([np.log10(i) for i in diffusivities])
        if diffusivity_errors is None:
            [slope, intercept], cov = curve_fit(linear, x, y)
            slope_sigma = np.sqrt(np.diag(cov))[0]
            intercept_sigma = np.sqrt(np.diag(cov))[1]
            y_error = None
        else:
            y_error = [np.log10(np.e) * diffusivity_errors[i] / diffusivities[i] for i in range(len(diffusivities))]
            [slope, intercept], cov = curve_fit(linear, x, y, sigma=y_error, absolute_sigma=True)
            slope_sigma = np.sqrt(np.diag(cov))[0]
            intercept_sigma = np.sqrt(np.diag(cov))[1]

        self.temperatures = temperatures
        self.diffusivities = diffusivities
        self.x = x  # 1000/T
        self.y = y  # log10(D)
        self.diffusivity_errors = diffusivity_errors
        self.y_error = y_error
        self.Ea = slope_to_eV * slope
        self.Ea_error = -1 * slope_to_eV * slope_sigma
        self.intercept = intercept
        self.intercept_sigma = intercept_sigma
        self.slope = slope
        self.slope_sigma = slope_sigma

    def predict_diffusivity(self, temperature):
        """

        :param temperature (float): target temperature
        :return: corresponding diffusivity, and diffusivity range based on error of linear fitting
        """
        logD = self.slope * (1000.0 / temperature) + self.intercept
        logD_sigma = np.sqrt(np.power(self.slope_sigma * (1000.0 / temperature), 2) + \
                             np.power(self.intercept_sigma, 2))
        D_min = np.power(10, logD - logD_sigma)
        D_max = np.power(10, logD + logD_sigma)
        return np.power(10, logD), [D_min, D_max]

    def predict_conductivity(self, temperature, structure, specie):
        """

        :param temperature (float): target temperature
        :param structure (Structure): one structure, used to calculate convertion factor from diffusivity
            to conductivity
        :param specie (str): string of diffusion specie, should contain oxidation state, such as Li+, O2-
        :return: corresponding conductivity, and conductivity range based on error of linear fitting
        """
        D, [D_min, D_max] = self.predict_diffusivity(temperature)
        factor = get_conversion_factor(structure, specie, temperature)
        return factor * D, [factor * D_min, factor * D_max]

    def get_arrhenius_plot(self):
        from pymatgen.util.plotting import pretty_plot
        plt = pretty_plot(12, 8)
        arr = np.power(10, self.slope * self.x + self.intercept)
        plt.plot(self.x, self.diffusivities, 'ko', self.x, arr, 'k--', markersize=10)
        plt.errorbar(self.x, self.diffusivities, yerr=self.diffusivity_errors,
                     fmt='ko', ecolor='k', capthick=2, linewidth=2)
        ax = plt.axes()
        ax.set_yscale('log')
        plt.text(0.6, 0.85, "E$_a$ = {:.3f} eV".format(self.Ea),
                 fontsize=30, transform=plt.axes().transAxes)
        plt.ylabel("D (cm$^2$/s)")
        plt.xlabel("1000/T (K$^{-1}$)")
        plt.tight_layout()
        return plt

    @classmethod
    def from_csv(cls, csv_file):
        with open(csv_file, 'r') as csvf:
            a = csv.reader(csvf, delimiter=str(","))
            data_keys = []
            data_list = []
            for row_index, row in enumerate(a):
                if row_index == 0:
                    data_keys = row
                    data_list = [[] for _ in range(len(row))]
                else:
                    for col_index, col in enumerate(row):
                        data_list[col_index].append(float(col))
        data_dict = {data_keys[i]: data_list[i] for i in range(len(data_keys))}
        if set(data_keys) >= set(['T', 'D', 'D_error']):
            return cls(data_dict['T'], data_dict['D'], data_dict['D_error'])
        elif set(data_keys) >= set(['T', 'D']):
            return cls(data_dict['T'], data_dict['D'])
        else:
            raise Exception("Please make sure the header name in csv file is T, D, D_error(optional)")

def get_conversion_factor(structure, specie, temperature):
    """
    Conversion factor to convert between cm^2/s diffusivity measurements and
    mS/cm conductivity measurements based on number of atoms of diffusing
    species.
    :param structure (Structure): Input structure.
    :param specie (string/specie): Diffusing species string, must contain oxidation state.
    :param temperature (float): Temperature of the diffusion run in Kelvin.
    :return: Conversion factor.
        Conductivity (in mS/cm) = Conversion Factor * Diffusivity (in cm^2/s)
    """
    if type(specie) is Specie:
        df_sp = specie
    else:
        try:
            df_sp = Specie.from_string(specie)
        except:
            raise Exception("Please provide oxidation decorated specie, like Li+, O2-")
    z = df_sp.oxi_state
    el, occu = structure.composition.items()[0]
    if isinstance(el, Specie):  # oxidation decorated structure
        n = structure.composition[specie]
    else:
        n = structure.composition[str(df_sp.element)]
    if n == 0:
        raise Exception("No specie {} in the structure composition: {}".format(specie, structure.composition))
    vol = structure.volume * 1e-24  # units cm^3
    N_A = 6.022140857e+23
    e = 1.6021766208e-19
    R = 8.3144598
    return 1000 * n / (vol * N_A) * z ** 2 * (N_A * e) ** 2 \
           / (R * temperature)

max_stars_repo_name: sivonxay/aimd

max_issues_repo_name: sivonxay/aimd

max_forks_repo_name: sivonxay/aimd


text:
import numpy as np
import pickle
from game import UNO

from agents.greedyAgent import GreedyAgent

all_actions = []
all_prev = []
all_hands = []

agent_1 = GreedyAgent()
agent_2 = GreedyAgent()

for i in range(10000):
  action_list = [[],[]]
  hands_list = [[],[]]
  prev_list = [[],[]]

  uno = UNO([agent_1, agent_2], "finite_deck", forced_draw_only=False)
  while True:
    while True:
      action = uno.get_action()
      current_player = uno.current_player
      hands_list[current_player].append(uno.hands[current_player].copy())
      prev_list[current_player].append(uno.previous_card)
      action_list[current_player].append(action)
      if uno.apply_action(action):
        break
    if uno.current_win():
      winner = uno.current_player
      break
    uno.next_player()
    uno.penalize()
  all_actions.append(action_list[0])
  all_actions.append(action_list[1])
  all_prev.append(prev_list[0])
  all_prev.append(prev_list[1])
  all_hands.append(hands_list[0])
  all_hands.append(hands_list[1])

with open('hands.pkl','wb') as f:
  pickle.dump(all_hands, f)

with open('prev.pkl','wb') as f:
  pickle.dump(all_prev, f)

with open('actions.pkl','wb') as f:
  pickle.dump(all_actions, f)

max_stars_repo_name: huiwy/UNO-Agents

max_issues_repo_name: huiwy/UNO-Agents

max_forks_repo_name: huiwy/UNO-Agents


text:
"""Various support methods for tests"""

import functools
import numpy as np
import sys
import unittest as ut
try:
    from unittest.case import _ExpectedFailure as ExpectedFailure
except ImportError:
    ExpectedFailure = None

def stochastic(func):
    """Decorator to mark test as stochastic - these tests are not fatal when they fail."""
    @functools.wraps(func)
    def wrapper(*args, **kwargs):
        try:
            return func(*args, **kwargs)
        except ut.TestCase.failureException:
            if ExpectedFailure is not None:  # added in Py 2.7
                raise ExpectedFailure(sys.exc_info())
    wrapper.__doc__ += ' (stochastic, failures ignored)'
    return wrapper

class PbTestCase(ut.TestCase):
    """Test case that adds some numeric assert functions"""

    def assertApproxEqual(self, X, Y):
        """Return true if X = Y to within machine precision

        Function for checking that different matrices from different
        computations are in some sense "equal" in the verification tests.
        """
        X = np.asarray(X)
        Y = np.asarray(Y)
        fuzz = 1.0e-8

        self.assertEqual(X.ndim, Y.ndim)
        self.assertEqual(X.shape, Y.shape)

        if np.all(X == Y):  # catches -inf == -inf etc.
            return
        if np.all(abs(X - Y) < fuzz):
            return
        self.fail("NumPy arrays {0} and {1} are not fuzzy equal (+- {2})".format(X, Y, fuzz))

    def assertArraysEqualNotSame(self, a, b):
        """Assert that numpy arrays a and b are equal, but are not the same instances"""
        self.assertNotEqual(id(a), id(b))
        self.assertApproxEqual(a, b)

    def assertArraysSame(self, a, b):
        # id(a) == id(b) doen't work for Cython memoryview arrays
        self.assertEqual(a.ndim, b.ndim)
        self.assertEqual(a.shape, b.shape)
        if a.ndim == 1:
            if a.shape[0] == 0:
                return  # hard to test in this case
            orig_a, orig_b = a[0], b[0]
            a[0] = b[0] = 0  # reset
            a[0] = 1  # does it propagate to b?
            self.assertEqual(a[0], b[0])
            a[0], b[0] = orig_a, orig_b  # set back
        elif a.ndim == 2:
            if a.shape[0] == 0 or a.shape[1] == 0:
                return  # hard to test
            orig_a, orig_b = a[0, 0], b[0, 0]
            a[0, 0] = b[0, 0] = 0  # reset
            a[0, 0] = 1  # does it propagate to b?
            self.assertEqual(a[0, 0], b[0, 0])
            a[0, 0], b[0, 0] = orig_a, orig_b  # set back
        else:
            self.fail("More than 2 or less than 1 dimensions not supported in this method")

    def assertRVsEqualNotSame(self, a, b):
        """Assert that :class:`~pybayes.pdfs.RV` objects a and b are equal, but
        are not the same instances neither shallow copies of themselves.

        RVs are special case during deepcopy - the RVComps should be referenced,
        not copied."""
        self.assertNotEqual(id(a), id(b))
        self.assertNotEqual(id(a.components), id(b.components))
        self.assertEqual(a.name, b.name)  # no need to test for id inequality - strings are immutable
        self.assertEqual(a.dimension, b.dimension)  # ditto
        for (a_comp, b_comp) in zip(a.components, b.components):
            # equality for rv comps is defined as object instance identity
            self.assertEqual(a_comp, b_comp)

max_stars_repo_name: strohel/PyBayes

max_issues_repo_name: strohel/PyBayes

max_forks_repo_name: strohel/PyBayes


text:
# -*- coding: utf-8 -*-

import sys
import math
import tflearn
import tensorflow as tf
from tensorflow.python.ops import rnn_cell
from tensorflow.python.ops import rnn
import chardet
import numpy as np
import struct

question_seqs = []
answer_seqs = []

max_w = 50
float_size = 4
word_vector_dict = {}
word_vec_dim = 200
max_seq_len = 8
word_set = {}

def load_word_set():
    file_object = open('./segment_result_lined.3000000.pair.less', 'r')
    while True:
        line = file_object.readline()
        if line:
            line_pair = line.split('|')
            line_question = line_pair[0]
            line_answer = line_pair[1]
            for word in line_question.decode('utf-8').split(' '):
                word_set[word] = 1
            for word in line_answer.decode('utf-8').split(' '):
                word_set[word] = 1
        else:
            break
    file_object.close()

def load_vectors(input):
    """从vectors.bin加载词向量，返回一个word_vector_dict的词典，key是词，value是200维的向量
    """
    print "begin load vectors"

    input_file = open(input, "rb")

    # 获取词表数目及向量维度
    words_and_size = input_file.readline()
    words_and_size = words_and_size.strip()
    words = long(words_and_size.split(' ')[0])
    size = long(words_and_size.split(' ')[1])
    print "words =", words
    print "size =", size

    for b in range(0, words):
        a = 0
        word = ''
        # 读取一个词
        while True:
            c = input_file.read(1)
            word = word + c
            if False == c or c == ' ':
                break
            if a < max_w and c != '\n':
                a = a + 1
        word = word.strip()

        vector = []
        for index in range(0, size):
            m = input_file.read(float_size)
            (weight,) = struct.unpack('f', m)
            vector.append(float(weight))

        # 将词及其对应的向量存到dict中

        if word_set.has_key(word.decode('utf-8')):
            word_vector_dict[word.decode('utf-8')] = vector[0:word_vec_dim]

    input_file.close()

    print "load vectors finish"

def init_seq(input_file):
    """读取切好词的文本文件，加载全部词序列
    """
    file_object = open(input_file, 'r')
    vocab_dict = {}
    while True:
        question_seq = []
        answer_seq = []
        line = file_object.readline()
        if line:
            line_pair = line.split('|')
            line_question = line_pair[0]
            line_answer = line_pair[1]
            for word in line_question.decode('utf-8').split(' '):
                if word_vector_dict.has_key(word):
                    question_seq.append(word_vector_dict[word])
            for word in line_answer.decode('utf-8').split(' '):
                if word_vector_dict.has_key(word):
                    answer_seq.append(word_vector_dict[word])
        else:
            break
        question_seqs.append(question_seq)
        answer_seqs.append(answer_seq)
    file_object.close()

def vector_sqrtlen(vector):
    len = 0
    for item in vector:
        len += item * item
    len = math.sqrt(len)
    return len

def vector_cosine(v1, v2):
    if len(v1) != len(v2):
        sys.exit(1)
    sqrtlen1 = vector_sqrtlen(v1)
    sqrtlen2 = vector_sqrtlen(v2)
    value = 0
    for item1, item2 in zip(v1, v2):
        value += item1 * item2
    return value / (sqrtlen1*sqrtlen2)

def vector2word(vector):
    max_cos = -10000
    match_word = ''
    for word in word_vector_dict:
        v = word_vector_dict[word]
        cosine = vector_cosine(vector, v)
        if cosine > max_cos:
            max_cos = cosine
            match_word = word
    return (match_word, max_cos)

class MySeq2Seq(object):
    """
    思路：输入输出序列一起作为input，然后通过slick和unpack切分
    完全按照论文说的编码器解码器来做
    输出的时候把解码器的输出按照词向量的200维展平，这样输出就是(?,seqlen*200)
    这样就可以通过regression来做回归计算了，输入的y也展平，保持一致
    """
    def __init__(self, max_seq_len = 16, word_vec_dim = 200, input_file='./segment_result_lined.3000000.pair.less'):
        self.max_seq_len = max_seq_len
        self.word_vec_dim = word_vec_dim
        self.input_file = input_file

    def generate_trainig_data(self):
        load_word_set()
        load_vectors("./vectors.bin")
        init_seq(self.input_file)
        xy_data = []
        y_data = []
        for i in range(len(question_seqs)):
        #for i in range(100):
            question_seq = question_seqs[i]
            answer_seq = answer_seqs[i]
            if len(question_seq) < self.max_seq_len and len(answer_seq) < self.max_seq_len:
                sequence_xy = [np.zeros(self.word_vec_dim)] * (self.max_seq_len-len(question_seq)) + list(reversed(question_seq))
                sequence_y = answer_seq + [np.zeros(self.word_vec_dim)] * (self.max_seq_len-len(answer_seq))
                sequence_xy = sequence_xy + sequence_y
                sequence_y = [np.ones(self.word_vec_dim)] + sequence_y
                xy_data.append(sequence_xy)
                y_data.append(sequence_y)

                #print "right answer"
                #for w in answer_seq:
                #    (match_word, max_cos) = vector2word(w)
                #    if len(match_word)>0:
                #        print match_word, vector_sqrtlen(w)

        return np.array(xy_data), np.array(y_data)

    def model(self, feed_previous=False):
        # 通过输入的XY生成encoder_inputs和带GO头的decoder_inputs
        input_data = tflearn.input_data(shape=[None, self.max_seq_len*2, self.word_vec_dim], dtype=tf.float32, name = "XY")
        encoder_inputs = tf.slice(input_data, [0, 0, 0], [-1, self.max_seq_len, self.word_vec_dim], name="enc_in")
        decoder_inputs_tmp = tf.slice(input_data, [0, self.max_seq_len, 0], [-1, self.max_seq_len-1, self.word_vec_dim], name="dec_in_tmp")
        go_inputs = tf.ones_like(decoder_inputs_tmp)
        go_inputs = tf.slice(go_inputs, [0, 0, 0], [-1, 1, self.word_vec_dim])
        decoder_inputs = tf.concat(1, [go_inputs, decoder_inputs_tmp], name="dec_in")

        # 编码器
        # 把encoder_inputs交给编码器，返回一个输出(预测序列的第一个值)和一个状态(传给解码器)
        (encoder_output_tensor, states) = tflearn.lstm(encoder_inputs, self.word_vec_dim, return_state=True, scope='encoder_lstm')
        encoder_output_sequence = tf.pack([encoder_output_tensor], axis=1)

        # 解码器
        # 预测过程用前一个时间序的输出作为下一个时间序的输入
        # 先用编码器的最后一个输出作为第一个输入
        if feed_previous:
            first_dec_input = go_inputs
        else:
            first_dec_input = tf.slice(decoder_inputs, [0, 0, 0], [-1, 1, self.word_vec_dim])
        decoder_output_tensor = tflearn.lstm(first_dec_input, self.word_vec_dim, initial_state=states, return_seq=False, reuse=False, scope='decoder_lstm')
        decoder_output_sequence_single = tf.pack([decoder_output_tensor], axis=1)
        decoder_output_sequence_list = [decoder_output_tensor]
        # 再用解码器的输出作为下一个时序的输入
        for i in range(self.max_seq_len-1):
            if feed_previous:
                next_dec_input = decoder_output_sequence_single
            else:
                next_dec_input = tf.slice(decoder_inputs, [0, i+1, 0], [-1, 1, self.word_vec_dim])
            decoder_output_tensor = tflearn.lstm(next_dec_input, self.word_vec_dim, return_seq=False, reuse=True, scope='decoder_lstm')
            decoder_output_sequence_single = tf.pack([decoder_output_tensor], axis=1)
            decoder_output_sequence_list.append(decoder_output_tensor)

        decoder_output_sequence = tf.pack(decoder_output_sequence_list, axis=1)
        real_output_sequence = tf.concat(1, [encoder_output_sequence, decoder_output_sequence])

        net = tflearn.regression(real_output_sequence, optimizer='sgd', learning_rate=0.1, loss='mean_square')
        model = tflearn.DNN(net)
        return model

    def train(self):
        trainXY, trainY = self.generate_trainig_data()
        model = self.model(feed_previous=False)
        model.fit(trainXY, trainY, n_epoch=1000, snapshot_epoch=False, batch_size=1)
        model.save('./model/model')
        return model

    def load(self):
        model = self.model(feed_previous=True)
        model.load('./model/model')
        return model

if __name__ == '__main__':
    phrase = sys.argv[1]
    if 3 == len(sys.argv):
        my_seq2seq = MySeq2Seq(word_vec_dim=word_vec_dim, max_seq_len=max_seq_len, input_file=sys.argv[2])
    else:
        my_seq2seq = MySeq2Seq(word_vec_dim=word_vec_dim, max_seq_len=max_seq_len)
    if phrase == 'train':
        my_seq2seq.train()
    else:
        model = my_seq2seq.load()
        trainXY, trainY = my_seq2seq.generate_trainig_data()
        predict = model.predict(trainXY)
        for sample in predict:
            print "predict answer"
            for w in sample[1:]:
                (match_word, max_cos) = vector2word(w)
                #if vector_sqrtlen(w) < 1:
                #    break
                print match_word, max_cos, vector_sqrtlen(w)

max_stars_repo_name: drpreetyrai/ChatBotCourse


text:
#PyQt imports
from PyQt5.QtWidgets import QApplication, QWidget, QPushButton, QLineEdit
from PyQt5.QtCore import QThread

import sys
import os
from time import sleep

import pyrealsense2 as rs;

import numpy as np;

import cv2 as cv;

import pickle

import gzip

import threading

from time import sleep

class WorkerThread(QThread):
    def __init__(self, containerObject):
        QThread.__init__(self)
        self.containerObject = containerObject

    def run(self):
        while(True):
            if (self.containerObject.captureButton.text() ==  'Stop Capturing'):
                self.containerObject.captureImage()
                sleep(3)
            else:
                sleep(.5)

class Container(QWidget):

    def __init__(self):
        super().__init__();
        print('Loading model')
        self.model = self.loadModel()
        print('Model loaded')
        self.initUI();

    def initUI(self):
        # Capture button to initialize capturing of data.
        captureButton = QPushButton('Capture Data', self);
        captureButton.setToolTip('Use this button to capture an image from camera');
        captureButton.resize(150, 100);
        captureButton.move(300, 0);
        captureButton.clicked.connect(self.onCaptureClick);

    def onCaptureClick(self):
        print(3)
        sleep(1)
        print(2)
        sleep(1)
        print(1)
        sleep(1)
        print('Capturing')
        self.captureImage()

    def loadModel(self):
        return pickle.load(gzip.open('trained_model.pklz'))

    def removeBackground(self, image):
        if (np.max(image) != 0):

            minVal = np.min(image[np.nonzero(image)])

            image[image > 1500 + minVal] = 0

        return image

    def classifyImage(self, image):
        print(self.model.predict(image.copy().flatten().reshape(1, -1)))

    def captureImage(self):
        # Create a pipeline
        pipeline = rs.pipeline()

        #Create a config and configure the pipeline to stream for depth at 30fps
        config = rs.config()
        config.enable_stream(rs.stream.depth, 640, 480, rs.format.z16, 30)

        # Start streaming pipeline
        profile = pipeline.start(config)

        # Getting depth sensor
        depth_sensor = profile.get_device().first_depth_sensor()
        depth_scale = depth_sensor.get_depth_scale()

        # Set a clipping distance of 1 meter
        clipping_distance_in_meters = 1 #1 meter
        clipping_distance = clipping_distance_in_meters / depth_scale

        # Align stream to depth
        align_to = rs.stream.color
        align = rs.align(align_to)

        # Get frame of depth
        frames = pipeline.wait_for_frames();
        aligned_frames = align.process(frames)
        aligned_depth_frame = aligned_frames.get_depth_frame()

        # Compile a depth image
        depth_image = np.asanyarray(aligned_depth_frame.get_data())

        # Remove background
        depth_image = self.removeBackground(depth_image)

        # Turn image into three channel image
        depth_image_3d = np.dstack((depth_image,depth_image,depth_image))

        # Create a colormap depth image
        depth_colormap = cv.applyColorMap(cv.convertScaleAbs(depth_image, alpha=0.03), cv.COLORMAP_JET)

        # Create gray image so we can find contours
        depth_gray = cv.cvtColor(depth_colormap, cv.COLOR_BGR2GRAY)
        canny_output = cv.Canny(depth_gray, 100, 200);
        contours, hierarchy = cv.findContours(canny_output, cv.RETR_TREE, cv.CHAIN_APPROX_SIMPLE)

        # Create copy to draw contours on
        contour_map = depth_colormap.copy()
        contour_map = np.zeros((canny_output.shape[0], canny_output.shape[1], 3), dtype=np.uint8)
        # Draw contours
        for i in range(len(contours)):
            color = (0, 255, 0)
            cv.drawContours(contour_map, contours, i, color, 2, cv.LINE_8, hierarchy, 0)

        self.classifyImage(contour_map)

        # Draw images side by side
        images = np.hstack((depth_colormap, contour_map))
        cv.namedWindow('Contour Example', cv.WINDOW_AUTOSIZE)
        cv.imshow('Contour Example', images)
        key = cv.waitKey(1)

def main():
    app = QApplication(sys.argv);
    window = Container();
    window.setGeometry(0, 0, 450, 100);
    window.setWindowTitle("Gesture Recognition Training GUI");
    window.show();
    sys.exit(app.exec_());

if __name__ == "__main__":
    main();

max_stars_repo_name: ShaneClancy/gesture-recognition

max_issues_repo_name: ShaneClancy/gesture-recognition

max_forks_repo_name: ShaneClancy/gesture-recognition