text:
'''
fetch_pdbs.py
Updated: 1/19/17

Script fetches pdbs from Protein Data Bank as defined in class.csv files in
data_folder path. class.csv files includes pdb identifier and chain identifier
pairs. PDBs are save into folders corresponding to class and PDB identifier.

Script parallelizes fetches over multiple cores using MPI.

To run: $python3 fetch_pdbs.py or $mpirun -n N python3 fetch_pdbs.py

'''
import os, wget
import numpy as np
from mpi4py import MPI

# Data folder path
data_folder = '../../data/PsiBlast/'

###############################################################################

seed = 1234

def fetch_PDB(path, pdb_id):
    '''
    Method fetches pdb file from Protein Data Bank repo and stores file in
    designated path.

    '''
    # Download PDB file from database.
    url = 'https://files.rcsb.org/download/' # URL used to fetch PDB files
    file_path = path + pdb_id.lower() + '.pdb'
    if not os.path.exists(file_path):
        file_path = wget.download(url + pdb_id + '.pdb', out=path, bar=None)

if __name__ == '__main__':

    # Set paths relative to this file
    os.chdir(os.path.dirname(os.path.realpath(__file__)))

    # MPI init
    comm = MPI.COMM_WORLD
    rank = comm.Get_rank()
    cores = comm.Get_size()

    # MPI task distribution
    if rank == 0:
        tasks = []

        # Iterate over PDB class .csv files
        for class_fn in sorted(os.listdir(data_folder)):
            if not class_fn.endswith('.csv'): continue

            # Make folder for PDB class
            class_ = class_fn.split('.')[0]
            if not os.path.exists(data_folder+class_): os.mkdir(data_folder+class_)

            # Iterate over PDB id and chain id pairs
            with open(data_folder+class_fn, 'r')as f:
                lines = f.readlines()
                for l in lines:

                    # Parse PDB IDs and chain IDs
                    l = l[:-1].split(',')
                    pdb_id = l[0].lower()
                    chain_id = l[1]
                    tasks.append([data_folder+class_, pdb_id, chain_id])

        # Shuffle for Random Distribution
        np.random.seed(seed)
        np.random.shuffle(tasks)

    else: tasks = None

    # Broadcast tasks to all nodes and select tasks according to rank
    tasks = comm.bcast(tasks, root=0)
    tasks = np.array_split(tasks, cores)[rank]

    # Fetch PDBs
    for t in tasks:

        # Task IDs
        folder_ = t[0]
        pdb_id = t[1]
        chain_id = t[2]

        # Fetch PDB file and rename with task IDs
        fetch_PDB(folder_, pdb_id)
        if not os.path.exists(folder_+'/'+ pdb_id +'_'+chain_id): os.mkdir(folder_+'/'+ pdb_id +'_'+chain_id)
        os.rename(folder_+'/'+pdb_id+'.pdb', folder_+'/'+pdb_id+'_'+chain_id+'/'+pdb_id+'_'+chain_id+'.pdb')

        print('Fetched: ' + pdb_id)

max_stars_repo_name: tjosc/Protein-Structure-Exploration

max_issues_repo_name: tjosc/Protein-Structure-Exploration

max_forks_repo_name: tjosc/Protein-Structure-Exploration


text:
#!/usr/bin/python3

import numpy as np
import matplotlib.pyplot as plt

anchors = np.array([
    [0, 0, 0],
    [1, 0, 0],
    [0, 1, 0],
    [0, 0, 1],
    [1, 1, 1]
])

N = 5
x = np.linspace(-5, 5, N)
y = np.linspace(-5, 5, N)
z = np.linspace(-5, 5, N)

X = np.zeros((3,np.power(N, 3)))

index = 0
while index < X.shape[1]:
    x_index = index % N
    y_index = int(index / N) % N
    z_index = int(index / (N*N))
    X[:, index] = [x[x_index], y[y_index], z[z_index]]
    index += 1

noise = np.random.normal(size=X.shape)
x_noisy = X + 0.1*noise

A = 2*np.array([
    [anchors[0, 0] - anchors[1, 0], anchors[0, 1] -
        anchors[1, 1], anchors[0, 2] - anchors[1, 2]],
    [anchors[0, 0] - anchors[2, 0], anchors[0, 1] -
        anchors[2, 1], anchors[0, 2] - anchors[2, 2]],
    [anchors[0, 0] - anchors[3, 0], anchors[0, 1] -
        anchors[3, 1], anchors[0, 2] - anchors[3, 2]],
    [anchors[0, 0] - anchors[4, 0], anchors[0, 1] -
        anchors[4, 1], anchors[0, 2] - anchors[4, 2]]
])

err = np.array([])
print(x_noisy.shape)

index = 0
while index < x_noisy.shape[1]:

    b = np.array([
    [
        np.square(np.linalg.norm(anchors[1] - x_noisy[:,index])) - np.square(np.linalg.norm(anchors[0] - x_noisy[:,index])) -
        np.square(anchors[1, 0]) + np.square(anchors[0, 0]) -
        np.square(anchors[1, 1]) + np.square(anchors[0, 1]) -
        np.square(anchors[1, 2]) + np.square(anchors[0, 2])
    ],
    [
        np.square(np.linalg.norm(anchors[2] - x_noisy[:,index])) - np.square(np.linalg.norm(anchors[0] - x_noisy[:,index])) -
        np.square(anchors[2, 0]) + np.square(anchors[0, 0]) -
        np.square(anchors[2, 1]) + np.square(anchors[0, 1]) -
        np.square(anchors[2, 2]) + np.square(anchors[0, 2])
    ],
    [
        np.square(np.linalg.norm(anchors[3] - x_noisy[:,index])) - np.square(np.linalg.norm(anchors[0] - x_noisy[:,index])) -
        np.square(anchors[3, 0]) + np.square(anchors[0, 0]) -
        np.square(anchors[3, 1]) + np.square(anchors[0, 1]) -
        np.square(anchors[3, 2]) + np.square(anchors[0, 2])
    ],
    [
    np.square(np.linalg.norm(anchors[4] - x_noisy[:,index])) - np.square(np.linalg.norm(anchors[0] - x_noisy[:,index])) -
        np.square(anchors[4, 0]) + np.square(anchors[0, 0]) -
        np.square(anchors[4, 1]) + np.square(anchors[0, 1]) -
        np.square(anchors[4, 2]) + np.square(anchors[0, 2])
    ]
    ])

    print(x_noisy[:,index])
    x_estimated = np.linalg.lstsq(A, b, rcond=None)[0]
    print(x_estimated.T)
    err = np.append(err, np.linalg.norm(X[:,index]-x_estimated.T))
    print(err[index])
    index += 1

print("avg: ", np.average(err))
print("min: ", np.min(err))
print("max: ", np.max(err))

ext: py

max_stars_repo_name: mguc/arts

max_issues_repo_name: mguc/arts

max_forks_repo_name: mguc/arts

reason: import numpy


text:
# Rainbow RL Algorithm class
# Based on: https://github.com/higgsfield/RL-Adventure

import numpy as np

import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F

from training.training import Trainer
from common.replay_buffer import ReplayBuffer
from models.models import RainbowDQN

USE_CUDA = torch.cuda.is_available()
def Variable(x): return x.cuda() if USE_CUDA else x

num_atoms = 51
Vmin = -10
Vmax = 10

class Rainbow(Trainer):

    def __init__(self, parameters):
        super(Rainbow, self).__init__(parameters)
        self.replay_buffer = ReplayBuffer(self.buffersize)

    def push_to_buffer(self, state, action, reward, next_state, done):
        self.replay_buffer.push(state, action, reward, next_state, done)

    def load_model(self):
        self.current_model = RainbowDQN(self.env.observation_space.shape[0],
                                   self.env.action_space.n, num_atoms, Vmin, Vmax)  # input:(1,84,84), output:6
        self.target_model = RainbowDQN(
            self.env.observation_space.shape[0], self.env.action_space.n, num_atoms, Vmin, Vmax)

        if USE_CUDA:
            self.current_model = self.current_model.cuda()
            self.target_model = self.target_model.cuda()

        self.update_target(self.current_model, self.target_model)  # sync nets

    def projection_distribution(self, next_state, rewards, dones):
        batch_size = next_state.size(0)

        delta_z = float(Vmax - Vmin) / (num_atoms - 1)
        support = torch.linspace(Vmin, Vmax, num_atoms)

        next_dist = self.target_model(next_state).data.cpu() * support
        next_action = next_dist.sum(2).max(1)[1]
        next_action = next_action.unsqueeze(1).unsqueeze(
            1).expand(next_dist.size(0), 1, next_dist.size(2))
        next_dist = next_dist.gather(1, next_action).squeeze(1)

        rewards = rewards.unsqueeze(1).expand_as(next_dist)
        dones = dones.unsqueeze(1).expand_as(next_dist)
        support = support.unsqueeze(0).expand_as(next_dist)

        Tz = rewards + (1 - dones) * 0.99 * support
        Tz = Tz.clamp(min=Vmin, max=Vmax)
        b = (Tz - Vmin) / delta_z
        l = b.floor().long()
        u = b.ceil().long()

        offset = torch.linspace(0, (batch_size - 1) * num_atoms, batch_size).long()\
            .unsqueeze(1).expand(batch_size, num_atoms)

        proj_dist = torch.zeros(next_dist.size())
        proj_dist.view(-1).index_add_(0, (l + offset).view(-1),
                                    (next_dist * (u.float() - b)).view(-1))
        proj_dist.view(-1).index_add_(0, (u + offset).view(-1),
                                    (next_dist * (b - l.float())).view(-1))

        return proj_dist

    def compute_td_loss(self, batch_size, *args):
        state, action, reward, next_state, done = self.replay_buffer.sample(batch_size)

        state = Variable(torch.FloatTensor(np.float32(state)))
        next_state = Variable(torch.FloatTensor(np.float32(next_state)))
        action = Variable(torch.LongTensor(action))
        reward = torch.FloatTensor(reward)
        done = torch.FloatTensor(np.float32(done))

        proj_dist = self.projection_distribution(next_state, reward, done)

        dist = self.current_model(state)
        action = action.unsqueeze(1).unsqueeze(1).expand(batch_size, 1, num_atoms)
        dist = dist.gather(1, action).squeeze(1)
        dist.data.clamp_(0.01, 0.99)
        loss = -(Variable(proj_dist) * dist.log()).sum(1)
        loss = loss.mean()

        self.optimizer.zero_grad()
        loss.backward()
        self.optimizer.step()

        self.current_model.reset_noise()
        self.target_model.reset_noise()

        return loss

max_stars_repo_name: huestisk/chess

max_issues_repo_name: huestisk/chess

max_forks_repo_name: huestisk/chess


text:
"""Use a generator with numpy loadtxt"""

import numpy as np

def file_stream(file_name):
    for line in open(file_name):
        yield line

def demo_use():
    file_name = 'census_abs2011_summary.csv'
    # Skip the header, return 3 columns: the age, weekly rental and total family income
    data = np.loadtxt(file_stream(file_name), delimiter=',', usecols=(0, 2, 3), skiprows=1)
    print('the loaded data')
    print(data.shape)
    print(data)

if __name__ == '__main__':
    demo_use()

max_stars_repo_name: lydiaknuefing/didbits

max_issues_repo_name: lydiaknuefing/didbits

max_forks_repo_name: lydiaknuefing/didbits


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

"""
Handy MCMC scripts.

Test for the different fit method (mcmc, ptmcmc, minimizer).

Author:
    Dimitri Misiak (misiak@ipnl.in2p3.fr)
"""

import numpy as np
import matplotlib.pyplot as plt
import sys
import scipy.signal as sgl
from os import path
import scipy.optimize as op

import mcmc_red as mcr

# close all plots
plt.close('all')

nsample = 100
#mu1, sigma1 = 1., 1.
mu1 = np.random.uniform(-10, 10)
sigma1 = np.random.uniform(0, 10)
#mu2, sigma2 = 2., 5.
mu2 = np.random.uniform(-10, 10)
sigma2 = np.random.uniform(0, 10)

print ("Generating blob at (mu1, mu2)=({0:.2f}, {1:.2f})"
       " and (sigma1, sigma2)=({2:.2f}, {3:.2f})".format(mu1, mu2, sigma1, sigma2))

blob = np.random.normal((mu1, mu2), (sigma1, sigma2), (nsample,2))

print("Checking")
print("mean =", np.mean(blob, axis=0))
print("std =", np.std(blob, axis=0))

def chi2(param):
    return mcr.chi2_simple(blob, param, (sigma1, sigma2))

# XXX MCMC
# save directory
sampler_path = 'mcmc_sampler/autosave'

# running the mcmc analysis
bounds = ((-20, 20),(-20, 20))
sampler = mcr.mcmc_sampler(chi2, bounds, nsteps=1000, path=sampler_path)

