text:
import random
import math
from .tools import rescale

# moved: between, odds, pick to here

def between(a, b):
    """
    Returns a random value between a and b (exclusive).

    An int is returned if both a and b are ints,
    otherwise a float is returned. The a value can be
    less than, equal to, or greater than b.

    Parameters
    ----------
    a : int | float
        The lower bound for the random selection.
    b : int | float
        The upper bound for the random selection.

    Returns
    -------
        A randomly chosen value between a and b.
    """
    if (a < b):
        if isinstance(a, int) and isinstance(b, int):
            return random.randrange(a, b)
        else:
            return a + (random.random() * (b - a))
    if (a > b):
        if isinstance(a, int) and isinstance(b, int):
            return a - random.randrange(b)
        else:
            return a - (random.random() * (a - b))
    return a

def odds(prob, true=True, false=False):
    """
    Returns true value if a random choice is less than prob else returns
    the false.

    Parameters
    ----------
    prob : float | int
        A value between 0 and 1 inclusive.
    true : value
        If odds is true the true value is returned.
    false : value
        If odds is false the false value is returned.

    Returns
    -------
    Either the true value or the false value based on probability.
    """
    if 0 <= prob <= 1.0:
        return true if random.random() < prob else false
    raise ValueError(f"{prob} not a number between 0 and 1.")

def pick(*args):
    """
    Randomly picks one of its arguments to return.

    Parameters
    ----------
    *args : variadic
        A series of values from which one will be randomly selected.

    Returns
    -------
    One value from the argument list.
    """
    return args[random.randrange(len(args))]

def vary(value, variance, shift=None):
    """
    Returns a random number that deviates from value or list of the same by
    up to variance (1=100%) according to shift.

    Parameters
    ----------
    value : number | list
        The value or list of values to vary.
    variance : number
        The maximum amount to vary the value by, e.g. 0.5 means the variance
        could be up to 50% of the input value.
    shift : None | "+" | "-"
        Specifies where the variance is placed. If None then then value is at the
        center of what could be returned. Shift "+" places value at the minimum
        of what could be returned and "-" means that value is the maximum possible
        value returned.
    """
    if isinstance(value, list):
        return [vary(v, variance, shift) for v in value]
    if value == 0 or variance == 0:
        return val
    r = abs(value * variance)
    v = random.random() * r
    if shift == None:
        return (value - (r * .5)) + v
    if shift == '+':
        return value + v
    if shift == '-':
        return value - v
    raise ValueError("shift should be None, '+', or '-' ({shift})")

def ransegs(num, mapto, rangen=None):
    """
    Returns a list of random number segments.

    Parameters
    ----------
    num : int
        The number of random values returned in the list.
    mapto : number | [low, high]
        If mapto is a number then the values retuned will sum to that number.
        If mapto is a list then the random values are sorted and scaled to
        lie proportionatly from low to high.
    rangen : None | function
        If function is provided it will be called to return random numbers,
        otherwise uniform random numbers between are generated.
    """
    #print(f"ransegs(num={num}, mapto={mapto})")
    segs = []
    rsum = 0
    if rangen == None:
        rangen = random.random
    for i in range(num):
        n = rangen()
        segs.append(n)
        rsum += n
    if isinstance(mapto, (list, tuple)):
        segs.sort()
        #print("raws", segs)
        for i in range(num):
            segs[i] = rescale(segs[i], segs[0], segs[-1], mapto[0], mapto[1])
    else:
        #print("raws", segs)
        for i in range(num):
            segs[i] = rescale(segs[i], 0, rsum, 0, mapto)
    #print("segs", segs)
    return segs

#-----------------------------------------------------------------------------
# Random number functions with different distributions

def uniran():
    """
    Returns a uniformly chosen floating point value between 0.0 and 1.0.
    See also: lowran(), midran(), highran().
    """
    return random.random()

def lowran():
    """
    Returns a floating point value between 0.0 and 1.0 with lower
    values more likely. See also: uniran(), midran(), highran().
    """
    return min(random.random(), random.random())

def midran():
    """
    Returns a floating point value between 0.0 and 1.0 with midrange
    values more likely. See also: uniran(), lowran(), highran().
    """
    return (random.random() + random.random()) / 2.0

def highran():
    """
    Returns a floating point value between 0.0 and 1.0 with higher
    values more likely. See also: uniran(), lowran(), midran().
    """
    return max(random.random(), random.random())

def beta(alpha=.5, beta=.5):
    """
    Returns value between 0 and 1 from the beta distribution (https://en.wikipedia.org/wiki/Beta_distribution).
    When a=b=1 the distribution is uniform. When alpha=beta, the distribution is
    symmetric around .5. When alpha<1 and beta<1 then the density of larger and
    smaller numbers increases. When alpha>1 and beta>1, density is similar to the
    gaussian distribution.
    """
    ra = 1.0 / alpha
    rb = 1.0 / beta
    while True:
        r1 = random.random()
        r2 = random.random()
        y1 = r1 ** ra
        y2 = r2 ** rb
        y3 = y1 + y2
        if y3 <= 1.0:
            return (y1 / y3)

def gamma(alpha=1):
    """
    Returns value greater than 0 from the gamma distribution
    (https://en.wikipedia.org/wiki/Gamma_distribution).
    Parameter alpha controls the distribution's shape and should be a positive integer
    (if a non-integer is provided, the value is rounded.) When alpha=1, the
    distribution is the same as exponential. As its value increases, the
    probability density function becomes a curve with mean=alpha and standard
    deviation = sqrt(alpha).
    """
    r = 1.0
    n = round(alpha)
    for i in range(n):
        r = r * (1 - random.random())
    return - math.log(r)

def poisson(alpha=1):
    """
    Returns value greater than 0 from the poisson distribution
    (https://en.wikipedia.org/wiki/Poisson_distribution).
    Returns positive integer values, theoretically unbounded but practically
    limited by the alpha parameter. Parameter alpha controls the distribution's
    shape and must be positive: the mean is a and the standard deviation is sqrt(a).
    """
    b = math.exp( - alpha)
    n = 0
    p = 1.0
    while True:
        p = p * random.random()
        n += 1
        if p < b:
            return n

def expran(alpha=1):
    """
    Returns a value greater than 0 from the exponential distribution
    (https://en.wikipedia.org/wiki/Exponential_distribution) with stretching factor alpha.
    Increasing alpha "prefers" smaller numbers. The distribution is unbounded but when
    alpha=1 then %99.9 of the time the value returned will be than 6.9077554, i.e. -log(.001).
    The distribution density is f(x)=(exp -x) with a mean of 1.0.
    """
    return (- math.log(1.0 - random.random())) / alpha

def gauss(sigma=1, mu=0):
    """
    Returns unbounded value from the normal distribution
    (https://en.wikipedia.org/wiki/Normal_distribution) with standard deviation
    sigma and mean mu.  The spread (standard deviation) is 1.0 centered at 0,
    so 68.26% of the results are between -1 and 1 inclusive and 99.74% of the
    results are between -3 and 3 inclusive.
    """
    while True:
        x = -1 + 2 * random.random()
        y = -1 + 2 * random.random()
        r2 = (x * x) + (y * y)
        if not ((r2 > 1.0) or (r2 == 0)):
            break
    return (sigma * y * math.sqrt(-2.0 * math.log(r2) / r2)) + mu

def cauchy(alpha=False):
    """
    Returns an unbounded value from the Cauchy distribution
    (https://en.wikipedia.org/wiki/Cauchy_distribution).
    The density function is a bell shaped curve centered at 0 similar to a
    normal distribution but with more values at the extremes. The mean and
    standard deviation of the Cauchy distribution are undefined. If parameter
    alpha is True then only positive values are returned. Density is f(x)=1/(pi(1 + x^2)).
    """
    #return math.tan( math.pi * 2 * (random.random() - 0.5))
    r = (1.5707964 if alpha else math.pi) * random.random()
    return sin(r) / cos(r)

#-----------------------------------------------------------------------------
# Noises white, brown and pink return sample values -1.0 to 1.0

def white(level=None):
    """
    Returns white (uniform) noise samples between -1.0 and 1.0.

    Parameters
    ----------
    level : None or int or float
        An optional scaler on the sample returned.
    """
    samp = random.random() * 2.0 - 1.0
    return samp * level if level else samp

# def _lowpass(current, previous, alpha):
#     """
#     The change from one filter output to the next is proportional to the
#     difference between the previous output and the next input.
#     """
#     return previous + (alpha * (current - previous))
#
# _brown = white()
# """ The previous brown sample, initialzed to a random number."""
#
# def brown(alpha=.02):
#     global _brown
#     samp = _lowpass(white(), _brown, alpha)
#     _brown = samp
#     return samp * 3.0  # boost the overall gain

_brown = 0.0

def brown():
    """
    Returns brownish (1/f**2) noise samples between -1 and 1.
    """
    # from http://vellocet.com/dsp/noise/VRand.h
    # but made to generate between -1 1
    global _brown
    while True:
        r = random.random() * 2 - 1
        _brown += r
        if (_brown < -16.0) or (_brown > 16.0):
            _brown -= r
        else:
            break
    # return interval -1 1.
    return _brown * 0.0625

_pow2 = 5
_pown = 32
_pinking = [0.0 for _ in range(_pown+1)]

def _one_over_f_aux(n, rarray, halfrange):
    s = 0.0
    for i in range(_pow2):
        p = 2.0 ** i
        if not ( (n / p) == ((n - 1) / p) ):
            rarray[i]=( (random.random() * 2 * halfrange) - halfrange)
        s += rarray[i]
    return s

_i = _pown

def pink():
    """
    Returns pinkish (1/f) noise samples between -1.0 and 1.0.  Based on Gardner
    (1978) and Dick Moore (1988?).
    """
    global _i
    if _i == _pown:
        r = [0.0 for _ in range(_pow2)]
        h = 1.0 / _pow2
        for n in range(_pown):
            _pinking[n]=_one_over_f_aux(n, r, h)
        _i = 0
    else:
        _i += 1
    return _pinking[_i]

if __name__ == '__main__':

    ransegs(num=5, mapto=10)
    ransegs(num=5, mapto=(0,4))
    ransegs(num=5, mapto=(-2,4))

    # x = 60
    # for i in range(20):
    #     print(x)
    #     x = drunk(x, width=2)
    # print(x)

    # for i in range(20):
    #     print(pink())

    # for i in range(20):
    #     print(brown())

    import numpy as np
    import matplotlib.pyplot as plt
    # beta(.3, .3)
    # x = [pink() for _ in range(10000)]
    # plt.hist(x, bins=20, facecolor="blue", alpha=0.5)
    # plt.show()

#     y = [expran(1.5) for _ in range(1000)]
#     plt.hist(y, bins=20, facecolor="blue", alpha=0.5)
#    #plt.plot(y)
#     plt.show()

    #ransegs(5, 2, 4)
    #ransegs(5, 0, 10)
    #ransegs(5, 0, 10)
    #ransegs(5, -10, 10)
    #ransegs(5, sumto=10)

    #x = ransegs(5, 10)
    #print("sum=", sum(x))

"""
(define (segs num sum)
  (let ((rgen ran)
        (args ()))
    (let* ((rsum (apply rgen args))
           (head (list rsum))
           (tail head))
      (print "rgen " rgen " rsum " rsum " args " args " head" head)
      (do ((i 1 (+ i 1))
           (n #f))
          ((not (< i num))
           (do ((tail head (cdr tail)))
               ((null? tail) head)
             (set-car! tail
                       (ffi_rescale (car tail) 0 rsum 0 sum 1))))
        (set! n (apply rgen args))
        (set-cdr! tail (list n))
        (set! tail (cdr tail))
        (set! rsum (+ rsum n))))
    ))
(let ((l (segs 7 2)))
  (print l)
  (apply + l))
"""

max_stars_repo_name: ricktaube/musx

max_issues_repo_name: musx-admin/musx

max_forks_repo_name: musx-admin/musx


text:
"""
POVM effect representation classes for the `densitymx_slow` evolution type.
"""
#***************************************************************************************************
# Copyright 2015, 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights
# in this software.
# 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 or in the LICENSE file in the root pyGSTi directory.
#***************************************************************************************************

import numpy as _np

# import functools as _functools
from .. import basereps as _basereps
from pygsti.baseobjs.statespace import StateSpace as _StateSpace
from ...tools import matrixtools as _mt

