text:
#!/usr/bin/python3

import sys
import numpy as np
import json
from scipy import stats
from nnweaver.validation import random_search, kfold_cross_validation
from nnweaver.nn import NN, Layer
from nnweaver.activations import Rectifier, TanH, Linear
from nnweaver.optimizers import ProximalBundleMethod
from nnweaver.losses import MEE
from nnweaver.regularizers import L1L2Regularizer

class LogScaleUniform(object):
    def __init__(self, low, high):
        self.log_low = np.log10(low)
        self.log_high = np.log10(high)
        self.dist = stats.uniform(0, 1)

    def rvs(self, random_state=None):
        rv = self.log_low + self.dist.rvs(random_state=random_state) * (
                self.log_high - self.log_low)
        return 10 ** rv

def nn_builder():
    nn = NN(10)
    nn.add_layer(Layer(97, Rectifier))
    nn.add_layer(Layer(19, TanH))
    nn.add_layer(Layer(2, Linear))
    return nn

train_args = {'mu': LogScaleUniform(0.1, 10),
              'm_L': stats.uniform(0.1, 0.4),
              'm_R': stats.uniform(0.5, 0.5),
              't_bar': stats.uniform(0.1, 0.9),
              'a_bar': LogScaleUniform(0.1, 100.),
              'regularizer': [L1L2Regularizer(0.00037205665611222174,
                                              0.0008170173137216012)],
              'gamma': LogScaleUniform(1e-8, 10.),
              'accuracy_tolerance': [1e-8],
              'max_iterations': [200]}

assert len(sys.argv) > 1, 'Dataset path not given'
dataset = np.genfromtxt(sys.argv[1] + '/train.csv', delimiter=',',
                        usecols=tuple(range(1, 13)))
x = dataset[:, :10]
y = dataset[:, 10:]
limit = int(len(x) * 0.8)
x_train = x[:limit]
y_train = y[:limit]
x_val = x[limit:]
y_val = y[limit:]

pbm = ProximalBundleMethod(MEE)
results = []
i = 0

while i < 10:
    try:
        _, loss, args, _ = random_search(nn_builder, pbm, x_train, y_train,
                                         train_args, {}, 1)
        args['regularizer'] = {"l1": args['regularizer'].l1,
                               "l2": args['regularizer'].l2}
        element = dict()
        element["loss"] = loss
        element["args"] = args
        results.append(element)
        print("\nResults\n"
              "=======\n"
              "Loss: {}\n"
              "Parameters: {}\n\n".format(str(loss), str(args)))
    except KeyboardInterrupt:
        break
    except ValueError:
        continue
    i += 1

file_path = '%s/results.json' % sys.argv[1]

try:
    with open(file_path, 'r') as file:
        experiments = json.load(file)
except FileNotFoundError:
    experiments = []

with open(file_path, 'w+') as file:
    experiments.extend(results)
    json.dump(experiments, file, indent=3)

max_stars_repo_name: gvinciguerra/nnweaver

max_issues_repo_name: gvinciguerra/nnweaver

max_forks_repo_name: gvinciguerra/nnweaver

reason: import numpy,from scipy


text:
import sys

import cv2
import numpy as np

from je_auto_control.utils.image.screenshot import pil_screenshot
from je_auto_control.utils.exception.exception_tag import screen_get_size
from je_auto_control.utils.exception.exception_tag import screen_screenshot
from je_auto_control.utils.exception.exceptions import AutoControlScreenException
from je_auto_control.wrapper.platform_wrapper import screen
from je_auto_control.utils.test_record.record_test_class import record_total

def size():
    """
    get screen size
    """
    try:
        try:
            record_total("size", None)
            return screen.size()
        except AutoControlScreenException:
            raise AutoControlScreenException(screen_get_size)
    except Exception as error:
        record_total("size", None, repr(error))
        print(repr(error), file=sys.stderr)

def screenshot(file_path: str = None, region: list = None):
    """
    use to capture current screen image
    :param file_path screenshot file save path
    :param region screenshot region
    """
    param = locals()
    try:
        try:
            record_total("screenshot", param)
            return cv2.cvtColor(np.array(pil_screenshot(file_path=file_path, region=region)), cv2.COLOR_RGB2BGR)
        except AutoControlScreenException as error:
            raise AutoControlScreenException(screen_screenshot + " " + repr(error))
    except Exception as error:
        record_total("screenshot", None, repr(error))
        print(repr(error), file=sys.stderr)

ext: py

max_stars_repo_name: JE-Chen/AutoControl

max_issues_repo_name: JE-Chen/AutoControl

max_forks_repo_name: JE-Chen/AutoControl

reason: import numpy


text:
import logging
import os
import shutil
import subprocess
import json
import requests
from retrying import retry
import numpy as np

KFP_PACKAGE = 'https://storage.googleapis.com/ml-pipeline/release/0.1.20/kfp.tar.gz'
def notebook_setup():
  # Install the SDK

  subprocess.check_call(["pip3", "install", "-r", "requirements.txt"])
  subprocess.check_call(["pip3", "install", KFP_PACKAGE, "--upgrade"])

  logging.basicConfig(format='%(message)s')
  logging.getLogger().setLevel(logging.INFO)

  subprocess.check_call(["gcloud", "auth", "configure-docker", "--quiet"])
  subprocess.check_call(["gcloud", "auth", "activate-service-account",
                         "--key-file=" + os.getenv("GOOGLE_APPLICATION_CREDENTIALS"),
                         "--quiet"])

def copy_data_to_nfs(nfs_path, model_dir):
  if not os.path.exists(nfs_path):
    shutil.copytree("ames_dataset", nfs_path)

  if not os.path.exists(model_dir):
    os.makedirs(model_dir)

@retry(wait_exponential_multiplier=1000, wait_exponential_max=5000,
       stop_max_delay=2*60*1000)
def predict_nparray(url, data, feature_names=None):
  pdata = {
      "data": {
          "names":feature_names,
          "tensor": {
              "shape": np.asarray(data.shape).tolist(),
              "values": data.flatten().tolist(),
          },
      }
  }
  serialized_data = json.dumps(pdata)
  r = requests.post(url, data={'json':serialized_data}, timeout=5)
  return r

max_stars_repo_name: rainerenglisch/examples

max_issues_repo_name: rainerenglisch/examples

max_forks_repo_name: rainerenglisch/examples


text:
import datetime
import numpy as np
from icarus_backend.flight.FlightModel import Flight
from users.models import IcarusUser as User
from icarus_backend.department.DepartmentModel import Department
from icarus_backend.department.tasks import DepartmentTasks
from django.utils import timezone
from django.contrib.gis.geos import Polygon
from django.db.models import Q

class DepartmentController:

    @staticmethod
    def create(name, owner_id, area, airbosses=[], watch_commanders=[]) -> (int, dict):
        numpy_area = np.asarray(area)
        if not (len(numpy_area.shape) == 3 and numpy_area.shape[2] == 2):
            return 400, {'message': 'Area is not properly formatted. Must have shape (-1, -1, 2).'}
        area = [[x[1],x[0]] for x in area[0]]
        area = Polygon(area)
        department = Department.objects.filter(name=name).first()
        if department:
            return 400, {'message': 'Department name already taken.'}
        owner = User.objects.filter(id=owner_id).first()
        department = Department(name=name, area=area, owner=owner)
        for airboss in airbosses:
            department.airbosses.add(airboss)
        for watch_commander in watch_commanders:
            department.watchCommanders.add(watch_commander)
        department.save()
        return 200, {'message': 'Department created.'}

    @staticmethod
    def edit(name, new_name, area) -> (int, dict):
        area = Polygon(area)
        department = Department.objects.filter(name=name).first()
        if not department:
            return 400, {'message': 'No department with this name exists.'}
        department.area = area
        department.name = new_name
        department.save()
        return 200, {'message': 'Department edited.'}

    @staticmethod
    def get(area=None) -> (int, dict):
        departments = Department.objects
        if area:
            area = Polygon(area)
            departments = departments.filter(area__intersects=area)
        departments = departments.all()
        response_array = []
        for department in departments:
            response_array += [department.as_dict()]
        return 200, response_array

    @staticmethod
    def add_airboss(name, airboss_id) -> (int, dict):
        department = Department.objects.filter(name=name).first()
        if not department:
            return 400, {'message': 'No department of that name exists.'}

        airboss = User.objects.filter(id=airboss_id).first()
        if not airboss:
            return 400, {'message': 'No user of that id exists.'}
        if department.airbosses.filter(id=airboss_id).exists():
            return 400, {'message': 'This user is already an airboss for this department.'}
        department.airbosses.add(airboss)
        return 200, {'message': 'Members successfully added.'}

    @staticmethod
    def remove_airboss(name, airboss_id) -> (int, dict):
        department = Department.objects.filter(name=name).first()
        if not department:
            return 400, {'message': 'No department of that name exists.'}

        airboss = User.objects.filter(id=airboss_id).first()
        if not airboss:
            return 400, {'message': 'No user of that id exists.'}
        if not department.airbosses.filter(id=airboss_id).exists():
            return 400, {'message': 'This user is not currently an airboss for this department.'}
        department.airbosses.remove(airboss)
        return 200, {'message': 'Members successfully removed.'}

    @staticmethod
    def add_watch_commander(name, watch_commander_id) -> (int, dict):
        department = Department.objects.filter(name=name).first()
        if not department:
            return 400, {'message': 'No department of that name exists.'}

        watch_commander = User.objects.filter(id=watch_commander_id).first()
        if not watch_commander:
            return 400, {'message': 'No user of that id exists.'}
        if department.airbosses.filter(id=watch_commander_id).exists():
            return 400, {'message': 'This user is already a watch commander for this department.'}
        department.watchCommanders.add(watch_commander)
        return 200, {'message': 'Members successfully added.'}

    @staticmethod
    def remove_watch_commander(name, watch_commander_id) -> (int, dict):
        department = Department.objects.filter(name=name).first()
        if not department:
            return 400, {'message': 'No department of that name exists.'}

        watch_commander = User.objects.filter(id=watch_commander_id).first()
        if not watch_commander:
            return 400, {'message': 'No user of that id exists.'}
        if not department.watchCommanders.filter(id=watch_commander_id).exists():
            return 400, {'message': 'This user is not currently an airboss for this department.'}
        department.watchCommanders.remove(watch_commander)
        return 200, {'message': 'Members successfully added.'}

    @staticmethod
    def flight_histogram(department_id, start_day, end_day, user) -> (int, dict):
        number_of_days = (end_day-start_day).days + 1
        histogram = []
        department = Department.objects.filter(id=department_id).first()
        if not department:
            return 400, {'message': 'No department with that id exists.'}
        for index in range(0, number_of_days):
            starts_at = start_day + datetime.timedelta(days=index)
            ends_at = start_day + datetime.timedelta(days=index+1)
            num_flights = Flight.objects.filter(area__intersects=department.area, starts_at__gt=starts_at, ends_at__lt=ends_at).count()
            histogram += [num_flights]
        return 200, histogram

    @staticmethod
    def message_pilots(department_name, message):
        _now = timezone.now()
        department = Department.objects.filter(name=department_name).first()
        local_pilots = User.objects.filter(flight__area__intersects=department.area, flight__ends_at__gt=_now)\
            .distinct()
        for pilot in local_pilots:
            DepartmentTasks.email_jurisdiction.delay(pilot.username, pilot.email, department.name, message)

    @staticmethod
    def get_user_department_query(user):
        return Department.objects.filter(Q(airbosses__in=[user]) | Q(watchCommanders__in=[user])).distinct()

    @staticmethod
    def get_user_department(user) -> (int, dict):
        departments_query = DepartmentController.get_user_department_query(user)
        departments = [x.as_dict() for x in departments_query]
        return 200, departments

    @staticmethod
    def is_department_member(user):
        return bool(DepartmentController.get_user_department_query(user))

    @staticmethod
    def can_edit(user, flight):
        can_edit = False
        for department in DepartmentController.get_user_department_query(user):
            if flight.area.intersects(department.area) and user in department.airbosses.all():
                can_edit = True
        return can_edit

    @staticmethod
    def info(department_id) -> (int, dict):
        department = Department.objects.filter(id=department_id).first()
        if not department:
            return 400, {'message': 'No department with this id exists.'}
        return 200, department.as_dict()

    @staticmethod
    def delete(department_id) -> (int, dict):
        department = Department.objects.filter(id=department_id).first()
        if not department:
            return 400, {'message': 'No department with this id exists.'}
        department.delete()
        return 200, {'message': 'Department successfully deleted.'}

