text:
import numpy as np
import tables
import glob

from scipy import ndimage,stats
from scipy.optimize import curve_fit

import matplotlib.cm as pltcm
import matplotlib.pyplot as plt
#from matplotlib.colors import LogNorm
from matplotlib.backends.backend_pdf import PdfPages

from optparse import OptionParser

import parula #matlab colormap used in previous foil plots
pltcm.register_cmap(name="parula",cmap=parula.parula_map);

from os import get_terminal_size

hg = [{},{}]; #histograms
ad = [{},{}]; #data arrays

gnames = ["inner","outer","blocked","defect","etching"];

convf = 4.34; #pixel to um

opt = OptionParser(usage="%prog [options]",version="%prog 0.1",description="Neo-TPCQA GEM foil plotter / database map generator.");
opt.option_list[0].help = "Print this help screen and exit";
opt.add_option("-S","--S-dir",default=None,metavar="DIR",action="store",type="string",dest="sdir",help="Source directory for the segmented side input files.");
opt.add_option("-U","--U-dir",default=None,metavar="DIR",action="store",type="string",dest="udir",help="Source directory for the unsegmented side input files.");
opt.add_option("-t","--type",default=None,metavar="TYPE",action="store",type="string",dest="type",help="Force foil type to TYPE=IROC/O1/O2/O3. By default, type is recognized from the source directory name.");
opt.add_option("-o","--output-dir",default="./",metavar="DIR",action="store",type="string",dest="outdir",help="Destination directory for the output files (report and maps) [default: %default]");
opt.add_option("-n","--name",default=None,metavar="NAME",action="store",type="string",dest="name",help="Name for the generated files. If unspecified, an attempt will be made to acquire the name from segmented side source path.");
opt.add_option("-c","--color-map",default="parula",metavar="CMAP",type="string",dest="colormap",help="Name of the colormap to use from the list of matplotlib colormaps. [default: %default]");
opt.add_option("-q","--quiet",action="store_true",dest="quiet",help="Quiet mode. Do not show plots.");
(options,args) = opt.parse_args();

if options.sdir is None:
	opt.error("S-side source unspecified.");
if options.udir is None:
	opt.error("U-side source unspecified.");

class ProgressBar():
	def __init__(self, l = None):
		if l is None:
			try:
				ts,_ = get_terminal_size();
				self.l = ts-12;
			except OSError:
				self.l = 50;
		else:
			self.l = l;
	def __enter__(self):
		return self;
	def update(self, r):
		a = np.ceil(r*self.l);
		b = np.floor(self.l-r*self.l);
		print("{:.2f}%\t|{}>{}|".format(100*r,int(a)*'=',int(b)*'-'),end='\r');
	def __exit__(self, type, val, tb):
		print("{:.2f}%\t|{}=|".format(100.0,self.l*'='));

srcdir = options.sdir.replace('\\','/'); #take the foil name from S-source
while srcdir[-1] == '/':
	srcdir = srcdir[:-1];
pathc = srcdir.split('/');

foiln = options.name;
if foiln is None:
	for p in reversed(pathc):
		if any([s in p for s in
			("I_","I-","O1_","O1-","O2_","O2-","O3_","O3-")]):
			foiln = p;
			break;
	else:
		opt.error("Unable to extract foil name. Use -n to manually specify the name.");

	for s in ["_S","_U","-s","-u","-"]:
		try:
			u = foiln.index(s);
			foiln = foiln[:u];
			break;
		except ValueError:
			pass;
if len(foiln) == 0:
	foiln = "output";

t = options.type;
if t is None:
	if any([s in foiln for s in ("I_","I-")]):
		t = "IROC";
	else:
		o = foiln.index('O');
		t = foiln[o:o+2];

try:
	nl = {
		"IROC":np.array([467.0,496.5]),
		"O1":np.array([595.8,353.0]),
		"O2":np.array([726.2,350.0]),
		"O3":np.array([867.0,379.0])
	}[t];
	el = {"IROC":88,"O1":63,"O2":63,"O3":63}[t];

	nb = [int(np.ceil(n)) for n in nl];
	print("Foil/{}: assuming type {}.".format(foiln,t));
except KeyError:
	opt.error("Unrecognized foil type / invalid option -t.");

def groupCounter(hf):
	for g in hf.walk_groups():
		yield g;

print("Loading sources...");
for u,srcdir in enumerate([options.sdir,options.udir]):
	print("[{}]".format(srcdir));
	with ProgressBar() as pb:
		hfls = glob.glob(srcdir+"/*.h5");
		if len(hfls) == 0:
			opt.error("No source files (*.h5) found.");

		for i,hfn in enumerate(hfls):
			hf = tables.open_file(hfn,"r");

			for N,_ in enumerate(groupCounter(hf)):
				pass;

			for j,g in enumerate(hf.walk_groups()):
				try:
					data = g.data.read();

					a = ad[u].get(g._v_name);
					if a is None:
						a = {"x":np.array([]),"y":np.array([]),"d":np.array([]),"fl":np.array([])};
						ad[u][g._v_name] = a;
					a["x"] = np.concatenate((a["x"],data[:,0]));
					a["y"] = np.concatenate((a["y"],data[:,1]));
					a["d"] = np.concatenate((a["d"],data[:,7]));
					a["fl"] = np.concatenate((a["fl"],data[:,24]));

				except tables.exceptions.NoSuchNodeError:
					pass;

				pb.update((i+j/N)/len(hfls));

			hf.close();

print("Booking...");
with ProgressBar() as pb:
	for u in range(0,2):
		for i,name in enumerate(gnames):
			a = ad[u][name];
			x = a["x"];
			y = a["y"];

			pb.update(0.5*i/len(gnames)+0.5*u);

			he = {};
			hg[u][name] = he;
			try:
				he["count"],_,_ = np.histogram2d(x,y,bins=[edgesx,edgesy]);
			except NameError:
				he["count"],edgesx,edgesy = np.histogram2d(x,y,bins=nb);

			he["diam"],_,_,_ = stats.binned_statistic_2d(x,y,a["d"]*convf,"mean",bins=[edgesx,edgesy]);
			he["diam"][np.isnan(he["diam"])] = 0;
			he["dstd"],_,_,_ = stats.binned_statistic_2d(x,y,a["d"]*convf,"std",bins=[edgesx,edgesy]);
			he["dstd"][np.isnan(he["dstd"])] = 0;

			he["fl"],_,_,_ = stats.binned_statistic_2d(x,y,a["fl"],"mean",bins=[edgesx,edgesy]);
			he["fl"][np.isnan(he["fl"])] = 0;
			he["flstd"],_,_,_ = stats.binned_statistic_2d(x,y,a["fl"],"std",bins=[edgesx,edgesy]);
			he["flstd"][np.isnan(he["flstd"])] = 0;

		he = {};
		hg[u]["rim"] = he;

		he["diam"] = 0.5*(hg[u]["outer"]["diam"]-hg[u]["inner"]["diam"]);
		he["diam"][np.isnan(he["diam"])] = 0.0;
		he["dstd"] = np.sqrt(0.25*hg[u]["outer"]["dstd"]**2+0.25*hg[u]["inner"]["dstd"]**2);
		he["dstd"][np.isnan(he["dstd"])] = 0.0;

def renderFoilPlot(fig, ax, h, title, cdfrange = True, zrnan = False):
	if cdfrange:
		hf = h.flatten();
		#estimate optimal plotting range by removing peaks
		ph,bins = np.histogram(hf[hf > 0],1000,density=True);
		cdf = np.cumsum(ph)*(bins[1]-bins[0]);
		a = np.argmax(cdf > 0.05);
		b = np.argmax(cdf > 0.95);
		umin,umax = bins[a],bins[b];
	else:
		umin,umax = None,None;

	if zrnan:
		h = h.copy();
		h[h < 1e-6] = np.nan;

	ax.set_title(title,fontsize=8);
	#ax.set_xlabel("x (mm)",fontsize=6);
	#ax.set_ylabel("y (mm)",fontsize=6);
	cmap = ax.imshow(h,interpolation="nearest",cmap=options.colormap,extent=(0,nl[0],0,nl[1]),
		vmin=umin,vmax=umax);
	ax.plot([0,el],[0,nl[1]],color="black");
	ax.plot([nl[0],nl[0]-el],[0,nl[1]],color="black");
	[t.label.set_fontsize(6) for t in ax.xaxis.get_major_ticks()];
	[t.label.set_fontsize(6) for t in ax.yaxis.get_major_ticks()];
	cbar = fig.colorbar(cmap,ax=ax);
	cbar.ax.tick_params(labelsize=6);

	ax.autoscale(False);

def renderPointOverlay(ax, x, y):
	x = nl[1]*(x-edgesx[0])/(edgesx[-1]-edgesx[0]);
	y = nl[0]*(y-edgesy[0])/(edgesy[-1]-edgesy[0]);
	ax.plot(y,nl[1]-x,marker="o",linestyle="none",color="red",mfc="none");

def gauss(x, *p):
	A,mu,sigma = p;
	return A*np.exp(-(x-mu)**2/(2.0*sigma*sigma));

options.outdir += "/";
print("Generating output... (dst: {})".format(options.outdir));
with PdfPages(options.outdir+foiln+".pdf") as pdf:

	for fi in [
		[("inner","diam","Inner Diameter (S)"),("inner","diam","Inner Diameter (U)"),
		("inner","dstd","Inner Diameter Std. Dev (S)"),("inner","dstd","Inner Diameter Std. Dev (U)"),
		("outer","diam","Outer Diameter (S)"),("outer","diam","Outer Diameter (U)"),
		("outer","dstd","Outer Diameter Std. Dev (S)"),("outer","dstd","Outer Diameter Std. Dev (U)")],

		[("inner","count","Inner N (S)"),("inner","count","Inner N (U)"),
		("defect","count","Defect N (S)",False),("defect","count","Defect N (U)",False),
		("etching","count","Etching N (S)",False),("etching","count","Etching N (U)",False),
		("blocked","count","Blocked N (S)",False),("blocked","count","Blocked N (U)",False)],

		[("rim","diam","Rim Diameter (S)"),("rim","diam","Rim Diameter (U)"),
		("rim","dstd","Rim Diameter Std. Dev (S)"),("rim","dstd","Rim Diameter Std. Dev (U)"),
		("inner","fl","Foreground Light (S)"),("inner","fl","Foreground Light (U)"),
		("inner","flstd","Foreground Light Std. Dev (S)"),("inner","flstd","Foreground Light Std. Dev (U)")]
		]:
		fig,ax = plt.subplots(4,2,figsize=(0.707*8,8));
		for i,r in enumerate(fi):
			h = hg[i%2][r[0]][r[1]];
			q = np.unravel_index(i,(4,2));
			b = r[3] if len(r) > 3 else True;
			renderFoilPlot(fig,ax[q],h,r[2],b,not b);
			if not b:
				renderPointOverlay(ax[q],ad[i%2][r[0]]["x"],ad[i%2][r[0]]["y"]);
		plt.savefig(pdf,format="pdf");

	#m = np.fromfunction(lambda i,j: np.vectorize(Get)(i,j),(x,y),dtype=float);
	#h = hg[name]["diam"];
	#labels,nlabels = ndimage.label(h > 1e-6);
	#labelr = ndimage.sum(h,labels,range(nlabels+1));

	#h[(labelr < 1000)[labels]] = 0.0;
	#h = h[~np.all(h == 0,axis=1)][:,~np.all(h == 0,axis=0)];

	#remove disconnected artifacts
	#h[(labelr < 100)[labels]] = np.nan;
	#q1 = np.bitwise_or(np.isnan(h),h == 0.0);
	#q2 = np.bitwise_and(~np.all(q1,axis=1),np.any(np.isnan(h),axis=1));
	#q3 = np.bitwise_and(~np.all(q1,axis=0),np.any(np.isnan(h),axis=0));
	#h = h[q2][:,q3];
	#h[np.isnan(h)] = 0.0;

	#remove connected artifacts
	#h = h[~np.bitwise_and(np.any(h > 15,axis=1),np.count_nonzero(h == 0.0,axis=1)/y > 0.9)];
	#h = h[:,~np.bitwise_or(
	#	np.bitwise_and(np.any(h > 15,axis=0),np.count_nonzero(h == 0.0,axis=0)/x > 0.9),
	#	np.all(h == 0.0,axis=0))];
	#while True:
	#	n = len(h);
	#	h = h[~np.bitwise_and(h[:,0] > 1e-6,np.count_nonzero(h == 0.0,axis=1)/y > 0.9)];
	#	if len(h) == n:
	#		break;
	#	h = h[:,~np.all(h == 0.0,axis=0)];

	yi,bins_i = [None]*2,[None]*2;
	yo,bins_o = [None]*2,[None]*2;
	yr,bins_r = [None]*2,[None]*2;
	fig,ax = plt.subplots(2,1,sharey=True,figsize=(0.707*8,8));
	#fig.subplots_adjust(wspace=0.0,hspace=0.0);
	for u in range(0,2):
		yi[u],bins_i[u],_ = ax[u].hist(ad[u]["inner"]["d"]*convf,300,range=(0,110),alpha=0.5);
		yo[u],bins_o[u],_ = ax[u].hist(ad[u]["outer"]["d"]*convf,300,range=(0,110),alpha=0.5);

		rh = 0.5*(hg[u]["outer"]["diam"]-hg[u]["inner"]["diam"]); #rim histogram
		rh[np.isnan(rh)] = 0.0;
		rh = rh.flatten();

		yr[u],bins_r[u] = np.histogram(rh,300,range=(1,110));
		w = min(np.max(yi[u]),np.max(yo[u]))/np.max(yr[u]);
		#ax.bar(0.5*(bins[:-1]+bins[1:]),y*w,align="center",color="black",alpha=0.3);
		ax[u].hist(rh,300,range=(1,110),weights=np.full(len(rh),w),alpha=0.5);

		ax[u].ticklabel_format(style='sci',axis='y',scilimits=(0,0))
		ax[u].set_title("Diameter Profile ("+["S","U"][u]+")");

	ax[1].set_xlabel("Diameter ($\\mathrm{\\mu m}$)");
	#ax.set_ylabel("Occurrence");
	plt.savefig(pdf,format="pdf");

	#binc = 0.5*(bins[:-1]+bins[1:]);
	#binc = binc[y > 10];
	#y = y[y > 10];
	#coeff,cov = curve_fit(gauss,binc,y,p0=[1,0,1000]);
	#ax.plot(binc,t);
	#ax.plot(binc,gauss(binc,*coeff));

