text:
import itertools
from tqdm import tqdm
from collections import Counter
# from arrs import *
import math
import numpy as np
from sympy import primefactors, sieve
import matplotlib.pyplot as plt
import plotly.graph_objects as go

import networkx as nx

from itertools import chain, combinations, combinations_with_replacement

def powerset(iterable):
    "powerset([1,2,3]) --> () (1,) (2,) (3,) (1,2) (1,3) (2,3) (1,2,3)"
    s = list(iterable)
    return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))

def bin_operation(a, b, mod):
    return (a*b)%mod

def check_idemp(a, mod):
    return bin_operation(a, a, mod) == a

def absorption_laws(a, b, mod1, mod2):
    return bin_operation(a, bin_operation(a, b, mod1), mod2 ) == a and bin_operation(a, bin_operation(a, b, mod2), mod1 ) == a

def distributivity(a, b, c, mod1, mod2):
    return bin_operation(a, bin_operation(b, c, mod2), mod1) == bin_operation( bin_operation(a, b, mod1), bin_operation(a, c, mod1), mod2)

def is_distributiv(setp, mod1, mod2):
    for i in combinations_with_replacement(setp, 3):
        if distributivity(*i, mod1, mod2) == 0:
            return 0
    return 1

def check_ordered(setp, mod):
    for i in combinations(setp, 2):
        if bin_operation(i[0], i[1], mod) not in i:
            return 0
    return 1

def check_ordered_double(setp, mod1, mod2):
    for i in combinations(setp, 2):

        if bin_operation(i[0], i[1], mod1) not in i:
            return 0

        if bin_operation(i[0], i[1], mod2) not in i:
            return 0

        if bin_operation(i[0], i[1], mod1) ==  bin_operation(i[0], i[1], mod2) :
            return 0

    return 1

def inters(a, b):
    return list(set(a) & set(b))

def is_a_semigroup(setp, mod):
    for i in combinations_with_replacement(setp, 2):
        if bin_operation(i[0], i[1], mod) not in setp:
            return 0
    return 1

def is_a_lattice(setp, mod1, mod2):
    for c in combinations_with_replacement(setp, 2):
        if absorption_laws(c[0], c[1], mod1, mod2) == 0:
            return 0
    return 1

def make_a_table(setp, mod):
    d = {setp[i]: i for i in range(len(setp))}

    l = len(setp)

    M = np.zeros((l, l))

    for i in range(l):
        for j in range(l):
            M[i, j] = d[bin_operation(setp[i], setp[j], mod)]

    return M

def isomorphism_two_table(M1, M2):

    if M1.shape != M2.shape:
        return 0

    if  (np.array_equal(M1, M2) or
         np.array_equal(np.transpose(M1), M2) or
         np.array_equal(M1, np.transpose(M2)) or
         np.array_equal(np.transpose(M1), np.transpose(M2)) ):

        return 1

def check_isomorphism(setp, mod1, mod2, p = 0):
    M1 = make_a_table(setp, mod1)
    M2 = make_a_table(setp, mod2)
    if p == 1:
        print(M1)
        print(M2)

    return isomorphism_two_table(M1, M2)

def characterise(setp, mod1, mod2):
    d1 = is_distributiv(setp, mod1, mod2)
    d2 = is_distributiv(setp, mod2, mod1)

    d4 = check_ordered(setp, mod1)
    d5 = check_ordered(setp, mod2)

    d7 = check_ordered_double(setp, mod1, mod2)

    d6 = check_isomorphism(setp, mod1, mod2)

    print("ordered 2: ", d7 ,'is distributiv ', d1, d2,' is  ordered ', d4, d5, ' is isom ', d6,  setp, mod1, mod2)

def make_a_matrix(setp, mod):
    l = len(setp)
    M = np.zeros((l, l))

    for i in range(l):
        for j in range(l):
            M[i, j] = ((setp[i]*setp[j])%mod == setp[i])

    return M

def find_one_nums_like_h():
    arr = []
    for i in tqdm(range(2, 4_000)):
        out = []
        for j in range(0, i):
            if check_idemp(j, i):
                out.append(j)
        t = 0

        r = (out, i)

        if len(out) > 2:
            arr.append(r)

    out = []

    for k in tqdm(combinations(arr, 2)):
        t = inters(k[0][0], k[1][0])
        if len(t) > 3:
            g = (sorted(t),  k[0][1], k[1][1] )
            out.append(g)
            #print(g)

    for p in out:

        for i in powerset(p[0]):
            if is_a_semigroup(i, p[1]) and is_a_semigroup(i, p[2]):
                if len(i) > 3 and is_a_lattice(i,  p[1],  p[2]):

                    characterise(i,  p[1],  p[2])

def find_sets_like_hilbert():
     for i in range(1, 31):
        t = 0

        #print(i)
        for s in powerset(range(0, i)):
            c = 1

            for k in combinations_with_replacement(s, 2):
                if k[0] == k[1]:
                    if (k[0] * k[0])%i != k[0]:
                        c = 0
                        break

                if (k[0] * k[1])%i not in s:
                    c = 0
                    break

            if c == 1 and s != () :
                #print([s, i], ',')
                t +=1

        print(i, t)

import matplotlib.pyplot as plt
import networkx as nx

def show_graph_with_labels(adjacency_matrix, mylabels):
    rows, cols = np.where(adjacency_matrix == 1)
    edges = zip(rows.tolist(), cols.tolist())
    gr = nx.Graph()
    gr.add_edges_from(edges)
    nx.draw(gr, node_size=500, labels=mylabels, with_labels=True)
    plt.show()

def visualize(arr, mod):
    M = make_a_matrix(arr, mod)
    print(M)
    mylabels = {i:str(arr[i]) for i in range(len(arr))}
    G = nx.from_numpy_matrix(M, create_using=nx.DiGraph)
    options = {
    'node_color': 'white',
    'node_size': 400,
    'width': 1,
    'arrowstyle': '-|>',
    'arrowsize': 22,
    }
    nx.draw_circular(G, labels=mylabels, arrows=True, **options)
    plt.draw()
    plt.show()

if __name__ == "__main__":
    #find_sets_like_hilbert()
    find_one_nums_like_h()
    #(0, 190, 855, 1065), 1330, 1995

    # visualize((0, 385, 595, 925), 1155)
    # visualize((0, 385, 595, 925), 2310)

    #check_isomorphism((0, 190, 855, 1065) ,1330, 1995, 1)

max_stars_repo_name: weightan/some-numerical-experiments

max_issues_repo_name: weightan/some-numerical-experiments

max_forks_repo_name: weightan/some-numerical-experiments

reason: import numpy,from sympy,import networkx


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

# In[1]:

import os
import random
import numpy as np

import torch
import torch.nn as nn
import torchvision

from tqdm import tqdm

from src import arch
from src import utils
from config import cfg
from src import pytorch_utils as ptu

import warnings
warnings.filterwarnings("ignore")

# assert torch.cuda.is_available(), "no CUDA"

# In[2]:

cfg.seed = 42
random.seed(cfg.seed)
torch.random.manual_seed(cfg.seed)
torch.manual_seed(cfg.seed)
torch.backends.cudnn.deterministic = True
cfg.num_workers = 0

# In[3]:

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
print(device)
if device.type == 'cuda':
    torch.backends.cudnn.benchmark = True
    print(torch.cuda.get_device_properties(device))

# In[4]:

print('clf.version\n', cfg.clf.version)

# In[5]:

cfg.clf.load = 'best'
cfg.preload_data = False
# cfg.prints = 'display'
# cfg.tqdm_bar = True

# In[6]:

if cfg.clf.load is not None and os.path.exists(os.path.join(cfg.models_dir, cfg.clf.version, ptu.naming_scheme(cfg.clf.version, epoch=cfg.clf.load)) + '.pth'):
    checkpoint = ptu.load_model(version=cfg.clf.version, models_dir=cfg.models_dir, epoch=cfg.clf.load)

# In[7]:

train_dataset = utils.Dataset(os.path.join(cfg.data_path, 'train'), cfg.clf.train_transforms, preload_data=cfg.preload_data, tqdm_bar=cfg.tqdm_bar)
val_dataset = utils.Dataset(os.path.join(cfg.data_path, 'val'), cfg.clf.val_transforms, preload_data=cfg.preload_data, tqdm_bar=cfg.tqdm_bar)

train_loader = torch.utils.data.DataLoader(train_dataset,
                                           batch_size=checkpoint.model.batch_size,
                                           num_workers=cfg.num_workers,
                                           drop_last=True, shuffle=True, pin_memory=True)

val_loader = torch.utils.data.DataLoader(val_dataset,
                                         batch_size=checkpoint.model.batch_size,
                                         num_workers=cfg.num_workers,
                                         drop_last=True, shuffle=False, pin_memory=True)

# In[8]:

train_loss, train_score, train_results = checkpoint.evaluate(train_loader,
                                                             device=device,
                                                             tqdm_bar=cfg.tqdm_bar)

# In[9]:

val_loss, val_score, val_results = checkpoint.evaluate(val_loader,
                                                       device=device,
                                                       tqdm_bar=cfg.tqdm_bar)

# In[10]:

print(f'train | loss: {train_loss:.4f} | top-1 acc: {train_score:.6f}')
print(f'val   | loss: {val_loss:.4f} | top-1 acc: {val_score:.6f}')

# In[ ]:

# In[ ]:

# In[ ]:

# In[ ]:

ext: py

reason: import numpy


text:
# src file of the pysfa class

import numpy as np
import cyipopt
import copy

from numpy             import exp, log, sqrt, pi
from numpy.linalg      import det, norm, solve
from scipy.special     import erf, erfc
from scipy.optimize    import bisect
from .bspline          import *
from sfa_utils.npufunc import log_erfc, special