max_issues_repo_name: ransbymich/FlyRight

max_forks_repo_name: ransbymich/FlyRight


text:
""" testdata for testcases """
import pandas as pd
import numpy as np
import os

testdata_path = os.path.join(__file__, "../../tests/data")

def linear_univariable(ss=True):
    predictor, response, subsample = __retrieve_linear_data(location="40.00km")
    if ss:
        return predictor, response, subsample
    else:
        return predictor, response

def __retrieve_linear_data(location=None):
    """ Testdata from intervention cases. Univariable (here: at one location)"""

    subsample = list(map(int, np.loadtxt('{}/intervention/subsample.csv'.format(testdata_path), delimiter=",")))
    predictor = pd.read_csv('{}/intervention/reference.csv'.format(testdata_path), index_col=0, header=0)
    response = predictor.copy()
    response.iloc[:] = np.nan
    response_data = pd.read_csv('{}/intervention/intervention.csv'.format(testdata_path), index_col=0, header=0)

    for column in response_data.columns:
        response.loc[:, column] = response_data.loc[:, column]

    # Now, select only one station
    if location is None:
        predictor = predictor.T
        response = response.T
    else:
        predictor = predictor.T[location]
        response = response.T[location]

    return predictor, response, subsample

max_stars_repo_name: kdberends/coral

max_issues_repo_name: kdberends/coral

max_forks_repo_name: kdberends/coral


text:
import h5py
from IceWaterPathMethod import IceWaterPathMethod
from GAN_generator import GAN_generator

import torch
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize

folder = './' # location of parameter file
file= 'network_parameters_CGAN.pt' # name of parameter file
checkpoint_parameter = torch.load(folder + file, map_location=torch.device('cpu'))
noise_parameter = checkpoint_parameter['noise_parameter']

print(noise_parameter)
H_gen=[576,16384, 256, 128, 64, 1]
device = torch.device("cuda:0" if (torch.cuda.is_available()) else "cpu")
epoch = checkpoint_parameter['epoch']
real_label = 1
fake_label = 0
beta1 = 0.5
criterion = torch.nn.BCELoss()
lr = 0.0002

f_size = 14

netG = GAN_generator(H_gen).float().to(device)
optimizerG = torch.optim.Adam(netG.parameters(), lr=lr, betas=(beta1, 0.999))
netG.load_state_dict(checkpoint_parameter['model_state_dict_gen'])

b_size = 1
D_in_gen = [b_size, 1,1,64]

def t_count(t, val):
    elements_equal_to_value = torch.eq(t, val)
    as_ints = elements_equal_to_value.numpy()
    count = np.sum(as_ints)
    return count

location = './modis_cloudsat_data/'
file_string = location + 'CGAN_test_data_with_temp_conc_ver2.h5' #for CGAN 2015
#file_string = location + 'modis_cloudsat_ElevLatLong_test_data_conc_2016_ver2.h5' #for CGAN 2016
hf = h5py.File(file_string, 'r')

cloudsat_scenes = torch.tensor(hf.get('cloudsat_scenes'))

modis_scenes = torch.tensor(hf.get('modis_scenes'))
print('modis scenes input size: ', modis_scenes.shape)
temp_modis_scenes = torch.cat([modis_scenes[:, :, :, 0:3], modis_scenes[:, :, :, 4:9]], 3)
modis_scenes = temp_modis_scenes
print('modis after change size: ', modis_scenes.shape)
temperature = torch.tensor(hf.get('temperature'))
dataset = torch.utils.data.TensorDataset(cloudsat_scenes, modis_scenes,temperature)

print('number of missing values in modis ', t_count(modis_scenes,(float(-1))))

xplot=np.linspace(0,64*1.1,64)
yplot=np.linspace(1,16.36,64)

norm = Normalize(-35, 20)
f, axs = plt.subplots(8,5, figsize=(13,18))
dataloader = torch.utils.data.DataLoader(dataset, batch_size=1, shuffle=True,
                                             num_workers=2)
file_number = 0
counter = 0
for i, data in enumerate(dataloader,0):
    if i > 4:
        break
    cs = data[0]

    print('cs first element: ',cs[0,0,0,0])
    string_ex = 'example' + str(i)
    modis = data[1]
    temperature_data = data[2]
    print(temperature_data.shape, 'shape of temperature data')

    modis = torch.transpose(modis, 1, 3)
    modis = torch.transpose(modis, 2, 3)
    print('modis shape', modis.shape)
    modis = modis.to(device)
    pl = axs[2,i].plot(xplot,np.transpose(np.array(modis[0,:,0,:])))
    axs[2,i].set_ylim([-1,1])
    axs[2,i].minorticks_on()
    axs[2, i].tick_params(axis='both', which='major', labelsize=f_size)
    axs[2, i].tick_params(axis='both', which='minor')
    axs[2,i].tick_params(labelbottom=False)
    if i == 0:
        axs[2, i].set_ylabel('Modis Input', fontsize=f_size)
    else:
        axs[2, i].tick_params(labelleft=False)

    cs = (cs + 1) * (55 / 2) - 35

    # Calculate cloud top height:
    altitudes_cs = np.zeros([1, 64])
    for k in range(0, 64):
        first_cloud_location_cs = 0
        for j in range(63, -1, -1):  # Loop backwards, starting from top of scene
            if cs[0,0, k, j] >= -15 and first_cloud_location_cs == 0:  # Set dBZ limit for cloud top detection
                altitudes_cs[0, k] = j  # save the index of the altitude where the cloud top is for each position
                first_cloud_location_cs = +1
    altitudes_cs = (altitudes_cs * 0.24) + 1  # Altitude of cloud top over sea level, [km]

    # Calculate indices for IWP
    indexes_zero = np.zeros([len(temperature_data), 1, 64, 1])

    for l in range(0, len(temperature_data)):
        for position in range(0, 64):
            for j in range(len(temperature_data[l, 0, position]) - 2, -1, -1):
                if (temperature_data[l, 0, position, j] - 273.15) * (temperature_data[l, 0, position, j + 1] - 273.15) <= 0:
                    indexes_zero[l, 0, position] = j + 1
                    break
                if j == 0:
                    if torch.max(temperature_data[l, 0, position, :]) <= 273.15:
                        indexes_zero[l, 0, position] = 63
                    else:
                        indexes_zero[l] = 0
                        print('All above zero: ', temperature_data[l, 0, position])

    indexes_zero_for_cs = np.ones([len(indexes_zero), 1, 64, 1]) * 63 - indexes_zero

    [IWP_cs_1, IWP_cs_2, IWP_cs_3] = IceWaterPathMethod(cs,indexes_zero_for_cs)
    cs = cs.to(device)
    pcm1 = axs[7, i].pcolormesh(xplot, yplot, np.transpose(cs.cpu().detach().numpy()[0,0]), norm=norm)
    title_str = 'Scene' + str(i)
    axs[7, i].minorticks_on()
    axs[7,i].set_xlabel('Position [km]', fontsize=f_size)
    if i == 0:
        axs[7, i].set_ylabel('Real \nAltitude [km]', fontsize=f_size)
    else:
        axs[7, i].tick_params(labelleft=False)
    axs[7, i].tick_params(axis='both', which='both', labelsize=f_size)

    color_array = ['blue','green','red','purple']
    for j in range(0,4):
        noise = (torch.randn(D_in_gen)).to(device)
        output = netG(noise, modis)
        output = (output + 1) * (55 / 2) - 35

        pcm = axs[7-(j+1), i].pcolormesh(xplot, yplot, output.detach().numpy()[0,0], norm=norm)
        title_str = 'Scene' + str(i)
        axs[7-(j+1), i].minorticks_on()
        axs[7-(j+1), i].tick_params(axis='both', which='major', labelsize=f_size)
        axs[7-(j+1), i].tick_params(axis='both', which='minor')
        axs[7-(j+1),i].tick_params(labelbottom=False)
        if i == 0:
            axs[7-(j+1), i].set_ylabel('Generated \nAltitude [km]', fontsize=f_size)
        else:
            axs[7 - (j + 1), i].tick_params(labelleft=False)

        with torch.no_grad():

            [IWP_generated_1, IWP_generated_2, IWP_generated_3] = IceWaterPathMethod(torch.transpose(output, 2, 3),indexes_zero_for_cs)
        pcm = axs[0, i].semilogy(xplot, (IWP_generated_1[0, 0]), color=color_array[j], label='IWP generated', basey=10, alpha=0.5)
        # Calculate cloud top height:
        altitudes_generated = np.zeros([1, 64])
        for k in range(0, 64):
            first_cloud_location_gen = 0
            for l in range(63, -1, -1):  # Loop backwards, starting from top of scene
                if output[0,0,l,k] >= -15 and first_cloud_location_gen == 0:  # Set dBZ limit for cloud top detection
                    altitudes_generated[0, k] = l  # save the index of the altitude where the cloud top is for each position
                    first_cloud_location_gen = +1
        altitudes_generated = (altitudes_generated * 0.24) + 1  # Altitude of cloud top over sea level, [km]

        alt = axs[1,i].plot(xplot, np.transpose(altitudes_generated), color = color_array[j], label ='CTH generated',alpha=0.3)
    alt_cs = axs[1,i].plot(xplot, np.transpose(altitudes_cs),color='black', label='CTH real')
    axs[1, i].minorticks_on()
    axs[1, i].tick_params(axis='both', which='major', labelsize=f_size)
    axs[1, i].tick_params(axis='both', which='minor')
    axs[1, i].tick_params(labelbottom=False)
    if i == 0:
        axs[1, i].set_ylabel('Cloud-Top \nHeight [km]', fontsize=f_size)
    else:
        axs[1, i].tick_params(labelleft=False)
    axs[1, i].set_ylim(1, 16.36)

    axs[0, i].minorticks_on()
    axs[0, i].tick_params(axis='both', which='major', labelsize=f_size)
    axs[0, i].tick_params(axis='both', which='minor')
    axs[0, i].set_ylim([10 ** (-2), 10 ** 1])

    pcm = axs[0, i].semilogy(xplot, (IWP_cs_1[0,0]), color='black', label='IWP real',basey=10)
    title_str = 'IWP generated' + str(i)
    axs[0, i].tick_params(labelbottom=False)
    if i == 0:
        axs[0, i].set_ylabel('IWP [g m⁻²]', fontsize=f_size)
    else:
        axs[0, i].tick_params(labelleft=False)
    #axs[0, i].set_aspect(1/4)

f.tight_layout()
f.subplots_adjust(right=0.88)
cbar_ax = f.add_axes([0.9, 0.0375, 0.025, 0.59]) #([distance from left, bottom height, width, height])
cbar1= f.colorbar(pcm1, cax=cbar_ax)
cbar1.set_label('Reflectivities [dBZ]', fontsize=f_size)
cbar1.ax.tick_params(labelsize=f_size)

print('image saved as: ', 'testepoch' + str(epoch) + '_CGAN')
plt.savefig('./early_epochs/cgan/testepoch' + str(epoch) + '_CGAN')


text:
import pandas as pd
import numpy as np
import math

#V0.1.1