#    # loading the mcmc results
logd, chain, lnprob, acc = mcr.get_mcmc_sampler(sampler_path)

#    LAB = ('$log\ a$', '$log\ t$', '$log\ s$')
#LAB = ('$log\ a1$', '$log\ a2$', '$log\ t1$', '$log\ t2$', '$log\ s$')
lab = ('$\mu1$','$\mu2$')
dim = int(logd['dim'])
xopt, inf, sup = mcr.mcmc_results(dim, chain, lnprob, acc, lab)

print(xopt, inf, sup)

max_stars_repo_name: DimitriMisiak/mcmc-red

max_issues_repo_name: DimitriMisiak/mcmc-red

max_forks_repo_name: DimitriMisiak/mcmc-red

reason: import numpy,import scipy


text:
# -*- coding: utf-8 -*-
r"""
Output functions

These are the output functions for latexing and ascii/unicode art versions of
partitions and tableaux.

AUTHORS:

- Mike Hansen (?): initial version
- Andrew Mathas (2013-02-14): Added support for displaying conventions and
  lines, and tableaux of skew partition, composition, and
  skew/composition/partition/tableaux tuple shape.
- Travis Scrimshaw (2020-08): Added support for ascii/unicode art
"""

from __future__ import absolute_import, print_function

from string import Template
from sage.combinat.tableau import Tableaux

# The tex macro used to latex individual cells in an array (as a template).
# When using bar should be replaced by '|' or ''.
lr_macro = Template(r'\def\lr#1{\multicolumn{1}{$bar@{\hspace{.6ex}}c@{\hspace{.6ex}}$bar}{\raisebox{-.3ex}{$$#1$$}}}')

def tex_from_array(array, with_lines=True):
    r"""
    Return a latex string for a two dimensional array of partition, composition or skew composition shape

    INPUT:

    - ``array`` -- a list of list
    - ``with_lines`` -- a boolean (default: ``True``)
       Whether to draw a line to separate the entries in the array.

    Empty rows are allowed; however, such rows should be given as
    ``[None]`` rather than ``[]``.

    The array is drawn using either the English or French convention
    following :meth:`Tableaux.options`.

    .. SEEALSO:: :meth:`tex_from_array_tuple`

    EXAMPLES::

        sage: from sage.combinat.output import tex_from_array
        sage: print(tex_from_array([[1,2,3],[4,5]]))
        {\def\lr#1{\multicolumn{1}{|@{\hspace{.6ex}}c@{\hspace{.6ex}}|}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[b]{*{3}c}\cline{1-3}
        \lr{1}&\lr{2}&\lr{3}\\\cline{1-3}
        \lr{4}&\lr{5}\\\cline{1-2}
        \end{array}$}
        }
        sage: print(tex_from_array([[1,2,3],[4,5]], with_lines=False))
        {\def\lr#1{\multicolumn{1}{@{\hspace{.6ex}}c@{\hspace{.6ex}}}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[b]{*{3}c}\\
        \lr{1}&\lr{2}&\lr{3}\\
        \lr{4}&\lr{5}\\
        \end{array}$}
        }
        sage: print(tex_from_array([[1,2,3],[4,5,6,7],[8]]))
        {\def\lr#1{\multicolumn{1}{|@{\hspace{.6ex}}c@{\hspace{.6ex}}|}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[b]{*{4}c}\cline{1-3}
        \lr{1}&\lr{2}&\lr{3}\\\cline{1-4}
        \lr{4}&\lr{5}&\lr{6}&\lr{7}\\\cline{1-4}
        \lr{8}\\\cline{1-1}
        \end{array}$}
        }
        sage: print(tex_from_array([[1,2,3],[4,5,6,7],[8]], with_lines=False))
        {\def\lr#1{\multicolumn{1}{@{\hspace{.6ex}}c@{\hspace{.6ex}}}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[b]{*{4}c}\\
        \lr{1}&\lr{2}&\lr{3}\\
        \lr{4}&\lr{5}&\lr{6}&\lr{7}\\
        \lr{8}\\
        \end{array}$}
        }
        sage: print(tex_from_array([[None,None,3],[None,5,6,7],[8]]))
        {\def\lr#1{\multicolumn{1}{|@{\hspace{.6ex}}c@{\hspace{.6ex}}|}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[b]{*{4}c}\cline{3-3}
        &&\lr{3}\\\cline{2-4}
        &\lr{5}&\lr{6}&\lr{7}\\\cline{1-4}
        \lr{8}\\\cline{1-1}
        \end{array}$}
        }
        sage: print(tex_from_array([[None,None,3],[None,5,6,7],[None,8]]))
        {\def\lr#1{\multicolumn{1}{|@{\hspace{.6ex}}c@{\hspace{.6ex}}|}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[b]{*{4}c}\cline{3-3}
        &&\lr{3}\\\cline{2-4}
        &\lr{5}&\lr{6}&\lr{7}\\\cline{2-4}
        &\lr{8}\\\cline{2-2}
        \end{array}$}
        }
        sage: print(tex_from_array([[None,None,3],[None,5,6,7],[8]], with_lines=False))
        {\def\lr#1{\multicolumn{1}{@{\hspace{.6ex}}c@{\hspace{.6ex}}}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[b]{*{4}c}\\
        &&\lr{3}\\
        &\lr{5}&\lr{6}&\lr{7}\\
        \lr{8}\\
        \end{array}$}
        }
        sage: print(tex_from_array([[None,None,3],[None,5,6,7],[None,8]], with_lines=False))
        {\def\lr#1{\multicolumn{1}{@{\hspace{.6ex}}c@{\hspace{.6ex}}}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[b]{*{4}c}\\
        &&\lr{3}\\
        &\lr{5}&\lr{6}&\lr{7}\\
        &\lr{8}\\
        \end{array}$}
        }
        sage: Tableaux.options.convention="french"
        sage: print(tex_from_array([[1,2,3],[4,5]]))
        {\def\lr#1{\multicolumn{1}{|@{\hspace{.6ex}}c@{\hspace{.6ex}}|}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[t]{*{3}c}\cline{1-2}
        \lr{4}&\lr{5}\\\cline{1-3}
        \lr{1}&\lr{2}&\lr{3}\\\cline{1-3}
        \end{array}$}
        }
        sage: print(tex_from_array([[1,2,3],[4,5]], with_lines=False))
        {\def\lr#1{\multicolumn{1}{@{\hspace{.6ex}}c@{\hspace{.6ex}}}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[t]{*{3}c}\\
        \lr{4}&\lr{5}\\
        \lr{1}&\lr{2}&\lr{3}\\
        \end{array}$}
        }
        sage: print(tex_from_array([[1,2,3],[4,5,6,7],[8]]))
        {\def\lr#1{\multicolumn{1}{|@{\hspace{.6ex}}c@{\hspace{.6ex}}|}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[t]{*{4}c}\cline{1-1}
        \lr{8}\\\cline{1-4}
        \lr{4}&\lr{5}&\lr{6}&\lr{7}\\\cline{1-4}
        \lr{1}&\lr{2}&\lr{3}\\\cline{1-3}
        \end{array}$}
        }
        sage: print(tex_from_array([[1,2,3],[4,5,6,7],[8]], with_lines=False))
        {\def\lr#1{\multicolumn{1}{@{\hspace{.6ex}}c@{\hspace{.6ex}}}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[t]{*{4}c}\\
        \lr{8}\\
        \lr{4}&\lr{5}&\lr{6}&\lr{7}\\
        \lr{1}&\lr{2}&\lr{3}\\
        \end{array}$}
        }
        sage: print(tex_from_array([[None,None,3],[None,5,6,7],[8]]))
        {\def\lr#1{\multicolumn{1}{|@{\hspace{.6ex}}c@{\hspace{.6ex}}|}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[t]{*{4}c}\cline{1-1}
        \lr{8}\\\cline{1-4}
        &\lr{5}&\lr{6}&\lr{7}\\\cline{2-4}
        &&\lr{3}\\\cline{3-3}
        \end{array}$}
        }
        sage: print(tex_from_array([[None,None,3],[None,5,6,7],[None,8]]))
        {\def\lr#1{\multicolumn{1}{|@{\hspace{.6ex}}c@{\hspace{.6ex}}|}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[t]{*{4}c}\cline{2-2}
        &\lr{8}\\\cline{2-4}
        &\lr{5}&\lr{6}&\lr{7}\\\cline{2-4}
        &&\lr{3}\\\cline{3-3}
        \end{array}$}
        }
        sage: print(tex_from_array([[None,None,3],[None,5,6,7],[8]], with_lines=False))
        {\def\lr#1{\multicolumn{1}{@{\hspace{.6ex}}c@{\hspace{.6ex}}}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[t]{*{4}c}\\
        \lr{8}\\
        &\lr{5}&\lr{6}&\lr{7}\\
        &&\lr{3}\\
        \end{array}$}
        }
        sage: print(tex_from_array([[None,None,3],[None,5,6,7],[None,8]], with_lines=False))
        {\def\lr#1{\multicolumn{1}{@{\hspace{.6ex}}c@{\hspace{.6ex}}}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[t]{*{4}c}\\
        &\lr{8}\\
        &\lr{5}&\lr{6}&\lr{7}\\
        &&\lr{3}\\
        \end{array}$}
        }
        sage: Tableaux.options._reset()
    """
    lr=lr_macro.substitute(bar='|' if with_lines else '')
    if Tableaux.options.convention == "English":
        return '{%s\n%s\n}' % (lr, tex_from_skew_array(array, with_lines))
    else:
        return '{%s\n%s\n}' % (lr, tex_from_skew_array(array[::-1], with_lines, align='t'))

def tex_from_array_tuple(a_tuple, with_lines=True):
    r"""
    Return a latex string for a tuple of two dimensional array of partition,
    composition or skew composition shape.

    INPUT:

    - ``a_tuple`` -- a tuple of lists of lists
    - ``with_lines`` -- a boolean (default: ``True``)
      Whether to draw lines to separate the entries in the components of ``a_tuple``.

    .. SEEALSO:: :meth:`tex_from_array` for the description of each array

    EXAMPLES::

        sage: from sage.combinat.output import tex_from_array_tuple
        sage: print(tex_from_array_tuple([[[1,2,3],[4,5]],[],[[None,6,7],[None,8],[9]]]))
        {\def\lr#1{\multicolumn{1}{|@{\hspace{.6ex}}c@{\hspace{.6ex}}|}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[b]{*{3}c}\cline{1-3}
        \lr{1}&\lr{2}&\lr{3}\\\cline{1-3}
        \lr{4}&\lr{5}\\\cline{1-2}
        \end{array}$},\emptyset,\raisebox{-.6ex}{$\begin{array}[b]{*{3}c}\cline{2-3}
        &\lr{6}&\lr{7}\\\cline{2-3}
        &\lr{8}\\\cline{1-2}
        \lr{9}\\\cline{1-1}
        \end{array}$}
        }
        sage: print(tex_from_array_tuple([[[1,2,3],[4,5]],[],[[None,6,7],[None,8],[9]]], with_lines=False))
        {\def\lr#1{\multicolumn{1}{@{\hspace{.6ex}}c@{\hspace{.6ex}}}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[b]{*{3}c}\\
        \lr{1}&\lr{2}&\lr{3}\\
        \lr{4}&\lr{5}\\
        \end{array}$},\emptyset,\raisebox{-.6ex}{$\begin{array}[b]{*{3}c}\\
        &\lr{6}&\lr{7}\\
        &\lr{8}\\
        \lr{9}\\
        \end{array}$}
        }
        sage: Tableaux.options.convention="french"
        sage: print(tex_from_array_tuple([[[1,2,3],[4,5]],[],[[None,6,7],[None,8],[9]]]))
        {\def\lr#1{\multicolumn{1}{|@{\hspace{.6ex}}c@{\hspace{.6ex}}|}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[t]{*{3}c}\cline{1-2}
        \lr{4}&\lr{5}\\\cline{1-3}
        \lr{1}&\lr{2}&\lr{3}\\\cline{1-3}
        \end{array}$},\emptyset,\raisebox{-.6ex}{$\begin{array}[t]{*{3}c}\cline{1-1}
        \lr{9}\\\cline{1-2}
        &\lr{8}\\\cline{2-3}
        &\lr{6}&\lr{7}\\\cline{2-3}
        \end{array}$}
        }
        sage: print(tex_from_array_tuple([[[1,2,3],[4,5]],[],[[None,6,7],[None,8],[9]]], with_lines=False))
        {\def\lr#1{\multicolumn{1}{@{\hspace{.6ex}}c@{\hspace{.6ex}}}{\raisebox{-.3ex}{$#1$}}}
        \raisebox{-.6ex}{$\begin{array}[t]{*{3}c}\\
        \lr{4}&\lr{5}\\
        \lr{1}&\lr{2}&\lr{3}\\
        \end{array}$},\emptyset,\raisebox{-.6ex}{$\begin{array}[t]{*{3}c}\\
        \lr{9}\\
        &\lr{8}\\
        &\lr{6}&\lr{7}\\
        \end{array}$}
        }
        sage: Tableaux.options._reset()
    """
    lr=lr_macro.substitute(bar='|' if with_lines else '')
    if Tableaux.options.convention == "English":
        return '{%s\n%s\n}' % (lr, ','.join(
            r'\emptyset' if comp==[] else tex_from_skew_array(comp, with_lines) for comp in a_tuple))
    else:
        return '{%s\n%s\n}' % (lr, ','.join(
            r'\emptyset' if comp==[] else tex_from_skew_array(comp[::-1], with_lines, align='t') for comp in a_tuple))