binc = 0.5*(bins_i[0][:-1]+bins_i[0][1:]);
with open(options.outdir+foiln+"_profile.txt","w") as pf:
	for i,_ in enumerate(binc):
		pf.write("\t".join(["{:4f}".format(f).rstrip('0').rstrip('.') for f in [
			binc[i],yr[0][i],yr[1][i],yi[0][i],yi[1][i],yo[0][i],yo[1][i]
		]])+"\n");

with open(options.outdir+foiln+"_map.txt","w") as mf:
	for (j,i),_ in np.ndenumerate(np.flip(hg[0]["inner"]["count"],axis=1)):
		mf.write("\t".join(["{:.4f}".format(f).rstrip('0').rstrip('.') for f in [
			float(i)*nl[0]/nb[0]+0.5/nl[0],float(j)*nl[1]/nb[1]+0.5/nl[1],
			hg[0]["inner"]["diam"][j][i],hg[0]["outer"]["diam"][j][i],
			hg[0]["blocked"]["diam"][j][i],hg[0]["defect"]["diam"][j][i],hg[0]["etching"]["diam"][j][i],
			hg[1]["inner"]["diam"][j][i],hg[1]["outer"]["diam"][j][i],
			hg[1]["blocked"]["diam"][j][i],hg[1]["defect"]["diam"][j][i],hg[1]["etching"]["diam"][j][i],
			hg[0]["inner"]["count"][j][i],hg[0]["outer"]["count"][j][i],
			hg[0]["blocked"]["count"][j][i],hg[0]["defect"]["count"][j][i],hg[0]["etching"]["count"][j][i],
			hg[1]["inner"]["count"][j][i],hg[1]["outer"]["count"][j][i],
			hg[1]["blocked"]["count"][j][i],hg[1]["defect"]["count"][j][i],hg[1]["etching"]["count"][j][i],
			hg[0]["inner"]["dstd"][j][i],hg[0]["outer"]["dstd"][j][i],hg[1]["inner"]["dstd"][j][i],hg[1]["outer"]["dstd"][j][i],
			hg[0]["inner"]["fl"][j][i],hg[1]["inner"]["fl"][j][i]
		]])+"\n");

if not options.quiet:
	plt.show();

max_stars_repo_name: jaelpark/neogem-TPCQA

max_issues_repo_name: jaelpark/neogem-TPCQA

max_forks_repo_name: jaelpark/neogem-TPCQA


text:
import gym
import numpy as np
from utils.Memory import Supervised_ExperienceReplay
from configs import StudentPongConfig as student_config
from Pacman.processor import process_observation, process_state_batch, init_state, append_frame
from copy import deepcopy

def accumulate_experience(teacher, exp_replay: Supervised_ExperienceReplay, config=student_config):
    """
    teacher feeds the Experience replay with new experiences
    :param teacher: teacher net, knows how to solve the problem
    :param exp_replay: the experience replay where the teacher saves its experiences
    :param config : holds customer variables such as OBSERVE
    :return: an experience replay filled with new experiences
    """
    env = gym.make('MsPacmanDeterministic-v4')
    steps = 0
    while 1:
        state = init_state()
        observation = env.reset()
        observation = process_observation(observation)
        state = np.asarray(state)
        done = False
        while not done:
            steps += 1
            state = append_frame(state, observation)
            teacher_q_value = teacher.get_q(state=np.expand_dims(state, axis=0))
            action = teacher.select_action(teacher_q_value)
            next_observation, reward, done, _ = env.step(action)
            next_observation = process_observation(next_observation)
            observation = next_observation
            exp_replay.add_memory(state, teacher_q_value, action)  # feeding the experience replay
        if steps > config.OBSERVE:  # we have OBSERVE  number of exp in exp_replay
            try:
                del env
            except ImportError:
                pass
            break

ext: py

max_stars_repo_name: dorlivne/PoPS

max_issues_repo_name: dorlivne/PoPS

max_forks_repo_name: dorlivne/PoPS

reason: import numpy


text:
import numpy as np
import matplotlib.pyplot as pl
import h5py
import platform
import os
from ipdb import set_trace as stop
from astropy.io import fits
import scipy.io as io
import time
import matplotlib.animation as manimation

os.environ["KERAS_BACKEND"] = "tensorflow"

if (platform.node() != 'vena'):
    os.environ["CUDA_VISIBLE_DEVICES"] = "0"

from keras.layers import Input, Convolution2D, merge, Activation, Lambda, BatchNormalization
from keras.callbacks import ModelCheckpoint, Callback
from keras.models import Model, model_from_json
import tensorflow as tf
import keras.backend.tensorflow_backend as ktf

import poppy
import pyfftw
import scipy.special as sp
from astropy import units as u
import congrid

def _even(x):
    return x%2 == 0

def _zernike_parity(j, jp):
    return _even(j-jp)

def _zernike_coeff_kolmogorov(D, r0, n_zernike):
    """
    Return Zernike coefficients in phase units
    """
    covariance = np.zeros((n_zernike,n_zernike))
    for j in range(n_zernike):
        n, m = poppy.zernike.noll_indices(j+1)
        for jpr in range(n_zernike):
            npr, mpr = poppy.zernike.noll_indices(jpr+1)

            deltaz = (m == mpr) and (_zernike_parity(j, jpr) or m == 0)

            if (deltaz):
                phase = (-1.0)**(0.5*(n+npr-2*m))
                t1 = np.sqrt((n+1)*(npr+1))
                t2 = sp.gamma(14./3.0) * sp.gamma(11./6.0)**2 * (24.0/5.0*sp.gamma(6.0/5.0))**(5.0/6.0) / (2.0*np.pi**2)

                Kzz = t2 * t1 * phase

                t1 = sp.gamma(0.5*(n+npr-5.0/3.0)) * (D / r0)**(5.0/3.0)
                t2 = sp.gamma(0.5*(n-npr+17.0/3.0)) * sp.gamma(0.5*(npr-n+17.0/3.0)) * sp.gamma(0.5*(n+npr+23.0/3.0))
                covariance[j,jpr] = Kzz * t1 / t2

    covariance[0,0] = 1.0

    out = np.random.multivariate_normal(np.zeros(n_zernike), covariance)

    out[0] = 0.0

    return out

class imax_degradation(object):

    def __init__(self, telescope_radius, pixel_size, fov, secondary_radius=None):
        self.telescope_radius = telescope_radius
        if (secondary_radius != None):
            self.secondary_radius = secondary_radius
        else:
            self.secondary_radius = 0.0
        self.pixel_size = pixel_size
        self.fov = fov
        self.zernike_max = 45
        self.r0 = 10.0 * u.cm

    def compute_psf(self, lambda0, rms_telescope=1.0/9.0, rms_atmosphere=1.0/9.0):
        self.lambda0 = lambda0

        osys = poppy.OpticalSystem()

        osys.add_pupil(poppy.CircularAperture(radius = self.telescope_radius))

        if (self.secondary_radius != 0):
            osys.add_pupil(poppy.SecondaryObscuration(secondary_radius = self.secondary_radius))

# Telescope diffraction + aberrations
        phase_telescope = np.random.randn(self.zernike_max)
        sigma = np.sqrt(np.sum(phase_telescope[4:]**2)) / 2.0 / np.pi
        phase_telescope *= rms_telescope / sigma
        phase_telescope[0:4] = 0.0

# Atmosphere
        phase_atmosphere = _zernike_coeff_kolmogorov(2.0 * self.telescope_radius.to(u.cm).value, self.r0.to(u.cm).value, self.zernike_max)
        sigma = np.sqrt(np.sum(phase_atmosphere[4:]**2)) / 2.0 / np.pi
        phase_atmosphere *= rms_atmosphere / sigma
        phase_atmosphere[0:4] = 0.0

        thinlens = poppy.ZernikeWFE(radius=self.telescope_radius.to(u.m).value, coefficients=(phase_telescope + phase_atmosphere) * lambda0.to(u.m).value / (2.0 * np.pi))
        osys.add_pupil(thinlens)

        osys.add_detector(pixelscale=self.pixel_size, fov_pixels=self.fov, oversample=1)

        psf = osys.calc_psf(wavelength=self.lambda0)

        self.psf = psf[0].data
        nx_psf, ny_psf = psf[0].data.shape

        psf = np.roll(self.psf, int(nx_psf/2), axis=0)
        psf = np.roll(psf, int(ny_psf/2), axis=1)

        self.psf_fft = pyfftw.interfaces.numpy_fft.fft2(psf)

    def apply_psf(self, image):
        image_fft = pyfftw.interfaces.numpy_fft.fft2(image)
        self.image = np.real(pyfftw.interfaces.numpy_fft.ifft2(self.psf_fft * image_fft))
        return self.image

    def rebin_image(self, nx, ny):
        return congrid.resample(self.image, (nx, ny))

def running_mean(x, N):
    cumsum = np.cumsum(np.insert(x, 0, 0))
    return (cumsum[N:] - cumsum[:-N]) / N

class trainDNNFull(object):

    def __init__(self, root, output, name_of_variable):

# Only allocate needed memory
        config = tf.ConfigProto()
        config.gpu_options.allow_growth=True
        session = tf.Session(config=config)
        ktf.set_session(session)

        self.root = root
        self.nx = 576
        self.ny = 576
        self.n_times = 2
        self.n_filters = 64
        self.batch_size = 1
        self.n_conv_layers = 20
        self.stride = 1
        self.skip_frequency = 2
        self.n_frames = 1
        self.output = output
        self.name_of_variable = name_of_variable

        telescope_radius = 0.5 * 0.965 * u.meter
        pixel_size = 0.02759 * u.arcsec / u.pixel
        fov = 1152 * u.pixel
        lambda0 = 500 * u.nm
        imax = imax_degradation(telescope_radius, pixel_size, fov)
        imax.compute_psf(lambda0)

        res = io.readsav('/net/viga/scratch1/deepLearning/opticalFlow/mancha/c3d_1152_cont4_4bin_012000_continuum.sav')['continuum']

        self.images = np.zeros((2,576,576), dtype='float32')

