text:
import numpy as np

def bond_hash_dualdirect(bond, natoms):
    """
    :param bond: bond data in hoomdxml format (name, id1, id2)
    :param natoms: total number of particles
    :return: hash table of with value in {bondname1: [idxes], bondname2:[idxes]...} for each particle (in dual direct)
    """
    bond_hash_nn = {}
    print('Building bond hash...')
    if not isinstance(bond, np.ndarray):
        return {}
    for i in range(natoms):
        bond_hash_nn[i] = []
    for b in bond:
        idx = b[1]
        jdx = b[2]
        bond_hash_nn[idx].append(jdx)
        bond_hash_nn[jdx].append(idx)
    # print('Done.')
    return bond_hash_nn

max_stars_repo_name: atomrq/simulab

max_issues_repo_name: atomrq/simulab

max_forks_repo_name: atomrq/simulab


text:
# -*- coding: utf-8 -*-
# ProDy: A Python Package for Protein Dynamics Analysis
#
# Copyright (C) 2010-2012 Ahmet Bakan
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>

"""This module defines miscellaneous functions dealing with protein data."""

__author__ = 'Ahmet Bakan'
__copyright__ = 'Copyright (C) 2010-2012 Ahmet Bakan'

import numpy as np

from prody.atomic import Atomic, Atom, AtomGroup, Selection, HierView
from prody.utilities import openFile, showFigure
from prody import SETTINGS

__all__ = ['showProtein', 'writePQR', ]

def writePQR(filename, atoms):
    """Write *atoms* in PQR format to a file with name *filename*.  Only
    current coordinate set is written.  Returns *filename* upon success.  If
    *filename* ends with :file:`.gz`, a compressed file will be written."""

    if not isinstance(atoms, Atomic):
        raise TypeError('atoms does not have a valid type')
    if isinstance(atoms, Atom):
        atoms = Selection(atoms.getAtomGroup(), [atoms.getIndex()],
                          atoms.getACSIndex(),
                          'index ' + str(atoms.getIndex()))
    stream = openFile(filename, 'w')
    n_atoms = atoms.numAtoms()
    atomnames = atoms.getNames()
    if atomnames is None:
        raise RuntimeError('atom names are not set')
    for i, an in enumerate(atomnames):
        lenan = len(an)
        if lenan < 4:
            atomnames[i] = ' ' + an
        elif lenan > 4:
            atomnames[i] = an[:4]

    s_or_u = np.array(['a']).dtype.char

    resnames = atoms._getResnames()
    if resnames is None:
        resnames = ['UNK'] * n_atoms
    resnums = atoms._getResnums()
    if resnums is None:
        resnums = np.ones(n_atoms, int)
    chainids = atoms._getChids()
    if chainids is None:
        chainids = np.zeros(n_atoms, s_or_u + '1')
    charges = atoms._getCharges()
    if charges is None:
        charges = np.zeros(n_atoms, float)
    radii = atoms._getRadii()
    if radii is None:
        radii = np.zeros(n_atoms, float)
    icodes = atoms._getIcodes()
    if icodes is None:
        icodes = np.zeros(n_atoms, s_or_u + '1')
    hetero = ['ATOM'] * n_atoms
    heteroflags = atoms._getFlags('hetatm')
    if heteroflags is None:
        heteroflags = atoms._getFlags('hetero')
    if heteroflags is not None:
        hetero = np.array(hetero, s_or_u + '6')
        hetero[heteroflags] = 'HETATM'
    altlocs = atoms._getAltlocs()
    if altlocs is None:
        altlocs = np.zeros(n_atoms, s_or_u + '1')

    format = ('{0:6s}{1:5d} {2:4s}{3:1s}' +
              '{4:4s}{5:1s}{6:4d}{7:1s}   ' +
              '{8:8.3f}{9:8.3f}{10:8.3f}' +
              '{11:8.4f}{12:7.4f}\n').format
    coords = atoms._getCoords()
    write = stream.write
    for i, xyz in enumerate(coords):
        write(format(hetero[i], i+1, atomnames[i], altlocs[i],
                     resnames[i], chainids[i], int(resnums[i]),
                     icodes[i], xyz[0], xyz[1], xyz[2], charges[i], radii[i]))
    write('TER\nEND')
    stream.close()
    return filename

def showProtein(*atoms, **kwargs):
    """Show protein representation using :meth:`~mpl_toolkits.mplot3d.Axes3D`.
    This function is designed for generating a quick view of the contents of a
    :class:`~.AtomGroup` or :class:`~.Selection`.

    Protein atoms matching ``"calpha"`` selection are displayed using solid
    lines by picking a random and unique color per chain.  Line with can
    be adjusted using *lw* argument, e.g. ``lw=12``. Default width is 4.
    Chain colors can be overwritten using chain identifier as in ``A='green'``.

    Water molecule oxygen atoms are represented by red colored circles.  Color
    can be changed using *water* keyword argument, e.g. ``water='aqua'``.
    Water marker and size can be changed using *wmarker* and *wsize* keywords,
    defaults values are ``wmarker='.', wsize=6``.

    Hetero atoms matching ``"hetero and noh"`` selection are represented by
    circles and unique colors are picked at random on a per residue basis.
    Colors can be customized using residue name as in ``NAH='purple'``.  Note
    that this will color all distinct residues with the same name in the same
    color.  Hetero atom marker and size can be changed using *hmarker* and
    *hsize* keywords, default values are ``hmarker='o', hsize=6``.

    ProDy will set the size of axis so the representation is not distorted when
    the shape of figure window is close to a square.  Colors are picked at
    random, except for water oxygens which will always be colored red."""

    alist = atoms
    for atoms in alist:
        if not isinstance(atoms, Atomic):
            raise TypeError('atoms must be an Atomic instance')
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    cf = plt.gcf()
    show = None
    for child in cf.get_children():
        if isinstance(child, Axes3D):
            show = child
            break
    if show is None:
        show = Axes3D(cf)
    from matplotlib import colors
    cnames = dict(colors.cnames)
    wcolor = kwargs.get('water', 'red').lower()
    avoid = np.array(colors.hex2color(cnames.pop(wcolor, cnames.pop('red'))))
    for cn, val in cnames.items():  # PY3K: OK
        clr = np.array(colors.hex2color(val))
        if clr.sum() > 2.4:
            cnames.pop(cn)
        elif np.abs(avoid - clr).sum() <= 0.6:
            cnames.pop(cn)
    cnames = list(cnames)
    import random
    random.shuffle(cnames)
    min_ = list()
    max_ = list()
    for atoms in alist:
        if isinstance(atoms, AtomGroup):
            title = atoms.getTitle()
        else:
            title = atoms.getAtomGroup().getTitle()
        calpha = atoms.select('calpha')
        if calpha:
            for ch in HierView(calpha, chain=True):
                xyz = ch._getCoords()
                chid = ch.getChid()
                show.plot(xyz[:, 0], xyz[:, 1], xyz[:, 2],
                          label=title + '_' + chid,
                          color=kwargs.get(chid, cnames.pop()).lower(),
                          lw=kwargs.get('lw', 4))
        water = atoms.select('water and noh')
        if water:
            xyz = atoms.select('water')._getCoords()
            show.plot(xyz[:, 0], xyz[:, 1], xyz[:, 2], label=title + '_water',
                      color=wcolor,
                      ls='None', marker=kwargs.get('wmarker', '.'),
                      ms=kwargs.get('wsize', 6))
        hetero = atoms.select('not protein and not nucleic and not water')
        if hetero:
            for res in HierView(hetero).iterResidues():
                xyz = res._getCoords()
                resname = res.getResname()
                resnum = str(res.getResnum())
                chid = res.getChid()
                show.plot(xyz[:, 0], xyz[:, 1], xyz[:, 2], ls='None',
                          color=kwargs.get(resname, cnames.pop()).lower(),
                          label=title + '_' + chid + '_' + resname + resnum,
                          marker=kwargs.get('hmarker', 'o'),
                          ms=kwargs.get('hsize', 6))
        xyz = atoms._getCoords()
        min_.append(xyz.min(0))
        max_.append(xyz.max(0))

    show.set_xlabel('x')
    show.set_ylabel('y')
    show.set_zlabel('z')
    min_ = np.array(min_).min(0)
    max_ = np.array(max_).max(0)
    center = (max_ + min_) / 2
    half = (max_ - min_).max() / 2
    show.set_xlim3d(center[0]-half, center[0]+half)
    show.set_ylim3d(center[1]-half, center[1]+half)
    show.set_zlim3d(center[2]-half, center[2]+half)
    if kwargs.get('legend', False):
        show.legend(prop={'size': 10})
    if SETTINGS['auto_show']:
        showFigure()
    return show

ext: py

max_stars_repo_name: npabon/ProDy

max_issues_repo_name: npabon/ProDy

max_forks_repo_name: npabon/ProDy

reason: import numpy


text:
# An alternative to Ridge for regularizing linear regression.
import mglearn

import numpy as np

import matplotlib.pyplot as plt

from sklearn.linear_model import Lasso, Ridge
from sklearn.model_selection import train_test_split

X, y = mglearn.datasets.load_extended_boston()
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

lasso = Lasso()
lasso.fit(X_train, y_train)

print("Training set score: {:.2f}".format(lasso.score(X_train, y_train)))
print("Test set score: {:.2f}".format(lasso.score(X_test, y_test)))
print("Number of features used: {0}".format(np.sum(lasso.coef_ != 0)))

print('---\n{0}'.format(lasso.coef_))

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

lasso001 = Lasso(alpha=0.01, max_iter=100000)
lasso001.fit(X_train, y_train)

print("Training set score: {:.2f}".format(lasso001.score(X_train, y_train)))
print("Test set score: {:.2f}".format(lasso001.score(X_test, y_test)))
print("# of features used in Lasso: {0}".format(np.sum(lasso001.coef_ != 0)))

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

lasso0001 = Lasso(alpha=0.0001, max_iter=100000)
lasso0001.fit(X_train, y_train)
print("Training set score: {:.2f}".format(lasso0001.score(X_train, y_train)))
print("Test set score: {:.2f}".format(lasso0001.score(X_test, y_test)))
print("# of features used in Lasso: {0}".format(np.sum(lasso0001.coef_ != 0)))

### Plotting ###
ridge01 = Ridge(alpha=0.1).fit(X_train, y_train)

plt.plot(lasso.coef_, 's', label="Lasso alpha=1")
plt.plot(lasso001.coef_, '^', label="Lasso alpha=0.01")
plt.plot(lasso0001.coef_, 'v', label="Lasso alpha=0.0001")
plt.plot(ridge01.coef_, 'o', label="Ridge alpha=0.1")
plt.legend(ncol=2, loc=(0, 1.05))
plt.ylim(-25, 25)
plt.xlabel("Coefficient index")
plt.ylabel("Coefficient magnitude")
plt.show()

max_stars_repo_name: honchardev/Fun

max_issues_repo_name: honchardev/Fun

max_forks_repo_name: honchardev/Fun


