text:
from collections import namedtuple

from astropy import convolution
import numpy as np
import math

import extinction

from .define_filters import define_filters, get_cenwave

#Like a light weight class with attributes wavelength and flux
#spectrum1d = namedtuple("spectrum1d", ['wave', 'flux'])
class spectrum1d(object):
    def __init__(self, wave, flux, error=None):
        #Could be np.asarray, but change fit_spectral_lines.py
        self.wave = np.float_(wave)
        self.flux = np.float_(flux)
        self.error = np.float_(error)

def degrade_resolution(flux, disp_in, disp_out):
    '''
    Convolve spectrum with a 1D gaussian kernel to degrade dispersion from disp_in to disp_out
    '''
    fwhm = math.sqrt(disp_out**2 - disp_in**2)
    stddev = fwhm/(2.*math.sqrt(2.*math.log(2.)))
    gauss_kernel=convolution.Gaussian1DKernel(stddev=stddev)
    degraded_flux = convolution.convolve(flux, gauss_kernel)
    return degraded_flux

def bin_data(data_arr, bin_size=None, bins=[], sum=False, median=False):
    '''
    Bin 1D array of data
    bin_size - this keyword can be set to bin data to a certain number of points. The last point
            may have fewer points
    bins - this keyword can be set to an array of bin edges. Values will always be treated as the left bin edge
    sum - if True, the returned result is the sum of all points in each bin
    median - if True, the returned result is the median of all points in each bin

    Note: one and only one of sum and median must be set to True
          one and only one of bin_size and bins must be set
    '''
    binned_array = []
    if bin_size:
        assert len(bins) == 0, 'cannot set both bin_size and bins'
        assert isinstance(bin_size, int), 'bin_size must be an integer'
        arr_len = len(data_arr)
        nbins = int(np.ceil(arr_len/bin_size))

        for bin in range(nbins):
            if sum:
                assert median is False, 'cannot set both median and sum to be True'
                binned_array.append(np.sum(data_arr[bin*bin_size:(bin+1)*bin_size]))
            if median:
                assert sum is False, 'cannot set both median and sum to be True'
                binned_array.append(np.median(data_arr[bin*bin_size:(bin+1)*bin_size]))
        return np.array(binned_array)
    if len(bins) != 0:
        assert bin_size is None, 'cannot set both bin_size and bins'
        for left_bin_edge, right_bin_edge in zip(bins[:-1], bins[1:]):
            assert isinstance(left_bin_edge, int) or isinstance(left_bin_edge, np.int64), 'bin edges must be integers'
            if sum:
                assert median is False, 'cannot set both median and sum to be True'
                binned_array.append(np.sum(data_arr[left_bin_edge:right_bin_edge]))
            if median:
                assert sum is False, 'cannot set both median and sum to be True'
                binned_array.append(np.median(data_arr[left_bin_edge:right_bin_edge]))

        return np.array(binned_array)

def calc_wavelength(header, pixels):
    if ('CTYPE1' in header.keys()) and (len(header['CTYPE1'])>1):
        if header['CTYPE1'] != 'LINEAR':
            keep_going = input('Attempting to apply a linear solution to non-linear parameters. Continue anyways? (y), n ')
            if keep_going == 'n':
                sys.exit()
    else:
        print('WARNING: No solution type specified, assuming linear')
    CRVAL1 = header['CRVAL1']
    if 'CRPIX1' in header.keys():
        CRPIX1 = header['CRPIX1']
    else:
        CRPIX1 = 0
    if 'CD1_1' in header.keys():
        CD1_1 = header['CD1_1']
    elif 'CDELT1' in header.keys():
        CD1_1 = header['CDELT1']
    elif 'PC1_1' in header.keys():
        CD1_1 = header['PC1_1']
    else:
        sys.exit('Could not identify linear wavelength term (tried CD1_1, CDELT1, PC1_1)')

    wavelength = CRVAL1 + CD1_1*(pixels - CRPIX1)
    return wavelength

def apply_redshift(wl, redshift, redden=False):
    '''
    Redden or deredden a wavelength (array) with a given redshift

    Inputs:
        wl: wavelength (can be single value or array)
        redshift: redshift value
        redden: If False, code will calculate rest wavelength assuming input is observed
                If True, code will calculate observed wavelength assuming input is rest
    '''
    if isinstance(wl, list):
        wl = np.array(wl)
    if redden is False:
        new_wl = wl/(1+redshift)
    elif redden is True:
        new_wl = wl * (1+redshift)
    return new_wl

def scale_spectra(spec1, spec2, wlmin = None, wlmax=None, scale_factor=False):
    '''
    Scale spec1 to spec2 by integrating from wlmin to wlmax and using
    the area as a scaling factor
    Inputs:
        spec1, spec2: two spectrum1d objects with attributes wave and flux
        wlmin: shortest wavelength to be used in scaling
        wlmax: longest wavelength to be used in scaling
    Outputs:
        spec1: spec1d object with scaled flux
        scale_factor (optional): if scale_factor = True, then the multiplicative factor is returned
    '''
    if wlmin is None:
        wlmin = max(spec1.wave.min(), spec2.wave.min())
    if wlmax is None:
        wlmax = min(spec1.wave.max(), spec2.wave.max())
    windx_1 = (spec1.wave >= wlmin) &  (spec1.wave <= wlmax)
    windx_2 = (spec2.wave >= wlmin) & (spec2.wave <= wlmax)
    delta_wl1 = np.mean(spec1.wave[1:] - spec1.wave[:-1])
    delta_wl2 = np.mean(spec2.wave[1:] - spec2.wave[:-1])
    area1 = np.sum(spec1.flux[windx_1])*delta_wl1
    area2 = np.sum(spec2.flux[windx_2])*delta_wl2
    spec1 = spectrum1d(spec1.wave, spec1.flux/area1*area2)
    if scale_factor is True:
        return spec1, area2/area1
    else:
        return spec1

def correct_for_galactic_extinction(spec, E_BV, R_V=3.1):
    '''
    Correct flux for galactic (Milky Way) extinction using a Cardelli law
    Inputs:
        spec: spectrum1d object to be corrected (has wave and flux attributes)
        E_BV: E(B-V) for correction
        R_V: default to 3.1 values
    Output:
        spectrum1d object with the dust corrected flux
    '''
    A_V = R_V*E_BV
    new_flux = extinction.apply(-extinction.ccm89(spec.wave, A_V, R_V), spec.flux)
    return spectrum1d(spec.wave, new_flux)

def calc_extinction(E_BV, band, Rv=3.1, E_BV_err=0):
    '''
    Calculate extinction using the Cardelli law
    E_BV: float
        E(B-V) value
    band: str
        filter corresponding to filter in define_filters dictionary
    Rv: float
        R(V) set to 3.1 by default
    E_BV_err: float
        1 sigma error in E(B-V), set to 0 by default
    '''
    Av = Rv*E_BV
    Av_err = Rv**2 * E_BV_err**2
    cenwave = np.array([float(get_cenwave(band))])
    A_band = extinction.ccm89(cenwave, Av, Rv)
    if Av == 0:
        A_band_err = 0.0
    else:
        A_band_err = np.sqrt(Av_err * A_band/Av)
    return A_band, A_band_err

############################
#Tests
############################
def test_full_bins():
    test_arr = np.arange(100)
    binned_array1 = bin_data(test_arr, bin_size = 10, sum=True)
    assert (binned_array1 == np.arange(45, 1000, 100)).all(), 'failed test on full bins'

def test_incomplete_last_bin():
    test_arr = np.arange(100)
    binned_array2 = bin_data(test_arr, bin_size = 15, sum=True)
    assert (binned_array2 == np.array([105, 330, 555, 780, 1005, 1230, 945])).all(), 'failed test on last bin not being full'

def test_premade_bins():
    test_arr = np.arange(100)
    bins = np.arange(0, 101, 10)
    binned_array1 = bin_data(test_arr, bin_size = 10)
    binned_array3 = bin_data(test_arr, bins = bins)
    assert (binned_array3 == binned_array1).all(), 'bins not functioning as intended'

def test_scale_spectra():
    spec1 = spec.spectrum1d(np.arange(6000, 8000, 1), np.ones(2000))
    spec2 = spec.spectrum1d(np.arange(6000, 8000, 1), np.ones(2000)*4.2)
    new_spec, scale_factor = spec.scale_spectra(spec1, spec2, scale_factor=True)
    plt.figure()
    plt.plot(spec1.wave, spec1.flux, ':', label='original flux')
    plt.plot(spec1.wave, spec1.flux*scale_factor, label='scaled flux')
    plt.plot(spec2.wave, spec2.flux, '--', label='flux to scale to')
    plt.legend(loc='best')
    assert np.isclose(scale_factor,4.2), 'scale_factor={}'.format(scale_factor)

def test_scale_spectra2():
    flux1 = np.ones(2000)
    flux1[1000:1500]+=1
    flux2 = np.ones(2000)+4.2
    flux2[1000:1500]+=3
    spec1 = spec.spectrum1d(np.arange(6000, 8000, 1), flux1)
    spec2 = spec.spectrum1d(np.arange(6000, 8000, 1), flux2)
    new_spec, scale_factor = spec.scale_spectra(spec1, spec2, scale_factor=True, wlmin=7000, wlmax=7500)
    plt.figure()
    plt.plot(spec1.wave, spec1.flux, ':', label='original flux')
    plt.plot(spec1.wave, spec1.flux*scale_factor, label='scaled flux')
    plt.plot(spec2.wave, spec2.flux, '--', label='flux to scale to')
    plt.legend(loc='best')
    assert np.isclose(scale_factor,4.101098901098902), 'scale_factor={}'.format(scale_factor)

def test_scale_spectra3():
    flux1 = np.ones(2000)
    flux1[1000:1500]+=1
    flux2 = np.ones(2000)+4.2
    flux2[1000:1500]+=3
    spec1 = spec.spectrum1d(np.arange(6000, 8000, 1), flux1)
    spec2 = spec.spectrum1d(np.arange(6000, 8000, 1), flux2)
    new_spec, scale_factor = spec.scale_spectra(spec1, spec2, scale_factor=True, wlmin=7200, wlmax=7250)
    plt.figure()
    plt.plot(spec1.wave, spec1.flux, ':', label='original flux')
    plt.plot(spec1.wave, spec1.flux*scale_factor, label='scaled flux')
    plt.plot(spec2.wave, spec2.flux, '--', label='flux to scale to')
    plt.legend(loc='best')
    assert np.isclose(scale_factor,4.101098901098902, rtol=0.1), 'scale_factor={}'.format(scale_factor)

max_stars_repo_name: abostroem/utilities

max_issues_repo_name: abostroem/utilities

max_forks_repo_name: abostroem/utilities

reason: import numpy,from astropy


text:
import sys
import pandas as pd
import numpy as np
import itertools
from sklearn.preprocessing import RobustScaler
from sklearn.ensemble import AdaBoostClassifier
from evaluate_model import evaluate_model

dataset = sys.argv[1]

pipeline_components = [RobustScaler, AdaBoostClassifier]
pipeline_parameters = {}

learning_rate_values = [0.01, 0.1, 0.5, 1.0, 10.0, 50.0, 100.0]
n_estimators_values = [10, 50, 100, 500]
random_state = [324089]

all_param_combinations = itertools.product(learning_rate_values, n_estimators_values, random_state)
pipeline_parameters[AdaBoostClassifier] = [{'learning_rate': learning_rate, 'n_estimators': n_estimators, 'random_state': random_state}
                                    for (learning_rate, n_estimators, random_state) in all_param_combinations]

evaluate_model(dataset, pipeline_components, pipeline_parameters)

ext: py

max_stars_repo_name: lacava/sklearn-benchmarks

max_issues_repo_name: LeMGarouani/AMLBID

max_forks_repo_name: arunsinghyadav/sklearn-benchmarks

reason: import numpy


text:
#
# utils.py
#

from torch.autograd.functional import jacobian
import torch
import numpy as np
from numpy import pi
import timeit
import random
import sys
try:
    import resource
except ImportError:
    resource = None

from contextlib import contextmanager

from scipy.interpolate import UnivariateSpline

from pykalman import KalmanFilter
from scipy.linalg import expm

from tqdm import tqdm

def _check_param_device(param, old_param_device):
    """This helper function is to check if the parameters are located
    in the same device. Currently, the conversion between model parameters
    and single vector form is not supported for multiple allocations,
    e.g. parameters in different GPUs, or mixture of CPU/GPU.
    Arguments:
        param ([Tensor]): a Tensor of a parameter of a model
        old_param_device (int): the device where the first parameter of a
                                model is allocated.
    Returns:
        old_param_device (int): report device for the first time
    """

    # Meet the first parameter
    if old_param_device is None:
        old_param_device = param.get_device() if param.is_cuda else -1
    else:
        warn = False
        if param.is_cuda:  # Check if in same GPU
            warn = (param.get_device() != old_param_device)
        else:  # Check if in CPU
            warn = (old_param_device != -1)
        if warn:
            raise TypeError('Found two parameters on different devices, '
                            'this is currently not supported.')
    return old_param_device