# SFA Main Class
# =============================================================================
class SFA:
	# constructor
	# -------------------------------------------------------------------------
	def __init__(self, x, z, d, s, vtype='half_normal', Y=None, n=None,
			add_intercept_to_x=False,
			add_intercept_to_z=False,
			add_intercept_to_d=False):
		# add intercept
		if add_intercept_to_x: x = np.insert(x, 0, 1.0, axis=1)
		if add_intercept_to_z: z = np.insert(z, 0, 1.0, axis=1)
		if add_intercept_to_d: d = np.insert(d, 0, 1.0, axis=1)
		# pass in the parameters
		self.m = x.shape[0]
		self.k_beta = x.shape[1]
		self.k_gama = z.shape[1]
		self.k_deta = d.shape[1]
		self.k = self.k_beta + self.k_gama + self.k_deta
		#
		self.id_beta = slice(0, self.k_beta)
		self.id_gama = slice(self.k_beta, self.k_beta + self.k_gama)
		self.id_deta = slice(self.k_beta + self.k_gama, self.k)
		#
		self.vtype = vtype
		#
		self.x = x
		self.z = z
		self.d = d
		self.s = s
		self.v = s**2
		self.Y = Y
		#
		if n is not None:
			self.n = n
			self.N = np.sum(n)
			self.X = np.repeat(self.x, n, axis=0)
			self.Z = np.repeat(self.z, n, axis=0)
			self.D = np.repeat(self.d, n, axis=0)
			self.S = np.repeat(self.s, n, axis=0)
			self.V = np.repeat(self.v, n, axis=0)
		else:
			self.n = np.ones(self.m)
			self.N = self.m
			self.X = self.x.copy()
			self.Z = self.z.copy()
			self.D = self.d.copy()
			self.S = self.s.copy()
			self.V = self.v.copy()
		#
		self.check()
		#
		self.defaultOptions()

	# check the data consistancy
	# -------------------------------------------------------------------------
	def check(self):
		# check input data dimension
		assert len(self.x)==self.m, 'x: wrong number of rows.'
		assert len(self.z)==self.m, 'z: wrong number of rows.'
		assert len(self.d)==self.m, 'z: wrong number of rows.'
		assert len(self.s)==self.m, 's: wrong number of elements.'
		assert len(self.n)==self.m, 'n: wrong number of elements.'
		#
		if self.Y is not None:
			assert len(self.Y)==self.N, 'Y: wrong number of measurements.'
		#
		# check the vtype
		assert self.vtype in ['half_normal', 'exponential'], \
			'vtype: only support half_normal and exponential vtype.'
		#
		# check if all elements in Z is non-negative
		assert np.all(self.Z>=0.0), 'Z: all elements must be non-negative.'
		assert np.all(self.D>=0.0), 'D: all elements must be non-negative.'
		#
		# check if the system is identifiable
		ux, sx, vx = np.linalg.svd(self.X)
		uz, sz, vz = np.linalg.svd(self.Z)
		ud, sd, vd = np.linalg.svd(self.D)
		#
		assert sx[-1] > 1e-10, 'X: model mis-specification, is singular.'
		assert sz[-1] > 1e-10, 'Z: model mis-specification, is singular.'
		assert sd[-1] > 1e-10, 'D: model mis-specification, is singular.'
		#
		cx = sx[0]/sx[-1]
		cz = sz[0]/sz[-1]
		cd = sd[0]/sd[-1]
		#
		# print out the data summary
		print('number of studies:     ', self.m)
		print('number of measurements:', self.N)
		print('dimension of beta:     ', self.k_beta)
		print('dimension of gama:     ', self.k_gama)
		print('dimension of deta:     ', self.k_deta)
		print('cond number of X cov:  ', cx)
		print('cond number of Z cov:  ', cz)
		print('cond number of D cov:  ', cd)

	def checkUPrior(self, v, vname, vsize=None):
		if v is None: return None
		assert v.shape[0]==2, vname+': num of rows must be 2.'
		assert np.all(v[0]<=v[1]),\
			vname+': ubounds must be greater or equal than lbounds.'
		if vsize is not None:
			assert v.shape[1]==vsize, vname+': wrong num of cols.'

	def checkGPrior(self, v, vname, vsize=None):
		if v is None: return None
		assert v.shape[0]==2, vname+': num of rows must be 2.'
		assert np.all(v[1]>0.0), vname+': variance must be positive.'
		if vsize is not None:
			assert v.shape[1]==vsize, vname+': wrong num of cols.'

	# default options
	# -------------------------------------------------------------------------
	def defaultOptions(self):
		# default not using trimming
		self.use_trimming = False
		#
		# default not using bspline
		self.bspline_knots = None
		self.bspline_degree = None
		self.bspline_col_id = None
		self.bspline_uprior = None
		self.bspline_gprior = None
		self.bspline_mono = None
		self.bspline_cvcv = None
		self.bspline_l_linear = None
		self.bspline_r_linear = None
		#
		# default not having constrains
		self.constraint_matrix = None
		self.constraint_values = None
		#
		# default solution
		self.soln = None
		self.info = 'no solution, haven\'t run the solver yet.'
		#
		# add default uprior (constrains)
		self.beta_uprior = self.defaultUPrior(self.k_beta)
		self.gama_uprior = self.positiveUPrior(self.k_gama)
		self.deta_uprior = self.positiveUPrior(self.k_deta)
		self.uprior = np.hstack((
			self.beta_uprior,
			self.gama_uprior,
			self.deta_uprior
			))
		#
		self.beta_gprior = self.defaultGPrior(self.k_beta)
		self.gama_gprior = self.defaultGPrior(self.k_gama)
		self.deta_gprior = self.defaultGPrior(self.k_deta)
		self.gprior = np.hstack((
			self.beta_gprior,
			self.gama_gprior,
			self.deta_gprior
			))

	def defaultUPrior(self, col_len):
		uprior = np.empty((2,col_len))
		uprior[0] = -np.inf
		uprior[1] =  np.inf
		#
		return uprior

	def defaultGPrior(self, col_len):
		gprior = np.empty((2,col_len))
		gprior[0] = 0.0
		gprior[1] = np.inf
		#
		return gprior

	def positiveUPrior(self, col_len):
		uprior = np.empty((2,col_len))
		uprior[0] = 1e-8
		uprior[1] = np.inf
		#
		return uprior

	def negativeUPrior(self, col_len):
		uprior = np.empty((2,col_len))
		uprior[0] = -np.inf
		uprior[1] = -1e-8
		#
		return uprior

	# add uniform prior
	# -------------------------------------------------------------------------
	def addUPrior(self, beta_uprior, gama_uprior, deta_uprior):
		# check input
		self.checkUPrior(beta_uprior, 'beta_uprior', vsize=self.k_beta)
		self.checkUPrior(gama_uprior, 'gama_uprior', vsize=self.k_gama)
		self.checkUPrior(deta_uprior, 'deta_uprior', vsize=self.k_deta)
		#
		assert np.all(gama_uprior[0]>=0.0), \
			'gama_uprior: lbounds for variance must be non-negative'
		assert np.all(deta_uprior[0]>=0.0), \
			'deta_uprior: lbounds for variance must be non-negative'
		#
		self.beta_uprior = beta_uprior
		self.gama_uprior = gama_uprior
		self.deta_uprior = deta_uprior
		#
		self.uprior = np.hstack((
			self.beta_uprior,
			self.gama_uprior,
			self.deta_uprior
			))

	# add gaussian prior
	# -------------------------------------------------------------------------
	def addGPrior(self, beta_gprior, gama_gprior):
		# check input
		self.checkGPrior(beta_gprior, 'beta_gprior', vsize=self.k_beta)
		self.checkGPrior(gama_gprior, 'gama_gprior', vsize=self.k_gama)
		self.checkGPrior(deta_gprior, 'deta_gprior', vsize=self.k_deta)
		#
		self.beta_gprior = beta_gprior
		self.gama_gprior = gama_gprior
		self.deta_gprior = deta_gprior
		#
		self.gprior = np.hstack((
			self.beta_gprior,
			self.gama_gprior,
			self.deta_gprior
			))

	# add trimming
	# -------------------------------------------------------------------------
	def addTrimming(self, h):
		# check if h is valid
		assert 0<h<self.N, 'h should be positive and strictly less than N.'
		#
		self.h = h
		self.w = np.repeat(self.h/self.N, self.N)
		#
		print('trim %0.3f of the data' % (1.0 - self.h/self.N))
		#
		self.use_trimming = True

	# add bspline
	# -------------------------------------------------------------------------
	def addBSpline(self, knots, degree, col_id=1,
			l_linear=False,
			r_linear=False,
			bspline_uprior=None,
			bspline_gprior=None,
			bspline_mono=None,
			bspline_cvcv=None):
		# record the spline settings
		self.bspline_knots = knots
		self.bspline_degree = degree
		self.bspline_col_id = col_id
		self.bspline_l_linear = l_linear
		self.bspline_r_linear = r_linear
		self.bspline_uprior = bspline_uprior
		self.bspline_gprior = bspline_gprior
		self.bspline_mono = bspline_mono
		self.bspline_cvcv = bspline_cvcv
		#
		# check if there enough columns in x cov
		assert self.k_beta>=2, 'no x cov for bspline.'
		#
		# check the priors of bspline
		self.checkUPrior(bspline_uprior, 'bspline_uprior', vsize=len(knots)-1)
		self.checkUPrior(bspline_gprior, 'bspline_gprior', vsize=len(knots)-1)
		self.bspline_uprior = bspline_uprior
		self.bspline_gprior = bspline_gprior
		#
		# create design matrix
		X_bs, H_bs = dmatrix(self.X[:,col_id], knots, degree,
					l_linear=l_linear, r_linear=r_linear)
		#
		self.k_beta_bs = X_bs.shape[1]
		#
		# constrains: monotonicity
		C_bs = np.array([]).reshape(0, self.k_beta_bs)
		c_bs = np.array([]).reshape(2, 0)
		#
		if bspline_mono is not None:
			C_bs = self.diff1Mat(self.k_beta_bs)
			#
			if bspline_mono == 'increasing':
				c_bs = self.positiveUPrior(self.k_beta_bs - 1)
			if bspline_mono == 'decreasing':
				c_bs = self.negativeUPrior(self.k_beta_bs - 1)
		#
		if bspline_cvcv is not None:
			C_bs = np.vstack((C_bs, self.diff2Mat(self.k_beta_bs)))
			#
			if bspline_cvcv == 'convex':
				c_bs = np.hstack(
					(c_bs, self.positiveUPrior(self.k_beta_bs - 2)))
			if bspline_cvcv == 'concave':
				c_bs = np.hstack(
					(c_bs, self.negativeUPrior(self.k_beta_bs - 2)))
		#
		# constrains: uniform prior on bspline
		if bspline_uprior is not None:
			C_bs = np.vstack((H_bs, C_bs))
			c_bs = np.hstack((bspline_uprior, c_bs))

		# extend x and priors
		self.bsplineExtend(X_bs, H_bs, C_bs, c_bs, col_id)

	def bsplineExtend(self, X_bs, H_bs, C_bs, c_bs, col_id):
		# decide the delete cols id
		del_id = [col_id]
		if np.all(self.X[:,0]==1.0): del_id.append(0)
		#
		# extend the current X cov matrix
		# ---------------------------------------------------------------------
		self.X = self.delThenStack(self.X, X_bs, del_id)
		self.k_beta = self.X.shape[1]
		self.k = self.k_beta + self.k_gama + self.k_deta
		#
		self.id_beta = slice(0, self.k_beta)
		self.id_gama = slice(self.k_beta, self.k_beta + self.k_gama)
		self.id_deta = slice(self.k_beta + self.k_gama, self.k)
		#
		# add to the current beta uprior
		# ---------------------------------------------------------------------
		beta_uprior_bs = self.defaultUPrior(self.k_beta_bs)
		self.beta_uprior =\
			self.delThenStack(self.beta_uprior, beta_uprior_bs, del_id)
		self.uprior = np.hstack((
			self.beta_uprior,
			self.gama_uprior,
			self.deta_uprior
			))
		#
		# add to the current beta gprior
		# ---------------------------------------------------------------------
		beta_gprior_bs = self.defaultGPrior(self.k_beta_bs)
		self.beta_gprior = \
			self.delThenStack(self.beta_gprior, beta_gprior_bs, del_id)
		self.gprior = np.hstack((
			self.beta_gprior,
			self.gama_gprior,
			self.deta_gprior
			))
		#
		# add bspline prior
		# ---------------------------------------------------------------------
		self.H = np.hstack(
			(H_bs, np.zeros((H_bs.shape[0], self.k_beta - self.k_beta_bs))))
		#
		# add bspline constraints
		# ---------------------------------------------------------------------
		if C_bs.size > 0:
			C = np.hstack(
				(C_bs, np.zeros((C_bs.shape[0], self.k - self.k_beta_bs))))
			self.constraint_matrix = C
			self.constraint_values = c_bs

	def diff1Mat(self, k):
		M = np.zeros((k-1,k))
		diag_id0 = np.diag_indices(M.shape[0])
		diag_id1 = (diag_id0[0],diag_id0[1]+1)
		M[diag_id0] = -1.0
		M[diag_id1] =  1.0
		#
		return M

	def diff2Mat(self, k):
		M = np.zeros((k-2,k))
		diag_id0 = np.diag_indices(M.shape[0])
		diag_id1 = (diag_id0[0],diag_id0[1]+1)
		diag_id2 = (diag_id1[0],diag_id1[1]+1)
		M[diag_id0] =  1.0
		M[diag_id1] = -2.0
		M[diag_id2] =  1.0
		#
		return M

	def delThenStack(self, X, Y, del_id):
		X = np.delete(X, del_id, 1)
		return np.hstack((Y, X))

	# forcast data
	# -------------------------------------------------------------------------
	def forcastData(self, X, v, add_intercept_to_x=False):
		assert X.ndim == 2
		assert X.shape[0] == v.size, 'x and v should have same size.'
		if add_intercept_to_x: X = np.insert(X, 0, 1.0, axis=1)

		X_bs, H_bs = dmatrix(X[:, self.bspline_col_id],
							 self.bspline_knots,
							 self.bspline_degree,
							 l_linear=self.bspline_l_linear,
							 r_linear=self.bspline_r_linear)

		del_id = [self.bspline_col_id]
		if np.all(X[:,0]==1.0): del_id.append(0)

		X = self.delThenStack(X, X_bs, del_id)

		return X.dot(self.beta_soln) - v

	# data simulation
	# -------------------------------------------------------------------------
	def simData(self, beta_t, gama_t, deta_t):
		# check input dimension
		assert len(beta_t)==self.k_beta, 'beta_t: inconsistant with k_beta'
		assert len(gama_t)==self.k_gama, 'gama_t: inconsistant with k_gama'
		assert len(deta_t)==self.k_deta, 'deta_t: inconsistant with k_deta'
		#
		self.E = np.random.randn(self.N)*self.S
		self.EU = np.random.randn(self.N)*sqrt(self.Z.dot(gama_t))
		if self.vtype == 'half_normal':
			self.EV = np.abs(np.random.randn(self.N)*\
				sqrt(self.D.dot(deta_t)))
		if self.vtype == 'exponential':
			self.EV = np.random.exponential(size=self.N)*\
				sqrt(self.D.dot(deta_t))
		self.Y = self.X.dot(beta_t) + self.EU - self.EV + self.E

	# optimization: the maximum likelihood
	# -------------------------------------------------------------------------
	def optimizeSFA(self, print_level=0, max_iter=50):
		# create problem
		if self.constraint_matrix is None:
			handle = cyipopt.Problem(
				n=self.k,
				m=0,
				problem_obj=sfaObj(self),
				lb=self.uprior[0],
				ub=self.uprior[1]
				)
		else:
			handle = cyipopt.Problem(
				n=self.k,
				m=self.constraint_matrix.shape[0],
				problem_obj=sfaObj(self),
				lb=self.uprior[0],
				ub=self.uprior[1],
				cl=self.constraint_values[0],
				cu=self.constraint_values[1]
				)
		# add options
		handle.add_option('print_level', print_level)
		if max_iter is not None: handle.add_option('max_iter', max_iter)
		# initial point
		if self.soln is None:
			beta0 = np.linalg.solve(self.X.T.dot(self.X), self.X.T.dot(self.Y))
			gama0 = np.repeat(0.01, self.k_gama)
			deta0 = np.repeat(0.01, self.k_deta)
			x0 = np.hstack((beta0, gama0, deta0))
		else:
			x0 = self.soln
		# solver the problem
		soln, info = handle.solve(x0)
		# extract the solution
		self.soln = soln
		self.info = info
		self.beta_soln = soln[self.id_beta]
		self.gama_soln = soln[self.id_gama]
		self.deta_soln = soln[self.id_deta]

	# optimization: trimming method on the maximum likelihood
	# -------------------------------------------------------------------------
	def optimizeSFAWithTrimming(self, h, stepsize=1.0, max_iter=100, tol=1e-6,
			verbose=False):
		# add the trimming parameters into object
		self.addTrimming(h)
		#
		# initialize the iteration
		self.optimizeSFA()
		g = self.wGrad(self.soln)
		#
		# start iteration
		err = tol + 1.0
		iter_count = 0
		#
		while err >= tol:
			# proximal gradient step
			w_new = self.wProj(self.w - stepsize*g)
			# update information
			err = np.linalg.norm(w_new - self.w)/stepsize
			iter_count += 1
			#
			np.copyto(self.w, w_new)
			#
			self.optimizeSFA()
			g = self.wGrad(self.soln)
			obj = self.w.dot(g)
			#
			if verbose:
				print('iter %4d, obj %8.2e, err %8.2e' %
					(iter_count, obj, err))
			#
			if iter_count >= max_iter:
				print('trimming reach maximum number of iterations')
				break

	def wGrad(self, x):
		beta = x[self.id_beta]
		gama = x[self.id_gama]
		deta = x[self.id_deta]
		deta = np.maximum(deta, 0.0)
		# residual and all variances and stds
		r  = self.Y - self.X.dot(beta)
		#
		vu = self.Z.dot(gama)
		vv = self.D.dot(deta)
		#
		su = sqrt(vu)
		sv = sqrt(vv)
		#
		v1 = self.V + vu
		v2 = self.V + vu + vv
		#
		if self.vtype == 'half_normal':
			g = 0.5*r**2/v2 + 0.5*log(2.0*pi*v2) - \
				log_erfc(r*sv/sqrt(2.0*v1*v2))
		if self.vtype == 'exponential':
			a = (v1 + r*sv)/sqrt(2.0*v1)
			b = a/sv
			abs_a = sqrt(a*a)
			g = 0.5*r**2/v1 + log(2.0*abs_a) - special(b)
		#
		return g

	def wProj(self, w):
		a = np.min(w) - 1.0
		b = np.max(w) - 0.0
		#
		f = lambda x:\
		    np.sum(np.maximum(np.minimum(w - x, 1.0), 0.0)) - self.h
		#
		x = bisect(f, a, b)
		#
		return np.maximum(np.minimum(w - x, 1.0), 0.0)

	# post analysis
	# -------------------------------------------------------------------------
	def estimateRE(self):
		r  = self.X.dot(self.beta_soln) - self.Y
		vu = self.Z.dot(self.gama_soln)
		vv = self.D.dot(self.deta_soln)
		#
		if self.vtype == 'half_normal':
			self.v_soln = np.maximum(0.0, vv*r/(self.V + vv + vu))
		if self.vtype == 'exponential':
			self.v_soln = np.maximum(0.0, r - (self.V + vu)/np.sqrt(vv))
		#
		self.u_soln = vu*(self.v_soln - r)/(self.V + vu)