text:
"""
Copyright (c) 2019 Yevheniia Soroka
Licensed under the MIT License
Author: Yevheniia Soroka
Email: ysoroka@cs.stonybrook.edu
Last modified: 18/12/2019

Usage:
Run this script to train GloVe model on Adobe tags.
"""

import os
import gluonnlp as nlp
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sb
from glove import Corpus, Glove
import glob
import pickle
from gensim.models import Word2Vec
from gensim.models import KeyedVectors

model_folder = "/nfs/bigfovea/add_disk0/eugenia/Emotion/wordembedding_models/"

#Creating a corpus object
corpus = Corpus()
#Training the corpus to generate the co occurence matrix which is used in GloVe
corpus.fit(lines, window=10)

# train the model
glove = Glove(no_components=5, learning_rate=0.05)
glove.fit(corpus.matrix, epochs=30, no_threads=4, verbose=True)
glove.add_dictionary(corpus.dictionary)

# save the model
glove.save(os.path.join(model_folder, 'glove_adobe.model'))

max_stars_repo_name: cvlab-stonybrook/Emotion-Prediction

max_issues_repo_name: cvlab-stonybrook/Emotion-Prediction

max_forks_repo_name: cvlab-stonybrook/Emotion-Prediction


text:
#from .optimizer import Optimizer, required
import torch
from torch.optim.optimizer import Optimizer, required
import numpy as np

from gutils import unit
from gutils import gproj
from gutils import clip_by_norm
from gutils import xTy
from gutils import gexp
from gutils import gpt
from gutils import gpt2

import pdb

class SGDG(Optimizer):
    r"""This optimizer updates variables with two different routines
        based on the boolean variable 'grassmann'.

        If grassmann is True, the variables will be updated by SGD-G proposed
        in 'Riemannian approach to batch normalization'.

        If grassmann is False, the variables will be updated by SGD.
        This routine was taken from https://github.com/pytorch/pytorch/blob/master/torch/optim/sgd.py.

        References:
           - Minhyung Cho and Jaehyung Lee, Riemannian approach to batch normalization
             (https://arxiv.org/abs/1709.09603)

    Args:
        params (iterable): iterable of parameters to optimize or dicts defining
            parameter groups

        -- common parameters
        lr (float): learning rate
        momentum (float, optional): momentum factor (default: 0)
        grassmann (bool, optional): whether to use SGD-G (default: False)

        -- parameters in case grassmann is False
        weight_decay (float, optional): weight decay (L2 penalty) (default: 0)
        dampening (float, optional): dampening for momentum (default: 0)
        nesterov (bool, optional): enables Nesterov momentum (default: False)

        -- parameters in case grassmann is True
        omega (float, optional): orthogonality regularization factor (default: 0)
        grad_clip (float, optional): threshold for gradient norm clipping (default: None)
    """

    def __init__(self, params, lr=required, momentum=0, dampening=0,
                 weight_decay=0, nesterov=False,
                 grassmann=False, omega=0, grad_clip=None):
        defaults = dict(lr=lr, momentum=momentum, dampening=dampening,
                        weight_decay=weight_decay, nesterov=nesterov,
                        grassmann=grassmann, omega=0, grad_clip=grad_clip)
        if nesterov and (momentum <= 0 or dampening != 0):
            raise ValueError("Nesterov momentum requires a momentum and zero dampening")
        super(SGDG, self).__init__(params, defaults)

    def __setstate__(self, state):
        super(SGDG, self).__setstate__(state)
        for group in self.param_groups:
            group.setdefault('nesterov', False)

    def step(self, closure=None):
        """Performs a single optimization step.

        Arguments:
            closure (callable, optional): A closure that reevaluates the model
                and returns the loss.
        """
        loss = None
        if closure is not None:
            loss = closure()

        for group in self.param_groups:
            momentum = group['momentum']
            grassmann = group['grassmann']

            if grassmann:
                grad_clip = group['grad_clip']
                omega = group['omega']

                for p in group['params']:
                    if p.grad is None:
                        continue

                    unity,_ = unit(p.data.view(p.size()[0],-1))
                    g = p.grad.data.view(p.size()[0],-1)

                    if omega != 0:
                      # L=|Y'Y-I|^2/2=|YY'-I|^2/2+c
                      # dL/dY=2(YY'Y-Y)
                      g.add_(2*omega, torch.mm(torch.mm(unity, unity.t()), unity) - unity)

                    h = gproj(unity, g)

                    if grad_clip is not None:
                        h_hat = clip_by_norm(h, grad_clip)
                    else:
                        h_hat = h

                    param_state = self.state[p]
                    if 'momentum_buffer' not in param_state:
                        param_state['momentum_buffer'] = torch.zeros(h_hat.size())
                        if p.is_cuda:
                          param_state['momentum_buffer'] = param_state['momentum_buffer'].cuda()

                    mom = param_state['momentum_buffer']
                    mom_new = momentum*mom - group['lr']*h_hat

                    p.data.copy_(gexp(unity, mom_new).view(p.size()))
                    mom.copy_(gpt(unity, mom_new))

            else:
                # This routine is from https://github.com/pytorch/pytorch/blob/master/torch/optim/sgd.py
                weight_decay = group['weight_decay']
                dampening = group['dampening']
                nesterov = group['nesterov']
                for p in group['params']:
                    if p.grad is None:
                        continue
                    d_p = p.grad.data
                    if weight_decay != 0:
                        d_p.add_(weight_decay, p.data)
                    if momentum != 0:
                        param_state = self.state[p]
                        if 'momentum_buffer' not in param_state:
                            buf = param_state['momentum_buffer'] = d_p.clone()
                        else:
                            buf = param_state['momentum_buffer']
                            buf.mul_(momentum).add_(1 - dampening, d_p)
                        if nesterov:
                            d_p = d_p.add(momentum, buf)
                        else:
                            d_p = buf

                    p.data.add_(-group['lr'], d_p)

        return loss

class AdamG(Optimizer):
    r"""This optimizer updates variables with two different routines
        based on the boolean variable 'grassmann'.

        If grassmann is True, the variables will be updated by Adam-G proposed
        in 'Riemannian approach to batch normalization'.

        If grassmann is False, the variables will be updated by SGD.
        This routine was taken from https://github.com/pytorch/pytorch/blob/master/torch/optim/sgd.py.

        References:
           - Minhyung Cho and Jaehyung Lee, Riemannian approach to batch normalization
             (https://arxiv.org/abs/1709.09603)

    Args:
        params (iterable): iterable of parameters to optimize or dicts defining
            parameter groups

        -- common parameters
        lr (float): learning rate
        momentum (float, optional): momentum factor (default: 0)
        grassmann (bool, optional): whether to use Adam-G (default: False)

        -- parameters in case grassmann is False
        weight_decay (float, optional): weight decay (L2 penalty) (default: 0)
        dampening (float, optional): dampening for momentum (default: 0)
        nesterov (bool, optional): enables Nesterov momentum (default: False)

        -- parameters in case grassmann is True
        beta2 (float, optional): the exponential decay rate for the second moment estimates (defulat: 0.99)
        epsilon (float, optional): a small constant for numerical stability (default: 1e-8)
        omega (float, optional): orthogonality regularization factor (default: 0)
        grad_clip (float, optional): threshold for gradient norm clipping (default: None)
    """

    def __init__(self, params, lr=required, momentum=0, dampening=0,
                 weight_decay=0, nesterov=False,
                 grassmann=False, beta2=0.99, epsilon=1e-8, omega=0, grad_clip=None):
        defaults = dict(lr=lr, momentum=momentum, dampening=dampening,
                        weight_decay=weight_decay, nesterov=nesterov,
                        grassmann=grassmann, beta2=beta2, epsilon=epsilon, omega=0, grad_clip=grad_clip)
        if nesterov and (momentum <= 0 or dampening != 0):
            raise ValueError("Nesterov momentum requires a momentum and zero dampening")
        super(AdamG, self).__init__(params, defaults)

    def __setstate__(self, state):
        super(AdamG, self).__setstate__(state)
        for group in self.param_groups:
            group.setdefault('nesterov', False)

    def step(self, closure=None):
        """Performs a single optimization step.

        Arguments:
            closure (callable, optional): A closure that reevaluates the model
                and returns the loss.
        """
        loss = None
        if closure is not None:
            loss = closure()

        for group in self.param_groups:
            grassmann = group['grassmann']

            if grassmann:
                beta1 = group['momentum']
                beta2 = group['beta2']
                epsilon = group['epsilon']
                grad_clip = group['grad_clip']
                omega = group['omega']

                for p in group['params']:
                    if p.grad is None:
                        continue

                    unity,_ = unit(p.data.view(p.size()[0],-1))
                    g = p.grad.data.view(p.size()[0],-1)

                    if omega != 0:
                      # L=|Y'Y-I|^2/2=|YY'-I|^2/2+c
                      # dL/dY=2(YY'Y-Y)
                      g.add_(2*omega, torch.mm(torch.mm(unity, unity.t()), unity) - unity)

                    h = gproj(unity, g)

                    if grad_clip is not None:
                        h_hat = clip_by_norm(h, grad_clip)
                    else:
                        h_hat = h

                    param_state = self.state[p]
                    if 'm_buffer' not in param_state:
                        size=p.size()
                        param_state['m_buffer'] = torch.zeros([size[0], int(np.prod(size[1:]))])
                        param_state['v_buffer'] = torch.zeros([size[0], 1])
                        if p.is_cuda:
                            param_state['m_buffer'] = param_state['m_buffer'].cuda()
                            param_state['v_buffer'] = param_state['v_buffer'].cuda()

                        param_state['beta1_power'] = beta1
                        param_state['beta2_power'] = beta2

                    m = param_state['m_buffer']
                    v = param_state['v_buffer']
                    beta1_power = param_state['beta1_power']
                    beta2_power = param_state['beta2_power']

                    mnew = beta1*m  + (1.0-beta1)*h_hat
                    vnew = beta2*v  + (1.0-beta2)*xTy(h_hat,h_hat)

                    alpha = np.sqrt(1.-beta2_power) / (1.-beta1_power)
                    deltas = mnew / vnew.add(epsilon).sqrt()
                    deltas.mul_(-alpha*group['lr'])

                    p.data.copy_(gexp(unity, deltas).view(p.size()))
                    m.copy_(gpt2(unity, mnew, deltas))
                    v.copy_(vnew)

                    param_state['beta1_power']*=beta1
                    param_state['beta2_power']*=beta2
            else:
                momentum = group['momentum']
                weight_decay = group['weight_decay']
                dampening = group['dampening']
                nesterov = group['nesterov']
                for p in group['params']:
                    if p.grad is None:
                        continue
                    d_p = p.grad.data
                    if weight_decay != 0:
                        d_p.add_(weight_decay, p.data)
                    if momentum != 0:
                        param_state = self.state[p]
                        if 'momentum_buffer' not in param_state:
                            buf = param_state['momentum_buffer'] = d_p.clone()
                        else:
                            buf = param_state['momentum_buffer']
                            buf.mul_(momentum).add_(1 - dampening, d_p)
                        if nesterov:
                            d_p = d_p.add(momentum, buf)
                        else:
                            d_p = buf

                    p.data.add_(-group['lr'], d_p)

        return loss

