text:
"""Based on LipNet cropping of the mouth ROI"""
import sys

if sys.version_info[0] < 3 and sys.version_info[1] < 2:
    raise Exception("Must be using >= Python 3.2")

from os import listdir, path

if not path.isfile('face_detection/detection/sfd/s3fd.pth'):
    raise FileNotFoundError('Save the s3fd model to face_detection/sfd/s3fd.pth \
							before running this script!')

import multiprocessing as mp
from concurrent.futures import ThreadPoolExecutor, as_completed
import numpy as np
import argparse, os, cv2, traceback, subprocess
from tqdm import tqdm
from glob import glob
from synthesizer import audio
from synthesizer.hparams import hparams as hp

import face_alignment
# import face_detection
import dlib
import shutil
import torch
import time

# from video_inference import VIDEO_TO_AUDIO_COMMAND

# fa = [face_alignment.FaceAlignment(face_alignment.LandmarksType._2D, flip_input=False,
#                                    device='cpu') for _ in range(args.ngpu)]  # S3FD face detectors CPU
# fa = [face_alignment.FaceAlignment(face_alignment.LandmarksType._2D, flip_input=False,
#                                    device='cuda:{}'.format(id)) for id in range(args.ngpu)]  # S3FD face detectors multiple GPU
fa = face_alignment.FaceAlignment(face_alignment.LandmarksType._2D, flip_input=False, device='cuda:{}'.format(0))  # S3FD face detectors 1 GPU

# fa = [face_detection.FaceAlignment(face_detection.LandmarksType._2D, flip_input=False,
#                                    device='cuda:{}'.format(id)) for id in range(args.ngpu)]  # S3FD face detectors GPU

# face_detectors = [dlib.get_frontal_face_detector() for _ in range(args.ngpu)]  # CPU HOG Face detectors
# shape_predictors = [dlib.shape_predictor('shape_predictor_68_face_landmarks.dat') for _ in range(args.ngpu)]  # CPU shape predictors based on landmarks
face_detector = dlib.get_frontal_face_detector()
shape_predictor = dlib.shape_predictor('shape_predictor_68_face_landmarks.dat')

"""3 ways of face and landmark detection:
1: DLIB detector (CPU) and DLIB landmarks (CPU) - very fast w/ multiple cores
2: S3FD face detector (GPU) and DLIB landmarks (CPU) - 1.3 videos per second, not as fast as 1
3: S3FD face detector (GPU) and Pytorch landmarks (GPU) - ~20FPS
"""

def get_frames_mouth(detector, predictor, frames):
    """
    1: DLIB face and landmark detector
    """
    MOUTH_WIDTH = 100
    MOUTH_HEIGHT = 50
    HORIZONTAL_PAD = 0.19
    normalize_ratio = None
    mouth_frames = []
    for frame in frames:
        dets = detector(frame, 1)
        # pred = detector.get_detections_for_batch(np.array([frame]))[0]
        # dets = [dlib.rectangle(*pred[1:])]
        shape = None
        for k, d in enumerate(dets):
            shape = predictor(frame, d)
            i = -1
        if shape is None: # Detector doesn't detect face, just return as is
            return frames
        mouth_points = []
        for part in shape.parts():
            i += 1
            if i < 48: # Only take mouth region
                continue
            mouth_points.append((part.x,part.y))
        np_mouth_points = np.array(mouth_points)

        mouth_centroid = np.mean(np_mouth_points[:, -2:], axis=0)

        if normalize_ratio is None:
            mouth_left = np.min(np_mouth_points[:, :-1]) * (1.0 - HORIZONTAL_PAD)
            mouth_right = np.max(np_mouth_points[:, :-1]) * (1.0 + HORIZONTAL_PAD)

            normalize_ratio = MOUTH_WIDTH / float(mouth_right - mouth_left)

        new_img_shape = (int(frame.shape[0] * normalize_ratio), int(frame.shape[1] * normalize_ratio))
        # resized_img = imresize(frame, new_img_shape)
        resized_img = cv2.resize(frame, new_img_shape)

        mouth_centroid_norm = mouth_centroid * normalize_ratio

        mouth_l = int(mouth_centroid_norm[0] - MOUTH_WIDTH / 2)
        mouth_r = int(mouth_centroid_norm[0] + MOUTH_WIDTH / 2)
        mouth_t = int(mouth_centroid_norm[1] - MOUTH_HEIGHT / 2)
        mouth_b = int(mouth_centroid_norm[1] + MOUTH_HEIGHT / 2)

        mouth_crop_image = resized_img[mouth_t:mouth_b, mouth_l:mouth_r]
        if mouth_crop_image.shape != (50, 100, 3):
            mouth_crop_image = cv2.resize(mouth_crop_image, (50, 100))

        mouth_frames.append(mouth_crop_image)
    return mouth_frames

def s3fd_detector_and_pytorch_landmarks(_frames):
        """
        # 3: S3FD detector and Pytorch landmarks
        # this works but is quite slow because only 1 GPU
        # takes approx 1 second per video
        # takes longer on CPU OMG!!!
        """
        all_landmarks = []
        for frame in _frames:
            landmarks = fa.get_landmarks_from_image(np.asarray(frame))
            if landmarks is not None:
                all_landmarks.append(landmarks[0])

        if len(_frames) != len(all_landmarks):
            return []

        MOUTH_WIDTH = 100
        MOUTH_HEIGHT = 50
        HORIZONTAL_PAD = 0.19
        normalize_ratio = None
        mouth_frames = []
        for frame, landmarks in zip(_frames, all_landmarks):
            mouth_points = []
            for i in range(len(landmarks)):
                if i < 48:
                    continue
                mouth_points.append((landmarks[i][0], landmarks[i][1]))

            np_mouth_points = np.array(mouth_points)
            mouth_centroid = np.mean(np_mouth_points[:, -2:], axis=0)

            if normalize_ratio is None:
                mouth_left = np.min(np_mouth_points[:, :-1]) * (1.0 - HORIZONTAL_PAD)
                mouth_right = np.max(np_mouth_points[:, :-1]) * (1.0 + HORIZONTAL_PAD)

                normalize_ratio = MOUTH_WIDTH / float(mouth_right - mouth_left)

            new_img_shape = (int(frame.shape[0] * normalize_ratio), int(frame.shape[1] * normalize_ratio))
            # resized_img = imresize(frame, new_img_shape)
            resized_img = cv2.resize(frame, new_img_shape)

            mouth_centroid_norm = mouth_centroid * normalize_ratio

            mouth_l = int(mouth_centroid_norm[0] - MOUTH_WIDTH / 2)
            mouth_r = int(mouth_centroid_norm[0] + MOUTH_WIDTH / 2)
            mouth_t = int(mouth_centroid_norm[1] - MOUTH_HEIGHT / 2)
            mouth_b = int(mouth_centroid_norm[1] + MOUTH_HEIGHT / 2)

            mouth_crop_image = resized_img[mouth_t:mouth_b, mouth_l:mouth_r]
            if mouth_crop_image.shape != (50, 100, 3):
                mouth_crop_image = cv2.resize(mouth_crop_image, (50, 100))
            mouth_frames.append(mouth_crop_image)

        return mouth_frames

def mouth_frames_bad(_mouth_frames):
    return len(_mouth_frames) == 0 or not all([mouth_frame.shape == (50, 100, 3) for mouth_frame in _mouth_frames])

def process_video_file(vfile, args, gpu_id, detector='dlib'):
    video_stream = cv2.VideoCapture(vfile)

    frames = []
    while 1:
        still_reading, frame = video_stream.read()
        if not still_reading:
            video_stream.release()
            break
        frames.append(frame)

    vidname = os.path.basename(vfile).split('.')[0]
    split = vfile.split('/')[-2]
    word = vfile.split('/')[-3]

    fulldir = path.join(args.preprocessed_root, word, split, vidname)
    os.makedirs(fulldir, exist_ok=True)
    #print (fulldir)

    wavpath = path.join(fulldir, 'audio.wav')
    specpath = path.join(fulldir, 'mels.npz')

    command = VIDEO_TO_AUDIO_COMMAND.format(
        input_video_path=vfile,
        output_audio_path=wavpath
    )
    subprocess.call(command, shell=True, stdout=subprocess.DEVNULL, stderr=subprocess.STDOUT)

    # this works and is super fast because of multiple CPU cores
    if detector == 'dlib':
        mouth_frames = get_frames_mouth(face_detector, shape_predictor, frames)
    else:
        mouth_frames = s3fd_detector_and_pytorch_landmarks(frames)

    if mouth_frames_bad(mouth_frames):
        shutil.rmtree(fulldir)
        return False

    for i, mouth_frame in enumerate(mouth_frames):
        cv2.imwrite(path.join(fulldir, '{}.jpg'.format(i)), mouth_frame)

    return True

def process_audio_file(vfile, args, gpu_id):
    vidname = os.path.basename(vfile).split('.')[0]
    split = vfile.split('/')[-2]
    word = vfile.split('/')[-3]

    fulldir = path.join(args.preprocessed_root, word, split, vidname)
    os.makedirs(fulldir, exist_ok=True)

    wavpath = path.join(fulldir, 'audio.wav')
    specpath = path.join(fulldir, 'mels.npz')

    wav = audio.load_wav(wavpath, hp.sample_rate)
    spec = audio.melspectrogram(wav, hp)
    lspec = audio.linearspectrogram(wav, hp)
    np.savez_compressed(specpath, spec=spec, lspec=lspec)

def mp_handler(job):
    vfile, args, gpu_id = job
    try:
        success = process_video_file(vfile, args, gpu_id)
        if success:
            process_audio_file(vfile, args, gpu_id)
    except Exception as e:
        print(e, flush=True)

def process(process_id, jobs):
    failed_videos = []
    for vfile in tqdm(jobs):
        try:
            success = process_video_file(vfile, args, process_id, detector='dlib')
            if success:
                process_audio_file(vfile, args, process_id)
            else:
                failed_videos.append(vfile)
        except Exception as e:
            print(e, flush=True)
            failed_videos.append(vfile)

    return failed_videos

def main(args):
    print('Started processing for {} with {} GPUs'.format(args.data_root, args.ngpu))

    filelist = list(glob(path.join(args.data_root, '*/{}/*.mp4'.format(args.split))))

    """
    jobs = [(vfile, args, i%args.ngpu) for i, vfile in enumerate(filelist)]
    p = ThreadPoolExecutor(args.ngpu)
    futures = [p.submit(mp_handler, j) for j in jobs]
    _ = [r.result() for r in tqdm(as_completed(futures), total=len(futures))]
    """

    """
    for i, vfile in tqdm(enumerate(filelist)):
        mp_handler((vfile, args, i%args.ngpu))
    """

    """
    # run dlib face and landmark detectors first on majority
    num_processes = args.ngpu
    num_tasks_per_process = len(filelist) // num_processes
    jobs = []
    for i in range(num_processes):
        start = i * num_tasks_per_process
        end = start + num_tasks_per_process
        jobs.append([i+1, filelist[start:end]])
    all_failed_videos = []
    with mp.Pool(processes=num_processes) as p:
        all_results = p.starmap(process, jobs)
        for failed_videos in all_results:
            all_failed_videos.extend(failed_videos)
    """

    # quickly grab the failed DLIB videos
    from pathlib import Path
    all_failed_videos = []
    for vfile in filelist:
        vidname = os.path.basename(vfile).split('.')[0]
        split = vfile.split('/')[-2]
        word = vfile.split('/')[-3]
        processed_path = Path(args.preprocessed_root).joinpath(word).joinpath(split).joinpath(vidname)
        if not processed_path.exists():
            all_failed_videos.append(vfile)

    print('Num failed so far:', len(all_failed_videos))

    # try GPU face/landmark detector on failed
    failed_videos = []
    final_fail_counter = 0
    for vfile in tqdm(all_failed_videos):
        try:
            success = process_video_file(vfile, args, 0, detector='gpu')
            if success:
                process_audio_file(vfile, args, 0)
            else:
                final_fail_counter += 1
        except Exception as e:
            print(e, flush=True)
            final_fail_counter += 1
    print('Num failed final:', final_fail_counter)

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

    parser.add_argument('--ngpu', help='Number of GPUs across which to run in parallel', default=1, type=int)
    parser.add_argument('--batch_size', help='Single GPU Face detection batch size', default=32, type=int)
    parser.add_argument("--data_root", help="Root folder of the LRW dataset", required=True)
    parser.add_argument("--preprocessed_root", help="Root folder of the preprocessed dataset", required=True)
    parser.add_argument("--split", help="* | train | val | test (* will preprocess all splits)", default="test",
                        choices=["*", "train", "val", "test"])

    args = parser.parse_args()

    main(args)


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