def tex_from_skew_array(array, with_lines=False, align='b'):
    r"""
    This function creates latex code for a "skew composition" ``array``.
    That is, for a two dimensional array in which each row can begin with
    an arbitrary number ``None``'s and the remaining entries could, in
    principle, be anything but probably should be strings or integers of similar
    width. A row consisting completely of ``None``'s is allowed.

    INPUT:

    - ``array`` -- The array

    - ``with_lines`` -- (Default: ``False``) If ``True`` lines are drawn, if
      ``False`` they are not

    - ``align`` -- (Default: ``'b'``) Determines the alignment on the latex
      array environments

    EXAMPLES::

        sage: array=[[None, 2,3,4],[None,None],[5,6,7,8]]
        sage: print(sage.combinat.output.tex_from_skew_array(array))
        \raisebox{-.6ex}{$\begin{array}[b]{*{4}c}\\
        &\lr{2}&\lr{3}&\lr{4}\\
        &\\
        \lr{5}&\lr{6}&\lr{7}&\lr{8}\\
        \end{array}$}

    TESTS::

        sage: sage.combinat.output.tex_from_skew_array([(1,2,3), (2,3,4)])
        '\\raisebox{-.6ex}{$\\begin{array}[b]{*{3}c}\\\\\n\\lr{1}&\\lr{2}&\\lr{3}\\\\\n\\lr{2}&\\lr{3}&\\lr{4}\\\\\n\\end{array}$}'
        sage: sage.combinat.output.tex_from_skew_array([((1,2,),)])
        '\\raisebox{-.6ex}{$\\begin{array}[b]{*{1}c}\\\\\n\\lr{(1, 2)}\\\\\n\\end{array}$}'
    """
    # first identify where the None's appear in ``array`` and define a
    # function end_line which puts in the required \cline's.
    if with_lines:
        # last position of None in each row
        nones=[1 if not None in row else 1+len(row)-row[::-1].index(None) for row in array]
        def end_line(r):
            # in a slightly unpythonic way, we label the lines as 0, 1, ..., len(array)
            if r==0:
                return r'\cline{%s-%s}'%(nones[0],len(array[0]))
            elif r==len(array):
                start=nones[r-1]
                finish=len(array[r-1])
            else:
                start=min(nones[r], nones[r-1])
                finish=max(len(array[r]), len(array[r-1]))
            return r'\\' if start>finish else r'\\\cline{%s-%s}'%(start, finish)
    else:
        end_line=lambda r: r'\\'

    # now we draw the array
    tex=r'\raisebox{-.6ex}{$\begin{array}[%s]{*{%s}c}'%(align,max(map(len,array)))
    tex+=end_line(0)+'\n'
    for r in range(len(array)):
        tex+='&'.join('' if c is None else r'\lr{%s}'%(c,) for c in array[r])
        tex+=end_line(r+1)+'\n'
    return tex+r'\end{array}$}'

def ascii_art_table(data, use_unicode=False, convention="English"):
    r"""
    Return an ascii art table of ``data``.

    EXAMPLES::

        sage: from sage.combinat.output import ascii_art_table

        sage: data = [[None, None, 1], [2, 2], [3,4,5], [None, None, 10], [], [6]]
        sage: print(ascii_art_table(data))
                +----+
                | 1  |
        +---+---+----+
        | 2 | 2 |
        +---+---+----+
        | 3 | 4 | 5  |
        +---+---+----+
                | 10 |
                +----+
        <BLANKLINE>
        +---+
        | 6 |
        +---+
        sage: print(ascii_art_table(data, use_unicode=True))
                ┌────┐
                │ 1  │
        ┌───┬───┼────┘
        │ 2 │ 2 │
        ├───┼───┼────┐
        │ 3 │ 4 │ 5  │
        └───┴───┼────┤
                │ 10 │
                └────┘
        <BLANKLINE>
        ┌───┐
        │ 6 │
        └───┘

        sage: data = [[1, None, 2], [None, 2]]
        sage: print(ascii_art_table(data))
        +---+   +---+
        | 1 |   | 2 |
        +---+---+---+
            | 2 |
            +---+
        sage: print(ascii_art_table(data, use_unicode=True))
        ┌───┐   ┌───┐
        │ 1 │   │ 2 │
        └───┼───┼───┘
            │ 2 │
            └───┘
    """
    if use_unicode:
        import unicodedata
        v  = unicodedata.lookup('BOX DRAWINGS LIGHT VERTICAL')
        h  = unicodedata.lookup('BOX DRAWINGS LIGHT HORIZONTAL')
        dl = unicodedata.lookup('BOX DRAWINGS LIGHT DOWN AND LEFT')
        dr = unicodedata.lookup('BOX DRAWINGS LIGHT DOWN AND RIGHT')
        ul = unicodedata.lookup('BOX DRAWINGS LIGHT UP AND LEFT')
        ur = unicodedata.lookup('BOX DRAWINGS LIGHT UP AND RIGHT')
        vr = unicodedata.lookup('BOX DRAWINGS LIGHT VERTICAL AND RIGHT')
        vl = unicodedata.lookup('BOX DRAWINGS LIGHT VERTICAL AND LEFT')
        uh = unicodedata.lookup('BOX DRAWINGS LIGHT UP AND HORIZONTAL')
        dh = unicodedata.lookup('BOX DRAWINGS LIGHT DOWN AND HORIZONTAL')
        vh = unicodedata.lookup('BOX DRAWINGS LIGHT VERTICAL AND HORIZONTAL')
        from sage.typeset.unicode_art import unicode_art as art
    else:
        v = '|'
        h = '-'
        dl = dr = ul = ur = vr = vl = uh = dh = vh = '+'
        from sage.typeset.ascii_art import ascii_art as art

    if not data:
        return dr + dl + '\n' + ur + ul

    # Convert the input into a rectangular array with the top and bottom row
    #   being all None's for ease later on.
    ncols = max(len(row) for row in data)
    str_tab = [[None]*ncols] + [[art(val) if val is not None else None for val in row] + [None]*(ncols-len(row))
                                for row in data]
    str_tab.append([None]*ncols)
    # Get the widths of the columns
    col_widths = [1]*len(str_tab[0])
    if use_unicode:
        # Special handling of overline not adding to printed length
        def get_len(e):
            if e is None:
                return 0
            return len(e) - list(str(e)).count(u"\u0304")
    else:
        def get_len(e):
            if e is None:
                return 0
            return len(e)
    for row in str_tab:
        for i,e in enumerate(row):
            col_widths[i] = max(col_widths[i], get_len(e))

    matr = []  # just the list of lines
    for nrow,row in enumerate(str_tab):
        if nrow == 0: # skip the first row
            continue

        l1 = ""
        l2 = ""
        for i,(e,w) in enumerate(zip(row, col_widths)):
            prev_row = str_tab[nrow-1]
            if i == 0:
                if e is None:
                    if prev_row[i] is None:
                        l1 += " "*(3+w)
                    else:
                        l1 += ur + h*(2+w)
                    l2 += " "*(3+w)
                else:
                    if prev_row[i] is None:
                        l1 += dr + h*(2+w)
                    else:
                        l1 += vr + h*(2+w)
                    l2 += "{} {:^{width}} ".format(v, e, width=w)
            else:
                if e is None:
                    if row[i-1] is None:
                        if prev_row[i-1] is None:
                            if prev_row[i] is None:
                                l1 += " "*(3+w)
                            else:
                                l1 += ur + h*(2+w)
                        else:
                            if prev_row[i] is None:
                                l1 += ul + " "*(2+w)
                            else:
                                l1 += uh + h*(2+w)
                        l2 += " "*(3+w)
                    else:
                        if prev_row[i-1] is None:
                            if prev_row[i] is None:
                                l1 += dl + " "*(2+w)
                            else:
                                l1 += vh + h*(2+w)
                        else:
                            if prev_row[i] is None:
                                l1 += vl + " "*(2+w)
                            else:
                                l1 += vh + h*(2+w)
                        l2 += v + " "*(2+w)
                else:
                    if row[i-1] is None:
                        if prev_row[i-1] is None:
                            if prev_row[i] is None:
                                l1 += dr + h*(2+w)
                            else:
                                l1 += vr + h*(2+w)
                        else:
                            l1 += vh + h*(2+w)
                    else:
                        if prev_row[i-1] is None and prev_row[i] is None:
                                l1 += dh + h*(2+w)
                        else:
                            l1 += vh + h*(2+w)
                    l2 += "{} {:^{width}} ".format(v, e, width=w)

        if row[-1] is None:
            if prev_row[-1] is None:
                l1 += " "
            else:
                l1 += ul
            l2 += " "
        else:
            if prev_row[-1] is None:
                l1 += dl
            else:
                l1 += vl
            l2 += v

        matr.append(l1)
        matr.append(l2)

    matr.pop()  # Remove the last row (which is blank)

    if convention == "English":
        return "\n".join(matr)
    else:
        output = "\n".join(reversed(matr))
        if use_unicode:
            tr = {
                ord(dl): ul, ord(dr): ur,
                ord(ul): dl, ord(ur): dr,
                ord(dh): uh, ord(uh): dh}
            return output.translate(tr)
        else:
            return output

max_forks_repo_name: dimpase/sage

reason: from sage


text:
"""
Reorders factors in a manual fashion.
"""

import numpy as np
import pandas as pd

from click import *
from logging import *
from sklearn.decomposition import NMF
from sklearn.externals import joblib

@command()
@option(
    '--model-input',
    required=True,
    help='the Pickle file to read the model from')
@option(
    '--basis-input',
    required=True,
    help='the CSV file to read the basis matrix from')
@option(
    '--score-input', required=True, help='the CSV file to read scores from')
@option(
    '--model-output',
    required=True,
    help='the Pickle file to write the model to')
@option(
    '--basis-output',
    required=True,
    help='the CSV file to write the basis matrix to')
@option(
    '--score-output', required=True, help='the CSV file to write scores to')
@option(
    '--at-end',
    type=int,
    multiple=True,
    help='the factors to move to the end (multiple permitted)')
def main(model_input, basis_input, score_input, model_output, basis_output,
         score_output, at_end):

    basicConfig(
        level=INFO,
        handlers=[
            StreamHandler(), FileHandler(
                '{}.log'.format(model_output), mode='w')
        ])

    # Load the model.

    info('Loading model')

    model = joblib.load(model_input)

    info('Loading basis matrix')

    basis = pd.read_csv(basis_input, index_col=0)

    info('Result: {}'.format(basis.shape))

    info('Loading scores')

    scores = pd.read_csv(score_input, index_col=0)

    info('Result: {}'.format(scores.shape))

    # Generate a new ordering for the factors.

    info('Reordering factors')

    factor_order = np.arange(model.components_.shape[0], dtype=int) + 1

    if at_end:

        at_end = np.array(at_end, dtype=int)

        factor_order = np.concatenate(
            [np.setdiff1d(factor_order, at_end), at_end])

    model.components_ = model.components_[factor_order - 1]

    basis = basis.iloc[:, factor_order - 1]

    basis.columns = factor_order

    scores = scores.iloc[:, factor_order - 1]

    scores.columns = factor_order

    # Write the output.

    info('Writing output')

    joblib.dump(model, model_output)

    basis.to_csv(basis_output)

    scores.to_csv(score_output)

if __name__ == '__main__':
    main()

max_stars_repo_name: morrislab/plos-medicine-joint-patterns

max_issues_repo_name: morrislab/plos-medicine-joint-patterns

max_forks_repo_name: morrislab/plos-medicine-joint-patterns


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

## Imports ##

import os
import pickle
import numpy as np
np.random.seed(2018)

import sys
sys.path.append(os.path.join(os.path.expanduser('~'), 'find-tune', 'train'))
import perf_utils as pu

from imblearn.over_sampling import SMOTE
from sklearn.preprocessing  import StandardScaler
from sklearn.model_selection import StratifiedKFold