max_stars_repo_name: taralloc/deeprkm

max_issues_repo_name: taralloc/deeprkm

max_forks_repo_name: taralloc/deeprkm


text:
import numpy as np
import cvxpy

E1 = np.array([[-1.0,0.0,0.0],[-1.0,0.0,-1.0]])
E2 = np.array([[0.0,0.0],[-1.0,0.0]])
E3 = np.array([[0.0,0.0],[1.0,0.0]])
E4 = np.array([[1.0,0.0,-1.0],[1.0,0.0,0.0]])

A1 = np.array([[0.9,-0.1,0.0],[0.1,1,-0.02],[0.0,0.0,1.0]])
A2 = np.array([[1.0,-0.02],[0.02,0.9]])
A3 = np.array([[1.0,-0.02],[0.02,0.9]])
A4 = np.array([[0.9,-0.1,0.0],[0.1,1,0.02],[0.0,0.0,1.0]])

R = np.concatenate((np.eye(A2.shape[0]),np.zeros((A2.shape[0],1))),axis=1)
L = R.T

E2_ = np.matmul(E2,R)
E3_ = np.matmul(E3,R)
A2_ = np.matmul(L,np.matmul(A2,R))
A3_ = np.matmul(L,np.matmul(A3,R))

I2 = np.eye(2)
I3 = np.eye(3)

P1 = cvxpy.Symmetric(3)
P2 = cvxpy.Symmetric(2)
P3 = cvxpy.Symmetric(2)
P4 = cvxpy.Symmetric(3)

U1 = cvxpy.Symmetric(2)
U2 = cvxpy.Symmetric(2)
U3 = cvxpy.Symmetric(2)
U4 = cvxpy.Symmetric(2)

W1 = cvxpy.Symmetric(2)
W2 = cvxpy.Symmetric(2)
W3 = cvxpy.Symmetric(2)
W4 = cvxpy.Symmetric(2)

Q12 = cvxpy.Symmetric(2)
Q21 = cvxpy.Symmetric(2)
Q23 = cvxpy.Symmetric(2)
Q32 = cvxpy.Symmetric(2)
Q34 = cvxpy.Symmetric(2)
Q43 = cvxpy.Symmetric(2)

c=20.0
constraints = []
constraints.append(c*I2 << P2 - (E2.T)*U2*E2)
constraints.append(c*I2 << P3 - (E3.T)*U3*E3)
constraints.append((A2.T)*P2*A2 - P2 + (E2.T)*W2*E2 << c*I2)
constraints.append((A3.T)*P3*A3 - P3 + (E3.T)*W3*E3 << c*I2)
constraints.append(c*I3 << P1 - (E1.T)*U1*E1)
constraints.append(c*I3 << P4 - (E4.T)*U4*E4)
constraints.append((A1.T)*P1*A1 - P1 + (E1.T)*W1*E1 << c*I3)
constraints.append((A4.T)*P4*A4 - P4 + (E4.T)*W4*E4 << c*I3)
constraints.append((A2.T)*P3*A2 - P2 + (E2.T)*Q23*E2 << c*I2)
constraints.append((A3.T)*P2*A3 - P3 + (E3.T)*Q32*E3 << c*I2)
constraints.append((A1.T)*(L*P2*R)*A1 - P1 + (E1.T)*Q12*E1 << c*I3)
constraints.append((A4.T)*(L*P3*R)*A4 - P4 + (E4.T)*Q43*E4 << c*I3)
constraints.append(A2_.T*P1*A2_ - L*P2*R + E2_.T*Q21*E2_ << c*I3)
constraints.append(A3_.T*P4*A3_ - L*P3*R + E3_.T*Q34*E3_ << c*I3)
constraints.append(U1 > 0)
constraints.append(U2 > 0)
constraints.append(U3 > 0)
constraints.append(U4 > 0)
constraints.append(W1 > 0)
constraints.append(W2 > 0)
constraints.append(W3 > 0)
constraints.append(W4 > 0)
constraints.append(Q12 > 0)
constraints.append(Q21 > 0)
constraints.append(Q23 > 0)
constraints.append(Q32 > 0)
constraints.append(Q34 > 0)
constraints.append(Q43 > 0)

obj = cvxpy.Minimize(0)
Prob = cvxpy.Problem(obj,constraints)
Prob.solve()

max_stars_repo_name: MrRubyRed/DLS

max_issues_repo_name: MrRubyRed/DLS

max_forks_repo_name: MrRubyRed/DLS

reason: import numpy,import cvxpy


text:
import json

import dash_bootstrap_components as dbc
import dash_core_components as dcc
import dash_html_components as html
from dash.dependencies import Input, Output, State
from dash.dash import no_update
from networkx.readwrite import json_graph
import pandas as pd

from dash_app.app import app, server
from dash_app.pages import home, menu, vis, user_guide, about

app.layout = html.Div([
    dcc.Location(id='url', refresh=False),
    dbc.Navbar(
        id='top-bar',
        color='#1a3775',
        style={'height': '65px'},
        children=[
            dbc.Nav(
                dbc.Row(
                    [
                        html.A(
                            html.Img(src=app.get_asset_url('PaintDB-logo-small.png'),
                                     id='legend',
                                     height='60px',
                                     ),
                            href='/'),
                        dbc.NavLink('Build Network', href='/menu', active='exact'),
                        dbc.NavLink('Explore Network', href='/vis', id='explore', disabled=True),
                        dbc.DropdownMenu(
                            children=[
                                dbc.DropdownMenuItem('User Guide', href='/user_guide'),
                                dbc.DropdownMenuItem('About', href='/about'),
                            ],
                            nav=True,
                            in_navbar=True,
                            label='More'
                        ),
                    ],
                    align='center')
            )
        ],
    ),
    # Different pages are shown here
    html.Div(id='page-content'),

    # Hidden divs to store and share JSON data across callbacks and pages
    html.Div(id='hidden-bionetwork', style={'display': 'none'}),
    html.Div(id='network-parameters', style={'display': 'none'}),
    html.Div(id='node-details-df', style={'display': 'none'}),
    html.Div(id='enrichment-results', style={'display': 'none'}),
    html.Div(id='cyto-network', style={'display': 'none'}),
    html.Div(id='genes-of-interest', style={'display': 'none'}),
])

@app.callback(
    Output('explore', 'disabled'),
    [Input('hidden-bionetwork', 'children'),
     Input('enrichment-results', 'children')]
)
def enable_explore_tab(bio_network, enrichment=None):
    """Disables Explore tab if there is no network to explore."""
    return False if bio_network and enrichment else True

@app.callback(
    [Output('page-content', 'children'),
     Output('cyto-network', 'children')],
    [Input('url', 'pathname')],
    [State('hidden-bionetwork', 'children'),
     State('node-details-df', 'children'),
     State('enrichment-results', 'children'),
     State('network-parameters', 'children')]
)
def display_page(pathname, bio_network, json_df, json_enrichment_results, network_params):
    """Navigates to the selected app page. Generates vis layout depending on BioNetwork type."""
    if pathname == '/':
        return home.layout, no_update
    elif pathname == '/menu':
        return menu.layout, no_update
    elif pathname == '/vis':
        if bio_network:
            # Load JSON data
            network = json_graph.node_link_graph(json.loads(bio_network))
            network_df = pd.read_json(json_df)
            enrichment_results = pd.read_json(json_enrichment_results)
            # Create cyto network for visualization
            cyto_network = vis.make_cyto_elements(network)
            # Generate layout using generated data
            vis_layout = vis.make_vis_layout(network_df, enrichment_results, cyto_network, network_params)
            return vis_layout, json.dumps(cyto_network)
        else:
            return dbc.Alert('Warning: You need to build a network first.',
                             color='warning',
                             style={'display': 'inline-block', 'margin': '10px'}), no_update
    elif pathname == '/user_guide':
        return user_guide.layout, no_update
    elif pathname == '/about':
        return about.layout, no_update
    else:
        return dbc.Alert('Page not found.',
                         color='danger',
                         style={'display': 'inline-block', 'margin': '10px'}), no_update

if __name__ == '__main__':
    app.run_server(debug=True)

max_stars_repo_name: yavyx/PaIntDB

max_issues_repo_name: yavyx/PaIntDB

max_forks_repo_name: yavyx/PaIntDB

reason: from networkx


text:
import jax
from jax import lax, vmap
import jax.numpy as jnp
import haiku as hk
from jax_md import space, partition, nn, util, energy, smap
from jax_md.energy import multiplicative_isotropic_cutoff, _sw_radial_interaction, _sw_angle_interaction
from DiffTRe import custom_nn, custom_interpolate
from functools import partial

from typing import Callable, Tuple, Dict, Any
# Types
f32 = util.f32
f64 = util.f64
Array = util.Array

PyTree = Any
Box = space.Box
DisplacementFn = space.DisplacementFn
DisplacementOrMetricFn = space.DisplacementOrMetricFn

NeighborFn = partition.NeighborFn
NeighborList = partition.NeighborList

def stillinger_weber_energy(dr,
                            dR,
                            mask=None,
                            A=7.049556277,
                            B=0.6022245584,
                            p=4,
                            lam=21.0,
                            epsilon=2.16826,
                            gamma=1.2,
                            sigma=2.0951,
                            cutoff=1.8*2.0951,
                            three_body_strength=1.0):
    """
    Stillinger-Weber (SW) potential [1] which is commonly used to model
    silicon and similar systems. This function uses the default SW parameters
    from the original paper. The SW potential was originally proposed to
    model diamond in the diamond crystal phase and the liquid phase, and is
    known to give unphysical amorphous configurations [2, 3]. For this reason,
    we provide a three_body_strength parameter. Changing this number to 1.5
    or 2.0 has been know to produce more physical amorphous phase, preventing
    most atoms from having more than four nearest neighbors. Note that this
    function currently assumes nearest-image-convention.

    [1] Stillinger, Frank H., and Thomas A. Weber. "Computer simulation of
    local order in condensed phases of silicon." Physical review B 31.8
    (1985): 5262.
    [2] Holender, J. M., and G. J. Morgan. "Generation of a large structure
    (105 atoms) of amorphous Si using molecular dynamics." Journal of
    Physics: Condensed Matter 3.38 (1991): 7241.
    [3] Barkema, G. T., and Normand Mousseau. "Event-based relaxation of
    continuous disordered systems." Physical review letters 77.21 (1996): 4358.

    Args:
        dr: A ndarray of pairwise distances between particles
        dR: An ndarray of pairwise displacements between particles
        A: A scalar that determines the scale of two-body term
        B: Factor for radial power term
        p: Power in radial interaction
        lam: A scalar that determines the scale of the three-body term
        epsilon: A scalar that sets the energy scale
        gamma: Exponential scale in three-body term
        sigma: A scalar that sets the length scale
        cutoff: Cut-off value defined as sigma * a
        three_body_strength: A scalar that determines the relative strength
                             of the angular interaction
        mask: ndarray of size dr masking non-existing neighbors in neighborlist (if applicable)
    Returns:
        The Stilinger-Weber energy for a snapshot.
    """

    # initialize
    if mask is None:
        N = dr.shape[0]
        mask = jnp.ones([N, N])
        angle_mask = jnp.ones([N, N, N])
    else:  # for neighborlist input
        max_neighbors = mask.shape[-1]
        angle_mask1 = jnp.tile(jnp.expand_dims(mask, 1), [1, max_neighbors, 1])  # first pair
        angle_mask2 = jnp.tile(jnp.expand_dims(mask, -1), [1, 1, max_neighbors])  # second pair
        angle_mask = angle_mask1 * angle_mask2
    sw_radial_interaction = partial(_sw_radial_interaction, sigma=sigma, p=p, B=B, cutoff=cutoff)
    sw_angle_interaction = partial(_sw_angle_interaction, gamma=gamma, sigma=sigma, cutoff=cutoff)
    sw_three_body_term = vmap(vmap(vmap(sw_angle_interaction, (0, None)), (None, 0)), 0)

    # compute SW energy
    radial_interactions = sw_radial_interaction(dr) * mask
    angular_interactions = sw_three_body_term(dR, dR) * angle_mask
    first_term = A * jnp.sum(radial_interactions) / 2.0
    second_term = lam * jnp.sum(angular_interactions) / 2.0
    return epsilon * (first_term + three_body_strength * second_term)