# cyipopt objective: maximum likelihood
# =============================================================================
class sfaObj:
	# constructor
	# -------------------------------------------------------------------------
	def __init__(self, sfa):
		self.sfa = sfa
		#
		if sfa.constraint_matrix is not None:
			C = sfa.constraint_matrix
			self.constraints = lambda x: C.dot(x)
			self.jacobian = lambda x: C.reshape(C.size)

	# objective
	# -------------------------------------------------------------------------
	def objective(self, x):
		sfa = self.sfa
		# likelihood
		g = sfa.wGrad(x)
		#
		# trimming
		# ---------------------------------------------------------------------
		if sfa.use_trimming:
			f = sfa.w.dot(g)
		else:
			f = np.sum(g)
		#
		# add priors
		# ---------------------------------------------------------------------
		f += np.sum( (x - sfa.gprior[0])**2/sfa.gprior[1] )
		#
		# add bspline priors
		# ---------------------------------------------------------------------
		if sfa.bspline_gprior is not None:
			f += np.sum( (mr.H.dot(beta) - mr.bspline_gprior[0])**2 /
					mr.bspline_gprior[1] )
		#
		return f

	# gradient
	# -------------------------------------------------------------------------
	def gradient(self, x, eps=1e-10):
		g = np.empty(x.size)
		z = x + 0j
		#
		for i in range(x.size):
			z[i] += 1j*eps
			f     = self.objective(z)
			g[i]  = f.imag/eps
			z[i] -= 1j*eps
		#
		return g