import numpy as np
from GT3.Core.Functions.DrawCoreLine import draw_core_line
from scipy.interpolate import UnivariateSpline
from shapely.geometry import LineString
from math import pi, sqrt

def create_surf_area_interp(rho2psinorm, psinorm2rho, psi_data, sep_pts, R0_a, a, sep_val):
    rho_vals = np.linspace(0, sep_val, 50, endpoint=True)
    r_vals = rho_vals * a
    psi_vals = rho2psinorm(rho_vals)

    surf_area = np.zeros(rho_vals.shape)
    # skip the first value. First value is zero.
    for i, psi_val in enumerate(psi_vals[1:]):
        if psi_val < 1:
            try:
                # try to get the actual length of the fs trace
                fs_line, fs_pts = draw_core_line(psi_data.R, psi_data.Z, psi_data.psi_norm, psi_val, sep_pts)
                surf_area[i + 1] = fs_line.length * 2 * pi * fs_pts.axis[0]
            except:
                # if too close to the seperatrix, it may return None. In this case, estimate perimeter using Ramanujan
                # TODO: Improve this estimation. Also, see comments in draw_core_line
                rho_val = psinorm2rho(psi_val)
                a_el = rho_val * a  # a of the ellipse
                b_el = rho_val * a * 1.5  # b of the ellipse, assumed kappa of 1.5. This needs to be improved

                surf_area[i + 1] = pi * (3*(a_el+b_el) - sqrt((3*a_el + b_el)*(a_el+3*b_el)))

        else:
            surf_area[i+1] = LineString(sep_pts).length * 2*pi*R0_a

    # For some reason, psi_vals will sometimes not be monotonically increasing, especially near the magnetic axis.
    # This prevents UnivariateSpline from working. To prevent this, we're just going to delete any points that are
    # lower than the previous point, prior to doing the fit. As long as the psi data isn't too wonky, this should be
    # fine.
    psinormvals_mi = []  # psivals that are monotonically increasing
    rvals_mi = []  # rvals corresponding to monotonically increasing psi
    rhovals_mi = []  # rhovals corresponding to monotonically increasing psi
    surfarea_mi = []  # surf_area corresponding to monotonically increasing psi
    for i, psi_val in enumerate(psi_vals):
        if i == 0:
            rvals_mi.append(0)
            rhovals_mi.append(0)
            psinormvals_mi.append(0)
            surfarea_mi.append(0)
        elif psi_val > psinormvals_mi[-1]:
            rvals_mi.append(r_vals[i])
            rhovals_mi.append(rho_vals[i])
            psinormvals_mi.append(psi_vals[i])
            surfarea_mi.append(surf_area[i])

    rvals_mi = np.asarray(rvals_mi)
    rhovals_mi = np.asarray(rhovals_mi)
    psinormvals_mi = np.asarray(psinormvals_mi)
    surfarea_mi = np.asarray(surfarea_mi)

    r2sa = UnivariateSpline(rvals_mi, surfarea_mi, k=3, s=0)
    rho2sa = UnivariateSpline(rhovals_mi, surfarea_mi, k=3, s=0)
    psinorm2sa = UnivariateSpline(psinormvals_mi, surfarea_mi, k=3, s=0)
    # r2sa = UnivariateSpline(r_vals, surf_area, k=3, s=0)
    # rho2sa = UnivariateSpline(rho_vals, surf_area, k=3, s=0)
    # psinorm2sa = UnivariateSpline(psi_vals, surf_area, k=3, s=0)

    return r2sa, rho2sa, psinorm2sa

ext: py

reason: import numpy,from scipy


text:
################################
########### Imports ############
################################

import numpy as np                                              # Basic operators
import matplotlib.pyplot as plt                                 # Plotting
import os                                                       # Listing filenames

from scipy.interpolate import InterpolatedUnivariateSpline      # Interpolation (spline) function
import lmfit as lm                                              # Fitting

from ipywidgets import interactive, fixed, FloatSlider, HBox, Layout, Button, Label, Output, VBox
from IPython.display import display

import warnings
warnings.filterwarnings("ignore")         # ignore warnings

################################
############ Data ##############
################################
# Import data from previously downloaded sparc website: http://astroweb.cwru.edu/SPARC/
# Lelli, Federico, Stacy S. McGaugh, and James M. Schombert. “SPARC: Mass Models for 175 Disk Galaxies with Spitzer Photometry and Accurate Rotation Curves.” The Astronomical Journal 152, no. 6 (2016): 157. https://doi.org/10.3847/0004-6256/152/6/157.

# Read the file where the chosen galaxy name is located (this was modified by the user)
textfile = open('python/chosengalaxy.txt', 'r')
galaxy = textfile.read()
textfile.close()

# Define SPARC directory
SPARC_file_directory='./data/sparc/'                       #note that '' means the string variable is blank

# Define file path for .dat files
SPARC_file_path = SPARC_file_directory + galaxy + '_rotmod.dat'

# Load the galaxy data
data = np.loadtxt(SPARC_file_path)

# Split columns into arrays
Rad,Vobs,errV,Vgas,Vdisk,Vbul,SBdisk,SBbul = data.T

# Check if it has all components
if np.sum(Vbul) == 0:     # Bulge
    nobulge = True
    warning_bulge = "There is no bulge component."
else:
    nobulge = False
    warning_bulge = ""

if np.sum(Vdisk) == 0:    # Disk
    nodisk = True
    warning_disk = "There is no disk component."
else:
    nodisk = False
    warning_disk = ""

if np.sum(Vgas) == 0:     # Gas
    nogas = True
    warning_gas = "There is no gas component."
else:
    nogas = False
    warning_gas = ""

# Define distance to galaxy in Mpc
firstline = open(SPARC_file_path).readline()
firstline = firstline.split()
distance = float(firstline[3])

##############################
### Interpolation function ###
##############################

def interpd(x,y):
    return InterpolatedUnivariateSpline(x,y,k=3)

################################
######### Components ###########
################################

# Bulge
def bulge(r,bpref):
    polynomial = interpd(Rad,bpref*Vbul)
    return polynomial(r)

# Disk
def disk(r,dpref):
    polynomial = interpd(Rad,dpref*Vdisk)
    return polynomial(r)

# Gas
def gas(r):
    polynomial = interpd(Rad,Vgas)
    return polynomial(r)

# Halo
## Set parameters
rho0 = 3.10e8       # Central mass density (in solar mass/kpc^3)
rc = 1.4            # Core radius (in kpc)
G = 4.300e-6        # Gravitational constant (kpc/solar mass*(km/s)^2)

## Equation for dark matter halo velocity
def halo(r,rc,rho0):
    v = np.sqrt(4*np.pi*G*rho0*rc**2*(1 - rc/r * np.arctan(r/rc)))
    return v

# Total velocity
def totalcurve(r,bpref,dpref,rc,rho0):
    total = np.sqrt(bulge(r,bpref)**2
                    + disk(r,dpref)**2
                    + halo(r,rc,rho0)**2
                    + gas(r)**2)
    return total

#################################
#### Find Fitting Parameters ####
#################################

# Setup
fit_mod = lm.Model(totalcurve)
fit_params = fit_mod.make_params()
weighdata = 1/errV

# Halo
fit_params.add('rc',    value=rc,    min=0.1)           # Core Radius (kpc)
fit_params.add('rho0',  value=rho0,  min=0)             # Halo Density (Solar mass/kpc^3)
# Bulge
fit_params.add('bpref', value=1,     min=0.5, max=100)  # Bulge Prefactor
# Disk
fit_params.add('dpref', value=1,     min=0.5, max=100)  # Disk Prefactor

# Do fit
fit = fit_mod.fit(Vobs,fit_params,r=Rad,weights=weighdata)

#################################
### Define fitting parameters ###
#################################

bestfit = fit.best_fit

fit_dict = fit.best_values
best_bpref = fit_dict['bpref']
best_dpref = fit_dict['dpref']
best_rc = fit_dict['rc']
best_rho0 = fit_dict['rho0']

#####################
###### Widget #######
#####################

# Define plotting function
def widgetfunction(bpref,dpref,rc,rho0):

    # Define radius
    r = np.linspace(np.min(Rad),np.max(Rad),1000)

    # Plot
    plt.figure(figsize=(11,7))
    plt.xlim(0,np.max(Rad)+0.2)
    plt.ylim(0,np.max(Vobs)+100)

    plt.errorbar(Rad,Vobs,yerr=errV,fmt='bo',label='Data')
    plt.plot(r,bulge(r,bpref),label=("Bulge"),color='orange')
    plt.plot(r,disk(r,dpref),label=("Disk"),color='purple')
    plt.plot(r,halo(r,rc,rho0),label=("Dark Matter Halo"),color='green')
    plt.plot(r,gas(r),label=("Gas"),color='blue')
    plt.plot(r,totalcurve(r,bpref,dpref,rc,rho0),label=("Total Curve"),color='red')
    plt.suptitle("Interactive Rotation Curve - Galaxy: {}".format(galaxy),fontsize='16')
    plt.title("Distance: {} Mpc".format(distance),fontsize='13')
    plt.xlabel("Radius (kpc)",fontsize='14')
    plt.ylabel("Velocity (km/s)",fontsize='14')

    # Chi squared and reduced chi squared
    # Residuals
    residuals = Vobs - totalcurve(Rad,bpref,dpref,rc,rho0)
    # Chi squared
    chisquared = np.sum(residuals**2/errV**2)
    dof = len(Rad) - 4                 # number of degrees of freedom = number of observed data - number of fitting parameters
    reducedchisquared = chisquared / dof

    if len(Rad) <= 4:
        warning_chi = "Number of measured points are too low to calculate the reduced chi-squared value."
    else:
        warning_chi = ""

    plt.annotate(r"$\chi^2$: {:.5f}".format(chisquared)+'\n'
                 +r"Reduced $\chi^2$: {:.5f}".format(reducedchisquared)+'\n'
                 +r"Note: {}".format(warning_chi)
                 + "{}".format(warning_bulge)
                 + "{}".format(warning_disk)
                 + "{}".format(warning_gas),
            xy=(0, 0), xytext=(0,0),
            xycoords=('axes fraction', 'figure fraction'),
            textcoords='offset points',
            size=12, ha='left', va='bottom')
    plt.tight_layout()
    plt.legend(bbox_to_anchor=(1,1), loc="upper left")
    plt.show()

###############################
######### Appearance ##########
###############################

style = {'description_width': 'initial'}
layout = {'width':'600px'}

################################
######## Define Sliders ########
################################

if nobulge == False:
    bpref = FloatSlider(min=0, max=5, step=0.1,
                        value=best_bpref,
                        description='Bulge Prefactor',
                        readout_format='.2f',
                        orientation='horizontal',
                        style=style, layout=layout)
if nobulge == True:
    bpref = FloatSlider(min=0, max=5, step=0.1,
                        value=best_bpref,
                        description='Bulge Prefactor',
                        readout_format='.2f',
                        orientation='horizontal',
                        style=style, layout=layout,
                        disabled=True)

if nodisk == False:
    dpref = FloatSlider(min=0, max=5, step=0.1,
                        value=best_dpref,
                        description='Disk Prefactor',
                        readout_format='.2f',
                        orientation='horizontal',
                        style=style, layout=layout)

if nodisk == True:
    dpref = FloatSlider(min=0, max=5, step=0.1,
                        value=best_dpref,
                        description='Disk Prefactor',
                        readout_format='.2f',
                        orientation='horizontal',
                        style=style, layout=layout,
                        disabled=True)

rc = FloatSlider(min=0, max=20, step=0.1,
                 value=best_rc,
                 description='Halo Core Radius [$kpc$]',
                 readout_format='.2f',
                 orientation='horizontal',
                 style=style, layout=layout)
#rc = fixed(best_rc)

rho0 = FloatSlider(min=0, max=1e11, step=1e5,
                    value=best_rho0,
                    description='Halo Central Mass Density [$M_{\odot} / kpc^3$]',
                    readout_format='.2e',
                    orientation='horizontal',
                    style=style, layout=layout)

################################
######### Interactive ##########
################################

