text:
# Analyze the distances of the positive and negative examples and see if the
# classification approaches can correctly classify concepts that further away
# from the centroid than a correctly identified negative example / closer to the
# centroid

from evaluation import load_results, get_truth
from load_models import get_cosine, load_model, get_nearest_neighbors
from utils import load_vecs, merge_wi_dicts, load_data, load_gold
#from nearest_neighbors import get_centroid
import pandas as pd
import glob
import numpy as np
import os
import sys

def load_cosines_centroid(model_name, feature):

    cosine_dict = dict()

    if os.path.isfile('../results/'+model_name+'/cosine_distances/'+feature+'.txt'):

        with open('../results/'+model_name+'/cosine_distances/'+feature+'.txt') as infile:
            lines = infile.read().strip().split('\n')

        for line in lines:
            word, cos, st = line.strip().split(',')

            if cos == 'nan':
                cosine_dict[word] = np.NaN
            else:
                cosine_dict[word] = float(cos)

        return cosine_dict
    else:
        return None

def collect_answers(feature, model_name):

    prediction_dict = dict()

    words_pos, words_neg = load_gold(feature)
    words = words_pos + words_neg
    truth = get_truth(feature)

    print(feature+','+str(len(words_pos))+','+str(len(words_neg)))

    # load all predictions
    results_paths  = glob.glob('../results/'+model_name+'/*/*/'+feature+'.txt')

    prediction_dict['gold'] = truth
    prediction_dict['words'] = words

    for res in results_paths:
        name = 'pred-'+'-'.join(res.split('/')[2:-1])

        if os.path.isfile(res):

            predictions = load_results(res)
            if predictions != None:
                prediction_dict[name] = predictions
            else:
                prediction_dict[name] = ['-' for w in words]

        else:
            prediction_dict[name] = ['-' for w in words]

    return prediction_dict

def create_cosine_prediction_overview(feature, model_name):

    predictions = collect_answers(feature, model_name)

    cosine_dict = load_cosines_centroid(model_name, feature)

    if cosine_dict != None:

        cosines = [str(cosine_dict[word]) for word  in predictions['words']]

        predictions['cosines'] = cosines

        df = pd.DataFrame.from_dict(predictions).set_index('words')

        sorted_df = df.sort_index(axis = 1)
        # sorted_df=unsorted_df.sort_index(axis=1)

        return sorted_df.sort_values('cosines', ascending = False)
    else:
        return '-'

def main():

    model_name = sys.argv[1]
    feature = sys.argv[2]

    print('feature, positive, negative')
    if not os.path.isdir('../analyses/'):
        os.mkdir('../analyses/')

    if feature == 'all':
        print('all')
        files = glob.glob('../gold/*-pos.txt')
        features = [f.split('/')[-1].split('.')[0].split('-')[0] for f in files]

    else:
        features = [feature]

    for feat in features:

        overview = create_cosine_prediction_overview(feat, model_name)

        if type(overview) != str:

            overview.to_csv('../analyses/'+model_name+'-'+feat+'.csv')

if __name__ == '__main__':

    main()

ext: py

reason: import numpy


text:
# Copyright 2019 United Kingdom Research and Innovation
# Author: Evgueni Ovtchinnikov (evgueni.ovtchinnikov@stfc.ac.uk)

"""MKL loader/wrapper.
"""

import ctypes
import multiprocessing
import numpy
import os
import site
import sys
from sys import platform

from . import verbosity

if platform == 'win32':
    mkl_name = 'mkl_rt.dll'
    ld_path = 'PATH'
else:
    mkl_name = 'libmkl_rt.so'
    ld_path = 'LD_LIBRARY_PATH'

def find(name, path):
    for root, dirs, files in os.walk(path):
        if name in files:
            return os.path.join(root, name)

def find_mkl():
    from . import env
    mkl_path = env.mkl_path
    if mkl_path is not None:
        return os.path.join(mkl_path, mkl_name)
    user_base = site.USER_BASE
    mkl = find(mkl_name, user_base)
    if mkl is not None:
        #print('MKL found in USER_BASE')
        return mkl
    v = sys.version_info
    major = repr(v[0])
    minor = repr(v[1])
    suff1 = 'ython' + major + minor
    suff2 = 'ython' + major + '.' + minor
    for path in sys.path:
        i = path.find(suff1)
        if i < 0:
            i = path.find(suff2)
            if i < 0:
                continue
        i = path.find('-packages')
        if i < 0:
            continue
        #print('searching in %s...' % path)
        mkl = find(mkl_name, path)
        return mkl
    return None

def _array_ptr(array, shift=0):
    return ctypes.c_void_p(array.ctypes.data + shift)

def mkl_version(mkl):
    str_v = numpy.ndarray((256,), dtype=numpy.uint8)
    ptr_v = _array_ptr(str_v)
    len_v = ctypes.c_int(256)
    mkl.mkl_get_version_string(ptr_v, len_v)
    return str_v.tostring().decode('ascii')

mkl_path = find_mkl()
if mkl_path is not None:
    if verbosity.level > 1:
        print('loading %s...' % mkl_path)
    mkl = ctypes.CDLL(mkl_path, mode=ctypes.RTLD_GLOBAL)
else:
    if verbosity.level > 1:
        print('trying to load %s...' % mkl_name)
    try:
        mkl = ctypes.CDLL(mkl_name, mode=ctypes.RTLD_GLOBAL)
    except:
        if verbosity.level > 0:
            print('\nCould not find ' + mkl_name + ' in folders listed in ' \
                  + ld_path + '.')
            print('\nIf you have MKL 10.3 or later, make sure the folder ' + \
                  'containing ' + mkl_name + ' is listed in ' + ld_path + \
                  ', or else do')
            print('    from raleigh.algebra import env')
            print('and assign' + \
                  ' the path to that folder to env.mkl_path.')
            print('\nIf you do not have MKL, you can install it using')
            print('    pip install --user mkl')
            print('(RALEIGH will find it, no need to edit ' + ld_path + ').\n')
        raise RuntimeError(mkl_name + ' not found')
if verbosity.level > 1:
    print('Loaded %s' % mkl_version(mkl).strip())

# find the total number of threads
num_cpus = multiprocessing.cpu_count()

# find the number of cores using MKL threading behaviour
mkl.MKL_Set_Dynamic(1)
# now mkl_get_max_threads() returns the number of cores!
num_cores = mkl.mkl_get_max_threads()

# if hyperthreading is used, increase the number of mkl threads
# to achieve slightly better performance
if num_cpus == 2*num_cores:
    num_threads = num_cpus - 1
    mkl.MKL_Set_Dynamic(0)
    mkl.MKL_Set_Num_Threads(num_threads)

if verbosity.level > 1:
    print('Using %d MKL threads' % mkl.mkl_get_max_threads())

class Cblas:
    '''MKL cblas wrapper.
    '''

    ColMajor = 102
    NoTrans = 111
    Trans = 112
    ConjTrans = 113

    def __init__(self, dt):
        if dt == numpy.float32:
            self.dsize = 4
            self.gemm = mkl.cblas_sgemm
            self.axpy = mkl.cblas_saxpy
            self.copy = mkl.cblas_scopy
            self.scal = mkl.cblas_sscal
            self.norm = mkl.cblas_snrm2
            self.norm.restype = ctypes.c_float
            self.inner = mkl.cblas_sdot
            self.inner.restype = ctypes.c_float
            self.mkl_one = ctypes.c_float(1.0)
            self.mkl_zero = ctypes.c_float(0.0)
            self.svd = mkl.LAPACKE_sgesvd
        elif dt == numpy.float64:
            self.dsize = 8
            self.gemm = mkl.cblas_dgemm
            self.axpy = mkl.cblas_daxpy
            self.copy = mkl.cblas_dcopy
            self.scal = mkl.cblas_dscal
            self.norm = mkl.cblas_dnrm2
            self.norm.restype = ctypes.c_double
            self.inner = mkl.cblas_ddot
            self.inner.restype = ctypes.c_double
            self.mkl_one = ctypes.c_double(1.0)
            self.mkl_zero = ctypes.c_double(0.0)
            self.svd = mkl.LAPACKE_dgesvd
        elif dt == numpy.complex64:
            self.dsize = 8
            self.gemm = mkl.cblas_cgemm
            self.axpy = mkl.cblas_caxpy
            self.copy = mkl.cblas_ccopy
            self.scal = mkl.cblas_cscal
            self.norm = mkl.cblas_scnrm2
            self.norm.restype = ctypes.c_float
            self.inner = mkl.cblas_cdotc_sub
            self.cmplx_val = numpy.zeros((2,), dtype=numpy.float32)
            self.cmplx_one = numpy.zeros((2,), dtype=numpy.float32)
            self.cmplx_one[0] = 1.0
            self.cmplx_zero = numpy.zeros((2,), dtype=numpy.float32)
            self.mkl_one = ctypes.c_void_p(self.cmplx_one.ctypes.data)
            self.mkl_zero = ctypes.c_void_p(self.cmplx_zero.ctypes.data)
            self.svd = mkl.LAPACKE_cgesvd
        elif dt == numpy.complex128:
            self.dsize = 16
            self.gemm = mkl.cblas_zgemm
            self.axpy = mkl.cblas_zaxpy
            self.copy = mkl.cblas_zcopy
            self.scal = mkl.cblas_zscal
            self.norm = mkl.cblas_dznrm2
            self.norm.restype = ctypes.c_double
            self.inner = mkl.cblas_zdotc_sub
            self.cmplx_val = numpy.zeros((2,), dtype=numpy.float64)
            self.cmplx_one = numpy.zeros((2,), dtype=numpy.float64)
            self.cmplx_one[0] = 1.0
            self.cmplx_zero = numpy.zeros((2,), dtype=numpy.float64)
            self.mkl_one = ctypes.c_void_p(self.cmplx_one.ctypes.data)
            self.mkl_zero = ctypes.c_void_p(self.cmplx_zero.ctypes.data)
            self.svd = mkl.LAPACKE_zgesvd
        else:
            raise ValueError('data type %s not supported' % repr(dt))
        self.__s = 'S'
        self.__o = 'O'
        self.small = ctypes.c_char(self.__s.encode('utf-8'))
        self.ovwrt = ctypes.c_char(self.__o.encode('utf-8'))