from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier

from keras.layers import Input, Dense
from keras.models import Model, load_model
from keras.optimizers import Adam
from keras.callbacks import ReduceLROnPlateau, EarlyStopping, ModelCheckpoint

from definitions import DATA_DIR

## Main ##
print('=== Preparing data ===')

# Load data from file
# X = embedding features (output from VGGish model)
# y = binary target variable (1 = target track, 0 = other sound)
# c = sound categories (e.g. 'Jingle Bells', 'target')
# s = seconds in the original YouTube track where this clip begins
# ids = YouTube unique identifiers for the track
# L = log-mel attributes (input to VGGish model)
data = np.load(os.path.join(DATA_DIR, 'embedding_data.npz'))
X, y, c, s, ids, L = data['X'], data['y'], data['c'], data['s'], data['i'], data['L']

# Train (80%) / val (20%) split
p = np.mean(y)
N = int(len(y)*0.2)
print('p: %.4f' % p) # ~ 7.5/1000
print('N: %d'   % N) # ~ 700

# Note that whilst the mean of Binomial(N, p) is ~5, there is a very real chance that
# we could have very few positive data points in the validation set if we're not careful
print('P(X < 3) = %.3f (where X ~ Bin(N, p))' % np.mean(np.random.binomial(N, p, size=100000) < 3))

# So, to ensure this doesn't happen, we will explicitly hold out 5 positive examples for testing
i_pos = np.where(y == 1)[0]
i_neg = np.where(y == 0)[0]

i_pos_va = np.random.choice(i_pos, size=5, replace=False) # 5 for validation
i_pos_tr = np.array(list(set(i_pos) - set(i_pos_va)))     # compliment for training

i_neg_tr = np.random.choice(i_neg, size=len(X)-N, replace=False) # Training
i_neg_va = np.array(list(set(i_neg) - set(i_neg_tr)))            # Validation

# Create split datasets
X_tr = X[np.append(i_neg_tr, i_pos_tr), :]
y_tr = y[np.append(i_neg_tr, i_pos_tr)]
c_tr = c[np.append(i_neg_tr, i_pos_tr)]
s_tr = s[np.append(i_neg_tr, i_pos_tr)]
ids_tr = ids[np.append(i_neg_tr, i_pos_tr)]
L_tr = L[np.append(i_neg_tr, i_pos_tr), :, :]

X_va = X[np.append(i_neg_va, i_pos_va), :]
y_va = y[np.append(i_neg_va, i_pos_va)]
c_va = c[np.append(i_neg_va, i_pos_va)]
s_va = s[np.append(i_neg_va, i_pos_va)]
ids_va = ids[np.append(i_neg_va, i_pos_va)]
L_va = L[np.append(i_neg_va, i_pos_va), :, :]

print('Positive labels in training set:   %d' % y_tr.sum())
print('Positive labels in validation set: %d' % y_va.sum())

print('Negative labels in training set:   %d' % (y_tr==0).sum())
print('Negative labels in validation set: %d' % (y_va==0).sum())

# As this is quite an imbalanced problem, we'll use SMOTE to over sample the positive class
sm = SMOTE(random_state=2018)
X_tr_bal, y_tr_bal = sm.fit_sample(X_tr, y_tr)

print('Positive labels in balanced training set: %d' % y_tr_bal.sum())
print('Negative labels in balanced training set: %d' % (y_tr_bal==0).sum())

# Centre and scale the features
ss = StandardScaler()
X_tr_bal = ss.fit_transform(X_tr_bal)
X_tr = ss.transform(X_tr)
X_va = ss.transform(X_va)

## Linear classifier ##
print('\n')
print('=== Linear classifier ===')
print()
# Yields accuracy and perfect recall on positive class
# Suffers from imperfect precision on positive class (gives false positives)

# Fit multiple classifiers with increasing regularization strength
# Print a report card for each model
for alpha in np.logspace(-4, -4, 1):
    print('alpha = %.e' % alpha)
    lr = LogisticRegression(C=1/alpha, solver='liblinear')
    lr.fit(X_tr_bal, y_tr_bal)
    y_pred_tr_bal = lr.predict(X_tr_bal)
    y_pred_tr = lr.predict(X_tr)
    y_pred_va = lr.predict(X_va)

    pu.print_scorecard(y_tr_bal, y_pred_tr_bal, title='TRAIN (BAL.)')
    pu.print_scorecard(y_tr, y_pred_tr, title='TRAIN')
    pu.print_scorecard(y_va, y_pred_va, title='VALIDATION')
    print()

# An analysis of incorrect predictions
pu.print_negatives(y_va, y_pred_va, c_va, ytids=ids_va, num_secs=s_va)

## Random Forest ##
print('\n')
print('=== Random forest classifier ===')
print()
# Performs quite well, good accuracy
# Perfect precision on positive class (desirable: don't want false positives)
# Note: overfits without cap on depth
rf = RandomForestClassifier(n_estimators = 50, max_depth = 5)
#rf.fit(X_tr_bal, y_tr_bal)
rf.fit(X_tr, y_tr)
y_pred_tr_bal = rf.predict(X_tr_bal)
y_pred_tr = rf.predict(X_tr)
y_pred_va = rf.predict(X_va)

pu.print_scorecard(y_tr_bal, y_pred_tr_bal, title='TRAIN (BAL.)')
pu.print_scorecard(y_tr, y_pred_tr, title='TRAIN')
pu.print_scorecard(y_va, y_pred_va, title='VALIDATION')
print()

pu.print_negatives(y_va, y_pred_va, c_va, ytids=ids_va, num_secs=s_va)

## Dense Neural Network ##
print('\n')
print('=== Neural network classifier ===')
print()
# Without any hyperparameter tuning, the network gets a perfect score on all metrics!

def fit_nn_model(lr0, h1, bsz, verbose=0, cp_path=None):
    """
    Generic model for fitting a 1-layer neural network
    """

    # Callbacks
    lr_red = ReduceLROnPlateau(monitor='val_loss', min_delta=0, factor=0.75, patience=2)
    e_stop = EarlyStopping(monitor='val_loss', min_delta=0, patience=15)
    cbacks = [lr_red, e_stop]
    if cp_path is not None:
        cbacks.append(ModelCheckpoint(cp_path, save_best_only=True))

    # Model architecture
    ix = Input((X_tr_bal.shape[1],), name='vggish_feat_input')
    x  = Dense(h1, activation='relu', name='fc_last')(ix)
    x  = Dense(1, activation='sigmoid', name='classify')(x)

    # Compile the Model
    model = Model(ix, x)
    optzr = Adam(lr=lr0)
    model.compile(optzr, loss='binary_crossentropy', metrics=['accuracy'])
    # Pre-training metrics
    ep0_tr = model.evaluate(X_tr_bal, y_tr_bal, verbose=verbose)
    ep0_va = model.evaluate(X_va, y_va, verbose=verbose)
    # Fit the model
    history = model.fit(X_tr_bal, y_tr_bal,
                        epochs=200, batch_size=bsz,
                        callbacks=cbacks,
                        validation_data=(X_va, y_va),
                        verbose=verbose,
                        )

    return (model, history, (ep0_tr, ep0_va))

# Hyperparameters
lr0 = 0.01
h1  = 128
bsz = 64

model, history, _ = fit_nn_model(lr0, h1, bsz, verbose=0)

# Make (probability) predictions
y_pred_tr_bal = model.predict(X_tr_bal)[:, 0]
y_pred_tr = model.predict(X_tr)[:, 0]
y_pred_va = model.predict(X_va)[:, 0]

p_thresh = 0.5
pu.print_scorecard(y_tr_bal, y_pred_tr_bal > p_thresh, title='TRAIN (BAL.)')
pu.print_scorecard(y_tr, y_pred_tr > p_thresh, title='TRAIN')
pu.print_scorecard(y_va, y_pred_va > p_thresh, title='VALIDATION')

# More rigorous test: use cross-validation, 5-folds
# Again, perfect scorecards for each fold - shuffling is important for diversity of training data!
print('\n')
print('=== 5-fold CV ===')
print()
histories = []
fold_cntr = 0
skf = StratifiedKFold(5, shuffle=True, random_state=2018)
for i_tr, i_va in skf.split(X, y):

    # Storage for data / models relating to this fold
    fold_dir = os.path.join(DATA_DIR, 'fold%d' % fold_cntr)
    os.makedirs(fold_dir, exist_ok=True)

    # Create split datasets
    X_tr = X[i_tr, :]
    y_tr = y[i_tr]
    c_tr = c[i_tr]
    s_tr = s[i_tr]
    ids_tr = ids[i_tr]
    L_tr = L[i_tr, :, :]

    X_va = X[i_va, :]
    y_va = y[i_va]
    c_va = c[i_va]
    s_va = s[i_va]
    ids_va = ids[i_va]
    L_va = L[i_va, :, :]

    print('FOLD %d' % fold_cntr)
    print('======')
    print()
    print('Positive labels in training set:   %d' % y_tr.sum())
    print('Positive labels in validation set: %d' % y_va.sum())

    # over sample with SMOTE
    sm = SMOTE(random_state=2018)
    X_tr_bal, y_tr_bal = sm.fit_sample(X_tr, y_tr)
    print('Positive labels in balanced set:   %d' % y_tr_bal.sum())

    print()
    print('=== Linear classifier ===')
    alpha = 1e-4
    lr = LogisticRegression(C=1/alpha, solver='liblinear')
    lr.fit(X_tr_bal, y_tr_bal)
    y_pred_va = lr.predict(X_va)
    pu.print_scorecard(y_va, y_pred_va, title='VALIDATION')
    pu.print_negatives(y_va, y_pred_va > p_thresh, c_va, ytids=ids_va, num_secs=s_va)

    # Centre and scale the features
    ss = StandardScaler()
    X_tr_bal = ss.fit_transform(X_tr_bal)
    X_tr = ss.transform(X_tr)
    X_va = ss.transform(X_va)
    with open(os.path.join(fold_dir, 'sc_fold%d.pkl' % fold_cntr), 'wb') as f:
        pickle.dump(ss, f)

    # Save data (with input scaled features) for later reference
    np.savez(os.path.join(fold_dir, 'foldwise_data_tr.npz'), X=X_tr, y=y_tr, c=c_tr, s=s_tr, i=ids_tr, L=L_tr)
    np.savez(os.path.join(fold_dir, 'foldwise_data_va.npz'), X=X_va, y=y_va, c=c_va, s=s_va, i=ids_va, L=L_va)
    np.savez(os.path.join(fold_dir, 'foldwise_data_bal.npz'), X=X_tr_bal, y=y_tr_bal)

    # Fit the model
    save_path = os.path.join(fold_dir, 'nn_fold%d.model' % fold_cntr)
    _, history, _ = fit_nn_model(lr0, h1, bsz, verbose=0, cp_path=save_path)
    histories.append(history)

    # Load the best model and predict
    print()
    print('=== Dense neural network ===')
    model = load_model(save_path)
    y_pred_va = model.predict(X_va)[:, 0]
    pu.print_scorecard(y_va, y_pred_va > p_thresh, title='VALIDATION')
    pu.print_negatives(y_va, y_pred_va > p_thresh, c_va, ytids=ids_va, num_secs=s_va)

    print()
    fold_cntr += 1

max_stars_repo_name: lukerm/find-tune

max_issues_repo_name: lukerm/find-tune

max_forks_repo_name: lukerm/find-tune


text:
import os
import csv
import sys
import argparse
import numpy as np
import tensorflow as tf
from tensorflow import GradientTape
from tensorflow_addons.optimizers import AdamW
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.losses import MSE