def qqdot_to_q(qqdot, dt):
    """
    Convert a timeseries of [q,qdot]_t to a T+1 timeseires of [q]_t
    Args:
        qqdot (torch.tensor): (B,T,qdim*2) [q,qdot]
        dt (float): timestemp
    Returns:
        q (torch.tensor): (B,T+1,qdim)
    Notes:
        Finds q_{1:T+1} s.t qqdot_t approx [0.5*(q_t+1 + q_t),
                                            (1./dt)*(q_t+1 - q_t)]
    """

    B,T,qdim2 = qqdot.shape
    qdim = qdim2//2

    # find a matrix D that maps q to qddot
    a = torch.randn(T+1,qdim).reshape(-1)
    a.requires_grad = True
    def lam(a):
        a_ = a.reshape(T+1,qdim)
        b1 = 0.5*(a_[1:] + a_[:-1])
        b2  = (1./dt) * (a_[1:] - a_[:-1])
        b = torch.cat([b1,b2],dim=-1)
        return b.reshape(-1)
    D = jacobian(lam, a)

    # solve least squares problem
    q = torch.zeros(B,(T+1)*qdim)
    qqdot = qqdot.reshape(B,-1,1)
    for b in range(B):
        sol = torch.lstsq(qqdot[b], D)[0]
        q[b,:] = sol[:(T+1)*qdim].reshape(-1)

    return q.reshape(B,T+1,qdim)

def bfill_uppertriangle(A, vec):
	ii, jj = np.triu_indices(A.size(-2), k=1, m=A.size(-1))
	A[..., ii, jj] = vec
	return A

def bfill_lowertriangle(A, vec):
    ii, jj = np.tril_indices(A.size(-2), k=-1, m=A.size(-1))
    A[..., ii, jj] = vec
    return A

def bfill_diagonal(A, vec):
    ii, jj = np.diag_indices(min(A.size(-2), A.size(-1)))
    A[..., ii, jj] = vec
    return A

def parameter_grads_to_vector(parameters):
    """Convert parameters to one vector
    Arguments:
        parameters (Iterable[Tensor]): an iterator of Tensors that are the
            parameters of a model.
    Returns:
        The parameters represented by a single vector
    """
    # Flag for the device where the parameter is located
    param_device = None

    vec = []
    for param in parameters:
        # Ensure the parameters are located in the same device
        param_device = _check_param_device(param, param_device)

        if param.grad is not None:
            vec.append(param.grad.view(-1))
        else:
            zrs = torch.zeros_like(param)
            vec.append(zrs.view(-1))
    return torch.cat(vec)

def vector_to_parameter_grads(vec, parameters):
    """Convert one vector to the parameters
    Arguments:
        vec (Tensor): a single vector represents the parameters of a model.
        parameters (Iterable[Tensor]): an iterator of Tensors that are the
            parameters of a model.
    """
    # Ensure vec of type Tensor
    if not isinstance(vec, torch.Tensor):
        raise TypeError('expected torch.Tensor, but got: {}'.format(torch.typename(vec)))
    # Flag for the device where the parameter is located
    param_device = None

    # Pointer for slicing the vector for each parameter
    pointer = 0
    for param in parameters:
        # Ensure the parameters are located in the same device
        param_device = _check_param_device(param, param_device)

        # ensure that param requires grad
        assert param.requires_grad

        # The length of the parameter
        num_param = param.numel()
        # Slice the vector, reshape it, and replace the old grad of the parameter
        param.grad.data = vec[pointer:pointer + num_param].view_as(param).data

        # Increment the pointer
        pointer += num_param

def smooth_and_diff(x, dt, nders=2, **kwargs):
    """
    Smooth a set of timeseries and differentiate them.
    Args:
        x (torch.tensor): (B,T,n) timeseries
        dt (float): time-step
        nders (int): number of derivatives to take
        **kwargs: arguments to UnivariateSpline
    Returns:
        smooth_series (list of torch.tensor): list of (B,T,n) smooth tensors (xsm, dxsmdt, ...)
    """

    retvals = [torch.zeros_like(x) for i in range(nders+1)]

    B, T, N = x.shape
    t = np.arange(T) * dt
    for b in range(B):
        for n in range(N):
            ser = x[b,:,n].detach().numpy()
            spl = UnivariateSpline(t, ser, **kwargs)

            retvals[0][b,:,n] = torch.tensor(spl(t))
            for d in range(nders):
                retvals[d+1][b,:,n] = torch.tensor(spl.derivative(d)(t))

    return retvals

def kalman_smooth_and_diff(x, dt, nders=2, em_Q=False):
    """
    Smooth a set of timeseries and differentiate them.
    Args:
        x (torch.tensor): (B,T,n) timeseries
        dt (float): time-step
        nders (int): number of derivatives to take
        em_Q (bool): learn transition cov?
    Returns:
        smooth_series (list of torch.tensor): list of (B,T,n) smooth tensors (xsm, dxsmdt, ...)
    """

    retvals = [torch.zeros_like(x) for i in range(nders+1)]

    em_vars = ['initial_state_mean', 'initial_state_covariance',
                        'observation_covariance']
    if em_Q:
        em_vars += ['transition_covariance']

    if nders != 2:
        raise NotImplementedError

    A = np.array([[0.,1.,0.],
                  [0.,0.,1.],
                  [0.,0.,0.]])
    Ad = expm(A * dt)
    Bd = np.array([[1.,0.,0.]])

    Q = np.diag([0.001,0.001,1.0])

    B, T, N = x.shape
    for b in tqdm(range(B)):
        for n in range(N):
            ser = x[b,:,n:n+1].detach().numpy()

            kf = KalmanFilter(transition_matrices=Ad, observation_matrices=Bd, transition_covariance=Q)
            kf.em(ser, em_vars=em_vars)
            sm_x, _ = kf.smooth(ser)
            sm_x = torch.tensor(sm_x)

            for i in range(3):
                retvals[i][b,:,n] = sm_x[:,i]

    return retvals

max_stars_repo_name: sisl/delsmm

max_issues_repo_name: albernsrya/delsmm

max_forks_repo_name: albernsrya/delsmm

reason: import numpy,from numpy,from scipy


text:
import numpy as np

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

dot_product = np.dot(a, b)
print(dot_product)

max_stars_repo_name: specbug/nnfs

max_issues_repo_name: specbug/nnfs

max_forks_repo_name: specbug/nnfs


text:
import argparse
import pickle
from collections import namedtuple
from itertools import count
import os
import time
import numpy as np
import matplotlib.pyplot as plt
import gym
import highway_env
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from torch.distributions import Normal, Categorical
from torch.utils.data.sampler import BatchSampler, SubsetRandomSampler
import matplotlib.pyplot as plt

class Actor(nn.Module):
    def __init__(self, num_inputs, num_actions):
        super(Actor, self).__init__()
        self.fc1 = nn.Linear(num_inputs, 100)
        self.action_head = nn.Linear(100, num_actions)

    def forward(self, x):
        x = torch.flatten(x, start_dim=1)
        x = F.relu(self.fc1(x))
        action_prob = F.softmax(self.action_head(x), dim=1)
        return action_prob

class Critic(nn.Module):
    def __init__(self, num_inputs, num_actions):
        super(Critic, self).__init__()
        self.fc1 = nn.Linear(num_inputs, 100)
        self.state_value = nn.Linear(100, 1)

    def forward(self, x):
        x = torch.flatten(x, start_dim=1)
        x = F.relu(self.fc1(x))
        value = self.state_value(x)
        return value

class PPO():
    def __init__(self):
        super(PPO, self).__init__()

    def init_hyp(self):
        self.buffer = []
        self.counter = 0
        self.training_step = 0
        self.clip_param = 0.2
        self.max_grad_norm = 0.5
        self.ppo_update_time = 10
        self.buffer_capacity = 1000
        self.batch_size = 32
        self.gamma = 0.99
        self.render = True
        self.seed = 1
        self.log_interval = 10
        self.checkpoint_path = "single-agent/env-highway/ppo/checkpoints/ppo-highway.pth"

    def init_env(self):
        self.env = gym.make('highway-v0')
        self.num_state = self.env.observation_space.shape[0]
        self.num_action = self.env.action_space.n
        torch.manual_seed(self.seed)
        self.env.seed(self.seed)

        print("obs space = {} action space = {} ".format(self.env.observation_space.shape,
              self.env.action_space.n))

    def create_model(self):
        self.actor_net = Actor(
            np.prod(self.env.observation_space.shape), self.env.action_space.n)
        self.critic_net = Critic(
            np.prod(self.env.observation_space.shape), self.env.action_space.n)
        self.actor_optimizer = optim.Adam(self.actor_net.parameters(), 1e-3)
        self.critic_net_optimizer = optim.Adam(
            self.critic_net.parameters(), 3e-3)

    def select_action(self, state):
        state = torch.from_numpy(state).float().unsqueeze(0)
        with torch.no_grad():
            action_prob = self.actor_net(state)
        c = Categorical(action_prob)
        action = c.sample()
        return action.item(), action_prob[:, action.item()].item()

    def get_value(self, state):
        state = torch.from_numpy(state)
        with torch.no_grad():
            value = self.critic_net(state)
        return value.item()

    def store_transition(self, transition):
        self.buffer.append(transition)
        self.counter += 1

    def update(self, i_ep):
        state = torch.tensor([t.state for t in self.buffer], dtype=torch.float)
        action = torch.tensor(
            [t.action for t in self.buffer], dtype=torch.long).view(-1, 1)
        reward = [t.reward for t in self.buffer]
        old_action_log_prob = torch.tensor(
            [t.a_log_prob for t in self.buffer], dtype=torch.float).view(-1, 1)

        R = 0
        Gt = []
        for r in reward[::-1]:
            R = r + self.gamma * R
            Gt.insert(0, R)
        Gt = torch.tensor(Gt, dtype=torch.float)

        for i in range(self.ppo_update_time):
            for index in BatchSampler(SubsetRandomSampler(range(len(self.buffer))), self.batch_size, False):
                if self.training_step % 1000 == 0:
                    print('I_ep {} ，train {} times'.format(
                        i_ep, self.training_step))

                Gt_index = Gt[index].view(-1, 1)
                V = self.critic_net(state[index])
                delta = Gt_index - V
                advantage = delta.detach()

                action_prob = self.actor_net(state[index]).gather(
                    1, action[index])

                ratio = (action_prob/old_action_log_prob[index])
                surr1 = ratio * advantage
                surr2 = torch.clamp(ratio, 1 - self.clip_param,
                                    1 + self.clip_param) * advantage

                action_loss = -torch.min(surr1, surr2).mean()

                print("loss= {}, step = {}".format(
                    action_loss.item(), self.training_step))

                self.actor_optimizer.zero_grad()
                action_loss.backward()
                nn.utils.clip_grad_norm_(
                    self.actor_net.parameters(), self.max_grad_norm)
                self.actor_optimizer.step()

                value_loss = F.mse_loss(Gt_index, V)
                self.critic_net_optimizer.zero_grad()
                value_loss.backward()
                nn.utils.clip_grad_norm_(
                    self.critic_net.parameters(), self.max_grad_norm)
                self.critic_net_optimizer.step()
                self.training_step += 1

        del self.buffer[:]

    def save_model(self):
        torch.save(self.critic_net, self.checkpoint_path)

def main():
    Transition = namedtuple(
        'Transition', ['state', 'action',  'a_log_prob', 'reward', 'next_state'])

    agent = PPO()
    agent.init_hyp()
    agent.init_env()
    agent.create_model()
    for i_epoch in range(1000):
        state = agent.env.reset()
        if agent.render:
            agent.env.render()

        for t in count():
            action, action_prob = agent.select_action(state)
            next_state, reward, done, _ = agent.env.step(action)
            trans = Transition(state, action, action_prob, reward, next_state)
            if agent.render:
                agent.env.render()
            agent.store_transition(trans)
            state = next_state

            if done:
                if len(agent.buffer) >= agent.batch_size:
                    agent.update(i_epoch)
                    agent.save_model()
                break

# standalone function evaluate model on the environment
def evaluate():
    checkpoint_path = "single-agent/env-highway/ppo/checkpoints/ppo-highway.pth"
    model = torch.load(checkpoint_path)
    env = gym.make("highway-v0")
    state = env.reset()
    done = False

    while done != True:
        state = torch.from_numpy(state).float().unsqueeze(0)
        with torch.no_grad():
            action_prob = model(state)
        c = Categorical(action_prob)
        action = c.sample()
        state, reward, done, info = env.step(action.item())
        env.render()

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


text:
import numpy as np
from scipy.io import wavfile

TEST_FILE_NIST = (
    "/mnt/shared_storage/datasets/csr_1/11-10.1/wsj0/sd_tr_s/00f/00fo040b.wv2"
)
TEST_FILE_WAV = "./bss_speech_dataset/data/channels2_room16_mix.wav"

_nist_type_decoder = {
    "-s": str,
    "-i": int,
    "-f": float,
}

def _read_nist_decode_header_field(field):
    key, t, value = field.split(maxsplit=2)

    if len(t) == 2:
        n_chars = len(value)
    else:
        n_chars = int(t[2:])
    value = _nist_type_decoder[t[:2]](value[:n_chars])

    return key, value

def _read_nist_dtype(n_bytes, byte_format_str):

    if n_bytes == 2:
        dt = np.dtype(np.int16)
        if byte_format_str == "01":
            return dt
        elif byte_format_str == "10":
            return dt.newbyteorder(">")
        else:
            raise NotImplementedError(
                f"n_bytes={n_bytes} byte_format_str={byte_format_str}"
            )
    else:
        raise NotImplementedError(
            f"n_bytes={n_bytes} byte_format_str={byte_format_str}"
        )