class SparseSymmetricMatrix:
    '''MKL csr matrix wrapper.
    '''

    def __init__(self, a, ia, ja):
        self.__a = a
        self.__ib = ia
        self.__ie = ia[1:]
        self.__ja = ja
        self.__n = ctypes.c_int(ia.shape[0] - 1)
        self.__dtype = a.dtype
        self.__oz = numpy.array([1.0, 0.0]).astype(a.dtype)
        self.__one = _array_ptr(self.__oz)
        self.__zero = _array_ptr(self.__oz[1:])
        self.__complex = False
        if a.dtype == numpy.float32:
            self.__csrmm = mkl.mkl_scsrmm
            self.__csrsymv = mkl.mkl_scsrsymv
            self.__descr = 'SUNF  '
        elif a.dtype == numpy.float64:
            self.__csrmm = mkl.mkl_dcsrmm
            self.__csrsymv = mkl.mkl_dcsrsymv
            self.__descr = 'SUNF  '
        elif a.dtype == numpy.complex64:
            self.__csrmm = mkl.mkl_ccsrmm
            self.__csrsymv = mkl.mkl_ccsrsymv
            self.__descr = 'HUNF  '
            self.__complex = True
        elif a.dtype == numpy.complex128:
            self.__csrmm = mkl.mkl_zcsrmm
            self.__csrsymv = mkl.mkl_zcsrsymv
            self.__descr = 'HUNF  '
            self.__complex = True
        else:
            raise ValueError('unsupported data type')
        self.__csrmm.restype = None
        self.__uplo = 'U'
        self.__trans = 'N'
        self.__u = ctypes.c_char_p(self.__uplo.encode('utf-8'))
        self.__t = ctypes.c_char_p(self.__trans.encode('utf-8'))
        self.__d = ctypes.c_char_p(self.__descr.encode('utf-8'))

    def dot(self, x, y):
        ptr_x = _array_ptr(x)
        ptr_y = _array_ptr(y)
        ptr_a = _array_ptr(self.__a)
        ptr_ib = _array_ptr(self.__ib)
        ptr_ie = _array_ptr(self.__ie)
        ptr_ja = _array_ptr(self.__ja)
        s = x.shape
        if (len(s) == 1 or s[0] == 1) and not self.__complex:
            #print('using *csrsymv...')
            if len(s) == 1:
                n = s[0]
            else:
                n = numpy.prod(s[1:])
            n = ctypes.c_int(n)
            self.__csrsymv(self.__u, ctypes.byref(n), \
                           ptr_a, ptr_ib, ptr_ja, ptr_x, ptr_y)
            return
        #print('using *csrmm...')
        if len(s) > 1:
            m, n = s[0], int(numpy.prod(s[1:]))
        else:
            m = 1
            n = self.__n.value
        mn = numpy.array([m, n])
        ptr_m = _array_ptr(mn)
        ptr_n = _array_ptr(mn[1:])
        n = self.__n
        self.__csrmm(self.__t, ctypes.byref(n), ptr_m, ctypes.byref(n), \
            self.__one, self.__d, ptr_a, ptr_ja, ptr_ib, ptr_ie, \
            ptr_x, ptr_n, self.__zero, ptr_y, ptr_n)

class ILUT:
    '''MKL dscrilut wrapper.
    '''

    def __init__(self, a, ia, ja):
        self.__a = a
        self.__ia = ia
        self.__ja = ja
        n = ia.shape[0] - 1
        nnz = ia[n] - ia[0]
        self.__rows = n
        self.__nnz = nnz//n
        self.__ipar = numpy.zeros((128,), dtype=numpy.int32)
        self.__dpar = numpy.zeros((128,))
        self.__w = numpy.zeros((n,))
        self.__ipar[0] = n
        self.__b = None
        self.__ib = None
        self.__jb = None
        self.__u = 'U'
        self.__l = 'L'
        self.__n = 'N'
        self.__U = ctypes.c_char_p(self.__u.encode('utf-8'))
        self.__L = ctypes.c_char_p(self.__l.encode('utf-8'))
        self.__N = ctypes.c_char_p(self.__n.encode('utf-8'))

    def factorize(self, tol=1e-6, max_fill_rel=1):
        max_fill_abs = min(self.__rows - 1, self.__nnz * max_fill_rel)
        self.__ipar[30] = 1
        self.__dpar[30] = tol
        n = self.__rows
        nnz = (2*max_fill_abs + 1)*n
        #print(n, nnz)
        self.__b = numpy.ndarray((nnz,))
        self.__ib = numpy.ndarray((n + 1,))
        self.__jb = numpy.ndarray((nnz,), dtype=numpy.int32)
        n_c = ctypes.c_int(n)
        mf_c = ctypes.c_int(max_fill_abs)
        ierr_c = ctypes.c_int()
        tol_c = ctypes.c_double(tol)
        ptr_a = _array_ptr(self.__a)
        ptr_ia = _array_ptr(self.__ia)
        ptr_ja = _array_ptr(self.__ja)
        ptr_b = _array_ptr(self.__b)
        ptr_ib = _array_ptr(self.__ib)
        ptr_jb = _array_ptr(self.__jb)
        ptr_ipar = _array_ptr(self.__ipar)
        ptr_dpar = _array_ptr(self.__dpar)
        mkl.dcsrilut(ctypes.byref(n_c), \
                     ptr_a, ptr_ia, ptr_ja, ptr_b, ptr_ib, ptr_jb, \
                     ctypes.byref(tol_c), ctypes.byref(mf_c), \
                     ptr_ipar, ptr_dpar, ctypes.byref(ierr_c))
        if ierr_c.value != 0: print(ierr_c.value)

    def solve(self, f, x):
        m, n = f.shape
        n = ctypes.c_int(self.__rows)
        #print(n)
        ptr_b = _array_ptr(self.__b)
        ptr_ib = _array_ptr(self.__ib)
        ptr_jb = _array_ptr(self.__jb)
        ptr_w = _array_ptr(self.__w)
        for i in range(m):
            ptr_f = _array_ptr(f[i,:])
            ptr_x = _array_ptr(x[i,:])
            mkl.mkl_dcsrtrsv(self.__L, self.__N, self.__U, ctypes.byref(n), \
                             ptr_b, ptr_ib, ptr_jb, ptr_f, ptr_w)
            mkl.mkl_dcsrtrsv(self.__U, self.__N, self.__N, ctypes.byref(n), \
                             ptr_b, ptr_ib, ptr_jb, ptr_w, ptr_x)

class ParDiSo:
    '''MKL pardiso wrapper.
    '''

    def __init__(self, dtype=numpy.float64, pos_def=False):
        self.__pardiso = mkl.pardiso
        self.__pt = numpy.ndarray((64,), dtype=numpy.int64)
        self.__iparm = numpy.ndarray((64,), dtype=numpy.int32)
        self.__handle = _array_ptr(self.__pt)
        if (dtype == numpy.float32 or dtype == numpy.float64):
            m = 2
            self.__real = True
        elif (dtype == numpy.complex64 or dtype == numpy.complex128):
            m = 4
            self.__real = False
        else:
            raise ValueError('ParDiSo constructor: wrong data type')
        if not pos_def:
            m = -m
        #print('matrix type: %d' % m)
        mtype = ctypes.c_int(m)
        self.__dtype = dtype
        self.__mtype = mtype
        self.__a = None
        self.__ia = None
        self.__ja = None
        self.__perm = None
        self.__dummy = numpy.ndarray((1,))
        self.__ptr = ctypes.c_void_p(self.__dummy.ctypes.data)
        self.__ptr_iparm = _array_ptr(self.__iparm)
        mkl.pardisoinit(self.__handle, ctypes.byref(mtype), self.__ptr_iparm)
        self.__iparm[4] = 2
        if dtype == numpy.float32 or dtype == numpy.complex64:
            self.__iparm[27] = 1

    def __del__(self):
        step = ctypes.c_int(-1)
        maxf = ctypes.c_int(1)
        mnum = ctypes.c_int(1)
        n = ctypes.c_int(1)
        m = ctypes.c_int(1)
        ptr = self.__ptr
        verb = ctypes.c_int(0)
        err = ctypes.c_long(0)
        self.__pardiso(self.__handle, \
            ctypes.byref(maxf), ctypes.byref(mnum), \
            ctypes.byref(self.__mtype), \
            ctypes.byref(step), ctypes.byref(n), ptr, ptr, ptr, \
            ptr, ctypes.byref(m), self.__ptr_iparm, ctypes.byref(verb), \
            ptr, ptr, ctypes.byref(err))
        if err.value != 0: print(err.value)

    def handle(self):
        return self.__handle

    def iparm(self):
        return self.__iparm

    def perm(self):
        return self.__perm

    def analyse(self, a, ia, ja):
        self.__a = a
        self.__ia = ia
        self.__ja = ja
        step = ctypes.c_int(11)
        maxf = ctypes.c_int(1)
        mnum = ctypes.c_int(1)
        rows = ia.shape[0] - 1
        n = ctypes.c_int(rows)
        m = ctypes.c_int(1)
        ptr_a = _array_ptr(self.__a)
        ptr_ia = _array_ptr(self.__ia)
        ptr_ja = _array_ptr(self.__ja)
        self.__perm = numpy.ndarray((rows,), dtype=numpy.int32)
        ptr_perm = _array_ptr(self.__perm)
        ptr_iparm = _array_ptr(self.__iparm)
        verb = ctypes.c_int(0)
        err = ctypes.c_long(0)
        self.__pardiso(self.__handle, ctypes.byref(maxf), ctypes.byref(mnum), \
                    ctypes.byref(self.__mtype), \
                    ctypes.byref(step), ctypes.byref(n), ptr_a, ptr_ia, ptr_ja, \
                    ptr_perm, ctypes.byref(m), ptr_iparm, ctypes.byref(verb), \
                    self.__ptr, self.__ptr, ctypes.byref(err))
        return err.value

    def factorize(self):
        step = ctypes.c_int(22)
        maxf = ctypes.c_int(1)
        mnum = ctypes.c_int(1)
        rows = self.__ia.shape[0] - 1
        n = ctypes.c_int(rows)
        m = ctypes.c_int(1)
        ptr_a = _array_ptr(self.__a)
        ptr_ia = _array_ptr(self.__ia)
        ptr_ja = _array_ptr(self.__ja)
        ptr_perm = _array_ptr(self.__perm)
        ptr_iparm = _array_ptr(self.__iparm)
        verb = ctypes.c_int(0)
        err = ctypes.c_long(0)
        self.__pardiso(self.__handle, ctypes.byref(maxf), ctypes.byref(mnum), \
                    ctypes.byref(self.__mtype), \
                    ctypes.byref(step), ctypes.byref(n), ptr_a, ptr_ia, ptr_ja, \
                    ptr_perm, ctypes.byref(m), ptr_iparm, ctypes.byref(verb), \
                    self.__ptr, self.__ptr, ctypes.byref(err))
        return err.value

    def solve(self, b, x, part = None):
        if len(b.shape) > 1:
            nrhs = b.shape[0]
        else:
            nrhs = 1
        if part == 'f':
            step = ctypes.c_int(331)
        elif part == 'd':
            step = ctypes.c_int(332)
        elif part == 'b':
            step = ctypes.c_int(333)
        else:
            step = ctypes.c_int(33)
        maxf = ctypes.c_int(1)
        mnum = ctypes.c_int(1)
        rows = self.__ia.shape[0] - 1
        n = ctypes.c_int(rows)
        m = ctypes.c_int(nrhs)
        ptr_a = _array_ptr(self.__a)
        ptr_ia = _array_ptr(self.__ia)
        ptr_ja = _array_ptr(self.__ja)
        ptr_perm = _array_ptr(self.__perm)
        ptr_iparm = _array_ptr(self.__iparm)
        verb = ctypes.c_int(0)
        err = ctypes.c_long(0)
        ptr_b = _array_ptr(b)
        ptr_x = _array_ptr(x)
        self.__pardiso(self.__handle, ctypes.byref(maxf), ctypes.byref(mnum), \
                    ctypes.byref(self.__mtype), \
                    ctypes.byref(step), ctypes.byref(n), ptr_a, ptr_ia, ptr_ja, \
                    ptr_perm, ctypes.byref(m), ptr_iparm, ctypes.byref(verb), \
                    ptr_b, ptr_x, ctypes.byref(err))
        return err.value

    def diag(self):
        rows = self.__ia.shape[0] - 1
        perm = self.__perm - 1
        f = numpy.zeros((4, rows), dtype=self.__dtype)
        w = numpy.zeros((4, rows), dtype=self.__dtype)
        for i in range(rows):
            j = i%4
            f[j, perm[i]] = 1
        self.solve(f, w, part = 'd')
        w = w[:,perm]
        d = numpy.zeros((2, rows), dtype=self.__dtype)
        i = 0
        while i < rows:
            if i > 0:
                d[1, i - 1] = w[0, i - 1]
            d[0, i] = w[0, i]
            d[1, i] = w[1, i]
            if i == rows - 1:
                break
            i += 1
            d[0, i] = w[1, i]
            d[1, i] = w[2, i]
            if i == rows - 1:
                break
            i += 1
            d[0, i] = w[2, i]
            d[1, i] = w[3, i]
            if i == rows - 1:
                break
            i += 1
            d[0, i] = w[3, i]
            i += 1
        return d

    def inertia(self):
        if self.__real:
            return self.__iparm[22], self.__iparm[21]
        rows = self.__ia.shape[0] - 1
        diags = self.diag()
        diag = numpy.real(diags[0, :])
        offd = diags[1, :]
        nneg = 0
        npos = 0
        i = 0
        while i < rows:
            a = diag[i]
            c = offd[i]
            if c == 0 or i == rows - 1:
                if a < 0:
                    nneg += 1
                elif a > 0:
                    npos += 1
                i += 1
            else:
                b = diag[i + 1]
                det = a*b - abs(c)**2
                s = a + b
                if det == 0:
                    if s < 0:
                        nneg += 1
                    elif s > 0:
                        npos += 1
                elif det < 0:
                    nneg += 1
                    npos += 1
                elif s < 0:
                    nneg += 2
                i += 2
        return nneg, npos