# Post analysis the result
# =============================================================================
def predictData(sfa, X=None, Z=None, D=None,
	add_intercept_to_x=False,
	add_intercept_to_z=False,
	add_intercept_to_d=False,
	use_bspline=False,
	include_random_effect=True,
	sample_size=1000,
	q=0.95,
	beta_sample=None):

	assert (X is None and Z is None and D is None) or \
		   (X is not None and Z is not None and D is not None), \
		   'X, Z, D: either provide all of them or none of them.'

	if X is None:
		X = sfa.X.copy()
		Z = sfa.Z.copy()
		D = sfa.D.copy()
	else:
		S = np.ones(X.size)
		sfa_pred = SFA(X, Z, D, S, vtype=sfa.vtype,
			add_intercept_to_x=add_intercept_to_x,
			add_intercept_to_z=add_intercept_to_z,
			add_intercept_to_d=add_intercept_to_d)
		if use_bspline:
			sfa_pred.addBSpline(
				sfa.bsplin_knots, sfa.bspline_degree,
				col_id=sfa.bspline_col_id,
				l_linear=sfa.bspline_l_linear,
				r_linear=sfa.bspline_r_linear,
				bspline_uprior=sfa.bspline_uprior,
				bspline_gprior=sfa.bspline_gprior,
				bspline_mono=sfa.bspline_mono,
				bspline_cvcv=sfa.bspline_cvcv)
		#
		X = sfa_pred.X.copy()
		Z = sfa_pred.Z.copy()
		D = sfa_pred.D.copy()

	if beta_sample is None:
		beta_sample = sampleBeta(sfa, sample_size)
	assert beta_sample.shape[0] == sample_size,\
		'beta_sample: length must be consistant with sample_size.'

	if include_random_effect:
		re_sample = sampleRE(sfa, Z, D, sample_size)
	else:
		re_sample = 0.0

	# sample observation
	N = X.shape[0]
	y_mean = np.zeros(N)
	y_negp = np.zeros(N)
	y_intv = np.zeros((N, 2))
	#
	y_sample = X.dot(beta_sample.T) + re_sample
	#
	for i in range(N):
		y_mean[i] = np.mean(y_sample[i])
		y_negp[i] = np.sum(y_sample[i] < 0.0)/sample_size
		y_intv[i] = qIntv(y_sample[i], q)
	#
	return beta_sample, re_sample, y_sample, y_mean, y_negp, y_intv

def qIntv(y_sample, q):
	assert 0.0<=q<=1.0, 'q: wrong percentage'
	#
	sample_size=y_sample.size
	#
	p = 0.5*(1.0 - q)
	head = int(p*sample_size)
	tail = head
	#
	y_sample_clip = np.sort(y_sample)[head:-tail]
	#
	return np.array([y_sample_clip[0], y_sample_clip[-1]])

def sampleRE(sfa, Z, D, sample_size):

	EU = np.random.randn(Z.shape[0], sample_size)*\
		sqrt(Z.dot(sfa.gama_soln)).reshape(Z.shape[0], 1)
	if sfa.vtype == 'half_normal':
		EV = np.abs(np.random.randn(D.shape[0], sample_size)*\
			sqrt(D.dot(sfa.deta_soln)).reshape(D.shape[0], 1))
	if sfa.vtype == 'exponential':
		EV = np.random.exponential(size=(D.shape[0], sample_size))*\
			sqrt(D.dot(sfa.deta_soln)).reshape(D.shape[0], 1)
	#
	re_sample = EU - EV
	#
	return re_sample

def sampleBeta(sfa, sample_size):
	beta_sample = np.zeros((sample_size, sfa.k_beta))
	sfa_copy = copy.deepcopy(sfa)

	for i in range(sample_size):
		sfa_copy.simData(sfa.beta_soln, sfa.gama_soln, sfa.deta_soln)
		sfa_copy.optimizeSFA()
		beta_sample[i] = sfa_copy.beta_soln
		print('sample beta progress %0.2f' % ((i+1.0)/sample_size), end='\r')

	return beta_sample

max_stars_repo_name: UW-AMO/StochasticFrontier

max_issues_repo_name: UW-AMO/StochasticFrontier

max_forks_repo_name: UW-AMO/StochasticFrontier

reason: import numpy,from numpy,from scipy


text:
# SPDX-FileCopyrightText: 2014-2020 Martin Roelfs
#
# SPDX-License-Identifier: MIT

from __future__ import division, print_function
import pytest
import pickle
from collections import OrderedDict

import numpy as np
from symfit import (
    Variable, Parameter, Fit, FitResults, Eq, Ge, CallableNumericalModel, Model
)
from symfit.distributions import BivariateGaussian
from symfit.core.minimizers import (
    BaseMinimizer, MINPACK, BFGS, NelderMead, ChainedMinimizer, BasinHopping
)
from symfit.core.objectives import (
    LogLikelihood, LeastSquares, VectorLeastSquares, MinimizeModel
)

def ge_constraint(a):  # Has to be in the global namespace for pickle.
    return a - 1