class TrainPose():
    def __init__(self, args):
        # initialize tensorflow
        gpu = tf.config.experimental.list_physical_devices('GPU')[0]
        tf.config.experimental.set_memory_growth(gpu, True)
        sys.setrecursionlimit(10000)

        print("Eager execution is " + str(tf.executing_eagerly()))
        print("Tensorflow version " + str(tf.__version__))
        # import classes
        sys.path.append("..")
        from models.feature_pose import FeaturePose
        from preprocessing.process_kaggle import ProcessKaggle
        # initialize training variables
        self.numEpochs = args.numEpochs
        self.verbose = args.verbose

        # self.optimizer = AdamW(1e-2 * args.lr, learning_rate=args.lr)
        self.optimizer = Adam(args.lr)

        self.loss = None

        self.restore = args.restore

        # initialize preprocessor and build/load data
        preprocessor = ProcessKaggle(args.rawDir, args.proDir,
                                     args.batchSize, args.threadNum,
                                     args.preFetch, args.valSize)

        self.trainData, self.valData = preprocessor.loadData()

        # initialize models

        self.model = FeaturePose()

        # initialize checkpoints and managers
        if not self.restore:
            for f in os.listdir('../data/checkpoints/feature_pose'):
                os.remove('../data/checkpoints/feature_pose/' + f)
            if os.path.isfile('../data/checkpoints/feature_pose_losses.csv'):
                os.remove('../data/checkpoints/feature_pose_losses.csv')

        self.checkpoint = tf.train.Checkpoint(optimizer=self.optimizer,
                                              model=self.model)
        self.manager = tf.train.CheckpointManager(self.checkpoint,
                                                  directory='../data/checkpoints/feature_pose',
                                                  max_to_keep=5, keep_checkpoint_every_n_hours=2)

    def _loss(self, x, y):
        return MSE(y, x)

    def _logTrain(self, epoch, batch):
        msg = 'Training Batch {} in Epoch {}, Model Loss: {}'
        print(msg.format(batch + 1, epoch + 1, self.loss), end='\r')

    def _logVal(self, epoch, batch):
        msg = 'Validating Batch {} in Epoch {}, Model Loss: {}'
        print(msg.format(batch + 1, epoch + 1, self.loss), end='\r')

    def _save(self):
        self.manager.save()
        with open('../data/checkpoints/feature_pose_losses.csv', 'a', newline='') as file:
            writer = csv.writer(file)
            writer.writerow([self.lossTrainAverage, self.lossValAverage])

    @tf.function
    def _validate(self, input, output):
        predicted = self.model.call(input)
        loss = self._loss(predicted, output)

        return tf.math.reduce_mean(loss)

    @tf.function
    def _update(self, input, output):
        with GradientTape() as tape:
            predicted = self.model.call(input, training=True)
            loss = self._loss(predicted, output)

            gradients = tape.gradient(loss, self.model.trainable_variables)
            self.optimizer.apply_gradients(zip(gradients, self.model.trainable_variables))
            return tf.math.reduce_mean(loss)

    def train(self):
        if self.restore:
            self.checkpoint.restore(self.manager.latest_checkpoint)

        for epoch in range(self.numEpochs):
            self.lossTrainAverage = list()
            self.lossValAverage = list()
            for ind, (input, output) in self.trainData.enumerate():
                self.loss = self._update(input, output).numpy()
                self.lossTrainAverage.append(self.loss)
                if self.verbose:
                    self._logTrain(epoch, ind)
            self.trainData = self.trainData.shuffle(10000)
            for ind, (input, output) in self.valData.enumerate():
                self.loss = self._validate(input, output).numpy()
                self.lossValAverage.append(self.loss)
                if self.verbose:
                    self._logVal(epoch, ind)

            self.lossTrainAverage = np.average(np.asarray(self.lossTrainAverage))
            self.lossValAverage = np.average(np.asarray(self.lossValAverage))
            print('Epoch {} of {} complete with training loss of {} and validation loss of {}'.format(
                epoch + 1, self.numEpochs, self.lossTrainAverage, self.lossValAverage))
            self.valData = self.valData.shuffle(10000)

            self._save()

if __name__ == '__main__':
    def str2bool(v):
        return v.lower() in ("yes", "true", "t", "1")

    parser = argparse.ArgumentParser()
    parser.add_argument('--numEpochs', default=100, type=int)
    parser.add_argument('--lr', default=1e-4, type=float)
    parser.add_argument('--rawDir', default='../data/raw/kaggle', type=str)
    parser.add_argument('--proDir', default='../data/processed/kaggle', type=str)
    parser.add_argument('--batchSize', default=32, type=int)
    parser.add_argument('--preFetch', default=1, type=int)
    parser.add_argument('--valSize', default=1000, type=int)
    parser.add_argument('--threadNum', default=4, type=int)
    parser.add_argument('--restore', default=True, type=str2bool)
    parser.add_argument('--verbose', default=True, type=str2bool)
    args = parser.parse_args()

    train = TrainPose(args)

    try:
        train.train()
    except KeyboardInterrupt:
        print('')
        print("Training stopped. Continue by running again with --restore True. Note: will reset weight decay for AdamW optimizer")

max_stars_repo_name: aaronlockhartdev/deepfakeGAN

max_issues_repo_name: aaronlockhartdev/deepfakeGAN

max_forks_repo_name: aaronlockhartdev/deepfakeGAN


text:
import time
import json
import numpy as np
import quote_pb2

symprice = quote_pb2.SymbolPrice()
symvolume = quote_pb2.SymbolVolume()
price_dict = {}
volume_dict = {}

symprice.code = 0
for i in range(50000):
    key = str(i).zfill(6)
    symprice.data[key] = price_dict[key] = np.random.random()

symvolume.code = 0
for i in range(50000):
    key = str(i).zfill(6)
    symvolume.data[key] = price_dict[key] = np.random.randint(1, 10000)

t1= time.time()
out1 = symprice.SerializeToString()
t2 = time.time()
print("dumps time", t2 - t1)

t1 = time.time()
out2 = symvolume.SerializeToString()
t2 = time.time()
print("dumps time", t2 - t1)

t1 = time.time()
data = quote_pb2.SymbolPrice()
data.ParseFromString(out1)
t2 = time.time()
print(data.data['000001'])
print("load time", t2 - t1, data.code)

t1 = time.time()
data = quote_pb2.SymbolVolume()
data.ParseFromString(out2)
t2 = time.time()
print(data.data['000001'], data.code)
print("load time", t2 - t1)

t1 = time.time()
d=json.dumps(price_dict)
t2 = time.time()

print("json dumps", t2 - t1)
t1 = time.time()
json.loads(d)
t2 = time.time()
print("json loads", t2 - t1)

max_stars_repo_name: colin-zhou/mrfs

max_issues_repo_name: colin-zhou/mrfs

max_forks_repo_name: colin-zhou/mrfs


text:
import os.path as op
import numpy as np
import pandas as pd
import itertools
import numpy.testing as npt
import pytest

import matplotlib
matplotlib.use('Agg')

import tcrdist as td

tempSkip = pytest.mark.skip(reason="Temporarily skipping for efficiency.")

@tempSkip
def exampleBetaChains():
    lines = """TRBV18*01 CASSPRHGISPLHF
    TRBV19*01,TRBV19*02 CASSPGGVTEAFF
    TRBV28*01 CASRGPGEPYEQYF
    TRBV4-1*01,TRBV4-1*02 CASSQAEGGLLSYEQYF
    TRBV6-6*01,TRBV6-6*02,TRBV6-6*03,TRBV6-6*04 CASMLGGPPQETQYF
    TRBV7-8*03 CASSHSWDSGTGELFF
    TRBV7-9*01,TRBV7-9*02,TRBV7-9*03 CASSSSGGAYNEQFF
    TRBV19*01,TRBV19*02 CASSIEGSGANVLTF
    TRBV7-8*01 CASSSSGGMNIQYF
    TRBV6-1*01 CASSEAPLLGGNEQYF
    TRBV29-1*01 CSVEARLF""".split('\n')

    tcrs = []
    for line in lines:
        l = line.split()
        tcr = td.objects.TCRChain()
        tcr.cdr3b = l[1]
        tcr.vb_genes = l[0].split(',')
        tcr.organism = 'human'
        tcr.vb_reps = ';'.join([td.util.get_rep(vb_gene, 'human') for vb_gene in tcr.vb_genes])
        tcrs.append(tcr)
    return tcrs

datasetsPath = op.join(td.__path__[0], 'datasets')

tempSkip = pytest.mark.skip(reason="Temporarily skipping for efficiency.")

@tempSkip
def test_human_vregion():
    """Getting lots of NA distances..."""
    pwDf = td.distances.computeVRegionDistances('human')
    assert np.all(pwDf.index == pwDf.columns)

@tempSkip
def test_mouse_vregion():
    pwDf = td.distances.computeVRegionDistances('mouse')
    assert np.all(pwDf.index == pwDf.columns)

@tempSkip
def test_chain():
    tcrs = exampleBetaChains()
    pwmat = np.zeros((len(tcrs), len(tcrs)))
    for i, j in itertools.product(list(range(len(tcrs))), list(range(len(tcrs)))):
        if j <= i:
            pwmat[i, j] = td.distances.basicDistance('B', tcrs[i], tcrs[j])
            pwmat[j, i] = pwmat[i, j]

@tempSkip
def test_compute_paired_chain_distance():
    psDf = td.datasets.loadPSData('test_human_pairseqs')
    d = td.distances.basicDistance('AB', psDf.iloc[0], psDf.iloc[-1])

@tempSkip
def test_compute_clones_distance():
    psDf = td.datasets.loadPSData('test_human_pairseqs')
    probDf = td.processing.computeProbs(psDf)
    psDf = psDf.join(probDf)
    clonesDf = td.processing.identifyClones(psDf)
    d = td.distances.basicDistance('AB', clonesDf.iloc[0], clonesDf.iloc[-1])

@tempSkip
def test_compute_all_distances():
    psDf = td.datasets.loadPSData('test_human_pairseqs')
    probDf = td.processing.computeProbs(psDf)
    psDf = psDf.join(probDf)
    clonesDf = td.processing.identifyClones(psDf)
    clonesDf = clonesDf.set_index('CLONEID')
    pwDf = td.distances.computeBasicPWDistances('AB', clonesDf)
    assert pwDf.shape[0] == pwDf.shape[1]
    assert pwDf.shape[0] == clonesDf.shape[0]

@tempSkip
def test_nbr_dist():
    psDf = td.datasets.loadPSData('test_human_pairseqs')
    psDf = psDf.set_index('TCRID')
    pwDf = td.distances.computeBasicPWDistances('AB', psDf)
    params = pd.objects.DistanceParams()
    VRegionDists = td.distance.computeVRegionDistances('human')

    d = nearestNeighborDistance(chains='AB',
                                tcr=psDf.iloc[0],
                                referenceTCRs=psDf,
                                neighborTopN=5,
                                pwDf=pwDf,
                                weighted=False,
                                VRegionDists=VRegionDists,
                                params=params)
    dw = nearestNeighborDistance(chains='AB',
                                tcr=psDf.iloc[0],
                                referenceTCRs=psDf,
                                neighborTopN=5,
                                pwDf=pwDf,
                                weighted=True,
                                VRegionDists=VRegionDists,
                                params=params)
    """Force recompute with pwDf = None"""
    d2 = nearestNeighborDistance(chains='AB',
                                tcr=psDf.iloc[0],
                                referenceTCRs=psDf,
                                neighborTopN=5,
                                pwDf=None,
                                weighted=False,
                                VRegionDists=VRegionDists,
                                params=params)
    assert d == d2

@tempSkip
def test_embedding():
    psDf = td.datasets.loadPSData('test_human_pairseqs')
    psDf = psDf.set_index('TCRID')
    pwDf = td.distances.computeBasicPWDistances('AB', psDf)

    # xyDf = embedDistanceMatrix(pwDf, method='kpca', n_components=2)
    xyDf = td.embedding.plotEmbedding(pwDf, labels=psDf.epitope, method='kpca')


text:
import numpy as np
import pandas as pd

class Scenes:
  def __init__(self, df, users):
    self.globalDf = df.loc[:, df.columns.intersection(['user','hash','emotion'])]
    self.data = self.globalDf.loc[self.globalDf.user.isin(users), :].to_numpy()
    hash = np.unique(self.data.transpose()[2])
    self.lookup = {ith:np.where(self.data == ith) for ith in hash}

  def getData(self, hash):
    # use as mask
    try:
      idx = self.lookup[hash][0]
      return self.data[idx].transpose(1,0)[:2]
    except:
      return np.array([ np.array([]),np.array([]) ])

  def updateUsers(self, users):
    self.data = self.globalDf.loc[self.globalDf.user.isin(users), :].to_numpy()
    hash = np.unique(self.data.transpose()[2])
    self.lookup = {ith:np.where(self.data == ith) for ith in hash}

  def __getitem__(self, attr):
    return self.getData(attr)

max_stars_repo_name: ryanapfel/clustering

max_issues_repo_name: ryanapfel/clustering

max_forks_repo_name: ryanapfel/clustering


text:
import torch
import torch.nn as nn
import torch.optim as optim
import torchvision.transforms as transforms
import os
from torch.autograd import Variable
import argparse
import numpy as np
from torch.optim.lr_scheduler import *
import csv

from model.resnet import resnet101
from data_pre.fashionTEST import fashiontest

parser = argparse.ArgumentParser()
parser.add_argument('--workers', type=int, default=2)
parser.add_argument('--batchSize', type=int, default=128)
parser.add_argument('--nepoch', type=int, default=1, help='number of epochs to train for')
parser.add_argument('--lr', type=float, default=0.0002, help='learning rate, default=0.0002')
parser.add_argument('--gpu', type=str, default='7', help='gpu ids: e.g. 0  0,1,2, 0,2. use -1 for CPU')
parser.add_argument('--attr', type=str,default='collar_design_labels')
opt = parser.parse_args()
print(opt)
os.environ["CUDA_VISIBLE_DEVICES"] = opt.gpu