max_stars_repo_name: evgueni-ovtchinnikov/raleigh

max_issues_repo_name: evgueni-ovtchinnikov/raleigh

max_forks_repo_name: evgueni-ovtchinnikov/raleigh


text:
import numpydoc
from distutils.version import LooseVersion as Version
from numpydoc.docscrape import FunctionDoc as _FunctionDoc
from numpydoc.docscrape import NumpyDocString  # noqa

try:
    if Version(numpydoc.__version__) < '0.7':
        raise ImportError('numpydoc mangling requires numpydoc 0.7.')
except AttributeError:
    raise ImportError('numpydoc mangling requires numpydoc 0.7.')

class FunctionDoc(_FunctionDoc):

    def _str_index(self):
        # This should not be necessary after
        # See https://github.com/numpy/numpydoc/pull/187/files
        idx = self['index']
        out = []
        output_index = False
        default_index = idx.get('default', '')
        if default_index:
            output_index = True
        out += ['.. index:: %s' % default_index]
        for section, references in idx.items():
            if section == 'default':
                continue
            output_index = True
            out += ['   :%s: %s' % (section, ', '.join(references))]
        if output_index:
            return out
        else:
            return ''

    def _str_see_also(self, func_role):
        # See https://github.com/numpy/numpydoc/pull/189
        if not self['See Also']:
            return []
        out = []
        out += self._str_header("See Also")
        last_had_desc = True
        for func, desc, role in self['See Also']:
            if role:
                link = ':%s:`%s`' % (role, func)
            elif func_role:
                link = ':%s:`%s`' % (func_role, func)
            else:
                link = "%s" % func
            if desc or last_had_desc:
                out += ['']
                out += [link]
            else:
                out[-1] += ", %s" % link
            if desc:
                out += self._str_indent([' '.join(desc)])
                last_had_desc = True
            else:
                last_had_desc = False
        out += ['']
        return out

    def __str__(self, func_role=''):
        out = []
        # See https://github.com/numpy/numpydoc/issues/190
        # out += self._str_signature()
        out += self._str_summary()
        out += self._str_extended_summary()
        for param_list in ('Parameters', 'Returns', 'Yields',
                           'Other Parameters', 'Raises', 'Warns'):
            out += self._str_param_list(param_list)
        out += self._str_section('Warnings')
        out += self._str_see_also(func_role)
        for s in ('Notes', 'References', 'Examples'):
            out += self._str_section(s)
        for param_list in ('Attributes', 'Methods'):
            out += self._str_param_list(param_list)
        out += self._str_index()
        return '\n'.join(out)

    def __add__(self, other):
        for section in self.sections:
            if section in ['Signature', 'Summary']:
                # Summary should be
                continue
            elif section == 'index':
                idx = {**other['index']}
                idx.update(self['index'])
                self['index'] = idx
            elif section == 'Extended Summary':
                self['Extended Summary'] += other['Summary']
                self['Extended Summary'] += other['Extended Summary']
            else:
                self[section] += other[section]
        return self

"""

The __str__, _str_index, and __str_see_also methods are slightly modified
versions of the methods that belong to the numpydoc developmenet team.

This copyright applies to those 3 methods listed above.

Copyright (C) 2008 Stefan van der Walt <stefan@mentat.za.net>,
Pauli Virtanen <pav@iki.fi>

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:

 1. Redistributions of source code must retain the above copyright
    notice, this list of conditions and the following disclaimer.
 2. Redistributions in binary form must reproduce the above copyright
    notice, this list of conditions and the following disclaimer in
    the documentation and/or other materials provided with the
    distribution.

THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.

"""

reason: import numpy,from numpy


text:
# coding=UTF-8

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

image = np.array([[[-1., -1., 1., -1., -1.],
                   [-1., 1., -1., 1., -1.],
                   [-1., 1., -1., 1., -1.],
                   [-1., 1., -1., 1., -1.],
                   [-1., -1., 1., -1., -1.]],

                  [[-1., -1., 1., -1., -1.],
                   [-1., 1., 1., -1., -1.],
                   [-1., -1., 1., -1., -1.],
                   [-1., -1., 1., -1., -1.],
                   [-1., -1., 1., -1., -1.]],

                  [[-1., 1., 1., 1., -1.],
                   [1., -1., -1., 1., -1.],
                   [-1., -1., -1., 1., -1.],
                   [-1., -1., 1., -1., -1.],
                   [1., 1., 1., 1., 1.]],

                  [[-1., 1., 1., 1., -1.],
                   [-1., -1., -1., -1., 1.],
                   [-1, -1., 1., 1., 1.],
                   [-1., -1., -1., -1., 1.],
                   [-1., 1., 1., 1., -1.]],

                  [[-1., -1., 1., 1., -1.],
                   [-1., 1., -1., 1., -1.],
                   [1., 1., 1., 1., 1.],
                   [-1., -1., -1., 1., -1.],
                   [-1., -1., -1., 1., -1.]],

                  [[-1., 1., 1., 1., -1.],
                   [-1., 1., -1., -1., -1.],
                   [-1., 1., 1., 1., -1.],
                   [-1., -1., -1., 1., -1],
                   [-1., 1., 1., 1., -1.]]])

image2 = np.array([[[1., -1., -1., -1., -1.],
                    [1., -1., -1., -1., -1.],
                    [1., -1., -1., -1., -1.],
                    [1., -1., -1., -1., -1.],
                    [1., -1., -1., -1., -1.]],

                   [[-1., 1., -1., -1., -1.],
                    [-1., 1., -1., -1., -1.],
                    [-1., 1., -1., -1., -1.],
                    [-1., 1., -1., -1., -1.],
                    [-1., 1., -1., -1., -1.]],

                   [[-1., -1., 1., -1., -1.],
                    [-1., -1., 1., -1., -1.],
                    [-1., -1., 1., -1., -1.],
                    [-1., -1., 1., -1., -1.],
                    [-1., -1., 1., -1., -1.]],

                   [[-1., -1., -1., 1., -1.],
                    [-1., -1., -1., 1., -1.],
                    [-1., -1., -1., 1., -1.],
                    [-1., -1., -1., 1., -1.],
                    [-1., -1., -1., 1., -1.]],

                   [[-1., -1., -1., -1., 1.],
                    [-1., -1., -1., -1., 1.],
                    [-1., -1., -1., -1., 1.],
                    [-1., -1., -1., -1., 1.],
                    [-1., -1., -1., -1., 1.]]])

image3 = np.array([[[1., 1., 1., 1., 1.],
                    [-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.]],

                   [[-1., -1., -1., -1., -1.],
                    [1., 1., 1., 1., 1.],
                    [-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.]],

                   [[-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.],
                    [1., 1., 1., 1., 1.],
                    [-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.]],

                   [[-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.],
                    [1., 1., 1., 1., 1.],
                    [-1., -1., -1., -1., -1.]],

                   [[-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.],
                    [-1., -1., -1., -1., -1.],
                    [1., 1., 1., 1., 1.]]])

if __name__ == "__main__":
    num_pattern = 5

    fig, axs = plt.subplots(nrows=1, ncols=num_pattern, figsize=(3 * num_pattern, 5))

    for i in range(num_pattern):
        plot = (image3[i] + 1.) / 2.

        axs[i].imshow(plot, cmap=cm.Greys_r)

    fig.savefig("image3.png")

max_stars_repo_name: saragai/hopfield

max_issues_repo_name: saragai/hopfield

max_forks_repo_name: saragai/hopfield


text:
#!/usr/bin/env python
# coding: utf8
#
# Copyright (c) 2021 Centre National d'Etudes Spatiales (CNES).
#
# This file is part of demcompare
# (see https://github.com/CNES/demcompare).
#
# 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.
#
"""
Nuth and Kaab universal co-registration
(Correcting elevation data for glacier change detection 2011).

Based on the work of geoutils project
https://github.com/GeoUtils/geoutils/blob/master/geoutils/dem_coregistration.py
Authors : Amaury Dehecq, Andrew Tedstone
Date : June 2015
License : MIT
"""

# Standard imports
import argparse
import os

# Third party imports
import matplotlib.pyplot as pl
import numpy as np
from scipy.interpolate import RectBivariateSpline
from scipy.optimize import leastsq

# DEMcompare imports
from .img_tools import load_dems, read_img_from_array, save_tif

def grad2d(dem):
    """

    :param dem:
    :return: slope (fast forward style) and aspect
    """
    grad2, grad1 = np.gradient(dem)  # in Python, x and y axis reversed

    slope = np.sqrt(grad1 ** 2 + grad2 ** 2)
    aspect = np.arctan2(-grad1, grad2)  # aspect=0 when slope facing north
    aspect = aspect + np.pi

    return slope, aspect

def nuth_kaab_single_iter(dh, slope, aspect, plot_file=None):
    """
    Compute the horizontal shift between 2 DEMs
    using the method presented in Nuth & Kaab 2011

    Inputs :
    - dh : array, elevation difference master_dem - slave_dem
    - slope/aspect : array, slope and aspect for the same locations as the dh
    - plot_file : file to where store plot.
        Set to None if plot is to be printed. Set to False for no plot at all.
    Returns :
    - east, north, c : f, estimated easting and northing of the shift,
        c is not used here but is related to the vertical shift
    """
    # The aim is to compute dh / tan(alpha) as a function of the aspect
    # -> hence we are going to be slice the aspect to average a value
    #     for dh / tan(alpha) on those sliced areas
    # -> then we are going to fit the values by the model a.cos(b-aspect)+c
    #    - a will be the magnitude of the horizontal shift
    #    - b will be its orientation
    #    - c will be a vertical mean shift

    # function to be correlated with terrain aspect
    # NB : target = dh / tan(alpha) (see Fig. 2 of Nuth & Kaab 2011)
    # Explicitely ignore divide by zero warning,
    #   as they will be processed as nan later.
    with np.errstate(divide="ignore", invalid="ignore"):
        target = dh / slope
    target = target[np.isfinite(dh)]
    aspect = aspect[np.isfinite(dh)]

    # compute median value of target for different aspect slices
    slice_bounds = np.arange(0, 2 * np.pi, np.pi / 36)
    mean = np.zeros([len(slice_bounds)])
    x_s = np.zeros([len(slice_bounds)])
    j = 0
    for i in slice_bounds:
        target_slice = target[
            (i < aspect) & (aspect < i + np.pi / 36)
        ]  # select target in the slice
        target_slice = target_slice[
            (target_slice < 200) & (target_slice > -200)
        ]  # avoid target>200 and target<-200
        mean[j] = np.median(target_slice)
        x_s[j] = i
        j = j + 1

    # function to fit according to Nuth & Kaab
    x = aspect.ravel()
    y = target.ravel()

    # remove non-finite values
    xf = x[(np.isfinite(x)) & (np.isfinite(y))]
    yf = y[(np.isfinite(x)) & (np.isfinite(y))]

    # remove outliers
    p1 = np.percentile(yf, 1)
    p99 = np.percentile(yf, 99)
    xf = xf[(p1 <= yf) & (yf <= p99)]
    yf = yf[(p1 <= yf) & (yf <= p99)]

    # set the first guess
    p0 = (3 * np.std(yf) / (2 ** 0.5), 0, np.mean(yf))

    # least square fit : peval defines the model chosen
    def peval(x, p):
        return p[0] * np.cos(p[1] - x) + p[2]

    def residuals(p, y, x):
        err = peval(x, p) - y
        return err

    # we run the least square fit
    # by minimizing the "distance" between y and peval (see residuals())
    plsq = leastsq(residuals, p0, args=(mean, x_s), full_output=1)
    yfit = peval(x_s, plsq[0])

    # plotting results
    if plot_file is not False:
        pl.figure(1, figsize=(7.0, 8.0))
        pl.plot(x_s, mean, "b.")
        pl.plot(x_s, yfit, "k-")
        # ax.set_ylim([np.min(mean),])
        pl.xlabel("Terrain aspect (rad)")
        pl.ylabel(r"dh/tan($\alpha$)")
        if plot_file:
            pl.savefig(plot_file, dpi=100, bbox_inches="tight")
        else:
            pl.show()
        pl.close()

    a, b, c = plsq[0]
    east = a * np.sin(b)  # with b=0 when north (origin=y-axis)
    north = a * np.cos(b)

    return east, north, c