def interactive_plot(widgetfunction):
    interact = interactive(widgetfunction,  bpref = bpref,
                               dpref = dpref,
                               rc = rc,
                               rho0 = rho0,
                               continuous_update=False)
    return interact

################################
########### Button #############
################################

button = Button(
    description="Best Fit",
    button_style='warning', # 'success', 'info', 'warning', 'danger' or ''
    icon='check')
out = Output()

def on_button_clicked(_):
    bpref.value = best_bpref
    dpref.value = best_dpref
    rc.value = best_rc
    rho0.value = best_rho0

button.on_click(on_button_clicked)

################################
######## Galaxy Image ##########
################################

from astroquery.skyview import SkyView
from astropy.wcs import WCS

def GalaxyImage(position=galaxy,survey=['DSS']):

    # DSS images of the target
    hdu = SkyView.get_images(galaxy, survey=survey)[0][0]
    gfilter = hdu.data

    # WCS
    plt.figure(figsize=(5,5))
    wcs = WCS(hdu.header)
    ax = plt.gca(projection=wcs)

    # Plot galaxy image
    ax.imshow(hdu.data, vmin=np.percentile(gfilter,0), vmax=np.percentile(gfilter,100), cmap='plasma')
    ax.set(xlabel="RA", ylabel="Dec")
    plt.title("{}".format(galaxy),fontsize='14')
    plt.show()

reason: import numpy,from scipy,from astropy


text:
#
# @Author: kuroitu (2020)
# @email: skuroitu@gmail.com
#

import numpy as np

from ..dual import *
from .const import *

##########
# 数学に関連するユーティリティ math related utilities
#######
def deg2rad(obj):
    obj = to_dual(obj)
    return obj * pi / 180

def rad2deg(obj):
    obj = to_dual(obj)
    return obj * 180 / pi

def floor(obj):
    obj = to_dual(obj)
    return Dual(np.floor(obj.re), np.floor(obj.im))

def trunc(obj):
    obj = to_dual(obj)
    return Dual(np.trunc(obj.re), np.trunc(obj.im))

def ceil(obj):
    obj = to_dual(obj)
    return Dual(np.ceil(obj.re), np.ceil(obj.im))

def round(obj, decimals=0):
    obj = to_dual(obj)
    return Dual(np.round(obj.re, decimals=decimals),
                np.round(obj.im, decimals=decimals))

def rint(obj):
    obj = to_dual(obj)
    return Dual(np.rint(obj.re), np.rint(obj.im), dtype=np.int)

def fix(obj):
    obj = to_dual(obj)
    return Dual(np.fix(obj.re), np.fix(obj.im), dtype=np.int)

def real(obj):
    obj = to_dual(obj)
    return obj.re

def imag(obj):
    obj = to_dual(obj)
    return obj.im

def conj(obj):
    obj = to_dual(obj)
    return Dual(obj.re, - obj.im)

def absolute(obj):
    return abs(obj)

def _compare(obj1, obj2, func):
    obj1, obj2 = to_dual(obj1), to_dual(obj2)

    result = Dual(np.zeros(obj1.shape), np.zeros(obj1.shape))
    result.re = func(obj1.re, obj2.re)
    indices = obj1.re == result.re
    result.im[indices] = obj1.im[indices]
    result.im[~indices] = obj2.im[~indices]

    eq_indices = obj1.re == obj2.re
    im_elements = func(obj1.im, obj2.im)
    im_indices = im_elements == obj1.im
    result.re[im_indices * eq_indices] = obj1.re[im_indices * eq_indices]
    result.im[im_indices * eq_indices] = obj1.im[im_indices * eq_indices]
    result.re[~im_indices * eq_indices] = obj2.re[~im_indices * eq_indices]
    result.im[~im_indices * eq_indices] = obj2.im[~im_indices * eq_indices]
    return result

def fmax(obj1, obj2):
    return _compare(obj1, obj2, np.fmax)

def fmin(obj1, obj2):
    return _compare(obj1, obj2, np.fmin)

def maximum(obj1, obj2):
    return _compare(obj1, obj2, np.maximum)

def minimum(obj1, obj2):
    return _compare(obj1, obj2, np.minimum)

def isnan(obj):
    obj = to_dual(obj)
    return np.isnan(obj.re) | np.isnan(obj.im)

def isfinite(obj):
    obj = to_dual(obj)
    return np.isfinite(obj.re) | np.isfinite(obj.im)

def isinf(obj):
    obj = to_dual(obj)
    return np.isinf(obj.re) | np.isinf(obj.im)

def sign(obj):
    obj = to_dual(obj)
    return np.sign(obj.re)

def sign_imag(obj):
    obj = to_dual(obj)
    return np.sign(obj.im)

def sign_dual(obj):
    obj = to_dual(obj)
    return Dual(np.sign(obj.re), np.sign(obj.im))

def signbit(obj):
    obj = to_dual(obj)
    return np.signbit(obj.re)

def signbit_imag(obj):
    obj = to_dual(obj)
    return np.signbit(obj.im)

def signbit_dual(obj):
    obj = to_dual(obj)
    return {"re": np.signbit(obj.re), \
            "im": np.signbit(obj.im)}

def copysign(obj1, obj2):
    obj1, obj2 = to_dual(obj1), to_dual(obj2)
    return Dual(np.copysign(obj1.re, obj2.re),
                np.copysign(obj1.im, obj2.im))

max_stars_repo_name: YumaShimomoto/dual

max_issues_repo_name: YumaShimomoto/dual

max_forks_repo_name: YumaShimomoto/dual

reason: import numpy


text:
import random
from collections import defaultdict
from dataclasses import dataclass
from pathlib import Path
from typing import Dict, List, Tuple, Union

import numpy as np
import tgt
import torch
from torch.utils.data import Dataset
from tqdm import tqdm

from src.data_process.config import VCTKDatasetParams

NUMBER = Union[int, float]
PHONES_TIER = "phones"
PAD_TOKEN = "<PAD>"

@dataclass
class VCTKSample:

    phonemes: List[int]
    num_phonemes: int
    speaker_id: int
    durations: np.ndarray
    mels: torch.Tensor

@dataclass
class VCTKInfo:

    text_path: Path
    mel_path: Path
    speaker_id: int
    phonemes_length: int

@dataclass
class VCTKBatch:

    phonemes: torch.Tensor
    num_phonemes: torch.Tensor
    speaker_ids: torch.Tensor
    durations: torch.Tensor
    mels: torch.Tensor

class VCTKDataset(Dataset[VCTKSample]):

    def __init__(
            self,
            sample_rate: int,
            hop_size: int,
            mels_mean: torch.Tensor,
            mels_std: torch.Tensor,
            phoneme_to_ids: Dict[str, int],
            data: List[VCTKInfo]
    ):
        self._phoneme_to_id = phoneme_to_ids
        self._dataset = data
        self._dataset.sort(key=lambda x: x.phonemes_length)
        self.sample_rate = sample_rate
        self.hop_size = hop_size
        self.mels_mean = mels_mean
        self.mels_std = mels_std

    def __len__(self) -> int:
        return len(self._dataset)

    def __getitem__(self, idx: int) -> VCTKSample:

        info = self._dataset[idx]
        text_grid = tgt.read_textgrid(info.text_path)
        phones_tier = text_grid.get_tier_by_name(PHONES_TIER)
        phoneme_ids = [self._phoneme_to_id[x.text] for x in phones_tier.get_copy_with_gaps_filled()]

        durations = np.array(
            [
                self.seconds_to_frame(x.duration())
                for x in phones_tier.get_copy_with_gaps_filled()
            ],
            dtype=np.float32
        )

        mels: torch.Tensor = torch.load(info.mel_path)
        mels = (mels - self.mels_mean) / self.mels_std

        pad_size = mels.shape[-1] - np.int64(durations.sum())
        if pad_size < 0:
            durations[-1] += pad_size
            assert durations[-1] >= 0
        if pad_size > 0:
            phoneme_ids.append(self._phoneme_to_id[PAD_TOKEN])
            np.append(durations, pad_size)

        return VCTKSample(
            phonemes=phoneme_ids,
            num_phonemes=len(phoneme_ids),
            speaker_id=info.speaker_id,
            mels=mels,
            durations=durations
        )

    def seconds_to_frame(self, seconds: float) -> float:
        return seconds * self.sample_rate / self.hop_size

class VCTKFactory:

    """Create VCTK Dataset

    Note:
        * All the speeches from speaker ``p315`` will be skipped due to the lack of the corresponding text files.
        * All the speeches from ``p280`` will be skipped for ``mic_id="mic2"`` due to the lack of the audio files.
        * Some of the speeches from speaker ``p362`` will be skipped due to the lack of  the audio files.
        * See Also: https://datashare.is.ed.ac.uk/handle/10283/3443
        * Make sure to put the files as the following structure:
            text
            ├── p225
            |   ├──p225_001.TextGrid
            |   ├──p225_002.TextGrid
            |   └──...
            └── pXXX
                ├──pXXX_YYY.TextGrid
                └──...
            mels
            ├── p225
            |   ├──p225_001.pkl
            |   ├──p225_002.pkl
            |   └──...
            └── pXXX
                ├──pXXX_YYY.pkl
                └──...
    """

    def __init__(
        self,
        sample_rate: int,
        hop_size: int,
        n_mels: int,
        config: VCTKDatasetParams,
        phonemes_to_id: Dict[str, int],
        speakers_to_id: Dict[str, int],
    ):
        self.sample_rate = sample_rate
        self.hop_size = hop_size
        self.n_mels = n_mels
        self._mels_dir = Path(config.mels_dir)
        self._text_dir = Path(config.text_dir)
        self._text_ext = config.text_ext
        self._mels_ext = config.mels_ext
        self.phoneme_to_id: Dict[str, int] = phonemes_to_id
        self.phoneme_to_id[PAD_TOKEN] = 0
        self.speaker_to_id: Dict[str, int] = speakers_to_id
        self._dataset: List[VCTKInfo] = self._build_dataset()
        self.mels_mean, self.mels_std = self._get_mean_and_std()

    @staticmethod
    def add_to_mapping(mapping: Dict[str, int], token: str) -> None:
        if token not in mapping:
            mapping[token] = len(mapping)

    def split_train_valid(
        self, test_fraction: float
    ) -> Tuple[VCTKDataset, VCTKDataset]:
        speakers_to_data_id: Dict[int, List[int]] = defaultdict(list)

        for i, sample in enumerate(self._dataset):
            speakers_to_data_id[sample.speaker_id].append(i)
        test_ids: List[int] = []
        for ids in speakers_to_data_id.values():
            test_size = int(len(ids) * test_fraction)
            if test_size > 0:
                test_indexes = random.choices(ids, k=test_size)
                test_ids.extend(test_indexes)

        train_data = []
        test_data = []
        for i in range(len(self._dataset)):
            if i in test_ids:
                test_data.append(self._dataset[i])
            else:
                train_data.append(self._dataset[i])
        train_dataset = VCTKDataset(
            sample_rate=self.sample_rate,
            hop_size=self.hop_size,
            mels_mean=self.mels_mean,
            mels_std=self.mels_std,
            phoneme_to_ids=self.phoneme_to_id,
            data=train_data
        )
        test_dataset = VCTKDataset(
            sample_rate=self.sample_rate,
            hop_size=self.hop_size,
            mels_mean=self.mels_mean,
            mels_std=self.mels_std,
            phoneme_to_ids=self.phoneme_to_id,
            data=test_data
        )
        return train_dataset, test_dataset

    def _build_dataset(self) -> List[VCTKInfo]:

        dataset: List[VCTKInfo] = []
        texts_set = {
            Path(x.parent.name) / x.stem
            for x in self._text_dir.rglob(f"*{self._text_ext}")
        }
        mels_set = {
            Path(x.parent.name) / x.stem
            for x in self._mels_dir.rglob(f"*{self._mels_ext}")
        }
        samples = list(mels_set & texts_set)
        for sample in tqdm(samples):
            tg_path = (self._text_dir / sample).with_suffix(self._text_ext)
            mels_path = (self._mels_dir / sample).with_suffix(self._mels_ext)

            text_grid = tgt.read_textgrid(tg_path)
            self.add_to_mapping(self.speaker_to_id, sample.parent.name)
            speaker_id = self.speaker_to_id[sample.parent.name]

            if PHONES_TIER in text_grid.get_tier_names():

                phones_tier = text_grid.get_tier_by_name(PHONES_TIER)
                phonemes = [x.text for x in phones_tier.get_copy_with_gaps_filled()]

                for phoneme in phonemes:
                    self.add_to_mapping(self.phoneme_to_id, phoneme)

                dataset.append(
                    VCTKInfo(
                        text_path=tg_path,
                        mel_path=mels_path,
                        phonemes_length=len(phonemes),
                        speaker_id=speaker_id
                    )
                )

        return dataset

    def _get_mean_and_std(self) -> Tuple[torch.Tensor, torch.Tensor]:
        mel_sum = torch.zeros((self.n_mels), dtype=torch.float64)
        mel_squared_sum = torch.zeros((self.n_mels), dtype=torch.float64)

        print(torch.load(self._dataset[0].mel_path).shape)

        for info in tqdm(self._dataset, desc="Computing mels mean and std"):
            mels: torch.Tensor = torch.load(info.mel_path)
            mels_mean: torch.Tensor = mels.mean(dim=-1).squeeze(0)
            mel_sum += mels_mean
            mel_squared_sum += mels_mean.pow(2)

        mels_mean: torch.Tensor = mel_sum / len(self._dataset)
        mels_std: torch.Tensor = (
            mel_squared_sum - mel_sum * mel_sum / len(self._dataset)
        ) / len(self._dataset)

        return mels_mean.view(-1, 1), mels_std.view(-1, 1)