transform_test = transforms.Compose([
    transforms.Scale((224, 224)),
    transforms.ToTensor(),
    transforms.Normalize((0.485, 0.456, 0.406), (0.229, 0.224, 0.225)),
])

AttrNum = {
    'coat_length_labels':8,
    'collar_design_labels':5,
    'lapel_design_labels':5,
    'neck_design_labels':5,
    'neckline_design_labels':10,
    'pant_length_labels':6,
    'skirt_length_labels':6,
    'sleeve_length_labels':9,
}

testset=fashiontest('/home/yhf/Challenge/FashionAI/ex_STL_FashionAI/data/2rank/Tests/question.csv',transform_test,opt.attr)
testloader=torch.utils.data.DataLoader(testset,batch_size=opt.batchSize,shuffle=False,num_workers=opt.workers)
model=resnet101(pretrained=True)
model.fc=nn.Linear(2048,AttrNum[opt.attr])
model.load_state_dict(torch.load('ckp/model_task_%s.pth'%opt.attr))
model.cuda()
model.eval()
results=[]

for image,addrs in testloader:
	image=Variable(image.cuda(),volatile=True)
	out=model(image)
	out=np.exp(out.cpu().data.numpy()).tolist()
	results.extend([[j,opt.attr,";".join([str(ii) for ii in i])] for (i,j) in zip(out,addrs)])

eval_csv=os.path.join(os.path.expanduser('.'),'deploy',opt.attr+'_eval.csv')
with open(eval_csv,'w',newline='') as f:
	writer=csv.writer(f,delimiter=',')
	for x in results:
		writer.writerow(x)

max_stars_repo_name: JackwithWilshere/FashionAI-

max_issues_repo_name: espectre/FashionAI

max_forks_repo_name: espectre/FashionAI


text:
# -*- coding: utf-8 -*-
# Authors: Victor Férat <victor.ferat@live.fr>
#
# License: BSD (3-clause)
import os
import pytest
import numpy as np

from ..io.inverse_solution import read_is
from ..source_estimate import read_ris, SourceEstimate
from ..regions_of_interest import RegionsOfInterest
from ..source_space import SourceSpace

dir_path = os.path.dirname(os.path.realpath(__file__))
data_path = os.path.join(dir_path, 'data')

def generate_source_space(size):
    """Generate Source space."""
    source_names = [('S' + str(i)) for i in range(0, size)]
    source_coords = np.random.rand(size, 3)
    source_space = SourceSpace(source_names, source_coords)
    return(source_space)

def generate_source_estimate(n_sources, n_times, sfreq):
    """Generate Source Estimate."""
    sources_tc = np.random.rand(n_sources, 3, n_times)
    source_estimate = SourceEstimate(sources_tc, sfreq)
    return(source_estimate)

def generate_rois(n_roi, n_sources):
    """Generate Rois."""
    rois_names = [('R' + str(i)) for i in range(0, n_roi)]
    roi_indices = []
    for i in range(0, n_roi):
        indices = np.random.randint(0, n_sources, size=i+1).tolist()
        roi_indices.append(indices)
    regions_of_interest = RegionsOfInterest(rois_names, roi_indices)
    return(regions_of_interest)

@pytest.mark.skip(reason='Requires large file')
def test_read_is():
    """Test read_is."""
    file_path = os.path.join('test')
    inverse_solution = read_is(file_path)
    if not inverse_solution['is_type'] == 'IS03':
        raise AssertionError()

def test_read_ris():
    """Test read_ris."""
    file_path = os.path.join(data_path, 'sample_test_ris.ris')
    read_ris(file_path)

def test_write_ris():
    """Test write_is."""
    file_path = os.path.join(data_path, 'sample_test_write_ris.ris')
    sources_tc = np.random.rand(5000, 3, 2048).astype('float32')
    sfreq = 512
    subject = 'test_subject'
    source_estimate = SourceEstimate(sources_tc, sfreq, subject=subject)
    source_estimate.save(file_path)
    read_source_estimate = read_ris(file_path)
    if not (read_source_estimate.sources_tc == sources_tc).all():
        raise AssertionError()
    if not read_source_estimate.sfreq == sfreq:
        raise AssertionError()
    if not read_source_estimate.filename == file_path:
        raise AssertionError()

def test_compute_tc():
    """Test SourceEstimate.compute_tc."""
    source_space = generate_source_space(5000)
    source_estimate = generate_source_estimate(5000, 2048, 512)
    source_estimate.source_space = source_space
    tc_mean = source_estimate.compute_tc(method='mean')
    tc_med = source_estimate.compute_tc(method='median')
    tc_svd = source_estimate.compute_tc(method='svd')
    print(tc_mean.shape)
    if not tc_mean.shape == (3, 2048):
        raise AssertionError()
    if not tc_med.shape == (3, 2048):
        raise AssertionError()
    if not tc_svd.shape == (1, 2048):
        raise AssertionError()

def test_compute_rois_tc():
    """Test SourceEstimate.compute_rois_tc."""
    source_space = generate_source_space(5000)
    source_estimate = generate_source_estimate(5000, 2048, 512)
    source_estimate.source_space = source_space
    regions_of_interest = generate_rois(32, 5000)
    regions_of_interest.source_space = source_space
    rois_estimates = source_estimate.compute_rois_tc(regions_of_interest)
    if not rois_estimates.sources_tc.shape == (32, 1, 2048):
        raise AssertionError()
    if not rois_estimates.sfreq == 512:
        raise AssertionError()

def test_per_roi():
    """Test SourceEstimate.per_roi"""
    source_space = generate_source_space(5000)
    source_estimate = generate_source_estimate(5000, 2048, 512)
    source_estimate.source_space = source_space
    regions_of_interest = generate_rois(32, 5000)
    regions_of_interest.source_space = source_space
    rois_estimates = source_estimate.per_roi(regions_of_interest)
    if not len(rois_estimates) == 32:
        raise AssertionError()

max_stars_repo_name: apadee/PyCartool

max_issues_repo_name: apadee/PyCartool

max_forks_repo_name: apadee/PyCartool


text:
import numpy as np
from mayavi import mlab

class ConnectedPoints:
    def __init__(
            self, color, line_color, mode, scale_factor, opacity_points, opacity_lines,
            point_names, high_quality=False):
        self.coords = []
        self.n_rendered_points = 0
        self.color = color
        self.mode = mode
        self.scale_factor = scale_factor
        self.line_color = line_color
        self.opacity_points = opacity_points
        self.opacity_lines = opacity_lines
        self.is_initialized = False
        self.points = None
        self.tube = None
        self.edges = []
        self.point_name_components = []
        self.point_names = point_names
        self.high_quality = high_quality

    def clear(self):
        self.coords.clear()
        self.edges.clear()

    def add_points(self, coords, edges, ignore_isolated=True):
        if ignore_isolated:
            coords, edges = self.get_nonisolated(coords, edges)
        n_points = len(self.coords)
        self.coords.extend(coords)
        self.edges.extend([[n_points + i, n_points + j] for i, j in edges])

    def get_nonisolated(self, coords, edges):
        selected_coords = []
        new_point_ids = {}
        selected_edges = []

        for (i, j) in edges:
            if i not in new_point_ids:
                new_point_ids[i] = len(selected_coords)
                selected_coords.append(coords[i])
            if j not in new_point_ids:
                new_point_ids[j] = len(selected_coords)
                selected_coords.append(coords[j])
            selected_edges.append([new_point_ids[i], new_point_ids[j]])
        return selected_coords, selected_edges

    def update(self):
        if not self.is_initialized:
            return self.initial_update()
        c = np.asarray(self.coords)

        if not self.coords:
            c = c.reshape((0, 3))

        if len(self.coords) == self.n_rendered_points:
            self.points.mlab_source.set(x=c[:, 0], y=c[:, 1], z=c[:, 2])
        else:
            self.n_rendered_points = len(self.coords)
            self.points.mlab_source.reset(x=c[:, 0], y=c[:, 1], z=c[:, 2])
            self.points.mlab_source.dataset.lines = self.edges

    def initial_update(self):
        if not self.coords:
            return
        c = np.asarray(self.coords)

        self.points = mlab.points3d(
            *c.T, scale_factor=self.scale_factor, color=self.color,
            opacity=self.opacity_points,
            mode=self.mode, resolution=8 if self.high_quality else 4, scale_mode='vector',
            reset_zoom=False)
        self.points.mlab_source.dataset.lines = self.edges
        tube = mlab.pipeline.tube(
            self.points, tube_radius=self.scale_factor / 3,
            tube_sides=6 if self.high_quality else 3)
        mlab.pipeline.surface(
            tube, color=self.line_color, opacity=self.opacity_lines, reset_zoom=False)

        # for i, (point, name, color) in enumerate(
        #         zip(c, self.point_names, colors.cycle_over_colors())):
        #     self.point_name_components.append(mlab.text3d(
        #         *point, f'{i}', scale=0.02, color=color))

        self.is_initialized = True

max_stars_repo_name: HY-Vision-Lab/metrabs

max_issues_repo_name: HY-Vision-Lab/metrabs

max_forks_repo_name: HY-Vision-Lab/metrabs


text:
from __future__ import print_function, division, absolute_import

import re
import time
import numpy as np

re_pos = re.compile('p=<(.+?)>')
re_vel = re.compile('v=<(.+?)>')
re_acc = re.compile('a=<(.+?)>')

np.set_printoptions(linewidth=1024, edgeitems=9)

def solve(inp, part=2):

    pmat = np.zeros([len(inp), 9], dtype=int)

    for id, line in enumerate(inp):
        pmatch = re_pos.search(line)
        vmatch = re_vel.search(line)
        amatch = re_acc.search(line)

        pos = np.array([int(s) for s in pmatch.groups()[0].split(',')], dtype=int)
        vel = np.array([int(s) for s in vmatch.groups()[0].split(',')], dtype=int)
        acc = np.array([int(s) for s in amatch.groups()[0].split(',')], dtype=int)

        pmat[id,:3] = pos
        pmat[id,3:6] = vel
        pmat[id,6:9] = acc

    # Simulate
    collided = set()
    for i in range(500):
        pmat[:, 3:6] += pmat[:, 6:9]  # update velocity
        pmat[:, :3] += pmat[:, 3:6]   # update position

        if part == 2:
            for row in range(pmat.shape[0]):
                if row in collided:
                    continue
                colliding_rows = np.where((pmat[:, :3] == pmat[row, :3]).all(axis=1))[0]
                if len(colliding_rows) > 1:
                    # Kill the collided rows
                    pmat[colliding_rows, :3] = 1000000000
                    pmat[colliding_rows, 3:] = 0
                    collided |= set(colliding_rows.tolist())

    # manhattan distance
    result = np.argmin(np.abs(pmat[:, :3]).sum(axis=1))

    if part == 1:
        print('min manhattan distance is particle #:', result)

    if part == 2:
        print('surviving particles:', pmat.shape[0] - len(collided))

if __name__ == '__main__':

    with open('input.txt', 'r') as f:
        puzzle_input = [line.strip() for line in f.readlines() if line]

    t0 = time.time()
    solve(puzzle_input, part=1)
    print('Time to solve part 1:', time.time()-t0, 'sec')

    t0 = time.time()
    solve(puzzle_input, part=2)
    print('Time to solve part 2:', time.time()-t0, 'sec')


text:
import numpy
import tensorflow

from google.cloud import storage
from tensorflow.keras.layers import Dense, Flatten, Conv2D
from tensorflow.keras import Model
from PIL import Image

# We keep model as global variable so we don't have to reload it in case of warm invocations
model = None

class CustomModel(Model):
    def __init__(self):
        super(CustomModel, self).__init__()
        self.conv1 = Conv2D(32, 3, activation="relu")
        self.flatten = Flatten()
        self.d1 = Dense(128, activation="relu")
        self.d2 = Dense(10, activation="softmax")

    def call(self, x):
        x = self.conv1(x)
        x = self.flatten(x)
        x = self.d1(x)
        return self.d2(x)

def download_blob(bucket_name, source_blob_name, destination_file_name):
    """Downloads a blob from the bucket."""
    storage_client = storage.Client()
    bucket = storage_client.get_bucket(bucket_name)
    blob = bucket.blob(source_blob_name)

    blob.download_to_filename(destination_file_name)

    print("Blob {} downloaded to {}.".format(source_blob_name, destination_file_name))