def stillinger_weber_pair(displacement,
                          A=7.049556277,
                          B=0.6022245584,
                          p=4,
                          lam=21.0,
                          epsilon=2.16826,
                          gamma=1.2,
                          sigma=2.0951,
                          cutoff=1.8*2.0951,
                          three_body_strength=1.0):
    """Convenience wrapper to compute stilinger-weber energy over a system with variable parameters."""

    def compute_fn(R, **dynamic_kwargs):
        d = partial(displacement, **dynamic_kwargs)
        dR = space.map_product(d)(R, R)  # N x N x3 displacement matrix
        dr = space.distance(dR)  # N x N distances
        return stillinger_weber_energy(dr, dR, None, A, B, p, lam, epsilon, gamma, sigma, cutoff, three_body_strength)

    return compute_fn

def stillinger_weber_neighborlist(displacement,
                                  box_size=None,
                                  A=7.049556277,
                                  B=0.6022245584,
                                  p=4,
                                  lam=21.0,
                                  epsilon=2.16826,
                                  gamma=1.2,
                                  sigma=2.0951,
                                  cutoff=1.8*2.0951,
                                  three_body_strength=1.0,
                                  dr_threshold=0.1,
                                  capacity_multiplier=1.25,
                                  initialize_neighbor_list=True):
    """Convenience wrapper to compute stilinger-weber energy using a neighbor list."""

    def energy_fn(R, neighbor, **dynamic_kwargs):
        d = partial(displacement, **dynamic_kwargs)
        N = R.shape[0]
        mask = neighbor.idx != N
        R_neigh = R[neighbor.idx]
        dR = space.map_neighbor(d)(R, R_neigh)
        dr = space.distance(dR)
        return stillinger_weber_energy(dr, dR, mask, A, B, p, lam, epsilon, gamma, sigma, cutoff, three_body_strength)

    if initialize_neighbor_list:
        assert box_size is not None
        neighbor_fn = partition.neighbor_list(displacement, box_size, cutoff, dr_threshold,
                                              capacity_multiplier=capacity_multiplier)
        return neighbor_fn, energy_fn

    return energy_fn

def generic_repulsion(dr: Array,
                      sigma: Array=1.,
                      epsilon: Array=1.,
                      exp: Array=12.,
                      **dynamic_kwargs) -> Array:
    """
    Repulsive interaction between soft sphere particles: U = epsilon * (sigma / r)**exp.

    Args:
      dr: An ndarray of pairwise distances between particles.
      sigma: Repulsion length scale
      epsilon: Interaction energy scale
      exp: Exponent specifying interaction stiffness

    Returns:
      Array of energies
    """

    dr = jnp.where(dr > 1.e-7, dr, 1.e7)  # save masks dividing by 0
    idr = (sigma / dr)
    U = epsilon * idr ** exp
    return U

def generic_repulsion_pair(displacement_or_metric: DisplacementOrMetricFn,
                     species: Array=None,
                     sigma: Array=1.0,
                     epsilon: Array=1.0,
                     exp: Array=12.,
                     r_onset: Array = 2.0,
                     r_cutoff: Array = 2.5,
                     per_particle: bool=False):
    """Convenience wrapper to compute generic repulsion energy over a system."""
    sigma = jnp.array(sigma, dtype=f32)
    epsilon = jnp.array(epsilon, dtype=f32)
    exp = jnp.array(exp, dtype=f32)
    r_onset = jnp.array(r_onset, dtype=f32)
    r_cutoff = jnp.array(r_cutoff, dtype=f32)

    return smap.pair(
        multiplicative_isotropic_cutoff(generic_repulsion, r_onset, r_cutoff),
        space.canonicalize_displacement_or_metric(displacement_or_metric),
        species=species,
        sigma=sigma,
        epsilon=epsilon,
        exp=exp,
        reduce_axis=(1,) if per_particle else None)

def generic_repulsion_neighborlist(displacement_or_metric: DisplacementOrMetricFn,
                     box_size: Box=None,
                     species: Array=None,
                     sigma: Array=1.0,
                     epsilon: Array=1.0,
                     exp: Array=12.,
                     r_onset: Array = 0.9,
                     r_cutoff: Array = 1.,
                     dr_threshold: float=0.2,
                     per_particle: bool=False,
                     capacity_multiplier: float=1.25,
                     initialize_neighbor_list: bool=True):
    """
    Convenience wrapper to compute generic repulsion energy over a system with neighborlist.

    Provides option not to initialize neighborlist. This is useful if energy function needs
    to be initialized within a jitted function.
    """

    sigma = jnp.array(sigma, dtype=f32)
    epsilon = jnp.array(epsilon, dtype=f32)
    exp = jnp.array(exp, dtype=f32)
    r_onset = jnp.array(r_onset, dtype=f32)
    r_cutoff = jnp.array(r_cutoff, dtype=f32)

    energy_fn = smap.pair_neighbor_list(
      multiplicative_isotropic_cutoff(generic_repulsion, r_onset, r_cutoff),
      space.canonicalize_displacement_or_metric(displacement_or_metric),
      species=species,
      sigma=sigma,
      epsilon=epsilon,
      exp=exp,
      reduce_axis=(1,) if per_particle else None)

    if initialize_neighbor_list:
        assert box_size is not None
        neighbor_fn = partition.neighbor_list(displacement_or_metric, box_size, r_cutoff, dr_threshold,
                                              capacity_multiplier=capacity_multiplier)
        return neighbor_fn, energy_fn

    return energy_fn

def tabulated(dr: Array, spline: Callable[[Array], Array], **unused_kwargs) -> Array:
    """
    Tabulated radial potential between particles given a spline function.

    Args:
        dr: An ndarray of pairwise distances between particles
        spline: A function computing the spline values at a given pairwise distance

    Returns:
        Array of energies
    """

    return spline(dr)

def tabulated_pair(displacement_or_metric: DisplacementOrMetricFn,
                   x_vals: Array,
                   y_vals: Array,
                   degree: int=3,
                   r_onset: Array=0.9,
                   r_cutoff: Array=1.,
                   species: Array = None,
                   per_particle: bool=False) -> Callable[[Array], Array]:
    """Convenience wrapper to compute tabulated energy over a system."""
    x_vals = jnp.array(x_vals, f32)
    y_vals = jnp.array(y_vals, f32)
    r_onset = jnp.array(r_onset, f32)
    r_cutoff = jnp.array(r_cutoff, f32)

    spline = custom_interpolate.MonotonicInterpolate(x_vals, y_vals)
    tabulated_partial = partial(tabulated, spline=spline)

    return smap.pair(
      multiplicative_isotropic_cutoff(tabulated_partial, r_onset, r_cutoff),
      space.canonicalize_displacement_or_metric(displacement_or_metric),
      species=species,
      reduce_axis=(1,) if per_particle else None)

def tabulated_neighbor_list(displacement_or_metric: DisplacementOrMetricFn,
                            x_vals: Array,
                            y_vals: Array,
                            box_size: Box,
                            degree: int=3,
                            r_onset: Array=0.9,
                            r_cutoff: Array=1.,
                            dr_threshold: Array=0.2,
                            species: Array = None,
                            capacity_multiplier: float=1.25,
                            initialize_neighbor_list: bool=True,
                            per_particle: bool=False) -> Callable[[Array], Array]:
    """
    Convenience wrapper to compute tabulated energy using a neighbor list.

    Provides option not to initialize neighborlist. This is useful if energy function needs
    to be initialized within a jitted function.
    """

    x_vals = jnp.array(x_vals, f32)
    y_vals = jnp.array(y_vals, f32)
    box_size = jnp.array(box_size, f32)
    r_onset = jnp.array(r_onset, f32)
    r_cutoff = jnp.array(r_cutoff, f32)
    dr_threshold = jnp.array(dr_threshold, f32)

    # Note: cannot provide the spline parameters via kwargs because only per-perticle parameters are supported
    spline = custom_interpolate.MonotonicInterpolate(x_vals, y_vals)
    tabulated_partial = partial(tabulated, spline=spline)

    energy_fn = smap.pair_neighbor_list(
      multiplicative_isotropic_cutoff(tabulated_partial, r_onset, r_cutoff),
      space.canonicalize_displacement_or_metric(displacement_or_metric),
      species=species,
      reduce_axis=(1,) if per_particle else None)

    if initialize_neighbor_list:
        neighbor_fn = partition.neighbor_list(displacement_or_metric, box_size, r_cutoff, dr_threshold,
                                              capacity_multiplier=capacity_multiplier)
        return neighbor_fn, energy_fn
    return energy_fn