class VCTKCollate:
    """
    Zero-pads model inputs and targets based on number of frames per setep
    """

    def __init__(self, n_frames_per_step: int = 1):
        self.n_frames_per_step = n_frames_per_step

    def __call__(self, batch: List[VCTKSample]) -> VCTKBatch:
        """Collate's training batch from normalized text and mel-spectrogram
        PARAMS
        ------
        batch: [{...}, {...}, ...]
        """
        # Right zero-pad all one-hot text sequences to max input length
        batch_size = len(batch)
        input_lengths, ids_sorted_decreasing = torch.sort(
            torch.LongTensor([len(x.phonemes) for x in batch]),
            dim=0,
            descending=True,
        )
        max_input_len = int(input_lengths[0])

        input_speaker_ids = torch.LongTensor(
            [batch[i].speaker_id for i in ids_sorted_decreasing]
        )

        text_padded = torch.zeros((batch_size, max_input_len), dtype=torch.long)
        durations_padded = torch.zeros((batch_size, max_input_len), dtype=torch.float)
        for i, idx in enumerate(ids_sorted_decreasing):
            text = batch[idx].phonemes
            text_padded[i, : len(text)] = torch.LongTensor(text)
            durations = batch[idx].durations
            durations_padded[i, : len(durations)] = torch.FloatTensor(durations)

        num_mels = batch[0].mels.squeeze(0).size(0)
        max_target_len = max([x.mels.squeeze(0).size(1) for x in batch])

        # include mel padded and gate padded
        mel_padded = torch.zeros(
            (batch_size, num_mels, max_target_len), dtype=torch.float
        )
        for i, idx in enumerate(ids_sorted_decreasing):
            mel: torch.Tensor = batch[idx].mels.squeeze(0)
            mel_padded[i, :, : mel.shape[1]] = mel
        mel_padded = mel_padded.permute(0, 2, 1)

        return VCTKBatch(
            phonemes=text_padded,
            num_phonemes=input_lengths,
            speaker_ids=input_speaker_ids,
            durations=durations_padded,
            mels=mel_padded
        )

max_stars_repo_name: IMDxD/NonAttentiveTacotron

max_issues_repo_name: IMDxD/NonAttentiveTacotron

max_forks_repo_name: IMDxD/NonAttentiveTacotron


text:
import collections

import numpy as np
import pytest

from project_code.classes import BranchTypeEnum
from project_code.topology_functions import create_coupler_mapping, \
    apply_couplers_on_branches_and_generators, remove_branches_with_loop_elements, \
    merge_tie_lines, get_most_connected_node, assign_nodes_to_ring_0, assign_nodes_to_other_rings, \
    remove_non_connected_nodes_and_branches, validate_topology, convert_couplers_to_lines

def test_branches_generators(branches_generators_nodes):
    if branches_generators_nodes.branches is None:
        print('Branches None, most likely because file not found')
        return
    print_examples(branches_generators_nodes.name,
                   branches_generators_nodes.branches, 'Branch')
    print_stats([line.impedance for line in branches_generators_nodes.branches],
                "Impedance")
    print_stats([line.PATL for line in branches_generators_nodes.branches], "PATL")
    print_examples(branches_generators_nodes.name,
                   branches_generators_nodes.generators, 'Generator')
    print_stats([gen.power for gen in branches_generators_nodes.generators], "Power")

def test_apply_couplers_on_branches(branches_generators_nodes):
    if branches_generators_nodes.branches is None:
        print('Branches None, most likely because file not found')
        return

    branches_generators_nodes.settings.do_merge_couplers = True
    branches = branches_generators_nodes.branches
    generators = branches_generators_nodes.generators
    nodes = branches_generators_nodes.nodes

    validate_topology(nodes, branches)

    couplers = [b for b in branches if b.type == BranchTypeEnum.Coupler]
    n_nodes = len(nodes)
    n_couplers = len(couplers)
    n_branches = len(branches)

    dict_couplers = create_coupler_mapping(couplers)
    apply_couplers_on_branches_and_generators(branches, generators, nodes)

    validate_topology(nodes, branches)

    n_nodes_after = len(nodes)
    n_branches_after = len(branches)
    set_replaced_branches = set(dict_couplers.keys()) - set(dict_couplers.values())
    for branch in branches:
        assert branch.name_from not in set_replaced_branches
        assert branch.name_to not in set_replaced_branches
    print(f'{n_couplers} couplers applied on a total of {n_nodes} nodes and {n_branches} branches')
    print(f'After application, left with {n_nodes_after} nodes and {n_branches_after} branches')
    print_examples(branches_generators_nodes.name, couplers, 'Bus bar coupler')
    print_stats([tr.impedance for tr in couplers], "Impedance")
    print_stats([tr.PATL for tr in couplers], "PATL")

def test_merge_tie_lines(branches_generators_nodes):
    if branches_generators_nodes.branches is None:
        print('Branches None, most likely because file not found')
        return
    branches = branches_generators_nodes.branches
    nodes = branches_generators_nodes.nodes

    validate_topology(nodes, branches)

    n_x_nodes = len([node for node in nodes if node.is_x_node()])
    # to account for x_nodes with more than 2 branches:
    n_redundant_branches = sum([max(len(node.branches)-2, 0) for node in nodes if node.is_x_node()])
    if n_x_nodes > 0:
        x_node = [node for node in nodes if node.is_x_node()][0]
        branch_names = (x_node.branches[0].name_branch, x_node.branches[1].name_branch)
        branch_node_0 = [n for n in [x_node.branches[0].node_from, x_node.branches[0].node_to]
                         if n != x_node]
        branch_node_1 = [n for n in [x_node.branches[1].node_from, x_node.branches[1].node_to]
                         if n != x_node]
        branch_imps = (x_node.branches[0].impedance, x_node.branches[1].impedance)
        branch_PATLs = (x_node.branches[0].PATL, x_node.branches[1].PATL)
    n_nodes = len(nodes)
    n_branches = len(branches)

    merge_tie_lines(branches, nodes)

    assert len([node for node in nodes if node.is_x_node()]) == 0
    assert len(nodes) == n_nodes - n_x_nodes
    # for each x-node, two branches replaced by 1, and correct for redundant branches:
    assert len(branches) == n_branches - n_x_nodes - n_redundant_branches
    if n_x_nodes > 0:
        print(f"Pre merge: {x_node.name} with branches {branch_names}, impedances {branch_imps}, "
              f"and PATLs {branch_PATLs} - between nodes {branch_node_0[0].name} and "
              f"{branch_node_1[0].name}.")
        # noinspection PyUnboundLocalVariable
        new_branch = [b for b in branches
                      if (branch_node_0[0] in (b.node_from, b.node_to))
                      and (branch_node_1[0] in (b.node_from, b.node_to))][0]
        print(f"Post merge: branch {new_branch.name_branch} with impedance "
              f"{new_branch.impedance:.5} and PATL {new_branch.PATL:.5} between "
              f"{new_branch.name_from} and {new_branch.name_to}.")

def test_remove_branches_with_loop_elements(branches_generators_nodes):
    if branches_generators_nodes.branches is None:
        print('Branches None, most likely because file not found')
        return
    branches = branches_generators_nodes.branches
    nodes = branches_generators_nodes.nodes
    nbranches_org = len(branches)
    remove_branches_with_loop_elements(branches, nodes)
    for branch in branches:
        assert branch.name_from != branch.name_to
    validate_topology(nodes, branches)
    print(f'Removed {nbranches_org - len(branches)} branches, {len(branches)} branches left.')

@pytest.mark.parametrize("do_merge_couplers", [True, False])
def test_set_branch_country_onoff_merge(do_merge_couplers, branches_generators_nodes):
    if branches_generators_nodes.branches is None:
        print('Branches None, most likely because file not found')
        return
    name, branches, generators, nodes, settings = branches_generators_nodes

    if do_merge_couplers:
        apply_couplers_on_branches_and_generators(branches, generators, nodes)
    else:
        convert_couplers_to_lines(branches)

    merge_tie_lines(branches, nodes)
    remove_branches_with_loop_elements(branches, nodes)

    assert all([b.country is not None for b in branches])
    assert all([n.country is not None for n in nodes])
    neg = ' ' if do_merge_couplers else ' not '
    nxnodes = len([n for n in nodes if n.country == 'X'])
    ntiebranches = len([b for b in branches if b.country == 'TIE'])
    print(f'{nxnodes} x-nodes, {ntiebranches} tie-lines ({nxnodes/len(nodes)*100:.1f} pct of nodes, '
          f'{ntiebranches/len(branches)*100:.1f} pct of branches)')
    print_list_occurrence(f'Nodes per country when couplers are{neg}merged',
                          [n.country for n in nodes])
    print_list_occurrence(f'Branches per country when couplers are{neg}merged',
                          [b.country for b in branches])

def test_get_most_connected_node(branches_generators_nodes):
    if branches_generators_nodes.branches is None:
        print('Branches None, most likely because file not found')
        return
    branches = branches_generators_nodes.branches
    nodes = branches_generators_nodes.nodes
    merge_tie_lines(branches, nodes)
    countries = branches_generators_nodes.settings.countries
    for country in countries:
        mc_node = get_most_connected_node(nodes, country)
        if mc_node is None:
            continue
        print(f"In {country}, the most connected node is {mc_node.name} with "
              f"{len(mc_node.branches)} branches.")

def test_assign_nodes_to_ring_0(branches_generators_nodes_preprocessed_merged,
                                branches_generators_nodes_preprocessed_nonmerged):

    for switch in range(2):
        if switch == 0:
            branches_generators_nodes = branches_generators_nodes_preprocessed_merged
        else:
            branches_generators_nodes = branches_generators_nodes_preprocessed_nonmerged

        if branches_generators_nodes.branches is None:
            print('Branches None, most likely because file not found')
            return
        branches = branches_generators_nodes.branches
        nodes = branches_generators_nodes.nodes
        countries = branches_generators_nodes.settings.countries

        validate_topology(nodes, branches)

        for country in countries:
            if country == 'XX':
                continue
            n_nodes_in_country = len([n for n in nodes if n.country == country])
            n_branches_in_country = len([b for b in branches if b.country == country])
            print(f'Connecting for country {country} with {n_nodes_in_country} nodes and '
                  f'{n_branches_in_country} branches:')

            assign_nodes_to_ring_0(nodes, branches, country)

            n_nodes_in_ring_0 = len([node for node in nodes if node.ring == 0])
            print(f"Ring 0 initialised with {n_nodes_in_ring_0} nodes.")
            cntry_nodes_non_conn = [n for n in nodes if n.country == country and n.ring != 0]
            for n in cntry_nodes_non_conn:
                print(f"Warning : {n.name} is not connected to {country}")
            conn_nodes_non_cntry = [n for n in nodes if n.ring == 0 and n.country != country]
            for n in conn_nodes_non_cntry:
                print(f"Warning : {n.name} is connected to {country} but belongs to {n.country}")
            assert len(conn_nodes_non_cntry) == 0

            reset_connectivity(branches, nodes)  # only for testing