class TestTestResult():

    @classmethod
    def setup_class(cls):
        xdata = np.linspace(1, 10, 10)
        ydata = 3 * xdata ** 2

        cls.a = Parameter('a')
        cls.b = Parameter('b')

        x = Variable('x')
        y = Variable('y')
        model = Model({y: cls.a * x ** cls.b})

        fit = Fit(model, x=xdata, y=ydata)
        cls.fit_result = fit.execute()
        fit = Fit(model, x=xdata, y=ydata, minimizer=MINPACK)
        cls.minpack_result = fit.execute()
        fit = Fit(model, x=xdata, objective=LogLikelihood)
        cls.likelihood_result = fit.execute()
        fit = Fit(model, x=xdata, y=ydata, minimizer=[BFGS, NelderMead])
        cls.chained_result = fit.execute()

        z = Variable('z')
        constraints = [
            Eq(cls.a, cls.b),
            CallableNumericalModel.as_constraint(
                {z: ge_constraint}, connectivity_mapping={z: {cls.a}},
                constraint_type=Ge, model=model
            )
        ]
        fit = Fit(model, x=xdata, y=ydata, constraints=constraints)
        cls.constrained_result = fit.execute()
        fit = Fit(model, x=xdata, y=ydata, constraints=constraints,
                  minimizer=BasinHopping)
        cls.constrained_basinhopping_result = fit.execute()

    def test_params_type(self):
        assert isinstance(self.fit_result.params, OrderedDict)

    def test_minimizer_output_type(self):
        assert isinstance(self.fit_result.minimizer_output, dict)
        assert isinstance(self.minpack_result.minimizer_output, dict)
        assert isinstance(self.likelihood_result.minimizer_output, dict)

    def test_fitting(self):
        """
        Test if the fitting worked in the first place.
        """
        assert isinstance(self.fit_result, FitResults)
        assert self.fit_result.value(self.a) == pytest.approx(3.0)
        assert self.fit_result.value(self.b) == pytest.approx(2.0)

        assert isinstance(self.fit_result.stdev(self.a), float)
        assert isinstance(self.fit_result.stdev(self.b), float)

        assert isinstance(self.fit_result.r_squared, float)
        # by definition since there's no fuzzyness
        assert self.fit_result.r_squared == 1.0

    def test_fitting_2(self):
        np.random.seed(43)
        mean = (0.62, 0.71)  # x, y mean 0.7, 0.7
        cov = [
            [0.102**2, 0],
            [0, 0.07**2]
        ]
        data_1 = np.random.multivariate_normal(mean, cov, 10**5)
        mean = (0.33, 0.28)  # x, y mean 0.3, 0.3
        cov = [  # rho = 0.25
            [0.05 ** 2, 0.25 * 0.05 * 0.101],
            [0.25 * 0.05 * 0.101, 0.101 ** 2]
        ]
        data_2 = np.random.multivariate_normal(mean, cov, 10**5)
        data = np.vstack((data_1, data_2))

        # Insert them as y,x here as np fucks up cartesian conventions.
        ydata, xedges, yedges = np.histogram2d(data[:, 1], data[:, 0], bins=200,
                                               range=[[0.0, 1.0], [0.0, 1.0]],
                                               density=True)
        xcentres = (xedges[:-1] + xedges[1:]) / 2
        ycentres = (yedges[:-1] + yedges[1:]) / 2

        # Make a valid grid to match ydata
        xx, yy = np.meshgrid(xcentres, ycentres, sparse=False)

        x = Variable('x')
        y = Variable('y')

        x0_1 = Parameter('x0_1', value=0.6, min=0.5, max=0.7)
        sig_x_1 = Parameter('sig_x_1', value=0.1, min=0.0, max=0.2)
        y0_1 = Parameter('y0_1', value=0.7, min=0.6, max=0.8)
        sig_y_1 = Parameter('sig_y_1', value=0.05, min=0.0, max=0.2)
        rho_1 = Parameter('rho_1', value=0.0, min=-0.5, max=0.5)
        A_1 = Parameter('A_1', value=0.5, min=0.3, max=0.7)
        g_1 = A_1 * BivariateGaussian(x=x, y=y, mu_x=x0_1, mu_y=y0_1,
                                      sig_x=sig_x_1, sig_y=sig_y_1, rho=rho_1)

        x0_2 = Parameter('x0_2', value=0.3, min=0.2, max=0.4)
        sig_x_2 = Parameter('sig_x_2', value=0.05, min=0.0, max=0.2)
        y0_2 = Parameter('y0_2', value=0.3, min=0.2, max=0.4)
        sig_y_2 = Parameter('sig_y_2', value=0.1, min=0.0, max=0.2)
        rho_2 = Parameter('rho_2', value=0.26, min=0.0, max=0.8)
        A_2 = Parameter('A_2', value=0.5, min=0.3, max=0.7)
        g_2 = A_2 * BivariateGaussian(x=x, y=y, mu_x=x0_2, mu_y=y0_2,
                                      sig_x=sig_x_2, sig_y=sig_y_2, rho=rho_2)

        model = g_1 + g_2
        fit = Fit(model, xx, yy, ydata)
        fit_result = fit.execute()

        assert fit_result.r_squared > 0.95
        for param in fit.model.params:
            try:
                assert fit_result.stdev(param)**2 == pytest.approx(fit_result.variance(param))
            except AssertionError:
                assert fit_result.variance(param) <= 0.0
                assert np.isnan(fit_result.stdev(param))

        # Covariance matrix should be symmetric
        for param_1 in fit.model.params:
            for param_2 in fit.model.params:
                assert fit_result.covariance(param_1, param_2) == pytest.approx(fit_result.covariance(param_2, param_1), rel=1e-3)

    def test_minimizer_included(self):
        """"The minimizer used should be included in the results."""
        assert isinstance(self.constrained_result.minimizer, BaseMinimizer)
        assert isinstance(self.constrained_basinhopping_result.minimizer, BaseMinimizer)
        assert isinstance(self.likelihood_result.minimizer, BaseMinimizer)
        assert isinstance(self.fit_result.minimizer, BaseMinimizer)
        assert isinstance(self.chained_result.minimizer, ChainedMinimizer)
        for minimizer, cls in zip(self.chained_result.minimizer.minimizers, [BFGS, NelderMead]):
            assert isinstance(minimizer, cls)

    def test_objective_included(self):
        """"The objective used should be included in the results."""
        assert isinstance(self.fit_result.objective, LeastSquares)
        assert isinstance(self.minpack_result.objective, VectorLeastSquares)
        assert isinstance(self.likelihood_result.objective, LogLikelihood)
        assert isinstance(self.constrained_result.objective, LeastSquares)
        assert isinstance(self.constrained_basinhopping_result.objective, LeastSquares)

    def test_constraints_included(self):
        """
        Test if the constraints have been properly fed to the results object so
        we can easily print their compliance.
        """
        # For a constrained fit we expect a list of MinimizeModel objectives.
        for constrained_result in [self.constrained_result, self.constrained_basinhopping_result]:
            assert isinstance(constrained_result.constraints, list)
            for constraint in self.constrained_result.constraints:
                assert isinstance(constraint, MinimizeModel)

    def test_message_included(self):
        """Status message should be included."""
        assert isinstance(self.fit_result.status_message, str)
        assert isinstance(self.minpack_result.status_message, str)
        assert isinstance(self.likelihood_result.status_message, str)
        assert isinstance(self.constrained_result.status_message, str)
        assert isinstance(self.constrained_basinhopping_result.status_message, str)

    def test_pickle(self):
        for fit_result in [self.fit_result, self.chained_result,
                           self.constrained_basinhopping_result,
                           self.constrained_result, self.likelihood_result]:
            dumped = pickle.dumps(fit_result)
            new_result = pickle.loads(dumped)
            assert sorted(fit_result.__dict__.keys()) == sorted(new_result.__dict__.keys())
            for k, v1 in fit_result.__dict__.items():
                v2 = new_result.__dict__[k]
                if k == 'minimizer':
                    assert type(v1) == type(v2)
                elif k != 'minimizer_output':  # Ignore minimizer_output
                    if isinstance(v1, np.ndarray):
                        assert v1 == pytest.approx(v2, nan_ok=True)

    def test_gof_presence(self):
        """
        Test if the expected goodness of fit estimators are present.
        """
        assert hasattr(self.fit_result, 'objective_value')
        assert hasattr(self.fit_result, 'r_squared')
        assert hasattr(self.fit_result, 'chi_squared')
        assert not hasattr(self.fit_result, 'log_likelihood')
        assert not hasattr(self.fit_result, 'likelihood')

        assert hasattr(self.minpack_result, 'objective_value')
        assert hasattr(self.minpack_result, 'r_squared')
        assert hasattr(self.minpack_result, 'chi_squared')
        assert not hasattr(self.minpack_result, 'log_likelihood')
        assert not hasattr(self.minpack_result, 'likelihood')

        assert hasattr(self.likelihood_result, 'objective_value')
        assert not hasattr(self.likelihood_result, 'r_squared')
        assert not hasattr(self.likelihood_result, 'chi_squared')
        assert hasattr(self.likelihood_result, 'log_likelihood')
        assert hasattr(self.likelihood_result, 'likelihood')

max_stars_repo_name: pckroon/symfit

max_issues_repo_name: pckroon/symfit

max_forks_repo_name: pckroon/symfit


text:
import taichi as ti
import numpy as np

def test_simple_array():
    ti.init(arch=ti.cpu, debug=True, print_ir=True, cfg_optimization=False)
    ci13 = ti.type_factory_.get_custom_int_type(13, True)
    cu19 = ti.type_factory_.get_custom_int_type(19, False)

    x = ti.field(dtype=ci13)
    y = ti.field(dtype=cu19)

    N = 12

    ti.root.dense(ti.i, N)._bit_struct(num_bits=32).place(x, y)

    ti.get_runtime().materialize()

    @ti.kernel
    def set_val():
        for i in range(N):
            x[i] = -2**i
            y[i] = 2**i - 1

    @ti.kernel
    def verify_val():
        for i in range(N):
            assert x[i] == -2**i
            assert y[i] == 2**i - 1

    set_val()
    verify_val()

def test_custom_int_load_and_store():
    ti.init(arch=ti.cpu, debug=True, print_ir=True, cfg_optimization=False)
    ci13 = ti.type_factory_.get_custom_int_type(13, True)
    cu14 = ti.type_factory_.get_custom_int_type(14, False)
    ci5 = ti.type_factory_.get_custom_int_type(5, True)

    x = ti.field(dtype=ci13)
    y = ti.field(dtype=cu14)
    z = ti.field(dtype=ci5)

    test_case_np = np.array(
        [[2**12 - 1, 2**14 - 1, -(2**3)], [2**11 - 1, 2**13 - 1, -(2**2)],
         [0, 0, 0], [123, 4567, 8], [10, 31, 11]],
        dtype=np.int32)

    ti.root._bit_struct(num_bits=32).place(x, y, z)
    test_case = ti.Vector.field(3, dtype=ti.i32, shape=len(test_case_np))
    test_case.from_numpy(test_case_np)

    @ti.kernel
    def set_val(idx: ti.i32):
        x[None] = test_case[idx][0]
        y[None] = test_case[idx][1]
        z[None] = test_case[idx][2]

    @ti.kernel
    def verify_val(idx: ti.i32):
        assert x[None] == test_case[idx][0]
        assert y[None] == test_case[idx][1]
        assert z[None] == test_case[idx][2]

    for idx in range(len(test_case_np)):
        set_val(idx)
        verify_val(idx)

max_stars_repo_name: zhiyaluo/taichi

max_issues_repo_name: zhiyaluo/taichi

max_forks_repo_name: zhiyaluo/taichi


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

import argparse
import prody
import os
import shutil
import subprocess
import numpy
from os.path import join

GMX_PATH = '/usr/local/gromacs/bin/'

mdp_string = '''
define = -DPOSRES
integrator = {integrator}
nsteps = 1000
emtol = 1
nstlist = 1
coulombtype = Cut-off
vdwtype = Cut-off
ns_type = simple
rlist = 1.8
rcoulomb = 1.8
rvdw = 1.8
pbc = xyz
implicit_solvent = GBSA
gb_algorithm = OBC
sa_algorithm = ACE-approximation
rgbradii = 1.8
;nstxout = 1
'''

def parse_args():
    parser = argparse.ArgumentParser(description='Generate trajectory with gaussian flucutations.')
    parser.add_argument('pdb_file', metavar='INPUT_PDB_FILE', help='path to input pdb file')
    parser.add_argument('trajectory', metavar='TRAJECTORY', help='path to input trajectory')
    parser.add_argument('out_file', metavar='OUTPUT_PDB_FILE', help='path to input pdb file')
    parser.add_argument('--pos_res_k', type=float, default=1000.)
    args = parser.parse_args()

    return (args.pdb_file, args.trajectory, args.out_file, args.pos_res_k)