class EffectRep(_basereps.EffectRep):
    def __init__(self, state_space):
        self.state_space = _StateSpace.cast(state_space)

    def probability(self, state):
        raise NotImplementedError()

class EffectRepConjugatedState(EffectRep):

    def __init__(self, state_rep):
        self.state_rep = state_rep
        super(EffectRepConjugatedState, self).__init__(state_rep.state_space)

    def __reduce__(self):
        return (EffectRepConjugatedState, (self.state_rep,))

    def probability(self, state):
        # can assume state is a StateRep and self.state_rep is
        return _np.dot(self.state_rep.data, state.data)  # not vdot b/c *real* data

    def to_dense(self, on_space):
        return self.state_rep.to_dense(on_space)

class EffectRepComputational(EffectRep):

    def __init__(self, zvals, basis, state_space):
        state_space = _StateSpace.cast(state_space)
        assert(basis.name == 'pp'), "Only Pauli-product computational effect vectors are currently supported"
        assert(state_space.num_qudits == len(zvals))
        assert(len(zvals) <= 64), "Cannot create a Computational basis rep with >64 qubits!"
        # Current storage of computational basis states converts zvals -> 64-bit integer

        base = 1
        self.zvals_int = 0
        for v in zvals:
            assert(v in (0, 1)), "zvals must contain only 0s and 1s"
            self.zvals_int += base * v
            base *= 2  # or left shift?

        self.zvals = zvals
        self.nfactors = len(zvals)  # (or nQubits)
        self.abs_elval = 1 / (_np.sqrt(2)**self.nfactors)
        self.basis = basis

        super(EffectRepComputational, self).__init__(state_space)

    def __reduce__(self):
        return (EffectRepComputational, (self.zvals, self.basis, self.state_space))

    def probability(self, state):
        scratch = _np.empty(self.state_space.dim, 'd')
        Edense = self.to_dense('HilbertSchmidt', scratch)
        return _np.dot(Edense, state.data)  # not vdot b/c data is *real*

    def to_dense(self, on_space, outvec=None):
        if on_space not in ('minimal', 'HilbertSchmidt'):
            raise ValueError("'densitymx' evotype cannot produce Hilbert-space ops!")
        return _mt.zvals_int64_to_dense(self.zvals_int, self.nfactors, outvec, False, self.abs_elval)

class EffectRepTensorProduct(EffectRep):

    def __init__(self, povm_factors, effect_labels, state_space):
        #Arrays for speeding up kron product in effect reps
        max_factor_dim = max(fct.state_space.dim for fct in povm_factors)
        kron_array = _np.ascontiguousarray(
            _np.empty((len(povm_factors), max_factor_dim), 'd'))
        factordims = _np.ascontiguousarray(
            _np.array([fct.state_space.dim for fct in povm_factors], _np.int64))

        self.povm_factors = povm_factors
        self.effect_labels = effect_labels
        self.kron_array = kron_array
        self.factor_dims = factordims
        self.max_factor_dim = max_factor_dim  # Unused
        state_space = _StateSpace.cast(state_space)
        assert(_np.product(factordims) == state_space.dim)
        super(EffectRepTensorProduct, self).__init__(state_space)
        self.factor_effects_have_changed()

    #TODO: fix this:
    #def __reduce__(self):
    #    return (EffectRepTensorProduct,
    #            (self.kron_array, self.factor_dims, self.nfactors, self.max_factor_dim, self.dim))

    def to_dense(self, on_space, outvec=None):

        if on_space not in ('minimal', 'HilbertSchmidt'):
            raise ValueError("'densitymx' evotype cannot produce Hilbert-space ops!")

        if outvec is None:
            outvec = _np.zeros(self.state_space.dim, 'd')

        N = self.state_space.dim
        nfactors = len(self.povm_factors)
        #Put last factor at end of outvec
        k = nfactors - 1  # last factor
        off = N - self.factor_dims[k]  # offset into outvec
        for i in range(self.factor_dims[k]):
            outvec[off + i] = self.kron_array[k, i]
        sz = self.factor_dims[k]

        #Repeatedly scale&copy last "sz" elements of outputvec forward
        # (as many times as there are elements in the current factor array)
        # - but multiply *in-place* the last "sz" elements.
        for k in range(nfactors - 2, -1, -1):  # for all but the last factor
            off = N - sz * self.factor_dims[k]
            endoff = N - sz

            #For all but the final element of self.kron_array[k,:],
            # mult&copy final sz elements of outvec into position
            for j in range(self.factor_dims[k] - 1):
                mult = self.kron_array[k, j]
                for i in range(sz):
                    outvec[off + i] = mult * outvec[endoff + i]
                off += sz

            #Last element: in-place mult
            #assert(off == endoff)
            mult = self.kron_array[k, self.factor_dims[k] - 1]
            for i in range(sz):
                outvec[endoff + i] *= mult
            sz *= self.factor_dims[k]

        return outvec

    def probability(self, state):  # allow scratch to be passed in?
        scratch = _np.empty(self.dim, 'd')
        Edense = self.to_dense('HilbertSchmidt', scratch)
        return _np.dot(Edense, state.data)  # not vdot b/c data is *real*

    def _fill_fast_kron(self):
        """ Fills in self._fast_kron_array based on current self.factors """
        for i, (factor_dim, Elbl) in enumerate(zip(self.factor_dims, self.effect_labels)):
            self.kron_array[i][0:factor_dim] = self.povm_factors[i][Elbl].to_dense('HilbertSchmidt')

    def factor_effects_have_changed(self):
        self._fill_fast_kron()  # updates effect reps

    #def to_dense(self):
    #    if len(self.factors) == 0: return _np.empty(0, complex if self._evotype == "statevec" else 'd')
    #        #NOTE: moved a fast version of to_dense to replib - could use that if need a fast to_dense call...
    #
    #        factorPOVMs = self.factors
    #        ret = factorPOVMs[0][self.effectLbls[0]].to_dense()
    #        for i in range(1, len(factorPOVMs)):
    #            ret = _np.kron(ret, factorPOVMs[i][self.effectLbls[i]].to_dense())
    #        return ret
    #    elif self._evotype == "stabilizer":
    #        # each factor is a StabilizerEffectVec
    #        raise ValueError("Cannot convert Stabilizer tensor product effect to an array!")
    #    # should be using effect.outcomes property...
    #    else:  # self._evotype in ("svterm","cterm")
    #        raise NotImplementedError("to_dense() not implemented for %s evolution type" % self._evotype)

class EffectRepComposed(EffectRep):
    def __init__(self, op_rep, effect_rep, op_id, state_space):
        self.op_rep = op_rep
        self.effect_rep = effect_rep
        self.op_id = op_id

        self.state_space = _StateSpace.cast(state_space)
        assert(self.state_space.is_compatible_with(effect_rep.state_space))

        super(EffectRepComposed, self).__init__(effect_rep.state_space)

    def __reduce__(self):
        return (EffectRepComposed, (self.op_rep, self.effect_rep, self.op_id, self.state_space))

    def probability(self, state):
        state = self.op_rep.acton(state)  # *not* acton_adjoint
        return self.effect_rep.probability(state)

ext: py

max_stars_repo_name: pyGSTi-Developers/pyGSTi

max_issues_repo_name: pyGSTi-Developers/pyGSTi

max_forks_repo_name: pyGSTi-Developers/pyGSTi

reason: import numpy


text:
from sympy import *
from sympy.matrices import matrices

class MatrixReduction():
    '''
    带步骤输出的矩阵消元类
    '''
    def __init__(self, matrix: MatrixBase) -> None:
        '''
        Parameter
            matrix  待消元的矩阵
        '''
        self.matrix = matrix
        self.course = [matrix.copy(), ]
        return

    def echelon_form(self) -> list:
        '''
        化矩阵为阶梯形，返回一个列表，存有每次初等变换后的矩阵
        '''
        if not self.matrix.is_echelon:
            for rowInd in range(self.matrix.shape[0]-1):
                self.__reduce_col(rowInd)
        return

    def __reduce_col(self, n: int) -> None:
        #用于逐行消元，n是用于消去后面行元素的行
        row = self.matrix[n, :]
        pivotInd = self.pivot_Ind(row)
        for rowInd in range(n+1, self.matrix.shape[0]):
            if pivotInd > self.pivot_Ind(self.matrix[rowInd, :]):
                self.matrix = self.matrix.elementary_row_op('n<->m', n, rowInd) #交换行，使得长的行放在最上面
                self.course.append(self.matrix.copy())
        pivotInd = self.pivot_Ind(self.matrix[n, :])    #因为可能交换行，所以重置了主元
        for rowInd in range(n+1, self.matrix.shape[0]):
            if self.matrix[rowInd, pivotInd] != 0:
                k = Rational(-self.matrix[rowInd, pivotInd],self.matrix[n, pivotInd])
                self.matrix = self.matrix.elementary_row_op('n->n+km', rowInd, k, n)    #用第n行消去后面的行

        return  self.course.append(self.matrix.copy())

    def pivot_Ind(self, row: MatrixBase) -> int:
        #找主元，主元是每行元素中从左到右第一个不为0的元素
        for colInd in range(row.shape[1]):
            if row[colInd] != 0:
                return colInd
        return row.shape[1] #如果一行全为0，主元序号设置为最大索引+1

def reduced_echelon_form(mat: MatrixBase):

    return

max_stars_repo_name: HelloMrGeorge/LinearAlgebraMachine

max_issues_repo_name: HelloMrGeorge/LinearAlgebraMachine

max_forks_repo_name: HelloMrGeorge/LinearAlgebraMachine

reason: from sympy


text:
''' Place within the full tomogram a subtomogram average in the original position
and orientation for each of the original subtomograms making up the average

Download to edit and run: :download:`sta_place_average.py <../../arachnid/snippets/sta_place_average.py>`

To run:

.. sourcecode:: sh

    $ python sta_place_average.py

.. note::

    Requires EMAN2 2.1

.. literalinclude:: ../../arachnid/snippets/sta_place_average.py
   :language: python
   :lines: 23-
   :linenos:

.. codeauthor:: Robert Langlois <rl2528@columbia.edu>
'''
import sys
import EMAN2
import numpy
#import os

if __name__ == '__main__':

    # Parameters
    tomogram = sys.argv[1]      # Decimated tomogram - 8 times
    alignment = sys.argv[2]     # Undecimated alignment parameters from EMAN2 e2spt_classaverage.py
    orig_stack = sys.argv[3]    # Undecimated volumes extracted with EMAN2 e2spt_boxer.py
    output = sys.argv[4]        # Output placement sans tomogram
    average = sys.argv[5]    # Undecimated volumes extracted with EMAN2 e2spt_boxer.py
    mask_file = sys.argv[6] if len(sys.argv)>6 else None
    bin_factor = 8

    tomo = EMAN2.EMData()
    tomo.read_image(tomogram)
    fulltomo = tomo.copy()
    tomo.to_zero()

    aligndict = EMAN2.js_open_dict(alignment) # @UndefinedVariable
    coords={}
    for key in aligndict.keys():
        off = key.find('_')+1
        id = int(key[off:])
        oem = EMAN2.EMData()
        oem.read_image(orig_stack, id)
        coords[id]=numpy.asarray(oem['ptcl_source_coord'])/8

    print "Tomogram dimensions", tomo.get_xsize(), tomo.get_ysize(), tomo.get_zsize()
    avg = EMAN2.EMData()
    avg.read_image(average)
    nx = avg.get_xsize()/bin_factor
    ny = avg.get_ysize()/bin_factor
    nz = avg.get_zsize()/bin_factor

    if mask_file is not None:
        emmask = EMAN2.EMData()
        emmask.read_image(mask_file)
        emmask = emmask.process('math.fft.resample', dict(n=bin_factor))
        npmask = EMAN2.EMNumPy.em2numpy(emmask)
        print npmask.max(), npmask.min(), npmask.sum(), numpy.prod(npmask.shape), (npmask==0).sum(), (npmask==1).sum()
        mask = numpy.nonzero(npmask>0.99999)
        nptomo = EMAN2.EMNumPy.em2numpy(tomo)
    else:
        mask = None

    print 'Place', len(aligndict), 'subtomos'
    for key in aligndict.keys():
        off = key.find('_')+1
        id = int(key[off:])
        trans = aligndict[key][0]
        cavg = avg.copy()
        itrans = trans.inverse()
        cavg.process_inplace("xform",{"transform":itrans})
        cavg = cavg.process('math.fft.resample', dict(n=bin_factor))
        x,y,z = coords[id]
        x = (2*int(x) - nx)/2
        y = (2*int(y) - ny)/2
        z = (2*int(z) - nz)/2
        if mask is not None:
            cmask = emmask.copy()
            sx,sy,sz = itrans.get_trans()
            itrans.set_trans((sx/bin_factor, sy/bin_factor, sz/bin_factor))
            cmask.process_inplace("xform",{"transform":itrans})
            npmask = EMAN2.EMNumPy.em2numpy(cmask)
            mask = numpy.nonzero(npmask>0.5)
            if id==0: print 'cavg', cavg.get_xsize()
            npcavg = EMAN2.EMNumPy.em2numpy(cavg)
            #print z,z+nz, y,y+ny, x,x+nx
            npsubset = nptomo[z:z+nz, y:y+ny, x:x+nx]
            npsubset[mask] = npcavg[mask]
        else:
            assert(False)
            tomo.insert_clip(cavg, (x, y, z))

    print 'done:', output
    tomo.write_image(output)