def test_assign_nodes_to_other_rings(branches_generators_nodes_preprocessed_merged,
                                     branches_generators_nodes_preprocessed_nonmerged):

    for switch in range(2):
        if switch == 0:
            branches_generators_nodes = branches_generators_nodes_preprocessed_merged
        else:
            branches_generators_nodes = branches_generators_nodes_preprocessed_nonmerged

        if branches_generators_nodes.branches is None:
            print('Branches None, most likely because file not found')
            return
        branches = branches_generators_nodes.branches
        nodes = branches_generators_nodes.nodes
        for country in branches_generators_nodes.settings.countries:
            if country == 'XX':
                continue
            n_nodes_in_country = len([n for n in nodes if n.country == country])
            n_branches_in_country = len([b for b in branches if b.country == country])
            print(f'\nConnecting for country {country} with {n_nodes_in_country} nodes and '
                  f'{n_branches_in_country} branches:')

            assign_nodes_to_ring_0(nodes, branches, country)
            assign_nodes_to_other_rings(nodes)

            # Test for consistency
            for node in [node.name for node in nodes if node.ring == 99 and node.connected]:
                print("Node " + node + " is connected but has no ring")
            assert len([node.name for node in nodes if node.ring == 99 and node.connected]) == 0
            for node in [node.name for node in nodes if node.ring < 99 and not node.connected]:
                print("Node " + node + " is in a ring but is not connected")
            assert len([node.name for node in nodes if node.ring < 99 and not node.connected]) == 0

            # Show results of first rings
            max_ring_nr_shown = 3
            for ring_idx in range(max_ring_nr_shown+1):
                data_nodes_in_ring = [(n.name, n.country) for n in nodes if n.ring == ring_idx]
                if ring_idx == 0:
                    print(f"{len(data_nodes_in_ring)} nodes in ring {ring_idx}")
                else:
                    print(f"{len(data_nodes_in_ring)} nodes in ring {ring_idx}: {data_nodes_in_ring}")

            reset_connectivity(branches, nodes)  # only for testing

def test_remove_non_connected_nodes_and_branches(branches_generators_nodes):
    if branches_generators_nodes.branches is None:
        print('Branches None, most likely because file not found')
        return
    branches = branches_generators_nodes.branches
    nodes = branches_generators_nodes.nodes
    merge_tie_lines(branches, nodes)

    # To limit time use of test, only sample first 5 countries
    countries = branches_generators_nodes.settings.countries[:5]

    for country in countries:
        if country == 'XX':
            continue
        assign_nodes_to_ring_0(nodes, branches, country)
        assign_nodes_to_other_rings(nodes)
        n_nodes = len(nodes)
        n_branches = len(branches)

        nodes_new, branches_new = remove_non_connected_nodes_and_branches(nodes)

        assert len([n for n in nodes_new if not n.connected]) == 0, \
            f"Non-connected nodes in node set for {country}."
        assert len([b for b in branches_new if not b.node_from.connected]) == 0, \
            f"Branches with non-connected nodes in node set for {country}."
        assert len([b for b in branches_new if not b.node_to.connected]) == 0, \
            f"Branches with non-connected nodes in node set for {country}."

        print(f"{(1 - len(nodes_new)/n_nodes)*100:.1f}% of nodes removed, ",
              f"{(1 - len(branches_new)/n_branches)*100:.1f}% of branches removed")

        # test for a too high (3 pct is a bit random treshold) number of elements lost
        assert (1 - len(nodes_new)/n_nodes)*100 < 3
        assert (1 - len(branches_new)/n_branches)*100 < 3
        reset_connectivity(branches, nodes)  # only for testing

def print_examples(name, lst, asset_type):
    print(f"\n {name}: {len(lst)} {asset_type}s")
    print(lst[0])
    print(lst[-1])

def print_stats(lst, title):
    print(f"{title } stats: mean {np.mean(lst)}, std {np.std(lst)}, "
          f"median {np.median(lst)}, min {np.min(lst)}, max {np.max(lst)}")

def print_list_occurrence(name, lst):
    print(f"{name}:")
    # noinspection PyArgumentList
    d = dict(collections.Counter(lst))
    for key in sorted(d, key=d.get, reverse=True):
        print(f"{key}\t{d[key]}")

def reset_connectivity(branches, nodes):
    for branch in branches:
        branch.ring = 99
    for node in nodes:
        node.ring = 99
        node.connected = False


text:
# -*- coding: utf-8 -*-
"""
Created on Fri Oct 22 2021

@author: Qirui Zhang

Multiple classes for controlling Sky130TempSensor I/Os through FT232H
"""

from pyftdi.gpio import GpioMpsseController
import time
import numpy as np

'''
GPIO Board USB addresses
'''
gpio_in_addr = "ftdi://ftdi:232h:00:fd/1" # {SEL_INST[1:0], SEL_GRP[2:0], SEL_DESIGN, SEL_CTR[3:0], RESET_REGn} = {C7:C0, D7:D5}
gpio_out1_addr = "ftdi://ftdi:232h:00:fe/1" # {C7:C0, D7:D4} = DOUT[23:12]
gpio_out0_addr = "ftdi://ftdi:232h:00:ff/1" # {C7:C0, D7:D3} = {DOUT[11:0], DONE}

'''
Conversion Time and Repeat-sampling LUT for different Design Groups
'''
CTR_LUT = [
        6, # design0, grp0
        6, # design0, grp1
        6, # design0, grp2
        6, # design0, grp3
        7, # design0, grp4
        7, # design0, grp5
        7, # design0, grp6
        7, # design0, grp7
        0, # design1, grp0
        0, # design1, grp1
        0, # design1, grp2
        0, # design1, grp3
        0, # design1, grp4
        0, # design1, grp5
        0, # design1, grp6
        0, # design1, grp7
    ]