def DimeNetPP_neighborlist(displacement: DisplacementFn,
                           R_test: Array,
                           neighbor_test,
                           r_cutoff: float,
                           embed_size: int = 32,
                           n_interaction_blocks: int = 4,
                           num_residual_before_skip: int = 1,
                           num_residual_after_skip: int = 2,
                           out_embed_size=None,
                           type_embed_size=None,
                           angle_int_embed_size=None,
                           basis_int_embed_size: int = 8,
                           num_dense_out: int = 3,
                           num_RBF: int = 6,
                           num_SBF: int = 7,
                           activation=jax.nn.swish,
                           envelope_p: int = 6,
                           init_kwargs: Dict[str, Any] = None,
                           n_species: int = 10,
                           max_angle_multiplier: float = 1.25,
                           max_edge_multiplier: float = 1.25
                           ) -> Tuple[nn.InitFn, Callable[[PyTree, Array], Array]]:
    """
    Convenience wrapper around DimeNetPPEnergy model using a neighbor list.

    The defaults correspond to the orinal values of DimeNet++, except for
    embedding sizes that are reduced by a factor of 4 for efficiency reasons.
    Computes molecular inputs: pairwise distances and angular triplets.
    For computational efficiency, a sparse edge and angle representation is
    implemented. A sample neighbor list to estimate max_edges and max_angles
    is therefore necessary.

    Args:
        displacement: Displacement function
        R_test: Sample state.position for initialization of network and estimate of max_edges / max_angles
        neighbor_test: Sample neighborlist for estimate of max_edges / max_angles
        r_cutoff: Cut-off of DimeNetPP and neighbor list
        embed_size: Size of embeddings; will scale interaction and output embedding sizes accordingly
        n_interaction_blocks: Number of interaction blocks
        num_residual_before_skip: Number of residual blocks before the skip connection
        num_residual_after_skip:Number of residual blocks after the skip connection
        out_embed_size: Embedding size of output block; If None is 2 * embed_size
        type_embed_size: Embedding size of type embedding; If None is 0.5 * embed_size
        angle_int_embed_size: Embedding size of angle interactions; If None is 0.5 * embed_size
        basis_int_embed_size: Embedding size of Bessel basis embedding in interaction block
        num_dense_out: Number of fully-connected layers in output block
        num_RBF: Number of radial Bessel embedding functions
        num_SBF: Number of spherical Bessel embedding functions
        activation: Activation function
        envelope_p: Power of envelope polynomial
        init_kwargs: Kwargs for initializaion of MLPs
        n_species: Number of species
        max_angle_multiplier: Multiplier for initial estimate of maximum triplets
        max_edge_multiplier: Multiplier for initial estimate of maximum edges

    Returns:
        A tuple of 2 functions: A init_fn that initializes the model parameters
        and an energy function that computes the energy for a particular state
        given model parameters.
    """
    r_cutoff = jnp.array(r_cutoff, dtype=f32)
    # Estimate maximum number of edges and angles for efficient computations during simulation:
    # neighbor.idx: an index j in row i encodes a directed edge from particle j to particle i. edge_idx[i, j]: j->i
    # if j == N: encodes no edge exists;  would index out-of-bounds, but in jax last element is returned in this case
    neighbor_displacement_fn = space.map_neighbor(displacement)
    R_neigh_test = R_test[neighbor_test.idx]
    test_displacement = neighbor_displacement_fn(R_test, R_neigh_test)
    pair_distances_test = space.distance(test_displacement)
    # adds all edges > cut-off to masked edges
    edge_idx_test = jnp.where(pair_distances_test < r_cutoff, neighbor_test.idx, R_test.shape[0])
    _, _, _, _, (n_edges_init, n_angles_init) = custom_nn.sparse_representation(pair_distances_test, edge_idx_test)
    max_angles = jnp.int32(jnp.ceil(n_angles_init * max_edge_multiplier))
    max_edges = jnp.int32(jnp.ceil(n_edges_init * max_angle_multiplier))

    if init_kwargs is None:
        init_kwargs = {
          'w_init': custom_nn.OrthogonalVarianceScalingInit(scale=1.),
          'b_init': hk.initializers.Constant(0.),
        }

    @hk.without_apply_rng
    @hk.transform
    def model(R: Array, neighbor, species=None, **dynamic_kwargs) -> Array:
        N, max_neighbors = neighbor.idx.shape
        if species is None:  # dummy species to allow streamlined use of different species
            species = jnp.zeros(N, dtype=jnp.int32)

        dynamic_displacement = partial(displacement, **dynamic_kwargs)  # necessary for pressure computation
        dyn_neighbor_displacement_fn = space.map_neighbor(dynamic_displacement)

        # compute pairwise distances
        R_neigh = R[neighbor.idx]
        pair_displacement = dyn_neighbor_displacement_fn(R, R_neigh)
        pair_distances = space.distance(pair_displacement)

        # compute adjacency matrix via neighbor_list, then build sparse representation to avoid part of padding overhead
        edge_idx_ji = jnp.where(pair_distances < r_cutoff, neighbor.idx, N)  # adds all edges > cut-off to masked edges
        pair_distances_sparse, pair_connections, angle_idxs, angular_connectivity, (n_edges, n_angles) = \
            custom_nn.sparse_representation(pair_distances, edge_idx_ji, max_edges, max_angles)
        too_many_edges_error_code = lax.cond(jnp.bitwise_or(n_edges > max_edges, n_angles > max_angles),
                                             lambda _: True, lambda _: False, n_edges)
        # TODO: return too_many_edges_error_code

        idx_i, idx_j, pair_mask = pair_connections
        # cutoff all non existing edges: are encoded as 0 by rbf envelope; non-existing angles will also be encoded as 0
        pair_distances_sparse = jnp.where(pair_mask[:, 0], pair_distances_sparse, 2. * r_cutoff)
        angles = custom_nn.angles_triplets(R, dynamic_displacement, angle_idxs, angular_connectivity)
        net = custom_nn.DimeNetPPEnergy(r_cutoff,
                                        n_particles=N,
                                        n_species=n_species,
                                        embed_size=embed_size,
                                        n_interaction_blocks=n_interaction_blocks,
                                        num_residual_before_skip=num_residual_before_skip,
                                        num_residual_after_skip=num_residual_after_skip,
                                        out_embed_size=out_embed_size,
                                        type_embed_size=type_embed_size,
                                        angle_int_embed_size=angle_int_embed_size,
                                        basis_int_embed_size=basis_int_embed_size,
                                        num_dense_out=num_dense_out,
                                        num_RBF=num_RBF,
                                        num_SBF=num_SBF,
                                        activation=activation,
                                        envelope_p=envelope_p,
                                        init_kwargs=init_kwargs)
        gnn_energy = net(pair_distances_sparse, angles, species, pair_connections, angular_connectivity)
        gnn_energy = gnn_energy[0]  # the net returns a 1D array as output, but grad needs a scalar for differentiation
        return gnn_energy

    return model.init, model.apply

max_stars_repo_name: moradza/difftre

max_issues_repo_name: moradza/difftre

max_forks_repo_name: moradza/difftre

reason: import jax,from jax


text:
import argparse
import os
import torch
from models import create_ensemble
from models.model_configs import config_from_string
import augmentations

import os
import argparse
import math
import numpy as np

import torch
import torch.autograd as autograd
import torch.utils.data as data
import torchvision
import torchvision.datasets.folder
import torchvision.transforms as transforms

from dataset import Dataset

parser = argparse.ArgumentParser(description='Defence')
parser.add_argument('--input_dir', metavar='DIR',
                    help='Input directory with images.')
parser.add_argument('--output_file', metavar='FILE', default='',
                    help='Output file to save labels.')

def main():
    args = parser.parse_args()

    tf = transforms.Compose([
        transforms.ToTensor(),
    ])

    dataset = Dataset(args.input_dir, transform=tf)
    loader = data.DataLoader(dataset, batch_size=8, shuffle=False)

    model_config_srings = ['Resnet18', 'InceptionResnetV2']
    cfgs = [config_from_string(s) for s in model_config_srings]

    ensemble = create_ensemble(cfgs, [1.0 for _ in cfgs])
    augment = augmentations.AugmentationComposer([
        augmentations.RandomCrop(269),
        augmentations.Mirror(0.5),
        augmentations.Blur(0.5, 0.5)
    ])

    for cfg, model in zip(cfgs, ensemble.models):
        checkpoint_path = os.path.join('/checkpoints/',cfg['checkpoint_file'])
        checkpoint = torch.load(checkpoint_path)
        if isinstance(checkpoint, dict) and 'state_dict' in checkpoint:
            model.get_core_model().load_state_dict(checkpoint['state_dict'])
        else:
            model.get_core_model().load_state_dict(checkpoint)
    ensemble.cuda()
    ensemble.eval()

    outputs = []
    for batch_idx, (input, _) in enumerate(loader):
        input = input.cuda()
        input_var = autograd.Variable(input, volatile=True)
        labels = ensemble(augment(input_var))
        labels = labels.max(1)[1]
        outputs.append(labels.data.cpu().numpy())
    outputs = np.concatenate(outputs, axis=0)

    with open(args.output_file, 'w') as out_file:
        filenames = dataset.filenames()
        for filename, label in zip(filenames, outputs):
            filename = os.path.basename(filename)
            out_file.write('{0},{1}\n'.format(filename, label))

if __name__ == '__main__':
    main()


text:
import cv2
import os
from PIL import Image
import numpy as np

video_path = 'Results/2.mp4'
data_path = 'data'

cam = cv2.VideoCapture(video_path)
try:
    if not os.path.exists(data_path):
        os.makedirs(data_path)
except OSError:
    print('Error')

currentframe = 0
while(True):
    ret, f = cam.read()
    if ret:
         name =data_path + '/frame' + str(currentframe) + '.png'
         cv2.imwrite(name, f)
         currentframe += 1
         #print('COUNT: ', currentframe)
    else:
         break
cam.release()
cv2.destroyAllWindows()


text:
# Copyright 2019 Cris Almodovar
#
# 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.

import grpc
from laservec.laservec_pb2_grpc import LaserGrpcApiStub
from laservec.laservec_pb2 import VectorizeRequest
import numpy as np

class LASER(object):
    def __init__(self, url: str):
        self.url = url
        self.channel = grpc.insecure_channel(url)
        self.stub = LaserGrpcApiStub(self.channel)

    def close(self):
        self.stub = None
        self.channel.close()
        self.channel = None

    def vectorize(self, text: str, lang: str = None):
        assert self.channel is not None
        assert self.stub is not None

        req = VectorizeRequest(text=text, lang=lang)
        res = self.stub.vectorize(req)
        if res:
            return np.array(res.embedding), res.lang
        return None, None

    def __enter__(self):
        return self

    def __exit__(self, exc_type, exc_val, exc_tb):
        self.close()

max_stars_repo_name: cris-almodovar/laservec

max_issues_repo_name: cris-almodovar/laservec

max_forks_repo_name: cris-almodovar/laservec


text:
import datetime
import multiprocessing
import picamera
import picamera.array
import cv2
import numpy as np

# Define
MOTION_DETECTION_TH = 0
BINARIZE_TH         = 30
PIXEL_VAL_MAX       = 255

RESOLUTION          = ( 1920, 1080 )
FRAMERATE           = 30