max_stars_repo_name: ezralanglois/subtomo

max_issues_repo_name: ezralanglois/subtomo

max_forks_repo_name: ezralanglois/subtomo


text:
from albumentations import augmentations
from scipy.fftpack import dct, idct, rfft, irfft
import tensorflow as tf
import numpy as np
from albumentations import *
from random import randint, uniform
import PIL
import PIL.Image
from io import BytesIO
import cv2
import random

# This file contains the defense methods compared in the paper.
# The FD algorithm's source code is from:
#   https://github.com/zihaoliu123/Feature-Distillation-DNN-Oriented-JPEG-Compression-Against-Adversarial-Examples/blob/master/utils/feature_distillation.py
# The FD algorithm is refer to the paper:
#   https://arxiv.org/pdf/1803.05787.pdf
# Some of the defense methods' code refering to Anish & Carlini's github: https://github.com/anishathalye/obfuscated-gradients

def defended(funcname, adv):
    func = eval('defend_'+funcname)
    defendadv = np.zeros((adv.shape[0],299,299,3))
    for i in range(adv.shape[0]):
        defendadv[i] = func(adv[i])
    return defendadv

def remap(img,index):
    n = img.shape[0]
    m = img.shape[1]
    if np.ndim(img)>2:
        outimgnp = np.zeros((n,m,img.shape[2]))
    else:
        outimgnp = np.zeros((n,m))
    for i in range(n):
        for j in range(m):
            outimgnp[j,i]=img[tuple(index[i][j])]
    return outimgnp

# GD algorithm that runs in the session (for EOT)
def tftensorGD(img,distort_limit = 0.25):
    num_steps = 10
    upflag = tf.round(tf.random_uniform((1,), 0, 1, dtype=tf.float32))
    leftflag = tf.round(tf.random_uniform((1,), 0, 1, dtype=tf.float32))
    xstep = tf.constant(1.0) + tf.random_uniform((num_steps + 1,),
                                                 -distort_limit,
                                                 distort_limit,
                                                 dtype=tf.float32)
    ystep = tf.constant(1.0) + tf.random_uniform((num_steps + 1,),
                                                 -distort_limit,
                                                 distort_limit,
                                                 dtype=tf.float32)
    img_shape = tf.shape(img)
    height, width = img_shape[0], img_shape[1]
    x_step = width // num_steps
    y_step = height // num_steps
    xs = tf.range(0.0, tf.dtypes.cast(width, tf.float32), delta=x_step)
    ys = tf.range(0.0, tf.dtypes.cast(height, tf.float32), delta=y_step)
    prev = tf.constant(0.0)
    listvec_x = tf.zeros((1, 1))
    for i in range(num_steps + 1):
        start = tf.cast(xs[i], tf.int32)
        end = tf.cast(xs[i], tf.int32) + x_step
        cur = tf.cond(end > width, lambda: tf.cast(width, tf.float32),
                      lambda: prev + tf.cast(x_step, tf.float32) * xstep[i])
        end = tf.cond(end > width, lambda: width, lambda: end)
        listvec_x = tf.concat([listvec_x, tf.reshape(tf.linspace(prev, cur, end - start), (1, -1))], -1)
        prev = cur
    xx = tf.cast(tf.clip_by_value(tf.round(listvec_x), 0, 298), tf.int32)
    map_x = tf.tile(xx[:, 1:], (299, 1))
    xx2 = tf.reverse((298 * tf.ones_like(xx, dtype=tf.int32) - xx), [1])
    map_x2 = tf.tile(xx2[:, :299], (299, 1))
    prev = tf.constant(0.0)
    listvec_y = tf.zeros((1, 1))
    for i in range(num_steps + 1):
        start = tf.cast(ys[i], tf.int32)
        end = tf.cast(ys[i], tf.int32) + y_step
        cur = tf.cond(end > width, lambda: tf.cast(height, tf.float32),
                      lambda: prev + tf.cast(y_step, tf.float32) * ystep[i])
        end = tf.cond(end > width, lambda: width, lambda: end)
        listvec_y = tf.concat([listvec_y, tf.reshape(tf.linspace(prev, cur, end - start), (1, -1))], -1)
        prev = cur
    yy = tf.cast(tf.clip_by_value(tf.round(listvec_y), 0, 298), tf.int32)
    map_y = tf.tile(tf.transpose(yy)[1:, :], (1, 299))
    yy2 = tf.reverse((298 * tf.ones_like(yy, dtype=tf.int32) - yy), [1])
    map_y2 = tf.tile(tf.transpose(yy2)[:299, :], (1, 299))
    index_x = tf.cond(leftflag[0] > 0.5, lambda: tf.identity(map_x), lambda: tf.identity(map_x2))
    index_y = tf.cond(upflag[0] > 0.5, lambda: tf.identity(map_y), lambda: tf.identity(map_y2))
    index = tf.stack([index_y, index_x], 2)
    x_gd = tf.gather_nd(img, index)
    return x_gd

# GD algorithm that runs outside the session (for BPDA)
def defend_GD(img,distort_limit = 0.25):
    num_steps = 10

    xsteps = [1 + random.uniform(-distort_limit, distort_limit) for i in range(num_steps + 1)]
    ysteps = [1 + random.uniform(-distort_limit, distort_limit) for i in range(num_steps + 1)]

    height, width = img.shape[:2]

    x_step = width // num_steps
    xx = np.zeros(width, np.float32)
    prev = 0
    for idx, x in enumerate(range(0, width, x_step)):
        start = x
        end = x + x_step
        if end > width:
            end = width
            cur = width
        else:
            cur = prev + x_step * xsteps[idx]

        xx[start:end] = np.linspace(prev, cur, end - start)
        prev = cur

    y_step = height // num_steps
    yy = np.zeros(height, np.float32)
    prev = 0
    for idx, y in enumerate(range(0, height, y_step)):
        start = y
        end = y + y_step
        if end > height:
            end = height
            cur = height
        else:
            cur = prev + y_step * ysteps[idx]

        yy[start:end] = np.linspace(prev, cur, end - start)
        prev = cur
    xx = np.round(xx).astype(int)
    yy = np.round(yy).astype(int)
    xx[xx >= 299] = 298
    yy[yy >= 299] = 298

    map_x, map_y = np.meshgrid(xx, yy)

    #     index=np.dstack((map_y,map_x))
    #     outimg = remap(img,index)
    #     if np.ndim(img)>2:
    #         outimg = outimg.transpose(1,0,2)
    #     else:
    #         outimg = outimg.T

    # to speed up the mapping procedure, OpenCV 2 is adopted
    map_x = map_x.astype(np.float32)
    map_y = map_y.astype(np.float32)
    outimg = cv2.remap(img, map1=map_x, map2=map_y, interpolation=1, borderMode=4, borderValue=None)
    return outimg

# FD algorithm
T = np.array([
        [0.3536, 0.3536, 0.3536, 0.3536, 0.3536, 0.3536, 0.3536, 0.3536],
        [0.4904, 0.4157, 0.2778, 0.0975, -0.0975, -0.2778, -0.4157, -0.4904],
        [0.4619, 0.1913, -0.1913, -0.4619, -0.4619, -0.1913, 0.1913, 0.4619],
        [0.4157, -0.0975, -0.4904, -0.2778, 0.2778, 0.4904, 0.0975, -0.4157],
        [0.3536, -0.3536, -0.3536, 0.3536, 0.3536, -0.3536, -0.3536, 0.3536],
        [0.2778, -0.4904, 0.0975, 0.4157, -0.4157, -0.0975, 0.4904, -0.2778],
        [0.1913, -0.4619, 0.4619, -0.1913, -0.1913, 0.4619, -0.4619, 0.1913],
        [0.0975, -0.2778, 0.4157, -0.4904, 0.4904, -0.4157, 0.2778, -0.0975]
    ])
num = 8
q_table = np.ones((num,num))*30
q_table[0:4,0:4] = 25

def dct2 (block):
    return dct(dct(block.T, norm = 'ortho').T, norm = 'ortho')
def idct2(block):
    return idct(idct(block.T, norm = 'ortho').T, norm = 'ortho')
def rfft2 (block):
    return rfft(rfft(block.T).T)
def irfft2(block):
    return irfft(irfft(block.T).T)

# Feature distillation for batch
def FD_fuction(input_matrix):
    output = []
    input_matrix = input_matrix*255

    n = input_matrix.shape[0]
    h = input_matrix.shape[1]
    w = input_matrix.shape[2]
    c = input_matrix.shape[3]
    horizontal_blocks_num = w / num
    output2=np.zeros((c,h, w))
    output3=np.zeros((n,3,h, w))
    vertical_blocks_num = h / num
    n_block = np.split(input_matrix,n,axis=0)
    for i in range(n):
        c_block = np.split(n_block[i],c,axis =3)
        j=0
        for ch_block in c_block:
            vertical_blocks = np.split(ch_block, vertical_blocks_num,axis = 1)
            k=0
            for block_ver in vertical_blocks:
                hor_blocks = np.split(block_ver,horizontal_blocks_num,axis = 2)
                m=0
                for block in hor_blocks:
                    block = np.reshape(block,(num,num))
                    block = dct2(block)
                    # quantization
                    table_quantized = np.matrix.round(np.divide(block, q_table))
                    table_quantized = np.squeeze(np.asarray(table_quantized))
                    # de-quantization
                    table_unquantized = table_quantized*q_table
                    IDCT_table = idct2(table_unquantized)
                    if m==0:
                        output=IDCT_table
                    else:
                        output = np.concatenate((output,IDCT_table),axis=1)
                    m=m+1
                if k==0:
                    output1=output
                else:
                    output1 = np.concatenate((output1,output),axis=0)
                k=k+1
            output2[j] = output1
            j=j+1
        output3[i] = output2
    output3 = np.transpose(output3,(0,2,1,3))
    output3 = np.transpose(output3,(0,1,3,2))
    output3 = output3/255
    output3 = np.clip(np.float32(output3),0.0,1.0)
    return output3
def padresult(cleandata):
    pad = augmentations.transforms.PadIfNeeded(min_height=304, min_width=304, border_mode=4)
    paddata = np.ones((cleandata.shape[0],304,304,3))
    for i in range(paddata.shape[0]):
        paddata[i] = pad(image = cleandata[i])['image']
    return paddata
def cropresult(paddata):
    crop = augmentations.transforms.Crop(0,0,299,299)
    resultdata = np.ones((paddata.shape[0],299,299,3))
    for i in range(resultdata.shape[0]):
        resultdata[i] = crop(image = paddata[i])['image']
    return resultdata

def defend_FD(data):
    paddata = padresult(data)
    defendresult = FD_fuction(paddata)
    resultdata = cropresult(defendresult)
    return resultdata