'''
The Temperature Sensor Chip Control Class
'''
class tempsensorIO():
    def __init__(self, input_addr_url, out1_addr_url, out0_addr_url, freq):
        self.input_url = input_addr_url
        self.out1_url = out1_addr_url
        self.out0_url = out0_addr_url
        self.freq = freq
        self.direction_in = 0xffff
        self.direction_out1 = 0x0000
        self.direction_out0 = 0x0000

        '''
        Input Board Pinout:
            C7-C6:  SEL_INST[1:0]. Default 0b00.
            C5-C3:  SEL_GRP[2:0]. Default 0b000.
            C2:     SEL_DESIGN. Default 0b0.
            C1-C0,
            D7-D6:  SEL_CTR[3:0]. Default 0b0000.
            D5:     RESET_REGn. Default 0b0.
        '''
        self.gpio_in = GpioMpsseController()
        self.gpio_in.configure(self.input_url, direction=self.direction_in, frequency=self.freq)
        self.gpio_in.write(0x0000)

        '''
        Output Board1 Pinout:
            C7-C0,
            D7-D4:  DOUT[23:12]

        Output Board0 Pinout:
            C7-C0,
            D7-D4:  DOUT[11:0]
            D3:     DONE
        '''
        self.gpio_out1 = GpioMpsseController()
        self.gpio_out1.configure(self.out1_url, direction=self.direction_out1, frequency=self.freq)

        self.gpio_out0 = GpioMpsseController()
        self.gpio_out0.configure(self.out0_url, direction=self.direction_out0, frequency=self.freq)

    def reset(self):
        self.gpio_in.write(0x0000)

    def close(self):
        self.gpio_in.close()
        self.gpio_out1.close()
        self.gpio_out0.close()

    def update_freq(self, freq):
        self.freq = freq
        self.gpio_in.set_frequency(self.freq)
        self.gpio_out1.set_frequency(self.freq)
        self.gpio_out0.set_frequency(self.freq)

    def get_state(self):
        return [f'{self.gpio_in.read()[0]:016b}', f'{self.gpio_out1.read()[0]:016b}', f'{self.gpio_out0.read()[0]:016b}']

    # active low
    def chip_reset(self, rst_val):
        cur_state = self.gpio_in.read()[0]
        cur_state &= ~(0b1 << 5)
        self.gpio_in.write(cur_state | rst_val << 5)

    def set_sel_ctr(self, sel_ctr_val):
        cur_state = self.gpio_in.read()[0]
        cur_state &= ~(0b1111 << 6)
        self.gpio_in.write(cur_state | sel_ctr_val << 6)

    def set_sel_design(self, sel_design_val):
        cur_state = self.gpio_in.read()[0]
        cur_state &= ~(0b1 << 10)
        self.gpio_in.write(cur_state | sel_design_val << 10)

    def set_sel_grp(self, sel_grp_val):
        cur_state = self.gpio_in.read()[0]
        cur_state &= ~(0b111 << 11)
        self.gpio_in.write(cur_state | sel_grp_val << 11)

    def set_sel_inst(self, sel_inst_val):
        cur_state = self.gpio_in.read()[0]
        cur_state &= ~(0b11 << 14)
        self.gpio_in.write(cur_state | sel_inst_val << 14)

    def get_done(self):
        return (self.gpio_out0.read()[0] >> 3) & 0b1

    def get_dout(self):
        return ((self.gpio_out1.read()[0] >> 4) << 12) + (self.gpio_out0.read()[0] >> 4)

    # Test RO frequency of all 64 designs on a chip
    def test_all_freqs(self, ctr_design0, ctr_design1, repeat0, repeat1, temp, freq_ref):
        # initialize columns in the measurement results table
        list_temp = []
        list_design = []
        list_grp = []
        list_inst = []
        list_ctr = []
        list_ref = []
        list_dout = []
        list_freq = []

        # Start looping all the 64 designs in a chip under a given temperature
        for sel_design in range(2):
            # select input ctr configurations according to design type
            if sel_design == 0:
                sel_ctr = ctr_design0
                repeat = repeat0
            else:
                sel_ctr = ctr_design1
                repeat = repeat1

            for sel_grp in range(8):
                for sel_inst in range(4):
                    print("Testing temperature sensor node " + str(sel_design) + '.' + str(sel_grp) + '.' + str(sel_inst))
                    # Set CTR
                    self.set_sel_ctr(sel_ctr)

                    # Select design
                    self.set_sel_design(sel_design)
                    self.set_sel_grp(sel_grp)
                    self.set_sel_inst(sel_inst)

                    freq_avg = 0
                    for i in range(repeat):
                        # reset chip
                        self.chip_reset(0)
                        time.sleep(1)

                        # release reset
                        self.chip_reset(1)

                        # Wait for done and read dout
                        while True:
                            done = self.get_done()
                            if done:
                                print("** DONE DETECTED **")
                                dout = self.get_dout()
                                print("DOUT result is " + str(dout))
                                break

                        freq = (dout/(16.0*(2**sel_ctr)))*freq_ref/2
                        freq_avg += freq

                    freq_avg /= repeat
                    print("Frequency result is " + str(freq_avg) + " kHz\n")

                    list_temp.append(temp)
                    list_design.append(sel_design)
                    list_grp.append(sel_grp)
                    list_inst.append(sel_inst)
                    list_ctr.append(sel_ctr)
                    list_ref.append(freq_ref)
                    list_dout.append(dout)
                    list_freq.append(freq_avg) # kHz

        dict_meas = {
            'temp':             list_temp,
            'design':           list_design,
            'group':            list_grp,
            'inst':             list_inst,
            'CTR':              list_ctr,
            'FreqRef (kHz)':    list_ref,
            'DOUT':             list_dout,
            'Freq (kHz)':       list_freq
        }

        return dict_meas

    # Test RO frequency of all 64 designs on a chip
    def test_all_freqs_wlut(self, repeat, temp, freq_ref):
        # initialize columns in the measurement results table
        list_temp = []
        list_design = []
        list_grp = []
        list_inst = []
        list_ctr = []
        list_ref = []
        list_dout = []
        list_freq = []

        # Start looping all the 64 designs in a chip under a given temperature
        for sel_design in range(2):
            # select input ctr configurations according to design type
            for sel_grp in range(8):
                sel_ctr = CTR_LUT[sel_design*8 + sel_grp]
                for sel_inst in range(4):
                    print("Testing temperature sensor node " + str(sel_design) + '.' + str(sel_grp) + '.' + str(sel_inst))
                    # print("Conversion Time is " + str(32.0*(2**sel_ctr)/freq_ref) + 'ms')
                    # Set CTR
                    self.set_sel_ctr(sel_ctr)

                    # Select design
                    self.set_sel_design(sel_design)
                    self.set_sel_grp(sel_grp)
                    self.set_sel_inst(sel_inst)
                    time.sleep(0.2)

                    dout_avg = 0
                    for i in range(repeat):
                        # reset chip
                        self.chip_reset(0)
                        time.sleep(0.1)

                        # release reset
                        self.chip_reset(1)

                        # Wait for done and read dout
                        while True:
                            done = self.get_done()
                            if done:
                                #print("** DONE DETECTED **")
                                dout = self.get_dout()
                                #print("DOUT result is " + str(dout))
                                break

                        dout_avg += dout

                    dout_avg /= repeat
                    freq = (dout_avg/(32.0*(2**sel_ctr)))*freq_ref

                    print("Frequency result is " + str(freq) + " kHz\n")

                    list_temp.append(temp)
                    list_design.append(sel_design)
                    list_grp.append(sel_grp)
                    list_inst.append(sel_inst)
                    list_ctr.append(sel_ctr)
                    list_ref.append(freq_ref)
                    list_dout.append(dout)
                    list_freq.append(freq) # kHz

        dict_meas = {
            'temp':             list_temp,
            'design':           list_design,
            'group':            list_grp,
            'inst':             list_inst,
            'CTR':              list_ctr,
            'FreqRef (kHz)':    list_ref,
            'DOUT':             list_dout,
            'Freq (kHz)':       list_freq
        }

        return dict_meas

    # Test RO frequency of all 64 designs on a chip
    def test_all_res_wlut(self, repeat, ctr2, temp, freq_ref):
        # initialize columns in the measurement results table
        list_ctr1 = []
        list_ctr2 = []
        list_ref = []
        list_res1 = [] # res1 uses the same conversion time for inaccuracy testing
        list_res2 = [] # res2 uses a common conversion time across all designs

        # Start looping all the 64 designs in a chip under a given temperature
        for sel_design in range(2):
            # select input ctr configurations according to design type
            for sel_grp in range(8):
                sel_ctr = CTR_LUT[sel_design*8 + sel_grp]
                for sel_inst in range(4):
                    print("Testing temperature sensor node " + str(sel_design) + '.' + str(sel_grp) + '.' + str(sel_inst))

                    # Select design
                    self.set_sel_design(sel_design)
                    self.set_sel_grp(sel_grp)
                    self.set_sel_inst(sel_inst)
                    time.sleep(0.1)

                    # Set CTR1
                    self.set_sel_ctr(sel_ctr)
                    time.sleep(0.1)
                    dout_list = []
                    for i in range(repeat):
                        # reset chip
                        self.chip_reset(0)
                        time.sleep(0.05)
                        # release reset
                        self.chip_reset(1)
                        # Wait for done and read dout
                        while True:
                            done = self.get_done()
                            if done:
                                dout = self.get_dout()
                                break
                        dout_list.append(dout)
                    dout_std = np.std(dout_list)
                    freq_std1 = (dout_std/(32.0*(2**sel_ctr)))*freq_ref
                    print("Frequency resolution-1 is " + str(freq_std1) + " kHz")

                    # Set CTR2
                    self.set_sel_ctr(ctr2)
                    time.sleep(0.1)
                    dout_list = []
                    for i in range(repeat):
                        # reset chip
                        self.chip_reset(0)
                        time.sleep(0.05)
                        # release reset
                        self.chip_reset(1)
                        # Wait for done and read dout
                        while True:
                            done = self.get_done()
                            if done:
                                dout = self.get_dout()
                                break
                        dout_list.append(dout)
                    dout_std = np.std(dout_list)
                    freq_std2 = (dout_std/(32.0*(2**ctr2)))*freq_ref
                    print("Frequency resolution-2 is " + str(freq_std2) + " kHz\n")

                    list_ctr1.append(sel_ctr)
                    list_ctr2.append(ctr2)
                    list_ref.append(freq_ref)
                    list_res1.append(freq_std1) # kHz
                    list_res2.append(freq_std2) # kHz

        dict_meas = {
            'CTR1':             list_ctr1,
            'CTR2':             list_ctr2,
            'FreqRef (kHz)':    list_ref,
            'Res1 (kHz)':       list_res1,
            'Res2 (kHz)':       list_res2
        }

        return dict_meas

    # Test chip current of all 64 designs on a chip
    def test_all_powers(self, ctr, meas_step, pmeas, temp, freq_ref, SMU):
        # initialize columns in the measurement results table
        list_Ivdd = []
        list_Ivdd1v8 = []

        # Set CTR
        self.set_sel_ctr(ctr)
        # Start looping all the 64 designs in a chip under a given temperature
        for sel_design in range(2):
            # Adaptive current measurement limits
            if sel_design == 0:
                I_VDD1v8_limit = 0.00001 # 10uA
            else:
                I_VDD1v8_limit = 0.0001 # 100uA
            SMU.write(':SENSe2:CURRent:DC:PROTection:LEVel %G' % (I_VDD1v8_limit))
            SMU.write(':INITiate:IMMediate:ALL (%s)' % ('@1,2'))
            time.sleep(10)

            # Dummy Test to heat-up SMU
            print("Dummy Test!")
                # Select design
            self.set_sel_design(sel_design)
            self.set_sel_grp(0)
            self.set_sel_inst(0)
                # reset chip
            self.chip_reset(0)
            time.sleep(1)
                # Initialize averaged Voltage, Current and Power
            Ivdd_rec = []
            Ivdd_avg = 0
            Ivdd1v8_rec = []
            Ivdd1v8_avg = 0
                # release reset
            self.chip_reset(1)
            time.sleep(1)
                # Wait for done and measure power
            while True:
                I_values = SMU.query_ascii_values(':MEASure:CURRent:DC? (%s)' % ('@1,2'))
                done = self.get_done()
                if done:
                    print("** DONE DETECTED **")
                    break
                else:
                    Ivdd_rec.append(I_values[0])
                    Ivdd1v8_rec.append(I_values[1])
                    time.sleep(meas_step)
            # Calculate average powers
            win_meas_len = len(Ivdd_rec)
            win_stable_len = int(len(Ivdd_rec)*pmeas)
            Ivdd_avg = sum(Ivdd_rec[(win_meas_len-win_stable_len):win_meas_len])/win_stable_len
            Ivdd1v8_avg = sum(Ivdd1v8_rec[(win_meas_len-win_stable_len):win_meas_len])/win_stable_len
            print("VDD1v8 dummy current is " + str(Ivdd1v8_avg*1e6) + " uA\n")

            # Formal Test
            for sel_grp in range(8):
                for sel_inst in range(4):
                    print("Testing temperature sensor node " + str(sel_design) + '.' + str(sel_grp) + '.' + str(sel_inst))
                    # Select design
                    self.set_sel_design(sel_design)
                    self.set_sel_grp(sel_grp)
                    self.set_sel_inst(sel_inst)

                    # reset chip
                    self.chip_reset(0)
                    time.sleep(1)

                    # Initialize averaged Voltage, Current and Power
                    Ivdd_rec = []
                    Ivdd_avg = 0
                    Ivdd1v8_rec = []
                    Ivdd1v8_avg = 0

                    # release reset
                    self.chip_reset(1)
                    time.sleep(1)

                    # Wait for done and measure power
                    while True:
                        I_values = SMU.query_ascii_values(':MEASure:CURRent:DC? (%s)' % ('@1,2'))
                        done = self.get_done()
                        if done:
                            print("** DONE DETECTED **")
                            break
                        else:
                            Ivdd_rec.append(I_values[0])
                            Ivdd1v8_rec.append(I_values[1])
                            time.sleep(meas_step)

                    # Calculate average powers
                    win_meas_len = len(Ivdd_rec)
                    win_stable_len = int(len(Ivdd_rec)*pmeas)
                    Ivdd_avg = sum(Ivdd_rec[(win_meas_len-win_stable_len):win_meas_len])/win_stable_len
                    Ivdd1v8_avg = sum(Ivdd1v8_rec[(win_meas_len-win_stable_len):win_meas_len])/win_stable_len
                    print("VDD1v8 current is " + str(Ivdd1v8_avg*1e6) + " uA\n")

                    list_Ivdd.append(Ivdd_avg)
                    list_Ivdd1v8.append(Ivdd1v8_avg)

        dict_meas = {
            'Ivdd (A)':     list_Ivdd,
            'Ivdd1v8 (A)':  list_Ivdd1v8
        }

        return dict_meas

'''
Some useful functions
'''

def convert_temp_read(temp):
    temp_bin = f'{temp:016b}'
    if temp_bin[0] == '1': # T < 0C
        temp = - (2**16 - temp)/10.0
    else:
        temp = temp/10.0
    return temp

def convert_temp_write(temp):
    if temp < 0: # T < 0C
        temp_set = 2**16 + temp * 10
    else:
        temp_set = temp * 10
    return temp_set

def twos_comp(val, bits):
    """compute the 2's complement of int value val"""
    if (val & (1 << (bits - 1))) != 0: # if sign bit is set e.g., 8bit: 128-255
        val = val - (1 << bits)        # compute negative value
    return val                         # return positive value as is

max_stars_repo_name: chenliyumich/OpenFASOC

max_issues_repo_name: chenliyumich/OpenFASOC

max_forks_repo_name: chenliyumich/OpenFASOC


text:
import numpy as np
import pandas as pd

def read_data(filename):
	'''
	Reads the input training data from filename and
	Returns the matrices X : [N X D] and Y : [N X 1] where D is number of features and N is the number of data points
	'''
	dataframe = pd.read_csv(filename, keep_default_na=False, na_values='')
	data = [np.array(dataframe[col]) for col in dataframe]
	for i, d in enumerate(data):
		data[i].shape = (data[i].shape[0], 1)
	data = np.concatenate(data, axis = 1)
	X = data[:,:-1]
	Y = data[:,-1]
	Y.shape = (Y.shape[0],1)
	return X, Y

def one_hot_encode(X, labels):
	'''
	X = input [N X 1] matrix data
	labels = list of all possible labels for current category
	Returns the matrix X : [N X len(labels)] in one hot encoded format
	'''
	X.shape = (X.shape[0], 1)
	newX = np.zeros((X.shape[0], len(labels)))
	label_encoding = {}
	for i, l in enumerate(labels):
		label_encoding[l] = i
	for i in range(X.shape[0]):
		newX[i, label_encoding[X[i,0]]] = 1
	return newX

def sse(X, Y, W):
	'''
	X = input feature matrix [N X D]
	Y = output values [N X 1]
	W = weight vector [D X 1]
	Returns the Sum Square Error for given set of weights W
	'''
	return np.linalg.norm(Y - X @ W) ** 2

def ridge_objective(X, Y, W, _lambda):
	'''
	X = input feature matrix [N X D]
	Y = output values [N X 1]
	W = weight vector [D X 1]
	_lambda = scalar parameter
	Returns the Ridge Loss function for given set of weights W
	'''
	return sse(X,Y,W) + _lambda * np.linalg.norm(W, ord=2)

def lasso_objective(X, Y, W, _lambda):
	'''
	X = input feature matrix [N X D]
	Y = output values [N X 1]
	W = weight vector [D X 1]
	_lambda = scalar parameter
	Returns the Lasso Loss function for given set of weights W
	'''
	return sse(X,Y,W) + _lambda * np.linalg.norm(W, ord=1)

def separate_data(X, Y):
	'''
	X = input feature matrix [N X D]
	Y = output values [N X 1]
	Segregate some part as train and some part as test
	Return the trainX, trainY, testX, testY
	'''
	trainX = X[0:1200, :]
	trainY = Y[0:1200, :]

	testX = X[1200:, :]
	testY = Y[1200:, :]

	return trainX, trainY, testX, testY

def plot_kfold(lambdas, scores):
	'''
	lambdas = list of scalar parameter lambda
	scores  = list of average SSE values for each lambda in lambdas
	Plots the kfold cross validation graph required b/w lambdas and scores
	'''
	import matplotlib.pyplot as plt
	plt.plot(lambdas, scores)
	plt.ylabel('Validation_SSE')
	plt.xlabel('Lambda')
	plt.show()

max_stars_repo_path: Lab 5/utils.py

max_issues_repo_path: Lab 5/utils.py

max_forks_repo_path: Lab 5/utils.py


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

x = np.random.randint(25,100,25)
y = np.random.randint(175,255,25)
z = np.hstack((x,y))
z = z.reshape((50,1))
z = np.float32(z)
plt.hist(z,256,[0,256]),plt.show()