def handler(request):
    global model
    class_names = [
        "T-shirt/top",
        "Trouser",
        "Pullover",
        "Dress",
        "Coat",
        "Sandal",
        "Shirt",
        "Sneaker",
        "Bag",
        "Ankle boot",
    ]

    # Model load which only happens during cold starts
    if model is None:
        download_blob(
            "tfserving1",
            "tensorflow/fashion_mnist_weights.index",
            "/tmp/fashion_mnist_weights.index",
        )
        download_blob(
            "tfserving1",
            "tensorflow/fashion_mnist_weights.data-00000-of-00001",
            "/tmp/fashion_mnist_weights.data-00000-of-00001",
        )
        model = CustomModel()
        model.load_weights("/tmp/fashion_mnist_weights")

    download_blob("tfserving1", "tensorflow/test.png", "/tmp/test.png")
    image = numpy.array(Image.open("/tmp/test.png"))
    input_np = (numpy.array(Image.open("/tmp/test.png")) / 255)[
        numpy.newaxis, :, :, numpy.newaxis
    ]
    predictions = model.call(input_np)
    print(predictions)
    print("Image is " + class_names[numpy.argmax(predictions)])

    return class_names[numpy.argmax(predictions)]

max_stars_repo_name: r-matsuzaka/cloud-function-example

max_issues_repo_name: r-matsuzaka/cloud-function-example

max_forks_repo_name: r-matsuzaka/cloud-function-example


text:
import sys, os, random, scipy
import numpy as np
import pandas as pd
from numba.typed import List
import statsmodels.api as sm
from statsmodels.stats.multitest import multipletests

from tqdm import tqdm
from sklearn.cluster import AgglomerativeClustering

from anndata import AnnData
from .base import lr, calc_neighbours, get_spot_lrs, get_lrs_scores, get_scores
from .merge import merge
from .perm_utils import get_lr_features, get_lr_bg

# Newest method #
def perform_spot_testing(
    adata: AnnData,
    lr_scores: np.ndarray,
    lrs: np.array,
    n_pairs: int,
    neighbours: List,
    het_vals: np.array,
    min_expr: float,
    adj_method: str = "fdr_bh",
    pval_adj_cutoff: float = 0.05,
    verbose: bool = True,
    save_bg=False,
    neg_binom=False,
    quantiles=(0.5, 0.75, 0.85, 0.9, 0.95, 0.97, 0.98, 0.99, 0.995, 0.9975, 0.999, 1),
):
    """Calls significant spots by creating random gene pairs with similar
    expression to given LR pair; only generate background for spots
    which have score for given LR.
    """
    quantiles = np.array(quantiles)

    lr_genes = np.unique([lr_.split("_") for lr_ in lrs])
    genes = np.array([gene for gene in adata.var_names if gene not in lr_genes])
    candidate_expr = adata[:, genes].to_df().values

    n_genes = round(np.sqrt(n_pairs) * 2)
    if len(genes) < n_genes:
        print(
            "Exiting since need atleast "
            f"{n_genes} genes to generate {n_pairs} pairs."
        )
        return

    if n_pairs < 100:
        print(
            "Exiting since n_pairs<100, need much larger number of pairs to "
            "get accurate backgrounds (e.g. 1000)."
        )
        return

    ####### Quantiles to select similar gene to LRs to gen. rand-pairs #######
    lr_expr = adata[:, lr_genes].to_df()
    lr_feats = get_lr_features(adata, lr_expr, lrs, quantiles)
    l_quants = lr_feats.loc[
        lrs, [col for col in lr_feats.columns if "L_" in col]
    ].values
    r_quants = lr_feats.loc[
        lrs, [col for col in lr_feats.columns if "R_" in col]
    ].values
    candidate_quants = np.apply_along_axis(
        np.quantile, 0, candidate_expr, q=quantiles, interpolation="nearest"
    )
    # Ensuring consistent typing to prevent numba errors #
    l_quants = l_quants.astype("<f4")
    r_quants = r_quants.astype("<f4")
    candidate_quants = candidate_quants.astype("<f4")

    ######## Background per LR, but only for spots where LR has a score ########
    # Determine the indices of the spots where each LR has a score #
    cols = ["n_spots", "n_spots_sig", "n_spots_sig_pval"]
    lr_summary = np.zeros((lr_scores.shape[1], 3), np.int)
    pvals = np.ones(lr_scores.shape, dtype=np.float64)
    pvals_adj = np.ones(lr_scores.shape, dtype=np.float64)
    log10pvals_adj = np.zeros(lr_scores.shape, dtype=np.float64)
    lr_sig_scores = lr_scores.copy()

    # If we are saving the backgrounds #
    if save_bg:
        adata.uns["lrs_to_bg"] = {}
        adata.uns["lr_spot_indices"] = {}

    with tqdm(
        total=lr_scores.shape[1],
        desc="Generating backgrounds & testing each LR pair...",
        bar_format="{l_bar}{bar} [ time left: {remaining} ]",
        disable=verbose == False,
    ) as pbar:

        gene_bg_genes = {}  # Keep track of genes which can be used to gen. rand-pairs.
        spot_lr_indices = [
            [] for i in range(lr_scores.shape[0])
        ]  # tracks the lrs tested in a given spot for MHT !!!!
        for lr_j in range(lr_scores.shape[1]):
            lr_ = lrs[lr_j]

            background, spot_indices = get_lr_bg(
                adata,
                neighbours,
                het_vals,
                min_expr,
                lr_,
                lr_scores[:, lr_j],
                l_quants[lr_j, :],
                r_quants[lr_j, :],
                genes,
                candidate_quants,
                gene_bg_genes,
                n_genes,
                n_pairs,
            )

            if save_bg:
                adata.uns["lrs_to_bg"][lr_] = background
                adata.uns["lr_spot_indices"][lr_] = spot_indices

            if not neg_binom:  # Calculate empirical p-values per-spot
                for spot_i, spot_index in enumerate(spot_indices):
                    n_greater = len(
                        np.where(background[spot_i, :] >= lr_scores[spot_index, lr_j])[
                            0
                        ]
                    )
                    n_greater = n_greater if n_greater != 0 else 1  # pseudocount
                    pvals[spot_index, lr_j] = n_greater / background.shape[1]
                    spot_lr_indices[spot_index].append(lr_j)
            else:  # Fitting NB per LR
                lr_j_scores = lr_scores[spot_indices, lr_j]
                bg_ = background.ravel()
                bg_wScore = np.array(list(lr_j_scores) + list(bg_))

                ##### 1) rounding discretisation
                # First multiple to get minimum value to be one before rounding #
                bg_1 = bg_wScore * (1 / min(bg_wScore[bg_wScore != 0]))
                bg_1 = np.round(bg_1)
                lr_j_scores_1 = bg_1[0 : len(lr_j_scores)]
                bg_1 = bg_1[len(lr_j_scores) : len(bg_1)]

                ###### Getting the pvalue from negative binomial approach
                round_pvals, _, _, _ = get_stats(
                    lr_j_scores_1,
                    bg_1,
                    len(bg_1),
                    neg_binom=True,
                    return_negbinom_params=False,
                )
                pvals[spot_indices, lr_j] = round_pvals
                for spot_index in spot_indices:
                    spot_lr_indices[spot_index].append(lr_j)

            pbar.update(1)

        # MHT correction # filling in other stats #
        for spot_i in range(lr_scores.shape[0]):
            lr_indices = spot_lr_indices[spot_i]
            if len(lr_indices) != 0:
                pvals_adj[spot_i, lr_indices] = multipletests(
                    pvals[spot_i, lr_indices], method=adj_method
                )[1]

            log10pvals_adj[spot_i, :] = -np.log10(pvals_adj[spot_i, :])

            # Recording per lr results for this LR #
            lrs_in_spot = lr_scores[spot_i] > min_expr
            sig_lrs_in_spot = pvals_adj[spot_i, :] < pval_adj_cutoff
            sigpval_lrs_in_spot = pvals[spot_i, :] < pval_adj_cutoff
            lr_summary[lrs_in_spot, 0] += 1
            lr_summary[sig_lrs_in_spot, 1] += 1
            lr_summary[sigpval_lrs_in_spot, 2] += 1

            lr_sig_scores[spot_i, sig_lrs_in_spot == False] = 0

    # Ordering the results according to number of significant spots per LR#
    order = np.argsort(-lr_summary[:, 1])
    lrs_ordered = lrs[order]
    lr_summary = lr_summary[order, :]
    lr_summary = pd.DataFrame(lr_summary, index=lrs_ordered, columns=cols)
    lr_scores = lr_scores[:, order]
    pvals = pvals[:, order]
    pvals_adj = pvals_adj[:, order]
    log10pvals_adj = log10pvals_adj[:, order]
    lr_sig_scores = lr_sig_scores[:, order]

    # Saving the results in AnnData #
    if verbose:
        print("\nStoring results:\n")

    adata.uns["lr_summary"] = lr_summary
    res_info = ["lr_scores", "p_vals", "p_adjs", "-log10(p_adjs)", "lr_sig_scores"]
    mats = [lr_scores, pvals, pvals_adj, log10pvals_adj, lr_sig_scores]
    for i, info_name in enumerate(res_info):
        adata.obsm[info_name] = mats[i]
        if verbose:
            print(f"{info_name} stored in adata.obsm['{info_name}'].")

    if verbose:
        print(
            "\nPer-spot results in adata.obsm have columns in same order as "
            "rows in adata.uns['lr_summary']."
        )
        print("Summary of LR results in adata.uns['lr_summary'].")

# Version 2, no longer in use, see above for newest method #
def perform_perm_testing(
    adata: AnnData,
    lr_scores: np.ndarray,
    n_pairs: int,
    lrs: np.array,
    lr_mid_dist: int,
    verbose: float,
    neighbours: List,
    het_vals: np.array,
    min_expr: float,
    neg_binom: bool,
    adj_method: str,
    pval_adj_cutoff: float,
):
    """Performs the grouped permutation testing when taking the stats approach."""
    if n_pairs != 0:  # Perform permutation testing
        # Grouping spots with similar mean expression point #
        genes = get_valid_genes(adata, n_pairs)
        means_ordered, genes_ordered = get_ordered(adata, genes)
        ims = np.array(
            [
                get_median_index(
                    lr_.split("_")[0],
                    lr_.split("_")[1],
                    means_ordered.values,
                    genes_ordered,
                )
                for lr_ in lrs
            ]
        ).reshape(-1, 1)

        if len(lrs) > 1:  # Multi-LR pair mode, group LRs to generate backgrounds
            clusterer = AgglomerativeClustering(
                n_clusters=None,
                distance_threshold=lr_mid_dist,
                affinity="manhattan",
                linkage="single",
            )
            lr_groups = clusterer.fit_predict(ims)
            lr_group_set = np.unique(lr_groups)
            if verbose:
                print(f"{len(lr_group_set)} lr groups with similar expression levels.")

        else:  # Single LR pair mode, generate background for the LR.
            lr_groups = np.array([0])
            lr_group_set = lr_groups

        res_info = ["lr_scores", "p_val", "p_adj", "-log10(p_adj)", "lr_sig_scores"]
        n_, n_sigs = np.array([0] * len(lrs)), np.array([0] * len(lrs))
        per_lr_results = {}
        with tqdm(
            total=len(lr_group_set),
            desc="Generating background distributions for the LR pair groups..",
            bar_format="{l_bar}{bar} [ time left: {remaining} ]",
        ) as pbar:
            for group in lr_group_set:
                # Determining common mid-point for each group #
                group_bool = lr_groups == group
                group_im = int(np.median(ims[group_bool, 0]))

                # Calculating the background #
                rand_pairs = get_rand_pairs(adata, genes, n_pairs, lrs=lrs, im=group_im)
                background = get_lrs_scores(
                    adata,
                    rand_pairs,
                    neighbours,
                    het_vals,
                    min_expr,
                    filter_pairs=False,
                ).ravel()
                total_bg = len(background)
                background = background[background != 0]  # Filtering for increase speed

                # Getting stats for each lr in group #
                group_lr_indices = np.where(group_bool)[0]
                for lr_i in group_lr_indices:
                    lr_ = lrs[lr_i]
                    lr_results = pd.DataFrame(index=adata.obs_names, columns=res_info)
                    scores = lr_scores[:, lr_i]
                    stats = get_stats(
                        scores,
                        background,
                        total_bg,
                        neg_binom,
                        adj_method,
                        pval_adj_cutoff=pval_adj_cutoff,
                    )
                    full_stats = [scores] + list(stats)
                    for vals, colname in zip(full_stats, res_info):
                        lr_results[colname] = vals

                    n_[lr_i] = len(np.where(scores > 0)[0])
                    n_sigs[lr_i] = len(
                        np.where(lr_results["p_adj"].values < pval_adj_cutoff)[0]
                    )
                    if n_sigs[lr_i] > 0:
                        per_lr_results[lr_] = lr_results
                pbar.update(1)

        print(f"{len(per_lr_results)} LR pairs with significant interactions.")

        lr_summary = pd.DataFrame(index=lrs, columns=["n_spots", "n_spots_sig"])
        lr_summary["n_spots"] = n_
        lr_summary["n_spots_sig"] = n_sigs
        lr_summary = lr_summary.iloc[np.argsort(-n_sigs)]

    else:  # Simply store the scores
        per_lr_results = {}
        lr_summary = pd.DataFrame(index=lrs, columns=["n_spots"])
        for i, lr_ in enumerate(lrs):
            lr_results = pd.DataFrame(index=adata.obs_names, columns=["lr_scores"])
            lr_results["lr_scores"] = lr_scores[:, i]
            per_lr_results[lr_] = lr_results
            lr_summary.loc[lr_, "n_spots"] = len(np.where(lr_scores[:, i] > 0)[0])
        lr_summary = lr_summary.iloc[np.argsort(-lr_summary.values[:, 0]), :]

    adata.uns["lr_summary"] = lr_summary
    adata.uns["per_lr_results"] = per_lr_results
    if verbose:
        print(
            "Summary of significant spots for each lr pair in adata.uns['lr_summary']."
        )
        print(
            "Spot enrichment statistics of LR interactions in adata.uns['per_lr_results']"
        )