def nuth_kaab_lib(dsm_dataset, ref_dataset, outdir_plot=None, nb_iters=6):
    """
    This is the lib api of nuth and kaab universal coregistration.
    It offers quite the same services as the classic main api
    but uses Xarray as input instead of raster files.
    This allows the user to pre-process the data
    and / or avoid unecessary reload of the data.

    Output coregister DSM might be saved.
    Plots might be saved as well (and then not printed) if outputPlot is set.

    NB : it is assumed that both dem3Draster have np.nan values inside
    the '.r' field as masked values

    :param dsm_dataset: xarray Dataset
    :param ref_dataset: xarray Dataset
    :param outdir_plot: path to output Plot directory
        (plots are printed if set to None)
    :param nb_iters: Nuth and Kaab method iterations number default: 6
    :return: x and y shifts (as 'dsm_dataset + (x,y) = ref_dataset')
    """

    # Set target dsm grid for interpolation purpose
    xgrid = np.arange(dsm_dataset["im"].data.shape[1])
    ygrid = np.arange(dsm_dataset["im"].data.shape[0])
    # TODO : not used, to clean
    # x_pixels, y_pixels = np.meshgrid(xgrid, ygrid)
    # trans = dsm_dataset["trans"].data
    # Xgeo = trans[0] + (x_pixels + 0.5) * trans[1]
    #       + (y_pixels + 0.5) * trans[2]
    # Ygeo = trans[3] + (x_pixels + 0.5) * trans[4]
    #       + (y_pixels + 0.5) * trans[5]

    initial_dh = ref_dataset["im"].data - dsm_dataset["im"].data
    coreg_ref = ref_dataset["im"].data

    # Display
    median = np.median(initial_dh[np.isfinite(initial_dh)])
    nmad_old = 1.4826 * np.median(
        np.abs(initial_dh[np.isfinite(initial_dh)] - median)
    )
    maxval = 3 * nmad_old
    pl.figure(1, figsize=(7.0, 8.0))
    pl.imshow(initial_dh, vmin=-maxval, vmax=maxval)
    color_bar = pl.colorbar()
    color_bar.set_label("Elevation difference (m)")
    if outdir_plot:
        pl.savefig(
            os.path.join(outdir_plot, "ElevationDiff_BeforeCoreg.png"),
            dpi=100,
            bbox_inches="tight",
        )
    else:
        pl.show()
    pl.close()

    # Since later interpolations will consider nodata values as normal values,
    # we need to keep track of nodata values
    # to get rid of them when the time comes
    nan_maskval = np.isnan(coreg_ref)
    dsm_from_filled = np.where(nan_maskval, -9999, coreg_ref)

    # Create spline function for interpolation
    spline_1 = RectBivariateSpline(ygrid, xgrid, dsm_from_filled, kx=1, ky=1)
    spline_2 = RectBivariateSpline(ygrid, xgrid, nan_maskval, kx=1, ky=1)
    xoff, yoff, zoff = 0, 0, 0

    print("Nuth & Kaab iterations: {}".format(nb_iters))
    coreg_dsm = dsm_dataset["im"].data
    for i in range(nb_iters):
        # remove bias
        coreg_ref -= median

        # Elevation difference
        dh = coreg_dsm - coreg_ref
        slope, aspect = grad2d(coreg_dsm)

        # compute offset
        if outdir_plot:
            plotfile = os.path.join(
                outdir_plot, "nuth_kaab_iter#{}.png".format(i)
            )
        else:
            plotfile = None
        east, north, z = nuth_kaab_single_iter(
            dh, slope, aspect, plot_file=plotfile
        )
        print(
            "# {} - Offset in pixels : "
            "({:.2f},{:.2f}), -bias : ({:.2f})".format(i + 1, east, north, z)
        )

        xoff += east
        yoff += north
        zoff += z

        # resample slave DEM in the new grid
        # spline 1 : positive y shift moves south
        znew = spline_1(ygrid - yoff, xgrid + xoff)
        nanval_new = spline_2(ygrid - yoff, xgrid + xoff)

        # we created nan_maskval so that non nan values are set to 0.
        # interpolation "creates" values
        # and the one not affected by nan are the one still equal to 0.
        # hence, all other values must be considered invalid ones
        znew[nanval_new != 0] = np.nan

        # update DEM
        if xoff >= 0:
            coreg_ref = znew[:, 0 : znew.shape[1] - int(np.ceil(xoff))]
            coreg_dsm = dsm_dataset["im"].data[
                :, 0 : dsm_dataset["im"].data.shape[1] - int(np.ceil(xoff))
            ]
        else:
            coreg_ref = znew[:, int(np.floor(-xoff)) : znew.shape[1]]
            coreg_dsm = dsm_dataset["im"].data[
                :, int(np.floor(-xoff)) : dsm_dataset["im"].data.shape[1]
            ]
        if -yoff >= 0:
            coreg_ref = coreg_ref[0 : znew.shape[0] - int(np.ceil(-yoff)), :]
            coreg_dsm = coreg_dsm[
                0 : dsm_dataset["im"].data.shape[0] - int(np.ceil(-yoff)), :
            ]
        else:
            coreg_ref = coreg_ref[int(np.floor(yoff)) : znew.shape[0], :]
            coreg_dsm = coreg_dsm[
                int(np.floor(yoff)) : dsm_dataset["im"].data.shape[0], :
            ]

        # print some statistics
        diff = coreg_ref - coreg_dsm
        diff = diff[np.isfinite(diff)]
        nmad_new = 1.4826 * np.median(np.abs(diff - np.median(diff)))
        median = np.median(diff)

        print(
            (
                "\tMedian : {0:.2f}, NMAD = {1:.2f}, Gain : {2:.2f}".format(
                    median, nmad_new, (nmad_new - nmad_old) / nmad_old * 100
                )
            )
        )
        nmad_old = nmad_new

    print(
        "Nuth & Kaab Final Offset in pixels (east, north):"
        "({:.2f},{:.2f})\n".format(xoff, yoff)
    )

    #
    # Get geo raster from coreg_ref array
    #
    coreg_dsm_dataset = read_img_from_array(
        coreg_dsm, from_dataset=dsm_dataset, no_data=-32768
    )
    coreg_ref_dataset = read_img_from_array(
        coreg_ref, from_dataset=dsm_dataset, no_data=-32768
    )
    initial_dh_dataset = read_img_from_array(
        initial_dh, from_dataset=dsm_dataset, no_data=-32768
    )
    final_dh_dataset = read_img_from_array(
        coreg_ref - coreg_dsm, from_dataset=dsm_dataset, no_data=-32768
    )

    return (
        xoff,
        yoff,
        zoff,
        coreg_dsm_dataset,
        coreg_ref_dataset,
        initial_dh_dataset,
        final_dh_dataset,
    )

def run(
    dsm_to,
    dsm_from,
    outfile=None,
    nb_iters=6,
    outdir_plot=None,
    nan_dsm_to=None,
    nan_dsm_from=None,
    save_diff=False,
):
    """
    Coregister dsm_from to dsm_to using Nuth & Kaab (2011) algorithm.

    Output coregister DSM might be saved.
    Plots might be saved as well (and then not printed) if outputPlot is set.

    If nan_dsm and/or nan_ref are not set,
    no data values are read from dsms metadata

    Both input dsm are projected on the same grid :
     - with dsm_to resolution
     - on the biggest common footprint

    TODO: not used in code, Refactor with nuth_kaab_lib function.

    :param dsm_to: path to dsm to coregister to
    :param dsm_from: path to dsm to coregister from
    :param outfile: path to dsm_from after coregistration to dsm_to
    :param nb_iters: Nuth and Kaab method iterations (default(6))
    :param outdir_plot: path to output Plot directory
        (plots are printed if set to None)
    :param nan_dsm_to:
    :param nan_dsm_from:
    :param save_diff: save ./initial_dh.tiff and ./final_dh.tiff
        with dsms diff before and after coregistration
    :return: x and y shifts (as 'dsm_from + (x,y) = dsm_to')
    """

    #
    # Create datasets
    #
    reproj_dem, reproj_ref = load_dems(
        dsm_to,
        dsm_from,
        ref_nodata=nan_dsm_from,
        dem_nodata=nan_dsm_to,
        load_data=False,
    )

    #
    # Coregister and compute diff
    #
    (
        xoff,
        yoff,
        zoff,
        coreg_dsm_dataset,
        coreg_ref_dataset,
        init_dh_dataset,
        final_dh_dataset,
    ) = nuth_kaab_lib(
        reproj_dem, reproj_ref, nb_iters=nb_iters, outdir_plot=outdir_plot
    )

    #
    # Save coreg dem
    #
    if outfile:
        save_tif(coreg_dsm_dataset, "./coreg_dsm.tif")
        save_tif(coreg_ref_dataset, "./coreg_ref.tif")

    #
    # Save diffs
    #
    if save_diff:
        save_tif(init_dh_dataset, "./initial_dh.tif")
        save_tif(final_dh_dataset, "./final_dh.tif")

    return xoff, yoff, zoff

def get_parser():
    """
    Parser of nuth kaab independent main
    TODO: To clean with main. Keep independent main ?
    """
    parser = argparse.ArgumentParser(
        os.path.basename(__file__),
        description="Universal co-registration method "
        "presented in Nuth & Kaab 2011."
        "NB : 1) It is supposed that both dsms share common reference"
        "            (whether it is geoid or ellipsoid)."
        "     2) DSMs must be georefenced.",
    )

    parser.add_argument("dsm_to", type=str, help="master dsm")
    parser.add_argument(
        "dsm_from", type=str, help="slave dsm you wish to coregister to dsm_to"
    )
    parser.add_argument(
        "-outfile", action="store_true", help="saves output coregistered DSM"
    )
    parser.add_argument(
        "-nb_iters",
        dest="nb_iters",
        type=int,
        default=6,
        help="number of iterations",
    )
    parser.add_argument(
        "-dirplot",
        dest="plot",
        type=str,
        default=None,
        help="path to output plot directory. "
        "Plots are printed if set to None (default)",
    )
    parser.add_argument(
        "-nodata1",
        dest="nodata1",
        type=str,
        default=None,
        help="no data value for DSM to compare "
        "(default value is read in metadata)",
    )
    parser.add_argument(
        "-nodata2",
        dest="nodata2",
        type=str,
        default=None,
        help="no data value for Reference DSM "
        "(default value is read in metadata)",
    )
    parser.add_argument(
        "-save_diff",
        action="store_true",
        help="store on file system a ./initial_dh.tiff and a ./final_dh.tiff "
        "with dsms differences before and after coregistration",
    )

    return parser