# 576 pixels are obtained by resampling 1152 pixels of 0.02759 "/px to 0.0545 "/px for IMaX
        self.images[0,:,:] = congrid.resample(imax.apply_psf(res[0,:,:]), (576, 576))
        self.images[1,:,:] = congrid.resample(imax.apply_psf(res[1,:,:]), (576, 576))

        res = io.readsav('/net/viga/scratch1/deepLearning/opticalFlow/mancha/c3d_1152_cont4_4bin_012000.isotau.sav')

        self.vx = np.zeros((3,576,576), dtype='float32')
        self.vy = np.zeros((3,576,576), dtype='float32')

        for i in range(3):
            self.vx[i,:,:] = congrid.resample(imax.apply_psf(res['vx'][i,:,:]), (576, 576))
            self.vy[i,:,:] = congrid.resample(imax.apply_psf(res['vy'][i,:,:]), (576, 576))

    def residual(self, inputs):
        x = Convolution2D(self.n_filters, 3, 3, border_mode='same', init='he_normal')(inputs)
        x = BatchNormalization()(x)
        x = Activation('relu')(x)
        x = Convolution2D(self.n_filters, 3, 3, border_mode='same', init='he_normal')(x)
        x = BatchNormalization()(x)
        x = merge([x, inputs], 'sum')

        return x

    def defineNetwork(self):
        print("Setting up network...")

        inputs = Input(shape=(self.nx, self.ny, self.n_times))
        conv = Convolution2D(self.n_filters, 3, 3, activation='relu', border_mode='same', init='he_normal')(inputs)

        x = self.residual(conv)
        for i in range(self.n_conv_layers):
            x = self.residual(x)

        x = Convolution2D(self.n_filters, 3, 3, border_mode='same', init='he_normal')(x)
        x = BatchNormalization()(x)
        x = merge([x, conv], 'sum')

        final = Convolution2D(6, 1, 1, activation='linear', border_mode='same', init='he_normal')(x)

        self.model = Model(input=inputs, output=final)

        print("Loading weights...")
        self.model.load_weights("{0}_weights.hdf5".format(self.root))

    def validation_generator(self):

        self.median_i = np.median(self.images[:,100:-100,100:-100])

        input_validation = np.zeros((self.batch_size,self.nx,self.ny,2), dtype='float32')

        while 1:
            for i in range(self.n_frames):

                print('{0}/{1}'.format(i,self.n_frames))

                input_validation[:,:,:,0] = self.images[i*self.batch_size:(i+1)*self.batch_size,:,:] / self.median_i
                input_validation[:,:,:,1] = self.images[i*self.batch_size+1:(i+1)*self.batch_size+1,:,:] / self.median_i

                yield input_validation

    def predict_validation(self):
        print("Predicting validation data...")

        tmp = np.load('/net/duna/scratch1/aasensio/deepLearning/opticalFlow/database/normalization.npz')
        min_i, max_i, min_v, max_v = tmp['arr_0'], tmp['arr_1'], tmp['arr_2'], tmp['arr_3']

        # ff = io.readsav(self.observations)
        # im = ff['cont']

        # x = np.arange(self.nx)
        # y = np.arange(self.ny)

        start = time.time()
        out = self.model.predict_generator(self.validation_generator(), self.n_frames, max_q_size=1)
        end = time.time()

        print("Prediction took {0} seconds...".format(end-start))

        ff = h5py.File(self.output, 'w')
        db = ff.create_dataset("velocity", (self.n_frames, self.nx, self.ny, 6), dtype='float32')

        for i in range(6):
            out[:,:,:,i] = out[:,:,:,i] * (max_v[i] - min_v[i]) + min_v[i]

        db[:] = out[:]

        db = ff.create_dataset("vx", (3, self.nx, self.ny), dtype='float32')
        db[:] = self.vx

        db = ff.create_dataset("vy", (3, self.nx, self.ny), dtype='float32')
        db[:] = self.vy
        ff.close()

    def compare_velocities(self):

        ff = h5py.File(self.output)
        prediction = ff.get("velocity")

        f, ax = pl.subplots(nrows=2, ncols=2, figsize=(10,10))

        im = ax[0,0].imshow(self.vx[0,200:300,200:300] / 1e5)
        pl.colorbar(im, ax=ax[0,0])

        im = ax[0,1].imshow(10*prediction[0,200:300,200:300,0])
        pl.colorbar(im, ax=ax[0,1])

        im = ax[1,0].imshow(self.vy[0,200:300,200:300] / 1e5)
        pl.colorbar(im, ax=ax[1,0])

        im = ax[1,1].imshow(10*prediction[0,200:300,200:300,1])
        pl.colorbar(im, ax=ax[1,1])

        pl.show()

        stop()

if (__name__ == '__main__'):

    # out = trainDNNFull('../training/cnns/resnet', 'cont.idl', 'imax_velocity.h5', 'cont')
    out = trainDNNFull('../../training/cnns/resnet2', 'mancha_velocity.h5', 'mov')
    out.defineNetwork()
    out.predict_validation()
    #out.compare_velocities()

max_stars_repo_name: aasensio/DeepLearning

max_issues_repo_name: aasensio/DeepLearning

max_forks_repo_name: aasensio/DeepLearning

reason: import numpy,import scipy,from astropy


text:
#!/usr/bin/env python
# RMS 2018

# A genetic algorithm to aid in feature selection in machine learning problems

import multiprocessing as mp

import numpy as np
from sklearn.model_selection import train_test_split

class GeneticAlgorithm(object):

    def __init__(self, X, Y, model, Niter=100, keep_fraction=0.5, mutation_rate="auto", nfeatures="auto",
                 test_size=0.3, njobs=1):

        '''

        A simple genetic algorithm designed to work with sklearn objects.

        The object of the algorithm is to select the optimal combination of columns of input datafame X that produce the
        best prediction of Y given the algorithm object provided.

        To do this, the algorithm starts by randomly selecting different
        combinations of columns. These become the first 'generation' of individuals

        For each combination, the sklearn algorithm is trained in the associated columns and tested. The test score
        becomes the 'fitness' of that combination/individual

        Individuals are then selected for 'breeding'. Those with the highest fitness scores have the highest chances
        of doing so. When two individuals 'breed'. Two children are produced by a simple crossover of their column IDs

        The child generation is then combined with the 'best' of the parent generation via the 'keep_fraction' argument

        The algorithm proceeds by testing the fitness of each generation. The result is an 'optimal' combination of
        features that correspond to the best fitness score.

        This should work with any supervised sklearn object

        Inputs:
        X - dataframe containing the predictors only
        Y - datafame or series containing the target only
        model - sklearn classifier or regression object, such as an RandomForestClassifier
        Niter - number of iterations of the genetic algorithm
        keep_fraction - the proportion of fittest parents to keep in each new generation
        mutation_rate - the probability of mutation in each child
        nfeatures - the maximum number of features that an output model can have
        test_size - test_size in the train_test_split that occurs in model fitness evaluation

        Main outputs:
        self.fitness_evolution - list of the best fitness value from each generation
        self.best_individual_evolution - list of arrays of the best individuals in each generation
        self.feature_selection - dataframe corresponding to the selected features
        self.best_fitness - fitness score correspondng to self.feature_selection
        self.best_individual - individual corresponding to self.feature_selection

        '''

        self.dataset = X
        self.response = Y
        self.model = model  # needs to be a sklearn object
        self.Niter = Niter  # number of iterations
        self.parent_keep = keep_fraction
        self.test_size = test_size
        self.nprocs = int(njobs)

        if self.nprocs > mp.cpu_count():
            raise ValueError("Entered number of processes > CPU count!")

        self.feature_columns = self.dataset.columns

        if nfeatures == 'auto':
            self.nfeatures = len(self.feature_columns)
        else:
            self.nfeatures = nfeatures

        self.P = 2 * int(np.ceil(self.nfeatures * 1.5 / 2))  # number of individuals in a given generation

        if mutation_rate == 'auto':

            self.mutation_rate = 1.0 / (self.P * np.sqrt(self.nfeatures))
        else:
            self.mutation_rate = mutation_rate

        self.fitness_evolution = []
        self.best_individual_evolution = []

        # These three things are typically the most desired output
        self.feature_selection = None
        self.best_fitness = None
        self.best_individual = None

    def fitness(self, generation):

        '''
        Assess the fitness of a generation of individuals

        This is the part that takes a long time because it must train a supervised ML algorithm on all individuals in
        a generation
        '''

        def determine_fitness(subgeneration, output, pos):

            fitness_array = np.zeros(np.shape(subgeneration)[0])

            for i in range(np.shape(subgeneration)[0]):
                individual = subgeneration[i, :]

                # Subset the columns based on this individual
                X_individual = self.dataset[
                    [self.dataset.columns[j] for j in range(len(individual)) if individual[j] == 1]]

                # Split into train-test datasets
                X_train, X_test, y_train, y_test = train_test_split(X_individual, self.response,
                                                                    test_size=self.test_size)

                # Fit the classifier
                self.model.fit(X_train, y_train)

                # Report fitness score (score in the testing dataset)
                fitness = self.model.score(X_test, y_test)

                # append to fitness array
                fitness_array[i] = fitness

            output.put((pos, fitness_array))

        process_output = mp.Queue()
        subarrays = np.array_split(generation, self.nprocs)
        processes = [mp.Process(target=determine_fitness, args=(subarrays[i], process_output, i)) for i in
                     range(self.nprocs)]

        for p in processes:
            p.start()

        for p in processes:
            p.join()

        results = [process_output.get() for p in processes]
        results.sort()
        rlist = []
        for element in results:
            r = element[1]
            for j in range(len(r)):
                rlist.append(r[j])
        rlist = np.array(rlist)

        return rlist

    def make_new_generation(self, old_generation, old_fitness_array):

        '''
        Make a new generation of individuals
        '''

        generation_size = len(old_fitness_array)

        # Vector describing the probability of reporduction of each individual in a generation
        prob_weights = 2 * np.argsort(old_fitness_array / (generation_size * (generation_size + 1)))[::-1]

        prob_reproduction = prob_weights / np.sum(prob_weights)

        # Make vector of indices to choose
        a = np.arange(generation_size)

        children = np.zeros([2 * generation_size, np.shape(old_generation)[1]])

        for i in range(generation_size):
            parent_index_pair = np.random.choice(a, size=2, replace=False, p=prob_reproduction)

            parent1 = old_generation[parent_index_pair[0]]
            parent2 = old_generation[parent_index_pair[1]]

            # Do cross over and apply mutation to generate two children for each parent pair
            child1 = parent1.copy()
            child2 = parent2.copy()

            # Generate locations of genetic information to swap
            pos = np.random.choice(len(parent1), size=int(len(parent1) / 2), replace=False)
            child1[pos] = parent2[pos]
            child2[pos] = parent1[pos]

            # Generate mutation vector
            mutate1 = np.random.binomial(1, self.mutation_rate, len(parent1))
            mutate2 = np.random.binomial(1, self.mutation_rate, len(parent1))

            # Generate children and fill child array
            child1 = (child1 + mutate1 >= 1).astype(int)
            child2 = (child2 + mutate2 >= 1).astype(int)

            children[i, :] = child1
            children[-(i + 1), :] = child2

        # shuffle and return only the same number of children as there were parents
        np.random.shuffle(children)

        new_generation = children[0:generation_size, :]

        # replace some fraction of the children with the fittest parents, if desired

        nparents_to_keep = int(self.parent_keep * generation_size)

        if nparents_to_keep > 0:
            parents_keep = np.argsort(old_fitness_array)[::-1][:nparents_to_keep]

            for i in range(len(parents_keep)):
                new_generation[i, :] = old_generation[parents_keep[i], :]

        np.random.shuffle(new_generation)

        return new_generation

    def fit(self):

        '''
        Run the genetic algorithm to obtain the optimal features for this problems
        This part takes a long time and could be parallelized
        '''

        # Make the first generation
        old_generation = np.zeros([self.P, self.nfeatures])
        for i in range(self.P):
            old_generation[i, :] = np.random.binomial(1, 0.5, self.nfeatures)

        old_fitness_array = self.fitness(old_generation)

        self.best_fitness = np.max(old_fitness_array)
        self.best_individual = old_generation[np.argmax(old_fitness_array), :]

        self.fitness_evolution.append(self.best_fitness)
        self.best_individual_evolution.append(self.best_individual)

        for n in range(1, self.Niter):
            print("GeneticAlgorithm: Testing generation %i" % n)

            # Make new generation
            new_generation = self.make_new_generation(old_generation, old_fitness_array)
            # Get fitness of new generation
            new_fitness_array = self.fitness(new_generation)

            # Locate and extract the best individual and its score
            self.best_fitness = np.max(new_fitness_array)
            self.best_individual = new_generation[np.argmax(new_fitness_array), :]
            self.fitness_evolution.append(self.best_fitness)
            self.best_individual_evolution.append(self.best_individual)

            old_fitness_array = new_fitness_array
            old_generation = new_generation

        # Get the features associated with the 'winning' individual

        self.feature_selection = self.dataset[
            [self.dataset.columns[j] for j in range(len(self.best_individual)) if self.best_individual[j] == 1]]

max_stars_repo_name: rmartinshort/GeneticAlgorithm

max_issues_repo_name: rmartinshort/GeneticAlgorithm