def do_snp_qc(snp_df, min_call_rate, min_maf, min_hwe_P):

    #Determine call rate.
    call_rate = 1-snp_df.isnull().sum()/len(snp_df.index)
    #Return call_rate
    selection = call_rate < min_call_rate
    #print("Pass CR: ");print(snp_df.isnull().sum())
    #print(call_rate)
    #print(call_rate < min_call_rate)
    call_rate = pd.DataFrame(data=call_rate.values, index=snp_df.columns, columns=['call_rate'])
    failed_snp_names  = list(snp_df.columns[selection])
    snp_df_c = snp_df.loc[:,list(snp_df.columns[~selection])]
    #print("Pass CR: ");print(snp_df_c.shape)
    if(len(snp_df_c.columns)==0):
        return snp_df_c.columns, failed_snp_names, call_rate, None, None
    #Determine MAF.
    genotypeCounter = np.zeros((len(snp_df_c.columns),3), dtype=np.int)
    for allele_index in [0,1,2]:
        genotypeCounter[:,allele_index] = np.nansum(np.around(snp_df_c.values)==allele_index,0)

    #Here we make sure that the major allele is temporarly 'coded' as 0 & directly calculate the MAF (based on allele counts and non NA samples)

    mac = np.zeros((len(snp_df_c.columns)), dtype=np.int)
    gc = np.zeros((len(snp_df_c.columns)), dtype=np.int)
    maf = np.zeros((len(snp_df_c.columns)), dtype=np.float)

    for snp in range(0, len(snp_df_c.columns)):
        if genotypeCounter[snp,0]<genotypeCounter[snp,2]:
            tmp = genotypeCounter[snp,0]
            genotypeCounter[snp,0] = genotypeCounter[snp,2]
            genotypeCounter[snp,2] = tmp
        mac[snp] = int((genotypeCounter[snp,2]*2)+genotypeCounter[snp,1])
        gc[snp] = int(genotypeCounter[snp,2]+genotypeCounter[snp,1]+genotypeCounter[snp,0])
        maf[snp] = mac[snp] / (float(2*gc[snp]))
    #Return maf
    selection = maf < min_maf
    maf = pd.DataFrame(data=maf, index=snp_df_c.columns, columns=['maf'])
    failed_snp_names.extend(list(snp_df_c.columns[selection]))
    snp_df_c = snp_df_c.iloc[:,~selection]
    #print("Pass MAF: ");print(snp_df_c.shape)
    if(len(snp_df_c.columns)==0):
        return snp_df_c.columns, failed_snp_names, call_rate, maf, None

    mac=mac[~selection]
    gc=gc[~selection]
    genotypeCounter = genotypeCounter[~selection,]

    #Determine HWE.
    hweP = np.zeros((len(snp_df_c.columns)), dtype=np.float)
    #print(len(snp_df_c.columns))
    #print(len(mac))
    #This can also be multi-threaded if we place it in an separate function. (And multi-threading is as easy as it seems)
    for snp in range(0, len(snp_df_c.columns)):
        rare_copies = mac[snp]
        genotypes = gc[snp]

        het_probs = np.zeros((rare_copies+1))
        # start at midpoint #

        mid = int(math.floor(rare_copies * (2 * genotypes - rare_copies) / (2 * genotypes)))

        # check to ensure that midpoint and rare alleles have same parity #
        if int(mid % 2) != int(rare_copies % 2) :
            mid+=1

        curr_homr = int(math.floor((rare_copies - mid) / 2))
        curr_homc = genotypes - mid - curr_homr

        #print(het_probs.shape)
        het_probs[mid] = 1.0
        sum_values = het_probs[mid]
        for curr_hets in range(mid, 0, -2) :
            het_probs[curr_hets - 2] = het_probs[curr_hets] * curr_hets * (curr_hets - 1.0) / (4.0 * (curr_homr + 1.0) * (curr_homc + 1.0))
            sum_values += het_probs[curr_hets - 2]
            ## 2 fewer heterozygotes for next iteration -> add one rare, one common homozygote ##
            curr_homr+=1
            curr_homc+=1

        curr_homr = int(math.floor((rare_copies - mid) / 2))
        curr_homc = genotypes - mid - curr_homr

        for curr_hets in range(mid, (rare_copies), 2) :
            het_probs[curr_hets + 2] = het_probs[curr_hets] * 4.0 * curr_homr * curr_homc / ((curr_hets + 2.0) * (curr_hets + 1.0))
            sum_values += het_probs[curr_hets + 2]
            curr_homr-=1
            curr_homc-=1

        p_hwe = 0.0
        het_probs[int(genotypeCounter[snp,1])]/= sum_values
        for index in range(0,rare_copies+1) :
            if index != int(genotypeCounter[snp,1]) :
                het_probs[index] /= sum_values
            if het_probs[index] <= het_probs[int(genotypeCounter[snp,1])] :
                p_hwe += het_probs[index];

        hweP[snp] = 1 if p_hwe > 1.0 else p_hwe;
    #Return hweP
    selection = hweP < min_hwe_P
    hweP = pd.DataFrame(data=hweP, index=snp_df_c.columns, columns=['hwe_p'])
    failed_snp_names.extend(list(snp_df_c.columns[selection]))
    snp_df_c = snp_df_c.loc[:,list(snp_df_c.columns[~selection])]
    #print("Pass HWE: ");print(snp_df_c.shape)
    return snp_df_c.columns, failed_snp_names, call_rate, maf, hweP

def do_snp_qc_stringent(snp_df, min_call_rate, min_maf, min_hwe_P, min_hmachR2):

    #Determine call rate.
    call_rate = 1-snp_df.isnull().sum()/len(snp_df.index)
    selection = call_rate < min_call_rate
    #print(call_rate)
    #print(call_rate < min_call_rate)
    failed_snp_names  = list(snp_df.columns[selection])
    snp_df_c = snp_df.loc[:,list(snp_df.columns[~selection])]
    if(len(snp_df_c.columns)==0):
        return snp_df_c.columns, failed_snp_names
    #Determine MAF.
    genotypeCounter = np.zeros((len(snp_df_c.columns),3), dtype=np.int)
    for allele_index in [0,1,2]:
        genotypeCounter[:,allele_index] = np.nansum(np.around(snp_df_c.values)==allele_index,0)

    #Here we make sure that the major allele is temporarly 'coded' as 0 & directly calculate the MAF (based on allele counts and non NA samples)
    mac = np.zeros((len(snp_df_c.columns)), dtype=np.int)
    gc = np.zeros((len(snp_df_c.columns)), dtype=np.int)
    maf = np.zeros((len(snp_df_c.columns)), dtype=np.float)
    for snp in range(0, len(snp_df_c.columns)):
        if genotypeCounter[snp,0]<genotypeCounter[snp,2]:
            tmp = genotypeCounter[snp,0]
            genotypeCounter[snp,0] = genotypeCounter[snp,2]
            genotypeCounter[snp,2] = tmp
        mac[snp] = int((genotypeCounter[snp,2]*2)+genotypeCounter[snp,1])
        gc[snp] = int(genotypeCounter[snp,2]+genotypeCounter[snp,1]+genotypeCounter[snp,0])
        maf[snp] = mac[snp] / (float(2*gc[snp]))

    selection = maf < min_maf

    failed_snp_names.extend(list(snp_df_c.columns[selection]))
    snp_df_c = snp_df_c.loc[:,list(snp_df_c.columns[~selection])]

    if(len(snp_df_c.columns)==0):
        return snp_df_c.columns, failed_snp_names

    mac=mac[~selection]
    gc=gc[~selection]
    genotypeCounter = genotypeCounter[~selection,]

    #Determine HWE.
    hweP = np.zeros((len(snp_df_c.columns)), dtype=np.float)

    for snp in range(0, len(snp_df_c.columns)):
        rare_copies = mac[snp]
        genotypes = gc[snp]

        het_probs = np.zeros((rare_copies+1))
        # start at midpoint #

        mid = int(math.floor(rare_copies * (2 * genotypes - rare_copies) / (2 * genotypes)))

        # check to ensure that midpoint and rare alleles have same parity #
        if int(mid % 2) != int(rare_copies % 2) :
            mid+=1

        curr_homr = int(math.floor((rare_copies - mid) / 2))
        curr_homc = genotypes - mid - curr_homr

        #print(het_probs.shape)
        het_probs[mid] = 1.0
        sum_values = het_probs[mid]
        for curr_hets in range(mid, 0, -2) :
            het_probs[curr_hets - 2] = het_probs[curr_hets] * curr_hets * (curr_hets - 1.0) / (4.0 * (curr_homr + 1.0) * (curr_homc + 1.0))
            sum_values += het_probs[curr_hets - 2]
            ## 2 fewer heterozygotes for next iteration -> add one rare, one common homozygote ##
            curr_homr+=1
            curr_homc+=1

        curr_homr = int(math.floor((rare_copies - mid) / 2))
        curr_homc = genotypes - mid - curr_homr

        for curr_hets in range(mid, (rare_copies), 2) :
            het_probs[curr_hets + 2] = het_probs[curr_hets] * 4.0 * curr_homr * curr_homc / ((curr_hets + 2.0) * (curr_hets + 1.0))
            sum_values += het_probs[curr_hets + 2]
            curr_homr-=1
            curr_homc-=1

        p_hwe = 0.0
        het_probs[int(genotypeCounter[snp,1])]/= sum_values
        for index in range(0,rare_copies+1) :
            if index != int(genotypeCounter[snp,1]) :
                het_probs[index] /= sum_values
            if het_probs[index] <= het_probs[int(genotypeCounter[snp,1])] :
                p_hwe += het_probs[index];

        hweP[snp] = 1 if p_hwe > 1.0 else p_hwe;
    selection = hweP < min_hwe_P
    failed_snp_names.extend(list(snp_df_c.columns[selection]))
    snp_df_c = snp_df_c.loc[:,list(snp_df_c.columns[~selection])]

    if(len(snp_df_c.columns)==0):
        return snp_df_c.columns, failed_snp_names

    gc = gc[~selection]
    genotypeCounter = genotypeCounter[~selection,]

    machR2 = np.zeros((len(snp_df_c.columns)), dtype=np.float)
    for snp in range(0, len(snp_df_c.columns)):
        dosageSum = 0.0
        dosageSqrSum = 0.0
        dosageSum += genotypeCounter[snp,1]
        dosageSum += genotypeCounter[snp,2]*2
        dosageSqrSum += math.pow(1, 2)*genotypeCounter[snp,1]
        dosageSqrSum += math.pow(2, 2)*genotypeCounter[snp,2]
        nonMissingCount = gc[snp]
        estimatedAlleleFrequency = dosageSum / (2 * nonMissingCount)
        if (estimatedAlleleFrequency <= 0 or estimatedAlleleFrequency >= 1) :
            machR2[snp] = 1
        tmpR2 = ((dosageSqrSum / nonMissingCount) - math.pow((dosageSum / nonMissingCount), 2)) / (2 * estimatedAlleleFrequency * (1 - estimatedAlleleFrequency))
        machR2[snp] = 1 if tmpR2 > 1.0 else tmpR2
    selection = machR2 < min_hmachR2

    failed_snp_names.extend(list(snp_df_c.columns[selection]))
    snp_df_c = snp_df_c.loc[:,list(snp_df_c.columns[~selection])]

    return snp_df_c.columns, failed_snp_names


text:
# Copyright (c) 2020 NVIDIA Corporation. All rights reserved.
# This work is licensed under the NVIDIA Source Code License - Non-commercial. Full
# text can be found in LICENSE.md

from ycb_render.ycb_renderer import *
from ycb_render.tless_renderer_tensor import *
from utils.poserbpf_utils import *
import numpy.ma as ma
import matplotlib.pyplot as plt
import time
import numpy

numpy.random.seed(0)