class Minimizer(object):
    def __init__(self, input_pdb_filename, trajectory_filename):
        self.input_pdb = self._load_pdb(input_pdb_filename)
        self.trajectory = self._load_pdb(trajectory_filename)

    def _load_pdb(self, in_file):
        protein = prody.parsePDB(in_file)
        return protein

    def _get_closest_frame(self):
        output = prody.AtomGroup('Cartesian average coordinates')
        output.setCoords( self.trajectory.getCoords() )
        output.setNames( self.trajectory.getNames() )
        output.setResnums( self.trajectory.getResnums() )
        output.setResnames( self.trajectory.getResnames() )

        ensemble = prody.PDBEnsemble(self.trajectory)
        ensemble.setCoords(self.input_pdb)
        ensemble.superpose()
        rmsds = ensemble.getRMSDs()
        min_index = numpy.argmin(rmsds)

        output.setCoords( ensemble.getCoordsets(min_index) )
        return output

    def _create_no_h_file(self, output_stream):
        # make the index file
        cmd = join(GMX_PATH, 'make_ndx')
        cmd += ' -f min_round_2.gro -o no_h.ndx'
        p1 = subprocess.Popen(cmd, shell=True, stdin=subprocess.PIPE,
                              stdout=output_stream, stderr=output_stream)
        p1.communicate('q\n')

        # run editconf
        edit_cmd = join(GMX_PATH, 'editconf')
        edit_cmd += ' -f min_round_2.gro -o no_h.gro -n no_h.ndx'
        p2 = subprocess.Popen(edit_cmd, shell=True, stdin=subprocess.PIPE,
                              stdout=output_stream, stderr=output_stream)
        p2.communicate('2\n')

    def _re_order(self, output_stream):
        # create a new index file
        lines = open('index.ndx').read().splitlines()
        header = lines[0]
        indices = []
        for line in lines[1:]:
            cols = line.split()
            for col in cols:
                indices.append( int(col) )
        resorted = [ indices.index(val)+1 for val in range( 1, max(indices)+1 ) ]
        with open('resort.ndx', 'w') as out:
            print >>out, header
            for val in resorted:
                print >>out, val

        # resort
        edit_cmd = join(GMX_PATH, 'editconf')
        edit_cmd += ' -f no_h.gro -o min.pdb -n resort.ndx'
        subprocess.check_call(edit_cmd, shell=True, stdout=output_stream,
                              stderr=output_stream)

    def run_minimization(self, posres_force_const=1000., output_stream=None):
        start = self._get_closest_frame()

        # create temp dir
        if os.path.isdir('Temp'):
            pass
        else:
            os.mkdir('Temp')
        os.chdir('Temp')

        # write the average file
        prody.writePDB('average.pdb', self.input_pdb)
        pdb_cmd = join(GMX_PATH, 'pdb2gmx')
        pdb_cmd += ' -f average.pdb -ff amber99sb-ildn -water none -n index.ndx -posrefc {} -o ref.gro -his'.format(
            posres_force_const)
        p = subprocess.Popen(pdb_cmd, shell=True, stdin=subprocess.PIPE,
                             stdout=output_stream, stderr=output_stream)
        p.communicate('0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n')
        # put it in a bigger box
        box_cmd = join(GMX_PATH, 'editconf')
        box_cmd += ' -f ref.gro -o ref_box.gro -c -box 999 999 999'
        subprocess.check_call(box_cmd, shell=True, stdout=output_stream,
                              stderr=output_stream)

        # write pdb file
        prody.writePDB('start.pdb', start)

        # pdb2gmx
        pdb_cmd = join(GMX_PATH, 'pdb2gmx')
        pdb_cmd += ' -f start.pdb -ff amber99sb-ildn -water none -n index.ndx -posrefc {} -his'.format(
            posres_force_const)
        p = subprocess.Popen(pdb_cmd, shell=True, stdin=subprocess.PIPE,
                             stdout=output_stream, stderr=output_stream)
        p.communicate('0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n0\n')

        # put it in a bigger box
        box_cmd = join(GMX_PATH, 'editconf')
        box_cmd += ' -f conf.gro -o box.gro -c -box 999 999 999'
        subprocess.check_call(box_cmd, shell=True, stdout=output_stream,
                              stderr=output_stream)

        #
        # Round 1
        #

        # write mdp file
        with open('min_round_1.mdp', 'w') as min_file:
            min_file.write( mdp_string.format(integrator='steep') )

        # run grompp
        grompp_cmd = join(GMX_PATH, 'grompp')
        grompp_cmd += ' -f min_round_1.mdp -c box.gro -p topol.top -o min_round_1 -r ref_box.gro'
        subprocess.check_call(grompp_cmd, shell=True, stdout=output_stream,
                              stderr=output_stream)

        # run mdrun
        md_cmd = join(GMX_PATH, 'mdrun')
        md_cmd += ' -deffnm min_round_1 -v -nt 1'
        subprocess.check_call(md_cmd, shell=True, stdout=output_stream,
                              stderr=output_stream)

        #
        # Round 2
        #

        # write mdp file
        with open('min_round_2.mdp', 'w') as min_file:
            min_file.write( mdp_string.format(integrator='l-bfgs') )

        # run grompp
        grompp_cmd = join(GMX_PATH, 'grompp')
        grompp_cmd += ' -f min_round_2.mdp -c min_round_1.gro -p topol.top -o min_round_2 -maxwarn 1 -r ref_box.gro'
        subprocess.check_call(grompp_cmd, shell=True, stdout=output_stream,
                              stderr=output_stream)

        # run mdrun
        md_cmd = join(GMX_PATH, 'mdrun')
        md_cmd += ' -deffnm min_round_2 -v -nt 1'
        subprocess.check_call(md_cmd, shell=True, stdout=output_stream,
                              stderr=output_stream)

        #
        # gather results
        #
        self._create_no_h_file(output_stream)
        self._re_order(output_stream)

        # load the pdb
        protein = prody.parsePDB('min.pdb').select('not hydrogen')

        # clean up
        os.chdir('..')
        shutil.rmtree('Temp')

        return protein

def main():
    r = parse_args()
    input_pdb_filename = r[0]
    trajectory_filename = r[1]
    output_pdb_filename = r[2]
    force_const = r[3]
    m = Minimizer(input_pdb_filename, trajectory_filename)
    minimized_protein = m.run_minimization(force_const)
    prody.writePDB(output_pdb_filename, minimized_protein)

if __name__ == '__main__':
    main()

max_stars_repo_name: grollins/netavg-django

max_issues_repo_name: grollins/netavg-django

max_forks_repo_name: grollins/netavg-django


text:
#!/usr/common/software/python/2.7-anaconda/bin/python
import os
import scipy
from scipy import array, linalg, dot
from csld.phonon import head
from numpy import *
import numpy as np
from numpy.linalg import *
import copy
import numpy.linalg
from get_matcov import get_matcov
import multiprocessing as mp
import time

def mkdir(dir):
    if os.path.isdir(dir):
        os.system(str("rm"+" -r "+dir))
        os.system(str("mkdir"+" "+dir))
    else:
        os.system(str("mkdir"+" "+dir))

def get_ifc(file):
    fin=None
    fin=open(file,'r')
    fifc=fin.readlines()
    fin.close()
    natom=int(fifc[0])
    ifc=array([[0.0 for i in range(3*natom)] for j in range(3*natom)])
    for i in range(natom):
        for j in range(natom):
            index=i*natom*4+j*4+1
            ifc[3*i+0, 3*j:3*j+3]=list(map(float, fifc[index+1].split()))
            ifc[3*i+1, 3*j:3*j+3]=list(map(float, fifc[index+2].split()))
            ifc[3*i+2, 3*j:3*j+3]=list(map(float, fifc[index+3].split()))
    return ifc

def get_dig(poscar,ifc):
    dyn=copy.deepcopy(ifc)
    #evc=copy.deepcopy(ifc)
    natom=poscar['natom']
    for i in range(natom):
        for j in range(natom):
            dyn[3*i:3*i+3, 3*j:3*j+3] = dyn[3*i:3*i+3, 3*j:3*j+3]/sqrt(poscar['mas'][i]*poscar['mas'][j])
    dyn=(dyn+transpose(dyn))/2.0
    #scipy.linalg.eig(dyn, evc, False, True, False,False,True)
    solution=numpy.linalg.eigh(dyn)
    return solution

def read_poscar(file,massall): #assuming SPOSCAR
    print(massall)
    fin=None
    fin=open(file, 'r')
    tmp=fin.readlines() #ftmp
    fin.close()
    poscar={}
    poscar['mas']=[] #done
    poscar['cor']=[] #done
    poscar['spe']=[]
    poscar['prefix']=tmp[0]
    poscar['scaling']=float(tmp[1])
    poscar['vec1']=list(map(float, tmp[2].split()))
    poscar['vec2']=list(map(float, tmp[3].split()))
    poscar['vec3']=list(map(float, tmp[4].split()))
    poscar['species']=tmp[5].split()
    poscar['specount']=list(map(int, tmp[6].split()))
    poscar['natom']=sum(poscar['specount'])
    poscar['coordinate']=tmp[7]
    for ia in range(poscar['natom']):
        poscar['cor'].append(array( list(map(float, (tmp[8+ia].split())[0:3] )) ))
    for ise in range(len(poscar['specount'])):
        iac=poscar['specount'][ise]
        for i in range(iac):
            poscar['mas'].append(massall[ise])
            poscar['spe'].append(poscar['species'][ise])
    return poscar

def write_poscar_cart(file, poscar, carposdis):
    #print("Writing to "+str(file))
    dataout="System name \n"
    dataout+="1.0 \n"
    dataout+=" ".join(map(str,poscar['vec1']))+"\n"
    dataout+=" ".join(map(str,poscar['vec2']))+"\n"
    dataout+=" ".join(map(str,poscar['vec3']))+"\n"
    dataout+=" ".join(map(str,poscar['species']))+"\n"
    dataout+=" ".join(map(str,poscar['specount']))+"\n"
    dataout+="C"+"\n"
    for pos in carposdis:
        dataout+=" ".join(map(str, pos))+"\n"
    ffout=None
    ffout=open(file, 'w')
    ffout.write(dataout)
    ffout.close()

def write_poscar_direct(file, poscar, carposdis):
    #print("Writing to "+str(file))
    dataout="System name \n"
    dataout+="1.0 \n"
    dataout+=" ".join(map(str,poscar['vec1']))+"\n"
    dataout+=" ".join(map(str,poscar['vec2']))+"\n"
    dataout+=" ".join(map(str,poscar['vec3']))+"\n"
    #PRN uses no species
    dataout+=" ".join(map(str,poscar['species']))+"\n"
    dataout+=" ".join(map(str,poscar['specount']))+"\n"
    dataout+="Direct"+"\n"
    for pos in carposdis:
        dataout+=" ".join(map(str, pos))+"\n"
    ffout=None
    ffout=open(file, 'w')
    ffout.write(dataout)
    ffout.close()

