text:
#!/usr/bin/env python
# Copyright 2020 Calico LLC

# 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

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# =========================================================================

from optparse import OptionParser
import glob
import os
import pickle
import shutil
import subprocess
import sys

import h5py
import numpy as np

import slurm

"""
basenji_bench_gtex.py

Compute SNP expression difference scores for variants in VCF files of
fine-mapped GTEx variants to benchmark as features in a classification
task.
"""

################################################################################
# main
################################################################################
def main():
  usage = 'usage: %prog [options] <params_file> <model_file>'
  parser = OptionParser(usage)

  # sad
  parser.add_option('-f', dest='genome_fasta',
      default='%s/data/hg38.fa' % os.environ['BASENJIDIR'],
      help='Genome FASTA for sequences [Default: %default]')
  parser.add_option('--local',dest='local',
      default=1024, type='int',
      help='Local SAD score [Default: %default]')
  parser.add_option('-n', dest='norm_file',
      default=None,
      help='Normalize SAD scores')
  parser.add_option('-o',dest='out_dir',
      default='sad_gtex',
      help='Output directory for tables and plots [Default: %default]')
  parser.add_option('--pseudo', dest='log_pseudo',
      default=1, type='float',
      help='Log2 pseudocount [Default: %default]')
  parser.add_option('--rc', dest='rc',
      default=False, action='store_true',
      help='Average forward and reverse complement predictions [Default: %default]')
  parser.add_option('--shifts', dest='shifts',
      default='0', type='str',
      help='Ensemble prediction shifts [Default: %default]')
  parser.add_option('--stats', dest='sad_stats',
      default='SAD',
      help='Comma-separated list of stats to save. [Default: %default]')
  parser.add_option('-t', dest='targets_file',
      default=None, type='str',
      help='File specifying target indexes and labels in table format')
  parser.add_option('--ti', dest='track_indexes',
      default=None, type='str',
      help='Comma-separated list of target indexes to output BigWig tracks')
  parser.add_option('--threads', dest='threads',
      default=False, action='store_true',
      help='Run CPU math and output in a separate thread [Default: %default]')
  parser.add_option('-u', dest='penultimate',
      default=False, action='store_true',
      help='Compute SED in the penultimate layer [Default: %default]')

  # multi
  parser.add_option('-e', dest='conda_env',
      default='tf2.2-gpu',
      help='Anaconda environment [Default: %default]')
  parser.add_option('-g', dest='gtex_vcf_dir',
      default='/home/drk/seqnn/data/gtex_fine/susie_pip90')
  parser.add_option('--name', dest='name',
      default='gtex', help='SLURM name prefix [Default: %default]')
  parser.add_option('--max_proc', dest='max_proc',
      default=None, type='int',
      help='Maximum concurrent processes [Default: %default]')
  parser.add_option('-p', dest='processes',
      default=None, type='int',
      help='Number of processes, passed by multi script. \
            (Unused, but needs to appear as dummy.)')
  parser.add_option('-q', dest='queue',
      default='gtx1080ti',
      help='SLURM queue on which to run the jobs [Default: %default]')
  parser.add_option('-r', dest='restart',
      default=False, action='store_true',
      help='Restart a partially completed job [Default: %default]')
  (options, args) = parser.parse_args()

  if len(args) != 2:
    parser.error('Must provide parameters and model files')
  else:
    params_file = args[0]
    model_file = args[1]

  #######################################################
  # prep work

  # output directory
  if not options.restart:
    if os.path.isdir(options.out_dir):
      print('Please remove %s' % options.out_dir, file=sys.stderr)
      exit(1)
    os.mkdir(options.out_dir)

  # pickle options
  options_pkl_file = '%s/options.pkl' % options.out_dir
  options_pkl = open(options_pkl_file, 'wb')
  pickle.dump(options, options_pkl)
  options_pkl.close()

  #######################################################
  # predict

  cmd_base = '. /home/drk/anaconda3/etc/profile.d/conda.sh;'
  cmd_base += ' conda activate %s;' % options.conda_env
  cmd_base += ' basenji_sad.py %s %s %s' % (options_pkl_file, params_file, model_file)

  jobs = []
  for gtex_pos_vcf in glob.glob('%s/*_pos.vcf' % options.gtex_vcf_dir):
    # positive job
    job_base = os.path.splitext(os.path.split(gtex_pos_vcf)[1])[0]
    out_dir = '%s/%s' % (options.out_dir, job_base)
    if not options.restart or not os.path.isfile('%s/sad.h5'%out_dir):
      cmd = '%s -o %s %s' % (cmd_base, out_dir, gtex_pos_vcf)
      name = '%s_%s' % (options.name, job_base)
      j = slurm.Job(cmd, name,
          '%s.out'%out_dir, '%s.err'%out_dir,
          queue=options.queue, gpu=1,
          mem=22000, time='1-0:0:0')
      jobs.append(j)

    # negative job
    gtex_neg_vcf = gtex_pos_vcf.replace('_pos.','_neg.')
    job_base = os.path.splitext(os.path.split(gtex_neg_vcf)[1])[0]
    out_dir = '%s/%s' % (options.out_dir, job_base)
    if not options.restart or not os.path.isfile('%s/sad.h5'%out_dir):
      cmd = '%s -o %s %s' % (cmd_base, out_dir, gtex_neg_vcf)
      name = '%s_%s' % (options.name, job_base)
      j = slurm.Job(cmd, name,
          '%s.out'%out_dir, '%s.err'%out_dir,
          queue=options.queue, gpu=1,
          mem=22000, time='1-0:0:0')
      jobs.append(j)

  slurm.multi_run(jobs, max_proc=options.max_proc, verbose=True,
                  launch_sleep=10, update_sleep=60)

  #######################################################
  # classify

  cmd_base = 'basenji_bench_classify.py -i 100 -p 2 -r 44 -s'

  jobs = []
  for gtex_pos_vcf in glob.glob('%s/*_pos.vcf' % options.gtex_vcf_dir):
    tissue = os.path.splitext(os.path.split(gtex_pos_vcf)[1])[0][:-4]
    sad_pos = '%s/%s_pos/sad.h5' % (options.out_dir, tissue)
    sad_neg = '%s/%s_neg/sad.h5' % (options.out_dir, tissue)
    out_dir = '%s/%s_class' % (options.out_dir, tissue)
    if not options.restart or not os.path.isfile('%s/stats.txt' % out_dir):
      cmd = '%s -o %s %s %s' % (cmd_base, out_dir, sad_pos, sad_neg)
      j = slurm.Job(cmd, tissue,
          '%s.out'%out_dir, '%s.err'%out_dir,
          queue='standard', cpu=2,
          mem=22000, time='1-0:0:0')
      jobs.append(j)

  slurm.multi_run(jobs, verbose=True)

def job_completed(options, pi):
  """Check whether a specific job has generated its
     output file."""
  if options.out_txt:
    out_file = '%s/job%d/sad_table.txt' % (options.out_dir, pi)
  elif options.out_zarr:
    out_file = '%s/job%d/sad.zarr' % (options.out_dir, pi)
  elif options.csv:
    out_file = '%s/job%d/sad_table.csv' % (options.out_dir, pi)
  else:
    out_file = '%s/job%d/sad.h5' % (options.out_dir, pi)
  return os.path.isfile(out_file) or os.path.isdir(out_file)

################################################################################
# __main__
################################################################################
if __name__ == '__main__':
  main()

max_stars_repo_name: shtoneyan/basenji

max_issues_repo_name: shtoneyan/basenji

max_forks_repo_name: shtoneyan/basenji


text:
import cv2
import numpy as np

cap = cv2.VideoCapture(0)

def time_warp_scan(rows, cols):
    alternate_frame = np.zeros((rows,cols,3))

    for i in range(rows):
        ret, frame = cap.read()
        if not ret:
            break

        alternate_frame[i, :] = frame[i, :]
        for j in range(i):
            frame[j, :] = alternate_frame[j, :]

        cv2.line(frame, (0, i-1), (cols, i-1), (212, 184, 0), 2)

        cv2.imshow('cam', frame)
        if cv2.waitKey(1) & 0xFF == ord('e'):
            break
    cv2.imshow('cam', frame)
    cv2.waitKey()
    return alternate_frame # in case we want to save or process the final img later, return it to main

def main():
    while True:
        ret, frame = cap.read()
        if not ret:
            break

        cv2.imshow('cam', frame)

        if cv2.waitKey(1) & 0xFF == ord('s'):
            rows, cols, _ = frame.shape
            time_warp_scan(rows, cols)

        if cv2.waitKey(1) & 0xFF == ord('e'):
            break

if __name__=='__main__':
    main()
    cap.release()
    cv2.destroyAllWindows()

ext: py

max_stars_repo_path: Assignment28/time warp scan/main.py

max_issues_repo_path: Assignment28/time warp scan/main.py

max_forks_repo_path: Assignment28/time warp scan/main.py

reason: import numpy


text:
# encoding: utf-8
# usando python 3
# Nota musical - Elaborar um programa que lê uma freqüência (em Hertz) e uma duração (em
# milissegundos) e emite um som na freqüência com a duração.

import pyaudio # folder spelt PyAudio
import numpy as np
import os

os.system('clear')

print("Este programa toca um som de acordo com a frequencia e tempo escolhido pelo usuário")
freq = int(input("Digite a frequencia: "))
tempo = int(input("Digite o tempo em millissegundos: "))

p = pyaudio.PyAudio()

volume = 0.5     # range [0.0, 1.0]
fs = 44100       # sampling rate, Hz, must be integer
duration = tempo/1000   # in seconds, may be float
f = freq        # sine frequency, Hz, may be float

# generate samples, note conversion to float32 array
samples = (np.sin(2*np.pi*np.arange(fs*duration)*f/fs)).astype(np.float32)

# for paFloat32 sample values must be in range [-1.0, 1.0]
stream = p.open(format=pyaudio.paFloat32,
                channels=1,
                rate=fs,
                output=True)

# play. May repeat with different volume values (if done interactively)
stream.write(volume*samples)

stream.stop_stream()
stream.close()

p.terminate()

max_stars_repo_path: 009 Sons/notamusical.py

max_issues_repo_path: 009 Sons/notamusical.py

max_forks_repo_path: 009 Sons/notamusical.py


text:
"""
consists of functions to get a and c weighting information
(DIN EN 61672-1) in different forms:
as weighting values for a given frequency numpy array
as a weighting vector (numpy array) for typical fft proecessing
as a weighting vector (numpy array) for a third-octave filterbank
as a weighting vector (numpy array) for an octave filterbank
as the complex transfer function
as filter coefficients (to be used by lfilter)
"""
# Author: J. Bitzer @ Jade Hochschule (copyright owner)
# License: BSD 3-clause license (https://opensource.org/licenses/BSD-3-Clause)
# Date: 12.12.2021
# used sources: DIN 61672-1, octave toolbox by Christophe COUVREUR
# version 1.0  12.12.2021 first build
# 1.1 14.12.2021 added lower() to allow A and a

import math
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sp

def __get_frequency_constants():
    """helper function with mathematical definition of the design frequencies (f1 to f4) for a and c-weighting curves"""
    fr = 1000 # section 5.4.6 in DIN 61672
    fl = 10**1.5 # section 5.4.6 in DIN 61672
    fh = 10**3.9 # section 5.4.6 in DIN 61672
    D = math.sqrt(1/2) # section 5.4.6 in DIN 61672
    b = (1/(1-D))*(fr**2 + fl**2 * fh**2/fr**2 - D*(fl**2 + fh**2) ) # eq 11 in DIN 61672
    c = fl**2 * fh**2 # eq 12 in DIN 61672
    f1 = ((-b - math.sqrt(b**2 - 4*c))*0.5)**0.5 # eq 9 in DIN 61672
    f4 = ((-b + math.sqrt(b**2 - 4*c))*0.5)**0.5 # eq 10 in DIN 61672
    fa = 10**2.45
    f2 = ((3-math.sqrt(5))*0.5)*fa # eq 13 in DIN 61672
    f3 = ((3+math.sqrt(5))*0.5)*fa # eq 14 in DIN 61672
    return f1,f2,f3,f4