class render_wrapper:
    def __init__(self, models, intrinsics, gpu_id, model_dir, model_ctg='ycb', im_w=640, im_h=480, initialize_render=True):

        obj_paths = []
        texture_paths = []
        assert model_ctg in ['ycb', 'tless'], "model {} is not supported".format(model_ctg)
        for item in models:
            if model_ctg == 'ycb':
                if os.path.exists('{}/ycb_models/{}/texture_map.png'.format(model_dir, item)):
                    obj_paths.append('{}/ycb_models/{}/textured_simple.obj'.format(model_dir, item))
                    texture_paths.append('{}/ycb_models/{}/texture_map.png'.format(model_dir, item))
                else:
                    obj_paths.append('{}/ycb_models/{}/textured_simple.ply'.format(model_dir, item))
                    texture_paths.append(''.format(model_dir, item))
            elif model_ctg == 'tless':
                obj_paths.append('{}/tless_models/{}.ply'.format(model_dir, item))
                texture_paths.append(''.format(model_dir, item))

        self.obj_paths = obj_paths
        self.texture_paths = texture_paths

        if initialize_render:
            if model_ctg == 'ycb':
                self.renderer = YCBRenderer(im_w, im_h, gpu_id=gpu_id)
            else:
                self.renderer = TLessTensorRenderer(im_w, im_h, gpu_id=gpu_id)
            self.renderer.load_objects(obj_paths, texture_paths)
            self.renderer.set_light_pos(cfg.TRAIN.TARGET_LIGHT1_POS)
            self.renderer.set_light_color([1.0, 1.0, 1.0])
            self.renderer_cam_pos = [0, 0, 0]
            self.renderer.set_camera_default()
            self.renderer.set_projection_matrix(im_w, im_h,
                                                intrinsics[0, 0], intrinsics[1, 1],
                                                intrinsics[0, 2], intrinsics[1, 2], 0.01, 10)
        else:
            self.renderer = None
        self.intrinsics = intrinsics
        self.im_w = im_w
        self.im_h = im_h

        # predefined parameters
        self.ymap_full = np.array([[j for i in range(int(self.im_w))] for j in range(self.im_h)])
        self.xmap_full = np.array([[i for i in range(int(self.im_w))] for j in range(self.im_h)])

        self.ymap_full_torch = torch.from_numpy(self.ymap_full).float().unsqueeze(2)
        self.xmap_full_torch = torch.from_numpy(self.xmap_full).float().unsqueeze(2)

        self.xmap_full_torch = self.xmap_full_torch.cuda()
        self.ymap_full_torch = self.ymap_full_torch.cuda()

        class_file = './datasets/ycb_video_classes.txt'
        if model_ctg == 'tless':
            class_file = './datasets/tless_classes.txt'

        with open(class_file, 'r') as class_name_file:
            self.class_names_all = class_name_file.read().split('\n')
        self.model_idx = []
        for model in models:
            self.model_idx.append(self.class_names_all.index(model) + 1)

    def set_intrinsics(self, intrinsics, im_w=640, im_h=480):
        self.renderer.set_camera_default()
        self.renderer.set_projection_matrix(im_w, im_h,
                                            intrinsics[0, 0], intrinsics[1, 1],
                                            intrinsics[0, 2], intrinsics[1, 2], 0.01, 10)
        self.intrinsics = intrinsics
        self.im_w = im_w
        self.im_h = im_h

    def render_pose(self, intrinsics, t, R, cls_id): # cls_id: target object index in the loaded classes
        self.set_intrinsics(intrinsics, im_w=self.im_w, im_h=self.im_h)
        t_v = t
        q_v = mat2quat(R)
        pose_v = np.zeros((7,))
        pose_v[:3] = t_v
        pose_v[3:] = q_v
        frame_cuda = torch.cuda.FloatTensor(int(self.im_h), int(self.im_w), 4)
        seg_cuda = torch.cuda.FloatTensor(int(self.im_h), int(self.im_w), 4)
        pc_cuda = torch.cuda.FloatTensor(int(self.im_h), int(self.im_w), 4)
        self.renderer.set_poses([pose_v])
        self.renderer.render([cls_id], frame_cuda, seg_cuda, pc2_tensor=pc_cuda)
        frame_cuda = frame_cuda.flip(0)
        frame_cuda = frame_cuda[:, :, :3].float()
        frame_cuda = frame_cuda.permute(2, 0, 1).unsqueeze(0)
        pc_cuda = pc_cuda.flip(0)
        depth_render = pc_cuda[:, :, 2]
        return frame_cuda, depth_render

    def render_pose_multiple(self, intrinsics, Tco_list, cls_id_list):
        self.set_intrinsics(intrinsics)

        # list of poses
        poses = []
        pose_v = np.zeros((7,))
        for Tco in Tco_list:
            t_v = Tco[:3, 3]
            q_v = mat2quat(Tco[:3, :3])
            pose_v[:3] = t_v
            pose_v[3:] = q_v
            poses.append(pose_v.copy())

        # render
        frame_cuda = torch.cuda.FloatTensor(int(self.im_h), int(self.im_w), 4)
        seg_cuda = torch.cuda.FloatTensor(int(self.im_h), int(self.im_w), 4)
        pc_cuda = torch.cuda.FloatTensor(int(self.im_h), int(self.im_w), 4)
        self.renderer.set_poses(poses)
        self.renderer.render(cls_id_list, frame_cuda, seg_cuda, pc2_tensor=pc_cuda)

        frame_cuda = frame_cuda.flip(0)
        frame_cuda = frame_cuda[:, :, :3].float()
        frame_cuda = frame_cuda.permute(2, 0, 1).unsqueeze(0)
        pc_cuda = pc_cuda.flip(0)
        depth_render = pc_cuda[:, :, 2]

        seg_cuda = seg_cuda.flip(0)
        seg_render_np = seg_cuda[:, :, 0].cpu().numpy()
        seg_est_np = np.zeros_like(seg_render_np)
        for i in cls_id_list:
            select = np.isclose(seg_render_np, float(self.renderer.colors[i][0]))
            seg_est_np[select] = self.model_idx[i]

        return frame_cuda, depth_render, seg_est_np

    def estimate_visibile_points(self, T_co_est, cls_id, depth, render_dist, intrinsics,
                                 delta=0.02, mask_input=None, output_cuda=True):

        t_v = T_co_est[:3, 3]
        R = T_co_est[:3, :3]

        # pose
        pose_v = np.zeros((7,))
        pose_v[:3] = t_v
        pose_v[3:] = mat2quat(R)

        # render
        frame_cuda = torch.cuda.FloatTensor(self.im_h, self.im_w, 4)
        seg_cuda = torch.cuda.FloatTensor(self.im_h, self.im_w, 4)
        pc_cuda = torch.cuda.FloatTensor(self.im_h, self.im_w, 4)

        self.renderer.set_poses([pose_v])
        self.renderer.render([cls_id], frame_cuda, seg_cuda, pc2_tensor=pc_cuda)

        uv = project(t_v, intrinsics)
        uv = np.expand_dims(uv, axis=0)
        z = np.array([[t_v[2]]])

        render_roi_cuda, _ = get_rois_cuda(pc_cuda.flip(0), uv, z,
                                                 intrinsics[0, 0], self.intrinsics[1, 1],
                                                 render_dist)
        if output_cuda:
            pc_cuda = render_roi_cuda[0, :3, :, :].permute(1, 2, 0)
        else:
            pc_cuda = render_roi_cuda[0, :3, :, :].cpu().permute(1, 2, 0)

        # get point cloud from depth
        pt2 = depth
        if output_cuda:
            pt2 = pt2.cuda()
        pt0 = (self.xmap_full_torch - intrinsics[0, 2]) * pt2 / intrinsics[0, 0]
        pt1 = (self.ymap_full_torch - intrinsics[1, 2]) * pt2 / intrinsics[1, 1]
        pc_c_torch = torch.cat((pt0, pt1, pt2), dim=2)

        if mask_input is not None:
            if np.sum(mask_input.numpy() == self.model_idx[cls_id]) == 0:
                mask_input = None

        if mask_input is None:
            pc_roi_cuda, _ = get_rois_cuda(pc_c_torch, uv, z, intrinsics[0, 0], intrinsics[1, 1],
                                                    render_dist)
            # compare the depth
            depth_meas_roi = pc_roi_cuda[0, 2, :, :].cpu().numpy()
            depth_render_roi = pc_cuda[:, :, 2].cpu().numpy()
            mask_depth_meas = ma.getmaskarray(ma.masked_not_equal(depth_meas_roi, 0))
            mask_depth_render = ma.getmaskarray(ma.masked_greater(depth_render_roi, 0))
            mask_depth_vis = ma.getmaskarray(ma.masked_less(np.abs(depth_render_roi - depth_meas_roi), delta))
            mask = mask_depth_meas * mask_depth_render * mask_depth_vis
            mask_disp = mask.copy()
            choose = mask.flatten().nonzero()[0]

            if np.sum(mask_depth_render) == 0:
                visible_ratio = 0
                ps_c = None
                mask_disp = None
                return ps_c, visible_ratio, mask_disp

            visible_ratio = np.sum(mask) * 1.0 / (np.sum(mask_depth_render) * 1.0)
            pt2_valid = depth_meas_roi.flatten()[choose][:, np.newaxis].astype(np.float32)
            pt0_valid = pc_roi_cuda[0, 0, :, :].cpu().numpy().flatten()[choose][:, np.newaxis].astype(np.float32)
            pt1_valid = pc_roi_cuda[0, 1, :, :].cpu().numpy().flatten()[choose][:, np.newaxis].astype(np.float32)
        else:
            pc_roi_cuda, _ = get_rois_cuda(pc_c_torch, uv, z, intrinsics[0, 0], intrinsics[1, 1],
                                                 render_dist)
            # compare the depth
            depth_meas_roi = pc_roi_cuda[0, 2, :, :].cpu().numpy()
            depth_render_roi = pc_cuda[:, :, 2].cpu().numpy()
            mask_depth_meas = ma.getmaskarray(ma.masked_not_equal(depth_meas_roi, 0))
            mask_depth_render = ma.getmaskarray(ma.masked_greater(depth_render_roi, 0))
            mask_depth_vis = ma.getmaskarray(ma.masked_less(np.abs(depth_render_roi - depth_meas_roi), delta))
            mask_prbpf = mask_depth_meas * mask_depth_render * mask_depth_vis
            visible_ratio_prbpf = np.sum(mask_prbpf) * 1.0 / (np.sum(mask_depth_render) * 1.0)

            mask_roi_cuda, _ = get_rois_cuda(mask_input, uv, z, intrinsics[0, 0], intrinsics[1, 1],
                                                   render_dist)
            mask_roi = mask_roi_cuda[0, 0, :, :].cpu().numpy()
            mask_input = ma.getmaskarray(ma.masked_equal(mask_roi, self.model_idx[cls_id]))
            mask_depth_vis2 = ma.getmaskarray(ma.masked_less(np.abs(depth_render_roi - depth_meas_roi), delta * 2.5))
            mask_fuse = mask_depth_meas * mask_input * mask_depth_vis2 * mask_depth_render
            visible_ratio_fuse = np.sum(mask_fuse) * 1.0 / (np.sum(mask_depth_render) * 1.0)

            if visible_ratio_fuse > 1.2 * visible_ratio_prbpf:
                choose = mask_fuse.flatten().nonzero()[0]
                visible_ratio = visible_ratio_fuse
                mask_disp = mask_fuse.copy()
            else:
                choose = mask_prbpf.flatten().nonzero()[0]
                visible_ratio = visible_ratio_prbpf
                mask_disp = mask_prbpf.copy()

            pt2_valid = depth_meas_roi.flatten()[choose][:, np.newaxis].astype(np.float32)
            pt0_valid = pc_roi_cuda[0, 0, :, :].cpu().numpy().flatten()[choose][:, np.newaxis].astype(np.float32)
            pt1_valid = pc_roi_cuda[0, 1, :, :].cpu().numpy().flatten()[choose][:, np.newaxis].astype(np.float32)

        ps_c = np.concatenate((pt0_valid, pt1_valid, pt2_valid, np.ones_like(pt0_valid)), axis=1)
        ps_c = torch.from_numpy(ps_c)
        if output_cuda:
            ps_c = ps_c.cuda()

        return ps_c, visible_ratio, mask_disp

    def estimate_visibile_points_multi(self, Tco_list, cls_id_list, depth, intrinsics, mask_input,
                                        delta=0.02, output_cuda=True):

        # render with current pose
        self.set_intrinsics(intrinsics)

        # list of poses
        poses = []
        pose_v = np.zeros((7,))
        for Tco in Tco_list:
            t_v = Tco[:3, 3]
            q_v = mat2quat(Tco[:3, :3])
            pose_v[:3] = t_v
            pose_v[3:] = q_v
            poses.append(pose_v.copy())

        # render
        frame_cuda = torch.cuda.FloatTensor(int(self.im_h), int(self.im_w), 4)
        seg_cuda = torch.cuda.FloatTensor(int(self.im_h), int(self.im_w), 4)
        pc_cuda = torch.cuda.FloatTensor(int(self.im_h), int(self.im_w), 4)
        self.renderer.set_poses(poses)
        self.renderer.render(cls_id_list, frame_cuda, seg_cuda, pc2_tensor=pc_cuda)

        frame_cuda = frame_cuda.flip(0)
        frame_cuda = frame_cuda[:, :, :3].float()
        frame_cuda = frame_cuda.permute(2, 0, 1).unsqueeze(0)
        pc_cuda = pc_cuda.flip(0)
        depth_render_np = pc_cuda[:, :, 2].cpu().numpy()
        seg_cuda = seg_cuda.flip(0)
        seg_render_np = seg_cuda[:, :, :3].cpu().numpy()

        seg_input_np = mask_input[:, :, 0].cpu().numpy()
        seg_est_np = np.zeros_like(seg_input_np)
        for i in cls_id_list:
            select = np.isclose(seg_render_np[:, :, 0], float(self.renderer.colors[i][0]))
            seg_est_np[select] = self.model_idx[i]

        # convert depth to point cloud
        pt2 = depth
        if output_cuda:
            pt2 = pt2.cuda()
        pt0 = (self.xmap_full_torch - intrinsics[0, 2]) * pt2 / intrinsics[0, 0]
        pt1 = (self.ymap_full_torch - intrinsics[1, 2]) * pt2 / intrinsics[1, 1]

        depth_meas_np = depth[:, :, 0].cpu().numpy()
        mask_depth_meas = ma.getmaskarray(ma.masked_not_equal(depth_meas_np, 0))

        if not mask_input is None:
            delta = delta * 2
        mask_depth_vis = ma.getmaskarray(ma.masked_less(np.abs(depth_render_np - depth_meas_np), delta))

        pt2_np = pt2.cpu().numpy()
        pt0_np = pt0.cpu().numpy()
        pt1_np = pt1.cpu().numpy()

        pc_list = []

        for i in cls_id_list:
            mask_seg_est = ma.getmaskarray(ma.masked_equal(seg_est_np, self.model_idx[i]))
            mask_seg_meas = ma.getmaskarray(ma.masked_equal(seg_input_np, self.model_idx[i]))

            mask_out_merge = mask_depth_meas * mask_seg_est * mask_seg_meas * mask_depth_meas * mask_depth_vis

            choose = mask_out_merge.flatten().nonzero()[0]

            pt2_valid = pt2_np.flatten()[choose][:, np.newaxis].astype(np.float32)
            pt0_valid = pt0_np.flatten()[choose][:, np.newaxis].astype(np.float32)
            pt1_valid = pt1_np.flatten()[choose][:, np.newaxis].astype(np.float32)

            ps_c = np.concatenate((pt0_valid, pt1_valid, pt2_valid, np.ones_like(pt0_valid)), axis=1)
            ps_c = torch.from_numpy(ps_c).clone()
            if output_cuda:
                ps_c = ps_c.cuda()

            pc_list.append(ps_c)

        return pc_list

    # evaluate particles according to depth measurements
    def evaluate_depths_init(self, cls_id, depth, uv, z, q_idx, intrinsics, render_dist, codepose, delta=0.03, tau=0.05, mask=None):

        score = np.zeros_like(z)

        # crop rois
        depth_roi_cuda, _ = get_rois_cuda(depth.detach(), uv, z, intrinsics[0, 0], intrinsics[1, 1],
                                                 render_dist)
        depth_roi_np = depth_roi_cuda.cpu().numpy()

        if mask is not None:
            mask_roi_cuda, _ = get_rois_cuda(mask.detach(), uv, z, intrinsics[0, 0], intrinsics[1, 1],
                                                    render_dist)
            mask_roi_np = mask_roi_cuda.cpu().numpy()

        # render
        pose_v = np.zeros((7,))
        frame_cuda = torch.cuda.FloatTensor(self.im_h, self.im_w, 4)
        seg_cuda = torch.cuda.FloatTensor(self.im_h, self.im_w, 4)
        pc_cuda = torch.cuda.FloatTensor(self.im_h, self.im_w, 4)

        q_idx_unique, idx_inv = np.unique(q_idx, return_inverse=True)
        pc_render_all = np.zeros((q_idx_unique.shape[0], 128, 128, 3), dtype=np.float32)
        q_render = codepose[q_idx_unique][:, 3:]
        for i in range(q_render.shape[0]):
            pose_v[:3] = [0, 0, render_dist]
            pose_v[3:] = q_render[i]
            self.renderer.set_poses([pose_v])
            self.renderer.render([cls_id], frame_cuda, seg_cuda, pc2_tensor=pc_cuda)
            render_roi_cuda, _ = get_rois_cuda(pc_cuda.flip(0),
                                                     np.array([[intrinsics[0, 2], intrinsics[1, 2], 1]]),
                                                     np.array([[render_dist]]),
                                                     intrinsics[0, 0], intrinsics[1, 1],
                                                     render_dist)
            pc_render_all[i] = render_roi_cuda[0, :3, :, :].permute(1, 2, 0).cpu().numpy()

        # evaluate every particle
        for i in range(uv.shape[0]):

            pc_render = pc_render_all[idx_inv[i]].copy()
            depth_mask = pc_render[:, :, 2] > 0
            pc_render[:, :, 2][depth_mask] = pc_render[:, :, 2][depth_mask] - render_dist + z[i]

            depth_render = pc_render[:, :, [2]]

            depth_meas_np = depth_roi_np[i, 0, :, :]
            depth_render_np = depth_render[:, :, 0]

            # compute visibility mask
            if mask is None:
                visibility_mask = estimate_visib_mask_numba(depth_meas_np, depth_render_np, delta=delta)
            else:
                visibility_mask = np.logical_and((mask_roi_np[i, 0, :, :] > 0), np.logical_and(depth_meas_np>0, depth_render_np>0))

            if np.sum(visibility_mask) == 0:
                continue

            # compute depth error
            depth_error = np.abs(depth_meas_np[visibility_mask] - depth_render_np[visibility_mask])
            depth_error /= tau
            depth_error[depth_error > 1] = 1

            # score computation
            total_pixels = np.sum((depth_render_np > 0).astype(np.float32))
            if total_pixels is not 0:
                vis_ratio = np.sum(visibility_mask.astype(np.float32)) / total_pixels
                score[i] = (1 - np.mean(depth_error)) * vis_ratio
            else:
                score[i] = 0

        return score

    # evaluate particles according to depth measurements
    def evaluate_depths_tracking(self, rbpf, cls_id, depth, uv, z, q_idx, intrinsics, render_dist, codepose, rbpf_ready,
                        delta=0.03, tau=0.05, mask=None):

        if mask is not None:
            depth = depth.float() * mask

        # crop rois
        depth_roi_cuda, _ = get_rois_cuda(depth.detach(), uv, z, intrinsics[0, 0], intrinsics[1, 1],
                                                 render_dist)
        depth_meas_all = depth_roi_cuda[:, 0, :, :]

        # render
        pose_v = np.zeros((7,))
        frame_cuda = torch.cuda.FloatTensor(self.im_h, self.im_w, 4)
        seg_cuda = torch.cuda.FloatTensor(self.im_h, self.im_w, 4)
        pc_cuda = torch.cuda.FloatTensor(self.im_h, self.im_w, 4)

        fast_rendering = False
        if rbpf_ready and np.linalg.norm(rbpf.trans_bar) > 0: # fast rendering for tracking
            pose_v[:3] = rbpf.trans_bar
            pose_v[3:] = mat2quat(rbpf.rot_bar)
            self.renderer.set_poses([pose_v])
            self.renderer.render([cls_id], frame_cuda, seg_cuda, pc2_tensor=pc_cuda)
            uv_crop = project(rbpf.trans_bar, intrinsics)
            uv_crop = np.repeat(np.expand_dims(uv_crop, axis=0), uv.shape[0], axis=0)
            z_crop = np.ones_like(z) * rbpf.trans_bar[2]
            render_roi_cuda, _ = get_rois_cuda(pc_cuda.flip(0), uv_crop, z_crop,
                                                      intrinsics[0, 0],
                                                      intrinsics[1, 1],
                                                      render_dist)
            depth_render_all = render_roi_cuda[:, 2, :, :]
            fast_rendering = True
        else:
            q_idx_unique, idx_inv = np.unique(q_idx, return_inverse=True)
            depth_render_uni = torch.zeros((q_idx_unique.shape[0], 128, 128), dtype=torch.float32).cuda()
            q_render = codepose[q_idx_unique][:, 3:]
            for i in range(q_render.shape[0]):
                pose_v[:3] = [0, 0, render_dist]
                pose_v[3:] = q_render[i]
                self.renderer.set_poses([pose_v])
                self.renderer.render([cls_id], frame_cuda, seg_cuda, pc_cuda)
                render_roi_cuda, _ = get_rois_cuda(pc_cuda.flip(0),
                                                         np.array([[intrinsics[0, 2], intrinsics[1, 2], 1]]),
                                                         np.array([[render_dist]]),
                                                         intrinsics[0, 0], intrinsics[1, 1],
                                                         render_dist)
                depth_render_uni[i] = render_roi_cuda[0, 2, :, :].clone()
            depth_render_all = torch.zeros((q_idx.shape[0], 128, 128), dtype=torch.float32).cuda()
            for i in range(q_idx.shape[0]):
                depth_render_all[i] = depth_render_uni[idx_inv[i]].clone()

        # eval particles
        # shift the rendered image
        depth_invalid_mask = depth_render_all == 0
        if fast_rendering:
            depth_shift_np = z[:, 0] - rbpf.trans_bar[2]
        else:
            depth_shift_np = z[:, 0] - render_dist
        depth_shift = torch.from_numpy(depth_shift_np).cuda().float().repeat(depth_render_all.size(1),
                                                                             depth_render_all.size(2),
                                                                             1).permute(2, 0, 1)
        depth_render_all += depth_shift
        depth_render_all[depth_invalid_mask] = 0

        # compute visibility mask
        visibility_mask_cuda = estimate_visib_mask_cuda(depth_meas_all, depth_render_all, delta=delta)

        # compute scores
        # depth errors
        depth_error = torch.ones_like(depth_render_all).cuda()
        depth_error[visibility_mask_cuda] = torch.abs(depth_meas_all[visibility_mask_cuda] -
                                                      depth_render_all[visibility_mask_cuda]) / tau
        depth_error = torch.clamp(depth_error, 0, 1.0)
        depth_error_mean = torch.mean(depth_error, (2, 1))

        # visible ratio
        total_pixels = torch.sum(depth_render_all > 0, (2, 1)).float()
        total_pixels[total_pixels == 0] = 10000
        vis_ratio = torch.sum(visibility_mask_cuda, (2, 1)).float() / total_pixels

        # scores
        score = (torch.ones_like(depth_error_mean) - depth_error_mean) * vis_ratio

        return score.unsqueeze(1), vis_ratio

    def transform_points(self, Transformation, points):
        return torch.matmul(Transformation, points.permute(1, 0)).permute(1, 0)