def main():
    """
    Main from Nuth Kaab API lib
    """
    parser = get_parser()
    args = parser.parse_args()

    run(
        args.dsm_to,
        args.dsm_from,
        outfile=args.outfile,
        nb_iters=args.nb_iters,
        outdir_plot=args.plot,
        nan_dsm_to=args.nodata1,
        nan_dsm_from=args.nodata2,
        save_diff=args.save_diff,
    )

if __name__ == "__main__":
    main()

max_stars_repo_name: fchapoton/demcompare

max_issues_repo_name: fchapoton/demcompare

max_forks_repo_name: fchapoton/demcompare

reason: import numpy,from scipy


text:
import pandas as pd
import numpy as np
import plotly.graph_objects as go
from optbinning import OptimalBinning
import math
from plotly.subplots import make_subplots
import plotly.figure_factory as ff
import plotly.express as px

def DistributionPlot(Df, PlotVar):
    '''
    Plots the distribution of a given variable in a dataframe
    '''
    Labels = [i for i in range(0, 100, 10)]

    BinSize = Df[PlotVar].describe().loc["std"] / 20

    fig = ff.create_distplot(
            hist_data = [Df[Df["Target"] == 0][PlotVar].values, Df[Df["Target"] == 1][PlotVar].values]
            , group_labels = [0, 1]
            , bin_size=BinSize
            , show_hist=True)

    fig.update_xaxes(
    zeroline = True
    , showgrid = True
    , title=PlotVar)

    fig.update_yaxes(
        zeroline=True
        , showgrid=True
        , title="Distribution")

    fig.update_layout(
        title = dict(text=str(PlotVar) + " Distribution"
                     , font=dict(color="Black", size=20))
        , font = dict(color="Black", size=10)
        , height = 700
        , width = 1100
        , legend_title='Target')

    fig.show(renderer='png', height=700, width=1100)

def Scatter(Df, PlotVar, Hue, Y, Title):

    '''
    Produces a plot of data pulled from specified dataframe split by a certain binary population
    PlotVars defines the independent variable
    Hue defines the population for which to split the plots
    Y is the dependent variable
    Title is title of the plot
    '''

    fig = go.Figure()

    fig.add_trace(
        go.Scatter(
            x = Df[Df[Hue] == 1][PlotVar]
            , y=Df[Df[Hue] == 1][Y]
            , legendgroup=Hue + " = 1"
            , name=Hue + " = 1"
            , mode='markers'
            , line=dict(color='red')
            , marker=dict(size=10, opacity=0.1)
            , showlegend= True))

    fig.add_trace(
        go.Scatter(
            x = Df[Df[Hue] == 0][PlotVar]
            , y=Df[Df[Hue] == 0][Y]
            , legendgroup=Hue + " = 0"
            , name=Hue + " = 0"
            , mode='markers'
            , line=dict(color='blue')
            , marker=dict(size=10, opacity=0.1)
            , showlegend= True))

    fig.update_xaxes(
        zeroline = True
        , showgrid = True
        , title = PlotVar
    )

    fig.update_yaxes(
        zeroline=True
        , showgrid=True
        , title = Y)

    fig.update_layout(
        title = dict(text=Title, font=dict(size=17)))

    fig.update_annotations(
        font = dict(size=14))

    fig.show(renderer="png", height=600, width=1000)

def Distribution(Df, Target, Variable):
    Graph = pd.pivot_table(Df, index=Variable, columns=Target, values="Track", aggfunc=len)
    Graph1 = pd.pivot_table(Df, index=Variable, values=Target, aggfunc="mean").sort_values(by="Target", ascending=False)

    Graph = Graph.reindex(Graph1.index)

    fig = make_subplots(specs=[[{"secondary_y": True}]])

    fig.add_trace(
        go.Scatter(
            y=Graph1[Target]*100
            , x=[Col.title() if type(Col) == "str" else Col for Col in Graph1.index.values]
            , name=Target
            , mode="lines"
            , showlegend= True)
        , secondary_y = True)

    fig.add_trace(
        go.Bar(
            y=Graph[0]
            , x=[Col.title() if type(Col) == "str" else Col for Col in Graph1.index.values]
            , name="Not"+str(Target)
            , showlegend= True)
        , secondary_y = False)

    fig.add_trace(
        go.Bar(
            y=Graph[1]
            , x=[Col.title() if type(Col) == "str" else Col for Col in Graph1.index.values]
            , name=Target
            , showlegend= True)
        , secondary_y = False)

    fig.update_xaxes(
        zeroline = True
        , showgrid = True
        , title = Variable
        , type="category")

    fig.update_yaxes(
        zeroline=True
        , showgrid=True
        , title="Frequency"
        , secondary_y = False)

    fig.update_yaxes(
        zeroline=True
        , showgrid=False
        , title=Target
        , ticksuffix="%"
        , range=[0, 100]
        , secondary_y = True)

    fig.update_layout(
        title = dict(text= str(Variable) +" Distribution vs. " + str(Target), font=dict(color="Black", size=20))
        , font = dict(color="Black", size=10)
        , height = 600
        , width = 900
        , barmode='stack')

    fig.update_annotations(
    font = dict(color="Black", size=14))

    fig.show(renderer="png", width=900, height=600)

class VariableBinning():
    '''
    Class to bin a variable according to a selection of metrics.
    :Attribute BinPlot: Plots Distribution vs. Event rate for a DataFrame with class Count columns and Event rate column.
    :Attribute BinVariable: Fits OptBinning algorithm and prints summary plot along with BinPlot (For visualising)
    :Attribute Transform: Fits and transforms variable, returns transformed series.

    '''

    def __init__(self, Df, Variable, Target, DType = "numerical"):
        self.Temp = Df.copy()[[Variable, Target, "Track"]]
        self.Variable = Variable
        self.Target = Target
        self.Mod = None
        self.DType = DType

    def BinPlot(self, Graph):
        '''
        :Param Graph: Dataframe containing Class count columns and event rate column
        '''

        fig = make_subplots(specs=[[{"secondary_y": True}]])

        fig.add_trace(
            go.Scatter(
                y=Graph["EventRate"]*100
                , x=[Col.title() if type(Col) == "str" else Col for Col in Graph.index.values]
                , name=self.Target
                , mode="lines+markers"
                , showlegend= True)
            , secondary_y = True)

        for Col in [str(self.Target)+" == 0", str(self.Target)+" == 1"]:
            fig.add_trace(
                go.Bar(
                    y=Graph[Col]
                    , x=[Col.title() if type(Col) == "str" else Col for Col in Graph.index.values]
                    , name=Col
                    , showlegend= True)
                , secondary_y = False)

        fig.update_xaxes(
            zeroline = True
            , showgrid = True
            , title = self.Variable
            , type='category' if self.DType == "categorical" else "-")

        fig.update_yaxes(
            zeroline=True
            , showgrid=True
            , title="Count"
            , secondary_y = False)

        fig.update_yaxes(
            zeroline=True
            , showgrid=False
            , title=self.Target
            , ticksuffix="%"
            , range=[0, 100]
            , secondary_y = True)

        fig.update_layout(
            title = dict(text= str(self.Variable) +" Distribution vs. " + str(self.Target), font=dict(color="Black", size=20))
            , font = dict(color="Black", size=10)
            , height = 600
            , width = 900
            , barmode='stack')

        fig.show()

    def BinVariable(self, Trend = "auto_asc_desc", Method = "bins", ShowTable = False):
        '''
        :Param Trend: Default = "auto_asc_desc", sets the assumed trend for binning
        :Param Method: Default = "bins", sets the desired transformation method.
        '''

        self.Mod = OptimalBinning(name=self.Variable, dtype=self.DType, solver="cp", max_n_prebins=100, monotonic_trend=Trend,
                                min_prebin_size=0.01, time_limit=30)

        #Fit and record
        self.Mod.fit(self.Temp[self.Variable], self.Temp[self.Target])

        BinningTable = self.Mod.binning_table
        Table = BinningTable.build()
        BinsValues = self.Mod.splits

        if ShowTable == True:
            #Print Status and Summary
            print(self.Mod.status)
            print(BinningTable.analysis())

        self.Temp["Transformed"] = self.Mod.transform(self.Temp[self.Variable], metric = Method)

        if ((self.DType == "numerical") & (Method != "woe")):
            self.Temp["Transformed"] = self.Temp["Transformed"].apply(lambda s: tuple(float(x) for x in s.replace('[', '').replace(')', '').split(',')))
            self.Temp["Transformed"] = self.Temp["Transformed"].apply(lambda x: x[0] + 1 if math.isinf(x[1]) else x[1])

        Graph = pd.pivot_table(self.Temp, index="Transformed", columns=self.Target, values = "Track", aggfunc="count")
        Graph = Graph.rename({0: str(self.Target)+" == 0", 1: str(self.Target)+" == 1"}, axis=1)

        Graph1 = pd.pivot_table(self.Temp, index="Transformed", values=self.Target, aggfunc="mean")

        Graph["EventRate"] = Graph1[self.Target]

        self.BinPlot(Graph)

    def Transform(self, Df = None, Method = 'woe'):
        '''
        :Param Trend: Default = "auto_asc_desc", sets the assumed trend for binning
        :Param Method: Default = "bins", sets the desired transformation method.
        '''

        if Df is not None:
            DataFrame = Df.copy()
            DataFrame = DataFrame[[self.Variable]]

            DataFrame["Transformed"] = self.Mod.transform(DataFrame[self.Variable], metric = Method)

            return DataFrame["Transformed"]

        else:
            self.Temp["Transformed"] = self.Mod.transform(self.Temp[self.Variable], metric = Method)

            Graph = pd.pivot_table(self.Temp, index="Transformed", columns=self.Target, values = "Track", aggfunc="count")
            Graph = Graph.rename({0: str(self.Target)+" == 0", 1: str(self.Target)+" == 1"}, axis=1)

            Graph1 = pd.pivot_table(self.Temp, index="Transformed", values=self.Target, aggfunc="mean").sort_values(by=self.Target, ascending=False)

            Graph = Graph.reindex(Graph1.index)
            Graph["EventRate"] = Graph1[self.Target]

            self.BinPlot(Graph)

            return self.Temp["Transformed"]

def Correlation(Df, PlotVars, Title):

    '''
    Provides a correlation matrix heatmap for data pulled from a specified dataframe
    PlotVars define all features in question
    Title is title of the plot
    '''

    Correlations = Df[PlotVars].corr()

    fig = go.Figure()

    fig.add_trace(
        go.Heatmap(
            z=Correlations
            , x=Correlations.index
            , y=Correlations.index
            , zmax=1
            , zmin=-1
            , hoverongaps = False
            , colorscale=[(0, "blue"), (0.5, "white"), (1, "red")]))

    fig.update_layout(
        title = dict(text=Title, font=dict(color="Black", size=20))
        , font = dict(color="Black", size=10)
        , height = 1000
        , width = 1000
        , legend_title='Period')

    fig.update_annotations(
    font = dict(color="Black", size=14))

    fig.show(renderer="png", height=900, width=900)