class CameraController():
    '''
    カメラ制御クラス
    PiCameraインスタンスとDetectionControllerインスタンスを持ち、
    動体の検出に応じて録画を開始/停止する
    '''

    def __init__( self ):
        self.camera            = picamera.PiCamera( resolution=RESOLUTION, framerate=FRAMERATE )
        self.detection_ctrlr   = DetectionController( self.camera )

    def control_recording( self ):
        '''
        動体の検出に応じて録画の開始/停止を制御する
        '''
        try:
            self.detection_ctrlr.set_detection( 'frame subtraction' )

            while True:
                if self.detection_ctrlr.detection_is_enabled == True:
                    is_detected = self.detection_ctrlr.detect()

                    if is_detected == True:
                        if self.camera.recording == False:
                            self.camera.start_recording( 'video'+datetime.datetime.now().strftime('%Y-%m-%d-%H%M%S')+'.h264')
                            print( 'Start recording' )
                    else:
                        if self.camera.recording == True:
                            self.camera.stop_recording()
                            print( 'Stop recording' )

        finally:
            if self.camera.recording == True:
                self.camera.stop_recording()

            self.detection_ctrlr.fin_detection()
            self.camera.close()
            print( 'Done closing' )

class DetectionController():
    ''' Detection Controllerクラス（Proxyパターン）'''

    def __init__( self, camera ):
        self.camera               = camera
        self.detection_is_enabled = True

    def set_detection( self, algorithm=None ):
        ''' 動体検出アルゴリズムを設定する '''

        if algorithm == None:
            print( 'Please select detection algorithm' )
            return

        # フレーム差分法
        if algorithm == 'frame subtraction':
            self.algorithm = FrameSubtraction( self.camera )
            print( 'Initialized detection' )

        print( 'Done setting up detection' )

    def fin_detection( self ):
        ''' Finalize detection '''

        self.algorithm.fin()
        print( 'Finalized detection' )

    def detect( self ):
        ''' Detect '''
        # return True/False

        return self.algorithm.func()

class FrameSubtraction():
    ''' Frame Subtraction algorithm '''

    def __init__( self, camera ):
        self.camera = camera
        self.frame1 = picamera.array.PiRGBArray( self.camera )
        self.frame2 = picamera.array.PiRGBArray( self.camera )
        self.frame3 = picamera.array.PiRGBArray( self.camera )

    def fin( self ):
        ''' Finalizer of frame subtraction '''

        self.frame1.close()
        self.frame2.close()
        self.frame3.close()

    def func( self ):
        ''' フレーム差分法による動体検出アルゴリズム '''
        # return True/False

        # フレームを３枚キャプチャ
        self.camera.capture( self.frame1, 'bgr', use_video_port=True )
        self.camera.capture( self.frame2, 'bgr', use_video_port=True )
        self.camera.capture( self.frame3, 'bgr', use_video_port=True )

        # RGBからGrayへ変換
        frame1_gray = cv2.cvtColor( self.frame1.array, cv2.COLOR_BGR2GRAY )
        frame2_gray = cv2.cvtColor( self.frame2.array, cv2.COLOR_BGR2GRAY )
        frame3_gray = cv2.cvtColor( self.frame3.array, cv2.COLOR_BGR2GRAY )

        # フレーム差分の計算
        dif_img1 = cv2.absdiff( frame1_gray, frame2_gray )
        dif_img2 = cv2.absdiff( frame2_gray, frame3_gray )

        dif_imga = cv2.bitwise_and( dif_img1, dif_img2 )

        # 2値化
        ret, bin_img = cv2.threshold( dif_imga, BINARIZE_TH, PIXEL_VAL_MAX, cv2.THRESH_BINARY)
        bin_img = cv2.medianBlur( bin_img, 7 )

        if np.sum(bin_img) > MOTION_DETECTION_TH:
            is_detected = True
            print( 'Detected' )
        else:
            is_detected = False
            print( 'Not Detected' )

        # PiCameraモジュールの決まりで、frame1, frame2, frame3を再利用する前に truncate(0)する必要があるらしい
        self.frame1.truncate( 0 )
        self.frame2.truncate( 0 )
        self.frame3.truncate( 0 )

        return is_detected

max_stars_repo_name: kzkhmd/cappi

max_issues_repo_name: kzkhmd/cappi

max_forks_repo_name: kzkhmd/cappi


text:
#!/usr/bin/env python3

import os
import time
import json
import logging
import argparse
import subprocess

import shutil

import xml.etree.ElementTree as ET
from xml.dom import minidom

from pathlib import Path

import numpy as np

from ase import Atoms
from ase.io import read, write
from ase.io.vasp import write_vasp
from ase.calculators.vasp import Vasp2
from ase.constraints import FixAtoms

# customised xsd reader
#import xsd_plus as xsd2
#from xsd_plus import read_xsd

vasp_params = {
    # INCAR
    "nwrite": 2,
    "istart": 0,
    "lcharg": False,
    "lwave": False,
    "lorbit": 10,
    "npar": 4,
    "xc": "pbe",
    "encut": 400,
    "prec": "Normal",
    "ediff": 1E-5,
    "nelm": 180,
    "nelmin": 6,
    "ispin": 1,
    "ismear": 1,
    "sigma": 0.2,
    "algo": "Fast",
    "lreal": "Auto",
    "isym": 0,
    "ediffg": -0.05,
    "nsw": 200,
    "ibrion": 2,
    "isif": 2,
    "potim": 0.02,
}

def read_xsd2(fd):
    tree = ET.parse(fd)
    root = tree.getroot()

    atomtreeroot = root.find('AtomisticTreeRoot')
    # if periodic system
    if atomtreeroot.find('SymmetrySystem') is not None:
        symmetrysystem = atomtreeroot.find('SymmetrySystem')
        mappingset = symmetrysystem.find('MappingSet')
        mappingfamily = mappingset.find('MappingFamily')
        system = mappingfamily.find('IdentityMapping')

        coords = list()
        cell = list()
        formula = str()

        names = list()
        restrictions = list()

        for atom in system:
            if atom.tag == 'Atom3d':
                symbol = atom.get('Components')
                formula += symbol

                xyz = atom.get('XYZ')
                if xyz:
                    coord = [float(coord) for coord in xyz.split(',')]
                else:
                    coord = [0.0, 0.0, 0.0]
                coords.append(coord)

                name = atom.get('Name')
                if name:
                    pass # find name
                else:
                    name = symbol + str(len(names)+1) # None due to copy atom
                names.append(name)

                restriction = atom.get('RestrictedProperties')
                if restriction:
                    if restriction == 'FractionalXYZ':
                        restrictions.append(True)
                    else:
                        raise ValueError('unknown RestrictedProperties')
                else:
                    restrictions.append(False)
            elif atom.tag == 'SpaceGroup':
                avec = [float(vec) for vec in atom.get('AVector').split(',')]
                bvec = [float(vec) for vec in atom.get('BVector').split(',')]
                cvec = [float(vec) for vec in atom.get('CVector').split(',')]

                cell.append(avec)
                cell.append(bvec)
                cell.append(cvec)

        atoms = Atoms(formula, cell=cell, pbc=True)
        atoms.set_scaled_positions(coords)

        # add constraints
        fixed_indices = [idx for idx, val in enumerate(restrictions) if val]
        if fixed_indices:
            atoms.set_constraint(FixAtoms(indices=fixed_indices))

        # add two atoms constrained optimisation
        constrained_indices = [
            idx for idx, name in enumerate(names) if name.endswith('_c')
        ]
        if constrained_indices:
            assert len(constrained_indices) == 2
            atoms.info['copt'] = constrained_indices

        return atoms
        # if non-periodic system
    elif atomtreeroot.find('Molecule') is not None:
        system = atomtreeroot.find('Molecule')

        coords = list()
        formula = str()

        for atom in system:
            if atom.tag == 'Atom3d':
                symbol = atom.get('Components')
                formula += symbol

                xyz = atom.get('XYZ')
                coord = [float(coord) for coord in xyz.split(',')]
                coords.append(coord)

        atoms = Atoms(formula, pbc=False)
        atoms.set_scaled_positions(coords)
        return atoms

def create_copt(atoms):
    """"""

    return content

def create_slurm(directory, partition='k2-hipri', time='3:00:00', ncpus='32'):
    """create job script"""
    content = "#!/bin/bash -l \n"
    content += "#SBATCH --partition=%s        # queue\n" %partition
    content += "#SBATCH --job-name=%s         # Job name\n" %Path(directory).name
    content += "#SBATCH --time=%s              # Time limit hrs:min:sec\n" %time
    content += "#SBATCH --nodes=1-2                  # Number of nodes\n"
    content += "#SBATCH --ntasks=%s                  # Number of cores\n" %ncpus
    content += "#SBATCH --cpus-per-task=1            # Number of cores per MPI task \n"
    content += "#SBATCH --mem=10G                    # Number of cores per MPI task \n"
    content += "#SBATCH --output=slurm.o%j           # Standard output and error log\n"
    content += "#SBATCH --error=slurm.e%j            # Standard output and error log\n"
    content += "\n"
    content += "module purge\n"
    content += "module load services/s3cmd\n"
    content += "module load libs/intel/2016u1\n"
    content += "module load libs/intel-mkl/2016u1/bin\n"
    content += "module load mpi/intel-mpi/2016u1/bin\n"
    content += "\n"
    content += 'echo `date "+%Y-%m-%d %H:%M:%S"` `pwd` >> $HOME/submitted\n'
    content += "mpirun -n 32 /mnt/scratch/chemistry-apps/dkb01416/vasp/installed/intel-2016/5.4.1-TS/vasp_std 2>&1 > vasp.out\n"
    content += 'echo `date "+%Y-%m-%d %H:%M:%S"` `pwd` >> $HOME/finished\n'

    with open(os.path.join(directory,'vasp.slurm'), 'w') as fopen:
        fopen.write(content)

    return