max_stars_repo_name: xEnVrE/PoseRBPF

max_issues_repo_name: xEnVrE/PoseRBPF

max_forks_repo_name: xEnVrE/PoseRBPF


text:
from typing import Optional, Union

import numpy as np

from ._math import euclidean_distance, manhattan_distance
from ..utils import check_is_fitted

class KNeighborsRegressor:
    """
    Performing regression using k-nearest neighbors model.

    k-NN is a non-parametric ML supervised learning model, thus unlike logistic regressions or neural network
    this model does NOT have a parameters to trained (i.e., weights and biases).

    Parameters
    ----------
    k:  int. The default is 2
        The number of nearest neighbors.

    p:  int. The default is 2.
        The method to calculate the data points distance
        p = 1 for l1 distance or manhattan distance
        p = 2 for l2 distance or euclidean distance

    Attributes
    ------

    """
    def __init__(self, k: int = 2, p: int = 2):
        min_k = 1
        if type(k) is not int or k <= min_k:
            raise ValueError(f"The number of nearest neighbors (k) has to be an integer "
                             f"that is larger than {min_k}! The input was k={k} with {type(k)} type.")
        if p not in [1, 2]:
            raise ValueError(f"The current model only supports either the l1, manhattan, distance (p=1) "
                             f"or l2, euclidean, distance (p=2)! The input was p={p}.")

        self.k = k
        self.p = p
        self.calc_distance = manhattan_distance if p == 1 else euclidean_distance

    def fit(self, X: np.array, y: np.array):
        """
        Fitting the K-NN regression model to the training dataset

        Parameters
        ----------
        X:  ndarray with (n_train_samples, n_features) shape. Training data with,
            n_train_samples number of samples and n_features number of features.
        y:  ndarray with (n_train_samples,) shape. The training data label/target.

        Return
        ------
        self:   KNeighborsRegressor
                The fitted k-NN regression model, based on X and y.
        """
        # Storing the training dataset
        self.X_ = X
        self.y_ = y

        # k-NN is a non-parametric model, it does NOT have parameters (i.e., weights and biases) to train!
        return self

    def predict(self, X: np.array) -> np.array:
        """
        Predict the target/label of the testing dataset,
        by performing the k-NN regression fitted by the training dataset.

        Parameters
        ----------
        X:  ndarray with (n_test_samples, n_features) shape. Testing data with,
            n_test_samples number of samples and n_features number of features.
            The target/label of this dataset will be predicted later.

        Return
        ------
        y_pred: ndarray with (n_test_samples,) shape.
                The predicted target/label of the testing set.
        """
        check_is_fitted(self)  # Check if self.fit is already performed before
        y_pred_list = []  # Initialize the prediction

        for i in range(len(X)):  # Using loop to avoid memory overflow if the array is too large
            # Computing the euclidean distance of each feature (X_i) to the testing set (X_ij)
            distance = self.calc_distance(self.X_, X[i])
            distance_sorted_id = np.argsort(distance)

            # Predicting the new label (y_pred) based on the average of the k-nearest y value
            y_pred_i = self.y_[distance_sorted_id][:self.k].mean()
            y_pred_list.append(y_pred_i)

        y_pred = np.array(y_pred_list)
        return y_pred

    def fit_predict(self, X_train, y_train, X_test) -> np.array:
        """
        Sequentially predict the target/label of the testing data after fitting the model
        to the training dataset.

        Parameters
        ----------
        X_train:    ndarray with (n_train_samples, n_features) shape. Training data with,
                    n_train_samples number of samples and n_features number of features.

        y_train:    ndarray with (n_train_samples,) shape. The training data label/target.

        X_test:     ndarray with (n_test_samples, n_features) shape. Testing data with,
                    n_test_samples number of samples and n_features number of features.
                    The target/label of this dataset will be predicted later.

        Return
        ------
        y_pred:     ndarray with (n_test_samples,) shape.
                    The predicted target/label of the testing set.
        """
        y_pred = self.fit(X_train, y_train).predict(X_test)
        return y_pred