def BarPlot(DataFrame, Title):

    fig = go.Figure()

    for Column in DataFrame.columns.values:

        fig.add_trace(
            go.Bar(
                y=DataFrame[Column]
                , x=DataFrame.columns.values
                , name=str(Column)
                , showlegend= True))

    fig.update_xaxes(
        zeroline = True
        , showgrid = True
        , title = "Features"
        , showticklabels=False)

    fig.update_yaxes(
        zeroline=True
        , showgrid=True
        , title="Importance")

    fig.update_layout(
        title = dict(text= Title, font=dict(color="Black", size=20))
        , font = dict(color="Black", size=10)
        , height = 600
        , width = 900)

    fig.update_annotations(
    font = dict(color="Black", size=14))

    fig.show(renderer="png", width=900, height=600)

def InformationValue(Df, Variable, Target):
    '''
    Computes the information for a given dataframe feature w.r.t a target/dependent variable.
    '''
    Pivot = pd.pivot_table(Df, index=Variable, values="Track", columns=Target, aggfunc="count").reset_index()
    Pivot = Pivot.rename({0:"Flops", 1:"Hits"}, axis=1)
    Pivot["Flops"] = Pivot["Flops"] / Pivot["Flops"].sum()
    Pivot["Hits"] = Pivot["Hits"] / Pivot["Hits"].sum()

    Pivot["IV"] = Pivot["Flops"] - Pivot["Hits"]
    Pivot["IV"] = Pivot["IV"]*Pivot[Variable]

    return Pivot["IV"].sum()


text:
import numpy as np

def infer(signal,logging,model):
    #START
    result = "cough"
    #END
    #return result as either "cough" or "non_cough"
    return result

max_stars_repo_name: MohamedAfham/cough-detector-app

max_issues_repo_name: MohamedAfham/cough-detector-app

max_forks_repo_name: MohamedAfham/cough-detector-app


text:
#!/usr/bin/env python
# coding: utf-8

# # COVID-19 Growth By State (US)
# > Growth of COVID-19 for the US by State.
#
# - comments: true
# - author: Avy Faingezicht
# - categories: [growth, US, states]
# - image: images/covid-growth-states.png
# - permalink: /growth-us-states/

# In[1]:

#hide
get_ipython().run_line_magic('matplotlib', 'inline')
import math
import requests
import pandas as pd
import numpy as np
from cycler import cycler
import matplotlib.pyplot as plt
from matplotlib.ticker import ScalarFormatter

SMALL_SIZE = 16
MEDIUM_SIZE = 18
BIGGER_SIZE = 20

plt.rc('font', size=SMALL_SIZE)          # controls default text sizes
plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title
plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels
plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize
plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title

states_url = "https://covidtracking.com/api/states/daily"
us_url = "https://covidtracking.com/api/us/daily"
case_threshold = 100

colors = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple',
          'tab:brown', 'tab:pink', 'tab:gray', 'tab:olive', 'tab:cyan']
default_cycler = cycler(linestyle=['-', '--', ':', '-.']) * cycler(color=colors)

r = requests.get(states_url)
states_df = pd.DataFrame(r.json())
states_df['date'] = pd.to_datetime(states_df.date, format="%Y%m%d")
states_df = states_df[['date', 'state', 'positive', 'death']].sort_values('date')
cols = {}
for state in states_df.state.unique():
    cases = states_df[(states_df.state == state) & (states_df.positive > case_threshold)]
    cases = cases.reset_index().positive.reset_index(drop=True)
    if len(cases) > 1:
        cols[state] = cases

r = requests.get(us_url)
us_df = pd.DataFrame(r.json())
us_df['date'] = pd.to_datetime(us_df.date, format="%Y%m%d")
us_df = us_df[['date', 'positive', 'death']].sort_values('date')
cols['US'] = us_df.positive.reset_index(drop=True)

# In[2]:

#hide
populations = {
    # Taken from US Census Bureau NST-EST2019-01: Table 1. Annual Estimates of the Resident Population
    # for the United States, Regions, States, and Puerto Rico: April 1, 2010 to July 1, 2019
    # https://www.census.gov/newsroom/press-kits/2019/national-state-estimates.html
    "AL": 4903185,
    "AK": 731545,
    "AZ": 7278717,
    "AR": 3017804,
    "CA": 39512223,
    "CO": 5758736,
    "CT": 3565287,
    "DE": 973764,
    "DC": 705749,
    "FL": 21477737,
    "GA": 10617423,
    "HI": 1415872,
    "ID": 1787065,
    "IL": 12671821,
    "IN": 6732219,
    "IA": 3155070,
    "KS": 2913314,
    "KY": 4467673,
    "LA": 4648794,
    "ME": 1344212,
    "MD": 6045680,
    "MA": 6892503,
    "MI": 9986857,
    "MN": 5639632,
    "MS": 2976149,
    "MO": 6137428,
    "MT": 1068778,
    "NE": 1934408,
    "NV": 3080156,
    "NH": 1359711,
    "NJ": 8882190,
    "NM": 2096829,
    "NY": 19453561,
    "NC": 10488084,
    "ND": 762062,
    "OH": 11689100,
    "OK": 3956971,
    "OR": 4217737,
    "PA": 12801989,
    "RI": 1059361,
    "SC": 5148714,
    "SD": 884659,
    "TN": 6829174,
    "TX": 28995881,
    "UT": 3205958,
    "VT": 623989,
    "VA": 8535519,
    "WA": 7614893,
    "WV": 1792147,
    "WI": 5822434,
    "WY": 578759,
    "PR": 3193694,
    "GU": 165_718,
    "VI": 104_914,
    "AS": 55_641,
    "MP": 55_194,
}

# In[3]:

#collapse-hide
fig = plt.figure(figsize=(16, 9))
ax = plt.axes()
ax.set_prop_cycle(default_cycler)
pd.concat(cols, axis=1).plot(ax=ax, marker='o')

plt.title('COVID19 Growth In The US As A Whole And By State', fontsize=BIGGER_SIZE)
plt.ylabel('Cumulative confirmed cases')
plt.xlabel(f'Number of days since {case_threshold}th case')
plt.annotate(f'Latest data from {states_df.date.max().strftime("%Y-%m-%d")}, varies by state\n'
             'Based on COVID Data Repository by the COVID Tracking Project\n'
             'Chart by Avy Faingezicht, @avyfain',
             (0.07, 0.01), xycoords='figure fraction', fontsize=10);

box = ax.get_position()

ax.set_position([box.x0, box.y0, box.width * 0.9, box.height])
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))

x = np.linspace(0, plt.xlim()[1])
plt.plot(x, 100 * (1.33) ** x, ls='--', color='k', label='33% daily growth', linewidth=1)
plt.annotate('33% Daily Growth',
             (0.7, 0.3), xycoords='figure fraction', fontsize=SMALL_SIZE);

plt.grid(axis='y')
formatter = ScalarFormatter()
formatter.set_scientific(False)
ax.yaxis.set_major_formatter(formatter)
ax.yaxis.set_minor_formatter(formatter)

# In[4]:

#collapse-hide
fig = plt.figure(figsize=(16, 8))
ax = plt.axes()
ax.set_prop_cycle(default_cycler)
pd.concat(cols, axis=1).plot(ax=ax, marker='o', logy=True)

plt.title('COVID19 Growth In US As A Whole And By State\nY-axis is on a log scale', fontsize=BIGGER_SIZE)
plt.ylabel('Cumulative confirmed cases (log scale)')
plt.xlabel(f'Number of days since {case_threshold}th case')
plt.annotate(f'Latest data from {states_df.date.max().strftime("%Y-%m-%d")}, varies by state\n'
             'Based on COVID Data Repository by The COVID Tracking Project\n'
             'Chart by Avy Faingezicht, @avyfain',
             (0.07, .01), xycoords='figure fraction', fontsize=10);

box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.9, box.height])

# Put a legend to the right of the current axis
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.plot(x, 100 * (1.33) ** x, ls='--', color='k', label='33% daily growth')
plt.annotate('33% Daily Growth',
             (0.7, 0.65), xycoords='figure fraction', fontsize=SMALL_SIZE);
plt.grid(axis='y')
ax.yaxis.set_major_formatter(formatter)

# In[5]:

#hide
fig.savefig('../images/covid-growth-states.png')

# In[6]:

#collapse-hide

# re-using the df's from above
death_threshold = 1
cols = {}

states_df['population'] = states_df['state'].map(populations)
states_df['deaths_per_m'] = states_df.death / (states_df.population / 1_000_000)

for state in states_df.state.unique():
    deaths = states_df[(states_df.state == state) & (states_df.death > death_threshold)]
    deaths = deaths.reset_index().deaths_per_m
    if len(deaths) > 1:
        cols[state] = deaths

fig = plt.figure(figsize=(16, 8))
ax = plt.axes()
ax.set_prop_cycle(default_cycler)
pd.concat(cols, axis=1).plot(ax=ax, marker='o')#, logy=True)
plt.grid(axis='y')
plt.title('Death Rate In The US By State', fontsize=BIGGER_SIZE)
plt.ylabel('Cumulative deaths per 1M people (log scale)')
plt.xlabel(f'Days since {death_threshold}st reported death')
plt.annotate(f'Latest data from {states_df.date.max().strftime("%Y-%m-%d")}, varies by state\n'
             'Based on COVID Data Repository by The COVID Tracking Project\n'
             'Chart by Avy Faingezicht, @avyfain',
             (0.07, .01), xycoords='figure fraction', fontsize=10);

box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.9, box.height])

ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
ax.yaxis.set_major_formatter(formatter)

# In[7]:

#hide
fig.savefig('../images/covid-death-rate-states.png')

# In[8]:

#collapse-hide

# re-using the df's from above
death_threshold = 1
cols = {}

for state in states_df.state.unique():
    deaths = states_df[(states_df.state == state) & (states_df.death > death_threshold)]
    deaths = deaths.reset_index().death
    if len(deaths) > 1:
        cols[state] = deaths

fig = plt.figure(figsize=(16, 8))
ax = plt.axes()
ax.set_prop_cycle(default_cycler)

pd.concat(cols, axis=1).plot(ax=ax, marker='o', logy=True)
plt.grid(axis='y')
plt.title('Cumalitive Deaths In The US By State', fontsize=BIGGER_SIZE)
plt.ylabel('Cumulative deaths (log scale)')
plt.xlabel(f'Days since {death_threshold}st reported death')
plt.annotate(f'Latest data from {states_df.date.max().strftime("%Y-%m-%d")}, varies by state\n'
             'Based on COVID Data Repository by The COVID Tracking Project\n'
             'Chart by Avy Faingezicht, @avyfain',
             (0.07, .01), xycoords='figure fraction', fontsize=10);

box = ax.get_position()
ax.set_position([box.x0, box.y0, box.width * 0.9, box.height])

ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
ax.yaxis.set_major_formatter(formatter)

# This visualization was made by [Avy Faingezicht](https://twitter.com/avyfain)[^1].
#
# [^1]:  Data sourced from ["The COVID Tracking Project"](https://covidtracking.com/). Link to [original notebook](https://github.com/avyfain/covid19/blob/master/covid19.ipynb).  Updated hourly by [GitHub Actions](https://github.com/features/actions).


text:
#!/usr/bin/env python3
import argparse
import fromage.utils.volume as vo
import fromage.io.read_file as rf
import numpy as np
import time
import sys