# No longer in use #
def permutation(
    adata: AnnData,
    n_pairs: int = 200,
    distance: int = None,
    use_lr: str = "cci_lr",
    use_het: str = None,
    neg_binom: bool = False,
    adj_method: str = "fdr",
    neighbours: list = None,
    run_fast: bool = True,
    bg_pairs: list = None,
    background: np.array = None,
    **kwargs,
) -> AnnData:

    """Permutation test for merged result
    Parameters
    ----------
    adata: AnnData          The data object including the cell types to count
    n_pairs: int            Number of gene pairs to run permutation test (default: 1000)
    distance: int           Distance between spots (default: 30)
    use_lr: str             LR cluster used for permutation test (default: 'lr_neighbours_louvain_max')
    use_het: str            cell type diversity counts used for permutation test (default 'het')
    neg_binom: bool         Whether to fit neg binomial paramaters to bg distribution for p-val est.
    adj_method: str         Method used by statsmodels.stats.multitest.multipletests for MHT correction.
    neighbours: list        List of the neighbours for each spot, if None then computed. Useful for speeding up function.
    **kwargs:               Extra arguments parsed to lr.
    Returns
    -------
    adata: AnnData          Data Frame of p-values from permutation test for each window stored in adata.uns['merged_pvalues']
                            Final significant merged scores stored in adata.uns['merged_sign']
    """

    # blockPrint()

    #  select n_pair*2 closely expressed genes from the data
    genes = get_valid_genes(adata, n_pairs)
    if len(adata.uns["lr"]) > 1:
        raise ValueError("Permutation test only supported for one LR pair scenario.")
    elif type(bg_pairs) == type(None):
        pairs = get_rand_pairs(adata, genes, n_pairs, lrs=adata.uns["lr"])
    else:
        pairs = bg_pairs

    """
    # generate random pairs
    lr1 = adata.uns['lr'][0].split('_')[0]
    lr2 = adata.uns['lr'][0].split('_')[1]
    genes = [item for item in adata.var_names.tolist() if not (item.startswith('MT-') or item.startswith('MT_') or item==lr1 or item==lr2)]
    random.shuffle(genes)
    pairs = [i + '_' + j for i, j in zip(genes[:n_pairs], genes[-n_pairs:])]
    """
    if use_het != None:
        scores = adata.obsm["merged"]
    else:
        scores = adata.obsm[use_lr]

    # for each randomly selected pair, run through cci_rank analysis and keep the scores
    query_pair = adata.uns["lr"]

    # If neighbours not inputted, then compute #
    if type(neighbours) == type(None):
        neighbours = calc_neighbours(adata, distance, index=run_fast)

    if not run_fast and type(background) == type(
        None
    ):  # Run original way if 'fast'=False argument inputted.
        background = []
        for item in pairs:
            adata.uns["lr"] = [item]
            lr(
                adata,
                use_lr=use_lr,
                distance=distance,
                verbose=False,
                neighbours=neighbours,
                **kwargs,
            )
            if use_het != None:
                merge(adata, use_lr=use_lr, use_het=use_het, verbose=False)
                background += adata.obsm["merged"].tolist()
            else:
                background += adata.obsm[use_lr].tolist()
        background = np.array(background)

    elif type(background) == type(None):  # Run fast if background not inputted
        spot_lr1s = get_spot_lrs(adata, pairs, lr_order=True, filter_pairs=False)
        spot_lr2s = get_spot_lrs(adata, pairs, lr_order=False, filter_pairs=False)

        het_vals = (
            np.array([1] * len(adata)) if use_het == None else adata.obsm[use_het]
        )
        background = get_scores(
            spot_lr1s.values, spot_lr2s.values, neighbours, het_vals
        ).ravel()

    # log back the original query
    adata.uns["lr"] = query_pair

    #### Negative Binomial fit
    pvals, pvals_adj, log10_pvals, lr_sign = get_stats(
        scores, background, neg_binom, adj_method
    )

    if use_het != None:
        adata.obsm["merged"] = scores
        adata.obsm["merged_pvalues"] = log10_pvals
        adata.obsm["merged_sign"] = lr_sign

        # enablePrint()
        print(
            "Results of permutation test has been kept in adata.obsm['merged_pvalues']"
        )
        print("Significant merged result has been kept in adata.obsm['merged_sign']")
    else:
        adata.obsm[use_lr] = scores
        adata.obsm["lr_pvalues"] = log10_pvals
        adata.obsm["lr_sign"] = lr_sign  # scores for spots with pval_adj < 0.05

        # enablePrint()
        print("Results of permutation test has been kept in adata.obsm['lr_pvalues']")
        print("Significant merged result has been kept in adata.obsm['lr_sign']")

    # return adata
    return background

def get_stats(
    scores: np.array,
    background: np.array,
    total_bg: int,
    neg_binom: bool = False,
    adj_method: str = "fdr_bh",
    pval_adj_cutoff: float = 0.01,
    return_negbinom_params: bool = False,
):
    """Retrieves valid candidate genes to be used for random gene pairs.
    Parameters
    ----------
    scores: np.array        Per spot scores for a particular LR pair.
    background: np.array    Background distribution for non-zero scores.
    total_bg: int           Total number of background values calculated.
    neg_binom: bool         Whether to use neg-binomial distribution to estimate p-values, NOT appropriate with log1p data, alternative is to use background distribution itself (recommend higher number of n_pairs for this).
    adj_method: str         Parsed to statsmodels.stats.multitest.multipletests for multiple hypothesis testing correction.
    Returns
    -------
    stats: tuple          Per spot pvalues, pvals_adj, log10_pvals_adj, lr_sign (the LR scores for significant spots).
    """
    ##### Negative Binomial fit
    if neg_binom:
        # Need to make full background for fitting !!!
        background = np.array(list(background) + [0] * (total_bg - len(background)))
        pmin, pmax = min(background), max(background)
        background2 = [item - pmin for item in background]
        x = np.linspace(pmin, pmax, 1000)
        res = sm.NegativeBinomial(
            background2, np.ones(len(background2)), loglike_method="nb2"
        ).fit(start_params=[0.1, 0.3], disp=0)
        mu = res.predict()  # use if not constant
        mu = np.exp(res.params[0])
        alpha = res.params[1]
        Q = 0
        size = 1.0 / alpha * mu**Q
        prob = size / (size + mu)

        if return_negbinom_params:  # For testing purposes #
            return size, prob

        # Calculate probability for all spots
        pvals = 1 - scipy.stats.nbinom.cdf(scores - pmin, size, prob)

    else:  ###### Using the actual values to estimate p-values
        pvals = np.zeros((1, len(scores)), dtype=np.float)[0, :]
        nonzero_score_bool = scores > 0
        nonzero_score_indices = np.where(nonzero_score_bool)[0]
        zero_score_indices = np.where(nonzero_score_bool == False)[0]
        pvals[zero_score_indices] = (total_bg - len(background)) / total_bg
        pvals[nonzero_score_indices] = [
            len(np.where(background >= scores[i])[0]) / total_bg
            for i in nonzero_score_indices
        ]

    pvals_adj = multipletests(pvals, method=adj_method)[1]
    log10_pvals_adj = -np.log10(pvals_adj)
    lr_sign = scores * (pvals_adj < pval_adj_cutoff)
    return pvals, pvals_adj, log10_pvals_adj, lr_sign

# @njit(parallel=True)
# def perm_pvals(scores, background):
#     """Determines the p-values based on the actual observed frequency of the \
#         indicated score or greater in the background.
#     """
#     pvals = np.zeros((1, len(scores)), np.float64)[0,:]
#     for i in prange(len(pvals)):

def get_valid_genes(adata: AnnData, n_pairs: int) -> np.array:
    """Retrieves valid candidate genes to be used for random gene pairs.
    Parameters
    ----------
    adata: AnnData          The data object including the cell types to count
    n_pairs: int            The number of random pairs will generate elsewhere.
    Returns
    -------
    genes: np.array          List of genes which could be valid pairs.
    """
    genes = np.array(
        [
            item
            for item in adata.var_names.tolist()
            if not (item.startswith("MT-") or item.startswith("MT_"))
        ]
    )
    if n_pairs >= len(genes) / 2:
        raise ValueError(
            "Too many genes pairs selected, please reduce to a smaller number."
        )
    return genes

def get_rand_pairs(
    adata: AnnData,
    genes: np.array,
    n_pairs: int,
    lrs: list = None,
    im: int = None,
):
    """Gets equivalent random gene pairs for the inputted lr pair.
    Parameters
    ----------
    adata: AnnData          The data object including the cell types to count
    lr: int            The lr pair string to get equivalent random pairs for (e.g. 'L_R')
    genes: np.array           Candidate genes to use as pairs.
    n_pairs: int             Number of random pairs to generate.
    Returns
    -------
    pairs: list          List of random gene pairs with equivalent mean expression (e.g. ['L_R'])
    """
    lr_genes = [lr.split("_")[0] for lr in lrs]
    lr_genes += [lr.split("_")[1] for lr in lrs]

    # get the position of the median of the means between the two genes
    means_ordered, genes_ordered = get_ordered(adata, genes)
    if type(im) == type(None):  # Single background per lr pair mode
        l, r = lrs[0].split("_")
        im = get_median_index(l, r, means_ordered.values, genes_ordered)

    # get n_pair genes sorted by distance to im
    selected = (
        abs(means_ordered - means_ordered[im])
        .sort_values()
        .drop(lr_genes)[: n_pairs * 2]
        .index.tolist()
    )
    selected = selected[0 : n_pairs * 2]
    adata.uns["selected"] = selected
    # form gene pairs from selected randomly
    random.shuffle(selected)
    pairs = [i + "_" + j for i, j in zip(selected[:n_pairs], selected[-n_pairs:])]

    return pairs

def get_ordered(adata, genes):
    means_ordered = adata.to_df()[genes].mean().sort_values()
    genes_ordered = means_ordered.index.values
    return means_ordered, genes_ordered

def get_median_index(l, r, means_ordered, genes_ordered):
    """ "Retrieves the index of the gene with a mean expression between the two genes in the lr pair.
    Parameters
    ----------
    X: np.ndarray          Spot*Gene expression.
    l: str                 Ligand gene.
    r: str                 Receptor gene.
    genes: np.array        Candidate genes to use as pairs.
    Returns
    -------
    pairs: list          List of random gene pairs with equivalent mean expression (e.g. ['L_R'])
    """
    # sort the mean of each gene expression
    i1 = np.where(genes_ordered == l)[0][0]
    i2 = np.where(genes_ordered == r)[0][0]
    if means_ordered[i1] > means_ordered[i2]:
        it = i1
        i1 = i2
        i2 = it

    im = np.argmin(np.abs(means_ordered - np.median(means_ordered[i1:i2])))
    return im
    # means = adata.to_df()[genes].mean().sort_values()
    # lr1 = lr[0].split("_")[0]
    # lr2 = lr[0].split("_")[1]
    # i1, i2 = means.index.get_loc(lr1), means.index.get_loc(lr2)
    # if means[lr1] > means[lr2]:
    #     it = i1
    #     i1 = i2
    #     i2 = it
    #
    # # get the position of the median of the means between the two genes
    # im = np.argmin(abs(means.values - means.iloc[i1:i2]))

# Disable printing
def blockPrint():
    sys.stdout = open(os.devnull, "w")

# Restore printing
def enablePrint():
    sys.stdout = sys.__stdout__

reason: import numpy,from numba,import statsmodels,from statsmodels