def create_vasp_inputs(atoms, incar=None, directory='vasp-test'):
    """only for molecules"""
    # pseudo
    pp_path = "/mnt/scratch/chemistry-apps/dkb01416/vasp/PseudoPotential"
    if 'VASP_PP_PATH' in os.environ.keys():
        os.environ.pop('VASP_PP_PATH')
    os.environ['VASP_PP_PATH'] = pp_path

    # vdw
    vdw_envname = 'ASE_VASP_VDW'
    vdw_path = "/mnt/scratch/chemistry-apps/dkb01416/vasp/pot"
    if vdw_envname in os.environ.keys():
        _ = os.environ.pop(vdw_envname)
    os.environ[vdw_envname] = vdw_path

    # ===== initialise =====
    calc = Vasp2(
        command = "vasp_std",
        directory = 'dummy',
        **vasp_params,
    )

    # convert to pathlib usage
    if directory != Path.cwd() and not directory.exists():
        directory.mkdir()

    calc.initialize(atoms)

    calc.string_params['system'] = Path(directory).name

    content = '\n>>>>> Modified ASE for VASP <<<<<\n'
    content += '    directory -> %s\n' %directory
    print(content)

    # ===== POSCAR =====
    poscar = os.path.join(directory, 'POSCAR')
    write_vasp(poscar, atoms, direct=True, vasp5=True)

    from collections import Counter
    cons = atoms.constraints
    assert len(cons) == 1
    cons_indices = cons[0].get_indices()
    symbols = atoms.get_chemical_symbols()
    all_atoms = Counter(symbols)
    fixed_symbols = [symbols[i] for i in cons_indices]
    fixed_atoms = Counter(fixed_symbols)

    natoms = len(atoms)
    nfixed = np.sum(list(fixed_atoms.values()))

    content = '\nPOSCAR -->\n'
    content += '    \033[4m Element       Numbers      \033[0m\n'
    for sym in all_atoms.keys():
        content += '         %2s \033[1;33m%4d\033[0m \033[32m%4d\033[0m(T)\033[31m%4d\033[0m(F)\n'\
                %(sym, all_atoms[sym], all_atoms[sym]-fixed_atoms[sym], fixed_atoms[sym])
    content += '    \033[4m                            \033[0m\n'
    content += '      %2s \033[1;33m%4d\033[0m \033[32m%4d\033[0m(T)\033[31m%4d\033[0m(F)\n'\
                %('Total', natoms, natoms-nfixed, nfixed)

    print(content)

    # ===== KPOINTS =====
    kpts = np.linalg.norm(atoms.cell, axis=1).tolist()
    kpts = [int(20./k)+1 for k in kpts]

    calc.input_params['gamma'] = True
    calc.input_params['kpts'] = kpts
    calc.write_kpoints(atoms=atoms, directory=directory)

    content = 'KPOINTS -->\n'
    content += "     Set k-point -> \033[1;35m{} {} {}\033[0m\n".format(*kpts)
    print(content)

    # ===== INCAR and 188 =====
    content = 'INCAR -->\n'

    if incar:
        calc.read_incar(incar)
        calc.string_params['system'] = Path(directory).name
        content += '    read template from %s\n' %incar

    copt = atoms.info.get('copt', None)
    if copt:
        # If constrained get distance and create fort.188
        symbols = atoms.get_chemical_symbols()
        ca, cb = atoms.info['copt'][0], atoms.info['copt'][1]
        pt1, pt2 = atoms.positions[ca], atoms.positions[cb]

        # Use Ax=b convert to cartisan coordinate
        distance = np.linalg.norm(pt1-pt2)

        # Create fort.188
        ts_content = '1\n3\n6\n4\n0.04\n%-5d%-5d%f\n0\n' % \
            (ca+1, cb+1, distance)

        with open(os.path.join(directory,'fort.188'), 'w') as f:
            f.write(ts_content)

        calc.int_params['ibrion'] = 1

        content += '\n'
        content += "     fort.188 has been created.\n"
        content += '     ' + '-'*20 + '\n'
        content += "     atom number: {:<5d}{:<5d}\n".format(ca+1, cb+1)
        content += "     atom name: {} {}\n".format(symbols[ca], symbols[cb])
        content += "     distance: {:f}\n".format(distance)
        content += '     ' + '-'*20 + '\n'

        # Set IBRION = 1
        content += '     Note: IBRION has been set to 1.\n'

    print(content)

    calc.write_incar(atoms, directory=directory)

    # ===== POTCAR =====
    calc.write_potcar(directory=directory)
    content = "POTCAR -->\n"
    content += "    Using POTCAR from %s\n" %pp_path
    print(content)

    # ===== VDW =====
    calc.copy_vdw_kernel(directory=directory)

    content = "VDW -->\n"
    content += "    Using vdW Functional as %s\n" %calc.string_params['gga']
    print(content)

    # ===== SLURM =====
    create_slurm(directory)

    content = "JOB -->\n"
    content += "    Write to %s\n" %str('vasp.slurm')
    print(content)

    return

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

    # Add optional argument.
    #parser.add_argument("-t", "--task", nargs='?', const='SC', help="set task type")
    #parser.add_argument("-k", "--kpoints", help="set k-points")
    #parser.add_argument("-q", "--queue", choices=('Gold', 'Test'), help="pbs queue type")
    #parser.add_argument("-n", "--ncpu", help="cpu number in total")

    parser.add_argument(
        "-d", "--cwd", default="./",
        help="current working directory"
    )
    parser.add_argument(
        "-f", "--file", help="structure file in any format"
    )
    parser.add_argument(
        "-i", "--incar",
        help="template incar file"
    )
    parser.add_argument(
        # "-c", "--copt", action='store_true',
        "-c", "--copt", type=int, nargs=2,
        help="use constrained optimisation for transition state search"
    )
    parser.add_argument(
        "-s", "--sort", action='store_false',
        help="sort atoms by elemental numbers and z-positions"
    )
    parser.add_argument(
        "--sub", action='store_true',
        help="submit the job after creating input files"
    )

    args = parser.parse_args()

    # Add all possible arguments in INCAR file.
    struct_path = Path(args.file)
    #if struct_path.suffix != '.xsd':
    #    raise ValueError('only support xsd format now...')

    cwd = Path(args.cwd)
    if cwd != Path.cwd():
        directory = cwd.resolve() / struct_path.stem
    else:
        directory = Path.cwd().resolve() / struct_path.stem

    if directory.exists():
        raise ValueError('targeted directory exists.')

    atoms = read_xsd2(struct_path)
    #atoms = read(struct_path)
    #atoms.set_constraint(FixAtoms(indices=range(len(atoms))))

    # sort atoms by symbols and z-positions especially for supercells
    if args.sort:
        numbers = atoms.numbers
        zposes = atoms.positions[:,2].tolist()
        sorted_indices = np.lexsort((zposes,numbers))
        atoms = atoms[sorted_indices]

        map_indices = dict(zip(sorted_indices, range(len(atoms))))
        copt = atoms.info.get('copt', None)
        if copt:
            new_copt = [map_indices[key] for key in atoms.info['copt']]
            atoms.info['copt'] = new_copt

    create_vasp_inputs(atoms, incar=args.incar, directory=directory)

    # submit job automatically
    if args.sub:
        command = 'sbatch vasp.slurm'
        proc = subprocess.Popen(
            command, shell=True, cwd=directory,
            stdout=subprocess.PIPE, stderr=subprocess.PIPE,
            encoding = 'utf-8'
        )
        errorcode = proc.wait(timeout=120) # 10 seconds
        if errorcode:
            raise ValueError('Error in generating random cells.')

        print(''.join(proc.stdout.readlines()))

    pass

max_stars_repo_name: hsulab/DailyScripts

max_issues_repo_name: hsulab/DailyScripts

max_forks_repo_name: hsulab/DailyScripts


text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import pycocotools.coco as coco
from pycocotools.cocoeval import COCOeval
import numpy as np
import json
import os

import torch.utils.data as data

class KITTIHP(data.Dataset):
    num_classes = 1
    num_joints = 9
    default_resolution = [384, 1280]

    mean = np.array([0.485, 0.456, 0.406], np.float32).reshape(1, 1, 3)
    std = np.array([0.229, 0.224, 0.225], np.float32).reshape(1, 1, 3)
    flip_idx = [[0, 1], [2, 3], [4, 5], [6, 7]]

    def __init__(self, opt, split,unlabel=False):
        super(KITTIHP, self).__init__()
        self.edges = [[0, 1], [0, 2], [1, 3], [2, 4],
                      [4, 6], [3, 5], [5, 6],
                      [5, 7]]

        self.acc_idxs = [1, 2, 3, 4, 5, 6, 7, 8]
        self.data_dir = os.path.join(opt.data_dir, 'kitti')
        self.img_dir= os.path.join(self.data_dir,'image')
        self.calib_dir = os.path.join(self.data_dir,'calib')
        if split == 'test':
            self.annot_path = os.path.join(
                self.data_dir, 'annotations',
                'image_info_test-dev2017.json').format(split)
        else:
            self.annot_path = os.path.join(
                self.data_dir, 'annotations',
                'kitti_{}_RTS3D.json').format(split)
        self.max_objs = opt.max_objs
        self.unlabel=unlabel
        self._data_rng = np.random.RandomState(123)
        self._eig_val = np.array([0.2141788, 0.01817699, 0.00341571],
                                 dtype=np.float32)
        self._eig_vec = np.array([
            [-0.58752847, -0.69563484, 0.41340352],
            [-0.5832747, 0.00994535, -0.81221408],
            [-0.56089297, 0.71832671, 0.41158938]
        ], dtype=np.float32)
        self.split = split
        self.opt = opt
        self.alpha_in_degree = False
        print('==> initializing kitti{} data.'.format(split))
        self.coco = coco.COCO(self.annot_path)
        image_ids = self.coco.getImgIds()

        if split == 'train':
            self.images = []
            for img_id in image_ids:
                idxs = self.coco.getAnnIds(imgIds=[img_id])
                if len(idxs) > 0:
                    self.images.append(img_id)
        else:
            self.images = image_ids
        self.num_samples = len(self.images)
        print('Loaded {} {} samples'.format(split, self.num_samples))

    def _to_float(self, x):
        return float("{:.2f}".format(x))

    def convert_eval_format(self, all_bboxes):
        # import pdb; pdb.set_trace()
        detections = []
        for image_id in all_bboxes:
            for cls_ind in all_bboxes[image_id]:
                category_id = 1
                for dets in all_bboxes[image_id][cls_ind]:
                    bbox = dets[:4]
                    bbox[2] -= bbox[0]
                    bbox[3] -= bbox[1]
                    score = dets[4]
                    bbox_out = list(map(self._to_float, bbox))
                    keypoints = np.concatenate([
                        np.array(dets[5:39], dtype=np.float32).reshape(-1, 2),
                        np.ones((17, 1), dtype=np.float32)], axis=1).reshape(51).tolist()
                    keypoints = list(map(self._to_float, keypoints))

                    detection = {
                        "image_id": int(image_id),
                        "category_id": int(category_id),
                        "bbox": bbox_out,
                        "score": float("{:.2f}".format(score)),
                        "keypoints": keypoints
                    }
                    detections.append(detection)
        return detections

    def __len__(self):
        return self.num_samples

    def save_results(self, results, save_dir):
        json.dump(self.convert_eval_format(results),
                  open('{}/results.json'.format(save_dir), 'w'))

    def run_eval(self, results, save_dir):
        # result_json = os.path.join(opt.save_dir, "results.json")
        # detections  = convert_eval_format(all_boxes)
        # json.dump(detections, open(result_json, "w"))
        self.save_results(results, save_dir)
        coco_dets = self.coco.loadRes('{}/results.json'.format(save_dir))
        coco_eval = COCOeval(self.coco, coco_dets, "keypoints")
        coco_eval.evaluate()
        coco_eval.accumulate()
        coco_eval.summarize()
        coco_eval = COCOeval(self.coco, coco_dets, "bbox")
        coco_eval.evaluate()
        coco_eval.accumulate()
        coco_eval.summarize()


text:
#! /usr/bin/env python3

import argparse
import pandas as pd
import numpy as np
import yaml
import os
from itertools import chain, zip_longest
import pickle
import json
import shutil