def parse_args():
    # parse the input
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "in_xyz", help="Input .xyz file with the unit cell geometry")
    parser.add_argument(
        "-l","--atom_label", help="Label of atom in the molecule to analyse", default=1, type=int)
    parser.add_argument(
        "-b", "--bonding", help="Type of bonding used for detecting full molecules. The options are dis, cov and vdw", default="dis", type=str)
    parser.add_argument(
        "-T", "--thresh", help="Threshold distance to select a bond. The default pairs are dis:1.8, cov:0.2, vdw:-0.3. To select default, just use a threshold of 999", default=999, type=float)
    parser.add_argument(
        "-bs", "--bonding_string", help="Alternate specification of bonding. Here the threshold and bonding are lumped up in one string like 'cov-0.1' or 12dis'", default="", type=str)
    parser.add_argument("-res", "--resolution",
                        help="The number of voxels per side of the box", default=100, type=int)
    parser.add_argument("-dim", "--dimensions",help="Dimensions of the box in Angstrom. Give x y and z like '30 30 30'",default=[25, 25, 25], type=float, nargs='*')

    user_input = sys.argv[1:]
    args = parser.parse_args(user_input)

    return args
def main(args):

    in_atoms = args.in_xyz
    atoms = rf.mol_from_file(in_atoms)

    out_file = open("volumes", "w")
    prox_grid = vo.CubeGrid()
    vdw_grid = vo.CubeGrid()

    mol = atoms.select(args.atom_label - 1)
    c_x, c_y, c_z = mol.centroid()

    x_dim, y_dim, z_dim = args.dimensions[0], args.dimensions[1], args.dimensions[2]

    prox_grid.grid_from_point(c_x, c_y, c_z, res=args.resolution, box=np.array(
        [[x_dim, 0.0, 0.0], [0.0, y_dim, 0.0], [0.0, 0.0, z_dim]]))
    vdw_grid.grid_from_point(c_x, c_y, c_z, res=100, box=np.array(
        [[x_dim, 0.0, 0.0], [0.0, y_dim, 0.0], [0.0, 0.0, z_dim]]))

    rest = []
    for atom in atoms:
        if atom not in mol:
            rest.append(atom)

    prox_grid.set_grid_coord()
    vdw_grid.set_grid_coord()

    prox_grid.proximity(mol, rest)
    vdw_grid.vdw_vol(mol)

    prox_grid.out_cube("voro.cube", atoms)
    out_file.write("Voronoi volume: " + str(prox_grid.volume()) + "\n")

    vdw_grid.out_cube("vdw.cube", atoms)
    out_file.write("VDW volume: " + str(vdw_grid.volume()) + "\n")

    prox_grid.add_grid(vdw_grid.grid)
    out_file.write("Union volume: " + str(prox_grid.volume()) + "\n")
    prox_grid.out_cube("add.cube", atoms)

    out_file.close()

if __name__ == '__main__':
    start = time.time()
    args = parse_args()
    main(args)
    end = time.time()
    print("\nTotal time: {}s".format(round((end - start), 1)))


text:
import importlib
import time
import numpy as np

from empyric.adapters import *
from empyric.collection.instrument import *

class Phidget1101(Instrument):
    """
    Phidgets 4x TC reader
    Many instrumet methods (setX, getY, etc.) are mapped by the adapter from the Phidgets device class
    """

    name = 'Phidget1101'

    supported_adapters = (
        (Phidget, {}),
    )

    # Available knobs
    knobs = ('type',)

    # Available meters
    meters = ('temperature',)

    def __init__(self, *args, **kwarhs):
        self.device_class = importlib.import_module('Phidget22.Devices.TemperatureSensor').TemperatureSensor
        Instrument.__init__(self, *args, **kwargs)

    @setter
    def set_type(self, type_):

        types = importlib.import_module('Phidget22.ThermocoupleType')

        type_dict = {
            'K': types.ThermocoupleType.THERMOCOUPLE_TYPE_K,
            'J': types.ThermocoupleType.THERMOCOUPLE_TYPE_J,
            'T': types.ThermocoupleType.THERMOCOUPLE_TYPE_T,
            'E': types.ThermocoupleType.THERMOCOUPLE_TYPE_E,
        }

        self.write('ThermocoupleType', type_dict[type_])

    @measurer
    def measure_temperature(self):
        return self.query('Temperature')

max_stars_repo_name: dmerthe/empyric

max_issues_repo_name: dmerthe/empyric

max_forks_repo_name: dmerthe/empyric


text:
# utils.py
"""Utilities for logging, timing, loading, and saving."""
import os
import sys
import h5py
import time
import logging
import numpy as np
import matplotlib.pyplot as plt

try:
    import rom_operator_inference as opinf
except ModuleNotFoundError:
    print("\nrom_operator_inference module not installed",
          "(python3 -m pip install --user -r requirements.txt)\n")
    raise
if opinf.__version__ != "1.2.1":
    raise ModuleNotFoundError("rom-operator-inference version 1.2.1 required "
                              "(python3 -m pip install --user "
                              "-r requirements.txt)")

import config

# Logging / timing tools ======================================================

def reset_logger(trainsize=None):
    """Switch to the log file within the folder for experiments with
    `trainsize` snapshots of training data. If `trainsize` is None,
    switch to the log file in the base folder.
    """
    # Remove all old logging handlers.
    logger = logging.getLogger()
    for hdlr in logger.handlers[:]:
        logger.removeHandler(hdlr)

    # Get the log filename and append a newline.
    log_filename = config.log_path(trainsize)
    with open(log_filename, 'a') as lf:
        lf.write('\n')

    # Get a new logging handler to the log file.
    handler = logging.FileHandler(log_filename, 'a')
    handler.setFormatter(logging.Formatter("%(levelname)s: %(message)s"))
    handler.setLevel(logging.INFO)
    logger.addHandler(handler)
    print(f"Logging to {log_filename}")

    # Log the session header.
    if hasattr(sys.modules["__main__"], "__file__"):
        _front = f"({os.path.basename(sys.modules['__main__'].__file__)})"
        _end = time.strftime('%Y-%m-%d %H:%M:%S')
        _mid = '-' * (79 - len(_front) - len(_end) - 20)
        header = f"NEW SESSION {_front} {_mid} {_end}"
    else:
        header = f"NEW SESSION {time.strftime(' %Y-%m-%d %H:%M:%S'):->61}"
    logging.info(header)

class timed_block:
    """Context manager for timing a block of code and reporting the timing.

    >>> with timed_operation("This is a test"):
    ...     # Code to be timed
    ...     time.sleep(2)
    ...
    This is a test...done in 2.00 s.
    """
    def __init__(self, message):
        self.message = message

    def __enter__(self):
        """Print the message and record the current time."""
        print(f"{self.message}...", end='', flush=True)
        self.start = time.time()

    def __exit__(self, exc_type, exc_value, exc_traceback):
        """Calculate and report the elapsed time."""
        self.end = time.time()
        elapsed = self.end - self.start
        if exc_type:    # Report an exception if present.
            print(f"{exc_type.__name__}: {exc_value}")
            logging.info(self.message.strip())
            logging.error(f"({exc_type.__name__}) {exc_value} "
                          f"(raised after {elapsed:.6f} s)")
        else:           # If no exception, report execution time.
            print(f"done in {elapsed:.2f} s.", flush=True)
            logging.info(f"{self.message.strip()}...done in {elapsed:.6f} s.")

# Data loaders ================================================================

class DataNotFoundError(FileNotFoundError):
    """Exception to be raised when attempting to load a missing data file for
    * GEMS simulation data
    * Scaled training data
    * POD basis
    * Projected training data
    * Trained ROMs
    """
    pass

def _checkexists(filename):
    """Check that the file `filename` exists; if not, raise an exception."""
    if not os.path.isfile(filename):
        raise DataNotFoundError(filename)
    return filename

def load_gems_data(rows=None, cols=None):
    """Load the indicated rows and colums of GEMS simulation data.
    This is a large file, so try to only load what is needed at the moment.

    Parameters
    ----------
    rows : int, slice, or (nrows,) ndarray of integer indices
        Which rows (spatial locations) to extract from the data (default all).
        If an integer, extract the first `rows` rows.

    cols : int or slice
        Which columns (temporal points) to extract from the data (default all).
        If an integer, extract the first `cols` columns.

    Returns
    -------
    gems_data : (nrows,ncols) ndarray
        The indicated rows / columns of the data.

    time_domain : (ncols,) ndarray
        The time (in seconds) associated with each column of extracted data.
    """
    # Locate the data.
    data_path = _checkexists(config.gems_data_path())

    # Ensure rows are loaded in ascending index order (HDF5 requirement).
    if isinstance(rows, (np.ndarray, list)):
        row_order = np.argsort(rows)
        rows = np.array(rows, copy=True)[row_order]
        old_row_order = np.argsort(row_order)
    elif np.isscalar(rows) or rows is None:
        rows = slice(None, rows)
    if np.isscalar(cols) or cols is None:
        cols = slice(None, cols)

    # Extract the data.
    NROWS = config.NUM_GEMSVARS * config.DOF
    with timed_block(f"Loading GEMS simulation data from {data_path}"):
        with h5py.File(data_path, 'r') as hf:
            # Check data shape.
            if hf["data"].shape[0] != NROWS:
                raise RuntimeError(f"data should have exactly {NROWS} rows")
            gems_data = hf["data"][rows,cols]
            time_domain = hf["time"][cols]

    # Restore row ordering if needed.
    if isinstance(rows, np.ndarray):
        gems_data = gems_data[old_row_order,:]

    return gems_data, time_domain

def load_time_domain(nsteps=None):
    """Load the time domain corresponding to the GEMS simulation data.

    Parameters
    ----------
    nsteps : int or None
        How many entries to extract from the time domain (default all).
    """
    # Locate the data.
    data_path = _checkexists(config.gems_data_path())

    # Extract the data.
    with timed_block(f"Loading time domain data from {data_path}"):
        with h5py.File(data_path, 'r') as hf:
            time_domain = hf["time"][:nsteps]

    # Check time spacing.
    if not np.allclose(np.diff(time_domain), config.DT):
        raise ValueError("time domain DT != config.DT")

    # If a larger number of time steps is requested, use np.linspace().
    if np.isscalar(nsteps) and time_domain.size < nsteps:
        t0 = time_domain[0]
        return np.linspace(t0, t0 + nsteps*config.DT, nsteps)

    return time_domain

def load_scaled_data(trainsize):
    """Load scaled snapshot data and the associated scaling factors.

    Parameters
    ----------
    trainsize : int
        Number of snapshots of scaled data to load. See step2a_transform.py.

    Returns
    -------
    Q : (NUM_ROMVARS*DOF,trainsize) ndarray
        Lifted, scaled, shifted data.

    time_domain : (trainsize) ndarray
        Time domain corresponding to the lifted, scaled data.

    qbar : (NUM_ROMVARS*DOF,) ndarray
        Mean snapshot of the scaled training data.

    scales : (NUM_ROMVARS,2) ndarray
        Factors used to scale the variables.
    """
    # Locate the data.
    data_path = _checkexists(config.scaled_data_path(trainsize))

    # Extract the data.
    with timed_block(f"Loading lifted, scaled snapshot data from {data_path}"):
        with h5py.File(data_path, 'r') as hf:
            # Check data shapes.
            if hf["data"].shape != (config.NUM_ROMVARS*config.DOF, trainsize):
                raise RuntimeError("data set 'data' has incorrect shape")
            if hf["time"].shape != (trainsize,):
                raise RuntimeError("data set 'time' has incorrect shape")
            if "mean" in hf:
                if hf["mean"].shape != (hf["data"].shape[0],):
                    raise RuntimeError("data set 'mean' has incorrect shape")
                mean = hf["mean"][:]
            else:
                mean = np.zeros(hf["data"].shape[0])
            if hf["scales"].shape != (config.NUM_ROMVARS, 2):
                raise RuntimeError("data set 'scales' has incorrect shape")

            # Load and return the data.
            return (hf["data"][:,:], hf["time"][:],
                    mean, hf["scales"][:,:])