def _read_nist_header(f):

    end_marker = "end_head"
    newline = b"\n"

    raw_header = b""

    # get read of first line
    c = ""
    while c != newline:
        c = f.read(1)
        raw_header += c

    # read number of bytes in header on second line
    n_bytes_header_str = b""
    c = ""
    while c != newline:
        c = f.read(1)
        n_bytes_header_str += c

    # concatenate to header
    raw_header += n_bytes_header_str

    # convert the byte string to an integer
    n_bytes_header = int(n_bytes_header_str)

    # read-in the rest of the header
    while len(raw_header) < n_bytes_header:
        raw_header += f.read(1)

    # convert to ASCII
    raw_header = raw_header.decode("ascii")

    lines = raw_header.split("\n")

    header = {}

    for line in lines[2:]:

        if line.startswith(end_marker):
            break
        else:
            key, value = _read_nist_decode_header_field(line)
            header[key] = value

    fmt_marker = lines[0][:-1]

    return fmt_marker, header

def read_audio_nist(filename):

    with open(filename, "rb") as f:
        fmt_marker, header = _read_nist_header(f)
        data = f.read()

    # sampling frequency
    fs = header["sample_rate"]

    # number of channels
    if "numChannels" in header and "channel_count" not in header:
        n_channels = header["numChannels"]
    elif "channel_count" in header:
        n_channels = header["channel_count"]

    if "channels_interleaved" in header:
        channels_interleaved = header["channels_interleaved"] == "TRUE"
    else:
        channels_interleaved = True

    # byte format
    n_bytes = header["sample_n_bytes"]
    byte_fmt_str = header["sample_byte_format"]
    dtype = _read_nist_dtype(n_bytes, byte_fmt_str)

    # total number of bytes
    n_samples = min(header["sample_count"], len(data)) // n_bytes

    n_samples = (n_samples // n_channels) * n_channels

    shape = (n_samples // n_channels, n_channels)
    if not channels_interleaved:
        order = "F"
    else:
        order = "C"

    audio = np.frombuffer(data, dtype=dtype, count=n_samples).reshape(
        shape, order=order
    )

    if n_channels == 1:
        audio = audio[:, 0]

    return fs, audio

reader = {
    "NIST": read_audio_nist,
    "RIFF": wavfile.read,
}

def read_audio(filename):

    with open(filename, "rb") as f:
        fmt = f.read(4).decode("ascii")

    return reader[fmt](filename)

if __name__ == "__main__":

    fs, audio = read_audio(TEST_FILE_NIST)
    # read_audio(TEST_FILE_WAV)

reason: import numpy,from scipy


text:
import math
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point, Polygon, LineString
import shapely.affinity

from Node import Node
from Person import Person
from A_star import a_star

from time import time

from CalculateWeight import calculateWeight
from A_star import checkValidPath

def main():
    # with open('test.txt','w') as file:
    #     file.write("==================\n")

    start = Node(position=[0,0])

    end = Node(position=[50,75])

    kevin = Person(x=20,y=25,velocity=0,direction=0)
    #john = Person(x=2,y=7,velocity=0,direction=0)
    #am = Person(x=80,y=20,velocity=0,direction=0)
    # mike = Person(x=10,y=-20,velocity=0,direction=0)

    adjacencyList, listOfShape = createAdjacencyList_new( nodes =[start,end], people=[kevin])

    for i,j in adjacencyList.items():
        i.cost += calculateWeight(end,i)

    # plotAllPath(adjacencyList, listOfShape)

    listOfPath = a_star( start, end ,adjacencyList, listOfShape)
    listOfPath.reverse()
    for i in listOfPath:
        print ("{} to {}".format(i[0],i[1]))
    plot(listOfPath,adjacencyList)

def plot(listOfPath,adjacencyList):
    x_values = []
    y_values = []

    for i in range(0,len(listOfPath)-1):
        plt.plot( [(listOfPath[i])[0],(listOfPath[i+1])[0]],[(listOfPath[i])[1],(listOfPath[i+1])[1]],'k-')

    for x,y in adjacencyList.items():
        x_values.append(x.position[0])
        y_values.append(x.position[1])

    # for p in y:
    #         plt.plot([x.position[0],p.position[0]],[x.position[1],p.position[1]],'k-')
    #         x_values.append(p.position[0])
    #         y_values.append(p.position[1])

    plt.scatter(x_values, y_values, s=5)
    plt.show()

def createAdjacencyList_new(nodes, people):
    listOfShape = []

    for p in people :
        nodes.extend(p.getKeyNodes())
        p.makeSocialZone()
        for circle in p.listOfCicle:
            listOfShape.append(circle)
            for node in circle:
                nodes.append(node)

    # for o in obstacles:
    #     listOfShape.append(obstacles.makeSocialZone)

    adjacencyList = {}
    for i in nodes:
        adjacencyList[i] = []
        for j in nodes :
            if ( i!=j ):
                adjacencyList[i].append(j)

    return adjacencyList, listOfShape

if __name__ == "__main__":
    main()


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

inter_min = 15
inter = '{}min'.format(inter_min)
DATA_PATH = 'langfang_data/'
SAFE_NAME = DATA_PATH + inter + '_data.h5'
ADJ_NAME = DATA_PATH + 'intersection.xlsx'
SAFE_ADJ_NAME = DATA_PATH + inter + '_adj_mx.pkl'
NUM_NODES = 65 if inter_min == 10 else 94
NUM_MONTH = 4
NUM_DIRECT = 4
CLEAN_THRES = 0.25
ADJ_THRES = 0.01

time_stamp = pd.date_range("2014-01-01 00:00","2014-0{}-30 23:59".format(NUM_MONTH),freq=inter)
TOTAL_TIME = time_stamp.size
flow = pd.DataFrame(0, columns=np.arange(1, NUM_NODES+1), index=time_stamp) # 2-D

listd = os.listdir(DATA_PATH+inter)

for each in listd:
    if 'F' in each:
        flow_node = pd.read_csv(os.path.join(DATA_PATH+inter, each),
                           sep=",", header=None)
        id_node, id_direct = each.split('_')
        id_node = int(id_node[2:])

        flow_node = flow_node.iloc[:,5:]
        assert flow_node.shape[1] == 24*60/inter_min
        flow_node = flow_node.values.flatten()
        flow_node = flow_node[:TOTAL_TIME].transpose()

        # Add to flow data
        flow[id_node] += flow_node

# Data cleaning
drop_cols = flow.columns[(flow == 0).sum() > CLEAN_THRES*flow.shape[1]]
flow.drop(drop_cols, axis = 1, inplace = True)
flow /= NUM_DIRECT

# Generate adj matrix
id_node_list = list(flow.columns)

def distance(origin, destination):
    lat1, lon1 = origin
    lat2, lon2 = destination
    radius = 6371 # km

    dlat = math.radians(lat2-lat1)
    dlon = math.radians(lon2-lon1)
    a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat1)) \
        * math.cos(math.radians(lat2)) * math.sin(dlon/2) * math.sin(dlon/2)
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
    d = radius * c

    return d

hdf = pd.HDFStore(SAFE_NAME)
hdf.put('data', flow, format='t')
hdf.close()

max_stars_repo_name: kangzf/DCRNN

max_issues_repo_name: kangzf/DCRNN

max_forks_repo_name: kangzf/DCRNN


text:
import numpy as np
from scipy.constants import pi
from numpy.fft import fftshift
from scipy.fftpack import fft, ifft
from six.moves import builtins
from cython_files.cython_integrand import *
import sys
assert_allclose = np.testing.assert_allclose
import numba
complex128 = numba.complex128

vectorize = numba.vectorize
autojit, jit = numba.autojit, numba.jit
cfunc = numba.cfunc
generated_jit = numba.generated_jit
guvectorize = numba.guvectorize

# Pass through the @profile decorator if line profiler (kernprof) is not in use
# Thanks Paul!
try:
    builtins.profile
except AttributeError:
    def profile(func):
        return func
from time import time
import pickle

#@profile
class Integrator(object):
    def __init__(self, int_fwm):
        if int_fwm.nm == 1:
            self.RK45mm = self.RK45CK_nm1
        elif int_fwm.nm == 2:
            self.RK45mm = self.RK45CK_nm2
        else:
            sys.exit('Too many modes!!')

        return None

    def RK45CK_nm1(self, dAdzmm, u1, dz, M1, M2, Q, tsh, dt, hf, w_tiled, gam_no_aeff):
        """
        Propagates the nonlinear operator for 1 step using a 5th order Runge
        Kutta method
        use: [A delta] = RK5mm(u1, dz)
        where u1 is the initial time vector
        hf is the Fourier transform of the Raman nonlinear response time
        dz is the step over which to propagate

        in output: A is new time vector
        delta is the norm of the maximum estimated error between a 5th
        order and a 4th order integration
        """
        (u1, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
        A1 = dz*dAdzmm(u1, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
        u2 = A2_temp(u1, A1)

        A2 = dz*dAdzmm(u1, Q, tsh, dt, hf, w_tiled,gam_no_aeff)

        u3 = A3_temp(u1, A1,A2)
        A3 = dz*dAdzmm(u1, Q, tsh, dt, hf, w_tiled,gam_no_aeff)

        u4 = A4_temp(u1, A1, A2, A3)
        A4 = dz*dAdzmm(u1, Q, tsh, dt, hf, w_tiled,gam_no_aeff)

        u5 = A5_temp(u1, A1, A2, A3, A4)
        A5 = dz*dAdzmm(u1, Q, tsh, dt, hf, w_tiled,gam_no_aeff)

        u6 = A6_temp(u1, A1, A2, A3, A4, A5)
        A6 = dz*dAdzmm(u1, Q, tsh, dt, hf, w_tiled,gam_no_aeff)

        A =  A_temp(u1, A1, A3, A4, A6) # Fifth order accuracy
        Afourth =  Afourth_temp(u1, A1, A3, A4,A5, A6) # Fourth order accuracy

        delta = np.linalg.norm(A - Afourth,2, axis = 1).max()
        return A, delta

    def RK45CK_nm2(self, dAdzmm, u1, dz, M1, M2,Q, tsh, dt, hf, w_tiled, gam_no_aeff):
        """
        Propagates the nonlinear operator for 1 step using a 5th order Runge
        Kutta method
        use: [A delta] = RK5mm(u1, dz)
        where u1 is the initial time vector
        hf is the Fourier transform of the Raman nonlinear response time
        dz is the step over which to propagate

        in output: A is new time vector
        delta is the norm of the maximum estimated error between a 5th
        order and a 4th order integration
        """

        A1 = dz*dAdzmm(u1,u1.conj(), M1, M2, Q, tsh, dt, hf, w_tiled, gam_no_aeff)
        u2 = A2_temp(u1, A1)

        A2 = dz*dAdzmm(u2,u2.conj(), M1, M2, Q, tsh, dt, hf, w_tiled, gam_no_aeff)

        u3 = A3_temp(u1, A1,A2)
        A3 = dz*dAdzmm(u3,u3.conj(), M1, M2, Q, tsh, dt, hf, w_tiled, gam_no_aeff)

        u4 = A4_temp(u1, A1, A2, A3)
        A4 = dz*dAdzmm(u4,u4.conj(), M1, M2, Q, tsh, dt, hf, w_tiled, gam_no_aeff)

        u5 = A5_temp(u1, A1, A2, A3, A4)
        A5 = dz*dAdzmm(u5,u5.conj(), M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff)

        u6 = A6_temp(u1, A1, A2, A3, A4, A5)
        A6 = dz*dAdzmm(u6,u6.conj(), M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff)

        A =  A_temp(u1, A1, A3, A4, A6) # Fifth order accuracy
        Afourth =  Afourth_temp(u1, A1, A3, A4,A5, A6) # Fourth order accuracy

        delta = np.linalg.norm(A - Afourth,2, axis = 1).max()
        return A, delta

trgt = 'cpu'
#trgt = 'parallel'
#trgt = 'cuda'

@jit(nopython=True,nogil = True)
def Afourth_temp(u1, A1, A3, A4, A5, A6):
    return u1 + (2825./27648)*A1 + (18575./48384)*A3 + (13525./55296) * \
        A4 + (277./14336)*A5 + (1./4)*A6

@jit(nopython=True,nogil = True)
def A_temp(u1, A1, A3, A4, A6):
    return u1 + (37./378)*A1 + (250./621)*A3 + (125./594) * \
        A4 + (512./1771)*A6

@jit(nopython=True,nogil = True)
def A2_temp(u1, A1):
    return u1 + (1./5)*A1

@jit(nopython=True,nogil = True)
def A3_temp(u1, A1, A2):
    return u1 + (3./40)*A1 + (9./40)*A2

@jit(nopython=True,nogil = True)
def A4_temp(u1, A1, A2, A3):
    return u1 + (3./10)*A1 - (9./10)*A2 + (6./5)*A3

@jit(nopython=True,nogil = True)
def A5_temp(u1, A1, A2, A3, A4):
    return u1 - (11./54)*A1 + (5./2)*A2 - (70./27)*A3 + (35./27)*A4

@jit(nopython=True,nogil = True)
def A6_temp(u1, A1, A2, A3, A4, A5):
    return u1 + (1631./55296)*A1 + (175./512)*A2 + (575./13824)*A3 +\
                   (44275./110592)*A4 + (253./4096)*A5

"""--------------------------Two modes-------------------------------------"""
#@jit(nogil = True)
def dAdzmm_roff_s0_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff):
    """
    calculates the nonlinear operator for a given field u0
    use: dA = dAdzmm(u0)
    """
    M3 = uabs_nm2(u0,u0_conj,M2)
    N = nonlin_kerr_nm2(M1, Q, u0, M3)
    N *= gam_no_aeff
    return N

#@jit(nogil = True)
def dAdzmm_roff_s1_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff):
    """
    calculates the nonlinear operator for a given field u0
    use: dA = dAdzmm(u0)
    """
    M3 = uabs_nm2(u0,u0_conj,M2)
    N = nonlin_kerr_nm2(M1, Q, u0, M3)
    N = gam_no_aeff * (N + tsh*ifft(w_tiled * fft(N)))
    return N