# Feature distillation for single imput
def FD_fuction_sig(input_matrix):
    output = []
    input_matrix = input_matrix * 255

    h = input_matrix.shape[0]
    w = input_matrix.shape[1]
    c = input_matrix.shape[2]
    horizontal_blocks_num = w / num
    output2 = np.zeros((c, h, w))
    vertical_blocks_num = h / num

    c_block = np.split(input_matrix, c, axis=2)
    j = 0
    for ch_block in c_block:
        vertical_blocks = np.split(ch_block, vertical_blocks_num, axis=0)
        k = 0
        for block_ver in vertical_blocks:
            hor_blocks = np.split(block_ver, horizontal_blocks_num, axis=1)
            m = 0
            for block in hor_blocks:
                block = np.reshape(block, (num, num))
                block = dct2(block)
                # quantization
                table_quantized = np.matrix.round(np.divide(block, q_table))
                table_quantized = np.squeeze(np.asarray(table_quantized))
                # de-quantization
                table_unquantized = table_quantized * q_table
                IDCT_table = idct2(table_unquantized)
                if m == 0:
                    output = IDCT_table
                else:
                    output = np.concatenate((output, IDCT_table), axis=1)
                m = m + 1
            if k == 0:
                output1 = output
            else:
                output1 = np.concatenate((output1, output), axis=0)
            k = k + 1
        output2[j] = output1
        j = j + 1

    output2 = np.transpose(output2, (1, 0, 2))
    output2 = np.transpose(output2, (0, 2, 1))
    output2 = output2 / 255
    output2 = np.clip(np.float32(output2), 0.0, 1.0)
    return output2
def padresult_sig(cleandata):
    pad = augmentations.transforms.PadIfNeeded(min_height=304, min_width=304, border_mode=4)
    paddata = pad(image=cleandata)['image']
    return paddata
def cropresult_sig(paddata):
    crop = augmentations.transforms.Crop(0, 0, 299, 299)
    resultdata = crop(image=paddata)['image']
    return resultdata

def defend_FD_sig(data):
    paddata = padresult_sig(data)
    defendresult = FD_fuction_sig(paddata)
    resultdata = cropresult_sig(defendresult)
    return resultdata

#####################################################################
# Defense methods used to compare with different indifferentiable functions that combined with GD algorithm
def rand_jpeg(img):
    jpeg = JpegCompression(quality_lower=10, quality_upper=80, p=1)
    augmented = jpeg(image=((img*255).astype(np.uint8)))
    auged = (augmented['image']/255).astype(float)
    return auged

def defend_RCDfense(img):
    # attach an random Jpeg compression procedure before the GD algorithm
    img = rand_jpeg(img)
    num_steps = 10
    distort_limit = 0.1

    xsteps = [1 + random.uniform(-distort_limit, distort_limit) for i in range(num_steps + 1)]
    ysteps = [1 + random.uniform(-distort_limit, distort_limit) for i in range(num_steps + 1)]

    height, width = img.shape[:2]

    x_step = width // num_steps
    xx = np.zeros(width, np.float32)
    prev = 0
    for idx, x in enumerate(range(0, width, x_step)):
        start = x
        end = x + x_step
        if end > width:
            end = width
            cur = width
        else:
            cur = prev + x_step * xsteps[idx]

        xx[start:end] = np.linspace(prev, cur, end - start)
        prev = cur

    y_step = height // num_steps
    yy = np.zeros(height, np.float32)
    prev = 0
    for idx, y in enumerate(range(0, height, y_step)):
        start = y
        end = y + y_step
        if end > height:
            end = height
            cur = height
        else:
            cur = prev + y_step * ysteps[idx]

        yy[start:end] = np.linspace(prev, cur, end - start)
        prev = cur
    xx = np.round(xx).astype(int)
    yy = np.round(yy).astype(int)
    xx[xx >= 299] = 298
    yy[yy >= 299] = 298

    map_x, map_y = np.meshgrid(xx, yy)
    #     index=np.dstack((map_y,map_x))

    #     outimg = remap(img,index)
    #     if np.ndim(img)>2:
    #         outimg = outimg.transpose(1,0,2)
    #     else:
    #         outimg = outimg.T

    map_x = map_x.astype(np.float32)
    map_y = map_y.astype(np.float32)
    outimg = cv2.remap(img, map1=map_x, map2=map_y, interpolation=1, borderMode=4, borderValue=None)

    return outimg

#####################################################################
# Defense methods used to compare the performance over the BPDA

def defend_BitReduct(arr, depth=3):
    arr = (arr * 255.0).astype(np.uint8)
    shift = 8 - depth
    arr = (arr >> shift) << shift
    arr = arr.astype(np.float32)/255.0
    return arr

def defend_FixedJpeg(input_array):
    pil_image = PIL.Image.fromarray((input_array*255.0).astype(np.uint8))
    f = BytesIO()
    pil_image.save(f, format='jpeg', quality=75) # quality level specified in paper
    jpeg_image = np.asarray(PIL.Image.open(f)).astype(np.float32)/255.0
    return jpeg_image

def bregman(image, mask, weight, eps=1e-3, max_iter=100):
    rows, cols, dims = image.shape
    rows2 = rows + 2
    cols2 = cols + 2
    total = rows * cols * dims
    shape_ext = (rows2, cols2, dims)

    u = np.zeros(shape_ext)
    dx = np.zeros(shape_ext)
    dy = np.zeros(shape_ext)
    bx = np.zeros(shape_ext)
    by = np.zeros(shape_ext)

    u[1:-1, 1:-1] = image
    # reflect image
    u[0, 1:-1] = image[1, :]
    u[1:-1, 0] = image[:, 1]
    u[-1, 1:-1] = image[-2, :]
    u[1:-1, -1] = image[:, -2]

    i = 0
    rmse = np.inf
    lam = 2 * weight
    norm = (weight + 4 * lam)

    while i < max_iter and rmse > eps:
        rmse = 0

        for k in range(dims):
            for r in range(1, rows + 1):
                for c in range(1, cols + 1):
                    uprev = u[r, c, k]

                    # forward derivatives
                    ux = u[r, c + 1, k] - uprev
                    uy = u[r + 1, c, k] - uprev

                    # Gauss-Seidel method
                    if mask[r - 1, c - 1]:
                        unew = (lam * (u[r + 1, c, k] +
                                       u[r - 1, c, k] +
                                       u[r, c + 1, k] +
                                       u[r, c - 1, k] +
                                       dx[r, c - 1, k] -
                                       dx[r, c, k] +
                                       dy[r - 1, c, k] -
                                       dy[r, c, k] -
                                       bx[r, c - 1, k] +
                                       bx[r, c, k] -
                                       by[r - 1, c, k] +
                                       by[r, c, k]
                                       ) + weight * image[r - 1, c - 1, k]
                                ) / norm
                    else:
                        # similar to the update step above, except we take
                        # lim_{weight->0} of the update step, effectively
                        # ignoring the l2 loss
                        unew = (u[r + 1, c, k] +
                                u[r - 1, c, k] +
                                u[r, c + 1, k] +
                                u[r, c - 1, k] +
                                dx[r, c - 1, k] -
                                dx[r, c, k] +
                                dy[r - 1, c, k] -
                                dy[r, c, k] -
                                bx[r, c - 1, k] +
                                bx[r, c, k] -
                                by[r - 1, c, k] +
                                by[r, c, k]
                                ) / 4.0
                    u[r, c, k] = unew

                    # update rms error
                    rmse += (unew - uprev) ** 2

                    bxx = bx[r, c, k]
                    byy = by[r, c, k]

                    # d_subproblem
                    s = ux + bxx
                    if s > 1 / lam:
                        dxx = s - 1 / lam
                    elif s < -1 / lam:
                        dxx = s + 1 / lam
                    else:
                        dxx = 0
                    s = uy + byy
                    if s > 1 / lam:
                        dyy = s - 1 / lam
                    elif s < -1 / lam:
                        dyy = s + 1 / lam
                    else:
                        dyy = 0

                    dx[r, c, k] = dxx
                    dy[r, c, k] = dyy

                    bx[r, c, k] += ux - dxx
                    by[r, c, k] += uy - dyy

        rmse = np.sqrt(rmse / total)
        i += 1

    return np.squeeze(np.asarray(u[1:-1, 1:-1]))
def defend_TotalVarience(input_array, keep_prob=0.5, lambda_tv=0.03):
    mask = np.random.uniform(size=input_array.shape[:2])
    mask = mask < keep_prob
    return bregman(input_array, mask, weight=2.0 / lambda_tv)