def load_basis(trainsize, r):
    """Load a POD basis and the associated scales.

    Parameters
    ----------
    trainsize : int
        Number of snapshots used when the SVD was computed.

    r : int
        Number of left singular vectors to load.

    Returns
    -------
    V : (NUM_ROMVARS*DOF,r) ndarray
        POD basis of rank `r`, i.e., the first `r` left singular vectors of
        the training data.

    qbar : (NUM_ROMVARS*DOF,) ndarray
        Mean snapshot that the training data was shifted by after scaling
        but before projection.

    scales : (NUM_ROMVARS,2) ndarray
        Factors used to scale the variables before projecting.
    """
    # Locate the data.
    data_path = _checkexists(config.basis_path(trainsize))

    # Secret! Return list of full singular values.
    if r == -1:
        data_path = data_path.replace(config.BASIS_FILE, "svdvals.h5")
        with h5py.File(data_path, 'r') as hf:
            return hf["svdvals"][:]

    # Extract the data.
    with timed_block(f"Loading POD basis from {data_path}"):
        with h5py.File(data_path, 'r') as hf:
            # Check data shapes.
            rmax = hf["basis"].shape[1]
            if r is not None and rmax < r:
                raise ValueError(f"basis only has {rmax} columns")
            if "mean" in hf:
                if hf["mean"].shape != (hf["basis"].shape[0],):
                    raise RuntimeError("basis and mean snapshot not aligned!")
                mean = hf["mean"][:]
            else:
                mean = np.zeros(hf["basis"].shape[0])

            # Load the data.
            return hf["basis"][:,:r], mean, hf["scales"][:]

def load_projected_data(trainsize, r):
    """Load snapshots that have been projected to a low-dimensional subspace.

    Parameters
    ----------
    trainsize : int
        Number of snapshots to load. This is also the number of
        snapshots that were used when the POD basis (SVD) was computed.

    r : int
        Number of retained POD modes used in the projection.

    Returns
    -------
    Q_ : (r,trainsize) ndarray
        Lifted, scaled, projected snapshots.

    Qdot_ : (r,trainsize) ndarray
        Velocity snapshots corresponding to Q_.

    time_domain : (trainsize) ndarray
        Time domain corresponding to the lifted, scaled data.
    """
    # Locate the data.
    data_path = _checkexists(config.projected_data_path(trainsize))

    # Extract the data.
    with timed_block(f"Loading projected training data from {data_path}"):
        with h5py.File(data_path, 'r') as hf:

            # Check data shapes.
            rmax = hf["data"].shape[0]
            if rmax < r:
                raise ValueError(f"basis only has {rmax} columns")
            if hf["data"].shape[1] != trainsize:
                raise RuntimeError("data set 'data' has incorrect shape")
            if hf["ddt"].shape != hf["data"].shape:
                raise RuntimeError("data sets 'data' and 'ddt' not aligned")
            if hf["time"].shape != (trainsize,):
                raise RuntimeError("data set 'time' has incorrect shape")

            # Get the correct rows of the saved projection data.
            return hf["data"][:r], hf["ddt"][:r], hf["time"][:]

def load_rom(trainsize, r, regs):
    """Load a single trained ROM.

    Parameters
    ----------
    trainsize : int
        Number of snapshots used to train the ROM. This is also the number
        of snapshots that were used when the POD basis (SVD) was computed.

    r : int
        Dimension of the ROM. Also the number of retained POD modes (left
        singular vectors) used to project the training data.

    regs : one, two, or three positive floats
        Regularization hyperparameters used in the Operator Inference
        least-squares problem for training the ROM.

    Returns
    -------
    rom : opinf.InferredContinuousROM
        The trained reduced-order model.
    """
    # Locate the data.
    data_path = _checkexists(config.rom_path(trainsize, r, regs))

    # Extract the trained ROM.
    try:
        rom = opinf.load_model(data_path)
    except FileNotFoundError as e:
        raise DataNotFoundError(f"could not locate ROM with {trainsize:d} "
                                f"training snapshots, r={r:d}, and "
                                f"{config.REGSTR(regs)}") from e
    # Check ROM dimension.
    if rom.r != r:
        raise RuntimeError(f"rom.r = {rom.r} != {r}")

    rom.trainsize, rom.regs = trainsize, regs
    return rom

def load_spatial_statistics(keys, k=None):
    """Load statistical features of the lifted data, computed over the
    spatial domain at each point in time.

    Parameters
    ----------
    keys : list(str)
        Which data set(s) to load. Options:
        * {var}_min : minimum of variable var
        * {var}_max : maximum of variable var
        * {var}_sum : sum (integral) of variable var
        * {var}_std : standard deviation of variable var
        * {var}_mean : mean of variable var
        Here var is a member of config.ROM_VARIABLES. Examples:
        * "T_mean" -> spatially averaged temperature
        * "vx_min" -> minimum x-velocity
        * "CH4_sum" -> methane molar concentration integral

    k : int, slice, or one-dimensional ndarray of sorted integer indices
        Number of time steps of data to load (default all).

    Returns
    -------
    features : dict(str -> (k,) ndarray) or (k,) ndarray
        Dictionary of statistical feature arrays with keys `keys`.
        If only one key is given, return the actual array, not a dict.

    t : (k,) ndarray
        Time domain corresponding to the statistical features.
    """
    # Locate the data.
    data_path = _checkexists(config.statistical_features_path())

    # Parse arguments.
    if isinstance(keys, str):
        keys = [keys]
    elif keys is None:
        keys = ["T_mean"] + [f"{spc}_int" for spc in config.SPECIES]
    if np.isscalar(k) or k is None:
        k = slice(None, k)

    # Extract the data.
    with timed_block(f"Loading statistical features from {data_path}"):
        with h5py.File(data_path, 'r') as hf:
            if len(keys) == 1:
                return hf[f"space/{keys[0]}"][k], hf["t"][k]
            return {key: hf[f"space/{key}"][k] for key in keys}, hf["t"][k]

def load_temporal_statistics(keys):
    """Load statistical features of the lifted data, computed over the
    temporal domain at each spatial point.

    Parameters
    ----------
    keys : list(str)
        Which data set(s) to load. Options:
        * {var}_min : minimum of variable var
        * {var}_max : maximum of variable var
        * {var}_sum : sum (integral) of variable var
        * {var}_std : standard deviation of variable var
        * {var}_mean : mean of variable var
        Here var is a member of config.ROM_VARIABLES. Examples:
        * "T_mean" -> time-averaged temperature
        * "vx_min" -> minimum x-velocity
        * "CH4_sum" -> methane molar concentration time integral

    Returns
    -------
    features : dict(str -> (N,) ndarray) or (N,) ndarray
        Dictionary of statistical feature arrays with keys `keys`.
        If only one key is given, return the actual array, not a dict.
    """
    # Locate the data.
    data_path = _checkexists(config.statistical_features_path())

    # Parse arguments.
    if isinstance(keys, str):
        keys = [keys]

    # Extract the data.
    with timed_block(f"Loading statistical features from {data_path}"):
        with h5py.File(data_path, 'r') as hf:
            if len(keys) == 1:
                return hf[f"time/{keys[0]}"][:]
            return {key: hf[f"time/{key}"][:] for key in keys}

# Figure saving ===============================================================

def save_figure(figname):
    """Save the current matplotlib figure to the figures folder."""
    save_path = os.path.join(config.figures_path(), figname)
    # plt.show() # Uncomment to display figure before saving.
    with timed_block(f"Saving {save_path}"):
        plt.savefig(save_path, bbox_inches="tight", dpi=1200)
        plt.close(plt.gcf())


text:
from os import listdir, path

import matplotlib.pyplot as plt
import numpy as np
from natsort import natsorted

data_files = listdir('data')

numatom_files = [
    file for file in data_files if file.startswith('coords') and '_n' in file
]

numatom_data = {'Base': {}, 'Harmonic': {}, 'Shake': {}}

for numatom_file in natsorted(numatom_files):
    simtype = numatom_file.split('coords_')[1].split('_n')[0]

    with open(path.join('data', numatom_file)) as file:
        num_atoms = int(file.name.split('_n')[1].split('.')[0])
        numatom_data[simtype][num_atoms] = []

        for line in file:
            if not line.startswith('MODEL'):
                break

            start_atom = file.readline()
            for i in range(num_atoms - 2):
                file.readline()
            end_atom = file.readline()
            file.readline()

            start_atom = start_atom.split('UNK     1     ')[1]
            start_atom = start_atom.split('  1.00  0.00')[0]
            end_atom = end_atom.split('UNK     1     ')[1]
            end_atom = end_atom.split('  1.00  0.00')[0]

            start_coords = np.asarray(start_atom.split(' '), dtype=float)
            end_coords = np.asarray(end_atom.split(' '), dtype=float)

            numatom_data[simtype][num_atoms].append(
                np.linalg.norm(start_coords - end_coords))

for simtype, data in numatom_data.items():
    fig, ax = plt.subplots()

    ax.set_title(simtype)
    ax.set_xlabel('NumMDSteps / TrajectoryOutputInterval')
    ax.set_ylabel('End-to-End Distance [$\AA$]')

    for temp, lengths in data.items():
        ax.plot(lengths, label=f'{temp} Atoms')

    fig.savefig(f'{simtype}_numatom.png', bbox_inches='tight')
    plt.close(fig)

max_stars_repo_name: anianruoss/StatPhys

max_issues_repo_name: anianruoss/StatPhys

max_forks_repo_name: anianruoss/StatPhys


text:
import numpy as np
import pandas as pd

df = pd.DataFrame({
    "group": ['cat', 'dog', 'cat', 'cat', 'dog'],
    "group2": ['a', 'a', 'b', 'b', 'b'],
    "A": ['a', 'b', 'c', 'd', 'x'],
    "B": ['e', 'f', 'g', 'h', 'y'],
    "C": np.random.randn(5),
})

def foo(group):
    # blödes Beispiel aber funktioniert auch nur, da Zeilen und spalten gleich groß sind
    # ansonsten könnte DF so nicht erstellt werden
    d = pd.DataFrame({
            'group': group.columns,
            'mean': group["C"] * group["C"].mean()
        }, index=[0, 1, 2, 3, 4])
    return d

print(df)

# gruppe erstellen, auch mit mehreren Spalten
group_df = df.groupby(["group", "group2"])
group_df = df.groupby("group")
print("===\n")
# gibt Index der Gruppen aus
print(group_df.groups)

for name, gruppe in group_df:
    print(f"Der Name ist: {name}")
    print("Und hier kommt die Gruppe:")
    print(gruppe)

print("spezielle Gruppe auslesen")
print(group_df.get_group("cat"))

print("Size der einzelnen Gruppen auslesen")
print(group_df.size())
print("Durchschnitt natürlich nur für Zahlen möglich")
print(group_df.mean())

# Funktionen auf der Gruppe ausführen
print(group_df.apply(foo))

max_stars_repo_name: Kreijeck/learning

max_issues_repo_name: Kreijeck/learning

max_forks_repo_name: Kreijeck/learning


text:
# coding: UTF-8

# works at around 90 FPS (97 FPS if disable imshow)

# Install XIMEA software package
# Copy 'XIMEA\API\Python\v3\ximea' to 'PythonXX\Lib'

from ximea import xiapi
import cv2
import numpy as np

# Connect to camera
cam = xiapi.Camera()
cam.open_device_by_SN('XXXXXXXX')  # Enter serial number of your Ximea camera

# Configuration
cam.set_exposure(5000)
cam.set_gain(24)
cam.set_gammaY(1)

# Allocate memory for ximea image
img = xiapi.Image()

# Start acquisition
cam.start_acquisition()

# Preview output from camera
key = -1
while key == -1:
    cam.get_image(img)
    cvimg = img.get_image_data_numpy()
    cv2.imshow('camera', cvimg)
    key = cv2.waitKey(1)

# Terminate
cam.stop_acquisition()
cam.close_device()