def write_zero_forces(file, poscar, carposdis):
    #print("Writing to "+str(file))
    dataout=""
    #dataout="System name \n"
    #dataout+="1.0 \n"
    #dataout+=" ".join(map(str,poscar['vec1']))+"\n"
    #dataout+=" ".join(map(str,poscar['vec2']))+"\n"
    #dataout+=" ".join(map(str,poscar['vec3']))+"\n"
    #PRN uses no species
    #dataout+=" ".join(map(str,poscar['species']))+"\n"
    #dataout+=" ".join(map(str,poscar['specount']))+"\n"
    #dataout+="Direct"+"\n"
    for pos in carposdis:
        dataout+="0.0  0.0  0.0\n"
    ffout=None
    ffout=open(file, 'w')
    ffout.write(dataout)
    ffout.close()

def dispWrap(dispAll,Lmatcov,natom,poscar,path,iconf):
    #dispAll = random.normal(mu, sigma, ndim)
    disp=(dot(Lmatcov, dispAll[iconf])).reshape((natom, 3))
    #print("---------------------------------------------------")
    #print(norm(disp)/sqrt(len(disp)))
    latvec=array([poscar['vec1'],poscar['vec2'],poscar['vec3']])
    dirpos=poscar['cor']+dot(disp, inv(latvec))
    mkdir(path+"disp-"+str(iconf+1))
    write_poscar_direct(path+"disp-"+str(iconf+1)+"/POSCAR", poscar, dirpos)
#    write_poscar_direct(path+"disp-"+str(iconf+1)+"/POSCAR.direct", poscar, dirpos)
#    write_zero_forces(path+"disp-"+str(iconf+1)+"/force.txt", poscar, dirpos)

#-----------------------------------------------------------------------
def get_qcv(masslist,temperature,path):
    readPOSCAR=0
    if readPOSCAR > 1:
        tmp=get_yaml("./mesh.yaml")
        poscar=tmp[0]
        matevc=tmp[-1]
        print("******Check evc normalization: "+str(dot(matevc[:,0], matevc[:,0])))
    else:
        poscar=read_poscar(path+"SPOSCAR",masslist)  # mass for each species

    latvec=array([poscar['vec1'],poscar['vec2'],poscar['vec3']]) * poscar['scaling']
    natom=poscar['natom']
    masses=poscar['mas']
    ndim=3*natom
    print("Number of atoms: "+str(natom))
    ifc=get_ifc(path+"FORCE_CONSTANTS_2ND")
    [eig, evc]=get_dig(poscar,ifc)
    print("Frequency-------------------------------------------------------")
    print("******Number of eigenfrequencies: "+str(len(eig)))
    freq=[]
    unit2thz=98.1761/2/pi # eV / A^2 / au to THz
    unit2mev=unit2thz*4.135665538536  #THZ to meV
    unit=unit2mev
    unit=unit2thz
    isImag=0
    for tmp in eig:
        if tmp < 0:
            freq.append(-1*sqrt(-tmp)*unit)
            if sqrt(-tmp)*unit > 0.005: # check if it the mode is imaginary
                isImag=isImag+1
        else:
            freq.append(sqrt(tmp)*unit)
#    head.writenumber(isImag, path+"isImag")
    idx = array(freq).argsort()[::1]
    sorteig=array(freq)[idx]
    sortevc=evc[:,idx]
    matcov=[[ 0.0 for i in range(ndim)] for j in range(ndim)]
    free_energy, matcov = get_matcov.pmatcov(temperature, natom, masses, path)
#    matcov=head.read2dmat(path+"matcov.dat")
#    print("Eigenvalues of covariance matrix: ")
#    print(eigvalsh(matcov))

    # Get lower-triangle matrix
    try: # Cholesky, but could fail
        Lmatcov = scipy.linalg.cholesky(matcov,lower=True)
    except: # Cholesky with added multiple of identity (Nocedal & Wright, p.51)
        if np.min(matcov) > 0 :
            tau = 0
        else:
            beta = 0.1
            tau = -np.min(matcov)+beta
        while True:
            matcov = matcov+tau*np.eye(ndim)
            try:
                Lmatcov = scipy.linalg.cholesky(matcov,lower=True)
                break
            except:
                tau = max(2*tau,beta)
#    head.write2dmat(Lmatcov.tolist(), path+"Lmatcov")

    return free_energy, Lmatcov, poscar

# ----------------------
def qcv_displace(Lmatcov,poscar,nconfig,numprocess,path):

#    print("Lattice vector: ")
#    print(latvec)
    #random seed to ensure the randoming sampling are same
#    random.seed(90001)
#    random.seed(100901)
    # pre-generate all random distributions
    mu = 0.0
    sigma = 1.0 # 0.5
    ndim = Lmatcov.shape[0]
    natom = int(ndim/3)
    dispAll=[random.normal(mu, sigma, ndim) for iconf in range(nconfig)] # Gaussian random sampling
    stime=time.time()
    if numprocess == 1:
        results=[ dispWrap(dispAll,Lmatcov,natom,poscar,path,iconf) for iconf in range(nconfig) ]
    elif numprocess > 1:
        pool = mp.Pool(processes=numprocess)
        results=[ pool.apply_async(dispWrap, args=(dispAll,Lmatcov,natom,poscar,path,iconf)) for iconf in range(nconfig) ]
        pool.close()
        print("Pool closed")
        pool.join()
    else:
        raise ValueError("numprocess should be natural number")
    print("Time for generating displacements: "+str(time.time()-stime)+"\n")

    return

#make sure path exist
#--here is path=./pbte-images/
#prepare FORCE_CONSTANT_CSLD, FORCE_CONSTANT_CSLD_OLD and POSCAR-org under path
#--FORCE_CONSTANT_CSLD and FORCE_CONSTANT_CSLD_OLD could be the same, if not, phonon frequency is the average one
#--FORCE_CONSTANT_CSLD and FORCE_CONSTANT_CSLD_OLD are all in Phonopy format
#--POSCAR-org is the supercell in line with FORCE_CONSTANT_CSLD

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


text:
from sympy import sin, cos, exp, E, series
from sympy.abc import x

def test_sin():
    e1 = sin(x).series(x, 0)
    e2 = series(sin(x), x, 0)
    assert e1 == e2

def test_cos():
    e1 = cos(x).series(x, 0)
    e2 = series(cos(x), x, 0)
    assert e1 == e2

def test_exp():
    e1 = exp(x).series(x, 0)
    e2 = series(exp(x), x, 0)
    assert e1 == e2

def test_exp2():
    e1 = exp(cos(x)).series(x, 0)
    e2 = series(exp(cos(x)),x,0)
    assert e1 == e2

max_stars_repo_name: matthew-brett/sympy

max_issues_repo_name: matthew-brett/sympy

max_forks_repo_name: matthew-brett/sympy

reason: from sympy


text:
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.tsa.api import acf

wave = pd.read_csv("data/wave.dat")

print(acf(wave["waveht"],fft=False))

plt.figure()
wave["waveht"].plot()
plt.figure()
wave["waveht"][0:60].plot()

plot_acf(wave["waveht"])
plt.show()

reason: from statsmodels


text:
from torch.utils.data import Dataset, DataLoader
from torch.utils.data.dataloader import default_collate
from collections import deque
from dataset.make_dataset import get_dataset
import torch
import numpy as np

import random

def bfs(start, adj, distance):
    """

    :param start:
    :param adj: a sparse adjacent matrix
    :param distance:
    :return:
    """
    num_nodes = adj.shape[0]
    visited = [False, ] * num_nodes
    q = deque()
    q.append((start, 0))
    visited[start] = True
    node_list = [start, ]
    while(len(q) != 0):
        cur_node, cur_dist = q.pop()
        node_list.append(cur_node)
        if(cur_dist + 1 > distance):
            break
        for next_node in adj[cur_node].nonzero()[1]:
            if not visited[next_node]:
                q.append((next_node, cur_dist + 1))
                visited[next_node] = True

    while(len(q) != 0):
        node, dist = q.pop()
        node_list.append(node)

    return node_list

def bfs_sample(start, adj, distance, sample_num):
    """

    :param start:
    :param adj:
    :param distance:
    :param sample_numbers: should be a list specific number of support node sampled at each distance
    :return:
    """
    assert distance == len(sample_num), "BFS distance should equal to length of sample_nums."
    num_nodes = adj.shape[0]
    visited = [False, ] * num_nodes
    nodelist = [start, ]
    curlist = [start, ]
    visited[start] = True
    for i in range(distance):
        nextlist = []
        for cur_node in curlist:
            downstream = []
            next_nodes = adj[cur_node].nonzero()[1]
            for node in next_nodes:
                if not visited[node]:
                    downstream.append(node)
            if len(downstream) > sample_num[i]:
                random.shuffle(downstream)
                downstream = downstream[:sample_num[i]]

            for node in downstream:
                visited[node] = True
            nextlist.extend(downstream)

        nodelist.extend(nextlist)
        curlist = nextlist

    return nodelist

class SubGraph(Dataset):
    def __init__(self, adj, features, labels, idx, num_samples):
        """

        :param adj: suppose adj a sparse adjacent matrix
        :param features: a numpy array in shape (num_nodes, num_features)
        :param labels: a numpy array in shape (num_nodes, 1) if not multi label task.
        """
        self.adj = adj
        self.features = features
        self.num_samples = num_samples
        self.num_layers = len(num_samples)
        self.labels = labels
        self.idx = idx
        self.num_nodes = len(idx)

    def __getitem__(self, item):
        node = self.idx[item]
        nodelist = bfs_sample(node, self.adj, self.num_layers, self.num_samples)
        min_adj = self.adj[nodelist][:, nodelist]
        min_features = self.features[nodelist]
        min_label = self.labels[node:node + 1]

        return (min_adj, min_features), min_label

    def __len__(self):
        return self.num_nodes

def custom_collate(batch):
    """

    :param batch: a list of tuple ((adj, features), label)
    :return:
    """
    max_nodes = max(map(lambda x: x[0][0].shape[0], batch))
    batch = [align(data, max_nodes) for data in batch]
    return default_collate(batch)