max_forks_repo_name: rmartinshort/GeneticAlgorithm


text:
import numpy as np

# V1 = np.random.uniform(low=-10, high=10, size=(5, 3))
# print(type(V1))
# print(V1)
#
# print("------------------------------")
# V2 = V1.tolist()
# print(type(V2))
# del V2[0]
# print(V2)
# V1 = np.array(V2)
# print(type(V1))
# print(V1)

v1 = np.random.uniform(low=-10, high=10, size=(5, 3))
print(v1)
wait = [1,2,3]
v2 = v1.tolist()
v2.append(wait)
print(v2)
v1 = np.array(v2)
print(v1)

max_stars_repo_name: gmjustforfun/code

max_issues_repo_name: gmjustforfun/code

max_forks_repo_name: gmjustforfun/code


text:
import os, glob
import pandas as pd
import numpy as np
import re

abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)

inputLoc = ".\\raw\\"
outputLoc = ".\\"

all_files = glob.glob(os.path.join(inputLoc, "*.csv"))
all_df = []
for f in all_files:
    df = pd.read_csv(f, sep=',', dtype=str)
    fileDate = re.findall('\\d\\d-\\d\\d-\\d\\d\\d\\d', f)[0]
    df['Date'] = fileDate
    df['Anzahl'] = df['Anzahl'].str.replace("\+|\*|\.", "", regex=True)
    df['Differenz zum Vortag'] = df['Differenz zum Vortag'].str.replace("\+|\*|\.", "", regex=True)
    df['Fälle/100.000 Einw.'] = df['Fälle/100.000 Einw.'].str.replace("\+|\*|\.", "", regex=True)
    df['Todesfälle'] = df['Todesfälle'].str.replace("\+|\*|\.|", "", regex=True)
    all_df.append(df)

combined_df = pd.concat(all_df, ignore_index=True, sort=True)
outputFileName = outputLoc + 'RKI_Covid19_ALL.csv'
# Sort by date
combined_df["Date"] = pd.to_datetime(combined_df["Date"], format="%d-%m-%Y")
combined_df = combined_df.sort_values(by="Date")
combined_df = combined_df[['Date', 'Bundesland', 'Besonders betroffene Gebiete', 'Anzahl', 'Differenz zum Vortag', 'Fälle/100.000 Einw.', 'Todesfälle']]

# Write to file
combined_df.to_csv(outputFileName, sep=',', encoding='utf-8', index=False)


text:
# Copyright (c) Microsoft Corporation. All rights reserved.
# Licensed under the MIT License.

import os
import time
import onnx
import logging
import argparse
from pathlib import Path
from onnx import helper
import numpy as np
from typing import List
import torch
from transformers import GPT2Config
from gpt2_helper import PRETRAINED_GPT2_MODELS
from convert_to_onnx import main as convert_gpt2_to_onnx
from benchmark_helper import Precision
from onnx import onnx_pb as onnx_proto
"""
This converts GPT2 model to onnx with beam search operator.

Examples:
   python convert_beam_search.py -m gpt2 --gpt2_onnx .\onnx_models\gpt2_past_fp32.onnx --output .\onnx_models\gpt2_beam_search.onnx --output_sequences_scores
"""

config: GPT2Config = None

logger = logging.getLogger('')

def parse_arguments(argv=None):
    parser = argparse.ArgumentParser()

    parser.add_argument('-m',
                        '--model_name_or_path',
                        required=True,
                        type=str,
                        help='Model path, or pretrained model name in the list: ' + ', '.join(PRETRAINED_GPT2_MODELS))

    parser.add_argument('--cache_dir',
                        required=False,
                        type=str,
                        default=os.path.join('.', 'cache_models'),
                        help='Directory to cache pre-trained models')

    parser.add_argument('--gpt2_onnx',
                        required=True,
                        type=str,
                        help='Output directory for GPT-2 onnx model, or model path ends with .onnx')

    parser.add_argument('--output',
                        required=False,
                        type=str,
                        help='Output directory for beam search model, or model path ends with .onnx')

    parser.add_argument("-p",
                        "--precision",
                        required=False,
                        type=Precision,
                        default=Precision.FLOAT32,
                        choices=[Precision.FLOAT32, Precision.FLOAT16],
                        help="Precision of model to run. fp32 for full precision, fp16 for half or mixed precision")

    parser.add_argument('--use_gpu', required=False, action='store_true', help="use GPU for inference")
    parser.set_defaults(use_gpu=False)

    parser.add_argument('-e', '--use_external_data_format', required=False, action='store_true')
    parser.set_defaults(use_external_data_format=False)

    parser.add_argument('--disable_parity', required=False, action='store_true', help="do not run parity test")
    parser.set_defaults(disable_parity=False)

    parser.add_argument('--torch_performance', required=False, action='store_true', help="test PyTorch performance")
    parser.set_defaults(torch_performance=False)

    parser.add_argument('--total_runs',
                        required=False,
                        type=int,
                        default=1,
                        help='Number of times of inference for latency measurement')

    beam_search_group = parser.add_argument_group("beam search options")

    beam_search_group.add_argument('--output_sequences_scores',
                                   required=False,
                                   action='store_true',
                                   help="output sequences scores")
    beam_search_group.set_defaults(output_sequences_scores=False)

    beam_search_group.add_argument('--output_token_scores',
                                   required=False,
                                   action='store_true',
                                   help="output token scores")
    beam_search_group.set_defaults(output_token_scores=False)

    beam_search_group.add_argument('--early_stopping', required=False, action='store_true')
    beam_search_group.set_defaults(early_stopping=False)

    beam_search_group.add_argument('--min_length', type=int, required=False, default=1, help='Min sequence length')

    beam_search_group.add_argument('--max_length', type=int, required=False, default=50, help='Max sequence length')

    beam_search_group.add_argument('--no_repeat_ngram_size',
                                   type=int,
                                   required=False,
                                   default=0,
                                   help='No repeat ngram size')

    beam_search_group.add_argument('--num_beams', type=int, required=False, default=4, help='Beam size')

    beam_search_group.add_argument('--num_return_sequences',
                                   type=int,
                                   required=False,
                                   default=1,
                                   help='Number of return sequence <= num_beams')

    beam_search_group.add_argument('--temperature',
                                   type=float,
                                   required=False,
                                   default=1,
                                   help='Softmax temperature for output logits.')

    beam_search_group.add_argument('--length_penalty',
                                   type=float,
                                   required=False,
                                   default=1,
                                   help='Positive. >1 to penalize and <1 to encorage short sentence.')

    beam_search_group.add_argument('--repetition_penalty',
                                   type=float,
                                   required=False,
                                   default=1,
                                   help='Positive. >1 to penalize and <1 to encorage.')

    beam_search_group.add_argument('--vocab_size',
                                   type=int,
                                   required=False,
                                   default=-1,
                                   help="Vocab_size of the underlying model")

    beam_search_group.add_argument(
        '--prefix_vocab_mask',
        required=False,
        action='store_true',
        help="This vocab mask applies only to first iteration, enable if last word in query might need auto complete")
    beam_search_group.set_defaults(prefix_vocab_mask=False)

    args = parser.parse_args(argv)

    return args

def gpt2_to_onnx(args):
    model_name = args.model_name_or_path

    print(f"use convert_to_onnx.py to convert model {model_name} to onnx {args.gpt2_onnx} ...")
    arguments = [
        '--model_name_or_path',
        model_name,
        '--output',
        args.gpt2_onnx,
        '--optimize_onnx',
        '--precision',
        'fp32' if args.precision == Precision.FLOAT32 else 'fp16',
        '--test_runs',
        '1',
        '--test_cases',
        '10',
        '--use_int32_inputs'  # BeamSearch requires to use int32 for input_ids, postion_ids and attention_mask
    ]
    if args.use_gpu:
        arguments.append('--use_gpu')
    if args.use_external_data_format:
        arguments.append('--use_external_data_format')

    if args.precision == Precision.FLOAT16:
        assert args.use_gpu, "fp16 or mixed precision model cannot run in CPU. Please add --use_gpu"
        # TODO: Use auto mixed precision for fp16 conversion: arguments.append('--auto_mixed_precision')
        #       Need change cuda kernel to support a combination of fp32 logits and fp16 past state.
        #       Currently logits and past state shall be same data type.
        arguments.extend(['--op_block_list', 'Add', 'LayerNormalization', 'FastGelu'])
    convert_gpt2_to_onnx(arguments)

def shape_inference(gpt2_onnx_path):
    # Run symbolic shape inference to walk around ORT shape inference issue for subgraph.
    from onnxruntime.tools.symbolic_shape_infer import SymbolicShapeInference
    out = SymbolicShapeInference.infer_shapes(onnx.load(gpt2_onnx_path), auto_merge=True, guess_output_rank=False)
    if out:
        # TODO: Use external format if input has extra data.
        onnx.save(out, gpt2_onnx_path)
    else:
        print("Failed to run symbolic shape inference on the model.")

def create_ort_session(model_path, use_gpu):
    from onnxruntime import SessionOptions, InferenceSession, __version__ as ort_version, GraphOptimizationLevel, get_available_providers
    sess_options = SessionOptions()
    sess_options.graph_optimization_level = GraphOptimizationLevel.ORT_DISABLE_ALL
    execution_providers = ['CUDAExecutionProvider', 'CPUExecutionProvider'] if use_gpu else ['CPUExecutionProvider']
    if use_gpu:
        if 'CUDAExecutionProvider' not in get_available_providers():
            raise RuntimeError("CUDAExecutionProvider is not avaiable for --use_gpu!")
        else:
            print("use CUDAExecutionProvider")

    ort_session = InferenceSession(model_path, sess_options, providers=execution_providers)
    return ort_session

def verify_gpt2_subgraph(graph, precision):
    is_float16 = (Precision.FLOAT16 == precision)

    input_count = len(graph.input)
    layer_count = input_count - 3

    expected_inputs = ['input_ids', 'position_ids', 'attention_mask'] + [f"past_{i}" for i in range(layer_count)]
    if len(graph.input) != len(expected_inputs):
        raise ValueError(f"Number of inputs expected to be {len(expected_inputs)}. Got {len(graph.input)}")

    for i, expected_input in enumerate(expected_inputs):
        if graph.input[i].name != expected_input:
            raise ValueError(f"Input {i} is expected to be {expected_input}. Got {graph.input[i].name}")

        expected_type = onnx_proto.TensorProto.INT32
        if i >= 3:
            expected_type = onnx_proto.TensorProto.FLOAT16 if is_float16 else onnx_proto.TensorProto.FLOAT32

        if graph.input[i].type.tensor_type.elem_type != expected_type:
            raise ValueError(
                f"Input {i} is expected to have onnx data type {expected_type}. Got {graph.input[i].type.tensor_type.elem_type}"
            )
    print("Verifying GPT-2 graph inputs: name and data type are good.")

    expected_outputs = ['logits'] + [f"present_{i}" for i in range(layer_count)]
    if len(graph.output) != len(expected_outputs):
        raise ValueError(f"Number of outputs expected to be {len(expected_outputs)}. Got {len(graph.output)}")

    for i, expected_output in enumerate(expected_outputs):
        if graph.output[i].name != expected_output:
            raise ValueError(f"Output {i} is expected to be {expected_output}. Got {graph.output[i].name}")

        expected_type = onnx_proto.TensorProto.FLOAT16 if is_float16 else onnx_proto.TensorProto.FLOAT32
        if graph.output[i].type.tensor_type.elem_type != expected_type:
            raise ValueError(
                f"Input {i} is expected to have onnx data type {expected_type}. Got {graph.output[i].type.tensor_type.elem_type}"
            )
    print("Verifying GPT-2 graph outputs: name and data type are good.")

    # TODO: verify shapes of inputs and outputs.
    return