def dAdzmm_ron_s0_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled, gam_no_aeff):
    """
    calculates the nonlinear operator for a given field u0
    use: dA = dAdzmm(u0)
    """
    M3 = uabs_nm2(u0,u0_conj,M2)
    M4 = dt*fftshift(ifft(fft(M3)*hf), axes = -1) # creates matrix M4
    N = nonlin_ram_nm2(M1, Q, u0, M3, M4)
    N *= gam_no_aeff
    return N

def dAdzmm_ron_s1_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff):
    """
    calculates the nonlinear operator for a given field u0
    use: dA = dAdzmm(u0)
    """
    M3 = uabs_nm2(u0,u0_conj,M2)
    M4 = dt*fftshift(ifft(multi(fft(M3),hf)), axes = -1) # creates matrix M4
    N = nonlin_ram_nm2(M1, Q, u0, M3, M4)

    N = gam_no_aeff * (N + tsh*ifft(multi(w_tiled,fft(N))))
    return N

@guvectorize(['void(complex128[:,:],complex128[:,:], int64[:,:], complex128[:,:])'],\
                 '(n,m),(n,m),(o,l)->(l,m)',target = trgt)
def uabs_nm2(u0,u0_conj,M2,M3):
    for ii in range(M2.shape[1]):
        M3[ii,:] = u0[M2[0,ii],:]*u0_conj[M2[1,ii],:]

@guvectorize(['void(int64[:,:], complex128[:,:], complex128[:,:],\
            complex128[:,:], complex128[:,:], complex128[:,:])'],\
            '(w,a),(i,a),(m,n),(l,n),(l,n)->(m,n)',target = trgt)
def nonlin_ram_nm2(M1, Q, u0, M3, M4, N):
    N[:,:] = 0
    for ii in range(M1.shape[1]):
        N[M1[0,ii],:] += u0[M1[1,ii],:]*(0.82*(2*Q[0,ii] + Q[1,ii]) \
                                *M3[M1[4,ii],:] + \
                               0.54*Q[0,ii]*M4[M1[4,ii],:])

@guvectorize(['void(int64[:,:], complex128[:,:], complex128[:,:],\
            complex128[:,:], complex128[:,:])'],\
            '(w,a),(i,a),(m,n),(l,n)->(m,n)',target = trgt)
def nonlin_kerr_nm2(M1, Q, u0, M3, N):
    N[:,:] = 0
    for ii in range(M1.shape[1]):
        N[M1[0,ii],:] += 0.82*(2*Q[0,ii] + Q[1,ii]) \
                                *u0[M1[1,ii],:]*M3[M1[4,ii],:]

"""------------------------------------------------------------------------"""

"""-----------------------------One mode-----------------------------------"""
#@jit(nogil = True)
def dAdzmm_roff_s0_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff):
    """
    calculates the nonlinear operator for a given field u0
    use: dA = dAdzmm(u0)
    """
    M3 = u0.real**2 + u0.imag**2
    #M3 = uabs_nm1(u0.real, u0.imag)
    N = nonlin_kerr_nm1(Q, u0, M3)
    N *= gam_no_aeff
    return N

#@jit(nogil = True)
def dAdzmm_roff_s1_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff):
    """
    calculates the nonlinear operator for a given field u0
    use: dA = dAdzmm(u0)
    """
    M3 = u0.real**2 + u0.imag**2
    #M3 = uabs_nm1(u0.real, u0.imag)
    N = nonlin_kerr_nm1(Q, u0, M3)
    N = gam_no_aeff * (N + tsh*ifft(w_tiled * fft(N)))
    return N

def dAdzmm_ron_s0_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff):
    """
    calculates the nonlinear operator for a given field u0
    use: dA = dAdzmm(u0)
    """
    M3 = u0.real**2 + u0.imag**2
    #M3 = uabs_nm1(u0.real, u0.imag)
    M4 = dt*fftshift(ifft(fft(M3)*hf), axes = -1) # creates matrix M4
    N = nonlin_ram_nm1(Q, u0, M3, M4)
    N *= gam_no_aeff
    return N

def dAdzmm_ron_s1_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff):
    """
    calculates the nonlinear operator for a given field u0
    use: dA = dAdzmm(u0)
    """
    M3 = u0.real**2 + u0.imag**2

    M3 = uabs_nm1(u0.real, u0.imag)
    M4 = dt *fftshift(ifft(fft(M3)*hf), axes = -1)
    N = nonlin_ram_nm1(Q, u0, M3, M4)
    N = gam_no_aeff * (N + tsh*ifft(w_tiled * fft(N)))
    return N

@vectorize(['float64(float64, float64)'], target=trgt)
def uabs_nm1(u0r, u0i):
    return u0r**2 + u0i**2

@vectorize(['complex128(complex128, complex128, float64, complex128)'], target=trgt)
def nonlin_ram_nm1(Q, u0, M3, temp):
    return Q*u0*(0.82*M3 + 0.18*temp)

@vectorize(['complex128(complex128, complex128, float64)'], target=trgt)
def nonlin_kerr_nm1(Q, u0, M3):
    return 0.82*Q*u0*M3
"""------------------------------------------------------------------------"""

@jit(nopython=True,nogil = True)
def multi(a,b):
    return a * b