def transfermagnitude2minimumphase(magH):
    '''Computes the minimum phase by a given magnitude of the transfer function by using the Hilbert transform'''
    second_half = magH[-1:0:-1]
    transfer_function = np.append(magH, second_half)
    tf_length = len(transfer_function)
    hilbert_in = np.log(np.abs(transfer_function))

    hilbert_out = -sp.hilbert(hilbert_in)
    phase = hilbert_out[:tf_length//2+1].imag
    return phase

def FDLS_design(order_of_numerator, order_of_denominator, mag_h, phase_h, freq_vek, fs = -1):

    if (len(freq_vek) <= 1):
        omega = np.linspace(start = 0, stop = np.pi, num = len(mag_h))
    else:
        omega = 2*freq_vek/fs*np.pi

    y = mag_h*np.cos(phase_h);
    X_D = -mag_h * np.cos(-1 * omega + phase_h);

    X_D = [X_D]
    for k in range(order_of_denominator-1):
        X_D.append(-mag_h * np.cos(-(k+2) * omega + phase_h))
    X_D = np.array(X_D).T

    # non recursive part
    X_N = np.cos(-1 * omega);

    X_N = [X_N]
    for k in range(order_of_numerator-1):
        X_N.append(np.cos(-(k+2) * omega))
    X_N = np.array(X_N).T

    #and define X as input matrix
    X = np.hstack([X_D, np.ones([len(mag_h), 1]), X_N])
    coeff = np.linalg.lstsq(X, y, rcond=None)
    a = [1, *coeff[0][:order_of_denominator]] # ANN: einfach eine Liste statt hstack (* fügt Elemente ein, anstatt die ganze Liste)
    b = coeff[0][order_of_denominator:order_of_denominator+order_of_numerator+1]
    return b,a

def get_complex_tf_weighting(f_Hz, weight_func = 'a'):
    """returns the complex transfer function for a given frequency vector f_Hz"""
    f1,f2,f3,f4 = __get_frequency_constants()

    om_vek = 2*np.pi*1j*f_Hz

    if weight_func.lower() == 'c':
        c1000 = -0.062
        cweight = ((4*np.pi**2*f4**2*om_vek**2)/((om_vek + 2*np.pi*f1)**2 *(om_vek + 2*np.pi*f4)**2))/(10**(c1000/20))
        return cweight

    if weight_func.lower() == 'a':
        a1000 = -2
        aweight = ((4*np.pi**2*f4**2 * om_vek**4)/((om_vek + 2*np.pi*f1)**2 *(om_vek + 2*np.pi*f2)*
                                                   (om_vek + 2*np.pi*f3)*(om_vek + 2*np.pi*f4)**2))/(10**(a1000/20))

        return aweight

def get_weight_value(f_Hz, weight_func = 'a', return_mode = 'log'):
    """returns the weighting values for a given frequency vector f_Hz in dB(log) or linear, specified by return_mode (default = 'log')"""
    f_Hz[f_Hz == 0] = 0.1 # prevent division by zero
    f1,f2,f3,f4 = __get_frequency_constants()
    if weight_func.lower() == 'c':
        c1000 = -0.062
        # eq 6 in DIN 61672
        if return_mode == 'log':
            cweight = 20*np.log10((f4**2*f_Hz**2)/((f_Hz**2 + f1**2)*(f_Hz**2 + f4**2))) - c1000
        else:
            cweight = ((f4**2*f_Hz**2)/((f_Hz**2 + f1**2)*(f_Hz**2 + f4**2)))/(10**(c1000/20))
        return cweight

    if weight_func.lower() == 'a':
        a1000 = -2
        if return_mode == 'log':
            # eq 7 in DIN 61672
            aweight = 20*np.log10((f4**2*f_Hz**4)/((f_Hz**2 + f1**2)*(f_Hz**2 + f2**2)**0.5*(f_Hz**2 + f3**2)**0.5*(f_Hz**2 + f4**2))) - a1000
        else:
            aweight = ((f4**2*f_Hz**4)/((f_Hz**2 + f1**2)*(f_Hz**2 + f2**2)**0.5*(f_Hz**2 + f3**2)**0.5*(f_Hz**2 + f4**2)))/(10**(a1000/20))
        return aweight

def get_fftweight_vector(fft_size, fs, weight_func = 'a', return_mode = 'log'):
    """ for a given fft_size the a or c weighting vector (fft_size/2 +1 elements from 0 to fs/2 Hz) is returned"""
    freq_vek = np.linspace(0, fs/2, num = int(fft_size/2+1))
    return get_weight_value(freq_vek, weight_func, return_mode), freq_vek

def get_onethirdweight_vector(weight_func = 'a'):
    """ the weights and frequencies for a 34 band one-third filterbank (startfreq = 10, endfreq = 20k) are returned"""
    a, freq_vek = get_spec(weight_func)
    return get_weight_value(freq_vek, weight_func), freq_vek

def get_octaveweight_vector(weight_func = 'a'):
    """the weights and frequencies for a 11 band octave filterbank (startfreq = 16, endfreq = 16k) are returned """
    freq_vek = np.array([16, 31.5, 63, 125, 250, 500, 1000, 2000, 4000, 8000, 16000])
    return get_weight_value(freq_vek, weight_func), freq_vek

def get_weight_coefficients(fs, weight_func = 'a'):
    """designs an IIR filter that implements the desired weighting function for a given sampling rate fs
    (fs > 40 kHz to fullfill class1 specification )
    design idea is based on the matlab octave toolbox. better solutions for lower fs should be possible
    see EOF for some ideas"""

    if fs < 40000:
        print("WARNING: the resulting filter coefficients will not fullfill the class1 constraint")

    f1,f2,f3,f4 = __get_frequency_constants()
    den = np.convolve([1., +4*np.pi * f4, (2*np.pi * f4)**2], [1., +4*np.pi * f1, (2*np.pi * f1)**2])
    if weight_func.lower() == 'c' :
        c1000 = -0.062
        num = [(2*np.pi*f4)**2 * (10**(-c1000 / 20.0)), 0., 0.]

    if weight_func.lower() == 'a':
        a1000 = -2
        den = np.convolve(np.convolve(den, [1., 2*np.pi * f3]),[1., 2*np.pi * f2])
        num = [(2*np.pi*f4)**2 * (10**(-a1000 / 20.0)), 0., 0., 0., 0.]

    b,a = sp.bilinear(num, den, fs)

    return b,a

def get_spec(weight_func = 'a'):
    """ returns the specification of a- and c-weighting function and the class1 limits
    parameter:
        weight_func:  'a' (default) or 'c'
    returns:
        tf_spec: a 34 element numpy-array  the specification for the transfer function of the weighting curve in dB
        f_reference: the corresponding vector with the reference frequencies
        class_1_upperlimit: corresponding vector of the allowed deviation in dB (positive)
        class_1_lowerlimit: corresponding vector of the allowed deviation in dB (negative dB, so you had to add this vector to tf_spec).
    """
    # all values from table 2 in  DIN 61672
    fref = np.array([10, 12.5, 16, 20, 25, 31.5, 40, 50, 63, 80, 100, 125, 160, 200, 250, 315, 400, 500, 630, 800,
                     1000, 1250, 1600, 2000, 2500, 3150, 4000, 5000, 6300, 8000, 10000, 12500, 16000, 20000 ])
    class1_upper_limit = np.array([3.5, 3.0,2.5,2.5, 2.5,2.0,1.5, 1.5,1.5,1.5, 1.5,1.5,1.5, 1.5,1.4,1.4, 1.4,1.4,1.4,
                                     1.4,1.1,1.4, 1.6,1.6,1.6, 1.6,1.6,2.1, 2.1,2.1,2.6, 3,3.5,4 ])
    class1_lower_limit = np.array([-np.inf, -np.inf,-4.5,-2.5, -2,-2,-1.5, -1.5,-1.5,-1.5, -1.5,-1.5,-1.5,  -1.5,-1.4,-1.4,
                                   -1.4,-1.4,-1.4, -1.4,-1.1,-1.4, -1.6,-1.6,-1.6, -1.6,-1.6,-2.1, -2.6,-3.1,-3.6, -6.0,-17.0,-np.inf  ])
    if (weight_func.lower() == 'a'):
        a_desired = np.array([-70.4, -63.4, -56.7, -50.5, -44.7, -39.4, -34.6, -30.2, -26.2, -22.5, -19.1, -16.1, -13.4,
                              -10.9, -8.6, -6.6, -4.8, -3.2, -1.9, -.8, 0, .6, 1, 1.2, 1.3, 1.2, 1, .5, -.1, -1.1, -2.5, -4.3, -6.6, -9.3 ])
        return a_desired, fref, class1_lower_limit,  class1_upper_limit

    if (weight_func.lower() == 'c'):
        c_desired = np.array([ -14.3, -11.2, -8.5, -6.2, -4.4, -3.0, -2.0, -1.3, -.8, -.5, -.3, -.2, -.1, 0, 0, 0, 0, 0, 0,
                              0, 0, 0, -.1, -.2, -.3, -.5, -.8, -1.3, -2.0, -3.0, -4.4, -6.2, -8.5, -11.2 ])
        return c_desired, fref, class1_lower_limit,  class1_upper_limit

if __name__ == '__main__':
    print('called as script')
    weight_type  = 'A'
    # print(get_weight_value(np.array([10,100,1000]),weight_type))
    a, f, ll, hl = get_spec(weight_type)
    fig,ax = plt.subplots()
    ax.plot(np.log(f),a, 'r')
    ax.plot(np.log(f),a+hl,'g:')
    ax.plot(np.log(f),a+ll, 'g:')
    fs = 48000
    afft,f_fft = get_fftweight_vector(1024,fs, weight_type)
    ax.plot(np.log(f_fft),afft, 'y+')
    b,a = get_weight_coefficients(fs,weight_type)
    #print (b)
    #print (a)
    w, H = sp.freqz(b, a, worN =2048, fs = fs)
    ax.plot(np.log(w[1:]),20*np.log10(np.abs(H[1:])),'k')
    ax.set_ylim([-70.0, 5.0])
    ax.set_xlim([np.log(10), np.log(20000)])
    ax.set_xticks(np.log(f[:-1:4]))
    ax.set_xticklabels((f[:-1:4]))
    plt.show()

"""
an idea to design the a and c filter by FDLS. it does not work,
another better solution would be, wo use SOS filter and correct the tf at Nyquist.
next idea: use a better IIR arbitrary designb routine

    #fft_size = 2*4096
    #tf,f = get_fftweight_vector(fft_size,fs,weight_func,'lin')
    #phase = transfermagnitude2minimumphase(np.sqrt(tf))

    # log weighting for LS Design
    #f_log = np.logspace(1.0,4.05, num = 150)
    f_log = np.linspace(100,10000, num = 200)
    # print(f_log)
    #index_f = (np.round( (2*f_log/fs*(fft_size/2+1)))).astype(int)
    #print(index_f)
    tf = get_complex_tf_weighting(f_log, weight_func)

    tf_mag = np.abs(tf)
    tf_phase = np.angle(tf)

    fig2,ax2 = plt.subplots(ncols=2)
    ax2[0].plot(f_log,20*np.log10(tf_mag))

    ax2[1].plot(f_log,np.unwrap(tf_phase))
    #plt.show()

    if weight_func == 'a':
        order_num = 8
        order_den = 8
    else:
        order_num = 4
        order_den = 4

    b,a = FDLS_design(order_num, order_den,tf_mag, np.unwrap(tf_phase), f_log, fs)

    return b,a

"""

"""Copyright <2021> <Joerg Bitzer>

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.

3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "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 COPYRIGHT HOLDER OR CONTRIBUTORS 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,import scipy


text:
###################################################################
###########  Setting up the workplace #############################
###################################################################
import os

import numpy as np

import pandas as pd

import pyreadr

pd.set_option('mode.chained_assignment', None)

''' Set the correct work directory '''

os.getcwd()

os.chdir('/Users/edgarakopyan/Desktop/Measuring Unfair Inequality/Data')

# We need to read the RData file using pyreadr package
data = pyreadr.read_r('RLMS-HSE_IND_1994_2018_STATA.RData')

# And finally to get the dataframe we extract it from the created dictionary
workdata = data["table"]

# Let's make a copy of the dataset and work on the copy for now and for coding convenience we reset index
samplework = workdata.copy()

# And select the variables we want to work with
samplework = samplework[['J8', 'J38', 'J59', 'I4', 'H5', 'idind', 'J217A', 'J217B', 'J216AC08', 'J216BC08', \
                           'age','J60', 'year', 'ID_H', 'origsm', 'inwgt']]

###################################################################
###########  Cleaning parent occupation and education #############
###################################################################

# Parental education is present only in two years. Let's fill in parent education with the maximum value
# (both mother and father variables)

samplework['J217A'] = samplework.groupby('idind')['J217A'].transform(lambda x: x.fillna(x.max()))

samplework['J217B'] = samplework.groupby('idind')['J217B'].transform(lambda x: x.fillna(x.max()))

# Same should be done with parental occupation (both mother and father variables)
samplework['J216AC08'] = samplework.groupby('idind')['J216AC08'].transform(lambda x: x.fillna(x.max()))

samplework['J216BC08'] = samplework.groupby('idind')['J216BC08'].transform(lambda x: x.fillna(x.max()))

# Let's get rid of people not of working age
samplework = samplework[samplework.age.isin(range(25,60))]

# Now we need to remove rows with missing parent occupation and education.
# We remove both NA values as well as answers like don't know and refused to answer.
samplework = samplework[samplework['J216AC08'].notna()]

samplework = samplework[samplework['J216BC08'].notna()]

samplework = samplework[samplework['J217A'].notna()]

samplework = samplework[samplework['J217B'].notna()]

samplework = samplework[samplework['J216AC08'] != 99999997]

samplework = samplework[samplework['J216AC08'] != 99999996]

samplework = samplework[samplework['J216AC08'] != 99999995]

samplework = samplework[samplework['J216AC08'] != 99999998]

samplework = samplework[samplework['J216AC08'] != 99999994]

samplework = samplework[samplework['J216AC08'] != 99999999]

samplework = samplework[samplework['J216AC08'] != 99999993]

samplework = samplework[samplework['J216BC08'] != 99999997]

samplework = samplework[samplework['J216BC08'] != 99999996]

samplework = samplework[samplework['J216BC08'] != 99999995]

samplework = samplework[samplework['J216BC08'] != 99999998]

samplework = samplework[samplework['J216BC08'] != 99999999]

samplework = samplework[samplework['J216BC08'] != 99999994]

samplework = samplework[samplework['J216BC08'] != 99999993]

samplework = samplework[samplework['J217A'] != 99999999]

samplework = samplework[samplework['J217A'] != 99999996]

samplework = samplework[samplework['J217A'] != 99999993]

samplework = samplework[samplework['J217A'] != 99999998]

samplework = samplework[samplework['J217A'] != 99999997]

samplework = samplework[samplework['J217B'] != 99999999]

samplework = samplework[samplework['J217B'] != 99999998]

samplework = samplework[samplework['J217B'] != 99999997]

samplework = samplework[samplework['J217B'] != 99999996]

samplework = samplework[samplework['J217B'] != 99999993]

samplework = samplework.reset_index(drop=True)

# Then we clear the occupation variable by taking only the first digit  and structure it into 3 distinct groups

for i in range(0, samplework.shape[0]):
    while samplework.J216AC08[i] >= 10:
        samplework.J216AC08[i] = samplework.J216AC08[i] // 10

for i in range(0, samplework.shape[0]):
    while samplework.J216BC08[i] >= 10:
        samplework.J216BC08[i] = samplework.J216BC08[i] // 10

# And finally we create a new variable that divides occupation into 3 groups following the original Measuring
# Unfair Inequality paper

samplework['fatherwork'] = None

for i in range(0, samplework.shape[0]):
    if samplework['J216AC08'][i] in [1, 2, 3]:
        samplework['fatherwork'][i] = '2'
    elif samplework['J216AC08'][i] in [4, 5, 7, 8]:
        samplework['fatherwork'][i] = '1'
    else:
        samplework['fatherwork'][i] = '0'

samplework['fatherwork'].value_counts()

samplework['motherwork'] = None

for i in range(0, samplework.shape[0]):
    if samplework['J216BC08'][i] in [1, 2, 3]:
        samplework['motherwork'][i] = '2'
    elif samplework['J216BC08'][i] in [4, 5, 7, 8]:
        samplework['motherwork'][i] = '1'
    else:
        samplework['motherwork'][i] = '0'

# And now we create a variable of maximum parent occupation for each individual
samplework['maxparentwork'] = None

for i in range(0, samplework.shape[0]):
    samplework['maxparentwork'][i] = max(samplework['motherwork'][i], samplework['fatherwork'][i])

# Now we need to do the same with education variable

samplework['fathereducation'] = None

for i in range(0, samplework.shape[0]):
    if samplework['J217A'][i] in [9, 10, 11]:
        samplework['fathereducation'][i] = 2
    elif samplework['J217A'][i] in [4, 5, 6, 7, 8]:
        samplework['fathereducation'][i] = 1
    else:
        samplework['fathereducation'][i] = 0

samplework['fathereducation'].value_counts()

samplework['mothereducation'] = None

for i in range(0, samplework.shape[0]):
    if samplework['J217B'][i] in [9, 10, 11]:
        samplework['mothereducation'][i] = 2
    elif samplework['J217B'][i] in [4, 5, 6, 7, 8]:
        samplework['mothereducation'][i] = 1
    else:
        samplework['mothereducation'][i] = 0

# And now we create a variable of maximum parent education
samplework['maxparented'] = None

for i in range(0, samplework.shape[0]):
    samplework['maxparented'][i] = max(samplework['mothereducation'][i], samplework['fathereducation'][i])

###################################################################
###########  Cleaning national (ethnic)  background ###############
###################################################################

# Now on to national background
# Let's look at the share of missing values
len(samplework[samplework['I4'] == 99999997]) + len(samplework[samplework['I4'] == 99999998]) +\
    len(samplework[samplework['I4'] == 99999999]) + len(samplework[samplework['I4'] == np.nan])

# We have 436 missing values. Let's see how many slavic peoples we have
(len(samplework[samplework['I4'] == 2]) + len(samplework[samplework['I4'] == 1]) +\
    len(samplework[samplework['I4'] == 4]))/samplework.shape[0]
samplework.I4.value_counts()
# 89% of our dataset are slavic peoples
# As we can see, Slavs heavily dominate this variable and for many nationalities have only 1 observation. This means
# that it will be impossible to construct sufficiently big types to compare them to each other.
# However, we can divided into Slavic and non-slavic groups. We first get rid of missing observations and then divide.
samplework = samplework[(samplework.I4 != 99999997) & (samplework.I4 != 99999998) & (samplework.I4 != 99999999)]

samplework = samplework[samplework.I4 != np.nan]

# And now we divide into the two groups. I will bundle Slavs and Europeans into one group and the rest into another.
samplework['nationality'] = None

for i in list(samplework.index):
    if samplework['I4'][i] in [1, 2, 4, 5, 7, 10, 13]:
        samplework['nationality'][i] = 1
    else:
        samplework['nationality'][i] = 0

###################################################################
###########  Cleaning outcome variable ############################
###################################################################

# Let's look at the outcome variable. No sudden jumps in dont knows
samplework[samplework.J60 == 99999997].groupby(samplework.year).J60.value_counts() / samplework.groupby(\
    samplework.year).size()

# No sudden jumps in refusals.
samplework[samplework.J60 == 99999998].groupby(samplework.year).J60.value_counts() / samplework.groupby(\
    samplework.year).size()

# Relatively high in 1994 and 1995 but otherwise low in other years
samplework[samplework.J60 == 99999999].groupby(samplework.year).J60.value_counts() / samplework.groupby(\
    samplework.year).size()

# More 0 in 1996 and 1998 but otherwise also consistent
samplework[samplework.J60 == 0].groupby(samplework.year).J60.value_counts() / samplework.groupby(\
    samplework.year).size()

# Hence we can remove unnecessary values
samplework = samplework[~((samplework.J60 == 99999997.0) | (samplework.J60 == 99999998.0) | (samplework.J60 == \
                                                                                99999999.0))]
# And we also need to remove observations where we have nonzero working hours but zero salary
samplework = samplework[~((samplework.J60 == 0) & ((samplework.J8 > 0) | (samplework.J38 > 0)| (samplework.J59 > 0)))]

samplework['income'] = samplework.J60
# We also need  to adjust for inflation. I use IMF inflation data as opposed to Penn Table because the latter goes only
# until 2017.

inflation = pd.read_excel('imf-dm-export-20200516.xls')

inflation = inflation.iloc[1, 15:40]  # Select the valid years

inflation = inflation/100 + 1
# We have to multiply each consequentive element by the previous one to see how much price levels have changed over the
# period
for i in range(1, inflation.shape[0]):
    inflation.iloc[i] = inflation.iloc[i] * inflation.iloc[i-1]

finalinflation = inflation.iloc[inflation.shape[0]-1]

for i in range(0, inflation.shape[0]):
    inflation.iloc[i] = (inflation.iloc[i])/finalinflation

inflation = inflation.drop([1997, 1999])  # Since we do not have those years anyway

# Now we adjust the income by inflation
samplework = samplework.reset_index(drop=True)

for i in range(1,samplework.shape[0]):
    samplework.income[i] = samplework.income[i] / inflation[int(samplework.year[i])]

# We also have to take into account the redonimination of ruble in 1998 with the rate 1000 to 1
samplework.income[samplework.year < 1998] = samplework.income[samplework.year < 1998]/1000

# Finally replace 0 incomes with 1 to avoid problems when applying logarithm
samplework.income[samplework.income == 0] = 1

###################################################################
###########  Create types #########################################
###################################################################

# And we make types
samplework['type'] = samplework['maxparentwork'].astype('str') + '_' + samplework['maxparented'].astype('str') + \
  '_' + samplework['H5'].astype('str') + '_' + samplework['nationality'].astype('str')

types = list((samplework.type.value_counts()).index)

###################################################################
###########  Save dataset #########################################
###################################################################

# Save this file to be used in the next part
samplework.to_pickle('Intermediarywork')

max_stars_repo_path: Data preparation.py

max_issues_repo_path: Data preparation.py

max_forks_repo_path: Data preparation.py


text:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# pos 读入的所有数据
# step 时间间隔
# i 第几列数据 (从1开始)
# name 图片名称，比如 acceleration.png
def picturePosition(data,step,i,name):
	i=i-1
	col = np.array([float(x.split(' ')[i]) for x in data])
	length = len(col)
	total = length*step
	xtic = np.arange(11)*(total/10)
	x = np.arange(length)*step
	plt.figure(figsize=(10,10))
	plt.plot(x,col)
	plt.xticks(xtic)
	plt.xlabel('time',fontsize=18)
	plt.ylabel('position',fontsize=18)
	plt.savefig(name)
	#plt.show()

def pictureVelocity(data,step,i,name):
	i=i-1
	col = np.array([float(x.split(' ')[i]) for x in data])
	length = len(col)
	total = length*step
	xtic = np.arange(11)*(total/10)
	x = np.arange(length)*step
	plt.figure(figsize=(10,10))
	plt.plot(x,col)
	plt.xticks(xtic)
	plt.xlabel('time',fontsize=18)
	plt.ylabel('velocity',fontsize=18)
	plt.savefig(name)
	#plt.show()

def pictureAcceleration(data,step,i,name):
	i=i-1
	col = np.array([float(x.split(' ')[i]) for x in data])
	length = len(col)
	total = length*step
	xtic = np.arange(11)*(total/10)
	x = np.arange(length)*step
	plt.figure(figsize=(10,10))
	plt.plot(x,col)
	plt.xticks(xtic)
	plt.xlabel('time',fontsize=18)
	plt.ylabel('acceleration',fontsize=18)
	plt.savefig(name)
	#plt.show()

def pictureContact(data,step,i,name):
	i=i-1
	col = np.array([float(x.split(' ')[i]) for x in data])
	length = len(col)
	total = length*step
	xtic = np.arange(11)*(total/10)
	x = np.arange(length)*step
	plt.figure(figsize=(10,10))
	plt.plot(x,col)
	plt.xticks(xtic)
	plt.xlabel('time',fontsize=18)
	plt.ylabel('contact force',fontsize=18)
	plt.savefig(name)
	#plt.show()

def pictureContactPosition(data,step,i,name):
	i=i-1
	col = np.array([float(x.split(' ')[i]) for x in data])
	length = len(col)
	total = length*step
	xtic = np.arange(11)*(total/10)
	x = np.arange(length)*step
	plt.figure(figsize=(10,10))
	plt.plot(x,col)
	plt.xticks(xtic)
	plt.xlabel('time',fontsize=18)
	plt.ylabel('contact position',fontsize=18)
	plt.savefig(name)
	#plt.show()

if __name__ == "__main__":
	with open('position.txt','r') as f:
		data1 = f.readlines()
	with open('velocity.txt','r') as f:
		data2 = f.readlines()
	with open('acceleration.txt','r') as f:
		data3 = f.readlines()
	with open('contact.txt','r') as f:
		data4 = f.readlines()
	with open('contactPosition.txt','r') as f:
		data5 = f.readlines()

	picturePosition(data1,0.005,3,'position.png')
	pictureVelocity(data2,0.005,3,'velocity.png')
	pictureAcceleration(data3,0.005,3,'acceleration.png')
	pictureContact(data4,0.005,1,'contact.png')
	pictureContactPosition(data5,0.005,3,'contactPosition.png')

max_stars_repo_name: yuanwentian/TestMujocoPhysics

max_issues_repo_name: yuanwentian/TestMujocoPhysics

max_forks_repo_name: yuanwentian/TestMujocoPhysics


text:
#!/usr/bin/env python
from nose.tools import *
from nose import SkipTest
import networkx as nx
from base_test import BaseTestAttributeMixing, BaseTestDegreeMixing
from networkx.algorithms.assortativity.correlation import attribute_ac

class TestDegreeMixingCorrelation(BaseTestDegreeMixing):
    @classmethod
    def setupClass(cls):
        global np
        global npt
        try:
            import numpy as np
            import numpy.testing as npt
        except ImportError:
            raise SkipTest('NumPy not available.')
        try:
            import scipy
            import scipy.stats
        except ImportError:
            raise SkipTest('SciPy not available.')

    def test_degree_assortativity_undirected(self):
        r = nx.degree_assortativity_coefficient(self.P4)
        npt.assert_almost_equal(r, -1.0 / 2, decimal=4)

    def test_degree_assortativity_directed(self):
        r = nx.degree_assortativity_coefficient(self.D)
        npt.assert_almost_equal(r, -0.57735, decimal=4)

    def test_degree_assortativity_multigraph(self):
        r = nx.degree_assortativity_coefficient(self.M)
        npt.assert_almost_equal(r, -1.0 / 7.0, decimal=4)

    def test_degree_assortativity_undirected(self):
        r = nx.degree_pearson_correlation_coefficient(self.P4)
        npt.assert_almost_equal(r, -1.0 / 2, decimal=4)

    def test_degree_assortativity_directed(self):
        r = nx.degree_pearson_correlation_coefficient(self.D)
        npt.assert_almost_equal(r, -0.57735, decimal=4)

    def test_degree_assortativity_multigraph(self):
        r = nx.degree_pearson_correlation_coefficient(self.M)
        npt.assert_almost_equal(r, -1.0 / 7.0, decimal=4)

class TestAttributeMixingCorrelation(BaseTestAttributeMixing):
    @classmethod
    def setupClass(cls):
        global np
        global npt
        try:
            import numpy as np
            import numpy.testing as npt

        except ImportError:
            raise SkipTest('NumPy not available.')

    def test_attribute_assortativity_undirected(self):
        r = nx.attribute_assortativity_coefficient(self.G, 'fish')
        assert_equal(r, 6.0 / 22.0)

    def test_attribute_assortativity_directed(self):
        r = nx.attribute_assortativity_coefficient(self.D, 'fish')
        assert_equal(r, 1.0 / 3.0)

    def test_attribute_assortativity_multigraph(self):
        r = nx.attribute_assortativity_coefficient(self.M, 'fish')
        assert_equal(r, 1.0)

    def test_attribute_assortativity_coefficient(self):
        # from "Mixing patterns in networks"
        a = np.array([[0.258, 0.016, 0.035, 0.013],
                      [0.012, 0.157, 0.058, 0.019],
                      [0.013, 0.023, 0.306, 0.035],
                      [0.005, 0.007, 0.024, 0.016]])
        r = attribute_ac(a)
        npt.assert_almost_equal(r, 0.623, decimal=3)

    def test_attribute_assortativity_coefficient2(self):
        a = np.array([[0.18, 0.02, 0.01, 0.03],
                      [0.02, 0.20, 0.03, 0.02],
                      [0.01, 0.03, 0.16, 0.01],
                      [0.03, 0.02, 0.01, 0.22]])

        r = attribute_ac(a)
        npt.assert_almost_equal(r, 0.68, decimal=2)

    def test_attribute_assortativity(self):
        a = np.array([[50, 50, 0], [50, 50, 0], [0, 0, 2]])
        r = attribute_ac(a)
        npt.assert_almost_equal(r, 0.029, decimal=3)

max_stars_repo_name: Richoor/HEFT

max_issues_repo_name: Richoor/HEFT

max_forks_repo_name: Richoor/HEFT

reason: import numpy,import scipy,import networkx,from networkx


text:
# -*- coding: utf-8 -*-
"""
window.py
======================================
Introduces the Window module that is used when splitting the path over:
    - Global
    - Sliding
    - Expanding
    - Dyadic
window types.
Code based on window code written by Patrick Kidger.
"""
import collections as co
import numpy as np

_Pair = co.namedtuple("Pair", ("start", "end"))

def _window_getter(
    window_name, window_depth=None, window_length=None, window_step=None
):
    """Gets the window method correspondent to the given string and initialises
    with specified parameters.

    Parameters
    ----------
    window_name: str, String from ['global', 'sliding', 'expanding', 'dyadic']
        used to access the window method.
    window_depth: int, The depth of the dyadic window. (Active only if
        `window_name == 'dyadic']`.
    window_length: int, The length of the sliding/expanding window. (Active
        only if `window_name in ['sliding, 'expanding'].
    window_step: int, The step of the sliding/expanding window. (Active
        only if `window_name in ['sliding, 'expanding'].

    Returns
    -------
    list:
        A list of lists where the inner lists are lists of tuples that
        denote the start and end indexes of each window.
    """
    # Setup all available windows here
    length_step = {"length": window_length, "step": window_step}
    window_dict = {
        "global": (_Global, {}),
        "sliding": (_Sliding, length_step),
        "expanding": (_Expanding, length_step),
        "dyadic": (_Dyadic, {"depth": window_depth}),
    }

    if window_name not in window_dict.keys():
        raise ValueError(
            "Window name must be one of: {}. Got: {}.".format(
                window_dict.keys(), window_name
            )
        )

    window_cls, window_kwargs = window_dict[window_name]
    return window_cls(**window_kwargs)

class _Window:
    """Abstract base class for windows.

    Each subclass must implement a __call__ method that returns a list of lists
    of 2-tuples. Each 2-tuple specifies the start and end of each window.

    These windows are grouped into a list that will (usually) cover the full
    time series. These lists are grouped into another list for situations
    where we consider windows of multiple scales.
    """

    def num_windows(self, length):
        """Method that returns the total number of windows in the set.
        Parameters
        ----------
        length: int, The length of the input path.
        Returns
        -------
        int: The number of windows.
        """
        return sum([len(w) for w in self(length)])

class _Global(_Window):
    """A single window over the full data."""

    def __call__(self, length=None):
        return [[_Pair(None, None)]]

class _ExpandingSliding(_Window):
    def __init__(self, initial_length, start_step, end_step):
        """
        Parameters
        ----------
        initial_length: int, Initial length of the input window.
        start_step: int, Initial step size.
        end_step: int, Final step size.
        """
        super(_ExpandingSliding, self).__init__()
        self.initial_length = initial_length
        self.start_step = start_step
        self.end_step = end_step

    def __call__(self, length):
        def _call():
            start = 0
            end = self.initial_length
            while end <= length:
                yield _Pair(start, end)
                start += self.start_step
                end += self.end_step

        windows = list(_call())
        if len(windows) == 0:
            raise ValueError(
                "Length {} too short for given window parameters.".format(length)
            )
        return [windows]

class _Sliding(_ExpandingSliding):
    """A window starting at zero and going to some point that increases
    between windows.
    """

    def __init__(self, length, step):
        """
        Parameters
        ----------
        length: int, The length of the window.
        step: int, The sliding step size.
        """
        super(_Sliding, self).__init__(
            initial_length=length, start_step=step, end_step=step
        )

class _Expanding(_ExpandingSliding):
    """A window of fixed length, slid along the dataset."""

    def __init__(self, length, step):
        """
        Parameters
        ----------
        length: int, The length of each window.
        step: int, The step size.
        """
        super(_Expanding, self).__init__(
            initial_length=length, start_step=0, end_step=step
        )

class _Dyadic(_Window):
    """Windows generated 'dyadically'.

    These are successive windows of increasing fineness. The windows are as
    follows:
        Depth 1: The global window over the full data.
        Depth 2: The windows of the first and second halves of the dataset.
        Depth 3: The dataset is split into quarters, and we take the windows of
            each quarter.
        ...
        Depth n: For a dataset of length L, we generate windows
            [0:L/(2^n), L/(2^n):(2L)/(2^n), ..., (2^(n-1))L/2^n:L].
    Each depth also contains all previous depths.

    Note: Ensure the depth, n, is chosen such that L/(2^n) >= 1, else it will
        be too high for the dataset.

    Parameters
    ----------
    depth: int, The depth of the dyadic window, explained in the class
        description.
    """

    def __init__(self, depth):
        super(_Dyadic, self).__init__()
        self.depth = depth

    def __call__(self, length):
        max_depth = int(np.floor(np.log2(length)))
        if self.depth > max_depth:
            raise ValueError(
                "Chosen dyadic depth is too high for the data length. "
                "We require depth <= {} for length {}. "
                "Depth given is: {}.".format(max_depth, length, self.depth)
            )
        return self.call(float(length))

    def call(self, length, _offset=0.0, _depth=0, _out=None):
        if _out is None:
            _out = [[] for _ in range(self.depth + 1)]
        _out[_depth].append(_Pair(int(_offset), int(_offset + length)))

        if _depth < self.depth:
            left = _Dyadic(self.depth)
            right = _Dyadic(self.depth)
            half_length = length / 2
            # The order of left then right is important, so that they add their
            # entries to _out in the correct order.
            left.call(half_length, _offset, _depth + 1, _out)
            right.call(half_length, _offset + half_length, _depth + 1, _out)

        return _out

max_stars_repo_name: mikofski/sktime

max_issues_repo_name: mikofski/sktime

max_forks_repo_name: mikofski/sktime


text:
import numpy as np
from astropy.visualization.wcsaxes.patches import _rotate_polygon
from matplotlib.patches import PathPatch
from matplotlib.path import Path
import astropy.units as u

class SphericalCircle(PathPatch):
    # created from the astropy.visualization.wcsaxes.patches.SphericalCircle class
    # changed to path from polygon to create disjointed parts
    """
    Create a patch representing a spherical circle - that is, a circle that is
    formed of all the points that are within a certain angle of the central
    coordinates on a sphere. Here we assume that latitude goes from -90 to +90

    This class is needed in cases where the user wants to add a circular patch
    to a celestial image, since otherwise the circle will be distorted, because
    a fixed interval in longitude corresponds to a different angle on the sky
    depending on the latitude.

    Parameters
    ----------
    center : tuple or `~astropy.units.Quantity`
        This can be either a tuple of two `~astropy.units.Quantity` objects, or
        a single `~astropy.units.Quantity` array with two elements.
    radius : `~astropy.units.Quantity`
        The radius of the circle
    resolution : int, optional
        The number of points that make up the circle - increase this to get a
        smoother circle.
    vertex_unit : `~astropy.units.Unit`
        The units in which the resulting polygon should be defined - this
        should match the unit that the transformation (e.g. the WCS
        transformation) expects as input.

    Notes
    -----
    Additional keyword arguments are passed to `~matplotlib.patches.Polygon`
    """

    def __init__(self, ra, dec, radius, resolution=100, vertex_unit=u.degree, **kwargs):

        # Extract longitude/latitude, either from a tuple of two quantities, or
        # a single 2-element Quantity.

        longitude, latitude = ra, dec

        # #longitude values restricted on domain of (-180,180]
        # if longitude.to_value(u.deg) > 180. :
        # 	longitude = -360. * u.deg + longitude.to(u.deg)

        # Start off by generating the circle around the North pole
        lon = np.linspace(0.0, 2 * np.pi, resolution + 1)[:-1] * u.radian
        lat = np.repeat(0.5 * np.pi - np.deg2rad(radius), resolution) * u.radian

        lon, lat = _rotate_polygon(lon, lat, longitude, latitude)

        # Extract new longitude/latitude in the requested units
        lon = lon.to_value(vertex_unit)
        lat = lat.to_value(vertex_unit)
        # Create polygon vertices
        vertices = np.array([lon, lat]).transpose()

        # split path into two sections if circle crosses -180, 180 bounds
        codes = []
        last = (4000.4 * u.degree).to_value(
            vertex_unit
        )  # 400.4 is a random number large enough so first element is "MOVETO"
        for v in vertices:
            if np.absolute(v[0] - last) > (300 * u.degree).to_value(vertex_unit):
                codes.append(Path.MOVETO)
            else:
                codes.append(Path.LINETO)
            last = v[0]

        circle_path = Path(vertices, codes)

        super().__init__(circle_path, **kwargs)

max_stars_repo_name: drJfunk/gbmgeometry

max_issues_repo_name: drJfunk/gbmgeometry

max_forks_repo_name: drJfunk/gbmgeometry

reason: import numpy,import astropy,from astropy


text:
import click
import glob
import numpy as np
import geopandas as gpd
import pandas as pd
import matplotlib.pyplot as plt
from copro import utils, evaluation
import os, sys

@click.command()
@click.option('-t0', '--start-year', type=int)
@click.option('-t1', '--end-year', type=int)
@click.option('-c', '--column', help='column name', default='chance_of_conflict', type=str)
@click.option('--geojson/--no-geojson', help='save output to geojson or not', default=False)
@click.option('--png/--no-png', help='save output to png or not', default=True)
@click.option('--verbose/--no-verbose', help='verbose on/off', default=False)
@click.argument('input_dir', type=click.Path())
@click.argument('output_dir',type=click.Path())
@click.argument('selected_polygons', type=click.Path())

def main(input_dir=None, output_dir=None, selected_polygons=None, start_year=None, end_year=None, geojson=None, png=None, column=None, verbose=None):
    """Post-processing script to calculate average model output over a user-specifeid period or all output geoJSON-files stored in input-dir.
    Computed average values can be outputted as geoJSON-file or png-file or both.

    Args:
        input_dir: path to input directory with geoJSON-files located per projection year.
        output_dir (str): path to directory where output will be stored.
        selected_polygons (str): path to a shp-file with all polygons used in a CoPro run.

    Output:
        geoJSON-file with average column value per polygon (if geojson is set).
        png-file with plot of average column value per polygon (if png is set)
    """

    # check if start/end time settings are consistent
    if ((start_year != None) and (end_year == None)) or ((end_year != None) and (start_year == None)):
        raise ValueError('ERROR: if start or end year is specified, the pendant must be specified too!')

    # read a shp-file with geometries of all selected polygons
    click.echo('\nreading shp-file with all polygons from {}'.format(os.path.abspath(selected_polygons)))
    selected_polygons_gdf = gpd.read_file(os.path.abspath(selected_polygons))
    # create dataframe
    global_df = utils.global_ID_geom_info(selected_polygons_gdf)

    # find all geojson-files in input-dir
    input_dir = os.path.abspath(input_dir)
    click.echo('getting geojson-files from {}'.format(input_dir))
    all_files = sorted(glob.glob(os.path.join(input_dir, '*.geojson')))

    # if specific start/end time is specified, find only those geojson-files for specified period
    if (start_year != None) and (end_year != None):
        # define period between start and ent time
        period = np.arange(start_year, end_year+1, 1)
        click.echo('using all geojson-files for years {} to {}'.format(period[0], period[-1]))
        # creating suffix for file saving later
        suffix = '{}_to_{}'.format(period[0], period[-1])
        # initiate empty list for selected geojson-files
        selected_files = []
        # select
        for fo in all_files:
            # if the year-suffix of geojson-file matches year in period, add to list
            year = int(str(str(os.path.basename(fo)).rsplit('.')[0]).rsplit('_')[-1])
            if year in period:
                if verbose: print('adding to selection')
                selected_files.append(fo)
    # if not end/start time is specified, use all geojson-files in input-dir
    else:
        click.echo('using all geojson-files in input-dir')
        # also here, create suffix for file saving laster
        suffix = 'all_years'
        selected_files = all_files

    # initiatie empyt dataframe for collecting all annual output
    y_df = pd.DataFrame()

    # go through all geojson-files left over after selection step
    for geojson in selected_files:
        # read files and convert to datatrame
        if verbose: click.echo('reading file {}'.format(geojson))
        gdf = gpd.read_file(geojson, driver='GeoJSON')
        df = pd.DataFrame(gdf)
        # append to dataframe
        y_df = y_df.append(df, ignore_index=True)

    # initiate dataframe for time-averaged output
    click.echo('creating one output dataframe from all geojson-files')
    y_out = pd.DataFrame()
    # get all unique IDs of polygons
    y_out['ID'] = y_df.ID.unique()
    click.echo('reading from column {}'.format(column))
    if column == 'chance_of_conflict':
        # add number of predictiosn made over all selected years
        y_out = pd.merge(y_out, y_df.nr_predictions.groupby(y_df.ID).sum().to_frame(), on='ID')
        # add number of predicted conflicts over all selected years
        y_out = pd.merge(y_out, y_df.nr_predicted_conflicts.groupby(y_df.ID).sum().to_frame(), on='ID')
        # determine chance of conflict over all selected years
        y_out[column] = y_out.nr_predicted_conflicts / y_out.nr_predictions
    elif column == 'avg_prob_1':
        y_out = pd.merge(y_out, pd.to_numeric(y_df[column]).groupby(y_df.ID).mean().to_frame(), on='ID')
    else:
        raise ValueError('ERROR: column {} is not yet supported'.format(column))
    # add geometry informatoin for each polygon
    y_out = pd.merge(y_out, global_df, on='ID', how='left')

    if not os.path.isdir(os.path.abspath(output_dir)):
        click.echo('creating output folder {}'.format(os.path.abspath(output_dir)))
        os.makedirs(os.path.abspath(output_dir))

    # convert to geo-dataframe
    gdf_out = gpd.GeoDataFrame(y_out, geometry=y_out.geometry)

    # if specified, save as geojson-file to output-dir
    if geojson:
        click.echo('saving to {}'.format(os.path.abspath(os.path.join(output_dir, '{}_merged_{}.geojson'.format(column, suffix)))))
        gdf_out.to_file(os.path.abspath(os.path.join(output_dir, '{}_merged_{}.geojson'.format(column, suffix))), driver='GeoJSON')

    # if specified, save as png-file to output-dir
    if png:
        fig, ax = plt.subplots(1, 1)
        gdf_out.plot(column=column, ax=ax,
                        cmap='Reds',
                        vmin=0, vmax=1,
                        legend=True,
                        legend_kwds={'label': column, 'orientation': "vertical"})

        click.echo('saving to {}'.format(os.path.abspath(os.path.join(output_dir, '{}_merged_{}.png'.format(column, suffix)))))
        plt.savefig(os.path.abspath(os.path.join(output_dir, '{}_merged_{}.png'.format(column, suffix))), dpi=300, bbox_inches='tight')
        plt.close()

if __name__ == '__main__':

    main()

max_issues_repo_name: JannisHoch/copro


text:
import numpy as np

class FeatureEncoder:
    def __init__(self):
        self.active = -1
        self.player_pos_x, self.player_pos_y = 0, 0

    def get_feature_dims(self):
        dims = {
            "player": 29,
            "ball": 18,
            "left_team": 7,
            "left_team_closest": 7,
            "right_team": 7,
            "right_team_closest": 7,
        }
        return dims

    def encode(self, obs):
        player_num = obs["active"]

        player_pos_x, player_pos_y = obs["left_team"][player_num]
        player_direction = np.array(obs["left_team_direction"][player_num])
        player_speed = np.linalg.norm(player_direction)
        player_role = obs["left_team_roles"][player_num]
        player_role_onehot = self._encode_role_onehot(player_role)
        player_tired = obs["left_team_tired_factor"][player_num]
        is_dribbling = obs["sticky_actions"][9]
        is_sprinting = obs["sticky_actions"][8]

        ball_x, ball_y, ball_z = obs["ball"]
        ball_x_relative = ball_x - player_pos_x
        ball_y_relative = ball_y - player_pos_y
        ball_x_speed, ball_y_speed, _ = obs["ball_direction"]
        ball_distance = np.linalg.norm([ball_x_relative, ball_y_relative])
        ball_speed = np.linalg.norm([ball_x_speed, ball_y_speed])
        ball_owned = 0.0
        if obs["ball_owned_team"] == -1:
            ball_owned = 0.0
        else:
            ball_owned = 1.0
        ball_owned_by_us = 0.0
        if obs["ball_owned_team"] == 0:
            ball_owned_by_us = 1.0
        elif obs["ball_owned_team"] == 1:
            ball_owned_by_us = 0.0
        else:
            ball_owned_by_us = 0.0

        ball_which_zone = self._encode_ball_which_zone(ball_x, ball_y)

        if ball_distance > 0.03:
            ball_far = 1.0
        else:
            ball_far = 0.0

        avail = self._get_avail(obs, ball_distance)
        player_state = np.concatenate(
            (
                avail[2:],
                obs["left_team"][player_num],
                player_direction * 100,
                [player_speed * 100],
                player_role_onehot,
                [ball_far, player_tired, is_dribbling, is_sprinting],
            )
        )

        ball_state = np.concatenate(
            (
                np.array(obs["ball"]),
                np.array(ball_which_zone),
                np.array([ball_x_relative, ball_y_relative]),
                np.array(obs["ball_direction"]) * 20,
                np.array(
                    [ball_speed * 20, ball_distance, ball_owned, ball_owned_by_us]
                ),
            )
        )

        obs_left_team = np.delete(obs["left_team"], player_num, axis=0)
        obs_left_team_direction = np.delete(
            obs["left_team_direction"], player_num, axis=0
        )
        left_team_relative = obs_left_team
        left_team_distance = np.linalg.norm(
            left_team_relative - obs["left_team"][player_num], axis=1, keepdims=True
        )
        left_team_speed = np.linalg.norm(obs_left_team_direction, axis=1, keepdims=True)
        left_team_tired = np.delete(
            obs["left_team_tired_factor"], player_num, axis=0
        ).reshape(-1, 1)
        left_team_state = np.concatenate(
            (
                left_team_relative * 2,
                obs_left_team_direction * 100,
                left_team_speed * 100,
                left_team_distance * 2,
                left_team_tired,
            ),
            axis=1,
        )
        left_closest_idx = np.argmin(left_team_distance)
        left_closest_state = left_team_state[left_closest_idx]

        obs_right_team = np.array(obs["right_team"])
        obs_right_team_direction = np.array(obs["right_team_direction"])
        right_team_distance = np.linalg.norm(
            obs_right_team - obs["left_team"][player_num], axis=1, keepdims=True
        )
        right_team_speed = np.linalg.norm(
            obs_right_team_direction, axis=1, keepdims=True
        )
        right_team_tired = np.array(obs["right_team_tired_factor"]).reshape(-1, 1)
        right_team_state = np.concatenate(
            (
                obs_right_team * 2,
                obs_right_team_direction * 100,
                right_team_speed * 100,
                right_team_distance * 2,
                right_team_tired,
            ),
            axis=1,
        )
        right_closest_idx = np.argmin(right_team_distance)
        right_closest_state = right_team_state[right_closest_idx]

        state_dict = {
            "player": player_state,
            "ball": ball_state,
            "left_team": left_team_state,
            "left_closest": left_closest_state,
            "right_team": right_team_state,
            "right_closest": right_closest_state,
            "avail": avail,
        }

        return state_dict

    def _get_avail(self, obs, ball_distance):
        avail = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
        (
            NO_OP,
            MOVE,
            LONG_PASS,
            HIGH_PASS,
            SHORT_PASS,
            SHOT,
            SPRINT,
            RELEASE_MOVE,
            RELEASE_SPRINT,
            SLIDE,
            DRIBBLE,
            RELEASE_DRIBBLE,
        ) = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11)

        ball_x, ball_y, _ = obs["ball"]
        if obs["ball_owned_team"] == 1:  # opponents owning ball
            (
                avail[LONG_PASS],
                avail[HIGH_PASS],
                avail[SHORT_PASS],
                avail[SHOT],
                avail[DRIBBLE],
            ) = (0, 0, 0, 0, 0)
        elif (
            obs["ball_owned_team"] == -1
            and ball_distance > 0.03
            and obs["game_mode"] == 0
        ):  # Ground ball  and far from me
            (
                avail[LONG_PASS],
                avail[HIGH_PASS],
                avail[SHORT_PASS],
                avail[SHOT],
                avail[DRIBBLE],
            ) = (0, 0, 0, 0, 0)
        else:  # my team owning ball
            avail[SLIDE] = 0
            if ball_x > 0.85 and (
                ball_y < -0.34 or ball_y > 0.34
            ):  # when the ball is near the opponent corner
                avail[LONG_PASS], avail[SHORT_PASS], avail[SHOT], avail[DRIBBLE] = (
                    0,
                    0,
                    0,
                    0,
                )

        # Dealing with sticky actions
        sticky_actions = obs["sticky_actions"]
        if sticky_actions[8] == 0:  # sprinting
            avail[RELEASE_SPRINT] = 0

        if sticky_actions[9] == 1:  # dribbling
            avail[SLIDE] = 0
        else:
            avail[RELEASE_DRIBBLE] = 0

        if np.sum(sticky_actions[:8]) == 0:
            avail[RELEASE_MOVE] = 0

        if ball_x < 0.64 or ball_y < -0.27 or 0.27 < ball_y:  # if too far, no shot
            avail[SHOT] = 0
        elif (0.64 <= ball_x and ball_x <= 1.0) and (
            -0.27 <= ball_y and ball_y <= 0.27
        ):  # In the penalty area,  no pass
            avail[HIGH_PASS], avail[LONG_PASS] = 0, 0

        if obs["game_mode"] == 2 and ball_x < -0.7:  # Our GoalKick
            avail = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
            avail[LONG_PASS], avail[HIGH_PASS], avail[SHORT_PASS] = 1, 1, 1
            return np.array(avail)

        elif obs["game_mode"] == 4 and ball_x > 0.9:  # Our CornerKick
            avail = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
            avail[LONG_PASS], avail[HIGH_PASS], avail[SHORT_PASS] = 1, 1, 1
            return np.array(avail)

        elif obs["game_mode"] == 6 and ball_x > 0.6:  # Our PenaltyKick
            avail = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
            avail[SHOT] = 1
            return np.array(avail)

        return np.array(avail)

    def _encode_ball_which_zone(self, ball_x, ball_y):
        MIDDLE_X, PENALTY_X, END_X = 0.2, 0.64, 1.0
        PENALTY_Y, END_Y = 0.27, 0.42
        if (-END_X <= ball_x and ball_x < -PENALTY_X) and (
            -PENALTY_Y < ball_y and ball_y < PENALTY_Y
        ):
            return [1.0, 0, 0, 0, 0, 0]
        elif (-END_X <= ball_x and ball_x < -MIDDLE_X) and (
            -END_Y < ball_y and ball_y < END_Y
        ):
            return [0, 1.0, 0, 0, 0, 0]
        elif (-MIDDLE_X <= ball_x and ball_x <= MIDDLE_X) and (
            -END_Y < ball_y and ball_y < END_Y
        ):
            return [0, 0, 1.0, 0, 0, 0]
        elif (PENALTY_X < ball_x and ball_x <= END_X) and (
            -PENALTY_Y < ball_y and ball_y < PENALTY_Y
        ):
            return [0, 0, 0, 1.0, 0, 0]
        elif (MIDDLE_X < ball_x and ball_x <= END_X) and (
            -END_Y < ball_y and ball_y < END_Y
        ):
            return [0, 0, 0, 0, 1.0, 0]
        else:
            return [0, 0, 0, 0, 0, 1.0]

    def _encode_role_onehot(self, role_num):
        result = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
        result[role_num] = 1.0
        return np.array(result)


text:
from unittest import TestCase
from nose_parameterized import parameterized
from collections import OrderedDict
import os
import gzip

from pandas import (
    Series,
    DataFrame,
    date_range,
    Timestamp,
    read_csv
)
from pandas.util.testing import assert_frame_equal

from numpy import (
    arange,
    zeros_like,
    nan,
)

import warnings

from pyfolio.utils import (to_utc, to_series, check_intraday,
                           detect_intraday, estimate_intraday)
from pyfolio.pos import (get_percent_alloc,
                         extract_pos,
                         get_sector_exposures,
                         get_max_median_position_concentration)

class PositionsTestCase(TestCase):
    dates = date_range(start='2015-01-01', freq='D', periods=20)

    def test_get_percent_alloc(self):
        raw_data = arange(15, dtype=float).reshape(5, 3)
        # Make the first column negative to test absolute magnitudes.
        raw_data[:, 0] *= -1

        frame = DataFrame(
            raw_data,
            index=date_range('01-01-2015', freq='D', periods=5),
            columns=['A', 'B', 'C']
        )

        result = get_percent_alloc(frame)
        expected_raw = zeros_like(raw_data)
        for idx, row in enumerate(raw_data):
            expected_raw[idx] = row / row.sum()

        expected = DataFrame(
            expected_raw,
            index=frame.index,
            columns=frame.columns,
        )

        assert_frame_equal(result, expected)

    def test_extract_pos(self):
        index_dup = [Timestamp('2015-06-08', tz='UTC'),
                     Timestamp('2015-06-08', tz='UTC'),
                     Timestamp('2015-06-09', tz='UTC'),
                     Timestamp('2015-06-09', tz='UTC')]
        index = [Timestamp('2015-06-08', tz='UTC'),
                 Timestamp('2015-06-09', tz='UTC')]

        positions = DataFrame(
            {'amount': [100., 200., 300., 400.],
             'last_sale_price': [10., 20., 30., 40.],
             'sid': [1, 2, 1, 2]},
            index=index_dup
        )
        cash = Series([100., 200.], index=index)

        result = extract_pos(positions, cash)

        expected = DataFrame(OrderedDict([
            (1, [100.*10., 300.*30.]),
            (2, [200.*20., 400.*40.]),
            ('cash', [100., 200.])]),
            index=index
        )
        expected.index.name = 'index'
        expected.columns.name = 'sid'

        assert_frame_equal(result, expected)

    @parameterized.expand([
        (DataFrame([[1.0, 2.0, 3.0, 10.0]]*len(dates),
                   columns=[0, 1, 2, 'cash'], index=dates),
         {0: 'A', 1: 'B', 2: 'A'},
         DataFrame([[4.0, 2.0, 10.0]]*len(dates),
                   columns=['A', 'B', 'cash'], index=dates),
         False),
        (DataFrame([[1.0, 2.0, 3.0, 10.0]]*len(dates),
                   columns=[0, 1, 2, 'cash'], index=dates),
         Series(index=[0, 1, 2], data=['A', 'B', 'A']),
         DataFrame([[4.0, 2.0, 10.0]]*len(dates),
                   columns=['A', 'B', 'cash'], index=dates),
         False),
        (DataFrame([[1.0, 2.0, 3.0, 10.0]]*len(dates),
                   columns=[0, 1, 2, 'cash'], index=dates),
         {0: 'A', 1: 'B'},
         DataFrame([[1.0, 2.0, 10.0]]*len(dates),
                   columns=['A', 'B', 'cash'], index=dates),
         True)
    ])
    def test_sector_exposure(self, positions, mapping,
                             expected_sector_exposure,
                             warning_expected):
        """
        Tests sector exposure mapping and rollup.

        """
        with warnings.catch_warnings(record=True) as w:
            result_sector_exposure = get_sector_exposures(positions,
                                                          mapping)

            assert_frame_equal(result_sector_exposure,
                               expected_sector_exposure)
            if warning_expected:
                self.assertEqual(len(w), 1)
            else:
                self.assertEqual(len(w), 0)

    @parameterized.expand([
        (DataFrame([[1.0, 2.0, 3.0, 14.0]]*len(dates),
                   columns=[0, 1, 2, 'cash'], index=dates),
         DataFrame([[0.15, 0.1, nan, nan]]*len(dates),
                   columns=['max_long', 'median_long',
                            'median_short', 'max_short'], index=dates)),
        (DataFrame([[1.0, -2.0, -13.0, 15.0]]*len(dates),
                   columns=[0, 1, 2, 'cash'], index=dates),
         DataFrame([[1.0, 1.0, -7.5, -13.0]]*len(dates),
                   columns=['max_long', 'median_long',
                            'median_short', 'max_short'], index=dates)),
        (DataFrame([[nan, 2.0, nan, 8.0]]*len(dates),
                   columns=[0, 1, 2, 'cash'], index=dates),
         DataFrame([[0.2, 0.2, nan, nan]]*len(dates),
                   columns=['max_long', 'median_long',
                            'median_short', 'max_short'], index=dates))
    ])
    def test_max_median_exposure(self, positions, expected):
        alloc_summary = get_max_median_position_concentration(positions)
        assert_frame_equal(expected, alloc_summary)

    __location__ = os.path.realpath(
        os.path.join(os.getcwd(), os.path.dirname(__file__)))

    test_returns = read_csv(
        gzip.open(
            __location__ + '/test_data/test_returns.csv.gz'),
        index_col=0, parse_dates=True)
    test_returns = to_series(to_utc(test_returns))
    test_txn = to_utc(read_csv(
        gzip.open(
            __location__ + '/test_data/test_txn.csv.gz'),
        index_col=0, parse_dates=True))
    test_pos = to_utc(read_csv(
        gzip.open(__location__ + '/test_data/test_pos.csv.gz'),
        index_col=0, parse_dates=True))

    @parameterized.expand([
        (test_pos, test_txn, False),
        (test_pos.resample('1W').last(), test_txn, True)
    ])
    def test_detect_intraday(self, positions, transactions, expected):
        detected = detect_intraday(positions, transactions, threshold=0.25)
        assert detected == expected

    @parameterized.expand([
        ('infer', test_returns, test_pos, test_txn, test_pos),
        (False, test_returns, test_pos, test_txn, test_pos)
    ])
    def test_check_intraday(self, estimate, returns,
                            positions, transactions, expected):
        detected = check_intraday(estimate, returns, positions, transactions)
        assert_frame_equal(detected, expected)

    @parameterized.expand([
        (test_returns, test_pos, test_txn, (1506, 8)),
        (test_returns, test_pos.resample('1W').last(), test_txn, (1819, 8))
    ])
    def test_estimate_intraday(self, returns, positions,
                               transactions, expected):
        intraday_pos = estimate_intraday(returns, positions, transactions)
        assert intraday_pos.shape == expected

max_stars_repo_name: tradingalgorithms/pyfolio

max_issues_repo_name: anojangra/pyfolio

max_forks_repo_name: anojangra/pyfolio

reason: from numpy


text:
import cirq
import numpy as np

####### For Hamiltonian ####
class Change_PauliWord_measurement_to_Z_basis(cirq.Gate):
    """
    Class to generate cirq circuit as gate... which generates CNOT entangling gates between non Idenity PauliWord
    qubits in order to perform PauliWord: e^(cofactor * theta * PauliWord)

    e.g.: QubitOperator('X0 Z2 Y3', 0.5)
        gives :
                0: ───H───────────

                3: ───Rx(-0.5π)───

    Args:
        PauliWord_QubitOp (openfermion.QubitOperator): QubitOperator of PauliWord

    Returns
        A cirq circuit object to be used by cirq.Circuit

    """

    def __init__(self, PauliWord_QubitOp):

        self.PauliWord = list(*PauliWord_QubitOp.terms.keys())

    def _decompose_(self, qubits):

        for qubitNo, PauliStr in (self.PauliWord):
            if PauliStr == 'X':
                # yield cirq.H(qubits[qubitNo])
                yield cirq.ry(-np.pi / 2)(qubits[qubitNo])
            elif PauliStr == 'Y':
                yield cirq.rx(np.pi / 2)(qubits[qubitNo])
            elif PauliStr == 'Z' or 'I':
                continue
            else:
                raise ValueError("Qubit Operation: {} is NOT a Pauli operation".format(PauliStr))

    def _circuit_diagram_info_(self, args):
        string_list = []
        for _ in range(self.num_qubits()):
            string_list.append('Changing pauliword to Z basis')
        return string_list

    def num_qubits(self):
        qubit_list, _ = zip(*(self.PauliWord))
        control_qubit = max(qubit_list)
        return control_qubit + 1  # index from 0

class Measure_PauliWord(cirq.Gate):
    """
    Class to generate cirq circuit as gate that measures PauliWord in Z BASIS (all non I terms).

    e.g.: QubitOperator('X0 Z2 Y3', 0.5)
        gives :
                0: ───M───
                      │
                2: ───M───
                      │
                3: ───M───
    Args:
        PauliWord_QubitOp (openfermion.QubitOperator): QubitOperator of PauliWord

    Returns
        A cirq circuit object to be used by cirq.Circuit

    """

    def __init__(self, PauliWord_QubitOp):

        self.PauliWord = list(*PauliWord_QubitOp.terms.keys())

    def _decompose_(self, qubits):

        # list of qubits to measure!
        qubit_list, _ = zip(*(self.PauliWord))

        qubits_to_measure = (qubits[q_No] for q_No in qubit_list)
        if qubit_list != []:
            yield cirq.measure(*qubits_to_measure)
        else:
            return None

    def _circuit_diagram_info_(self, args):
        string_list = []
        for _ in range(self.num_qubits()):
            string_list.append('measuring pauliword')
        return string_list

    def num_qubits(self):
        qubit_list, _ = zip(*(self.PauliWord))
        control_qubit = max(qubit_list)
        return control_qubit + 1  # index from 0

class change_pauliword_to_Z_basis_then_measure(cirq.Gate):
    """
    Class to generate cirq circuit as gate that changes basis and measures PauliWord in Z BASIS (all non I terms).

    e.g.: QubitOperator('X0 Z2 Y3', 0.5)
        gives :
                0: ───H───────────M───
                                  │
                2: ───────────────M───
                                  │
                3: ───Rx(-0.5π)───M───

    Args:
        PauliWord_QubitOp (openfermion.QubitOperator): QubitOperator of PauliWord

    Returns
        A cirq circuit object to be used by cirq.Circuit

    """

    def __init__(self, PauliWord_QubitOp):
        self.PauliWord_QubitOp = PauliWord_QubitOp

    def _decompose_(self, qubits):
        change_basis_to_measure_circuit = Change_PauliWord_measurement_to_Z_basis(self.PauliWord_QubitOp)
        Measure_circuit = Measure_PauliWord(self.PauliWord_QubitOp)

        change_basis_to_measure_circuit_gen = change_basis_to_measure_circuit._decompose_(qubits)
        Measure_circuit_gen = Measure_circuit._decompose_(qubits)

        list_gen = [change_basis_to_measure_circuit_gen, Measure_circuit_gen]
        yield list_gen

    def _circuit_diagram_info_(self, args):
        string_list = []
        for _ in range(self.num_qubits()):
            string_list.append('measuring pauliword: {}'.format(self.PauliWord_QubitOp))
        return string_list

    def num_qubits(self):
        qubit_list, _ = zip(*(list(*self.PauliWord_QubitOp.terms.keys())))
        control_qubit = max(qubit_list)
        return control_qubit + 1  # index from 0

def Generate_Full_Q_Circuit(Full_Ansatz_Q_Circuit, PauliWord_QubitOp):
    """
     Function that appends Ansatz Quantum Circuit to Pauli perform and measure circuit instance.

    Args:
        Full_Ansatz_Q_Circuit (cirq.circuits.circuit.Circuit): Full cirq Ansatz Q circuit
        PauliWord_QubitOp (openfermion.QubitOperator): QubitOperator of PauliWord

    Returns:
        full_circuit (cirq.circuits.circuit.Circuit): Full cirq VQE circuit

    """

    measure_PauliString_in_Z_basis = change_pauliword_to_Z_basis_then_measure(PauliWord_QubitOp)
    measure_PauliString_in_Z_basis_Q_circ = cirq.Circuit(cirq.decompose_once(
        (measure_PauliString_in_Z_basis(*cirq.LineQubit.range(measure_PauliString_in_Z_basis.num_qubits())))))
    full_circuit = cirq.Circuit(
       [
           Full_Ansatz_Q_Circuit.all_operations(),
           *measure_PauliString_in_Z_basis_Q_circ.all_operations(),
       ]
    )
    return full_circuit

def Generate_Full_Q_Circuit_of_Molecular_Hamiltonian(Full_Ansatz_Q_Circuit, QubitOperator_Hamiltonian):
    """
     Function that appends Ansatz Quantum Circuit to Pauli perform and measure circuit instance.

    Args:
        Full_Ansatz_Q_Circuit (cirq.circuits.circuit.Circuit): Full cirq Ansatz Q circuit
        PauliWord_str_list (list): List of tuples containing PauliWord to perform and constant -> (PauliWord, constant)
        n_qubits (int): number of qubits

    Returns:
        dic_holder (dict): Returns a dictionary of each quantum circuit, with cofactor, PauliWord and cirq Q Circuit

    """
    dic_holder = {}
    for index, PauliWord_QubitOp in enumerate(QubitOperator_Hamiltonian):
        for QubitOp_str, const in PauliWord_QubitOp.terms.items():
            temp_d = {}
            if QubitOp_str:
                temp_d['circuit'] = Generate_Full_Q_Circuit(Full_Ansatz_Q_Circuit, PauliWord_QubitOp)
                temp_d['PauliWord'] = PauliWord_QubitOp
            else:
                temp_d['circuit'] = None
                temp_d['PauliWord'] = PauliWord_QubitOp

            dic_holder[index] = temp_d

    return dic_holder

max_stars_repo_name: AlexisRalli/VQE-code

max_issues_repo_name: AlexisRalli/VQE-code

max_forks_repo_name: AlexisRalli/VQE-code


text:
import numpy as np
import time
import cv2
from pygame import mixer
import turtle

def state_machine(sumation,sound):

	# Check if blue color object present in the ROI
	yes = (sumation) > Hatt_thickness[0]*Hatt_thickness[1]*0.8

	# If present play the respective instrument.
	if yes and sound==1:
		drum_clap.play()

	elif yes and sound==2:
		drum_snare.play()
		time.sleep(0.001)
	elif yes and sound == 3:
		drum_bass.play()
		time.sleep(0.1)
	elif yes and sound == 4:
		drum_symbol.play()
		time.sleep(0.1)

def ROI_analysis(frame,sound, color_selected):
	# converting the image into HSV
	hsv = cv2.cvtColor(frame, cv2.COLOR_BGR2HSV)
	# generating mask for
	if color_selected == 'b':
		mask = cv2.inRange(hsv, blueLower, blueUpper)
	elif color_selected == 'g':
		mask = cv2.inRange(hsv, greenLower, greenUpper)
	else:
		mask = cv2.inRange(hsv, redLower, redUpper)

	# Calculating the nuber of white pixels depecting the blue color pixels in the ROI
	sumation = np.sum(mask)
	# Function that decides to play the instrument or not.
	state_machine(sumation,sound)
	return mask

#shows region of intrest when true, shows drums when false
Verbsoe = False

#Import audio files here
#Potentially can test shorter and longer sounds
mixer.init()
drum_clap = mixer.Sound('batterrm.wav')
drum_snare = mixer.Sound('button-2.ogg')
drum_bass = mixer.Sound('BassKick.wav')
drum_symbol = mixer.Sound('symbol_sound.wav')

# HSV range for detecting blue color
#Change these values to test other colors and see result
blueLower = (80, 150, 50)
blueUpper = (120, 255, 255)

redLower = (170, 80, 70)
redUpper = (255, 120, 230)

greenLower = (40, 100, 100)
greenUpper = (80, 255, 255)

#frame accusation from webcam/ usb camera
camera = cv2.VideoCapture(0)
ret, frame = camera.read()

#Shows the images of the drums on the screen
# Also can change this for testing
#to read an image the syntax is cv2.imread("path of the image") to reshape the image as per
#desired shape we use cv2.resize(image,(width,heigh),interpolation=INTERPOLATION_METHOD_FLAG_VALUE)
#The size for augmenting the objects is decided based on the ROI.

Hatt = cv2.resize(cv2.imread('Hatt.png'), (200,100), interpolation=cv2.INTER_CUBIC)
Snare = cv2.resize(cv2.imread('Snare.png'), (200,100), interpolation=cv2.INTER_CUBIC)
Symbol = cv2.resize(cv2.imread('cymbal.png'), (200,100), interpolation=cv2.INTER_CUBIC)
Bass = cv2.resize(cv2.imread('BDrum.png'), (200,100), interpolation=cv2.INTER_CUBIC)

# Setting the corners of ROI for blue colour detection
#To detect blue colour we need to perform certain operations on each captured frame.
# These operations need some computations to be performed by the processor.
# Since our instruments are fixed in this application and we want to play the sound only if
# the blue colour object hits the instrument (detected inside the ROI) it is a good idea to
# perform all these operations only inside the ROI.

#4 drum set up bottom 2
Hatt_center_4 = [np.shape(frame)[1]*3//9,np.shape(frame)[0]*8//9]
Snare_center_4 = [np.shape(frame)[1]*6//9,np.shape(frame)[0]*8//9]
Bass_center = [np.shape(frame)[1]*2//10,np.shape(frame)[0]*3//8]
Symbol_center = [np.shape(frame)[1]*8//10,np.shape(frame)[0]*3//8]

# 3 Drum set up
Hatt_center_3 = [np.shape(frame)[1]*8//10,np.shape(frame)[0]*4//8]
Snare_center_3 = [np.shape(frame)[1]*2//10,np.shape(frame)[0]*4//8]
Bass_center_3 = [np.shape(frame)[1]*5//10,np.shape(frame)[0]*8//9]

# 2 drum set up
Hatt_center = [np.shape(frame)[1]*2//8,np.shape(frame)[0]*6//8]
Snare_center = [np.shape(frame)[1]*6//8,np.shape(frame)[0]*6//8]

# DRUMSET OF 4
Bass_thickness = [200,100]
Bass_top = [Bass_center[0]-Bass_thickness[0]//2,Bass_center[1]-Bass_thickness[1]//2]
Bass_btm = [Bass_center[0]+Bass_thickness[0]//2,Bass_center[1]+Bass_thickness[1]//2]

Symbol_thickness = [200,100]
Symbol_top = [Symbol_center[0]-Symbol_thickness[0]//2,Symbol_center[1]-Symbol_thickness[1]//2]
Symbol_btm = [Symbol_center[0]+Symbol_thickness[0]//2,Symbol_center[1]+Symbol_thickness[1]//2]

Hatt_thickness = [200,100]
Hatt_top_4 = [Hatt_center_4[0]-Hatt_thickness[0]//2,Hatt_center_4[1]-Hatt_thickness[1]//2]
Hatt_btm_4 = [Hatt_center_4[0]+Hatt_thickness[0]//2,Hatt_center_4[1]+Hatt_thickness[1]//2]

Snare_thickness = [200,100]
Snare_top_4 = [Snare_center_4[0]-Snare_thickness[0]//2,Snare_center_4[1]-Snare_thickness[1]//2]
Snare_btm_4 = [Snare_center_4[0]+Snare_thickness[0]//2,Snare_center_4[1]+Snare_thickness[1]//2]

# Drumset of 3
Hatt_top_3 = [Hatt_center_3[0]-Hatt_thickness[0]//2,Hatt_center_3[1]-Hatt_thickness[1]//2]
Hatt_btm_3 = [Hatt_center_3[0]+Hatt_thickness[0]//2,Hatt_center_3[1]+Hatt_thickness[1]//2]

Snare_top_3 = [Snare_center_3[0]-Snare_thickness[0]//2,Snare_center_3[1]-Snare_thickness[1]//2]
Snare_btm_3 = [Snare_center_3[0]+Snare_thickness[0]//2,Snare_center_3[1]+Snare_thickness[1]//2]

Bass_top_3 = [Bass_center_3[0]-Bass_thickness[0]//2,Bass_center_3[1]-Bass_thickness[1]//2]
Bass_btm_3 = [Bass_center_3[0]+Bass_thickness[0]//2,Bass_center_3[1]+Bass_thickness[1]//2]
#DRUMSET OF 2

Hatt_top = [Hatt_center[0]-Hatt_thickness[0]//2,Hatt_center[1]-Hatt_thickness[1]//2]
Hatt_btm = [Hatt_center[0]+Hatt_thickness[0]//2,Hatt_center[1]+Hatt_thickness[1]//2]

Snare_top = [Snare_center[0]-Snare_thickness[0]//2,Snare_center[1]-Snare_thickness[1]//2]
Snare_btm = [Snare_center[0]+Snare_thickness[0]//2,Snare_center[1]+Snare_thickness[1]//2]

time.sleep(1)
color = "b"
numDrums = 2
print("Press 'b' for blue drumsticks, 'r' for red, or 'g' for green. Press '2' for two drums, '3' for three drums, and '4' for four drums. When you are done press 'q' to quit.")

while True:

	# grab the current frame
	ret, frame = camera.read()
	frame = cv2.flip(frame,1)

	if not(ret):
		break
	# Selecting ROI corresponding to drums selecteed
	if numDrums ==2:
		snare_ROI = np.copy(frame[Snare_top[1]:Snare_btm[1],Snare_top[0]:Snare_btm[0]])
		hatt_ROI = np.copy(frame[Hatt_top[1]:Hatt_btm[1],Hatt_top[0]:Hatt_btm[0]])

	elif numDrums == 3:
		snare_ROI = np.copy(frame[Snare_top_3[1]:Snare_btm_3[1],Snare_top_3[0]:Snare_btm_3[0]])
		hatt_ROI = np.copy(frame[Hatt_top_3[1]:Hatt_btm_3[1],Hatt_top_3[0]:Hatt_btm_3[0]])
		bass_ROI = np.copy(frame[Bass_top_3[1]:Bass_btm_3[1],Bass_top_3[0]:Bass_btm_3[0]])

	elif numDrums == 4:
		snare_ROI = np.copy(frame[Snare_top_4[1]:Snare_btm_4[1],Snare_top_4[0]:Snare_btm_4[0]])
		hatt_ROI = np.copy(frame[Hatt_top_4[1]:Hatt_btm_4[1],Hatt_top_4[0]:Hatt_btm_4[0]])
		bass_ROI = np.copy(frame[Bass_top[1]:Bass_btm[1],Bass_top[0]:Bass_btm[0]])
		symbol_ROI = np.copy(frame[Symbol_top[1]:Symbol_btm[1],Symbol_top[0]:Symbol_btm[0]])
		# Selecting ROI corresponding to Hatt

	mask = ROI_analysis(snare_ROI,1, color)
	mask = ROI_analysis(hatt_ROI,2,color)

 	# Selecting ROI corresponding to bass
	if numDrums == 4:
		mask = ROI_analysis(bass_ROI,3,color)
		mask = ROI_analysis(symbol_ROI,4,color)
	elif numDrums == 3:
		mask = ROI_analysis(bass_ROI,3,color)

	# A writing text on an image.
	cv2.putText(frame,'Project: AR DrumKit',(10,30),2,1,(20,20,20),2)

	# Display the ROI to view the blue colour being detected
	if Verbsoe:
		# Displaying the ROI in the Image
		if numDrums == 2:
			frame[Snare_top[1]:Snare_btm[1],Snare_top[0]:Snare_btm[0]] = cv2.bitwise_and(frame[Snare_top[1]:Snare_btm[1],Snare_top[0]:Snare_btm[0]],frame[Snare_top[1]:Snare_btm[1],Snare_top[0]:Snare_btm[0]], mask=mask[Snare_top[1]:Snare_btm[1],Snare_top[0]:Snare_btm[0]])
			frame[Hatt_top[1]:Hatt_btm[1],Hatt_top[0]:Hatt_btm[0]] = cv2.bitwise_and(frame[Hatt_top[1]:Hatt_btm[1],Hatt_top[0]:Hatt_btm[0]],frame[Hatt_top[1]:Hatt_btm[1],Hatt_top[0]:Hatt_btm[0]],mask=mask[Hatt_top[1]:Hatt_btm[1],Hatt_top[0]:Hatt_btm[0]])
		#Allows for 4 drums
		elif numDrums == 4:
			frame[Snare_top_4[1]:Snare_btm_4[1],Snare_top_4[0]:Snare_btm_4[0]] = cv2.bitwise_and(frame[Snare_top_4[1]:Snare_btm_4[1],Snare_top_4[0]:Snare_btm_4[0]],frame[Snare_top_4[1]:Snare_btm_4[1],Snare_top_4[0]:Snare_btm_4[0]], mask=mask[Snare_top_4[1]:Snare_btm_4[1],Snare_top_4[0]:Snare_btm_4[0]])
			frame[Hatt_top_4[1]:Hatt_btm_4[1],Hatt_top_4[0]:Hatt_btm_4[0]] = cv2.bitwise_and(frame[Hatt_top_4[1]:Hatt_btm_4[1],Hatt_top_4[0]:Hatt_btm_4[0]],frame[Hatt_top_4[1]:Hatt_btm_4[1],Hatt_top_4[0]:Hatt_btm_4[0]],mask=mask[Hatt_top_4[1]:Hatt_btm_4[1],Hatt_top_4[0]:Hatt_btm_4[0]])
			frame[Bass_top[1]:Bass_btm[1],Bass_top[0]:Bass_btm[0]] = cv2.bitwise_and(frame[Bass_top[1]:Bass_btm[1],Bass_top[0]:Bass_btm[0]],frame[Bass_top[1]:Bass_btm[1],Bass_top[0]:Bass_btm[0]], mask=mask[Bass_top[1]:Bass_btm[1],Bass_top[0]:Bass_btm[0]])
			frame[Symbol_top[1]:Symbol_btm[1],Symbol_top[0]:Symbol_btm[0]] = cv2.bitwise_and(frame[Symbol_top[1]:Symbol_btm[1],Symbol_top[0]:Symbol_btm[0]],frame[Symbol_top[1]:Symbol_btm[1],Symbol_top[0]:Symbol_btm[0]],mask=mask[Symbol_top[1]:Symbol_btm[1],Symbol_top[0]:Symbol_btm[0]])
		elif numDrums == 3:
			frame[Snare_top_3[1]:Snare_btm_3[1],Snare_top_3[0]:Snare_btm_3[0]] = cv2.bitwise_and(frame[Snare_top_3[1]:Snare_btm_3[1],Snare_top_3[0]:Snare_btm_3[0]],frame[Snare_top_3[1]:Snare_btm_3[1],Snare_top_3[0]:Snare_btm_3[0]], mask=mask[Snare_top_3[1]:Snare_btm_3[1],Snare_top_3[0]:Snare_btm_3[0]])
			frame[Hatt_top_3[1]:Hatt_btm_3[1],Hatt_top_3[0]:Hatt_btm_3[0]] = cv2.bitwise_and(frame[Hatt_top_3[1]:Hatt_btm_3[1],Hatt_top_3[0]:Hatt_btm_3[0]],frame[Hatt_top_3[1]:Hatt_btm_3[1],Hatt_top_3[0]:Hatt_btm_3[0]],mask=mask[Hatt_top_3[1]:Hatt_btm_3[1],Hatt_top_3[0]:Hatt_btm_3[0]])
			frame[Bass_top_3[1]:Bass_btm_3[1],Bass_top_3[0]:Bass_btm_3[0]] = cv2.bitwise_and(frame[Bass_top_3[1]:Bass_btm_3[1],Bass_top_3[0]:Bass_btm_3[0]],frame[Bass_top_3[1]:Bass_btm_3[1],Bass_top_3[0]:Bass_btm_3[0]], mask=mask[Bass_top_3[1]:Bass_btm_3[1],Bass_top_3[0]:Bass_btm_3[0]])
	# Augmenting the instruments in the output frame.
	else:
		# Augmenting the image of the instruments on the frame.
		if numDrums == 2:
			frame[Snare_top[1]:Snare_btm[1],Snare_top[0]:Snare_btm[0]] = cv2.addWeighted(Snare, 1, frame[Snare_top[1]:Snare_btm[1],Snare_top[0]:Snare_btm[0]], 1, 0)
			frame[Hatt_top[1]:Hatt_btm[1],Hatt_top[0]:Hatt_btm[0]] = cv2.addWeighted(Hatt, 1, frame[Hatt_top[1]:Hatt_btm[1],Hatt_top[0]:Hatt_btm[0]], 1, 0)
		#allows for 4 drums
		elif numDrums == 4:
			frame[Snare_top_4[1]:Snare_btm_4[1],Snare_top_4[0]:Snare_btm_4[0]] = cv2.addWeighted(Snare, 1, frame[Snare_top_4[1]:Snare_btm_4[1],Snare_top_4[0]:Snare_btm_4[0]], 1, 0)
			frame[Hatt_top_4[1]:Hatt_btm_4[1],Hatt_top_4[0]:Hatt_btm_4[0]] = cv2.addWeighted(Hatt, 1, frame[Hatt_top_4[1]:Hatt_btm_4[1],Hatt_top_4[0]:Hatt_btm_4[0]], 1, 0)
			frame[Bass_top[1]:Bass_btm[1],Bass_top[0]:Bass_btm[0]] = cv2.addWeighted(Bass, 1, frame[Bass_top[1]:Bass_btm[1],Bass_top[0]:Bass_btm[0]], 1, 0)
			frame[Symbol_top[1]:Symbol_btm[1],Symbol_top[0]:Symbol_btm[0]] = cv2.addWeighted(Symbol, 1, frame[Symbol_top[1]:Symbol_btm[1],Symbol_top[0]:Symbol_btm[0]], 1, 0)
		elif numDrums == 3:
			frame[Snare_top_3[1]:Snare_btm_3[1],Snare_top_3[0]:Snare_btm_3[0]] = cv2.addWeighted(Snare, 1, frame[Snare_top_3[1]:Snare_btm_3[1],Snare_top_3[0]:Snare_btm_3[0]], 1, 0)
			frame[Hatt_top_3[1]:Hatt_btm_3[1],Hatt_top_3[0]:Hatt_btm_3[0]] = cv2.addWeighted(Hatt, 1, frame[Hatt_top_3[1]:Hatt_btm_3[1],Hatt_top_3[0]:Hatt_btm_3[0]], 1, 0)
			frame[Bass_top_3[1]:Bass_btm_3[1],Bass_top_3[0]:Bass_btm_3[0]] = cv2.addWeighted(Bass, 1, frame[Bass_top_3[1]:Bass_btm_3[1],Bass_top_3[0]:Bass_btm_3[0]], 1, 0)

	cv2.imshow('Output',frame)
	key = cv2.waitKey(1) & 0xFF
	# if the 'q' key is pressed, stop the loop
	if key == ord("q"):
		break
	elif key == ord("b"):
		color = "b"
	elif key == ord("g"):
		color = "g"
	elif key == ord("r"):
		color = "r"
	elif key == ord("2"):
		numDrums = 2
	elif key == ord("3"):
		numDrums = 3
	elif key == ord("4"):
		numDrums = 4

# cleanup the camera and close any open windows
camera.release()
cv2.destroyAllWindows()

max_stars_repo_name: csu-hci-projects/AR-Drum-Kit

max_issues_repo_name: csu-hci-projects/AR-Drum-Kit

max_forks_repo_name: csu-hci-projects/AR-Drum-Kit