def convert_model(args):
    if os.path.exists(args.gpt2_onnx):
        print(f"skip convert_to_onnx since path existed: {args.gpt2_onnx}")
    else:
        gpt2_to_onnx(args)

    print(f"Run symbolic shape inference on {args.gpt2_onnx}. The file will be overwritten.")
    shape_inference(args.gpt2_onnx)

    global config
    config = GPT2Config.from_pretrained(args.model_name_or_path, cache_dir=args.cache_dir)
    print(config)

    eos_token_id = config.eos_token_id
    pad_token_id = config.eos_token_id
    vocab_size = config.vocab_size

    # if vocab_size is given in parameters use that.
    if args.vocab_size != -1:
        vocab_size = args.vocab_size

    model = onnx.load(args.gpt2_onnx)
    verify_gpt2_subgraph(model.graph, args.precision)

    model.graph.name = "gpt2 subgraph"
    inputs = [
        "input_ids", "max_length", "min_length", "num_beams", "num_return_sequences", "temperature", "length_penalty",
        "repetition_penalty", "vocab_mask"
    ]
    if args.prefix_vocab_mask:
        inputs.append("prefix_vocab_mask")

    outputs = ["sequences"]
    if args.output_sequences_scores:
        outputs.append("sequences_scores")

    if args.output_token_scores:
        assert args.output_sequences_scores, "--output_token_scores requires --output_sequences_scores"
        outputs.append("scores")

    node = helper.make_node('BeamSearch', inputs=inputs, outputs=outputs, name='BeamSearch_GPT2')
    node.domain = "com.microsoft"
    node.attribute.extend([
        helper.make_attribute("eos_token_id", eos_token_id),
        helper.make_attribute("pad_token_id", pad_token_id),
        helper.make_attribute("no_repeat_ngram_size", args.no_repeat_ngram_size),
        helper.make_attribute("early_stopping", 1 if args.early_stopping else 0),
        helper.make_attribute("body", model.graph),
    ])

    from onnx import TensorProto

    # graph inputs
    input_ids = helper.make_tensor_value_info('input_ids', TensorProto.INT32, ['batch_size', 'sequence_length'])
    max_length = helper.make_tensor_value_info('max_length', TensorProto.INT32, [1])
    min_length = helper.make_tensor_value_info('min_length', TensorProto.INT32, [1])
    num_beams = helper.make_tensor_value_info('num_beams', TensorProto.INT32, [1])
    num_return_sequences = helper.make_tensor_value_info('num_return_sequences', TensorProto.INT32, [1])
    temperature = helper.make_tensor_value_info('temperature', TensorProto.FLOAT, [1])
    length_penalty = helper.make_tensor_value_info('length_penalty', TensorProto.FLOAT, [1])
    repetition_penalty = helper.make_tensor_value_info('repetition_penalty', TensorProto.FLOAT, [1])
    vocab_mask = helper.make_tensor_value_info('vocab_mask', TensorProto.INT32, [vocab_size])

    graph_inputs = [
        input_ids, max_length, min_length, num_beams, num_return_sequences, temperature, length_penalty,
        repetition_penalty, vocab_mask
    ]

    if args.prefix_vocab_mask:
        prefix_vocab_mask = helper.make_tensor_value_info('prefix_vocab_mask', TensorProto.INT32,
                                                          ['batch_size', vocab_size])
        graph_inputs.append(prefix_vocab_mask)

    # graph outputs
    sequences = helper.make_tensor_value_info('sequences', TensorProto.INT32,
                                              ['batch_size', 'num_return_sequences', 'max_length'])

    sequences_scores = helper.make_tensor_value_info('sequences_scores', TensorProto.FLOAT,
                                                     ['batch_size', 'num_return_sequences'])

    scores = helper.make_tensor_value_info('scores', TensorProto.FLOAT,
                                           ['max_length - sequence_length', 'batch_size', 'num_beams', vocab_size])

    initializers = []

    graph_outputs = [sequences]

    if args.output_sequences_scores:
        graph_outputs.append(sequences_scores)

    if args.output_token_scores:
        graph_outputs.append(scores)

    new_graph = helper.make_graph([node], 'gpt2-beam-search', graph_inputs, graph_outputs, initializers)

    # Create the model
    new_model = helper.make_model(new_graph, producer_name='onnxruntime.transformers', opset_imports=model.opset_import)
    onnx.save(new_model, args.output)

def test_torch_performance(args, model, input_ids, attention_mask, eos_token_id, pad_token_id, bad_words_ids):
    if args.use_gpu and not torch.cuda.is_available():
        logger.error("Please install PyTorch with Cuda, and use a machine with GPU for testing gpu performance.")
        return None

    if args.precision == Precision.FLOAT16:
        model.half()

    device = torch.device("cuda:0" if args.use_gpu else "cpu")
    model.to(device)

    torch.set_grad_enabled(False)
    input_ids = input_ids.to(device)
    attention_mask = attention_mask.to(device)

    torch_latency = []
    for _ in range(args.total_runs):
        start = time.time()
        _ = model.generate(input_ids=input_ids,
                           attention_mask=attention_mask,
                           max_length=args.max_length,
                           min_length=args.min_length,
                           num_beams=args.num_beams,
                           early_stopping=args.early_stopping,
                           no_repeat_ngram_size=args.no_repeat_ngram_size,
                           eos_token_id=eos_token_id,
                           pad_token_id=pad_token_id,
                           num_return_sequences=args.num_return_sequences,
                           temperature=args.temperature,
                           length_penalty=args.length_penalty,
                           repetition_penalty=args.repetition_penalty,
                           bad_words_ids=bad_words_ids,
                           return_dict_in_generate=True,
                           output_scores=args.output_sequences_scores or args.output_token_scores)
        torch_latency.append(time.time() - start)
    batch_size = input_ids.shape[0]
    from benchmark_helper import get_latency_result
    return get_latency_result(torch_latency, batch_size)

def test_model(args, use_vocab_mask: bool = False, sentences: List[str] = None):

    if args.prefix_vocab_mask:
        print("Skipping parity test as prefix vocab mask is not implemented by Hugging Face")
        return

    from transformers import GPT2Tokenizer, GPT2LMHeadModel

    tokenizer = GPT2Tokenizer.from_pretrained(args.model_name_or_path, cache_dir=args.cache_dir)
    tokenizer.padding_side = "left"
    tokenizer.pad_token = tokenizer.eos_token

    model = GPT2LMHeadModel.from_pretrained(args.model_name_or_path,
                                            cache_dir=args.cache_dir,
                                            pad_token_id=tokenizer.eos_token_id)

    # Use different length sentences to test batching
    if sentences is None:
        sentences = ["The product is released", "I enjoy walking in the park", "Test best way to invest"]

    inputs = tokenizer(sentences, return_tensors='pt', padding=True)
    input_ids = inputs["input_ids"]
    attention_mask = inputs["attention_mask"]

    bad_words = "walk in park"
    bad_words_ids = tokenizer.encode(bad_words, add_prefix_space=True)
    bad_words_ids = [[word_id] for word_id in bad_words_ids]  # Convert to list of list
    if use_vocab_mask:
        print("bad_words_ids", bad_words_ids)
    else:
        bad_words_ids = None

    global config
    config = model.config
    eos_token_id = config.eos_token_id
    pad_token_id = config.eos_token_id
    vocab_size = config.vocab_size

    torch_decoded_sequences = []
    if not args.disable_parity:
        print('-' * 50)
        print("Test PyTorch model and beam search with huggingface transformers...")
        beam_outputs = model.generate(input_ids=input_ids,
                                      attention_mask=attention_mask,
                                      max_length=args.max_length,
                                      min_length=args.min_length,
                                      num_beams=args.num_beams,
                                      early_stopping=args.early_stopping,
                                      no_repeat_ngram_size=args.no_repeat_ngram_size,
                                      eos_token_id=eos_token_id,
                                      pad_token_id=pad_token_id,
                                      num_return_sequences=args.num_return_sequences,
                                      temperature=args.temperature,
                                      length_penalty=args.length_penalty,
                                      repetition_penalty=args.repetition_penalty,
                                      bad_words_ids=bad_words_ids,
                                      return_dict_in_generate=True,
                                      output_scores=args.output_sequences_scores or args.output_token_scores)
        print("input_ids", input_ids)
        print("huggingface transformers outputs:")
        print("sequences", beam_outputs.sequences)
        if args.output_sequences_scores:
            print("sequences_scores", beam_outputs.sequences_scores)
        if args.output_token_scores:
            print("scores", beam_outputs.scores)
        for i, sequence in enumerate(beam_outputs.sequences):
            decoded_sequence = tokenizer.decode(sequence, skip_special_tokens=True)
            torch_decoded_sequences.append(decoded_sequence)
            print("{}: {}".format(i, decoded_sequence))

    print('-' * 50)
    print("Test ONNX model and bream search with onnxruntime...")

    ort_session = create_ort_session(args.output, args.use_gpu)

    vocab_mask = np.ones((vocab_size), dtype=np.int32)
    if use_vocab_mask:
        for bad_word_id in bad_words_ids:
            vocab_mask[bad_word_id] = 0

    inputs = {
        "input_ids": input_ids.cpu().numpy().astype(np.int32),
        "max_length": np.array([args.max_length], dtype=np.int32),
        "min_length": np.array([args.min_length], dtype=np.int32),
        "num_beams": np.array([args.num_beams], dtype=np.int32),
        "num_return_sequences": np.array([args.num_return_sequences], dtype=np.int32),
        "temperature": np.array([args.temperature], dtype=np.float32),
        "length_penalty": np.array([args.length_penalty], dtype=np.float32),
        "repetition_penalty": np.array([args.repetition_penalty], dtype=np.float32),
        "vocab_mask": vocab_mask
    }

    test_data_dir = Path(args.output).parent.as_posix()
    print("test_data_dir", test_data_dir)
    from bert_test_data import output_test_data
    all_inputs = [inputs]
    for i, inputs in enumerate(all_inputs):
        dir = os.path.join(test_data_dir, 'test_data_set_' + str(i))
        output_test_data(dir, inputs)

    print("inputs", inputs)

    # Test performance
    latency = []
    for _ in range(args.total_runs):
        start = time.time()
        result = ort_session.run(None, inputs)
        latency.append(time.time() - start)
    batch_size = input_ids.shape[0]
    from benchmark_helper import get_latency_result
    output = get_latency_result(latency, batch_size)

    print("ORT outputs:")
    sequences = result[0]
    print("sequences", sequences)
    if args.output_sequences_scores:
        print("sequences_scores", result[1])
    if args.output_token_scores:
        print("scores", result[2])

    (batch_size, num_sequences, max_length) = sequences.shape
    ort_decoded_sequences = []
    for i in range(batch_size):
        for j in range(num_sequences):
            decoded_sequence = tokenizer.decode(sequences[i][j], skip_special_tokens=True)
            ort_decoded_sequences.append(decoded_sequence)
            print(f"batch {i} sequence {j}: {decoded_sequence}")

    if not args.disable_parity:
        torch_sequences = beam_outputs.sequences.reshape(batch_size, args.num_return_sequences, -1)
        ort_sequences = torch.LongTensor(sequences)
        print("-" * 50)
        print("Torch Sequences:")
        print(torch_sequences)
        print(torch_decoded_sequences)
        print("-" * 50)
        print("ORT Sequences:")
        print(ort_sequences)
        print(ort_decoded_sequences)
        print("-" * 50)
        # Compare the generated text instead of word IDs since ORT pads to max sequence length but Torch not.
        is_same = (torch_decoded_sequences == ort_decoded_sequences)
        print("Torch and ORT result is ", "same" if is_same else "different")
        output["parity"] = is_same

    if args.torch_performance:
        torch_latency_output = test_torch_performance(args, model, input_ids, attention_mask, eos_token_id,
                                                      pad_token_id, bad_words_ids)
        print("Torch Latency", torch_latency_output)

    print("ORT", output)
    return output

def main(argv=None, sentences=None):
    args = parse_arguments(argv)

    if os.path.exists(args.output):
        print(f"skip conversion since path existed: {args.output}")
    else:
        convert_model(args)

    return test_model(args, use_vocab_mask=True, sentences=sentences)

if __name__ == '__main__':
    main()

max_stars_repo_name: jambayk/onnxruntime

max_issues_repo_name: jambayk/onnxruntime

max_forks_repo_name: jambayk/onnxruntime


text:
import numpy as np
from matplotlib import pyplot as plt
import pandas as pd
from pathlib import Path

wd = f'{Path(__file__).resolve().parent}'

for obj in ['flood', 'pumping']:
  for gcms in ['allgcms', 'real5', 'priority4']:

    df = pd.read_csv(f'{wd}/sobol_{gcms}_si_{obj}.csv', index_col=0, parse_dates=True)
    #df = pd.read_csv('sobol_%s_si_%s.csv' % (gcms,obj), index_col=0, parse_dates=True)
    cond = df.columns.str.contains('_T')

    for s in ['si', 'st']:
      df_s = df.loc[:, ~cond] if s == 'si' else df.loc[:, cond]

      for a in ['raw', 'avg']:

        if a == 'raw':
          df_s.plot.area()
        else:
          df_s.rolling(10).mean().plot.area()
          # plt.xlim(['2030','2100'])

        plt.savefig(f'{wd}/figs/sobol/svg/{obj}_{gcms}_{s}_{a}.svg')
        plt.savefig(f'{wd}/figs/sobol/png/{obj}_{gcms}_{s}_{a}.png', dpi=300)
        # plt.savefig('figs/sobol/svg/%s_%s_%s_%s.svg' % (obj,gcms,s,a))
        # plt.savefig('figs/sobol/png/%s_%s_%s_%s.png' % (obj,gcms,s,a), dpi=300)
        plt.close()