max_stars_repo_name: abrosua/ml-with-numpy

max_issues_repo_name: abrosua/ml-with-numpy

max_forks_repo_name: abrosua/ml-with-numpy


text:
###############
#
# Transform R to Python Copyright (c) 2019 Masahiro Imai Released under the MIT license
#
###############

import os

import numpy as np
import pystan
import pandas
import pickle
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn.preprocessing import LabelEncoder

fish_num_climate_4 = pandas.read_csv('4-3-1-fish-num-4.csv')
print(fish_num_climate_4.head())
print(fish_num_climate_4.describe())

sns.scatterplot(
    x='temperature',
    y='fish_num',
    hue='human',
    data=fish_num_climate_4
)
plt.show()

fish_num_climate_4_d = pandas.get_dummies(fish_num_climate_4, columns=["human"])
print(fish_num_climate_4_d.head())

fish_num = fish_num_climate_4_d['fish_num']
sample_num = len(fish_num)
temperature = fish_num_climate_4_d['temperature']

# creating teamID
le = LabelEncoder()
le = le.fit(fish_num_climate_4['human'])
fish_num_climate_4['human'] = le.transform(fish_num_climate_4['human'])

sns.scatterplot(
    x='temperature',
    y='fish_num',
    hue='human',
    legend="full",
    data=fish_num_climate_4
)
plt.show()

human_id = fish_num_climate_4['human'].values
human_num = len(np.unique(human_id))
human_id = (human_id + 1)
print(human_id)

stan_data = {
    'N': sample_num,
    'fish_num': fish_num,
    'temp': temperature,
    'human_id': human_id
}

if os.path.exists('4-3-1-poisson-glmm.pkl'):
    sm = pickle.load(open('4-3-1-poisson-glmm.pkl', 'rb'))
    # sm = pystan.StanModel(file='4-3-1-poisson-glmm.stan')
else:
    # a model using prior for mu and sigma.
    sm = pystan.StanModel(file='4-3-1-poisson-glmm.stan')

mcmc_result = sm.sampling(
    data=stan_data,
    seed=1,
    chains=4,
    iter=2000,
    warmup=1000,
    thin=1
)

print(mcmc_result)
mcmc_result.plot()
plt.show()

# saving compiled model
if not os.path.exists('4-3-1-poisson-glmm.pkl'):
    with open('4-3-1-poisson-glmm.pkl', 'wb') as f:
        pickle.dump(sm, f)

mcmc_sample = mcmc_result.extract()
print(mcmc_sample)

# visualization
label_temp = np.arange(10,21)
df = pandas.DataFrame(mcmc_sample, columns=['Intercept', 'b_temp', 'b_human', 'b_temp_human'])

for h in np.arange(human_num):
    y = []
    for i in label_temp:
        y.append(np.exp(df['Intercept'].mean() +
                        df['b_temp'].mean() * i +
                        df['b_human'].mean() * (h+1) +
                        df['b_temp_human'].mean() * (h+1) * i))

    plt.plot(label_temp, np.array(y), 'red', label='%d' % (h + 1))
    plt.scatter(fish_num_climate_4.query('human == %d' % h)["temperature"],
                fish_num_climate_4.query('human == %d' % h)["fish_num"],
                c='r')
    plt.legend()
    plt.show()

max_stars_repo_name: MasazI/python-r-stan-bayesian-model

max_issues_repo_name: MasazI/python-r-stan-bayesian-model

max_forks_repo_name: MasazI/python-r-stan-bayesian-model


text:
# ----------------------------------------------
# -  Input deck for interpolation performance  -
# ----------------------------------------------

import numpy as np
from input.utils import *
from input.misc import *

####################### Preprocessing ##########################

run_name = "test_performance_remap"

## electron and trajectory
gamma=1000
# Ct=2
delta_alpha_B = 0.2  # assume the bending section spreads over \Delta\alpha=0.2
lbeam = 3000        # um
dbeam = 50          # um
delta_alpha_tilde = 10000 #gamma * delta_alpha_B / 2
delta_chi_tilde = dbeam * 1e-6 * gamma**2.0

if (mpi_rank==0) : print("delta_alpha_tilde={}, delta_chi_tilde={}".format(delta_alpha_tilde, delta_chi_tilde))

mesh_alpha_npt = 4000 #int(delta_alpha_tilde/1.0+1.0)
mesh_chi_npt = 4000 #int(delta_chi_tilde/1.0+1.0)

if not mesh_alpha_npt%2:
    mesh_alpha_npt += 1
if not mesh_chi_npt%2:
    mesh_chi_npt += 1

if (mpi_rank==0) : print("mesh_alpha_npt={}, mesh_chi_npt={}".format(mesh_alpha_npt, mesh_chi_npt))

Nw = 2000

####################### Main setup ##########################

# wavelet emission
num_wavefronts = 200               # number of wavefronts
num_dirs = 0                       # number of field lines
num_step = num_wavefronts          # number of steps
dt = delta_alpha_B/num_wavefronts  # time step in electron rest frame
emission_interval = num_step - 1   # only emit wavefronts at simulation end (test purpose)
# dt=Ct/(gamma*Nw)

## remap
remap_interval = num_step - 1      # interval of doing remapping (in time steps)
remap_scatter = False              # use scatter weights form for remap
remap_adaptive = False             # use adaptive smoothing length for remap
remap_scaling[0] = 0.2             # support/smoothing length scaling factor
remap_scaling[1] = 0.2             # support/smoothing length scaling factor
remap_verbose = False              # print remap statistics

# --- eletron beam info ---
beam_charge = 0.01                 # nC
num_particles = 10                 # number of particles
trajectory_type = 2                # 1: straight line, 2: circular, 3: sinusoidal
parameters[0] = gamma              # parameters for describing the motion of electron
parameters[1] = 100.0

beam = init_beam(num_particles, gamma, lbeam, dbeam, mpi_rank)
del beam                           # init a single particle instead

# --- comoving mesh info ---
num_gridpt_hor = mesh_alpha_npt                # number of points in x-axis
num_gridpt_ver = mesh_chi_npt                  # number of points in y-axis
mesh_span_angle = delta_alpha_tilde/gamma**3   # in radians
mesh_width = delta_chi_tilde/gamma**2          # in unit of radius

# --- wavelets ---
#------- generate regular mesh points and used as wavelets
x = np.linspace(-delta_alpha_tilde/2.0, delta_alpha_tilde/2.0, np.int(mesh_alpha_npt*1))
y = np.linspace(-delta_chi_tilde/2.0, delta_chi_tilde/2.0, np.int(mesh_chi_npt*1))
ymax = y.max()
ymin = y.min()
wavelet_x, wavelet_y = np.meshgrid(x,y)
wavelet_x = wavelet_x.flatten()
wavelet_y = wavelet_y.flatten()
wavelet_field = np.cos(5.0*wavelet_y*2*np.pi/(ymax-ymin))
num_wavelet_fields = 1

# unscale if necessary
wavelet_x /= gamma**3.0
wavelet_y /= gamma**2.0
wavelet_field *= gamma**4.0

min_emit_angle = 10.0

wavelet_type = 1   # 0: use global (x,y) coordinate; 1: use local (x',y') coordinate; 2: (TODO) use local cylindrical coordinate
if (wavelet_type == 0) :
   rotation_angle = 0.1
   wavelet_x, wavelet_y = convert2global(wavelet_x, wavelet_y, rotation_angle)

use_wavelet_for_subcycle = False # True: loaded wavelets will be repeatedly emitted at each step and copied into internal wavelets array, otherwise only used when interpolation is done and not copied into internal wavelets array
num_wavelet_fields = 1

####################### Diagnostics ##########################

print_interval = 100        # interval for printing simulation steps
beam_output_interval = num_step - 1
mesh_output_interval = num_step - 1
wavelet_output_interval = num_step - 1
beam_output = False

if (mpi_rank==0):
  make_output_dirs(run_name+"/beam", num_step, beam_output_start, beam_output_interval)
  mesh_output = True
  make_output_dirs(run_name+"/mesh", num_step, mesh_output_start, mesh_output_interval)
  wavelet_output = False
  make_output_dirs(run_name+"/wavelet", num_step, wavelet_output_start, wavelet_output_interval)
  make_output_dirs(run_name+"/traj", num_step, wavelet_output_start, wavelet_output_interval)

max_stars_repo_name: jsdomine/cosyr

max_issues_repo_name: jsdomine/cosyr

max_forks_repo_name: jsdomine/cosyr


text:
import os
import sys
import numpy as np
import pickle

datafile = sys.argv[1]
dir = os.path.dirname(datafile)

embs_bootleg = {'ent':{'inname':'ctx_embeddings.npy', 'outname':'ent_embedding.npy', 'outvocab':'ent_vocab.pkl'}}

for k, names_lst in embs_bootleg.items():
    datafile = dir+"/"+names_lst['inname']

    if os.path.isfile(datafile):
        emb_data = np.load(datafile)

        new_embs = np.zeros((emb_data.shape[0] + 2, emb_data.shape[1]), dtype='float')
        new_embs[2:] = emb_data
        np.save('{DIR}/{out}'.format(DIR=dir,out=names_lst['outname']), new_embs)

        v = [i-2 for i in range(new_embs.shape[0])]
        vocab_file = '{DIR}/{outvoc}'.format(DIR=dir,outvoc=names_lst['outvocab'])

        with open(vocab_file, 'wb') as outfile:
            pickle.dump(v, outfile)
        print("Saved {}".format(k))

max_stars_repo_name: simran-arora/tacred-relation

max_issues_repo_name: simran-arora/tacred-relation

max_forks_repo_name: simran-arora/tacred-relation


text:
import os
import h5py
import numpy as np
import logging

logger = logging.getLogger(__name__)

class HDF5BaseWriter:

    def __init__(self, dbPath, estNumSamples=500, maxBufferSize=50000, dbResizeFactor=2):
        """Class constructor.

        Arguments:
            dbPath {str} -- Path to where the features database will be stored.

        Keyword Arguments:
            estNumSamples {int} -- Approximate number of samples in the dataset. (default: {500})
            maxBufferSize {int} -- Maximum number of feature vectors to be stored in memory until the buffer is flushed to HDF5. (default: {50000})
            dbResizeFactor {int} -- Factor to resize the database when needed to expand it. (default: {2})
        """

        self.dbPath = dbPath
        self.estNumSamples = estNumSamples
        self.maxBufferSize = maxBufferSize
        self.dbResizeFactor = dbResizeFactor

        # Indexes dictionary to control datasets append
        self.idxs = {}

    def _writeBuffers(self):
        raise RuntimeError('Must be implemented by subclasses.')

    def _writeBuffer(self, dataset, datasetName, buf, sparse=False):
        """Write the data in the buffers to disk.

        Arguments:
            dataset {[type]} -- [description]
            datasetName {[type]} -- [description]
            buf {[type]} -- [description]

        Keyword Arguments:
            sparse {bool} -- [description] (default: {False})
        """

		# If the buffer is a list, then compute the ending index based on
		# the lists length
        if type(buf) is list:
            end = self.idxs[datasetName] + len(buf)

		# Otherwise, assume that the buffer is a NumPy/SciPy array, so
		# compute the ending index based on the array shape
        else:
            end = self.idxs[datasetName] + buf.shape[0]

		# Check to see if the dataset needs to be resized
        if end > dataset.shape[0]:
            logger.debug('Triggering `{}` db resize'.format(datasetName))
            self._resizeDataset(dataset, datasetName, baseSize=end)

		# if this is a sparse matrix, then convert the sparse matrix to a
		# dense one so it can be written to file
        if sparse:
            buf = buf.toarray()

		# dump the buffer to file
        logger.debug('Writing `{}` buffer'.format(datasetName))
        dataset[self.idxs[datasetName]:end] = buf

    def _resizeDataset(self, dataset, dbName, baseSize=0, finished=0):
        """[summary]

        Arguments:
            dataset {[type]} -- [description]
            dbName {[type]} -- [description]

        Keyword Arguments:
            baseSize {int} -- [description] (default: {0})
            finished {int} -- [description] (default: {0})
        """

		# grab the original size of the dataset
        origSize = dataset.shape[0]

		# check to see if we are finished writing rows to the dataset, and if
		# so, make the new size the current index
        if finished > 0:
            newSize = finished

		# otherwise, we are enlarging the dataset so calculate the new size
		# of the dataset
        else:
            newSize = baseSize * self.dbResizeFactor

		# determine the shape of (to be) the resized dataset
        shape = list(dataset.shape)
        shape[0] = newSize

		# show the old versus new size of the dataset
        dataset.resize(tuple(shape))
        logger.debug('old size of `{}`: {:,}; new size: {:,}'.format(dbName, origSize, newSize))