class Integrand(object):
    def __init__(self,nm,ram, ss, cython = True, timing = False):
        print('number of modes: ', nm)
        if nm == 2:
            if cython:
                if ss == 0 and ram == 'off':
                    self.dAdzmm = dAdzmm_roff_s0_cython_nm2
                elif ss == 0 and ram == 'on':
                    self.dAdzmm = dAdzmm_ron_s0_cython_nm2
                elif ss == 1 and ram == 'off':
                    self.dAdzmm = dAdzmm_roff_s1_cython_nm2
                else:
                    self.dAdzmm = dAdzmm_ron_s1_cython_nm2
            else:
                if ss == 0 and ram == 'off':
                    self.dAdzmm = dAdzmm_roff_s0_nm2
                elif ss == 0 and ram == 'on':
                    self.dAdzmm = dAdzmm_ron_s0_nm2
                elif ss == 1 and ram == 'off':
                    self.dAdzmm = dAdzmm_roff_s1_nm2
                else:
                    self.dAdzmm = dAdzmm_ron_s1_nm2
            if timing:
                self.dAdzmm = self.timer_nm2
        elif nm ==1:
            if cython:
                if ss == 0 and ram == 'off':
                    self.dAdzmm = dAdzmm_roff_s0_cython_nm1
                elif ss == 0 and ram == 'on':
                    self.dAdzmm = dAdzmm_ron_s0_cython_nm1
                elif ss == 1 and ram == 'off':
                    self.dAdzmm = dAdzmm_roff_s1_cython_nm1
                else:
                    self.dAdzmm = dAdzmm_ron_s1_cython_nm1
            else:
                if ss == 0 and ram == 'off':
                    self.dAdzmm = dAdzmm_roff_s0_nm1
                elif ss == 0 and ram == 'on':
                    self.dAdzmm = dAdzmm_ron_s0_nm1
                elif ss == 1 and ram == 'off':
                    self.dAdzmm = dAdzmm_roff_s1_nm1
                else:
                    self.dAdzmm = dAdzmm_ron_s1_nm1
            if timing:
                self.dAdzmm = self.timer_nm1
        else:
            sys.exit('Too many modes!!!')
    def timer_nm2(self,u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff):
        """
        Times the functions of python, cython etc.
        """
        dt1, dt2, dt3, dt4, dt5, dt6, dt7, dt8 = [], [], [], [],\
                                                [], [], [], []

        NN = 100
        for i in range(NN):
            '------No ram, no ss--------'
            t = time()
            N1 = dAdzmm_roff_s0_cython_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt1.append(time() - t)

            t = time()
            N2 = dAdzmm_roff_s0_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt2.append(time() - t)
            assert_allclose(N1, N2)

            '------ ram, no ss--------'
            t = time()
            N1 = dAdzmm_ron_s0_cython_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt3.append(time() - t)

            t = time()
            N2 = dAdzmm_ron_s0_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt4.append(time() - t)
            assert_allclose(N1, N2)

            '------ no ram, ss--------'
            t = time()
            N1 = dAdzmm_roff_s1_cython_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt5.append(time() - t)

            t = time()
            N2 = dAdzmm_roff_s1_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt6.append(time() - t)
            assert_allclose(N1, N2)

            '------ ram, ss--------'
            t = time()
            N1 = dAdzmm_ron_s1_cython_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt7.append(time() - t)

            t = time()
            N2 = dAdzmm_ron_s1_nm2(u0,u0_conj, M1, M2, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt8.append(time() - t)
            assert_allclose(N1, N2)

        print('cython_ram(off)_s0: {} +/- {}'.format(np.average(dt1),np.std(dt1)))
        print('python_ram(off)_s0: {} +/- {}'.format(np.average(dt2),np.std(dt2)))
        print('Cython is {} times faster'.format(np.average(dt2)/np.average(dt1)))
        print('--------------------------------------------------------')
        print('cython_ram(on)_s0: {} +/- {}'.format(np.average(dt3),np.std(dt3)))
        print('python_ram(on)_s0: {} +/- {}'.format(np.average(dt4),np.std(dt4)))
        print('Cython is {} times faster'.format(np.average(dt4)/np.average(dt3)))
        print('--------------------------------------------------------')

        print('cython_ram(off)_s1: {} +/- {}'.format(np.average(dt5),np.std(dt5)))
        print('python_ram(off)_s1: {} +/- {}'.format(np.average(dt6),np.std(dt6)))
        print('Cython is {} times faster'.format(np.average(dt6)/np.average(dt5)))
        print('--------------------------------------------------------')

        print('cython_ram(on)_s1: {} +/- {}'.format(np.average(dt7),np.std(dt7)))
        print('python_ram(on)_s1: {} +/- {}'.format(np.average(dt8),np.std(dt8)))
        print('Cython is {} times faster'.format(np.average(dt8)/np.average(dt7)))
        print('--------------------------------------------------------')
        sys.exit()
        return N

    def timer_nm1(self,u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff):
        """
        Times the functions of python, cython etc.
        """
        dt1, dt2, dt3, dt4, dt5, dt6, dt7, dt8 = [], [], [], [],\
                                                [], [], [], []

        NN = 100
        for i in range(NN):

            '------No ram, no ss--------'
            t = time()
            N1 = dAdzmm_roff_s0_cython_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt1.append(time() - t)

            t = time()
            N2 = dAdzmm_roff_s0_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt2.append(time() - t)
            assert_allclose(N1, N2)

            '------ ram, no ss--------'
            t = time()
            N1 = dAdzmm_ron_s0_cython_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt3.append(time() - t)

            t = time()
            N2 = dAdzmm_ron_s0_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt4.append(time() - t)
            assert_allclose(N1, N2)

            '------ no ram, ss--------'
            t = time()
            N1 = dAdzmm_roff_s1_cython_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt5.append(time() - t)

            t = time()
            N2 = dAdzmm_roff_s1_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt6.append(time() - t)
            assert_allclose(N1, N2)

            '------ ram, ss--------'
            t = time()
            N1 = dAdzmm_ron_s1_cython_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt7.append(time() - t)

            t = time()
            N2 = dAdzmm_ron_s1_nm1(u0, Q, tsh, dt, hf, w_tiled,gam_no_aeff)
            dt8.append(time() - t)
            assert_allclose(N1, N2)

        print('cython_ram(off)_s0: {} +/- {}'.format(np.average(dt1),np.std(dt1)))
        print('python_ram(off)_s0: {} +/- {}'.format(np.average(dt2),np.std(dt2)))
        print('Cython is {} times faster'.format(np.average(dt2)/np.average(dt1)))
        print('--------------------------------------------------------')
        print('cython_ram(on)_s0: {} +/- {}'.format(np.average(dt3),np.std(dt3)))
        print('python_ram(on)_s0: {} +/- {}'.format(np.average(dt4),np.std(dt4)))
        print('Cython is {} times faster'.format(np.average(dt4)/np.average(dt3)))
        print('--------------------------------------------------------')

        print('cython_ram(off)_s1: {} +/- {}'.format(np.average(dt5),np.std(dt5)))
        print('python_ram(off)_s1: {} +/- {}'.format(np.average(dt6),np.std(dt6)))
        print('Cython is {} times faster'.format(np.average(dt6)/np.average(dt5)))
        print('--------------------------------------------------------')

        print('cython_ram(on)_s1: {} +/- {}'.format(np.average(dt7),np.std(dt7)))
        print('python_ram(on)_s1: {} +/- {}'.format(np.average(dt8),np.std(dt8)))
        print('Cython is {} times faster'.format(np.average(dt8)/np.average(dt7)))
        print('--------------------------------------------------------')
        sys.exit()
        return N

max_stars_repo_name: ibegleris/w-fopo

max_issues_repo_name: ibegleris/w-fopo

max_forks_repo_name: ibegleris/w-fopo

reason: import numpy,from numpy,from scipy,import numba


text:
import copy as cp
from inspect import signature

import numpy as np

from skmultiflow.core import BaseSKMObject, ClassifierMixin, MetaEstimatorMixin
from sklearn.tree import DecisionTreeClassifier

import warnings

def BatchIncremental(nominal_attributes=None):     # pragma: no cover
    warnings.warn(
        "'BatchIncremental' has been renamed to 'BatchIncrementalClassifier' in v0.5.0.\n"
        "The old name will be removed in v0.7.0",
        category=FutureWarning)
    return BatchIncrementalClassifier(nominal_attributes=nominal_attributes)

class BatchIncrementalClassifier(BaseSKMObject, ClassifierMixin, MetaEstimatorMixin):
    """ Batch Incremental ensemble classifier.

    This is a wrapper that allows the application of any batch model to a
    stream by incrementally building an ensemble of instances of the batch model.
    A window of examples is collected, then used to train a new model, which is
    added to the ensemble. A maximum number of models ensures memory use is finite
    (the oldest model is deleted when this number is exceeded).

    Parameters
    ----------
    base_estimator: skmultiflow.core.BaseSKMObject or sklearn.BaseEstimator
        (default=DecisionTreeClassifier) Each member of the ensemble is
        an instance of the base estimator

    window_size: int (default=100)
        The size of the training window (batch), in other words,
        how many instances are kept for training.

    n_estimators: int (default=100)
        Number of estimators in the ensemble.

    Notes
    -----
    Not yet multi-label capable.

    Examples
    --------
    >>> # Imports
    >>> from skmultiflow.data import SEAGenerator
    >>> from skmultiflow.meta import BatchIncrementalClassifier
    >>>
    >>> # Setup a data stream
    >>> stream = SEAGenerator(random_state=1)
    >>>
    >>> # Pre-training the classifier with 200 samples
    >>> X, y = stream.next_sample(200)
    >>> batch_incremental_cfier = BatchIncrementalClassifier()
    >>> batch_incremental_cfier.partial_fit(X, y)
    >>>
    >>> # Preparing the processing of 5000 samples and correct prediction count
    >>> n_samples = 0
    >>> correct_cnt = 0
    >>> while n_samples < 5000 and stream.has_more_samples():
    >>>     X, y = stream.next_sample()
    >>>     y_pred = batch_incremental_cfier.predict(X)
    >>>     if y[0] == y_pred[0]:
    >>>         correct_cnt += 1
    >>>     batch_incremental_cfier.partial_fit(X, y)
    >>>     n_samples += 1
    >>>
    >>> # Display results
    >>> print('Batch Incremental ensemble classifier example')
    >>> print('{} samples analyzed'.format(n_samples))
    >>> print('Performance: {}'.format(correct_cnt / n_samples))
    """

    def __init__(self, base_estimator=DecisionTreeClassifier(), window_size=100, n_estimators=100):
        self.window_size = window_size
        self.n_estimators = n_estimators
        self.base_estimator = base_estimator
        # The ensemble
        self.ensemble = []
        self.i = -1
        self.X_batch = None
        self.y_batch = None
        self.sample_weight = None

    def partial_fit(self, X, y=None, classes=None, sample_weight=None):
        """ Partially (incrementally) fit the model.

        Parameters
        ----------
        X : numpy.ndarray of shape (n_samples, n_features)
            The features to train the model.

        y: numpy.ndarray of shape (n_samples)
            An array-like with the labels of all samples in X.

        classes: Not used (default=None)

        sample_weight: numpy.ndarray of shape (n_samples), optional (default=None)
            Samples weight. If not provided, uniform weights are assumed.

        Returns
        -------
            self

        """
        N, D = X.shape

        if self.i < 0:
            # No models yet -- initialize
            self.X_batch = np.zeros((self.window_size, D))
            self.y_batch = np.zeros(self.window_size)
            self.sample_weight = np.zeros(self.window_size)
            self.i = 0

        for n in range(N):
            # For each instance ...
            # TODO not very pythonic at the moment
            self.X_batch[self.i] = X[n]
            self.y_batch[self.i] = y[n]
            self.sample_weight[self.i] = sample_weight[n] if sample_weight else 1.0

            self.i = self.i + 1
            if self.i == self.window_size:
                # Get rid of the oldest model
                if len(self.ensemble) >= self.n_estimators:
                    self.ensemble.pop(0)
                # A new model
                h = cp.deepcopy(self.base_estimator)
                # Train it
                if 'sample_weight' in signature(h.fit).parameters:
                    h.fit(X=self.X_batch, y=self.y_batch.astype(int), sample_weight=sample_weight)
                else:
                    h.fit(X=self.X_batch, y=self.y_batch.astype(int))
                # Add it
                self.ensemble.append(h)
                # Reset the window
                self.i = 0

        return self

    def predict_proba(self, X):
        """ Estimates the probability of each sample in X belonging to each of the class-labels.

        Parameters
        ----------
        X : numpy.ndarray of shape (n_samples, n_features)
            The matrix of samples one wants to predict the class probabilities for.

        Returns
        -------
        A numpy.ndarray of shape (n_samples, n_labels), in which each outer entry is
        associated with the X entry of the same index. And where the list in index [i] contains
        len(self.target_values) elements, each of which represents the probability that
        the i-th sample of X belongs to a certain class-label.

        """
        N, D = X.shape
        votes = np.zeros(N)
        if len(self.ensemble) <= 0:
            # No models yet, just predict zeros
            return votes
        for h_i in self.ensemble:
            # Add vote (normalized by number of models)
            votes = votes + 1. / len(self.ensemble) * h_i.predict(X)
        return votes

    def predict(self, X):
        """ Predict classes for the passed data.

        Parameters
        ----------
        X : numpy.ndarray of shape (n_samples, n_features)
            The set of data samples to predict the labels for.

        Returns
        -------
        A numpy.ndarray with all the predictions for the samples in X.

        """
        votes = self.predict_proba(X)
        # Suppose a threshold of 0.5
        return (votes >= 0.5) * 1.

    def reset(self):
        """ Resets the estimator to its initial state.

        Returns
        -------
            self

        """
        self.ensemble = []
        self.i = -1
        self.X_batch = None
        self.y_batch = None
        return self

max_stars_repo_name: denisesato/scikit-multiflow

max_issues_repo_name: denisesato/scikit-multiflow

max_forks_repo_name: denisesato/scikit-multiflow


text:
from sklearn.mixture import GaussianMixture
from gmm_tf_debug import GaussianMixture as GMM
import numpy as np
import time
from sklearn.metrics import homogeneity_completeness_v_measure

np.random.seed(42)

n_feature = 32
n_components = 5
sampe_size = [100, 200, 300, 200, 100]
scale = 1.1

X = []

for i in range(n_components):
    A = np.random.rand(n_feature, n_feature)
    B = np.dot(A, A.transpose())
    X.append(np.random.multivariate_normal(mean = np.random.normal(scale = scale, size = n_feature), cov = B, size = sampe_size[i]))

X = np.concatenate(X)

Y = np.repeat(range(n_components), sampe_size)

it = time.time()
gmm_sklearn = GaussianMixture(n_components=n_components, random_state=0, n_init = 50, tol=1e-3).fit(X)
ft = time.time()
print('sklearn time:', ft - it)
print('sklearn accuracy:', homogeneity_completeness_v_measure(Y, gmm_sklearn.predict(X)))

gmm_tf = GMM(n_components=n_components, random_state=0, n_init = 50, eps=1e-3, batch=800)
it2 = time.time()
gmm_tf.fit(X)
ft2 = time.time()
print('TF time:', ft2 - it2)
print('TF accuracy:', homogeneity_completeness_v_measure(Y, gmm_tf.predict(X)))


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

if __name__ == "__main__":

    # # Initialization

    # matplotlib.rcParams.update({'font.size': 12, 'text.usetex': True})

    # gParams

    (Lx, Ly, Lz) = (21, 21, 21)
    (dx, dy, dz) = (0.375, 0.375, 0.375)

    # (Lx, Ly, Lz) = (21, 21, 21)
    # (dx, dy, dz) = (0.25, 0.25, 0.25)

    NGridPoints_cart = (1 + 2 * Lx / dx) * (1 + 2 * Ly / dy) * (1 + 2 * Lz / dz)

    aIBi = -2
    datapath = '/media/kis/Storage/Dropbox/VariationalResearch/HarvardOdyssey/genPol_data/NGridPoints_{:.2E}/spherical'.format(NGridPoints_cart)
    outputdatapath = '/media/kis/Storage/Dropbox/VariationalResearch/DataAnalysis/newdata/NGridPoints_{:.2E}/dynamics_spherical/aIBi_{:.2f}'.format(NGridPoints_cart, aIBi)

    # datapath = '/home/kis/Dropbox/VariationalResearch/HarvardOdyssey/genPol_data/NGridPoints_{:.2E}/spherical'.format(NGridPoints_cart)
    # outputdatapath = '/home/kis/Dropbox/VariationalResearch/DataAnalysis/newdata/NGridPoints_{:.2E}/dynamics_spherical/aIBi_{:.2f}'.format(NGridPoints_cart, aIBi)

    if os.path.isdir(outputdatapath) is False:
        os.mkdir(outputdatapath)

    # # Individual Datasets

    for ind, filename in enumerate(os.listdir(datapath)):
        if filename == 'quench_Dataset_sph.nc':
            continue
        ds = xr.open_dataset(datapath + '/' + filename)

        if(ds.attrs['aIBi'] != aIBi):
            continue

        print(filename)
        P = ds.attrs['P']
        tGrid = ds.coords['t'].values
        PVec = P * np.ones(tGrid.size)
        CSAmp = (ds['Real_CSAmp'] + 1j * ds['Imag_CSAmp']).values
        CSAmp = CSAmp.reshape(CSAmp.size)

        # generates data file with columns representing P, t, Phase, Phonon Momentum, !Momentum Dispersion, Phonon Number, Re(Dynamical Overlap), Im(Dynamical Overlap)
        data = np.concatenate((aIBiVec[:, np.newaxis], PVec[:, np.newaxis], tGrid[:, np.newaxis], Real_DynOv_Vec[:, np.newaxis], Imag_DynOv_Vec[:, np.newaxis]), axis=1)
        np.savetxt(outputdatapath + '/quench_P_%.2f.dat' % P, data)

    # # steady state

    # aIBi = -1
    # datapath = '/home/kis/Dropbox/VariationalResearch/HarvardOdyssey/genPol_data/NGridPoints_{:.2E}/steadystate_spherical/quench_st_Dataset_sph.nc'.format(NGridPoints_cart)
    # outputdatapath = '/home/kis/Dropbox/VariationalResearch/DataAnalysis/newdata/NGridPoints_{:.2E}/steadystate_spherical'.format(NGridPoints_cart)

    # ds = xr.open_dataset(datapath)
    # dsa = ds.sel(aIBi=aIBi)

    # data = np.concatenate((dsa.coords['P'].values[:, np.newaxis], dsa['PB'].values[:, np.newaxis], dsa['NB'].values[:, np.newaxis], dsa['Z_factor'].values[:, np.newaxis], dsa['Energy'].values[:, np.newaxis], dsa['effMass'].values[:, np.newaxis], dsa['Pcrit'].values[:, np.newaxis]), axis=1)
    # np.savetxt(outputdatapath + '/quench_aIBi_{:.2f}.dat'.format(aIBi), data)

max_stars_repo_name: kseetharam/genPolaron

max_issues_repo_name: kseetharam/genPolaron

max_forks_repo_name: kseetharam/genPolaron


text:
import colorlover as cl
import plotly.graph_objs as go
import plotly.express as px
import numpy as np
import utils.bmex as bmex
from sklearn import metrics
from pickle import dump, load
import tensorflow as tf
from tensorflow.keras import backend as K

def serve_prediction_plot(
    model, X_train, X_test, y_train, y_test, Z, xx, yy, mesh_step, threshold
):
    # Get train and test score from model
    y_pred_train = (model.decision_function(X_train) > threshold).astype(int)
    y_pred_test = (model.decision_function(X_test) > threshold).astype(int)
    train_score = metrics.accuracy_score(y_true=y_train, y_pred=y_pred_train)
    test_score = metrics.accuracy_score(y_true=y_test, y_pred=y_pred_test)

    # Compute threshold
    scaled_threshold = threshold * (Z.max() - Z.min()) + Z.min()
    range = max(abs(scaled_threshold - Z.min()), abs(scaled_threshold - Z.max()))

    # Colorscale
    bright_cscale = [[0, "#ff3700"], [1, "#0b8bff"]]
    cscale = [
        [0.0000000, "#ff744c"],
        [0.1428571, "#ff916d"],
        [0.2857143, "#ffc0a8"],
        [0.4285714, "#ffe7dc"],
        [0.5714286, "#e5fcff"],
        [0.7142857, "#c8feff"],
        [0.8571429, "#9af8ff"],
        [1.0000000, "#20e6ff"],
    ]

    # Create the plot
    # Plot the prediction contour of the SVM
    trace0 = go.Contour(
        x=np.arange(xx.min(), xx.max(), mesh_step),
        y=np.arange(yy.min(), yy.max(), mesh_step),
        z=Z.reshape(xx.shape),
        zmin=scaled_threshold - range,
        zmax=scaled_threshold + range,
        hoverinfo="none",
        showscale=False,
        contours=dict(showlines=False),
        colorscale=cscale,
        opacity=0.9,
    )

    # Plot the threshold
    trace1 = go.Contour(
        x=np.arange(xx.min(), xx.max(), mesh_step),
        y=np.arange(yy.min(), yy.max(), mesh_step),
        z=Z.reshape(xx.shape),
        showscale=False,
        hoverinfo="none",
        contours=dict(
            showlines=False, type="constraint", operation="=", value=scaled_threshold
        ),
        name=f"Threshold ({scaled_threshold:.3f})",
        line=dict(color="#708090"),
    )

    # Plot Training Data
    trace2 = go.Scatter(
        x=X_train[:, 0],
        y=X_train[:, 1],
        mode="markers",
        name=f"Training Data (accuracy={train_score:.3f})",
        marker=dict(size=10, color=y_train, colorscale=bright_cscale),
    )

    # Plot Test Data
    trace3 = go.Scatter(
        x=X_test[:, 0],
        y=X_test[:, 1],
        mode="markers",
        name=f"Test Data (accuracy={test_score:.3f})",
        marker=dict(
            size=10, symbol="triangle-up", color=y_test, colorscale=bright_cscale
        ),
    )

    layout = go.Layout(
        xaxis=dict(ticks="", showticklabels=False, showgrid=False, zeroline=False),
        yaxis=dict(ticks="", showticklabels=False, showgrid=False, zeroline=False),
        hovermode="closest",
        legend=dict(x=0, y=-0.01, orientation="h"),
        margin=dict(l=0, r=0, t=0, b=0),
        plot_bgcolor="#282b38",
        paper_bgcolor="#282b38",
        font={"color": "#a5b1cd"},
    )

    data = [trace0, trace1, trace2, trace3]
    figure = go.Figure(data=data, layout=layout)

    return figure

def isotope_chain_go(Z, NRange, model, axis_label, func):

    output = []

    for i in NRange:
        output.append(func(i, Z, model))

    trace0 = go.Scatter(
        x=NRange, y=output, mode="lines", name="Test Data", marker={"color": "#13c6e9"}
    )

    layout = go.Layout(
        #title=f"ROC Curve (AUC = {auc_score:.3f})",
        title=f"Isotopic Chain",
        xaxis=dict(title="Neutrons", gridcolor="#2f3445"),
        yaxis=dict(title=axis_label, gridcolor="#2f3445"),
        #legend=dict(x=0, y=1.05, orientation="h"),
        #margin=dict(l=100, r=10, t=25, b=40),
        plot_bgcolor="#282b38",
        paper_bgcolor="#282b38",
        font={"color": "#a5b1cd"},
    )

    data = [trace0]
    figure = go.Figure(data=data, layout=layout)

    return figure

def isotope_chain(Z, NRange, model, axis_label, func):

    output = []

    layout = go.Layout(
        #title=f"ROC Curve (AUC = {auc_score:.3f})",
        title=f"Isotopic Chain",
        xaxis=dict(title="Neutrons", gridcolor="#2f3445",title_font_size=14),
        yaxis=dict(title=axis_label, gridcolor="#2f3445",title_font_size=14),
        #legend=dict(x=0, y=1.05, orientation="h"),
        #margin=dict(l=100, r=10, t=25, b=40),
        plot_bgcolor="#282b38",
        paper_bgcolor="#282b38",
        font={"color": "#a5b1cd", "size": 14},
    )

    N = bmex.df[(bmex.df["Z"]==Z) & (bmex.df["Model"]==model) & (bmex.df["N"] >= NRange[0]) & (bmex.df["N"] <= NRange[1])]["N"]
    for i in N:
        output.append(func(i, Z, model))

    trace0 = go.Scatter(
        x=N, y=output, mode="lines+markers", name="Test Data", marker=\
            {
                "color": "#13c6e9",
                #"size": 20,
            }
    )

    #figure = px.line(x=Z, y=output, markers=True)
    data = [trace0]
    figure = go.Figure(data=data, layout=layout)
    figure.update_xaxes(title_font_size=20)
    figure.update_yaxes(title_font_size=20)
    figure.update_layout(title_font_size=24)
    return figure

def R2(y_true, y_pred):
    SS_res =  K.sum(K.square( y_true-y_pred ))
    SS_tot = K.sum(K.square( y_true - K.mean(y_true) ) )
    return ( 1 - SS_res/(SS_tot + K.epsilon()) )

def PESpredict_single(A, Z, Q20, Q30, model, scaler):

    X = np.zeros(4)
    X[0] = A
    X[1] = Z
    X[2] = Q20
    X[3] = Q30
    xx = scaler.transform(X.reshape(1,-1))
    pes=model.predict(xx)

    return pes.T.squeeze()

def PESpredict(X, model, scaler):

    xx = scaler.transform(X)
    pes=model.predict(xx)

    return pes.T.squeeze()

def pesnet_surface(N, Z):

    model = tf.keras.models.load_model('utils/PESNet.model',custom_objects={'R2':R2})
    xscaler = load( open( 'utils/xscaler.pkl', "rb" ) )

    An = N+Z
    q20num = 121
    q30num = 21
    Q20 = np.linspace(0,250,q20num)
    Q30 = np.linspace(0,60,q30num)
    Q20v, Q30v = np.meshgrid(Q20,Q30, indexing='ij')
    Q20v = Q20v.flatten()
    Q30v = Q30v.flatten()
    Av = np.zeros(Q20v.shape)
    Zv = np.zeros(Q20v.shape)
    Av[:] = An
    Zv[:] = Z
    #xx = (Q20,Q30)
    #PESv = np.vectorize(PESpredict)
    #pes = PESv(A, Z, Q20v, Q30v, model, xscaler)
    X=(Av,Zv,Q20v,Q30v)
    X = np.asarray(X)
    X = X.T

    pes = np.reshape(PESpredict(X,model,xscaler),(q20num,q30num)).T

    layout = go.Layout(
        #title=f"ROC Curve (AUC = {auc_score:.3f})",
        title=f"PESnet Prediction",
        xaxis=dict(title="Q20", gridcolor="#2f3445",title_font_size=14),
        yaxis=dict(title="Q30", gridcolor="#2f3445",title_font_size=14),
        #legend=dict(x=0, y=1.05, orientation="h"),
        #margin=dict(l=100, r=10, t=25, b=40),
        plot_bgcolor="#282b38",
        paper_bgcolor="#282b38",
        font={"color": "#a5b1cd", "size": 28},
    )

    fig = go.Figure(data =
    go.Contour(
        z=pes,
        x=Q20, # horizontal axis
        y=Q30, # vertical axis
        colorscale="Spectral_r",
    ),layout=layout)
    return fig

def true_surface(df, N, Z):

    An = N+Z
    q20num = 126
    q30num = 21
    Q20 = np.linspace(0,250,q20num)
    Q30 = np.linspace(0,60,q30num)
    Q20v, Q30v = np.meshgrid(Q20,Q30, indexing='ij')
    Q20v = Q20v.flatten()
    Q30v = Q30v.flatten()
    Av = np.zeros(Q20v.shape)
    Zv = np.zeros(Q20v.shape)
    Av[:] = An
    Zv[:] = Z
    #xx = (Q20,Q30)
    #PESv = np.vectorize(PESpredict)
    #pes = PESv(A, Z, Q20v, Q30v, model, xscaler)
    PES_sub = df[(df['A'] == An) & (df['Z'] == Z)]

    A = PES_sub['A']
    Z = PES_sub['Z']
    #Q20 = PES_sub['Q20']
    #Q30 = PES_sub['Q30']
    E = PES_sub['HFB_cubic']

    layout = go.Layout(
        #title=f"ROC Curve (AUC = {auc_score:.3f})",
        title=f"True PES",
        xaxis=dict(title="Q20", gridcolor="#2f3445",title_font_size=14),
        yaxis=dict(title="Q30", gridcolor="#2f3445",title_font_size=14),
        #legend=dict(x=0, y=1.05, orientation="h"),
        #margin=dict(l=100, r=10, t=25, b=40),
        plot_bgcolor="#282b38",
        paper_bgcolor="#282b38",
        font={"color": "#a5b1cd", "size": 28},
    )

    fig = go.Figure(data =
    go.Contour(
        z=np.reshape(E.to_numpy(),(q20num,q30num)).T,
        x=Q20, # horizontal axis
        y=Q30, # vertical axis
        colorscale="Spectral_r",
    ),layout=layout)
    return fig

def isotope_chain_gp(Z, NRange, model, axis_label, func):

    output = []

    layout = go.Layout(
        #title=f"ROC Curve (AUC = {auc_score:.3f})",
        title=f"Isotopic Chain",
        xaxis=dict(title="Neutrons", gridcolor="#2f3445",title_font_size=14),
        yaxis=dict(title=axis_label, gridcolor="#2f3445",title_font_size=14),
        #legend=dict(x=0, y=1.05, orientation="h"),
        #margin=dict(l=100, r=10, t=25, b=40),
        plot_bgcolor="#282b38",
        paper_bgcolor="#282b38",
        font={"color": "#a5b1cd", "size": 14},
    )

    N = bmex.df[(bmex.df["Z"]==Z) & (bmex.df["Model"]=="Exp") & (bmex.df["N"] >= NRange[0]) & (bmex.df["N"] <= NRange[1])]["N"]
    plotFRDMGP,plotFRDMGPBandplus,plotFRDMGPBandminus = func(model[0],model[1],model[2],Z)

    trace0 = go.Scatter(
        x=N, y=plotFRDMGP, mode="lines+markers", name="Test Data", marker=\
            {
                "color": "#13c6e9",
                #"size": 20,
            }
    )

    #figure = px.line(x=Z, y=output, markers=True)
    data = [trace0]
    figure = go.Figure([
        go.Scatter(
            x=plotFRDMGP.T[0],
            y=plotFRDMGP.T[1],
            line=dict(color='rgb(100,100,80)'),
            mode='lines'
        ),
        go.Scatter(
            x=plotFRDMGP.T[0]+plotFRDMGP.T[0][::-1], # x, then x reversed
            y=plotFRDMGPBandplus.T[0]+plotFRDMGPBandminus.T[1][::-1], # upper, then lower reversed
            fill='toself',
            fillcolor='rgba(255,0,0)',
            line=dict(color='rgba(255,255,255,0)'),
            hoverinfo="skip",
            #showlegend=False
        )
    ], layout=layout)
    figure.update_xaxes(title_font_size=20)
    figure.update_yaxes(title_font_size=20)
    figure.update_layout(title_font_size=24)
    return figure

def isotone_chain_go(N, ZRange, model, axis_label, func):

    output = []

    for i in ZRange:
        output.append(func(N, i, model))

    trace0 = go.Scatter(
        x=ZRange, y=output, mode="lines", name="Test Data", marker={"color": "#13c6e9"}
    )

    layout = go.Layout(
        #title=f"ROC Curve (AUC = {auc_score:.3f})",
        title=f"Isotonic Chain",
        xaxis=dict(title="Protons", gridcolor="#2f3445"),
        yaxis=dict(title=axis_label, gridcolor="#2f3445"),
        #legend=dict(x=0, y=1.05, orientation="h"),
        #margin=dict(l=100, r=10, t=25, b=40),
        plot_bgcolor="#282b38",
        paper_bgcolor="#282b38",
        font={"color": "#a5b1cd"},
    )

    data = [trace0]
    figure = go.Figure(data=data, layout=layout)

    return figure

def isotone_chain(N, ZRange, model, axis_label, func):

    output = []

    layout = go.Layout(
        #title=f"ROC Curve (AUC = {auc_score:.3f})",
        title=f"Isotonic Chain",
        xaxis=dict(title="Protons", gridcolor="#2f3445",title_font_size=14),
        yaxis=dict(title=axis_label, gridcolor="#2f3445",title_font_size=14),
        #legend=dict(x=0, y=1.05, orientation="h"),
        #margin=dict(l=100, r=10, t=25, b=40),
        plot_bgcolor="#282b38",
        paper_bgcolor="#282b38",
        font={"color": "#a5b1cd", "size": 14},
    )

    Z = bmex.df[(bmex.df["N"]==N) & (bmex.df["Model"]==model) & (bmex.df["Z"] >= ZRange[0]) & (bmex.df["Z"] <= ZRange[1])]["Z"]
    for i in Z:
        output.append(func(N, i, model))

    trace0 = go.Scatter(
        x=Z, y=output, mode="lines+markers", name="Test Data", marker=\
            {
                "color": "#13c6e9",
                #"size": 20,
            }
    )

    #figure = px.line(x=Z, y=output, markers=True)
    data = [trace0]
    figure = go.Figure(data=data, layout=layout)
    figure.update_xaxes(title_font_size=20)
    figure.update_yaxes(title_font_size=20)
    figure.update_layout(title_font_size=24)
    return figure

def serve_pie_confusion_matrix(model, X_test, y_test, Z, threshold):
    # Compute threshold
    scaled_threshold = threshold * (Z.max() - Z.min()) + Z.min()
    y_pred_test = (model.decision_function(X_test) > scaled_threshold).astype(int)

    matrix = metrics.confusion_matrix(y_true=y_test, y_pred=y_pred_test)
    tn, fp, fn, tp = matrix.ravel()

    values = [tp, fn, fp, tn]
    label_text = ["True Positive", "False Negative", "False Positive", "True Negative"]
    labels = ["TP", "FN", "FP", "TN"]
    blue = cl.flipper()["seq"]["9"]["Blues"]
    red = cl.flipper()["seq"]["9"]["Reds"]
    colors = ["#13c6e9", blue[1], "#ff916d", "#ff744c"]

    trace0 = go.Pie(
        labels=label_text,
        values=values,
        hoverinfo="label+value+percent",
        textinfo="text+value",
        text=labels,
        sort=False,
        marker=dict(colors=colors),
        insidetextfont={"color": "white"},
        rotation=90,
    )

    layout = go.Layout(
        title="Confusion Matrix",
        margin=dict(l=50, r=50, t=100, b=10),
        legend=dict(bgcolor="#282b38", font={"color": "#a5b1cd"}, orientation="h"),
        plot_bgcolor="#282b38",
        paper_bgcolor="#282b38",
        font={"color": "#a5b1cd"},
    )

    data = [trace0]
    figure = go.Figure(data=data, layout=layout)

    return figure

max_stars_repo_name: kylegodbey/bmex-web

max_issues_repo_name: kylegodbey/bmex-web

max_forks_repo_name: kylegodbey/bmex-web


text:
import os
from random import random
import sys

from math import factorial as f
import numpy as np
import torch
import torchvision

import common.config as cfg

class AddGaussianNoise():
    def __init__(self, mean=0., std=1.):
        self.std = std
        self.mean = mean

    def __call__(self, tensor):
        return tensor + torch.randn(tensor.size()) * self.std + self.mean

    def __repr__(self):
        return self.__class__.__name__ + '(mean={0}, std={1})'.format(self.mean, self.std)

def add_gaussian_noise(args, tensor):
    noise = AddGaussianNoise(0.0, args.noise)
    return noise(tensor)

def booltype(arg):
    return bool(int(arg))

def decimal_format(num, places=4):
    return round(num, places)

def decor_print(msg):
    print('+'*80)
    print(msg)

def flip(p):
    return True if random() < p else False

def get_device(args):
    USE_CUDA = not args.no_cuda and torch.cuda.is_available()

    torch.manual_seed(args.seed)
    return torch.device("cuda:{}".format(
        args.device_id[0]) if USE_CUDA else "cpu")

def get_paths(args, distributed=False):
    ckpt_path = cfg.ckpt_dir
    tb_path = cfg.tb_dir
    if args.num_workers:
        folder = '{}_{}'.format(args.dataset, args.num_workers)
    else:
        folder = args.dataset
    path = {}

    model_name = 'clf_{}_optim_{}_uniform_{}_non_iid_{}' \
                 '_num_workers_{}_lr_{}_decay_{}_batch_{}'.format(
                     args.clf, args.optim, args.uniform_data,
                     args.non_iid if args.non_iid else 0,
                     args.num_workers if args.num_workers else 0,
                     args.lr, args.decay, args.batch_size
                 )
    pca_ = args.ncomponent if args.ncomponent else args.pca_var
    dga_ = 'dga_{}_bs_{}'.format(args.num_dga, args.dga_bs)
    path['pca_path'] = os.path.join(
        ckpt_path, folder, 'models',  model_name + '.pca_{}'.format(pca_))
    path['dga_path'] = os.path.join(
        ckpt_path, folder, 'models', 'clf_{}_{}_non_iid_{}.pkl'.format(
            args.clf, dga_, args.non_iid
        ))

    if args.dry_run:
        print(model_name)
        model_name = 'debug'
    if args.noise:
        model_name += '_noise_{}'.format(args.noise)
    if args.paradigm:
        if 'atomo' in args.paradigm:
            model_name += '_atomo_{}'.format(args.atomo_r)
        if 'signsgd' in args.paradigm:
            model_name += '_signsgd'
        if 'topk' in args.paradigm:
            model_name += '_topk_{}'.format(args.topk)
        if 'pca' in args.paradigm:
            model_name += '_pca_{}'.format(pca_)
        if 'lbgm' in args.paradigm:
            model_name += '_lbgm_{}'.format(args.error_tol)
        elif 'rp' in args.paradigm:
            model_name += '_rp_{}'.format(
                args.rp_eps if not args.ncomponent else args.ncomponent)
        elif 'kgrad' in args.paradigm:
            model_name += '_kgrad_{}_tol_{}'.format(
                args.kgrads, args.error_tol)
        if args.residual:
            model_name += '_residual'
        if args.sdir_full:
            model_name += '_full'
    if args.start_epoch != 1:
        model_name += '_start_from_{}'.format(args.start_epoch)
    if distributed:
        model_name += '_distributed'

    path['model_name'] = model_name
    path['log_file'] = '{}/{}/logs/{}.log'.format(
        ckpt_path, folder, model_name)
    path['tb_path'] = '{}/{}_{}'.format(
        tb_path, folder, model_name)
    path['init_path'] = '{}/{}/{}_{}.init'.format(ckpt_path, 'init',
                                                  args.dataset, args.clf)
    path['best_path'] = os.path.join(
        ckpt_path, folder, 'models',  model_name + '.best')
    path['stop_path'] = os.path.join(
        ckpt_path, folder, 'models',  model_name + '.stop')
    path['data_path'] = '../ckpts/{}_{}/data/n_classes_per_node_{}' \
                        '_stratify_{}_uniform_{}_repeat_{}.pkl'.format(
                            args.dataset, args.num_workers, args.non_iid,
                            args.stratify, args.uniform_data, args.repeat)
    path['hist_path'] = '{}/{}/history/{}.pkl'.format(
        ckpt_path, folder, model_name)
    path['plot_path'] = '{}/{}/plots/{}.jpg'.format(
        ckpt_path, folder, model_name)

    return Struct(**path)

def get_random_vector(mean, std, size):
    return torch.normal(mean, std, size)

def in_range(elem, upper, lower):
    return (elem >= lower) and (elem <= upper)

def init_logger(log_file, dry_run, append):
    print("Logging: ", log_file)
    std_out = sys.stdout
    if not dry_run:
        log_file = open(
            log_file,
            'a' if os.path.exists(log_file) and append else 'w'
        )
        sys.stdout = log_file

    return log_file, std_out

def is_approx(args):
    return args.paradigm and (
        'kgrad' in args.paradigm or
        'rp' in args.paradigm or
        'pca' in args.paradigm or
        'dga' in args.paradigm
    )

def nCr(n, r):
    return f(n)//f(r)//f(n-r)

class Struct:
    def __init__(self, **entries):
        self.__dict__.update(entries)

def tb_model_summary(model, test_loader, tb, device):
    for data, _ in test_loader:
        grid = torchvision.utils.make_grid(data[:64])
        tb.add_image('test images', grid)
        tb.add_graph(model.module, data.to(device))
        break

def vec_angle(firsts, currents):
    angles = []
    for a, b in zip(firsts, currents):
        dot = torch.dot(a.flatten(), b.flatten())
        a_norm = torch.norm(a)
        b_norm = torch.norm(b)
        cos = dot/(2*a_norm*b_norm)
        angles.append(torch.acos(cos).item() * 180/3.14)
    angles = np.array(angles)

    return angles.mean(), angles.std()

def vec_unit_dot(a, b):
    dot = torch.dot(a.flatten(), b.flatten())
    a_norm = torch.norm(a)
    b_norm = torch.norm(b)
    cos = dot/(a_norm*b_norm)

    return cos.item()


text:
"""
>>> error()
Traceback (most recent call last):
    ...
NumbaError: 17:0: object of type int cannot be sliced
"""

import sys
import ctypes

from numba import *
import numpy as np

@autojit(backend='ast')
def error():
    i = 10
    i[:]

if __name__ == "__main__":
    import doctest
    doctest.testmod()

max_stars_repo_name: aseyboldt/numba

max_issues_repo_name: aseyboldt/numba

max_forks_repo_name: aseyboldt/numba

reason: import numpy,from numba


text:
from .. cupy_utils import to_numpy, trapz, xp
from ..utils import powerlaw

import numpy as np

from astropy.cosmology import Planck15

class PowerLawRedshift(object):
    """
    Redshift model from Fishbach+ https://arxiv.org/abs/1805.10270
    Note that this is deliberately off by a factor of dVc/dz
    """

    def __init__(self):
        self.zs_ = np.linspace(1e-3, 1, 1000)
        self.zs = xp.asarray(self.zs_)
        self.dvc_dz_ = (
            Planck15.differential_comoving_volume(self.zs_).value * 4 * np.pi)
        self.dvc_dz = xp.asarray(self.dvc_dz_)
        self.cached_dvc_dz = None

    def __call__(self, dataset, lamb):
        p_z = powerlaw(1 + dataset['redshift'], alpha=(lamb - 1),
                       high=(1 + self.zs_[-1]), low=1)
        try:
            p_z *= self.cached_dvc_dz
        except (TypeError, ValueError):
            self._cache_dvc_dz(dataset['redshift'])
            p_z *= self.cached_dvc_dz
        p_z /= self.normalisation(lamb)
        return p_z

    def normalisation(self, lamb):
        p_z_ = powerlaw(1 + self.zs, alpha=(lamb - 1),
                        high=(1 + self.zs_[-1]), low=1)
        norm = trapz(p_z_ * self.dvc_dz, self.zs)
        return norm

    def _cache_dvc_dz(self, redshifts):
        self.cached_dvc_dz = xp.asarray(np.interp(
            to_numpy(redshifts), self.zs_, self.dvc_dz_))

power_law_redshift = PowerLawRedshift()

max_stars_repo_name: chase-kimball/gwpopulation

max_issues_repo_name: chase-kimball/gwpopulation

max_forks_repo_name: chase-kimball/gwpopulation


text:
import numpy as np

class Cox(object):
    def __init__(self):
        pass

    def fit(self):
        pass

    def predict(self):
        pass

    def plot(self):
        pass


text:
import datetime as dt
import numpy as np
import pandas as pd

import sqlalchemy
from sqlalchemy.ext.automap import automap_base
from sqlalchemy.orm import Session
from sqlalchemy import create_engine, func

from flask import Flask, jsonify

# Database Setup
engine = create_engine("sqlite:///hawaii.sqlite")
# reflect an existing database into a new model
Base = automap_base()
# reflect the tables
Base.prepare(engine, reflect=True)

Measurement = Base.classes.measurement
Station = Base.classes.station
session = Session(engine)
#################################################
# Flask Setup
#################################################
app = Flask(__name__)
@app.route("/")
def welcome():
 return(
        f"Available Routes:
"
        f"
"
        f"/api/v1.0/precipitation"
         f"
"
        f"/api/v1.0/stations"
         f"
"
        f"/api/v1.0/tobs"
         f"
"
        f"/api/v1.0/<start>"
         f"
"
        f"/api/v1.0/<start>/<end>")

@app.route("/api/v1.0/precipitation")
def precipitation():
           # Create our session (link) from Python to the DB
    session = Session(engine)

    precipitation_score = session.query(Measurement.date,Measurement.prcp).filter(Measurement.date >='2016-08-23').filter(Measurement.date <='2017-08-24').order_by(Measurement.date).all()
    session.close()
    percp_total = []

    for date, prcp in precipitation_score:
        perc_date = {}
        perc_date["date"] = date
        perc_date["prcp"] = prcp
        percp_total.append(perc_date)
    return jsonify(percp_total)

@app.route("/api/v1.0/stations")
def stations ():
    session = Session(engine)

    station_count=session.query(Station.name,Station.station)
    session.close()
    stations = []
    for name,station in station_count:
        station_total={}
        station_total["name"] = name
        station_total["station"] = station
        stations.append(station_total)
    return jsonify(stations)

@app.route("/api/v1.0/tobs")
def tobs():
    session = Session(engine)
    last_date_entered=session.query(func.max(func.strftime("%Y-%m-%d", Measurement.date))).all()
    query_date=dt.date(2017,8,23) - dt.timedelta(days=365)
    precipitation_score = session.query(Measurement.date,Measurement.prcp).filter(Measurement.date >='2016-08-23').filter(Measurement.date <='2017-08-24').order_by(Measurement.date).all()
    session.close()
    tobs = []
    for date, prcp in precipitation_score:
        perc_total = {}
        perc_total["date"] = date
        perc_total["prcp"] = prcp
        tobs.append(perc_total)
    return jsonify(tobs)

@app.route ("/api/v1.0/<start>")
def start1(start):
    start = dt.datetime.strptime(start, '%Y-%m-%d')
    one_year = dt.timedelta(days=365)
    input_date = start - one_year
    end_date = dt.date(2017,8,23)
    trip_data=session.query(func.min(Measurement.tobs), func.avg(Measurement.tobs), func.max(Measurement.tobs)).\
        filter(Measurement.date >= input_date).filter(Measurement.date <= end_date).all()
    trip = list(trip_data)
    return jsonify(trip)

@app.route ("/api/v1.0/<start>/<end>")
def startend(start,end):
    start = dt.datetime.strptime(start, '%Y-%m-%d')
    end=dt.datetime.strptime(end, '%Y-%m-%d')
    one_year = dt.timedelta(days=365)
    input_start_date = start - one_year
    input_end_date = end - one_year
    trip_data=session.query(func.min(Measurement.tobs), func.avg(Measurement.tobs), func.max(Measurement.tobs)).\
        filter(Measurement.date >= input_start_date).filter(Measurement.date <= input_end_date).all()
    trip = list(trip_data)
    return jsonify(trip)
if __name__ == "__main__":
    app.run(debug=True)

max_stars_repo_name: stewart-lacy/sqlalchemy-challenge

max_issues_repo_name: stewart-lacy/sqlalchemy-challenge

max_forks_repo_name: stewart-lacy/sqlalchemy-challenge


text:
# -*- coding: utf-8 -*-
"""
Created on Fri May 21 15:12:22 2021

@author: ryans
"""
# Initialize GPU
import pyclesperanto_prototype as cle
import beetlesafari as bs
import numpy as np
from skimage.io import imread
import tribolium_clustering as tc
from skimage import io

'''
============================ Wolgast ========================================
'''

import pandas as pd
cle.select_device('GTX')

cc_dataset = bs.ClearControlDataset(tc.get_folder_path('Choose Wolgast Dataset Directory'))
savedata_location = tc.get_folder_path('Save Location')

print("Duration in seconds:", cc_dataset.get_duration_in_seconds())
print("Duration in hours:", bs.seconds_to_hours(cc_dataset.get_duration_in_seconds()))

timeslice_duration = bs.hours_to_seconds(3)/32
start_time = bs.hours_to_seconds(11)
gastrulation_times = [start_time + (i+1) * timeslice_duration for i in range(32)]
print('beginning: {} hours'.format(bs.seconds_to_hours(gastrulation_times[0])))
print('ending   : {} hours'.format(bs.seconds_to_hours(gastrulation_times[-1])))
voxel_size = cc_dataset.get_voxel_size_zyx()

test_time = [gastrulation_times[0]]
for i in gastrulation_times:
    rescaled = tc.resample_isotropic_nodownsample(cc_dataset.get_image(time_in_seconds=i)
                                                  , voxel_size)
    io.imsave(savedata_location + 'rescaled t = {}s.tif'.format(i), rescaled)
    print('dopamine')

rescaled = None
gpu_output = None

np.save(savedata_location + 'timepoints.npy', gastrulation_times)
#savedata_location = tc.get_folder_path('Save Location')
#np.save(savedata_location + 'rotation.npy', (-30,180,85))

'''
=============================== Finsterwalde =================================
'''
'''
import pandas as pd
cle.select_device('GTX')

cc_dataset = bs.ClearControlDataset('D://2019-10-28-17-22-59-23-Finsterwalde_Tribolium_nGFP')
savedata_location = 'C://Users//ryans//OneDrive//Documents//Master Thesis//output data//'

print("Duration in seconds:", cc_dataset.get_duration_in_seconds())
print("Duration in hours:", bs.seconds_to_hours(cc_dataset.get_duration_in_seconds()))

timeslice_duration = bs.hours_to_seconds(3)/32
start_time = bs.hours_to_seconds(13)
gastrulation_times = [start_time + (i+1) * timeslice_duration for i in range(32)]
print('beginning: {} hours'.format(bs.seconds_to_hours(gastrulation_times[0])))
print('ending   : {} hours'.format(bs.seconds_to_hours(gastrulation_times[-1])))
voxel_size = cc_dataset.get_voxel_size_zyx()

test_time = [gastrulation_times[0]]
for i in gastrulation_times:
    rescaled = cle.pull(cc_dataset.get_resampled_image(time_in_seconds=i))
    io.imsave(savedata_location + 'Finsterwalde Gastrulation Labels//rescaled t = {}s.tif'.format(i), rescaled)
    print('satisfaction')
'''

max_stars_repo_path: Code/Image Analysis and Feature Extraction/saving resclaed images to avoid later processing for montefano and finsterwalde.py

max_stars_repo_name: Cryaaa/Master-Thesis-Repository

max_issues_repo_name: Cryaaa/Master-Thesis-Repository

max_forks_repo_name: Cryaaa/Master-Thesis-Repository


text:
# based on https://github.com/tensorflow/models/tree/master/resnet
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np
import tensorflow as tf

class Model(object):
    """ResNet model."""

    def __init__(self, mode):
        """ResNet constructor.

        Args:
          mode: One of 'train' and 'eval'.
        """
        self.mode = mode
        self._build_model()

    def add_internal_summaries(self):
        pass

    def _stride_arr(self, stride):
        """Map a stride scalar to the stride array for tf.nn.conv2d."""
        return [1, stride, stride, 1]

    def _build_model(self):
        assert self.mode == 'train' or self.mode == 'eval'
        """Build the core model within the graph."""
        with tf.variable_scope('input'):

            self.x_input = tf.placeholder(
                tf.float32,
                shape=[None, 32, 32, 3])

            self.y_input = tf.placeholder(tf.int64, shape=None)

            input_standardized = tf.map_fn(lambda img: tf.image.per_image_standardization(img),
                                           self.x_input)
            x = self._conv('init_conv', input_standardized, 3, 3, 16, self._stride_arr(1))

        strides = [1, 2, 2]
        activate_before_residual = [True, False, False]
        res_func = self._residual

        # Uncomment the following codes to use w28-10 wide residual network.
        # It is more memory efficient than very deep residual network and has
        # comparably good performance.
        # https://arxiv.org/pdf/1605.07146v1.pdf
        filters = [16, 160, 320, 640]

        # Update hps.num_residual_units to 9

        with tf.variable_scope('unit_1_0'):
            x = res_func(x, filters[0], filters[1], self._stride_arr(strides[0]), activate_before_residual[0])
        for i in range(1, 5):
            with tf.variable_scope('unit_1_%d' % i):
                x = res_func(x, filters[1], filters[1], self._stride_arr(1), False)

        with tf.variable_scope('unit_2_0'):
            x = res_func(x, filters[1], filters[2], self._stride_arr(strides[1]), activate_before_residual[1])
        for i in range(1, 5):
            with tf.variable_scope('unit_2_%d' % i):
                x = res_func(x, filters[2], filters[2], self._stride_arr(1), False)

        with tf.variable_scope('unit_3_0'):
            x = res_func(x, filters[2], filters[3], self._stride_arr(strides[2]), activate_before_residual[2])
        for i in range(1, 5):
            with tf.variable_scope('unit_3_%d' % i):
                x = res_func(x, filters[3], filters[3], self._stride_arr(1), False)

        with tf.variable_scope('unit_last'):
            x = self._batch_norm('final_bn', x)
            x = self._relu(x, 0.1)
            x = self._global_avg_pool(x)

        with tf.variable_scope('logit'):
            self.pre_softmax = self._fully_connected(x, 10)

        self.predictions = tf.argmax(self.pre_softmax, 1)
        self.correct_prediction = tf.equal(self.predictions, self.y_input)
        self.num_correct = tf.reduce_sum(
            tf.cast(self.correct_prediction, tf.int64))
        self.accuracy = tf.reduce_mean(
            tf.cast(self.correct_prediction, tf.float32))

        with tf.variable_scope('costs'):
            self.y_xent = tf.nn.sparse_softmax_cross_entropy_with_logits(logits=self.pre_softmax, labels=self.y_input)
            self.xent = tf.reduce_sum(self.y_xent, name='y_xent')
            self.mean_xent = tf.reduce_mean(self.y_xent)
            self.weight_decay_loss = self._decay()

    def _batch_norm(self, name, x):
        """Batch normalization."""
        with tf.name_scope(name):
            return tf.contrib.layers.batch_norm(
                inputs=x,
                decay=.9,
                center=True,
                scale=True,
                activation_fn=None,
                updates_collections=None,
                is_training=(self.mode == 'train'))

    def _residual(self, x, in_filter, out_filter, stride,
                  activate_before_residual=False):
        """Residual unit with 2 sub layers."""
        if activate_before_residual:
            with tf.variable_scope('shared_activation'):
                x = self._batch_norm('init_bn', x)
                x = self._relu(x, 0.1)
                orig_x = x
        else:
            with tf.variable_scope('residual_only_activation'):
                orig_x = x
                x = self._batch_norm('init_bn', x)
                x = self._relu(x, 0.1)

        with tf.variable_scope('sub1'):
            x = self._conv('conv1', x, 3, in_filter, out_filter, stride)

        with tf.variable_scope('sub2'):
            x = self._batch_norm('bn2', x)
            x = self._relu(x, 0.1)
            x = self._conv('conv2', x, 3, out_filter, out_filter, [1, 1, 1, 1])

        with tf.variable_scope('sub_add'):
            if in_filter != out_filter:
                orig_x = tf.nn.avg_pool(orig_x, stride, stride, 'VALID')
                orig_x = tf.pad(
                    orig_x, [[0, 0], [0, 0], [0, 0],
                             [(out_filter - in_filter) // 2, (out_filter - in_filter) // 2]])
            x += orig_x

        tf.logging.debug('image after unit %s', x.get_shape())
        return x

    def _decay(self):
        """L2 weight decay loss."""
        costs = []
        for var in tf.trainable_variables():
            if var.op.name.find('DW') > 0:
                costs.append(tf.nn.l2_loss(var))
        return tf.add_n(costs)

    def _conv(self, name, x, filter_size, in_filters, out_filters, strides):
        """Convolution."""
        with tf.variable_scope(name):
            n = filter_size * filter_size * out_filters
            kernel = tf.get_variable(
                'DW', [filter_size, filter_size, in_filters, out_filters],
                tf.float32, initializer=tf.random_normal_initializer(
                    stddev=np.sqrt(2.0 / n)))
            return tf.nn.conv2d(x, kernel, strides, padding='SAME')

    def _relu(self, x, leakiness=0.0):
        """Relu, with optional leaky support."""
        return tf.where(tf.less(x, 0.0), leakiness * x, x, name='leaky_relu')

    def _fully_connected(self, x, out_dim):
        """FullyConnected layer for final output."""
        num_non_batch_dimensions = len(x.shape)
        prod_non_batch_dimensions = 1
        for ii in range(num_non_batch_dimensions - 1):
            prod_non_batch_dimensions *= int(x.shape[ii + 1])
        x = tf.reshape(x, [tf.shape(x)[0], -1])
        w = tf.get_variable(
            'DW', [prod_non_batch_dimensions, out_dim],
            initializer=tf.uniform_unit_scaling_initializer(factor=1.0))
        b = tf.get_variable('biases', [out_dim],
                            initializer=tf.constant_initializer())
        return tf.nn.xw_plus_b(x, w, b)

    def _global_avg_pool(self, x):
        assert x.get_shape().ndims == 4
        return tf.reduce_mean(x, [1, 2])


text:
from os import PathLike

import click
import numpy as np

from .data_munging import read_tlt_params, angle2rotm

def emanjson2imodxf(input_json_file: PathLike):
    # Get useful info from EMAN json file
    # shifts in px, theta in radians, y tilt in degrees
    dxy_eman, theta, tilt_angles, tilt_azimuth = read_tlt_params(
        json_file=input_json_file
    )
    rotation_matrices = angle2rotm(theta + tilt_azimuth)

    # EMAN translates then rotates, IMOD rotates then translates
    # R(x+dx, y+dy) = R(x, y) + R(dx, dy)
    # IMOD (dx, dy) is the rotated EMAN (dx, dy)
    dxy_imod = rotation_matrices @ dxy_eman

    # get lines for IMOD xf file
    A11 = rotation_matrices[:, 0, 0]
    A12 = rotation_matrices[:, 0, 1]
    A21 = rotation_matrices[:, 1, 0]
    A22 = rotation_matrices[:, 1, 1]
    dx_imod = dxy_imod[:, 0].squeeze()
    dy_imod = dxy_imod[:, 1].squeeze()

    xf = np.stack((A11, A12, A21, A22, dx_imod, dy_imod), axis=-1)
    return xf, tilt_angles

@click.command()
@click.option(
    '-i', '--input', 'input_json_file',
    help='EMAN2 json file',
    type=click.Path(exists=True),
    required=True,
)
@click.option(
    '-o', '--output-basename', 'output_basename',
    help='basename for IMOD .xf and .tlt files',
    type=click.Path(),
    required=True,
)
def cli(input_json_file, output_basename):
    """Convert an EMAN2 json file into an IMOD xf file.
    """
    output_xf_file = f'{output_basename}.xf'
    output_tlt_file = f'{output_basename}.tlt'

    click.echo(f'Input file: {input_json_file}')
    click.echo(f'Output files: {output_xf_file}, {output_tlt_file}')

    xf, tilt_angles = emanjson2imodxf(input_json_file)
    np.savetxt(output_xf_file, xf, fmt='%.6f', delimiter=' ')
    np.savetxt(output_tlt_file, tilt_angles, fmt='%.2f')
    click.echo(f'Success!')