max_stars_repo_name: JohnRushKucharski/systemuncertainty

max_issues_repo_name: JohnRushKucharski/systemuncertainty

max_forks_repo_name: JohnRushKucharski/systemuncertainty


text:
# Copyright 2021 The TensorFlow Authors. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ==============================================================================
"""Functional tests for fused batch-norm related to determinism."""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np

from tensorflow.python.eager import backprop
from tensorflow.python.framework import config
from tensorflow.python.framework import constant_op
from tensorflow.python.framework import dtypes
from tensorflow.python.framework import errors_impl
from tensorflow.python.framework import test_util
from tensorflow.python.ops import nn_impl
# The following import is required to register the gradient function
from tensorflow.python.ops.nn_grad import _FusedBatchNormV3Grad  # pylint: disable=unused-import
from tensorflow.python.platform import test

@test_util.run_all_in_graph_and_eager_modes
class FusedBatchNormalizationDeterministicTest(test.TestCase):
  """Test determinsitic functionality and exceptions for FusedBatchNorm.

  Test that tf.errors.UnimplementedError is thrown, as
  appropriate, by the GPU code-path through FusedBatchNormFreezeGrad when
  deterministic ops are enabled. This test assumes that
  nn_fused_batchnorm_test.py runs equivalent test cases when deterministic ops
  are not enabled and will therefore detect erroneous exception throwing in
  those cases.

  Also test that the other code-paths, running on both CPU and GPU, operate
  deterministically.
  """

  def _genParams(self, data_format, x_dtype, large_batch):
    if large_batch:
      batch_size = 5000
      height = width = 4
    else:
      batch_size = 10
      height = 5
      width = 5000
    channel_count = 3
    if data_format == 'NHWC':
      x_shape = (batch_size, height, width, channel_count)
    else:  # 'NCHW'
      x_shape = (batch_size, channel_count, height, width)
    # Using random_ops.random_normal would produce different values on each run
    x = constant_op.constant(np.random.normal(size=x_shape), dtype=x_dtype)
    scale_shape = (channel_count,)
    scale = constant_op.constant(
        np.random.normal(size=scale_shape), dtype=dtypes.float32)
    offset = constant_op.constant(
        np.random.normal(size=scale_shape), dtype=dtypes.float32)
    mean = np.random.normal(size=scale_shape)
    variance = np.random.normal(size=scale_shape)
    y_shape = x_shape
    y_dtype = x_dtype
    upstream_gradients = constant_op.constant(
        np.random.normal(size=y_shape), dtype=y_dtype)
    return x, scale, offset, mean, variance, upstream_gradients

  def testForward(self):
    with self.cached_session():
      for data_format in ['NHWC', 'NCHW']:
        for large_batch in [False, True]:
          for x_dtype in [dtypes.float16, dtypes.float32]:  # skipping bfloat16
            x, scale, offset, mean, variance, _ = self._genParams(
                data_format, x_dtype, large_batch)
            for is_training in [False, True]:
              op_output = nn_impl.fused_batch_norm(
                  x,
                  scale,
                  offset,
                  mean,
                  variance,
                  data_format=data_format,
                  is_training=is_training,
                  exponential_avg_factor=1.01)
              y_a, running_mean_a, running_var_a = op_output
              y_a = self.evaluate(y_a)
              if is_training:
                running_mean_a = self.evaluate(running_mean_a)
                running_var_a = self.evaluate(running_var_a)
              for _ in range(5):
                op_output_b = nn_impl.fused_batch_norm(
                    x,
                    scale,
                    offset,
                    mean,
                    variance,
                    data_format=data_format,
                    is_training=is_training,
                    exponential_avg_factor=1.01)
                y_b, running_mean_b, running_var_b = op_output_b
                y_b = self.evaluate(y_b)
                self.assertAllEqual(y_a, y_b)
                if is_training:
                  running_mean_b = self.evaluate(running_mean_b)
                  running_var_b = self.evaluate(running_var_b)
                  self.assertAllEqual(running_mean_a, running_mean_b)
                  self.assertAllEqual(running_var_a, running_var_b)

  @test_util.disable_xla('XLA is deterministic')
  def testBackward(self):
    with self.cached_session():
      for data_format in ['NHWC', 'NCHW']:
        for large_batch in [False, True]:
          for x_dtype in [dtypes.float16, dtypes.float32]:  # skipping bfloat16
            params = self._genParams(data_format, x_dtype, large_batch)
            x, scale, offset, mean, variance, upstream_gradients = params
            for is_training in [False, True]:
              for backprop_to in [x, scale, offset]:
                with backprop.GradientTape(persistent=True) as tape:
                  tape.watch(backprop_to)
                  op_output = nn_impl.fused_batch_norm(
                      x,
                      scale,
                      offset,
                      mean,
                      variance,
                      data_format=data_format,
                      is_training=is_training,
                      exponential_avg_factor=0.99)
                  gradient_injector_output = op_output[0] * upstream_gradients
                if (len(config.list_physical_devices('GPU')) and
                    not is_training):
                  # Only backprop to offset is nondeterministic (on GPU, when
                  # is_training=False), but backprop to the other parameters is
                  # calculated using the same kernel.
                  with self.assertRaisesRegex(
                      errors_impl.UnimplementedError,
                      'A deterministic GPU implementation of fused batch-norm' +
                      ' backprop, when training is disabled, is not currently' +
                      ' available.'):
                    grad = tape.gradient(gradient_injector_output, backprop_to)
                    self.evaluate(grad)
                else:
                  grad_a = tape.gradient(gradient_injector_output, backprop_to)
                  grad_a = self.evaluate(grad_a)
                  for _ in range(5):
                    grad_b = tape.gradient(gradient_injector_output,
                                           backprop_to)
                    grad_b = self.evaluate(grad_b)
                    self.assertAllEqual(grad_a, grad_b)

if __name__ == '__main__':
  # TODO(reedwm): Merge this file with nn_fused_batchnorm_test.py
  config.enable_deterministic_ops(True)
  test.main()

max_stars_repo_name: kalosisz/tensorflow

max_issues_repo_name: kalosisz/tensorflow

max_forks_repo_name: kalosisz/tensorflow


text:
import numpy as np
from scipy.spatial import cKDTree
from tqdm import tqdm
import btrack
from btrack.constants import BayesianUpdates
import pandas as pd
import logging

logger = logging.getLogger()

tracker_config = {
    "TrackerConfig":
        {
            "MotionModel":
                {
                    "name": "cell_motion",
                    "dt": 1.0,
                    "measurements": 3,
                    "states": 6,
                    "accuracy": 7.5,
                    "prob_not_assign": 0.001,
                    "max_lost": 3,
                    "A": {
                        "matrix": [1, 0, 0, 1, 0, 0,
                                   0, 1, 0, 0, 1, 0,
                                   0, 0, 1, 0, 0, 1,
                                   0, 0, 0, 1, 0, 0,
                                   0, 0, 0, 0, 1, 0,
                                   0, 0, 0, 0, 0, 1]
                    },
                    "H": {
                        "matrix": [1, 0, 0, 0, 0, 0,
                                   0, 1, 0, 0, 0, 0,
                                   0, 0, 1, 0, 0, 0]
                    },
                    "P": {
                        "sigma": 150.0,
                        "matrix": [0.1, 0, 0, 0, 0, 0,
                                   0, 0.1, 0, 0, 0, 0,
                                   0, 0, 0.1, 0, 0, 0,
                                   0, 0, 0, 1, 0, 0,
                                   0, 0, 0, 0, 1, 0,
                                   0, 0, 0, 0, 0, 1]
                    },
                    "G": {
                        "sigma": 15.0,
                        "matrix": [0.5, 0.5, 0.5, 1, 1, 1]

                    },
                    "R": {
                        "sigma": 5.0,
                        "matrix": [1, 0, 0,
                                   0, 1, 0,
                                   0, 0, 1]
                    }
                },
            "ObjectModel":
                {},
            "HypothesisModel":
                {
                    "name": "cell_hypothesis",
                    "hypotheses": ["P_FP", "P_init", "P_term", "P_link", "P_branch", "P_dead"],
                    "lambda_time": 5.0,
                    "lambda_dist": 3.0,
                    "lambda_link": 10.0,
                    "lambda_branch": 50.0,
                    "eta": 1e-10,
                    "theta_dist": 20.0,
                    "theta_time": 5.0,
                    "dist_thresh": 40,
                    "time_thresh": 2,
                    "apop_thresh": 5,
                    "segmentation_miss_rate": 0.1,
                    "apoptosis_rate": 0.001,
                    "relax": True
                }
        }
}

def track(adata,
          time_var="Metadata_Time",
          group_vars=("Metadata_Well", "Metadata_Field"),
          x_loc="Cells_Location_Center_X",
          y_loc="Cells_Location_Center_Y",
          parent_id="Metadata_Trace_Parent",
          tree_id="Metadata_Trace_Tree",
          do_filter=True,
          allowed_dummies=0,
          min_track_len=10,
          max_search_radius=100,
          field_size=(2048, 2048),
          approx=False
          ):
    """Wrapper for Bayesian Tracker.

    Tracks objects per field of view defined by group_vars.
    Objects can be filtered for trajectory length and dummy objects.

    Reference:
        Automated deep lineage tree analysis using a Bayesian single cell tracking approach
        Ulicna K, Vallardi G, Charras G and Lowe AR.
        bioRxiv (2020)

    Args:
        adata (anndata.AnnData): Morphological data from cells with different time stamps.
        time_var (str): Variable name for time point in annotations.
        group_vars (iterable): Variables in annotations.
            Should point to specific wells/ or fields on plates that can be compared over time.
        x_loc (str): Identifier for x location in annotations.
        y_loc (str): Identifier for y location in annotations.
        parent_id (str): Variable name used to store index of parent cell.
        tree_id (str): Variable name used to store unique branch number for a certain field.
        do_filter (bool): Filter tracks with dummies and minimum length if True.
        allowed_dummies (int): Max allowed dummies in track.
        min_track_len (int): Min track length.
        max_search_radius (int): Local spatial search radius (pixels)
        field_size (tuple of ints): height and width in of field.
        approx (bool): Speed up processing on very large datasets.

    Returns:
        adata (anndata.AnnData)
    """
    # check variables
    if isinstance(group_vars, str):
        group_vars = [group_vars]
    elif isinstance(group_vars, tuple):
        group_vars = list(group_vars)

    if isinstance(group_vars, list):
        assert all(gv in adata.obs.columns for gv in group_vars), \
            f"One or all group_vars not in .obs.columns: {group_vars}"
    else:
        raise KeyError(f"Expected type(list) or type(str) for group_vars, "
                       f"instead got {type(group_vars)}")

    assert time_var in adata.obs.columns, f"time_var not in .obs.columns: {time_var}"
    assert x_loc in adata.obs.columns, f"x_loc not in .obs.columns: {x_loc}"
    assert x_loc in adata.obs.columns, f"x_loc not in .obs.columns: {x_loc}"
    assert y_loc in adata.obs.columns, f"y_loc not in .obs.columns: {y_loc}"
    assert x_loc != y_loc, f"x_loc expected to be different from y_loc, " \
                           f"x_loc: {x_loc}, y_loc: {y_loc}"

    field_height, field_width = field_size

    # create new column to store index of parent object
    adata.obs[parent_id] = np.nan
    # create new column and store id for every trace tree
    adata.obs[tree_id] = np.nan

    # load config settings
    config = tracker_config['TrackerConfig']
    t_config = {
        "MotionModel": btrack.utils.read_motion_model(config),
        "ObjectModel": btrack.utils.read_object_model(config),
        "HypothesisModel": btrack.optimise.hypothesis.read_hypothesis_model(config),
    }

    tree_start_id = 0

    # iterate over every field and get field at different times
    total_its = np.prod([len(adata.obs[gv].unique()) for gv in group_vars])
    for ix, (groups, field_df) in tqdm(enumerate(adata.obs.groupby(list(group_vars))),
                                       desc='Tracking cells in every field',
                                       total=int(total_its)):
        # create objects for btracker
        z = np.zeros((len(field_df),))
        _id = field_df.index

        objects_raw = {'x': field_df[x_loc].to_numpy(),
                       'y': field_df[y_loc].to_numpy(),
                       'z': z,
                       't': field_df[time_var].to_numpy(),
                       'id': _id}
        objects = btrack.dataio.objects_from_dict(objects_raw)

        # start tracker
        logger.disabled = True
        with btrack.BayesianTracker(verbose=False) as tracker:
            tracker.configure(t_config)

            tracker.max_search_radius = max_search_radius

            if approx:
                tracker.update_method = BayesianUpdates.APPROXIMATE

            tracker.append(objects)

            # set the volume (Z axis volume is set very large for 2D data)
            tracker.volume = ((0, field_height), (0, field_width), (-1e5, 1e5))

            # track them
            tracker.track()

            # generate hypotheses and run the global optimizer
            tracker.optimize()

            tracks = tracker.tracks

        logger.disabled = False

        if do_filter:
            filtered_tracks = filter_tracks(tracks, track_len=min_track_len, allowed_dummies=allowed_dummies)
        else:
            filtered_tracks = tracks

        t_header = ["ID", "t"] + ["z", "y", "x"][-2:]
        p_header = ["t", "state", "generation", "root", "parent"]

        # ensure lexicographic ordering of tracks
        ordered = sorted(list(filtered_tracks), key=lambda t: t.ID)
        header = t_header + p_header
        dict_tracks = tracks_as_dict(ordered, header)

        # add tree and parent information
        dict_tracks = add_parent_id(dict_tracks)
        dict_tracks = add_tree_id(dict_tracks, start_id=tree_start_id)
        cat_tracks = cat_track_dicts(dict_tracks)
        output = pd.DataFrame(cat_tracks)

        tree_start_id += len(dict_tracks)

        # update adata
        if len(output) > 0:
            adata.obs.loc[output['id'], [tree_id, parent_id]] = output[['tree_id', 'parent_id']].to_numpy()

    return adata