def make_defend_quilt(sess):
    # setup for quilting
    quilt_db = np.load('data/quilt_db.npy')
    quilt_db_reshaped = quilt_db.reshape(1000000, -1)
    TILE_SIZE = 5
    TILE_OVERLAP = 2
    tile_skip = TILE_SIZE - TILE_OVERLAP
    K = 10
    db_tensor = tf.placeholder(tf.float32, quilt_db_reshaped.shape)
    query_imgs = tf.placeholder(tf.float32, (TILE_SIZE * TILE_SIZE * 3, None))
    norms = tf.reduce_sum(tf.square(db_tensor), axis=1)[:, tf.newaxis] \
            - 2 * tf.matmul(db_tensor, query_imgs)
    _, topk_indices = tf.nn.top_k(-tf.transpose(norms), k=K, sorted=False)

    def min_error_table(arr, direction):
        assert direction in ('horizontal', 'vertical')
        y, x = arr.shape
        cum = np.zeros_like(arr)
        if direction == 'horizontal':
            cum[:, -1] = arr[:, -1]
            for ix in range(x - 2, -1, -1):
                for iy in range(y):
                    m = arr[iy, ix + 1]
                    if iy > 0:
                        m = min(m, arr[iy - 1, ix + 1])
                    if iy < y - 1:
                        m = min(m, arr[iy + 1, ix + 1])
                    cum[iy, ix] = arr[iy, ix] + m
        elif direction == 'vertical':
            cum[-1, :] = arr[-1, :]
            for iy in range(y - 2, -1, -1):
                for ix in range(x):
                    m = arr[iy + 1, ix]
                    if ix > 0:
                        m = min(m, arr[iy + 1, ix - 1])
                    if ix < x - 1:
                        m = min(m, arr[iy + 1, ix + 1])
                    cum[iy, ix] = arr[iy, ix] + m
        return cum

    def index_exists(arr, index):
        if arr.ndim != len(index):
            return False
        return all(i > 0 for i in index) and all(index[i] < arr.shape[i] for i in range(arr.ndim))

    def assign_block(ix, iy, tile, synth):
        posx = tile_skip * ix
        posy = tile_skip * iy

        if ix == 0 and iy == 0:
            synth[posy:posy + TILE_SIZE, posx:posx + TILE_SIZE, :] = tile
        elif iy == 0:
            # first row, only have horizontal overlap of the block
            tile_left = tile[:, :TILE_OVERLAP, :]
            synth_right = synth[:TILE_SIZE, posx:posx + TILE_OVERLAP, :]
            errors = np.sum(np.square(tile_left - synth_right), axis=2)
            table = min_error_table(errors, direction='vertical')
            # copy row by row into synth
            xoff = np.argmin(table[0, :])
            synth[posy, posx + xoff:posx + TILE_SIZE] = tile[0, xoff:]
            for yoff in range(1, TILE_SIZE):
                # explore nearby xoffs
                candidates = [(yoff, xoff), (yoff, xoff - 1), (yoff, xoff + 1)]
                index = min((i for i in candidates if index_exists(table, i)), key=lambda i: table[i])
                xoff = index[1]
                synth[posy + yoff, posx + xoff:posx + TILE_SIZE] = tile[yoff, xoff:]
        elif ix == 0:
            # first column, only have vertical overlap of the block
            tile_up = tile[:TILE_OVERLAP, :, :]
            synth_bottom = synth[posy:posy + TILE_OVERLAP, :TILE_SIZE, :]
            errors = np.sum(np.square(tile_up - synth_bottom), axis=2)
            table = min_error_table(errors, direction='horizontal')
            # copy column by column into synth
            yoff = np.argmin(table[:, 0])
            synth[posy + yoff:posy + TILE_SIZE, posx] = tile[yoff:, 0]
            for xoff in range(1, TILE_SIZE):
                # explore nearby yoffs
                candidates = [(yoff, xoff), (yoff - 1, xoff), (yoff + 1, xoff)]
                index = min((i for i in candidates if index_exists(table, i)), key=lambda i: table[i])
                yoff = index[0]
                synth[posy + yoff:posy + TILE_SIZE, posx + xoff] = tile[yoff:, xoff]
        else:
            # glue cuts along diagonal
            tile_up = tile[:TILE_OVERLAP, :, :]
            synth_bottom = synth[posy:posy + TILE_OVERLAP, :TILE_SIZE, :]
            errors_up = np.sum(np.square(tile_up - synth_bottom), axis=2)
            table_up = min_error_table(errors_up, direction='horizontal')
            tile_left = tile[:, :TILE_OVERLAP, :]
            synth_right = synth[:TILE_SIZE, posx:posx + TILE_OVERLAP, :]
            errors_left = np.sum(np.square(tile_left - synth_right), axis=2)
            table_left = min_error_table(errors_left, direction='vertical')
            glue_index = -1
            glue_value = np.inf
            for i in range(TILE_OVERLAP):
                e = table_up[i, i] + table_left[i, i]
                if e < glue_value:
                    glue_value = e
                    glue_index = i
            # copy left part first, up to the overlap column
            xoff = glue_index
            synth[posy + glue_index, posx + xoff:posx + TILE_OVERLAP] = tile[glue_index, xoff:TILE_OVERLAP]
            for yoff in range(glue_index + 1, TILE_SIZE):
                # explore nearby xoffs
                candidates = [(yoff, xoff), (yoff, xoff - 1), (yoff, xoff + 1)]
                index = min((i for i in candidates if index_exists(table_left, i)), key=lambda i: table_left[i])
                xoff = index[1]
                synth[posy + yoff, posx + xoff:posx + TILE_OVERLAP] = tile[yoff, xoff:TILE_OVERLAP]
            # copy right part, down to overlap row
            yoff = glue_index
            synth[posy + yoff:posy + TILE_OVERLAP, posx + glue_index] = tile[yoff:TILE_OVERLAP, glue_index]
            for xoff in range(glue_index + 1, TILE_SIZE):
                # explore nearby yoffs
                candidates = [(yoff, xoff), (yoff - 1, xoff), (yoff + 1, xoff)]
                index = min((i for i in candidates if index_exists(table_up, i)), key=lambda i: table_up[i])
                yoff = index[0]
                synth[posy + yoff:posy + TILE_OVERLAP, posx + xoff] = tile[yoff:TILE_OVERLAP, xoff]
            # copy rest of image
            synth[posy + TILE_OVERLAP:posy + TILE_SIZE, posx + TILE_OVERLAP:posx + TILE_SIZE] = tile[TILE_OVERLAP:,
                                                                                TILE_OVERLAP:]

    KNN_MAX_BATCH = 1000

    def quilt(arr, graphcut=True):
        h, w, c = arr.shape
        assert (h - TILE_SIZE) % tile_skip == 0
        assert (w - TILE_SIZE) % tile_skip == 0
        horiz_blocks = (w - TILE_SIZE) // tile_skip + 1
        vert_blocks = (h - TILE_SIZE) // tile_skip + 1
        num_patches = horiz_blocks * vert_blocks
        patches = np.zeros((TILE_SIZE * TILE_SIZE * 3, num_patches))
        idx = 0
        for iy in range(vert_blocks):
            for ix in range(horiz_blocks):
                posx = tile_skip * ix
                posy = tile_skip * iy
                patches[:, idx] = arr[posy:posy + TILE_SIZE, posx:posx + TILE_SIZE, :].ravel()
                idx += 1

        ind = []
        for chunk in range(num_patches // KNN_MAX_BATCH + (1 if num_patches % KNN_MAX_BATCH != 0 else 0)):
            start = KNN_MAX_BATCH * chunk
            end = start + KNN_MAX_BATCH
            # for some reason, the code below is 10x slower when run in a Jupyter notebook
            # not sure why...
            indices_ = sess.run(topk_indices, {db_tensor: quilt_db_reshaped, query_imgs: patches[:, start:end]})
            for i in indices_:
                ind.append(np.random.choice(i))

        synth = np.zeros((299, 299, 3))

        idx = 0
        for iy in range(vert_blocks):
            for ix in range(horiz_blocks):
                posx = tile_skip * ix
                posy = tile_skip * iy
                tile = quilt_db[ind[idx]]
                if not graphcut:
                    synth[posy:posy + TILE_SIZE, posx:posx + TILE_SIZE, :] = tile
                else:
                    assign_block(ix, iy, tile, synth)
                idx += 1
        return synth

    return quilt

def nearest_neighbour_scaling(label, new_h, new_w, patch_size, patch_n, patch_m):
    if len(label.shape) == 2:
        label_new = np.zeros([new_h, new_w])
    else:
        label_new = np.zeros([new_h, new_w, label.shape[2]])
    n_pos = np.arange(patch_n)
    m_pos = np.arange(patch_m)
    n_pos = n_pos.repeat(patch_size)[:299]
    m_pos = m_pos.repeat(patch_size)[:299]
    n_pos = n_pos.reshape(n_pos.shape[0], 1)
    n_pos = np.tile(n_pos, (1, new_w))
    m_pos = np.tile(m_pos, (new_h, 1))
    assert n_pos.shape == m_pos.shape
    label_new[:, :] = label[n_pos[:, :], m_pos[:, :]]
    return label_new
def jpeg(input_array, quali):
    pil_image = PIL.Image.fromarray((input_array * 255.0).astype(np.uint8))
    f = BytesIO()
    pil_image.save(f, format='jpeg', quality=quali)  # quality level specified in paper
    jpeg_image = np.asarray(PIL.Image.open(f)).astype(np.float32) / 255.0
    return jpeg_image
def defend_SHIELD(x, qualities=(20, 40, 60, 80), patch_size=8):
    n = x.shape[0]
    m = x.shape[1]
    patch_n = n / patch_size
    patch_m = m / patch_size
    num_qualities = len(qualities)
    n = x.shape[0]
    m = x.shape[1]
    if n % patch_size > 0:
        patch_n = np.int(n / patch_size) + 1
        delete_n = 1
    if m % patch_size > 0:
        patch_m = np.int(m / patch_size) + 1
        delet_m = 1

    R = np.tile(np.reshape(np.arange(n), (n, 1)), [1, m])
    C = np.reshape(np.tile(np.arange(m), [n]), (n, m))
    mini_Z = (np.random.rand(patch_n, patch_m) * num_qualities).astype(int)
    Z = (nearest_neighbour_scaling(mini_Z, n, m, patch_size, patch_n, patch_m)).astype(int)
    indices = np.transpose(np.stack((Z, R, C)), (1, 2, 0))
    # x = img_as_ubyte(x)
    x_compressed_stack = []

    for quali in qualities:
        processed = jpeg(x, quali)
        x_compressed_stack.append(processed)

    x_compressed_stack = np.asarray(x_compressed_stack)
    x_slq = np.zeros((n, m, 3))
    for i in range(n):
        for j in range(m):
            x_slq[i, j] = x_compressed_stack[tuple(indices[i][j])]
    return x_slq

def defend_pixel_deflection(img, deflections=600, window=10):
    img = np.copy(img)
    H, W, C = img.shape
    while deflections > 0:
        #for consistency, when we deflect the given pixel from all the three channels.
        for c in range(C):
            x,y = randint(0,H-1), randint(0,W-1)
            while True: #this is to ensure that PD pixel lies inside the image
                a,b = randint(-1*window,window), randint(-1*window,window)
                if x+a < H and x+a > 0 and y+b < W and y+b > 0: break
            # calling pixel deflection as pixel swap would be a misnomer,
            # as we can see below, it is one way copy
            img[x,y,c] = img[x+a,y+b,c]
        deflections -= 1
    return img

#####################################################################
# Defense methods used to compare the performance over the BPDA + EOT

# x is a square image (3-tensor)
def defend_crop(x, crop_size=90, ensemble_size=30):
    x_size = tf.to_float(x.shape[1])
    frac = crop_size/x_size
    start_fraction_max = (x_size - crop_size)/x_size
    def randomizing_crop(x):
        start_x = tf.random_uniform((), 0, start_fraction_max)
        start_y = tf.random_uniform((), 0, start_fraction_max)
        return tf.image.crop_and_resize([x], boxes=[[start_y, start_x, start_y+frac, start_x+frac]],
                                 box_ind=[0], crop_size=[crop_size, crop_size])

    return tf.concat([randomizing_crop(x) for _ in range(ensemble_size)], axis=0)

# for single sample
def randomizing_crop(x):
    x_size = tf.to_float(x.shape[1])
    crop_size=90
    frac = crop_size/x_size
    start_fraction_max = (x_size - crop_size)/x_size
    start_x = tf.random_uniform((), 0, start_fraction_max)
    start_y = tf.random_uniform((), 0, start_fraction_max)
    return tf.image.crop_and_resize([x], boxes=[[start_y, start_x, start_y+frac, start_x+frac]],
                             box_ind=[0], crop_size=[crop_size, crop_size])

PAD_VALUE = 0.5
def defend_randomization(input_tensor):
    rnd = tf.random_uniform((), 299, 400, dtype=tf.int32)
    rescaled = tf.image.crop_and_resize(input_tensor, [[0, 0, 1, 1]], [0], [rnd, rnd])
    h_rem = 400 - rnd
    w_rem = 400 - rnd
    pad_left = tf.random_uniform((), 0, w_rem, dtype=tf.int32)
    pad_right = w_rem - pad_left
    pad_top = tf.random_uniform((), 0, h_rem, dtype=tf.int32)
    pad_bottom = h_rem - pad_top
    padded = tf.pad(rescaled, [[0, 0], [pad_top, pad_bottom], [pad_left, pad_right], [0, 0]], constant_values=PAD_VALUE)
    padded.set_shape((1, 400, 400, 3))
    return padded

from skimage import transform
def defend_onlyrand(img):
    rnd = np.random.randint(299,400,(1,))[0]
    rescaled = transform.resize(img,(rnd,rnd))
    h_rem = 400 - rnd
    w_rem = 400 - rnd
    pad_left = np.random.randint(0,w_rem,(1,))[0]
    pad_right = w_rem - pad_left
    pad_top = np.random.randint(0,h_rem,(1,))[0]
    pad_bottom = h_rem - pad_top
    padded = np.pad(rescaled,((pad_top,pad_bottom),(pad_left,pad_right),(0,0)),'constant',constant_values = 0.5)
    padded = transform.resize(padded,(299,299))
    return padded

reason: import numpy,from scipy


text:
import pytest
import numpy as np

from profiler import Cylinder

def test_cyclinder_class_should_accept_diameter_and_length():
    c = Cylinder(diameter=10, length=2)

    assert 10 == c.diameter
    assert 2 == c.length

def test_cylinder_class_should_fail_if_diameter_in_not_number():
    with pytest.raises(ValueError):
        c = Cylinder(diameter='a')

def test_cylinder_class_should_fail_if_length_in_not_number():
    with pytest.raises(ValueError):
        c = Cylinder(length='a')

def test_cylinder_class_should_calculate_area_correctly():
    c = Cylinder(diameter=10,length=2)

    expected_area = 78.5398

    assert np.isclose(expected_area, c.area, rtol=1e-4, atol=1e-6, equal_nan= False)

def test_cylinder_class_should_calculate_volume_correctly():
    c = Cylinder(diameter=10, length=2)

    expected_volume = 157.0796

    assert np.isclose(expected_volume, c.volume, rtol=1e-4, atol=1e-6,equal_nan=False)

def test_two_cylinder_class_should_add_together():
    c0 = Cylinder(diameter=10, length=4)
    c1 = Cylinder(diameter=10, length=2)
    c = c0 - c1
    expected_volume = c0.volume - c1.volume

    assert expected_volume == c.volume

def test_calculate_cylinder_volume_works():
    c = Cylinder(diameter=10, length=2)
    expected_volume = 157.0796

    assert np.isclose(expected_volume, c.calculate_volume(10,2), rtol=1e-4, atol=1e-6,equal_nan=False)

def test_calculate_cylinder_area_works():

    diameter = 3

    expected_area = np.pi * (diameter/2)**2

    assert expected_area == Cylinder.calculate_area(diameter)

max_stars_repo_name: NOAA-PMEL/EDD-DepthSimulator

max_issues_repo_name: NOAA-PMEL/EDD-DepthSimulator

max_forks_repo_name: NOAA-PMEL/EDD-DepthSimulator


text:
#!/usr/bin/python

import matplotlib
import matplotlib.pyplot as plt
from matplotlib import gridspec
import numpy as np
import datetime
import json
import pandas
import matplotlib
import argparse
import os
matplotlib.style.use('ggplot')

# sacrificial plot for single legend
matplotlib.rcParams['figure.figsize'] = 1, 1
randFig = plt.figure()
randAx = plt.subplot()
randAx.plot(0, 0, color='green', label='good', marker='+')
randAx.plot(0, 0, color='red', label='failed', marker='x')
randAx.plot(0, 0, color='black', label='sent', linestyle='--')
randAx.plot(0, 0, color='green', label='50th quantile')
randAx.plot(0, 0, color='orange', label='90th quantile')
randAx.plot(0, 0, color='red', label='99th quantile')
handles, labels = randAx.get_legend_handles_labels()

# big ol' plotting method
def plot_section(all_data, title, outputPath):
    # group calls by the endpoint/method
    actions = all_data.groupby('action')
    h = len(actions.groups.keys())
    matplotlib.rcParams['figure.figsize'] = 20, 3 * h

    fig = plt.figure()
    fig.legend(handles, labels, ncol=6, fontsize=16, framealpha=0, loc='upper center')
    if title is not None:
        fig.suptitle(title, fontsize=20, y=0.93)
    gs = gridspec.GridSpec(h, 3)

    # figure out left and right datetime bounds
    started = all_data['sent'].min()
    stopped = all_data['finished'].max()

    i = 0
    # plot one row of charts for each endpoint/method combination
    for section in actions.groups.keys():
        # setup the tree charts
        ax = fig.add_subplot(gs[i, 0])
        ax.set_title(section)
        ax.set_xlim(started, stopped)
        ax2 = fig.add_subplot(gs[i, 2])
        ax2.set_xlim(started, stopped)
        ax3 = fig.add_subplot(gs[i, 1])
        ax3.set_xlim(started, stopped)

        # find the maximum y value and set it across all three charts
        calls = actions.get_group(section)
        tookMax = calls['took'].max()
        ax.set_ylim(0, tookMax+tookMax*0.1)
        ax2.set_ylim(0, tookMax+tookMax*0.1)
        ax3.set_ylim(0, tookMax+tookMax*0.1)

        groups = calls.groupby('type')
        if groups.groups.get('error', False) is not False:
            bad = groups.get_group('error')
            ax.plot_date(bad['finished'], bad['took'], color='red', marker='x', label='error')

            bad_rate = bad.set_index('finished')
            bad_rate['rate'] = [0] * len(bad_rate.index)
            bad_rate = bad_rate.resample('5S').count()
            bad_rate['rate'] = bad_rate['rate'].divide(5)
            rateMax = bad_rate['rate'].max()
            ax2.plot_date(bad_rate.index, bad_rate['rate'], linestyle='-', marker='', color='red', label='error')
        if groups.groups.get('good', False) is not False:
            good = groups.get_group('good')
            ax.plot_date(good['finished'], good['took'], color='green', marker='+', label='good')

            good_rate = good.set_index('finished')
            good_rate['rate'] = [0] * len(good_rate.index)
            good_rate = good_rate.resample('5S').count()
            good_rate['rate'] = good_rate['rate'].divide(5)
            rateMax = good_rate['rate'].max()
            ax2.plot_date(good_rate.index, good_rate['rate'], linestyle='-', marker='', color='green', label='good')
        ax.set_ylabel('Latency (ms)')

        # calculate the request rate
        sent_rate = pandas.DataFrame(calls['sent'])
        sent_rate = sent_rate.set_index('sent')
        sent_rate['rate'] = [0] * len(sent_rate.index)
        sent_rate = sent_rate.resample('5S').count()
        sent_rate['rate'] = sent_rate['rate'].divide(5)
        if sent_rate['rate'].max() > rateMax:
            rateMax = sent_rate['rate'].max()
        ax2.plot_date(sent_rate.index, sent_rate['rate'], linestyle='--', marker='', color='black', label='sent')
        ax2.set_ylim(0, rateMax+rateMax*0.1)
        ax2.set_ylabel('Rate (per second)')

        # calculate and plot latency quantiles
        calls = calls.set_index('finished')
        calls = calls.sort_index()
        quan = pandas.DataFrame(calls['took'])
        for q, c in [[.5, 'green'], [.9, 'orange'], [.99, 'red']]:
            quanN = quan.rolling(500, center=True).quantile(q)
            ax3.plot(quanN['took'].index, quanN['took'], color=c)

        ax3.set_ylabel('Latency quantiles (ms)')

        i += 1

    # format x axes
    for ax in fig.axes:
        matplotlib.pyplot.sca(ax)
        plt.xticks(rotation=30, ha='right')
        majorFormatter = matplotlib.dates.DateFormatter('%H:%M:%S')
        ax.xaxis.set_major_formatter(majorFormatter)

    # save image
    gs.update(wspace=0.275, hspace=0.5)
    fig.savefig(outputPath, bbox_inches='tight')

# and the main event
parser = argparse.ArgumentParser()
parser.add_argument('chartData', type=str, help='Path to file containing JSON chart output from load-generator')
parser.add_argument('--output', type=str, help='Path to save output to', default='latency-chart.png')
parser.add_argument('--title', type=str, help='Chart title')
args = parser.parse_args()

with open(args.chartData) as data_file:
    stuff = []
    for l in data_file.readlines():
        stuff.append(json.loads(l))

df = pandas.DataFrame(stuff)
df['finished'] = pandas.to_datetime(df['finished']).astype(datetime.datetime)
df['sent'] = pandas.to_datetime(df['sent']).astype(datetime.datetime)
df['took'] = df['took'].divide(1000000)

plot_section(df, args.title, args.output)

max_stars_repo_name: CyberEncryptCA/boulder

max_issues_repo_name: CyberEncryptCA/boulder

max_forks_repo_name: CyberEncryptCA/boulder


text:
from collections import OrderedDict

import numpy as np
from astropy.table import Table, vstack, Column

Zsun = 0.0152  # this value from CMD website
Zmin = 0.0001
Zmax = 0.07
logtmax = 10.13
logtmin = 1.

class Isochrone(Table):
    """ This class accepts astropy.table.Table objects with 4 columns:
            _lgmini:    log10(M_ini/M_sun)
            _feh:       [Fe/H]
            _lgage:     log10(age/yr)
            _eep:       EEP number

    """

    meta = OrderedDict()
    model = None

    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)

    def add_column_safely(self, col):
        if col.name in self.colnames:
            self[col.name] = col.data
        else:
            self.add_column(col)

    def post_proc(self, model=None):
        """ post-processing of isochrone

        Parameters
        ----------
        model:
            {None | "parsec" | "mist"}

        """
        if model is None:
            pass

        elif model == "parsec":
            _lgmini = Column(np.log10(self["Mini"]), name="_lgmini")
            _fehini = Column(np.log10(self["Zini"]/Zsun), name="_fehini") #Z or M ???
            _lgage  = Column(np.log10(self["Age"]), name="_lgage")
            for col in [_lgmini, _lgage, _fehini]:
                self.add_column_safely(col)
            self.model = model

        elif model == "mist":
            _lgmini = Column(np.log10(self["initial_mass"]), name="_lgmini")
            _fehini = Column(self["Fe/H]_init"], name="_fehini")
            _lgage = Column(np.log10(self["log10_isochrone_age_yr"]), name="_lgage")
            for col in [_lgmini, _lgage, _fehini]:
                self.add_column_safely(col)
            self.model = model

        else:
            raise ValueError("@Isochrone: model *{}* unknown!".format(model))

    @staticmethod
    def vstack(*args, **kwargs):
        return vstack(*args, **kwargs)

    @staticmethod
    def from_table(tbl):
        return Isochrone(tbl.columns, meta=tbl.meta)

    def to_table(self):
        return Table(self.columns)

    # info of isochrone
    @property
    def Neep(self):
        # N_points of isochrone
        return len(self)

    @property
    def Ncol(self):
        # N_columns of isochrone
        return len(self.colnames)

    # parameters of this isochrone
    @property
    def lgage(self):
        self._age = self["_lgage"].data[0]
        return self._age

    @property
    def feh(self):
        return self["_feh"].data[0]

    @property
    def maxlgmass(self):
        return np.max(self["_lgmass"])

    @property
    def minlgmass(self):
        return np.min(self["_lgmass"])

    @property
    def maxeep(self):
        return np.max(self["_eep"])

    @property
    def mineep(self):
        return np.min(self["_eep"])

    @staticmethod
    def combine(isolists):
        result = []
        for isolist in isolists:
            for iso in isolist:
                result.append(iso)
        return result