def align(data_tuple, max_nodes):
    """
    consider adj as sparse matrix and features as dense ndarray.
    :param data_tuple:
    :param max_nodes:
    :return:
    """
    (adj, features), label = data_tuple
    adj.resize(max_nodes, max_nodes)
    adj_fill = adj.toarray()
    features_fill = np.zeros((max_nodes, features.shape[1]), )
    features_fill[:features.shape[0]] = features

    # convert numpy/scipy to torch tensor
    adj_fill = torch.FloatTensor(adj_fill)
    features_fill = torch.FloatTensor(features_fill)
    label = torch.LongTensor(label)
    return (adj_fill, features_fill), label

if __name__ == "__main__":
    # data_name = "citeseer"
    # adj, features, labels = get_dataset(data_name, './data/npz/{}.npz'.format(data_name), standardize=True,
    #                                     train_examples_per_class=40, val_examples_per_class=100)
    # features = features.toarray()
    from utils import load_reddit
    adj, features, labels, idx_train, idx_val, idx_test = load_reddit()

    # %%
    subgraph = SubGraph(adj, features, labels, np.arange(adj.shape[0]), 2, [25, 10])
    dataloader = DataLoader(subgraph, batch_size = 2, num_workers= 4, collate_fn= custom_collate)
    # %%
    for i, data in enumerate(dataloader):
        (min_adj, min_feat), label = data
        print(min_adj.shape, min_feat.shape, label.shape)
        if i > 20:
            break

max_stars_repo_name: Yindong-Zhang/myGAT

max_issues_repo_name: Yindong-Zhang/myGAT

max_forks_repo_name: Yindong-Zhang/myGAT


text:
"""Collection of tests for unified reduction functions."""

# global
import numpy as np
from hypothesis import given, strategies as st

# local
import ivy
import ivy.functional.backends.numpy as ivy_np
import ivy_tests.test_ivy.helpers as helpers

# random_uniform
@given(
    data=st.data(),
    shape=helpers.get_shape(),
    dtype=st.sampled_from(ivy_np.valid_float_dtypes),
    as_variable=st.booleans(),
)
def test_random_uniform(data, shape, dtype, as_variable, device, call):
    low, high = data.draw(helpers.get_bounds(dtype))
    # smoke test
    if as_variable and call is helpers.mx_call:
        # mxnet does not support 0-dimensional variables
        return
    low_tnsr = ivy.array(low, dtype=dtype, device=device) if low is not None else None
    high_tnsr = (
        ivy.array(high, dtype=dtype, device=device) if high is not None else None
    )
    if as_variable and (low is not None):
        low_tnsr = ivy.variable(low_tnsr)
    if as_variable and (high is not None):
        high_tnsr = ivy.variable(high_tnsr)
    kwargs = {
        k: v for k, v in zip(["low", "high"], [low_tnsr, high_tnsr]) if v is not None
    }
    if shape is not None:
        kwargs["shape"] = shape
    ret = ivy.random_uniform(**kwargs, device=device)
    # type test
    assert ivy.is_ivy_array(ret)
    # cardinality test
    if shape is None:
        assert ret.shape == ()
    else:
        assert ret.shape == shape
    # value test
    ret_np = call(ivy.random_uniform, **kwargs, device=device)
    assert np.min((ret_np < (high if high else 1.0)).astype(np.int32)) == 1
    assert np.min((ret_np >= (low if low else 0.0)).astype(np.int32)) == 1

# random_normal
@given(
    data=st.data(),
    shape=helpers.get_shape(),
    dtype=st.sampled_from(ivy_np.valid_float_dtypes),
    as_variable=st.booleans(),
)
def test_random_normal(data, shape, dtype, as_variable, device, call):
    mean, std = data.draw(helpers.get_mean_std(dtype))
    ivy.seed(0)
    # smoke test
    if as_variable and call is helpers.mx_call:
        # mxnet does not support 0-dimensional variables
        return
    mean_tnsr = (
        ivy.array(mean, dtype=dtype, device=device) if mean is not None else None
    )
    std_tnsr = ivy.array(std, dtype=dtype, device=device) if std is not None else None
    if as_variable and (mean is not None):
        mean_tnsr = ivy.variable(mean_tnsr)
    if as_variable and (std is not None):
        std_tnsr = ivy.variable(std_tnsr)
    kwargs = {
        k: v for k, v in zip(["mean", "std"], [mean_tnsr, std_tnsr]) if v is not None
    }
    if shape is not None:
        kwargs["shape"] = shape
    ret = ivy.random_normal(**kwargs, device=device)
    # type test
    assert ivy.is_ivy_array(ret)
    # cardinality test
    if shape is None:
        assert ret.shape == ()
    else:
        assert ret.shape == shape

# multinomial
@given(
    data=st.data(),
    num_samples=st.integers(1, 2),
    replace=st.booleans(),
    dtype=st.sampled_from(ivy_np.valid_float_dtypes),
    tensor_fn=st.sampled_from([ivy.array]),
)
def test_multinomial(data, num_samples, replace, dtype, tensor_fn, device, call):
    probs, population_size = data.draw(helpers.get_probs(dtype))
    if (
        call in [helpers.mx_call, helpers.tf_call, helpers.tf_graph_call]
        and not replace
        or dtype == "float64"
    ):
        # mxnet and tenosorflow do not support multinomial without replacement
        return
    # smoke test
    probs = tensor_fn(probs, dtype=dtype, device=device) if probs is not None else probs
    batch_size = probs.shape[0] if probs is not None else 2
    ret = ivy.multinomial(population_size, num_samples, batch_size, probs, replace)
    # type test
    assert ivy.is_ivy_array(ret)
    # cardinality test
    assert ret.shape == tuple([batch_size] + [num_samples])

# randint
@given(
    data=st.data(),
    shape=helpers.get_shape(allow_none=False),
    dtype=st.sampled_from(ivy_np.valid_int_dtypes),
    as_variable=st.booleans(),
)
def test_randint(data, shape, dtype, as_variable, device, call):
    # smoke test
    low, high = data.draw(helpers.get_bounds(dtype))
    if (
        call in [helpers.mx_call, helpers.torch_call]
        and as_variable
        or dtype == "uint64"
        or call == helpers.torch_call
        and dtype[0] == "u"
    ):
        # PyTorch and MXNet do not support non-float variables
        return
    low_tnsr = ivy.array(low, dtype=dtype, device=device)
    high_tnsr = ivy.array(high, dtype=dtype, device=device)
    if as_variable:
        low_tnsr, high_tnsr = ivy.variable(low_tnsr), ivy.variable(high_tnsr)
    kwargs = {
        k: v for k, v in zip(["low", "high"], [low_tnsr, high_tnsr]) if v is not None
    }
    kwargs["shape"] = shape
    ret = ivy.randint(**kwargs, device=device)
    # type test
    assert ivy.is_ivy_array(ret)
    # cardinality test
    assert ret.shape == shape
    # value test
    ret_np = call(ivy.randint, **kwargs, device=device)
    assert np.min((ret_np < high).astype(np.int64)) == 1
    assert np.min((ret_np >= low).astype(np.int64)) == 1

# seed
@given(
    seed_val=st.integers(min_value=0, max_value=2147483647),
)
def test_seed(seed_val):
    # smoke test
    seed_val = ivy.array(seed_val)
    ivy.seed(seed_val)

# shuffle
@given(
    data=st.data(),
    dtype=st.sampled_from(ivy_np.valid_float_dtypes),
    as_variable=st.booleans(),
)
def test_shuffle(data, dtype, as_variable, device, call):
    # smoke test
    shape = data.draw(helpers.get_shape(min_num_dims=1))
    x = data.draw(helpers.array_values(dtype, shape))
    x = ivy.array(x, dtype=dtype, device=device)
    if as_variable:
        x = ivy.variable(x)
    ret = ivy.shuffle(x)
    # type test
    assert ivy.is_ivy_array(ret)
    # cardinality test
    assert ret.shape == x.shape
    # value test
    ivy.seed(0)
    first_shuffle = call(ivy.shuffle, x)
    ivy.seed(0)
    second_shuffle = call(ivy.shuffle, x)
    assert np.array_equal(first_shuffle, second_shuffle)


text:
# Copyright 2020 Ian Rankin
#
# Permission is hereby granted, free of charge, to any person obtaining a copy of this
# software and associated documentation files (the "Software"), to deal in the Software
# without restriction, including without limitation the rights to use, copy, modify, merge,
# publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons
# to whom the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all copies or
# substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
# INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE
# FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
# OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.
#
## @package HSignature.py
# Written Ian Rankin - January 2020
#
# A basic structure to handle HSignatures
# Each HSignature is stored as the set of each non-signature obstacles
# or some partial or complete list.

import numpy as np
import copy
import sys

# for checking python version (required for hashing function)
import sys

class HSignature(object):

    ## @overide
    # @param edge - the HEdge that should be added to the H signature to update
    #               to a new state
    # @return - True if valid edge crossing, false if the crossing is invalid (loop)
    def edge_cross(self, edge):
        raise NotImplementedError()

    ## cross
    # A function to add a crossing to the HSignature
    # @param id - the id of the feature
    # @param value - the sign of the crossing (+1, 0, -1) 0, makes no sense to be given
    def cross(self, id, value):
        raise NotImplementedError()

    ## compute_line_segment
    # This function turns the current HSignature into the h signature for a
    # line-segment
    # @param pt_a - the first point of the line segment (numpy)
    # @param pt_a - the second point of the line segment (numpy)
    def compute_line_segment(self, pt_a, pt_b, features, ray_angle=np.pi/2):
        raise NotImplementedError()

    ## ensure, correct copying of the signatures
    def copy(self):
        raise NotImplementedError()

    ## concatination
    # Note concationation for HSignatures may not be commutative (a+b)!=(b+a)
    # For Homotopy signatures (a+b)!=(b+a)
    # For Homology signatures a+b=b+a
    def __add__(self, other):
        raise NotImplementedError()

    ## inverse concatination
    def __sub__(self, other):
        raise NotImplementedError()

    def __hash__(self):
        raise NotImplementedError()

#