max_stars_repo_name: kleysonr/snsdl

max_issues_repo_name: kleysonr/snsdl

max_forks_repo_name: kleysonr/snsdl


text:
import numpy as np

from ..local_interpolation import ThirdOrderHermitePolynomialInterpolation
from .base import AbstractAdaptiveSDESolver, AbstractStratonovichSolver
from .runge_kutta import AbstractERK, ButcherTableau

_midpoint_tableau = ButcherTableau(
    a_lower=(np.array([0.5]),),
    b_sol=np.array([0.0, 1.0]),
    b_error=np.array([1.0, -1.0]),
    c=np.array([0.5]),
)

class Midpoint(AbstractERK, AbstractStratonovichSolver, AbstractAdaptiveSDESolver):
    """Midpoint method.

    2nd order explicit Runge--Kutta method. Has an embedded Euler method for adaptive
    step sizing.

    Also sometimes known as the "modified Euler method".

    When used to solve SDEs, converges to the Stratonovich solution.
    """

    tableau = _midpoint_tableau
    interpolation_cls = ThirdOrderHermitePolynomialInterpolation.from_k

    def order(self, terms):
        return 2

    def strong_order(self, terms):
        return 0.5

max_stars_repo_name: nicholasjng/diffrax

max_issues_repo_name: nicholasjng/diffrax

max_forks_repo_name: nicholasjng/diffrax


text:
import os
import numpy as np
import cv2
from torch.utils.data.dataset import Dataset
from utils.common_util import *
import pandas as pd
from config.config import *
from datasets.tool import *
from utils.augment_util import *
from tqdm import tqdm
import shutil

class ProteinMLDataset(Dataset):
    def __init__(self,
                 split_file,
                 img_size=512,
                 transform=None,
                 return_label=True,
                 is_trainset=True,
                 in_channels=3,
                 ):
        self.is_trainset = is_trainset
        self.img_size = img_size
        self.return_label = return_label
        self.in_channels = in_channels
        self.transform = transform
        data_type = 'train' if is_trainset else 'test'
        split_df = pd.read_csv(split_file)
        base_dir = DATA_DIR

        if EXTERNAL not in split_df.columns:
            split_df[EXTERNAL] = False

        self.is_trainset = is_trainset

        self.split_df = split_df
        self.is_external = self.split_df[EXTERNAL].values
        self.img_ids = self.split_df[ID].values

        if self.return_label:
            self.labels = self.split_df[ANTIBODY_CODE].values.astype(int)

        self.num = len(self.img_ids)

        self.set_img_dir(base_dir, data_type)

        if ANTIBODY in self.split_df.columns:
            print(self.split_df[ANTIBODY].nunique())

    def set_img_dir(self, base_dir, data_type):
        if self.img_size<=512:
            self.img_dir = opj(base_dir, data_type, 'images')
            self.external_img_dir = opj(base_dir, 'train', 'external_v18_512')
        elif self.img_size<=768:
            self.img_dir = opj(base_dir, data_type, 'images_768')
            self.external_img_dir = opj(base_dir, 'train', 'external_v18_768')
        else:
            self.img_dir = opj(base_dir, data_type, 'images_1536')
            self.external_img_dir = opj(base_dir, 'train', 'external_v18_1536')
        print(self.img_dir)
        print(self.external_img_dir)

    def read_rgby(self, img_dir, img_id, index):
        if self.is_external[index]:
            img_is_external = True
        else:
            img_is_external = False

        suffix = '.jpg' if img_is_external else '.png'
        if self.in_channels == 3:
            colors = ['red', 'green', 'blue']
        else:
            colors = ['red', 'green', 'blue', 'yellow']

        flags = cv2.IMREAD_GRAYSCALE
        img = [cv2.imread(opj(img_dir, img_id + '_' + color + suffix), flags)
               for color in colors]
        img = np.stack(img, axis=-1)

        img = cv2.resize(img, (self.img_size, self.img_size), interpolation=cv2.INTER_LINEAR)
        img = img / 255.0
        return img

    def __getitem__(self, index):
        img_id = self.img_ids[index]
        if self.is_external[index]:
            img_dir = self.external_img_dir
        else:
            img_dir = self.img_dir

        image = self.read_rgby(img_dir, img_id, index)
        if image is None:
            print(img_dir, img_id)

        if self.transform is not None:
            image = self.transform(image)

        image = image_to_tensor(image)
        if self.return_label:
            label = self.labels[index]
            return image, label, index
        else:
            return image, index

    def __len__(self):
        return self.num

max_stars_repo_name: guitarmind/HPA-competition-solutions

max_issues_repo_name: guitarmind/HPA-competition-solutions

max_forks_repo_name: guitarmind/HPA-competition-solutions


text:
import numpy as np
from keras.layers import LSTM, Input
from keras.models import Model
from keras.utils.np_utils import to_categorical

from generate_data import generate_sort_data, generate_x_y_for_inference
from PointerLSTM import PointerLSTM
import argparse
from sklearn.metrics import accuracy_score

# Instantiate the parser
parser = argparse.ArgumentParser(description='Train LSTM_encoder+PointerLSTM_decoder for sorting numbers')

# Required positional argument
parser.add_argument('n_steps', type=int,
                    help='Sequence length (recommended: 5)')

# Optional argument
parser.add_argument('--n_examples', type=int, default=10000,
                    help='n_examples (recommended: 10000)')

# Optional argument
parser.add_argument('--upper_limit', type=int, default=5,
                    help='upper_limit of the input data (recommended: 5)')

# Optional argument
parser.add_argument('--epochs', type=int, default=10,
                    help='no_of_epochs to be trained for (recommended: 10)')

# Optional argument
parser.add_argument('--test_sequence', type=int, nargs='+', default=[3, 2, 1, 0, 4],
                    help='test_sequence to view predicted output sequence')

args = parser.parse_args()

n_steps = args.n_steps
upper_limit = args.upper_limit
n_examples = args.n_examples
epochs = args.epochs
split_at = int(.9*n_examples)
batch_size = 100

hidden_size = 64
# weights_file = 'model_weights/achbogga_model_weights_{}_steps_{}.hdf5'.format(n_steps, hidden_size)

x, y = generate_sort_data(n_steps, n_examples, upper_limit)
x = np.expand_dims(x, axis=2)
# prep position indicators

n = np.arange(n_steps)
nn = np.tile(n, (x.shape[0], 1))
nn = np.expand_dims(nn, axis=2)

xx = np.dstack((x, nn))

YY = []
for y_ in y:
    YY.append(to_categorical(y_))
YY = np.asarray(YY)

x_train = xx[:split_at]
x_test = xx[split_at:]

y_test = y[split_at:]
YY_train = YY[:split_at]
YY_test = YY[split_at:]

assert (n_steps == x_train.shape[1])
n_features = x_train.shape[2]

print("building model...")
main_input = Input(shape=(n_steps, n_features), name='main_input')

encoder = LSTM(output_dim=hidden_size, return_sequences=True, name="encoder")(main_input)
decoder = PointerLSTM(hidden_size, output_dim=hidden_size, name="decoder")(encoder)

model = Model(input=main_input, output=decoder)

# print(("loading weights from {}...".format(weights_file)))
# try:
#     model.load_weights(weights_file)
# except IOError:
#     print("no weights file, starting anew.")

model.compile(optimizer='rmsprop',
              loss='categorical_crossentropy',
              metrics=['accuracy'])

print('training and saving model weights each epoch...')

validation_data = (x_test, YY_test)

history = model.fit(x_train, YY_train, epochs=epochs, batch_size=batch_size,
                    validation_data=validation_data)

# p = model.predict(x_test)
# print("test accuracy: ", accuracy_score(YY_test, p))

x, y = generate_x_y_for_inference(args.test_sequence, n_steps)
x = np.expand_dims(x, axis=2)
assert (n_steps == x.shape[1])

n = np.arange(n_steps)
nn = np.tile(n, (x.shape[0], 1))
nn = np.expand_dims(nn, axis=2)

xx = np.dstack((x, nn))

p = model.predict(xx)

print ("input_sequence: ", xx[0][:,0].flatten())

for i in range(len(xx)):
    y_true = xx[i][:,0][y].flatten()
    y_pred = xx[i][:,0][p.argmax(axis=1)].flatten()
    print(("ground_truth:", y_true))
    print(("predicted_sequence:     ", y_pred))
    print("accuracy: ", accuracy_score(y_true, y_pred))

max_stars_repo_name: achbogga/pointer-networks-experiments

max_issues_repo_name: achbogga/pointer-networks-experiments

max_forks_repo_name: achbogga/pointer-networks-experiments


text:
import cv2
import operator
import numpy as np
from imageio import imread
from skimage.transform import resize
from PIL import Image

def load_and_align_images(filepaths, margin,image_size):

    aligned_images = []
    for filepath in filepaths:
        img = imread(filepath)
        aligned = resize(img, (image_size, image_size), mode='reflect')
        aligned_images.append(aligned)

    return np.array(aligned_images)

def read_image(image_path):
    return Image.open(image_path)

def read_image_to_arr(image_path):
    return np.array(Image.open(image_path))

def prewhiten(x):
    if x.ndim == 4:
        axis = (1, 2, 3)
        size = x[0].size
    elif x.ndim == 3:
        axis = (0, 1, 2)
        size = x.size
    else:
        raise ValueError('Dimension should be 3 or 4')

    mean = np.mean(x, axis=axis, keepdims=True)
    std = np.std(x, axis=axis, keepdims=True)
    std_adj = np.maximum(std, 1.0/np.sqrt(size))
    y = (x - mean) / std_adj
    return y


text:
# Author: Jintao Huang
# Email: hjt_study@qq.com
# Date:

import torch
from torch import Tensor
import numpy as np

def gather(input: Tensor, dim: int, index: Tensor) -> Tensor:
    """只支持input, index, return的shape一致

    :param input: shape[...]
    :param dim:
    :param index: Long. shape[...]
    :return: shape[...]
    """
    idxs = [range(length) for length in input.shape]
    idxs = np.ix_(*idxs)
    idxs = [torch.tensor(arr, dtype=torch.long) for arr in idxs]  # 转Tensor
    idxs[dim] = index
    #
    out = input[idxs]
    return out


text:
import setup_path
import cv2
import airsim
import time
import numpy as np
from PIL import Image
from controller import PS4Controller

MAX_SPEED = 1
MAX_THETA = 1

if __name__ == '__main__':
    client = airsim.CarClient()
    client.confirmConnection()
    client.enableApiControl(True)
    car_controls = airsim.CarControls()

    controller = PS4Controller()
    controller.start()

    """PS4 Controller's axes map

                 Axis 1                  Axis 4

                  -1                      -1
                   ^                       ^
                   |                       |
    Axis 0  -1 <---0---> 1   Axis 3 -1 <---0---> 1
                   |                       |
                   v                       v
                   1                       1
    """

    while True:
        _, axis_data, _ = controller.read()
        raw_speed = axis_data[1]
        raw_theta = axis_data[3]
        speed = np.interp(-raw_speed, (-1, 1), (-MAX_SPEED, MAX_SPEED))
        theta = np.interp(raw_theta, (-1, 1), (-MAX_THETA, MAX_THETA))

        car_controls.throttle = speed
        car_controls.steering = theta
        client.setCarControls(car_controls)

        car_state = client.getCarState()
        print(f'Speed {car_state.speed}, Gear {car_state.gear}')

        # Getting images
        responses = client.simGetImages([
            # png format
            airsim.ImageRequest("0", airsim.ImageType.Scene),
            # uncompressed RGB array bytes
            airsim.ImageRequest("1", airsim.ImageType.Scene, False, False),
            # floating point uncompressed image
            airsim.ImageRequest("1", airsim.ImageType.DepthPlanner, True)])

        response = responses[1]
        # get numpy array
        img1d = np.fromstring(response.image_data_uint8, dtype=np.uint8)
        # reshape array to 3 channel image array H X W X 3
        img_rgb = img1d.reshape(response.height, response.width, 3)

        cv2.imshow('image', img_rgb)
        if cv2.waitKey(1) & 0xFF == ord('q'):
            cv2.destroyAllWindows()
            break

    # restore to original state
    client.reset()

    client.enableApiControl(False)