# Define criteria = ( type, max_iter = 10 , epsilon = 1.0 )
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 10, 1.0)

# Set flags (Just to avoid line break in the code)
flags = cv2.KMEANS_RANDOM_CENTERS

# Apply KMeans
compactness,labels,centers = cv2.kmeans(z,2,None,criteria,10,flags)

A = z[labels==0]
B = z[labels==1]

# Now plot 'A' in red, 'B' in blue, 'centers' in yellow
plt.hist(A,256,[0,256],color = 'r')
plt.hist(B,256,[0,256],color = 'b')
plt.hist(centers,32,[0,256],color = 'y')
plt.show()

max_stars_repo_path: OpenCV Python/8. Machine Learning/3. K-Means Clusterig/code.py


text:
import numpy as np
import pickle
import json
import torch
import pickle5
from external.ldif.representation.structured_implicit_function import StructuredImplicit

def read_pkl(pkl_file, protocol=4):
    with open(pkl_file, 'rb') as file:
        if protocol == 4:
            data = pickle.load(file)
        elif protocol == 5:
            data = pickle5.load(file)
        else:
            raise NotImplementedError()
    return data

def write_pkl(obj, pkl_file):
    with open(pkl_file, 'wb') as file:
        pickle.dump(obj, file)

def read_json(json_file):
    with open(json_file, 'r') as file:
        json_data = json.load(file)
    return json_data

def write_json(obj, json_file):
    with open(json_file, 'w') as file:
        json.dump(obj, file, indent=4)

def get_flag(flags, k):
    if isinstance(flags, dict):
        flag = flags.get(k, None)
    elif isinstance(flags, (list, tuple, set)):
        flag = k in flags
    else:
        flag = flags
    return flag

def list_of_dict_to_dict_of_array(old_lis, force_list=None, force_cat=None, stop_recursion=None, to_tensor=None):
    """
    From a list of dict to a dict of array (or list)
    """
    if isinstance(old_lis, list):
        if not old_lis:
            return {}
        elif not isinstance(old_lis[0], dict):
            return old_lis
    else:
        return old_lis

    new_dic = {}
    if len(old_lis) == 0:
        return new_dic

    key_type = {}
    for d in old_lis:
        key_type.update({k: type(v) for k, v in d.items()})

    for k, t in key_type.items():
        if_force_list = get_flag(force_list, k)
        if_force_cat = get_flag(force_cat, k)
        if_stop_recursion = get_flag(stop_recursion, k)
        if_to_tensor = get_flag(to_tensor, k)

        new_value = [dic[k] for dic in old_lis if k in dic]

        if t == list and if_force_cat == True:
            new_value = [v for l in new_value for v in l]
            if_force_cat = False

        if if_stop_recursion != True:
            new_value = list_of_dict_to_dict_of_array(
                new_value, if_force_list, if_force_cat, if_stop_recursion, if_to_tensor)

        if not (isinstance(new_value, dict) or if_force_list == True):
            if t == torch.Tensor:
                new_value = torch.cat(new_value) if if_force_cat else torch.stack(new_value)
            elif t == list and if_force_cat:
                new_value = [v for l in new_value for v in l]
            elif t == StructuredImplicit:
                new_value = StructuredImplicit.cat(new_value)
            elif t not in (str, dict):
                new_value = np.concatenate(new_value) if if_force_cat else np.stack(new_value)

        if if_to_tensor == True:
            new_value = recursively_to(new_value, dtype='tensor')
        new_dic[k] = new_value

    return new_dic

def dict_of_array_to_list_of_dict(old_dic, split=None, keep_dim=False):
    """
    From a dict of array (or list) to a list of dict
    """
    new_list = []
    if len(old_dic) > 0:
        for key, value in old_dic.items():
            if isinstance(value, dict):
                value = dict_of_array_to_list_of_dict(value)
                if isinstance(split, dict) and split and key in split:
                    value = [value[start:end] for start, end in split[key]]
            if new_list:
                assert len(new_list) == len(value)
            for i_obj in range(len(value)):
                obj = value[i_obj:i_obj+1] if keep_dim else value[i_obj]
                if i_obj >= len(new_list):
                    new_list.append({})
                obj = recursively_to(obj, dtype='numpy')
                if isinstance(obj, np.ndarray) and obj.size == 1:
                    if obj.dtype in (np.float32, np.float64):
                        obj = float(obj)
                    elif obj.dtype == np.bool:
                        obj = bool(obj)
                    else:
                        obj = int(obj)
                new_list[i_obj][key] = obj
        if isinstance(split, torch.Tensor):
            new_list = [new_list[start:end] for start, end in split]
    return new_list

def recursively_to(data, dtype=None, device=None):
    if isinstance(data, dict):
        new_data = {}
        for k, v in data.items():
            new_data[k] = recursively_to(v, dtype, device)
    elif isinstance(data, (list, tuple, set)):
        new_data = []
        for v in data:
            new_data.append(recursively_to(v, dtype, device))
        new_data = type(data)(new_data)
    elif isinstance(data, torch.Tensor):
        if dtype in ('numpy', 'list'):
            new_data = data.detach().cpu().numpy()
            if dtype == 'list':
                new_data = new_data.tolist()
        elif dtype == 'cuda':
            new_data = data.cuda()
        else:
            new_data = data
        if device is not None:
            new_data = new_data.to(device)
    elif isinstance(data, np.ndarray):
        if dtype == 'list':
            new_data = data.tolist()
        elif dtype in ('tensor', 'cuda'):
            new_data = torch.from_numpy(data)
            if new_data.type() == 'torch.DoubleTensor':
                new_data = new_data.float()
            if dtype == 'cuda':
                new_data = new_data.cuda()
            if device is not None:
                new_data = new_data.to(device)
        else:
            new_data = data
    else:
        new_data = data
    return new_data

def recursively_update(dict1, dict2):
    ''' Update two config dictionaries recursively.

    Args:
        dict1 (dict): first dictionary to be updated
        dict2 (dict): second dictionary which entries should be used

    '''
    for k, v in dict2.items():
        if k not in dict1:
            dict1[k] = dict()
        if isinstance(v, dict):
            recursively_update(dict1[k], v)
        else:
            dict1[k] = v

def get_any_array(dic):
    if isinstance(dic, dict):
        vs = dic.values()
    elif isinstance(dic, (list, tuple, set)):
        vs = dic
    else:
        return dic

    for v in vs:
        t = get_any_array(v)
        if isinstance(t, (np.ndarray, torch.Tensor)):
            return t

    return dic

def recursively_ignore(dic, ignore_keys):
    if not isinstance(dic, dict):
        return dic
    if isinstance(ignore_keys, dict):
        return {k: recursively_ignore(v, ignore_keys.get(k, {})) for k, v in dic.items() if ignore_keys.get(k) != True}
    else:
        return {k: recursively_ignore(v, {}) for k, v in dic.items() if k not in ignore_keys}

max_stars_repo_name: chengzhag/DeepPanoContext

max_issues_repo_name: chengzhag/DeepPanoContext

max_forks_repo_name: chengzhag/DeepPanoContext


text:
"""Prepare train and valid for V2E MVSEC recording.

This script also does some plotting and calculating
statistics.

Author: Yuhuang Hu
Email : yuhuang.hu@ini.uzh.ch
"""

import argparse
import glob
import os
import pickle

import cv2
import h5py
import numpy as np
import matplotlib.pyplot as plt

from v2ecore.v2e_utils import video_writer
from v2e_exps.utils import expandpath
from v2e_exps.utils import find_nearest
from v2e_exps.utils import make_dvs_frame
from v2e_exps.compute_voxel import events_to_voxel_grid

parser = argparse.ArgumentParser()

parser.add_argument("--v2e_mvsec_data_path", type=expandpath)
parser.add_argument("--ori_mvsec_data_path", type=expandpath)
parser.add_argument("--mvsec_frame_data_path", type=expandpath)
parser.add_argument("--v2e_inds_save_path", type=expandpath)
parser.add_argument("--show", action="store_true")
parser.add_argument("--out_video_path", type=expandpath, default="")
parser.add_argument("--frame_rate", type=int, default=43)

# data generation params
parser.add_argument("--data_output_root", type=expandpath)
parser.add_argument("--num_events", type=int, default=25000)
parser.add_argument("--num_slices", type=int, default=10)

args = parser.parse_args()

# open v2e mvsec data
v2e_data = h5py.File(args.v2e_mvsec_data_path, "r")
v2e_events = v2e_data["events"]

v2e_time = (v2e_events[-1, 0]-v2e_events[0, 0])/1e6

print("V2E MVSEC Number of Events: {}".format(v2e_events.shape[0]))
print("V2E MVSEC Duration: {}s".format(v2e_time))

# open mvsec data
mvsec_data = h5py.File(args.ori_mvsec_data_path, "r")
mvsec_events = mvsec_data["davis"]["left"]["events"]
mvsec_img_inds = mvsec_data["davis"]["left"][
    "image_raw_event_inds"][()]

mvsec_time = (mvsec_events[-1, 2]-mvsec_events[0, 2])

print("Original MVSEC Number of Events: {}".format(mvsec_events.shape[0]))
print("Original MVSEC Duration: {}s".format(mvsec_time))

# mvsec frame time
frame_ts = mvsec_data["davis"]["left"]["image_raw_ts"][()]
frame_ts -= frame_ts[0]  # convert into ms
frame_ts *= 1e6

# select time, ignore the ones that are larger than total duration
# TODO: WARNING: FOR NIGHT ONLY
frame_ts = np.array([ts for ts in frame_ts if ts < v2e_time*1e6])

frame_list = sorted(
    glob.glob(
        os.path.join(args.mvsec_frame_data_path, "*.png")))

print("Number of MVSEC Frames: {}".format(len(frame_list)))
print("Number of MVSEC Frame Ts: {}".format(len(frame_ts)))

# build up a event index array
v2e_event_inds = []
out_writer = None

start_idx = 0
# skipping at least the first frame since it does not represent
# any events
skip_start = 1

for idx, frame_t in enumerate(frame_ts[skip_start:-1]):
    # TODO: WARNING: FOR NIGHT ONLY
    # skip for adapting frame rate
    if idx % 2 != 0:
        continue

    if not os.path.isfile(args.v2e_inds_save_path) or args.show:
        nearest_value_idx = find_nearest(
            v2e_events, start_idx, frame_t)

        start_idx = nearest_value_idx

        frame_path = frame_list[idx+skip_start]
        v2e_event_inds.append([frame_path, nearest_value_idx, frame_t])

        print("Frame {}:".format(idx+skip_start))
        print("Frame Ts: {}".format(frame_t))
        print("Event TS: {}".format(v2e_events[nearest_value_idx, 0]))
    else:
        print("Index calculated before.")
        break

    if args.show:

        aps_img = cv2.imread(frame_path)

        # V2E Events
        dvs_events = v2e_events[
            max(0, nearest_value_idx-10000):
            min(nearest_value_idx+10000, v2e_events.shape[0])][()]

        dvs_img = make_dvs_frame(
            dvs_events, height=260, width=346, color=False, clip=3)

        # original MVSEC events
        ori_ind = mvsec_img_inds[idx+skip_start]
        ori_events = mvsec_events[
            max(0, ori_ind-10000):
            min(ori_ind+10000, mvsec_events.shape[0])][()]
        ori_events = ori_events[:, [2, 0, 1, 3]]

        ori_dvs_img = make_dvs_frame(
            ori_events, height=260, width=346, color=False, clip=3)

        fig = plt.figure(figsize=(15, 5))
        fig.add_subplot(3, 3, (1, 4))
        plt.imshow(aps_img)
        plt.title("APS Frame")
        plt.axis("off")
        fig.add_subplot(3, 3, (2, 5))
        plt.imshow(dvs_img, cmap="gray")
        plt.title("V2E MVSEC DVS Frame")
        plt.axis("off")
        fig.add_subplot(3, 3, (3, 6))
        plt.imshow(ori_dvs_img, cmap="gray")
        plt.title("Original MVSEC DVS Frame")
        plt.axis("off")

        fig.add_subplot(3, 3, (7, 9))
        ts = frame_ts[:idx+skip_start+1]/1e6
        y = np.zeros_like(ts)
        plt.plot(ts, y, color="blue",
                 linewidth=5)
        plt.plot(frame_ts[idx+skip_start]/1e6, 0, color="red",
                 linewidth=5, marker="o", markersize=10)
        plt.xlabel("Time (s)")
        plt.xlim([frame_ts[skip_start]/1e6, frame_ts[-1]/1e6])
        plt.yticks([])
        plt.tight_layout(pad=0.5)

        fig.canvas.draw()

        data = np.frombuffer(fig.canvas.tostring_rgb(), dtype=np.uint8)
        data = data.reshape(fig.canvas.get_width_height()[::-1] + (3,))

        if out_writer is None:
            out_writer = video_writer(
                args.out_video_path, data.shape[0], data.shape[1],
                frame_rate=args.frame_rate)

        out_writer.write(data)

        #  cv2.imshow("V2E MVSEC Show", data)
        #  cv2.waitKey(1)

        plt.cla()
        fig.clear()
        plt.close()
        del fig