max_stars_repo_name: hypergravity/berliner

max_issues_repo_name: hypergravity/berliner

max_forks_repo_name: hypergravity/berliner

reason: import numpy,from astropy


text:
import numpy as np
import os

from chainercv.chainer_experimental.datasets.sliceable import GetterDataset
from chainercv.utils import read_image

linemod_object_diameters = {
    'ape': 0.103,
    'benchvise': 0.286908,
    'cam': 0.173,
    'can': 0.202,
    'cat': 0.155,
    'driller': 0.262,
    'duck': 0.109,
    'eggbox': 0.176364,
    'glue': 0.176,
    'holepuncher': 0.162,
    'iron': 0.303153,
    'lamp': 0.285155,
    'phone': 0.213}

class LinemodDataset(GetterDataset):

    def __init__(self, base_dir, obj_name='ape', split='train',
                 return_msk=False):
        super(LinemodDataset, self).__init__()
        split_path = os.path.join(
            base_dir, 'LINEMOD', obj_name, '{}.txt'.format(split))
        self.base_dir = base_dir
        with open(split_path, 'r') as f:
            self.img_paths = f.readlines()

        self.add_getter(('img', 'point', 'label'), self._get_example)
        if return_msk:
            self.add_getter('msk', self._get_msk)

    def __len__(self):
        return len(self.img_paths)

    def _get_example(self, i):
        img_path = os.path.join(self.base_dir, self.img_paths[i].rstrip())
        img = read_image(img_path)

        anno_path = img_path.replace(
            'images', 'labels').replace(
                'JPEGImages', 'labels').replace(
                    '.jpg', '.txt').replace('.png', '.txt')

        anno = np.zeros(50*21)
        if os.path.getsize(anno_path):
            _, H, W = img.shape
            tmp = read_truths_args(anno_path, 8.0/W)
            size = tmp.size
            if size > 50*21:
                anno = tmp[0:50*21]
            elif size > 0:
                anno[0:size] = tmp
        anno = anno.reshape(-1, 21)
        anno = anno[:truths_length(anno)]
        point = anno[:, 1:19].reshape(-1, 9, 2).astype(np.float32)
        point[:, :, 0] *= W
        point[:, :, 1] *= H
        label = anno[:, 0].astype(np.int32)
        return img, point, label

    def _get_msk(self, i):
        img_path = os.path.join(self.base_dir, self.img_paths[i].rstrip())
        mskpath = img_path.replace('JPEGImages', 'mask').replace(
            '/00', '/').replace('.jpg', '.png')
        msk = read_image(mskpath, color=False)[0]
        return msk > 0

def truths_length(truths):
    for i in range(50):
        if truths[i][1] == 0:
            return i