def add_parent_id(tracks):
    new_tracks = []
    for tr in tracks:
        sort_ix = np.argsort(tr['t'])
        parent_id = np.zeros(sort_ix.shape)
        parent_id[0] = np.nan
        parent_id[1:] = tr['id'][sort_ix[:-1]]
        tr['parent_id'] = parent_id
        new_tracks.append(tr)

    return new_tracks

def add_tree_id(tracks, start_id=None):
    n_tracks = len(tracks)
    if start_id is None:
        start_id = 0

    any_key = list(tracks[0].keys())[0] if len(tracks) > 0 else None
    ids = list(range(start_id, start_id + n_tracks))
    new_tracks = []
    for _id, tr in zip(ids, tracks):
        tr['tree_id'] = [_id] * len(tr[any_key])
        new_tracks.append(tr)

    return new_tracks

def tracks_as_dict(tracks: list, properties: list):
    dict_tracks = []

    for ix, tr in enumerate(tracks):
        trk = tr.to_dict(properties)

        data = {}

        for key in trk.keys():
            prop = trk[key]
            if not isinstance(prop, (list, np.ndarray)):
                prop = [prop] * len(tr)

            assert len(prop) == len(tr)
            data[key] = prop

        dict_tracks.append(data)

    return dict_tracks

def cat_track_dicts(tracks: list):
    data = {}

    for ix, trk in enumerate(tracks):

        if not data:
            data = {k: [] for k in trk.keys()}

        for key in data.keys():
            prop = trk[key]
            assert (
                isinstance(prop, (list, np.ndarray))
            ), "all track properties must be lists"

            data[key].append(prop)

    for key in data.keys():
        data[key] = np.concatenate(data[key])

    return data

def filter_tracks(tracks, allowed_dummies=0, track_len=10):
    dummy_filtered = []
    if allowed_dummies is not None:
        for tr in tracks:
            if np.sum(tr['dummy']) <= allowed_dummies:
                dummy_filtered.append(tr)
    else:
        dummy_filtered = tracks

    len_filtered = []
    if track_len is not None:
        for tr in dummy_filtered:
            if len(tr) >= track_len:
                len_filtered.append(tr)
    else:
        len_filtered = dummy_filtered

    return len_filtered

def track_nn(adata,
             time_var="Metadata_Time",
             group_vars=("Metadata_Well", "Metadata_Field"),
             x_loc="Cells_Location_Center_X",
             y_loc="Cells_Location_Center_Y",
             parent_id="Metadata_Trace_Parent",
             tree_id="Metadata_Trace_Tree",
             start_tp=0):
    """Neares neighbor tracking.

    Iterates over every well or field that can be traced over time.
    Fore every object and time point calculate the nearest neighbor from
    the previous time point.

    Args:
        adata (anndata.AnnData): Morphological data from cells with different time stamps.
        time_var (str): Variable name for time point in annotations.
        group_vars (iterable): Variables in annotations.
            Should point to specific wells/ or fields on plates that can be compared over time.
        x_loc (str): Identifier for x location in annotations.
        y_loc (str): Identifier for y location in annotations.
        parent_id (str): Variable name used to store index of parent cell.
        tree_id (str): Variable name used to store unique branch number for a certain field.
        start_tp (int): Start time point.

    Returns:
        adata (anndata.AnnData)
    """

    # check variables
    if isinstance(group_vars, str):
        group_vars = [group_vars]
    elif isinstance(group_vars, tuple):
        group_vars = list(group_vars)

    if isinstance(group_vars, list):
        assert all(gv in adata.obs.columns for gv in group_vars), \
            f"One or all group_vars not in .obs.columns: {group_vars}"
    else:
        raise KeyError(f"Expected type(list) or type(str) for group_vars, "
                       f"instead got {type(group_vars)}")

    assert time_var in adata.obs.columns, f"time_var not in .obs.columns: {time_var}"
    assert x_loc in adata.obs.columns, f"x_loc not in .obs.columns: {x_loc}"
    assert x_loc in adata.obs.columns, f"x_loc not in .obs.columns: {x_loc}"
    assert y_loc in adata.obs.columns, f"y_loc not in .obs.columns: {y_loc}"
    assert x_loc != y_loc, f"x_loc expected to be different from y_loc, " \
                           f"x_loc: {x_loc}, y_loc: {y_loc}"

    # create new column to store index of parent object
    adata.obs[parent_id] = np.nan
    # create new column and store id for every trace tree
    adata.obs[tree_id] = np.nan
    n_start_tp = len(adata[adata.obs[time_var] == start_tp])
    if n_start_tp > 1:
        adata.obs.loc[adata.obs[time_var] == start_tp, tree_id] = range(n_start_tp)
    else:
        raise ValueError(f"No observation with time_var {time_var} and start_tp {start_tp}")

    # iterate over every field and get field at different times
    total_its = np.prod([len(adata.obs[gv].unique()) for gv in group_vars])
    for ix, (groups, field_df) in tqdm(enumerate(adata.obs.groupby(list(group_vars))),
                                       desc='Tracking cells in every field',
                                       total=int(total_its)):

        # cache lagged values
        lagged = None

        # iterate over timepoints
        for t, t_df in field_df.groupby(time_var):

            # find closest object in lagged objects
            if lagged is not None:
                # get locations of objects and lagged objects
                t_loc = t_df[[x_loc, y_loc]].to_numpy()
                t_loc_lagged = lagged[[x_loc, y_loc]].to_numpy()

                # get nearest object in lagged objects for every object
                tree = cKDTree(t_loc_lagged)
                _, parent_ix = tree.query(t_loc, k=1)

                # assign lagged trace ids to objects
                adata.obs.loc[t_df.index, tree_id] = lagged.iloc[parent_ix][tree_id].tolist()
                # assign trace parents to objects
                adata.obs.loc[t_df.index, parent_id] = lagged.iloc[parent_ix].index

            # cache field_df
            lagged = adata.obs.loc[t_df.index, [x_loc, y_loc, tree_id]]

    return adata

max_stars_repo_name: marx-alex/Morphelia

max_issues_repo_name: marx-alex/Morphelia

max_forks_repo_name: marx-alex/Morphelia

reason: import numpy,from scipy


text:
from math import atan, cos, sqrt
from os import path
from matplotlib.colors import PowerNorm
from pyproj import Proj
from math import pi, sin, cos
import math
import numpy as np
import matplotlib.pyplot as plt

DEG_TO_RAD = math.pi / 180

class KinematicModel(object):
    def __init__(self, x, y, psi, v, f_len, r_len):
        self.x = x
        self.y = y
        self.psi = psi
        self.v = v

        self.f_len = f_len
        self.r_len = r_len

    def get_state(self):
        return self.x, self.y, self.psi, self.v

    def update_state(self, a, delta, dt):
        beta = math.atan((self.r_len / (self.r_len + self.f_len)) * math.tan(delta))

        self.x = self.x + self.v * math.cos(self.psi + beta) * dt
        self.y = self.y + self.v * math.sin(self.psi + beta) * dt
        self.psi = self.psi + (self.v / (self.f_len + self.r_len)) * math.sin(beta) * dt
        self.v = self.v + a * dt
        return self.x, self.y, self.psi, self.v

def read_global_path(file_log):
    """[summary]
    Args:
        file_log ([str]): log path

    Returns:
        [list]: [ [lng, lat, alt, yaw], ... ]  in deg
    """
    res = []
    vx = 0
    rf = open(file_log, "r")
    for line in rf:
        line = line.strip()
        fields = line.split(",")
        if "$GINS" in fields[0] and len(fields) > 10:
            lat, lng, alt = float(fields[3]), float(fields[4]), float(fields[5])
            yaw = (float(fields[11]) + 360) % 360
            res.append([lng, lat, alt, vx, yaw])

        if "$GNRMC" in fields[0] and len(fields) > 10:
            vx = float(fields[7]) * 0.774

    return res

def test(log_path):
    proj_text = "+proj=utm +zone=50 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
    trsfm = Proj(proj_text)
    path_point = read_global_path(log_path)

    x, y = [], []
    for point in path_point:
        utm_x, utm_y = trsfm(point[0], point[1])
        x.append(utm_x)
        y.append(utm_y)
    x = np.array(x)
    y = np.array(y)

    km = KinematicModel(0, 0, 0, 0, 0, 0)

    yaw_list = []
    lat_erros = []

    model_x = []
    model_y = []

    # 113.745098853 22.737577497
    # km.x, km.y = trsfm(113.745098853, 22.737577497)
    km.x, km.y = x[0], y[0]
    km.v =  10 / 3.6 # path_point[0][-2]  m/s
    km.psi = pi/2 - path_point[0][-1] * DEG_TO_RAD

    km.r_len = 1.5
    km.f_len = 1.5

    target_pts_idx = 0
    iter = 0

    kp, ki, kd = 0.01, 10, 0
    g_rate = 0.1

    last_error = -1
    last_last_error = -1

    while iter < 100000:
        if target_pts_idx >= len(x):
            break
        # km.v = path_point[target_pts_idx][-2]
        iter += 1
        ld = max(km.v * g_rate, 2)
        dx = km.x - x[target_pts_idx]
        dy = km.y - y[target_pts_idx]
        euler_dist = math.sqrt(dx**2 + dy**2)

        if euler_dist < ld:
            target_pts_idx += 1
            continue

        print(target_pts_idx, ld, euler_dist, dx, dy, km.x, km.y)

        cos_target_heading = cos(path_point[target_pts_idx][-1] * DEG_TO_RAD)
        sin_target_heading = sin(path_point[target_pts_idx][-1] * DEG_TO_RAD)

        # lateral_error = cos_target_heading * dy - sin_target_heading * dx
        # print("lateral_error", lateral_error)
        lateral_error = cos_target_heading * dx - sin_target_heading * dy
        longtidual_error = sin_target_heading * dx + cos_target_heading * dy

        print("lateral_error", lateral_error)
        # print("longtidual_error", longtidual_error)
        # print("euler_dist:", math.sqrt(dx**2 + dy**2))

        alpha = atan(2 * 2.9 * lateral_error / euler_dist**2) + kp * lateral_error + ki * (lateral_error - last_error) + kd * (lateral_error - 2*(last_error) + last_last_error)
        alpha = min(alpha, 30 * DEG_TO_RAD)
        alpha = max(alpha, -30 * DEG_TO_RAD)
        last_last_error = last_error
        last_error = lateral_error
        # print("ld:", ld)
        print("alpha:", alpha)
        km.update_state(0, alpha, 0.01)
        model_x.append(km.x)
        model_y.append(km.y)

    # plt.plot(x, lat_erros)
    # plt.plot(x, yaw_list)
    plt.scatter(model_x, model_y, label="model", s=1)
    plt.scatter(x, y, label="true", s=1)
    plt.legend()
    plt.show()