text:
import numpy as np

lines = open('output.txt').read().strip().splitlines()
data = np.stack([np.array(list(line)) for line in lines])

scaffolds = data == '#'
intersections = np.zeros_like(scaffolds)
Nx, Ny = scaffolds.shape
for i in range(1, Nx-1):
    for j in range(1, Ny-1):
        if scaffolds[i, j]:
            is_int = (scaffolds[i-1,j] and
                      scaffolds[i+1,j] and
                      scaffolds[i,j-1] and
                      scaffolds[i,j+1])
            intersections[i,j] = is_int

xs, ys = np.where(intersections)
print((xs * ys).sum())

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

trans = Pipeline(steps=[
    ('scl', StandardScaler(with_mean=True, with_std=True, copy=True)),
    ('pca', PCA(
        n_components='mle', svd_solver='full',
        whiten=False, copy=True))
])

meta = {
    'id': 'dim2',
    'name': 'PCA Minka MLE',
    'description': "Minka's MLE to guess n_components",
    'keywords': [
        'dimensionality reduction', 'principal component anlysis',
        'StandardScaler', 'PCA', 'Minka MLE'],
    'feature_names_prefix': 'dim_minka'
}

"""Example

from verto.dim2 import trans, meta
from datasets.demo1 import X_train
from seasalt import create_feature_names
import pandas as pd
import numpy as np

X_new = trans.fit_transform(X_train).astype(np.float16)
names = create_feature_names(meta['feature_names_prefix'], X_new.shape[1])
df = pd.DataFrame(data=X_new, columns=names)
df
"""

max_stars_repo_name: kmedian/potpourri

max_issues_repo_name: kmedian/potpourri

max_forks_repo_name: kmedian/potpourri


text:
from taurex.contributions import Contribution
from functools import lru_cache
import pycuda.autoinit
from pycuda.compiler import SourceModule
from pycuda.gpuarray import GPUArray, zeros
import pycuda.driver as drv
import numpy as np
import math
@lru_cache(maxsize=400)
def _sum_kernal(nlayers, ngrid):

    code = f"""

    __global__ void sum_sigma(double* dest, const double* __restrict__ arr)
    {{
        unsigned int i = (blockIdx.x * blockDim.x) + threadIdx.x;
        unsigned int j = (blockIdx.y * blockDim.y) + threadIdx.y;

        if (i >= {ngrid})
            return;

        if (j >= {nlayers})
            return;

        dest[j*{ngrid} + i] += arr[j*{ngrid} + i];
    }}
    """
    mod = SourceModule(code)
    return mod.get_function('sum_sigma')

@lru_cache(maxsize=400)
def _contribute_tau_kernal(nlayers, grid_size, with_sigma_offset=False, start_layer=0):
    extra = '+layer'
    if with_sigma_offset:
        extra = ''

    code = f"""
    __global__ void contribute_tau(double* dest, const double* __restrict__ sigma,
                                   const double* __restrict__ density, const double* __restrict__ path,
                                   const int* __restrict__ startK, const int* __restrict__ endK,
                                   const int* __restrict__ density_offset, const int total_layers)
    {{
        unsigned int i = (blockIdx.x * blockDim.x) + threadIdx.x;

        if ( i >= {grid_size} )
            return;

        for (unsigned int layer={start_layer}; layer<{nlayers}; layer++)
        {{
            unsigned int _startK = startK[layer];
            unsigned int _endK = endK[layer];
            for (unsigned int k = _startK; k < _endK; k++)
            {{
                double _path = path[layer*{nlayers} + k];
                double _density = density[k+density_offset[layer]];
                dest[layer*{grid_size} + i] += sigma[(k{extra})*{grid_size} + i]*_path*_density;
            }}

        }}

    }}

    """
    mod = SourceModule(code)
    func = mod.get_function('contribute_tau')
    func.prepare('PPPPPPPi')
    return func

@lru_cache(maxsize=400)
def _contribute_tau_kernal_II(nlayers, grid_size, with_sigma_offset=False, start_layer=0):
    extra = '+layer'
    if with_sigma_offset:
        extra = ''

    code = f"""
    __global__ void contribute_tau(double* dest, const double* __restrict__ sigma,
                                   const double* __restrict__ density, const double* __restrict__ path,
                                   const int* __restrict__ startK, const int* __restrict__ endK,
                                   const int* __restrict__ density_offset, const int total_layers)
    {{
        __shared__ double path_cache[{nlayers}];
        __shared__ double dens_cache[{nlayers}];

        const unsigned int grid = (blockIdx.x * blockDim.x) + threadIdx.x;
        const unsigned int layer = (blockIdx.y * blockDim.y) + threadIdx.y; + {start_layer};

        unsigned int _s_layer = threadIdx.y*blockDim.x + threadIdx.x;

        if (_s_layer < {nlayers})
        {{
            path_cache[_s_layer] = path[layer*{nlayers} + _s_layer];
            dens_cache[_s_layer] = density[_s_layer];
        }}

        __syncthreads();

        if ( grid >= {grid_size} )
            return;
        if (layer >= {nlayers})
            return;

        const unsigned int _startK = startK[layer];
        const unsigned int _endK = endK[layer];
        const unsigned int _offset = density_offset[layer];
        double _result = 0.0;
        for (unsigned int k = _startK; k < _endK; k++)
        {{
            double _path = path_cache[k];
            double _density = dens_cache[k+_offset];
            _result += sigma[(k{extra})*{grid_size} + grid]*_path*_density;
        }}
        dest[layer*{grid_size} + grid] += _result;

    }}

    """
    mod = SourceModule(code,options=['-O3'])
    func = mod.get_function('contribute_tau')
    func.prepare('PPPPPPPi')
    return func

def cuda_contribute_tau_old(startK, endK, density_offset, sigma, density, path,
                        nlayers, ngrid, tau=None, with_sigma_offset=False, start_layer=0,total_layers=None,
                        stream=None):

    kernal = _contribute_tau_kernal(nlayers, ngrid, with_sigma_offset=with_sigma_offset, start_layer=start_layer)
    my_tau = tau
    if total_layers is None:
        total_layers = nlayers
    if my_tau is None:
        my_tau = GPUArray(shape=(nlayers,ngrid),dtype=np.float64)

    #THREAD_PER_BLOCK_X = 32
    #THREAD_PER_BLOCK_Y = 32

    THREAD_PER_BLOCK_X = 128
    THREAD_PER_BLOCK_Y = 1
    #NUM_BLOCK_Y = int(math.ceil((total_layers)/THREAD_PER_BLOCK_Y))
    NUM_BLOCK_X = int(math.ceil((ngrid)/THREAD_PER_BLOCK_X))
    NUM_BLOCK_Y = 1

    kernal.prepared_call(
        (NUM_BLOCK_X, NUM_BLOCK_Y, 1),
        (THREAD_PER_BLOCK_X, THREAD_PER_BLOCK_Y, 1),
        my_tau.gpudata, sigma.gpudata, density.gpudata, path.gpudata, startK.gpudata, endK.gpudata, density_offset.gpudata,np.int32(total_layers))
    if tau is None:
        return my_tau

def cuda_contribute_tau(startK, endK, density_offset, sigma, density, path,
                        nlayers, ngrid, tau=None, with_sigma_offset=False, start_layer=0,total_layers=None,
                        stream=None):

    kernal = _contribute_tau_kernal_II(nlayers, ngrid, with_sigma_offset=with_sigma_offset, start_layer=start_layer)
    my_tau = tau
    if total_layers is None:
        total_layers = nlayers
    if my_tau is None:
        my_tau = GPUArray(shape=(nlayers,ngrid),dtype=np.float64)

    THREAD_PER_BLOCK_X = 32
    THREAD_PER_BLOCK_Y = 32

    #THREAD_PER_BLOCK_X = 128
    #THREAD_PER_BLOCK_Y = 1
    NUM_BLOCK_Y = int(math.ceil((total_layers)/THREAD_PER_BLOCK_Y))
    NUM_BLOCK_X = int(math.ceil((ngrid)/THREAD_PER_BLOCK_X))
    #NUM_BLOCK_Y = 1

    kernal.prepared_call(
        (NUM_BLOCK_X, NUM_BLOCK_Y, 1),
        (THREAD_PER_BLOCK_X, THREAD_PER_BLOCK_Y, 1),
        my_tau.gpudata, sigma.gpudata, density.gpudata, path.gpudata, startK.gpudata, endK.gpudata, density_offset.gpudata,np.int32(total_layers))
    if tau is None:
        return my_tau

class CudaContribution(Contribution):

    def __init__(self, name):
        super().__init__(name)
        self._is_cuda_model = False

    def contribute(self, model, start_layer, end_layer,
                   density_offset, layer, density, tau, path_length=None, with_sigma_offset=False, streams=None):
        """
        Computes an integral for a single layer for the optical depth.

        Parameters
        ----------
        model: :class:`~taurex.model.model.ForwardModel`
            A forward model

        start_layer: int
            Lowest layer limit for integration

        end_layer: int
            Upper layer limit of integration

        density_offset: int
            offset in density layer

        layer: int
            atmospheric layer being computed

        density: :obj:`array`
            density profile of atmosphere

        tau: :obj:`array`
            optical depth to store result

        path_length: :obj:`array`
            integration length

        """

        streams = None

        self.debug(' %s %s %s %s %s %s %s', start_layer, end_layer,
                   density_offset, layer, density, tau, self._ngrid)

        if streams is None:
            cuda_contribute_tau(start_layer, end_layer, density_offset,
                                self.sigma_xsec, density, path_length,
                                self._nlayers, self._ngrid, tau,with_sigma_offset)
        else:
            num_streams=len(streams)
            split = self._nlayers//num_streams
            stream_vals = [(x*split, (x+1)*split, streams[x]) for x in range(num_streams-1)] + [((num_streams-1)*split,self._nlayers,streams[-1])]

            for start, end, stream in stream_vals:
                cuda_contribute_tau(start_layer, end_layer, density_offset,
                                    self.sigma_xsec, density, path_length,
                                    end, self._ngrid, tau,with_sigma_offset,start_layer=start,total_layers=self._nlayers,stream=stream)
        self.debug('DONE')

    def prepare(self, model, wngrid):
        """

        Used to prepare the contribution for the calculation.
        Called before the forward model performs the main optical depth
        calculation. Default behaviour is to loop through :func:`prepare_each`
        and sum all results into a single cross-section.

        Parameters
        ----------
        model: :class:`~taurex.model.model.ForwardModel`
            Forward model

        wngrid: :obj:`array`
            Wavenumber grid
        """

        self._ngrid = wngrid.shape[0]
        self._nlayers = model.nLayers

        sigma_xsec = zeros(shape=(self._nlayers, self._ngrid), dtype=np.float64)

        sum_kernal = _sum_kernal(self._nlayers, self._ngrid)
        THREAD_PER_BLOCK_X = 16
        THREAD_PER_BLOCK_Y = 16

        NUM_BLOCK_X = int(math.ceil(self._ngrid/THREAD_PER_BLOCK_X))
        NUM_BLOCK_Y = int(math.ceil(self._nlayers/THREAD_PER_BLOCK_Y))

        for gas, sigma in self.prepare_each(model, wngrid):
            self.debug('Gas %s', gas)
            self.debug('Sigma %s', sigma)
            sum_kernal(sigma_xsec, sigma, block=(THREAD_PER_BLOCK_X, THREAD_PER_BLOCK_Y,1), grid=(NUM_BLOCK_X, NUM_BLOCK_Y,1))

        self.sigma_xsec = sigma_xsec
        self.debug('Final sigma is %s', self.sigma_xsec)
        self.info('Done')

reason: import numpy,import pycuda,from pycuda