def read_truths(lab_path):
    if os.path.getsize(lab_path):
        truths = np.loadtxt(lab_path)
        # to avoid single truth problem
        truths = truths.reshape(truths.size//21, 21)
        return truths
    else:
        return np.array([])

def read_truths_args(lab_path, min_box_scale):
    truths = read_truths(lab_path)
    new_truths = []
    for i in range(truths.shape[0]):
        new_truths.append(
            [truths[i][0], truths[i][1], truths[i][2],
             truths[i][3], truths[i][4], truths[i][5],
             truths[i][6], truths[i][7], truths[i][8],
             truths[i][9], truths[i][10], truths[i][11],
             truths[i][12], truths[i][13], truths[i][14],
             truths[i][15], truths[i][16], truths[i][17],
             truths[i][18]])
    return np.array(new_truths)

max_stars_repo_name: take-cheeze/models

max_issues_repo_name: take-cheeze/models

max_forks_repo_name: take-cheeze/models


text:
# Copyright 2020 Tensorforce Team. All Rights Reserved.
#
# 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 numpy as np

from tensorforce.environments import Environment

class ArcadeLearningEnvironment(Environment):
    """
    [Arcade Learning Environment](https://github.com/mgbellemare/Arcade-Learning-Environment)
    adapter (specification key: `ale`, `arcade_learning_environment`).

    May require:
    ```bash
    sudo apt-get install libsdl1.2-dev libsdl-gfx1.2-dev libsdl-image1.2-dev cmake
    ```

    Args:
        level (string): ALE rom file
            (required).
        loss_of_life_termination: Signals a terminal state on loss of life
            (default: false).
        loss_of_life_reward (float): Reward/Penalty on loss of life (negative values are a penalty)
            (default: 0.0).
        repeat_action_probability (float): Repeats last action with given probability
            (default: 0.0).
        visualize (bool): Whether to visualize interaction
            (default: false).
        frame_skip (int > 0): Number of times to repeat an action without observing
            (default: 1).
        seed (int): Random seed
            (default: none).
    """

    def __init__(
        self, level, life_loss_terminal=False, life_loss_punishment=0.0,
        repeat_action_probability=0.0, visualize=False, frame_skip=1, seed=None
    ):
        super().__init__()

        from ale_py import ALEInterface

        self.environment = ALEInterface()
        self.rom_file = level

        self.life_loss_terminal = life_loss_terminal
        self.life_loss_punishment = life_loss_punishment

        self.environment.setFloat(b'repeat_action_probability', repeat_action_probability)
        self.environment.setBool(b'display_screen', visualize)
        self.environment.setInt(b'frame_skip', frame_skip)
        if seed is not None:
            self.environment.setInt(b'random_seed', seed)

        # All set commands must be done before loading the ROM.
        self.environment.loadROM(rom_file=self.rom_file.encode())
        self.available_actions = tuple(self.environment.getLegalActionSet())

        # Full list of actions:
        # No-Op, Fire, Up, Right, Left, Down, Up Right, Up Left, Down Right, Down Left, Up Fire,
        # Right Fire, Left Fire, Down Fire, Up Right Fire, Up Left Fire, Down Right Fire, Down Left
        # Fire

    def __str__(self):
        return super().__str__() + '({})'.format(self.rom_file)

    def states(self):
        width, height = self.environment.getScreenDims()
        return dict(type='float', shape=(height, width, 3), min_value=0.0, max_value=1.0)

    def actions(self):
        return dict(type='int', num_values=len(self.available_actions))

    def close(self):
        # Causes "Segmentation fault (core dumped)"
        # self.environment.__del__()
        self.environment = None

    def get_states(self):
        screen = np.copy(self.environment.getScreenRGB(screen_data=self.screen))
        screen = screen.astype(dtype=np.float32) / 255.0
        return screen

    def reset(self):
        self.environment.reset_game()
        width, height = self.environment.getScreenDims()
        self.screen = np.empty((height, width, 3), dtype=np.uint8)
        self.lives = self.environment.lives()
        return self.get_states()

    def execute(self, actions):
        reward = self.environment.act(action=self.available_actions[actions])
        terminal = self.environment.game_over()
        states = self.get_states()

        next_lives = self.environment.lives()
        if next_lives < self.lives:
            if self.life_loss_terminal:
                terminal = True
            elif self.life_loss_punishment > 0.0:
                reward -= self.life_loss_punishment
            self.lives = next_lives

        return states, terminal, reward

max_stars_repo_name: DLPerf/tensorforce

max_issues_repo_name: DLPerf/tensorforce

max_forks_repo_name: DLPerf/tensorforce


text:
from typing import *
import logging
import warnings
import numpy as np
from .defaults import load_config, ReconstructionConfig, connect_to_data, DataConnection
from .core import build_Design_Matrix, prepare_design_masked, prepare_design_symmetry, prepare_design_symmetry_masked, prepare_observations, prepare_observations_symmetry
from .optimization import ReconstructorFastScipyNB

class Tomographer:
    """Main object to work with tomographic reconstruction problems

    The analysis steps are:

    #. First load the defaults

        * tg = Tomographer()
        * tg.load_cfg("config_file.hdf5")  # or set this manually
        * tg.prepare_design()

    #. Then load the data

        * tg.connect_data("data.hdf5")

    #. Finally, reconstruct one gene

        * tg.reconstruct(gene="Lmx1a")
    """

    def __init__(self) -> None:
        self.cfg = ReconstructionConfig()  # type: ReconstructionConfig
        self.data = None  # type: DataConnection
        self.reconstructor = None  # type: Reconstructor
        self.ready_flag = False

    def load_cfg(self, cfg: Union[str, ReconstructionConfig]) -> None:
        """Load defaults from file
        Args
        ----
        defaults: filename (str) or ReconstructionDefaults (object)
        """
        if type(cfg) == str:
            self.cfg = load_config(cfg)
        else:
            assert type(cfg) == ReconstructionConfig
            self.cfg = cfg

    def connect_data(self, datafile: str) -> None:
        """Load hdf5 dataset
        """
        self.data = connect_to_data(datafile)
        self.data.bind(self.cfg)

    def prepare_design(self, symmetry: bool=None, masked_formulation: bool=None, kfolds: int = 5) -> None:
        if symmetry is not None:
            self.cfg.symmetry = symmetry
        if masked_formulation is not None:
            self.cfg.masked_formulation = masked_formulation

        D, proj_len = build_Design_Matrix(self.cfg.angles_values,
                                          self.cfg.widths,
                                          self.cfg.mask_g,
                                          self.cfg.mask_thrs,
                                          notation_inverse=True,
                                          return_projlen=True)
        self.cfg.proj_len = proj_len

        if self.cfg.symmetry:
            self.prepare = prepare_observations_symmetry
            if self.cfg.masked_formulation:
                self.cfg.A = prepare_design_symmetry_masked(D, self.cfg.mask_bw)
                self.reconstructor = ReconstructorFastScipyNB(config=self.cfg)
            else:
                self.cfg.A = prepare_design_symmetry(D)
                self.reconstructor = Reconstructor(config=self.cfg)
        else:
            self.prepare = prepare_observations
            if self.cfg.masked_formulation:
                self.cfg.A = prepare_design_masked(D, self.cfg.mask_bw)
                self.reconstructor = ReconstructorFastScipyNB(config=self.cfg, kfolds=kfolds)
            else:
                self.cfg.A = D
                self.reconstructor = Reconstructor(config=self.cfg)

        self.ready_flag = True

    def reconstruct(self, gene: str, alpha_beta: Union[str, Tuple[float, float]]="auto",return_heuristic=False,
                    warm_start: bool=False,
                    crossval_kwargs: Dict[str, Any]={}, nb_r = 0.3, kfolds=3) -> np.ndarray:
        """Perform reconstruction of a gene

        Args
        ----
        gene: str
            gene name as in datafile

        alpha_beta: str or tuple
            options are:
                "auto": it will try to use the alpha and beta parameters present in the data file, if none is provided it will perform crossvalidation
                "crossvalidation": it will force crossvalidation
                (alpha, beta): passing a tuple will perform the reconstruction using this parameters

        return_heuristic: bool
            whether to output the reconstruction along with the heuristic score calculated by 2log(sparsity) + log(expression)
            returns tuple(reconstruction, heuristic)

        warm_start:
            whether to use the previous reconstruction result as warm start.
            NOTE what the previous reconstruction is might change depending the alpha_beta parameter

        Return
        ------
        reconstructed: np.ndarray (2d)
            the reconstructed signal already reshaped

        NOTE you can also access the recosntruction parameters using the attribute `.reconstructor`
        """
        if not self.ready_flag:
            warnings.warn("prepare_design should be called before attempting reconstruction otherwise it is called automatically with default arguments!", UserWarning)
            self.prepare_design(kfolds=kfolds)

        logging.info("Preparing Reconstruction")
        b = self.prepare(projections=self.data[gene], xs=self.data.xs,
                         first_points=self.cfg.first_points, projs_len=self.cfg.proj_len,
                         interpolation="linear")

        # Normalize b taking some precaoution against high outlier and against division by 0
        # More naive form: b / b.max()
        self.unnormalized_b = np.copy(b)

        if type(self.reconstructor) is ReconstructorFastScipyNB:
            self.norm_factor = 1
            self.reconstructor.r = nb_r
        else:
            self.norm_factor = np.clip(np.percentile(b, 98), np.max(b) / 3., np.inf)
            b = self.unnormalized_b / self.norm_factor

        self.reconstructor.norm_factor = self.norm_factor

        if alpha_beta == "auto":
            try:
                _alpha, _beta = self.data.alphas_betas[gene]
                logging.info("alpha and beta found in the data file")
                self.reconstructor.change_par(alpha=_alpha, beta=_beta)
            except KeyError:
                logging.debug("alpha and beta not found in the data file, performing crossvalidation")
                self.reconstructor.optimize(b, **crossval_kwargs)
        elif alpha_beta == "crossvalidation":
            self.reconstructor.optimize(b, **crossval_kwargs)
        elif isinstance(alpha_beta, tuple) or isinstance(alpha_beta, list) or isinstance(alpha_beta, np.array):
            _alpha, _beta = alpha_beta
            self.reconstructor.change_par(alpha=_alpha, beta=_beta)
        logging.info("Reconstructing")
        result = self.reconstructor.fit_predict(b, warm_start=warm_start)
        logging.info("Finished Reconstructing %s" % gene)

        if return_heuristic:
            sparsity_measure = 1. / np.sum(result > 0.4 * result)
            expression = np.mean(b[b > 0])
            reliability_heuristic = 2 * np.log(sparsity_measure) + np.log(expression)
            return result, reliability_heuristic

        return result

class TomographerDebug(Tomographer):
    def __init__(self) -> None:
        super().__init__()

    def prepare_design(self, symmetry: bool=None, masked_formulation: bool=None) -> None:
        if symmetry is not None:
            self.cfg.symmetry = symmetry
        if masked_formulation is not None:
            self.cfg.masked_formulation = masked_formulation

        D, proj_len = build_Design_Matrix(self.cfg.angles_values,
                                          self.cfg.widths,
                                          self.cfg.mask_g,
                                          self.cfg.mask_thrs,
                                          notation_inverse=True,
                                          return_projlen=True)
        self.cfg.proj_len = proj_len

        if self.cfg.symmetry:
            self.prepare = prepare_observations_symmetry
            if self.cfg.masked_formulation:
                self.cfg.A = prepare_design_symmetry_masked(D, self.cfg.mask_bw)
                self.reconstructor = ReconstructorFastTest(config=self.cfg)
            else:
                self.cfg.A = prepare_design_symmetry(D)
                self.reconstructor = Reconstructor(config=self.cfg)
        else:
            self.prepare = prepare_observations
            if self.cfg.masked_formulation:
                self.cfg.A = prepare_design_masked(D, self.cfg.mask_bw)
                self.reconstructor = ReconstructorFastTest(config=self.cfg)
            else:
                self.cfg.A = D
                self.reconstructor = Reconstructor(config=self.cfg)

        self.ready_flag = True

    def reconstruct(self, gene: str, alpha_beta: Union[str, Tuple[float, float]]="auto",
                    warm_start: bool=False,
                    crossval_kwargs: Dict[str, Any]={}) -> np.ndarray:
        """Perform reconstruction of a gene

        Args
        ----
        gene: str
            gene name as in datafile

        alpha_beta: str or tuple
            options are:
                "auto": it will try to use the alpha and beta parameters present in the data file, if none is provided it will perform crossvalidation
                "crossvalidation": it will force crossvalidation
                (alpha, beta): passing a tuple will perform the reconstruction using this parameters

        warm_start:
            whether to use the previous reconstruction result as warm start.
            NOTE what the previous reconstruction is might change depending the alpha_beta parameter

        Return
        ------
        reconstructed: np.ndarray (2d)
            the reconstructed signal already reshaped

        NOTE you can also access the recosntruction parameters using the attribute `.reconstructor`
        """
        if not self.ready_flag:
            warnings.warn("prepare_design should be called before attempting reconstruction otherwise it is called automatically with default arguments!", UserWarning)
            self.prepare_design()

        logging.info("Preparing Reconstruction")
        b = self.prepare(projections=self.data[gene], xs=self.data.xs,
                         first_points=self.cfg.first_points, projs_len=self.cfg.proj_len,
                         interpolation="linear")

        # Normalize b taking some precaoution against high outlier and against division by 0
        # More naive form: b / b.max()
        self.unnormalized_b = np.copy(b)
        self.norm_factor = np.maximum(np.percentile(b, 98), np.max(b) / 3.)
        # hopefully this is the count data as i comment ourt the self.norm_factor lines
        b = self.unnormalized_b / self.norm_factor
        self.reconstructor.norm_factor = self.norm_factor

        if alpha_beta == "auto":
            try:
                _alpha, _beta = self.data.alphas_betas[gene]
                logging.info("alpha and beta found in the data file")
                self.reconstructor.change_par(alpha=_alpha, beta=_beta)
            except KeyError:
                logging.debug("alpha and beta not found in the data file, performing crossvalidation")
                self.reconstructor.optimize(b, **crossval_kwargs)
        elif alpha_beta == "crossvalidation":
            self.reconstructor.optimize(b, **crossval_kwargs)
        elif isinstance(alpha_beta, tuple) or isinstance(alpha_beta, list) or isinstance(alpha_beta, np.array):
            _alpha, _beta = alpha_beta
            self.reconstructor.change_par(alpha=_alpha, beta=_beta)
        logging.info("Reconstructing")
        result = self.reconstructor.fit_predict(b, warm_start=warm_start)
        logging.info("Finished Reconstructing %s" % gene)
        return result

class TomographerManualDesignMatrix(Tomographer):
    def __init__(self, A: np.ndarray, projection_lengths: list) -> None:
        super().__init__()
        self.A_manual = A
        self.projection_lengths = projection_lengths

       # self.reconstructor = ReconstructorFastScipyNB()

    def prepare_design(self, symmetry: bool=None, masked_formulation: bool=None) -> None:
        if symmetry is not None:
            self.cfg.symmetry = symmetry
        if masked_formulation is not None:
            self.cfg.masked_formulation = masked_formulation

        self.cfg.A = self.A_manual
        D = self.A_manual
        self.cfg.proj_len =  self.projection_lengths #proj_len

        if self.cfg.symmetry:
            self.prepare = prepare_observations_symmetry
            if self.cfg.masked_formulation:
                #self.cfg.A = prepare_design_symmetry_masked(D, self.cfg.mask_bw)
                self.reconstructor = ReconstructorFastScipyNB(config=self.cfg)
            else:
                #self.cfg.A = prepare_design_symmetry(D)
                self.reconstructor = ReconstructorFastScipyNB(config=self.cfg)
        else:
            self.prepare = prepare_observations
            if self.cfg.masked_formulation:
                #self.cfg.A = prepare_design_masked(D, self.cfg.mask_bw)
                self.reconstructor = ReconstructorFastScipyNB(config=self.cfg)
            else:
                #self.cfg.A = D
                self.reconstructor = ReconstructorFastScipyNB(config=self.cfg)

        self.cfg.A = self.A_manual
        self.ready_flag = True

    def reconstruct(self, b: np.ndarray, alpha_beta: Union[str, Tuple[float, float]]="auto",
                    warm_start: bool=False,
                    crossval_kwargs: Dict[str, Any]={}) -> np.ndarray:
        """Perform reconstruction of a gene

        Args
        ----
        gene: str
            gene name as in datafile

        alpha_beta: str or tuple
            options are:
                "auto": it will try to use the alpha and beta parameters present in the data file, if none is provided it will perform crossvalidation
                "crossvalidation": it will force crossvalidation
                (alpha, beta): passing a tuple will perform the reconstruction using this parameters

        warm_start:
            whether to use the previous reconstruction result as warm start.
            NOTE what the previous reconstruction is might change depending the alpha_beta parameter

        Return
        ------
        reconstructed: np.ndarray (2d)
            the reconstructed signal already reshaped

        NOTE you can also access the recosntruction parameters using the attribute `.reconstructor`
        """

        #b = self.observations[gene] # select observation vector

        # Normalize b taking some precaoution against high outlier and against division by 0
        # More naive form: b / b.max()
        self.unnormalized_b = np.copy(b)
        self.norm_factor = np.maximum(np.percentile(b, 98), np.max(b) / 3.)
        # hopefully this is the count data as i comment ourt the self.norm_factor lines
        b = self.unnormalized_b / self.norm_factor

        #self.reconstructor.norm_factor = self.norm_factor

        if alpha_beta == "auto":
            try:
                _alpha, _beta = self.data.alphas_betas[gene]
                logging.info("alpha and beta found in the data file")
                self.reconstructor.change_par(alpha=_alpha, beta=_beta)
            except KeyError:
                logging.debug("alpha and beta not found in the data file, performing crossvalidation")
                self.reconstructor.optimize(b, **crossval_kwargs)
        elif alpha_beta == "crossvalidation":
            self.reconstructor.optimize(b, **crossval_kwargs)
        elif isinstance(alpha_beta, tuple) or isinstance(alpha_beta, list) or isinstance(alpha_beta, np.array):
            _alpha, _beta = alpha_beta
            self.reconstructor.change_par(alpha=_alpha, beta=_beta)
        logging.info("Reconstructing")
        result = self.reconstructor.fit_predict(b, warm_start=warm_start)
        logging.info("Finished Reconstructing %")
        return result

max_stars_repo_name: lamanno-epfl/tomographer

max_issues_repo_name: lamanno-epfl/tomographer

max_forks_repo_name: lamanno-epfl/tomographer


text:
# Author: Chuan Wang
import json
import pandas as pd
from rnn.biLSTM_inference import biLSTM_inference
from torch import from_numpy
from numpy import load, copy
from flask import Flask, jsonify, request
from flask_cors import CORS
from rnn.parameter import FEATURES
app = Flask(__name__)
CORS(app)

time = '20200115-194901'
best_epoch = 46
best_accuracy = 0.96
filepath = './data/mimic/'

dataset = load(filepath + 'unique_data/uniqDataLabelsIds')

def jsonData2numpy():
    return None

def numpyData2Json(input, featureIdx):
    T = input.shape[0]
    features = input[:, featureIdx]
    featuresDf = pd.DataFrame(
        features, index=['T' + str(t + 1) for t in range(T)])

    featuresDf.columns = featureIdx
    featuresDf['time'] = featuresDf.index

    return json.loads(featuresDf.round(2).to_json(orient='records'))

@app.route("/", methods=["POST", "GET"])
def index():
    return jsonify({'abc': 'ddd', 'bbd': 'ccc'})

@app.route("/counterfactual", methods=["POST"])
def counterfactual():
    d = request.get_json()
    instanceId = int(d['instanceId'])

    arr = load(filepath + 'counterfactuals/all_counterfactuals20')

    # get counterfactuals for selected data
    jsonObjArr = []
    for ctfct in arr[instanceId]:
        jsonObjArr.append({
            'idx': ctfct[0],
            'similarity': 1 - ctfct[1]
        })

    return jsonify(jsonObjArr)

@app.route("/pdplot", methods=["POST"])
def pdplot():
    d = request.get_json()
    instanceId = int(d['instanceId'])

    # with open(filepath + 'pdplot/pdp-negs-30-46-47_clusters_ci_0.999.json') as jsonfile:
    with open(filepath + 'pdplot/pdp-negs-30-46-47_clusters_ci_new.json') as jsonfile:
        jsonObjArr = json.load(jsonfile)

    return jsonify(jsonObjArr)

@app.route("/predict", methods=["POST"])
def predict():
    d = request.get_json()
    instanceId = int(d['instanceId'])
    featureIdx = d['featureIdx']

    print('\n ---- predict instanceId', instanceId)
    # print('\n ---', d['updatedData'])

    # print(dataset['data'].shape)  # (14165, 48, 37)
    # print(dataset['labels'].shape)  # (14165, 2)
    # print(dataset['ids'].shape)  # (14165,)
    data = dataset['data']
    labels = dataset['labels']

    selectedDatum = data[instanceId]
    selectedLabel = labels[instanceId]

    # get updated input from frontend

    # input = from_numpy(selectedData[instanceId])
    input = from_numpy(selectedDatum)
    # counterfactual = from_numpy(data[9009]) # 6712
    # counterfactual = from_numpy(data[5875])  # 6712
    # counterfactual = from_numpy(data[5947])  # 4895
    counterfactual = from_numpy(data[12457])  # 4895
    # print(numpyData2Json(input.numpy(), featureIdx))

    model = biLSTM_inference(filepath, time, best_epoch, best_accuracy)
    result = model.predict(input)

    # updated input
    updatedInput = copy(input)
    # d['data2predict'] is empty [] at load and when inputing text id
    # d['data2predict'] has values when updating barchart values
    updatedDataDf = pd.DataFrame(d['data2predict'])
    for column in updatedDataDf.head():
        if (column != "time"):
            updatedInput[:, int(column)] = updatedDataDf[column]
    # prediction result for updated input
    updatedResult = model.predict(from_numpy(updatedInput))

    print('prediction result\n', result)
    print('updated prediction result\n', updatedResult)

    obj = {
        'xName': 'time',
        'yName': [str(i) for i in featureIdx],
        'yLabel': [FEATURES[i] for i in featureIdx],
        'instanceCnt': 7,
        'result':
            [{
                'class': 'dead',
                'original': round(result[0][0], 4),
                'predict':round(updatedResult[0][0], 4)
            },
            {
                'class': 'alive',
                'original': round(result[0][1], 4),
                'predict':round(updatedResult[0][1], 4)
            }],
        'original': numpyData2Json(input.numpy(), featureIdx),
        'counterfactual': numpyData2Json(counterfactual.numpy(), featureIdx),
        'updated': numpyData2Json(updatedInput, featureIdx)
    }
    return jsonify(obj)

@app.route("/tsne", methods=["POST"])
def tsne():
    tsneFile = './data/mimic/tsne4vis.json'
    with open(tsneFile) as json_file:
        tsneData = json.load(json_file)
        return jsonify(tsneData)

if __name__ == "__main__":
    app.run(debug=True)

max_stars_repo_name: nauhc/hyppersteer

max_issues_repo_name: nauhc/hyppersteer

max_forks_repo_name: nauhc/hyppersteer

reason: from numpy


text:
from pylayers.antprop.antenna import *
from pylayers.antprop.antvsh import *
import matplotlib.pylab as plt
from numpy import *
import pdb
"""
This test :

    1 : loads a measured antenna
    2 : applies an electrical delay obtained from data with getdelay method
    3 : evaluate the antenna vsh coefficient with a downsampling factor of 2
    4 : display the 16 first
"""
filename = 'S1R1.mat'

A = Antenna(filename,directory='ant/UWBAN/Matfile')

#plot(freq,angle(A.Ftheta[:,maxPowerInd[1],maxPowerInd[2]]*exp(2j*pi*freq.reshape(len(freq))*electricalDelay)))
freq = A.fa.reshape(104,1,1)
delayCandidates = arange(-10,10,0.001)
electricalDelay = A.getdelay(freq,delayCandidates)
disp('Electrical Delay = ' + str(electricalDelay)+' ns')

A.Ftheta = A.Ftheta*exp(2*1j*pi*freq*electricalDelay)
A.Fphi   = A.Fphi*exp(2*1j*pi*freq*electricalDelay)

dsf = 2
#
# Calculate Vector Spherical Harmonics
#
A = vsh(A,dsf)
A.C.s1tos2(15)
EBr,EBi,ECr,ECi= A.Fsynth2s()
plt.figure()
plt.subplot(221)
plt.plot(EBr)
plt.subplot(222)
plt.plot(EBi)
plt.subplot(223)
plt.plot(ECr)
plt.subplot(224)
plt.plot(ECi)
plt.subplot(224)

max_stars_repo_name: usmanwardag/pylayers

max_issues_repo_name: usmanwardag/pylayers

max_forks_repo_name: usmanwardag/pylayers


text:
import numpy as np
from Common import Activation, Diff , SGD , LossFun
from TwoLayerNetmy import Twolayernetmy
import os
import sys
from matplotlib import pyplot as plt

sys.path.append(os.pardir)
sys.path.append(os.pardir+"dataset")

from dataset.mnist import load_mnist

(x_train,y_train) , (x_test,y_test) = load_mnist(True,True)

trainLossList = []

# Setting
iter = 10
trainSize = x_train.shape[0]
batchSize = 100
lr = 0.1

inputSize = 784
hiddenSize = 50
outputSize = 10

net = Twolayernetmy(
            inputSize ,
            hiddenSize ,
            outputSize ,
            isOneHot=  False)

for i in range(iter):
    print( "Iter : {0}".format(i) )
    batchMask = np.random.choice(trainSize,batchSize)
    xBatch = x_train[batchMask]
    tBatch = y_train[batchMask]

    #Gradient Calculate
    grad = net.NGradient(xBatch,tBatch)

    for key in ('W1','b1','W2','b2'):
        # Update Parameter with SGD
        net.params[key] -= lr * grad[key]

    loss = net.Loss(xBatch,tBatch)
    trainLossList.append(loss)

xlabels = range(0,len(trainLossList))
plt.plot(xlabels,trainLossList  )

plt.show()