if out_writer is not None:
    out_writer.release()

# save if necessary
if not os.path.isfile(args.v2e_inds_save_path):
    with open(args.v2e_inds_save_path, "wb") as f:
        pickle.dump(v2e_event_inds, f)

if not args.show:
    with open(args.v2e_inds_save_path, "rb") as f:
        v2e_event_inds = pickle.load(f)

    # prepare the training pairs, no need testing pairs
    num_v2e_events = v2e_events.shape[0]
    num_frame_generated = 0

    if not os.path.isdir(args.data_output_root):
        os.makedirs(args.data_output_root)

    for frame_idx, frame_item in enumerate(v2e_event_inds[1:-1]):
        # adding this because we skipped the first frame.
        true_frame_idx = frame_idx+1
        frame_path, frame_ind, frame_t = frame_item

        ev_vol_start_idx = max(frame_ind-args.num_events//2, 0)
        ev_vol_end_idx = min(frame_ind+args.num_events//2, num_v2e_events)

        # condition on too few events
        if (ev_vol_end_idx - ev_vol_start_idx) < args.num_events:
            print("Bad frame {}, not enough events".format(
                  frame_idx))
            continue

        # if the start and the end reached another frame, then we dump this
        # means too few events here
        ev_start_ts = v2e_events[ev_vol_start_idx, 0]
        ev_end_ts = v2e_events[ev_vol_end_idx, 0]

        curr_frame_ts = frame_t
        pre_frame_ts = v2e_event_inds[true_frame_idx-1][2]
        next_frame_ts = v2e_event_inds[true_frame_idx+1][2]

        pre_frame_ind = v2e_event_inds[true_frame_idx-1][1]
        next_frame_ind = v2e_event_inds[true_frame_idx+1][1]

        #  print("Pre: ", frame_ind-pre_frame_ind)
        #  print("Next: ", next_frame_ind-frame_ind)
        #  print(pre_frame_ind, frame_ind, next_frame_ind)
        #  if (next_frame_ind-pre_frame_ind) < args.num_events:
        if ev_start_ts < pre_frame_ts or ev_end_ts > next_frame_ts:
            print("Bad frame {}, too long".format(
                frame_idx))
            print("Pre ts {} {} {}".format(
                pre_frame_ts, ev_start_ts, ev_start_ts < pre_frame_ts))
            print("Next ts {} {} {}".format(
                next_frame_ts, ev_end_ts, ev_end_ts > pre_frame_ts))
            print("Num events: {}".format(next_frame_ind-pre_frame_ind))
            continue

        # select events
        candidate_events = v2e_events[
            ev_vol_start_idx:ev_vol_end_idx][()].astype("float32")
        candidate_events[:, 0] /= 1e6

        # get voxel grid
        voxel = events_to_voxel_grid(
            candidate_events.copy(), num_bins=args.num_slices,
            width=346, height=260)

        # normalization
        nonzero_ev = (voxel != 0)
        num_nonzeros = nonzero_ev.sum()

        # has events
        if num_nonzeros > 0:
            mean = voxel.sum()/num_nonzeros
            stddev = np.sqrt((voxel**2).sum()/num_nonzeros-mean**2)
            mask = nonzero_ev.astype("float32")
            voxel = mask*(voxel-mean)/stddev

        # output path
        output_path = os.path.join(
            args.data_output_root, "v2e_frame_ev_pair_{:05d}".format(
                num_frame_generated+1))

        # frame
        frame = cv2.imread(frame_path, cv2.IMREAD_GRAYSCALE)

        # save data
        np.savez(output_path+".npz", ev_img=voxel,
                 img=frame)

        num_frame_generated += 1

        print("Generated data for frame {} at {}".format(
            frame_idx, num_frame_generated))


text:
import math

import numpy as np
import pytest
from qcodes.utils.validators import Ints

from .conftest import AClass, a_func

ints = [0, 1, 10, -1, 100, 1000000, int(-1e15), int(1e15),
        # warning: True==1 and False==0 - we *could* prohibit these, using
        # isinstance(v, bool)
        True, False,
        # np scalars
        np.int64(3)]
not_ints = [0.1, -0.1, 1.0, 3.5, -2.3e6, 5.5e15, 1.34e-10, -2.5e-5,
            math.pi, math.e, '', None, float("nan"), float("inf"),
            -float("inf"), '1', [], {}, [1, 2], {1: 1}, b'good',
            AClass, AClass(), a_func]

def test_unlimited():
    n = Ints()
    n.validate(n.valid_values[0])

    for v in ints:
        n.validate(v)

    for v in not_ints:
        with pytest.raises(TypeError):
            n.validate(v)

def test_min():
    for min_val in ints:
        n = Ints(min_value=min_val)
        for v in ints:
            if v >= min_val:
                n.validate(v)
            else:
                with pytest.raises(ValueError):
                    n.validate(v)

    for v in not_ints:
        with pytest.raises(TypeError):
            n.validate(v)

def test_max():
    for max_val in ints:
        n = Ints(max_value=max_val)
        for v in ints:
            if v <= max_val:
                n.validate(v)
            else:
                with pytest.raises(ValueError):
                    n.validate(v)

    for v in not_ints:
        with pytest.raises(TypeError):
            n.validate(v)

def test_range():
    n = Ints(0, 10)

    for v in ints:
        if 0 <= v <= 10:
            n.validate(v)
        else:
            with pytest.raises(ValueError):
                n.validate(v)

    for v in not_ints:
        with pytest.raises(TypeError):
            n.validate(v)

    assert repr(n) == '<Ints 0<=v<=10>'
    assert n.is_numeric

def test_failed_numbers():
    with pytest.raises(TypeError):
        Ints(1, 2, 3)

    with pytest.raises(TypeError):
        Ints(1, 1)  # min >= max

    for val in not_ints:
        with pytest.raises((TypeError, OverflowError)):
            Ints(max_value=val)

        with pytest.raises((TypeError, OverflowError)):
            Ints(min_value=val)

def test_valid_values():
    val = Ints()
    for vval in val.valid_values:
        val.validate(vval)

max_forks_repo_name: nikhartman/Qcodes


text:
from typing import Union

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd

from . import stats

def read_csv(filename: str) -> dict:
    raw_results = {"Experimental": {}, "Calculated": []}
    expt_block = False
    calc_block = False
    with open(filename, "r") as f:
        for line in f:
            if "Experiment" in line:
                expt_block = True
            if "Calculate" in line or "Relative" in line:
                expt_block = False
                calc_block = True
            if expt_block and len(line.split(",")) == 3 and line[0] != "#":
                expt = ExperimentalResult(*line.split(","))
                raw_results["Experimental"][expt.ligand] = expt
            if calc_block and len(line.split(",")) == 5 and line[0] != "#":
                calc = RelativeResult(*line.split(","))
                raw_results["Calculated"].append(calc)
    return raw_results

class RelativeResult(object):
    def __init__(self, ligandA, ligandB, calc_DDG, mbar_error, other_error):
        self.ligandA = str(ligandA).strip()
        self.ligandB = str(ligandB).strip()
        # scope for an experimental dDDG?
        self.calc_DDG = float(calc_DDG)
        self.mbar_dDDG = float(mbar_error)
        self.other_dDDG = float(other_error)

        self.calc_dDDG = self.mbar_dDDG + self.other_dDDG  # is this definitely always additive?

    def toDF(self):
        # TODO - can we do the handling of the dataframe in a different way?
        # Or inside the plotting function that needs it?
        return pd.DataFrame(
            {
                "ligandA": self.ligandA,
                "ligandB": self.ligandB,
                "calc_DDG": self.calc_DDG,
                "mbar_dDDG": self.mbar_dDDG,
                "other_dDDG": self.other_dDDG,
                "calc_dDDG": self.calc_dDDG,
            },
            index=[f"{self.ligandA}_{self.ligandB}"],
        )

class ExperimentalResult(object):
    def __init__(self, ligand, expt_DG, expt_dDG):
        self.ligand = ligand
        self.DG = float(expt_DG)
        self.dDG = float(expt_dDG.strip("\n"))

class FEMap(object):
    def __init__(self, csv):
        self.results = read_csv(csv)
        self.graph = nx.DiGraph()
        self.n_edges = len(self.results)

        self.generate_graph_from_results()

    def generate_graph_from_results(self):
        self._name_to_id = {}
        self._id_to_name = {}
        id = 0
        for result in self.results["Calculated"]:
            if result.ligandA not in self._name_to_id.keys():
                self._name_to_id[result.ligandA] = id
                self._id_to_name[id] = result.ligandA
                id += 1
            if result.ligandB not in self._name_to_id.keys():
                self._name_to_id[result.ligandB] = id
                self._id_to_name[id] = result.ligandB
                id += 1
            # # TODO need some exp error for mle to converge for exp... this is a horrible hack
            self.graph.add_edge(
                self._name_to_id[result.ligandA],
                self._name_to_id[result.ligandB],
                calc_DDG=result.calc_DDG,
                calc_dDDG=result.calc_dDDG,
            )

        for node in self.graph.nodes(data=True):
            name = self._id_to_name[node[0]]
            node[1]["name"] = name
            node[1]["exp_DG"] = self.results["Experimental"][name].DG
            node[1]["exp_dDG"] = self.results["Experimental"][name].dDG

        for edge in self.graph.edges(data=True):
            DG_A = self.graph.nodes[edge[0]]["exp_DG"]
            DG_B = self.graph.nodes[edge[1]]["exp_DG"]
            edge[2]["exp_DDG"] = DG_B - DG_A
            dDG_A = self.graph.nodes[edge[0]]["exp_dDG"]
            dDG_B = self.graph.nodes[edge[1]]["exp_dDG"]
            edge[2]["exp_dDDG"] = (dDG_A ** 2 + dDG_B ** 2) ** 0.5

        self.n_ligands = self.graph.number_of_nodes()
        self.degree = self.graph.number_of_edges() / self.n_ligands

        # check the graph has minimal connectivity
        self.check_weakly_connected()
        if not self.weakly_connected:
            print("Graph is not connected enough to compute absolute values")
        else:
            self.generate_absolute_values()

    def check_weakly_connected(self):
        undirected_graph = self.graph.to_undirected()
        self.weakly_connected = nx.is_connected(undirected_graph)
        return nx.is_connected(undirected_graph)

    def generate_absolute_values(self):
        # TODO this could work if either relative or absolute expt values are provided
        if self.weakly_connected:
            f_i_calc, C_calc = stats.mle(self.graph, factor="calc_DDG")
            variance = np.diagonal(C_calc)
            for i, (f_i, df_i) in enumerate(zip(f_i_calc, variance ** 0.5)):
                self.graph.nodes[i]["calc_DG"] = f_i
                self.graph.nodes[i]["calc_dDG"] = df_i

    def draw_graph(self, title: str = "", filename: Union[str, None] = None):
        plt.figure(figsize=(10, 10))
        self._id_to_name = {}
        for i, j in self._name_to_id.items():
            self._id_to_name[j] = i
        nx.draw_circular(self.graph, labels=self._id_to_name, node_color="hotpink", node_size=250)
        long_title = f"{title} \n Nedges={self.n_edges} \n Nligands={self.n_ligands} \n Degree={self.degree:.2f}"
        plt.title(long_title)
        if filename is None:
            plt.show()
        else:
            plt.savefig(filename, bbox_inches="tight")

max_stars_repo_name: glass-w/openff-arsenic

max_issues_repo_name: glass-w/openff-arsenic

max_forks_repo_name: glass-w/openff-arsenic

reason: import numpy,import networkx