def main():
    # zone mean utm zone , beijing is 50
    proj_text = "+proj=utm +zone=50 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
    P = Proj(proj_text)

    point_init = (22.534624368,113.949931864, 5.720, -173.04 * DEG_TO_RAD)
    point_current = (22.526802456, 113.949383126, 6.151, -173.04 * DEG_TO_RAD)
    point_target = (22.529278272, 113.949620569, 5.518, pi/3)

    utm_init = P(point_init[1], point_init[0],)
    utm_current = P(point_current[1], point_current[0])
    utm_target = P(point_target[1], point_target[0])

    print("utm_init:", utm_init)
    print("utm_current:", utm_current)
    print("utm_target:", utm_target)

    x_diff_map = utm_current[0] - utm_init[0]
    y_diff_map = utm_current[1] - utm_init[1]
    theta_diff = point_current[-1] - point_init[-1]

    cos_map_veh = cos(point_init[-1])
    sin_map_veh = sin(point_init[-1]);

    x_diff_veh = cos_map_veh * x_diff_map + sin_map_veh * y_diff_map;
    y_diff_veh = -sin_map_veh * x_diff_map + cos_map_veh * y_diff_map;

    print("current_veh_x:", x_diff_veh)
    print("current_veh_y:", y_diff_veh)

    cos_theta_diff = cos(-theta_diff)
    sin_theta_diff = sin(-theta_diff)

    tx = -(cos_theta_diff * x_diff_veh - sin_theta_diff * y_diff_veh)
    ty = -(sin_theta_diff * x_diff_veh + cos_theta_diff * y_diff_veh)

    print(tx, ty)

    x = utm_target[0]
    y = utm_target[1]
    theta = point_target[-1]

    x_new = cos_theta_diff * x - sin_theta_diff * y + tx
    y_new = sin_theta_diff * x + cos_theta_diff * y + ty
    theta_new = theta - theta_diff

    print("target_new_x:", x_new)
    print("target_new_y:", y_new)

    dx = utm_current[0] - utm_target[0]
    dy = utm_current[1] - utm_target[1]

    print("dx:", dx, "dy:", dy)

    cos_target_heading = cos(point_target[-1])
    sin_target_heading = sin(point_target[-1])

    lateral_error = cos_target_heading * dy - sin_target_heading * dx
    print("lateral_error", lateral_error)
    lateral_error = cos_target_heading * dx - sin_target_heading * dy
    print("lateral_error", lateral_error)

if __name__ == "__main__":

    # x = []
    # y = []
    # km = KinematicModel(0, 0, pi/3, 5, 0.5 ,0.5)
    # for i in range(200):
    #     km.update_state(0, 0, 0.1)
    #     x.append(km.x)
    #     y.append(km.y)
    #     plt.plot(x, y)
    #     plt.pause(0.1)
    test("cutecom.log")
    # main()

max_stars_repo_name: fulincao/algorithm-review

max_issues_repo_name: fulincao/algorithm-review

max_forks_repo_name: fulincao/algorithm-review


text:
#!/usr/bin/env python

import numpy as np
from sklearn.datasets import load_digits
import pytest
import sys
sys.path.append('../') # not os agnostic

from mymlimp.classification import knn

X, y = load_digits(return_X_y=True)

def test_fit():
    """
    Tests that the fir method loads data correctly
    """
    model = knn(k=3)
    model.fit(X,y)
    assert model.X_train.shape == X.shape
    assert model.y_train.shape == y.shape

def test_internal_prediction():
    """
    Test the return of the interval single observation predict @staticmethod
    """
    obs = X[0,:]
    k = 3
    prediction = knn._pred(obs, y, X, k)
    assert prediction.shape == ()

def test_external_prediction():
    """
    Tests the shape and contents of external prediction method
    """
    model = knn(k=3)
    model.fit(X,y)
    prediction = model.predict(X)
    assert prediction.shape == y.shape
    assert len(np.unique(prediction)) <= len(np.unique(y))

max_stars_repo_name: bradleypick/ML-implementations

max_issues_repo_name: bradleypick/ML-implementations

max_forks_repo_name: bradleypick/ML-implementations


text:
from collections import defaultdict
import pickle
import astropy.coordinates as coord
import astropy.units as u
import numpy as np
from scipy.interpolate import interp1d
from scipy.spatial import cKDTree
from scipy.stats import binned_statistic

__all__ = ['AbundanceTorusMaschine', 'run_bootstrap_coeffs',
           'get_cos2th_zerocross']

class AbundanceTorusMaschine:

    def __init__(self, aaf, tree_K=64, sinusoid_K=2):
        """
        Parameters
        ----------
        aaf : `astropy.table.Table`
            Table of actions, angles, and abundances.
        tree_K : int (optional)
            The number of neighbors used to estimate the action-local
            mean abundances.
        sinusoid_K : int (optional)
            The number of cos/sin terms in the sinusoid fit to the
            abundance anomaly variations with angle.
        """

        self.aaf = aaf

        # config
        self.tree_K = int(tree_K)
        self.sinusoid_K = int(sinusoid_K)

#     def get_theta_z_anomaly(self, elem_name, action_unit=30*u.km/u.s*u.kpc):
#         action_unit = u.Quantity(action_unit)

#         # Actions without units:
#         X = self.aaf['actions'].to_value(action_unit)
#         angz = coord.Angle(self.aaf['angles'][:, 2]).wrap_at(360*u.deg).radian

#         # element abundance
#         elem = self.aaf[elem_name]
#         elem_errs = self.aaf[f"{elem_name}_ERR"]
#         ivar = 1 / elem_errs**2

#         tree = cKDTree(X)
#         dists, idx = tree.query(X, k=self.tree_K+1)

#         # compute action-local abundance anomaly
#         errs = np.sqrt(1 / np.sum(ivar[idx[:, 1:]], axis=1))
#         means = np.sum(elem[idx[:, 1:]] * ivar[idx[:, 1:]], axis=1) * errs**2

#         d_elem = elem - means
#         d_elem_errs = np.sqrt(elem_errs**2 + errs**2)
# #         d_elem_errs = np.full_like(d_elem, 0.04)  # MAGIC NUMBER

#         return angz, d_elem, d_elem_errs

    def get_theta_z_anomaly(self, elem_name, action_unit=30*u.km/u.s*u.kpc):
        action_unit = u.Quantity(action_unit)

        # Actions without units:
        X = self.aaf['actions'].to_value(action_unit)
        angz = coord.Angle(self.aaf['angles'][:, 2]).wrap_at(360*u.deg).radian

        # element abundance
        elem = self.aaf[elem_name]
        elem_errs = self.aaf[f"{elem_name}_ERR"]

        tree = cKDTree(X)
        dists, idx = tree.query(X, k=self.tree_K+1)

        xhat = np.mean(X[idx[:, 1:]], axis=1) - X
        dx = X[idx[:, 1:]] - X[:, None]
        x = np.einsum('nij,nj->ni', dx, xhat)
        y = elem[idx[:, 1:]]

        w = np.sum(x**2, axis=1)[:, None] - x * np.sum(x, axis=1)[:, None]
        means = np.sum(y * w, axis=1) / np.sum(w, axis=1)

        d_elem = np.array(elem - means)
        d_elem_errs = np.array(elem_errs)

        return angz, d_elem, d_elem_errs

    def get_M(self, x):
        M = np.full((len(x), 1 + 2*self.sinusoid_K), np.nan)
        M[:, 0] = 1.

        for n in range(self.sinusoid_K):
            M[:, 1 + 2*n] = np.cos((n+1) * x)
            M[:, 2 + 2*n] = np.sin((n+1) * x)

        return M

    def get_coeffs(self, M, y, yerr):
        Cinv_diag = 1 / yerr**2
        MT_Cinv = M.T * Cinv_diag[None]
        MT_Cinv_M = MT_Cinv @ M
        coeffs = np.linalg.solve(MT_Cinv_M, MT_Cinv @ y)
        coeffs_cov = np.linalg.inv(MT_Cinv_M)
        return coeffs, coeffs_cov

    def get_coeffs_for_elem(self, elem_name):
        tz, d_elem, d_elem_errs = self.get_theta_z_anomaly(elem_name)
        M = self.get_M(tz)
        return self.get_coeffs(M, d_elem, d_elem_errs)

    def get_binned_anomaly(self, elem_name, theta_z_step=5*u.deg):
        """
        theta_z_step : `astropy.units.Quantity` [angle] (optional)
            The bin step size for the vertical angle bins. This is only
            used in methods if `statistic != 'mean'`.
        """
        theta_z_step = coord.Angle(theta_z_step)
        angz_bins = np.arange(0, 2*np.pi+1e-4,
                              theta_z_step.to_value(u.rad))
        theta_z, d_elem, d_elem_errs = self.get_theta_z_anomaly(elem_name)
        d_elem_ivar = 1 / d_elem_errs**2
#         d_elem_ivar = np.full_like(d_elem, 1 / 0.04**2)  # MAGIC NUMBER

        stat1 = binned_statistic(theta_z, d_elem * d_elem_ivar,
                                 bins=angz_bins,
                                 statistic='sum')
        stat2 = binned_statistic(theta_z, d_elem_ivar,
                                 bins=angz_bins,
                                 statistic='sum')

        binx = 0.5 * (angz_bins[:-1] + angz_bins[1:])
        means = stat1.statistic / stat2.statistic
        errs = np.sqrt(1 / stat2.statistic)

        return binx, means, errs

def run_bootstrap_coeffs(aafs, elem_name, bootstrap_K=128, seed=42,
                         overwrite=False, cache_path=None):
    if cache_path is not None:
        cache_filename = f'coeffs-bootstrap{bootstrap_K}-{elem_name}.pkl'
        coeffs_cache = cache_path / cache_filename

        if coeffs_cache.exists() and not overwrite:
            with open(coeffs_cache, 'rb') as f:
                all_bs_coeffs = pickle.load(f)
            return all_bs_coeffs

    all_bs_coeffs = {}
    for name in aafs:
        aaf = aafs[name]

        if seed is not None:
            np.random.seed(seed)

        bs_coeffs = []
        for k in range(bootstrap_K):
            bootstrap_idx = np.random.choice(len(aaf), size=len(aaf))
            atm = AbundanceTorusMaschine(aaf[bootstrap_idx])
            coeffs, _ = atm.get_coeffs_for_elem(elem_name)
            bs_coeffs.append(coeffs)

        all_bs_coeffs[name] = np.array(bs_coeffs)

    if cache_path is not None:
        with open(coeffs_cache, 'wb') as f:
            pickle.dump(all_bs_coeffs, f)

    return all_bs_coeffs

def get_cos2th_zerocross(coeffs):
    summary = defaultdict(lambda *args: defaultdict(list))
    for i in range(5):
        for k in coeffs:
            summary[i]['mdisk'].append(float(k))
            summary[i]['y'].append(np.mean(coeffs[k][:, i]))
            summary[i]['y_err'].append(np.std(coeffs[k][:, i]))

        summary[i]['mdisk'] = np.array(summary[i]['mdisk'])
        summary[i]['y'] = np.array(summary[i]['y'])
        summary[i]['y_err'] = np.array(summary[i]['y_err'])

        idx = summary[i]['mdisk'].argsort()
        for key in summary[i].keys():
            summary[i][key] = summary[i][key][idx]

    # cos2theta term:
    s = summary[3]
    zero_cross = interp1d(s['y'], s['mdisk'], fill_value="extrapolate")(0.)
    zero_cross1 = interp1d(np.array(s['y']) - np.array(s['y_err']),
                           s['mdisk'], fill_value="extrapolate")(0.)
    zero_cross2 = interp1d(np.array(s['y']) + np.array(s['y_err']),
                           s['mdisk'], fill_value="extrapolate")(0.)

    zero_cross_err = (np.abs(zero_cross2 - zero_cross),
                      np.abs(zero_cross1 - zero_cross))
    return summary, zero_cross, zero_cross_err

class ZeroCrossWorker:

    def __init__(self, aafs, cache_path=None, bootstrap_K=128):
        self.aafs = aafs
        self.cache_path = cache_path
        self.bootstrap_K = int(bootstrap_K)

    def __call__(self, elem_name):
        clean_aafs = {}
        for k in self.aafs:
            clean_aafs[k] = self.aafs[k][self.aafs[k][elem_name] > -3]

        bs_coeffs = run_bootstrap_coeffs(clean_aafs, elem_name,
                                         bootstrap_K=self.bootstrap_K,
                                         cache_path=self.cache_path)
        s, zc, zc_err = get_cos2th_zerocross(bs_coeffs)
        return elem_name, [zc, zc_err]

max_stars_repo_name: adrn/thrift-shop

max_issues_repo_name: adrn/thrift-shop

max_forks_repo_name: adrn/thrift-shop

reason: import numpy,from scipy,import astropy