def plyify(r, t, ply_thickness, reverse=False):
    """Fill thickness distribution using plies."""
    active = []
    done = []
    np = len(r)
    for i in range(np):
        while t[i] > len(active) * ply_thickness:
            active.append([r[i], -1])
        while t[i] <= (len(active) - 1) * ply_thickness:
            pop = active.pop()
            pop[-1] = r[i]
            done.append(pop)
    for j in active:
        j[-1] = r[i]
        done.append(j)
    for i in done:
        i[0] = round(i[0] * 1e3, -1) * 1e-3
        i[1] = round(i[1] * 1e3, -1) * 1e-3
    # reverse the stack so that longest plies get the lowest ply numbers
    if reverse:
        return done
    else:
        return [i for i in reversed(done)]

def ply_stack(r, t, t_ply=1.0, reverse=False, subdivisions=5000, material=11):
    """
    Taking in a r,t thickness over radius distribution, split it up in
    plies, note that it rounds down in thickness (if the last ply does not fit
    in the thickness it does not go on), so if you have thick plies,
    make sure that the mm thickness is slightly up from a whole number of plies.
    """
    x = np.linspace(min(r), max(r), subdivisions)
    y = np.interp(x, list(r), t)
    st = plyify(x, y, t_ply, reverse)
    return [[material, t_ply] + i for i in st]

def coreblock(r, t, subdivisions=200, material=11):
    """Make a thickness distribution into a block divisions."""

    x = np.array(sorted(list(np.linspace(min(r), max(r), int(subdivisions))) + list(r)))
    y = np.interp(0.5 * (x[:-1] + x[1:]), list(r), t)

    if type(material) == str and material.startswith(
        "["
    ):  # deal with multi-material cores here, parse the multi-material distribution
        mat = list(zip(*eval(material)))
        mm = np.interp(np.linspace(0, 1, len(x)), mat[1], mat[0])
        mm = [min(mat[0], key=lambda x: abs(i - x)) for i in mm]
    elif type(material) == int or type(eval(material)) == int:
        mm = np.ones(len(x)) * int(material)

    stack = [i for i in list(zip(mm, y, x[:-1], x[1:])) if i[1] != 0]
    return stack

def number_stack(stack, splitstack, key, increment):
    """Create ply numbering for the plies in the stack, depends on the start key and increment as well as the above/below ratio (splitstack)."""
    sp = np.nan_to_num(splitstack.astype(float))
    if sp.sum() != 1.0:
        exit("sum %s not 1." % splitstack)

    ky = np.nan_to_num(key.astype(float)).astype(int)
    inc = np.nan_to_num(increment.astype(float), nan=1).astype(int)

    ss = (len(stack) * sp).round().astype(int)
    st, sb = np.array(range(ss[0])), np.array(range(ss[1]))
    stt = ky[0] + st * inc[0]
    sbt = ky[1] + sb * inc[1]
    seq = np.array(
        [i for i in chain.from_iterable(zip_longest(stt, sbt)) if i is not None]
    ).astype(int)
    return seq

def get_coverage(coverage, parameters, rr):
    # replace the parameter names by array
    for p in parameters:
        if coverage.find(p) != -1:
            coverage = coverage.replace(p, repr(parameters[p]))

    o = []
    for i in coverage:
        out = [
            i,
            np.ones(len(rr)) * coverage[i][0],
            np.ones(len(rr)) * coverage[i][1],
            np.ones(len(rr)) * coverage[i][2],
        ] + [rr]
        o.append(out)

    return o

def lamplan2plies(blade):
    root_radius = blade["planform"]["z"][0][1]
    tip_radius = blade["planform"]["z"][-1][1]

    slabs = blade["laminates"]["slabs"]

    allstacks = []
    parameters = {}  # parameters defined in the laminate plan csv file

    # use a multiple of lamplan length as radius grid to interpolate geometric variables to
    n_s = round(tip_radius * 4)

    rr = np.linspace(root_radius, tip_radius, n_s)

    material_map = {}

    for i in slabs:
        name = i
        material = "glass_biax" if "material" not in slabs[i] else slabs[i]["material"]
        if material not in material_map:
            material_map[material] = len(material_map) + 1
        grid = "lamplan" if "grid" not in slabs[i] else slabs[i]["grid"]
        coverage = (
            slabs[i]["cover"]
            if "cover" in slabs[i]
            else exit("no cover defined for slab %s" % i)
        )
        draping = "plies" if "draping" not in slabs[i] else slabs[i]["draping"]
        splitstack = (
            np.array([1, 0])
            if "splitstack" not in slabs[i]
            else np.array(slabs[i]["splitstack"])
        )
        key = (
            np.array([0, 4000]) if "key" not in slabs[i] else np.array(slabs[i]["key"])
        )
        increment = (
            np.array([1, -1])
            if "increment" not in slabs[i]
            else np.array(slabs[i]["increment"])
        )
        scale = (
            tip_radius - root_radius if "rscale" not in slabs[i] else slabs[i]["rscale"]
        )
        ply_thickness = (
            1.0 if "ply_thickness" not in slabs[i] else slabs[i]["ply_thickness"]
        )
        r, t = np.array(slabs[i]["slab"]).T
        r *= scale
        t *= ply_thickness

        cover = get_coverage(coverage, parameters, rr)
        if draping == "blocks":
            stack = coreblock(r, t, material=material_map[material])
        elif draping == "plies":
            stack = ply_stack(
                r, t, float(ply_thickness), material=material_map[material]
            )

        # assigns keys to the plies in the stack depending on how the stack is split up
        # what increment the ply keys are stacked with (non-1 increment allowing interleaving of plies)
        # print(splitstack, key)
        stack_numbering = number_stack(stack, splitstack, key, increment)
        allstacks.append((name.strip(), grid.strip(), cover, stack_numbering, stack))

    print("material map ", material_map)
    if "materials" in blade:
        mdb = blade["materials"]
        assert (os.path.isfile(mdb), "mdb {mdb} not a file")
        material_map["matdb"] = mdb
        # copy the material db over to the working directory, this means that rerunning
        # a model in a workdir is not affected by changes in the source matdb, this is
        # the desired behaviour
        shutil.copyfile(
            blade["materials"], os.path.join(blade["general"]["workdir"], mdb)
        )
    else:
        print("no material db defined in blade file")
    matmap = os.path.join(blade["general"]["workdir"], "material_map.json")
    json.dump(material_map, open(matmap, "w"))
    print("written material map to %s" % matmap)
    return allstacks

def main():
    parser = argparse.ArgumentParser(
        description="Split up the slab based laminate plan into plies (for laminae) and blocks (for core materials), write to a .pck file for use in draping"
    )
    parser.add_argument("yaml", help="Laminate plan yaml file")
    parser.add_argument("--out", default="__lamplan.pck", help="Output file name")
    args = parser.parse_args()

    blade = yaml.load(open(args.yaml, "r"), Loader=yaml.CLoader)

    allstacks = lamplan2plies(blade)

    of = args.out
    pickle.dump(allstacks, open(of, "wb"))

    print("written plydrape to %s" % of)

if __name__ == "__main__":
    main()


text:
import numpy as np
import pytest

import nengo
from nengo.exceptions import BuildError, ValidationError

@pytest.mark.parametrize("use_dist", (False, True))
@pytest.mark.parametrize("use_scipy", (False, True))
def test_sparse(use_dist, use_scipy, Simulator, rng, seed, plt, monkeypatch, allclose):
    if use_scipy:
        scipy_sparse = pytest.importorskip("scipy.sparse")
    else:
        monkeypatch.setattr(nengo.transforms, "scipy_sparse", None)
        monkeypatch.setattr(nengo.utils.numpy, "scipy_sparse", None)
        monkeypatch.setattr(nengo.utils.numpy, "is_spmatrix", lambda obj: False)

    input_d = 4
    output_d = 2
    shape = (output_d, input_d)

    inds = np.asarray([[0, 0], [1, 1], [0, 2], [1, 3]])
    weights = rng.uniform(0.25, 0.75, size=4)
    if use_dist:
        init = nengo.dists.Uniform(0.25, 0.75)
        indices = inds
    elif use_scipy:
        init = scipy_sparse.csr_matrix((weights, inds.T), shape=shape)
        indices = None
    else:
        init = weights
        indices = inds

    transform = nengo.transforms.Sparse(shape, indices=indices, init=init)

    sim_time = 1.0
    with nengo.Network(seed=seed) as net:
        x = nengo.processes.WhiteSignal(period=sim_time, high=10, seed=seed + 1)
        u = nengo.Node(x, size_out=4)
        a = nengo.Ensemble(100, 2)
        conn = nengo.Connection(u, a, synapse=None, transform=transform)
        ap = nengo.Probe(a, synapse=0.03)

    def run_sim():
        with Simulator(net) as sim:
            sim.run(sim_time)
        return sim

    if use_scipy:
        sim = run_sim()
    else:
        with pytest.warns(UserWarning, match="require Scipy"):
            sim = run_sim()

    actual_weights = sim.data[conn].weights

    full_transform = np.zeros(shape)
    full_transform[inds[:, 0], inds[:, 1]] = weights
    if use_dist:
        actual_weights = actual_weights.toarray()
        assert np.array_equal(actual_weights != 0, full_transform != 0)
        full_transform[:] = actual_weights

    conn.transform = full_transform
    with Simulator(net) as ref_sim:
        ref_sim.run(sim_time)

    plt.plot(ref_sim.trange(), ref_sim.data[ap], ":")
    plt.plot(sim.trange(), sim.data[ap])

    assert allclose(sim.data[ap], ref_sim.data[ap])

@pytest.mark.parametrize("encoders", (True, False))
@pytest.mark.parametrize("decoders", (True, False))
def test_sparse_nef(encoders, decoders, Simulator):
    """Sparse transforms currently don't work with NEF connections."""

    with nengo.Network() as net:
        transform = nengo.transforms.Sparse((2, 2), indices=[[0, 1], [1, 0]])
        a = nengo.Ensemble(2, 2)
        b = nengo.Ensemble(2, 2)
        nengo.Connection(
            a if decoders else a.neurons,
            b if encoders else b.neurons,
            transform=transform,
        )

    if decoders:
        # error if decoders
        with pytest.raises(BuildError, match="decoded connection"):
            with Simulator(net):
                pass
    else:
        # no error
        with Simulator(net):
            pass

def test_sparse_validation_errors():
    with pytest.raises(ValidationError, match="Either `init` must be a `scipy.sparse"):
        nengo.Sparse((1, 1))

    nengo.Sparse((3, 3), indices=[(0, 0), (1, 2)], init=[1, 2])
    with pytest.raises(ValidationError, match="Must be a vector.*length as `ind"):
        nengo.Sparse((3, 3), indices=[(0, 0), (1, 2)], init=[1, 2, 3])

def test_sparseinitparam_errors():
    class TestClass:
        sparse = nengo.transforms.SparseInitParam("sparse")

    test = TestClass()
    with pytest.raises(ValidationError, match="Must be `.*SparseMatrix` or .*spmatrix"):
        test.sparse = "a"

def test_notransform():
    transform = nengo.transforms.NoTransform(3)
    assert transform.size_in == transform.size_out == 3

    with pytest.raises(TypeError, match="Cannot sample a NoTransform"):
        transform.sample()