text:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.tri import (
    Triangulation, UniformTriRefiner, CubicTriInterpolator)
from copy import deepcopy

FONTSIZE = 12

def plot_dynamics(history_pi_0, history_pi_1, pi_alpha_gradient_history, pi_beta_gradient_history, title=''):
    cmap = plt.get_cmap('viridis')
    colors = range(len(history_pi_1))
    fig = plt.figure(figsize=(6, 5))
    ax = fig.add_subplot(111)

    scatter = ax.scatter(history_pi_0, history_pi_1, c=colors, s=1)
    ax.scatter(0.5, 0.5, c='r', s=10., marker='*')
    colorbar = fig.colorbar(scatter, ax=ax)
    colorbar.set_label('Iterations', rotation=270, fontsize=FONTSIZE)

    skip = slice(0, len(history_pi_0), 50)
    ax.quiver(history_pi_0[skip],
              history_pi_1[skip],
              pi_alpha_gradient_history[skip],
              pi_beta_gradient_history[skip],
              units='xy', scale=10., zorder=3, color='blue',
              width=0.007, headwidth=3., headlength=4.)

    ax.set_ylabel("Policy of Player 2", fontsize=FONTSIZE)
    ax.set_xlabel("Policy of Player 1", fontsize=FONTSIZE)
    ax.set_ylim(0, 1)
    ax.set_xlim(0, 1)
    ax.set_title(title, fontsize=FONTSIZE+8)
    plt.tight_layout()
    plt.show()

def IGA(pi_alpha,
        pi_beta,
        payoff_0,
        payoff_1,
        u_alpha,
        u_beta,
        iteration=1000,
        lr=0.01):
    pi_alpha_history = [pi_alpha]
    pi_beta_history = [pi_beta]
    pi_alpha_gradient_history = [0.]
    pi_beta_gradient_history = [0.]
    for i in range(iteration):
        pi_alpha_gradient = (pi_beta * u_alpha + payoff_0[(0, 1)] - payoff_0[(1, 1)])
        pi_beta_gradient = (pi_alpha * u_beta + payoff_1[(1, 0)] - payoff_1[(1, 1)])
        pi_alpha_gradient_history.append(pi_alpha_gradient)
        pi_beta_gradient_history.append(pi_beta_gradient)
        pi_alpha_next = pi_alpha + lr * pi_alpha_gradient
        pi_beta_next = pi_beta + lr * pi_beta_gradient
        pi_alpha = max(0., min(1., pi_alpha_next))
        pi_beta = max(0., min(1., pi_beta_next))
        pi_alpha_history.append(pi_alpha)
        pi_beta_history.append(pi_beta)
    return pi_alpha_history, \
           pi_beta_history, \
           pi_alpha_gradient_history, \
           pi_beta_gradient_history

def WoLF_IGA(pi_alpha,
             pi_beta,
             payoff_0,
             payoff_1,
             u_alpha,
             u_beta,
             iteration=1000,
             pi_alpha_nash=0.5,
             pi_beta_nash=0.5,
             lr_min=0.01,
             lr_max=0.04):
    pi_alpha_history = [pi_alpha]
    pi_beta_history = [pi_beta]
    pi_alpha_gradient_history = [0.]
    pi_beta_gradient_history = [0.]
    # V_nash_alpha = V(pi_alpha_nash, pi_beta_nash, payoff_0)
    # V_nash_beta = V(pi_alpha_nash, pi_beta_nash, payoff_1)
    # print(V_nash_alpha, V_nash_beta)
    for i in range(iteration):
        lr_alpha = lr_max
        lr_beta = lr_max
        if V(pi_alpha, pi_beta, payoff_0) > V(pi_alpha_nash, pi_beta, payoff_0):
            lr_alpha = lr_min
        if V(pi_alpha, pi_beta, payoff_1) > V(pi_alpha, pi_beta_nash, payoff_0):
            lr_beta = lr_min

        pi_alpha_gradient = (pi_beta * u_alpha + payoff_0[(0, 1)] - payoff_0[(1, 1)])
        pi_beta_gradient = (pi_alpha * u_beta + payoff_1[(1, 0)] - payoff_1[(1, 1)])
        pi_alpha_gradient_history.append(pi_alpha_gradient)
        pi_beta_gradient_history.append(pi_beta_gradient)
        pi_alpha_next = pi_alpha + lr_alpha * pi_alpha_gradient
        pi_beta_next = pi_beta + lr_beta * pi_beta_gradient
        pi_alpha = max(0., min(1., pi_alpha_next))
        pi_beta = max(0., min(1., pi_beta_next))
        pi_alpha_history.append(pi_alpha)
        pi_beta_history.append(pi_beta)
    return pi_alpha_history, \
           pi_beta_history, \
           pi_alpha_gradient_history, \
           pi_beta_gradient_history

def IGA_PP(pi_alpha,
           pi_beta,
           payoff_0,
           payoff_1,
           u_alpha,
           u_beta,
           iteration=10000,
           lr=0.01,
           gamma=0.01,
           single=False):
    pi_alpha_history = [pi_alpha]
    pi_beta_history = [pi_beta]
    pi_alpha_gradient_history = [0.]
    pi_beta_gradient_history = [0.]
    for i in range(iteration):
        pi_beta_pp = pi_beta + gamma * (pi_alpha * u_beta + payoff_1[(1, 0)] - payoff_1[(1, 1)])
        pi_alpha_gradient = (pi_beta_pp * u_alpha + payoff_0[(0, 1)] - payoff_0[(1, 1)])
        pi_alpha_next = pi_alpha + lr * pi_alpha_gradient
        if not single:
            pi_alpha_pp = pi_alpha + gamma * (pi_beta * u_alpha + payoff_0[(0, 1)] - payoff_0[(1, 1)])
            pi_beta_gradient = (pi_alpha_pp * u_beta + payoff_1[(1, 0)] - payoff_1[(1, 1)])
            pi_beta_next = pi_beta + lr * pi_beta_gradient
        else:
            pi_beta_gradient = (pi_alpha * u_beta + payoff_1[(1, 0)] - payoff_1[(1, 1)])
            pi_beta_next = pi_beta + lr * pi_beta_gradient
        pi_alpha_gradient_history.append(pi_alpha_gradient)
        pi_beta_gradient_history.append(pi_beta_gradient)

        pi_alpha = max(0., min(1., pi_alpha_next))
        pi_beta = max(0., min(1., pi_beta_next))
        pi_alpha_history.append(pi_alpha)
        pi_beta_history.append(pi_beta)
    return pi_alpha_history, \
           pi_beta_history, \
           pi_alpha_gradient_history, \
           pi_beta_gradient_history

def V(alpha, beta, payoff):
    u = payoff[(0, 0)] - payoff[(0, 1)] - payoff[(1, 0)] + payoff[(1, 1)]
    return alpha * beta * u + alpha * (payoff[(0, 1)] - payoff[(1, 1)]) + beta * (payoff[(1, 0)] - payoff[(1, 1)]) + payoff[(1, 1)]

if __name__ == '__main__':

    game_name = 'wolf_05_05'

    pi_alpha = 0.
    pi_beta = 0.9
    pi_alpha_nash = 0.5
    pi_beta_nash = 0.5
    payoff_0 = np.array([[0, 3], [1, 2]])
    payoff_1 = np.array([[3, 2], [0, 1]])
    # \alpha, \beta = 0, 0.9, nash is 0.5 0.5

    u_alpha = payoff_0[(0, 0)] - payoff_0[(0, 1)] - payoff_0[(1, 0)] + payoff_0[(1, 1)]
    u_beta = payoff_1[(0, 0)] - payoff_1[(0, 1)] - payoff_1[(1, 0)] + payoff_1[(1, 1)]
    print(u_alpha, u_beta)

    agent = 'WPL'

    if agent == 'IGA':
        pi_alpha_history, \
        pi_beta_history, \
        pi_alpha_gradient_history, \
        pi_beta_gradient_history = IGA(pi_alpha,
                                       pi_beta,
                                       payoff_0,
                                       payoff_1,
                                       u_alpha,
                                       u_beta)
    elif agent == 'WoLF-IGA':
        pi_alpha_history, \
        pi_beta_history, \
        pi_alpha_gradient_history, \
        pi_beta_gradient_history = WoLF_IGA(pi_alpha,
                                            pi_beta,
                                            payoff_0,
                                            payoff_1,
                                            u_alpha,
                                            u_beta)

    elif agent == 'IGA-PP':
        pi_alpha_history, \
        pi_beta_history, \
        pi_alpha_gradient_history, \
        pi_beta_gradient_history = IGA_PP(pi_alpha,
                                          pi_beta,
                                          payoff_0,
                                          payoff_1,
                                          u_alpha,
                                          u_beta,
                                          single=False)

    plot_dynamics(pi_alpha_history,
                  pi_beta_history,
                  pi_alpha_gradient_history,
                  pi_beta_gradient_history,
                  agent)
    print('Done')


text:
from sklearn.metrics import accuracy_score
from sklearn.svm import LinearSVC
import numpy as np
from numpy import linalg
import matplotlib.pyplot as plt
from time import time
import DataReader as DR

def find_anchors_from_class_0(X,y):
    y_u=np.unique(y)
    points_of_class1=X[y==y_u[1]]
    distances= []
    #find all distances
    for a in points_of_class1:
        distances.append(linalg.norm(X[y== y_u[0]] -a, axis=1, ord=1))
    #get minimum distances
    distances = np.array(distances)
    min_distance= np.min(distances,axis=0)
    min_distance = min_distance / np.std(min_distance)

    inds=[]
    # get the indicies of the points whose minimum distance is greater than the lower bound
    inds.append(min_distance >= np.min(min_distance)+np.std(min_distance)*0.43)
    inds.append(min_distance <= np.max(min_distance))
    all_inds= np.all(inds,axis=0)
    tempX=X[y==y_u[0]]
    return  tempX[all_inds]

def find_anchors_from_class_1(X, y):
    y_u = np.unique(y)
    points_of_class0 = X[y == y_u[0]]
    distances = []
    # find all distances
    for a in points_of_class0:
        distances.append(linalg.norm(X[y == y_u[1]] - a, axis=1, ord=1))
    # get minimum distances
    distances = np.array(distances)
    min_distance = np.min(distances, axis=0)
    min_distance = min_distance / np.std(min_distance)
    inds = []
    # get the indicies of the points whose minimum distance is greater than the lower bound
    inds.append(min_distance >= np.min(min_distance) + np.std(min_distance) * 0.7)
    inds.append(min_distance <= np.max(min_distance))
    all_inds = np.all(inds, axis=0)
    tempX = X[y == y_u[1]]
    return tempX[all_inds]

def map_min_1_2(X, anchors):
    temp = []
    # begin: argmin
    for a in anchors[0]:
        temp.append(linalg.norm(X - a, axis=1, ord=1))
    temp = np.array(temp)
    mins1 = np.min(temp, axis=0)
    mins1 = mins1 / np.std(mins1)
    temp = []
    for a in anchors[1]:
        temp.append(linalg.norm(X - a, axis=1, ord=1))
    temp = np.array(temp)
    mins2 = np.min(temp, axis=0)
    mins2 = mins2 / np.std(mins2)
    nz = [mins1 != 0, mins2 != 0]
    d = np.all(nz, axis=0)
    mins1[mins1 == 0] = np.mean(mins1[d])
    mins2[mins2 == 0] = np.mean(mins2[d])
    # end: argmin
    X1 = np.hstack((X, mins1.reshape((len(X), 1))))
    return np.hstack((X1, mins2.reshape((len(X), 1))))

X_train, X_test, y_train, y_test = DR.Phoneme()
y_unique = np.unique(y_train)
sets_of_anchors = []
sets_of_anchors.append(find_anchors_from_class_0(X_train,y_train))
sets_of_anchors.append(find_anchors_from_class_1(X_train, y_train))

XD= map_min_1_2(X_train,sets_of_anchors)
XT= map_min_1_2(X_test,sets_of_anchors)
s=time()
clf = LinearSVC(C=1, dual=False).fit(XD,y_train)
train_time=round(time()-s,4)
score_training= accuracy_score(y_train, clf.predict(XD))
score_test=round(100*accuracy_score(y_test, clf.predict(XT)),2)
score_training= round(100*score_training,2)

plt.scatter(XD[y_train==y_unique[0]][:, -2], XD[y_train==y_unique[0]][:, -1], c='red', edgecolors='k')
plt.scatter(XD[y_train == y_unique[1]][:, -2], XD[y_train == y_unique[1]][:, -1], c='blue', edgecolors='k')
plt.xlabel(r'$||x_i-a_i^{(1)}||$')
plt.ylabel(r'$||x_i-a_i^{(2)}||$')
plt.text(XD[:,-2].max() - 0.5,XD[:,-1].max() - 2, "Training Acc.   :  %.2f\nTest Acc.         :  %.2f\nTraining T.       :  %.3f"% (score_training,score_test,train_time),  size=10,
                     va="baseline", ha="right", multialignment="left",
                     bbox=dict(fc="none"))

plt.title('Phonome')

plt.show()

ext: py

max_stars_repo_name: sibirbil/SimpleKernels

max_issues_repo_name: sibirbil/SimpleKernels

max_forks_repo_name: sibirbil/SimpleKernels

reason: import numpy,from numpy


text:
#
# Stephen Vondenstein
# 10/03/2018
#
import pandas as pd
import numpy as np
import os
from keras.preprocessing import image
from tqdm import tqdm
from src.utils.utility import get_max_unused_filename

def rle(img, order='F', format=True):
    bytes = img.reshape(img.shape[0] * img.shape[1], order=order)
    runs = []  # list of run lengths
    r = 0  # the current run length
    pos = 1  # count starts from 1 per WK
    for c in bytes:
        if (c == 0):
            if r != 0:
                runs.append((pos, r))
                pos += r
                r = 0
            pos += 1
        else:
            r += 1

    # if last run is unsaved (i.e. data ends with 1)
    if r != 0:
        runs.append((pos, r))
        pos += r
        r = 0

    if format:
        z = ''

        for rr in runs:
            z += '{} {} '.format(rr[0], rr[1])
        return z[:-1]
    else:
        return runs

def prepare_submission(source_dir, output_path, sub_prefix):
    pred_ids = next(os.walk(source_dir))[2]
    preds = []

    tt = tqdm(range(len(pred_ids)), total=len(pred_ids),
              desc="Loading ")

    for t in tt:
        img = image.load_img(source_dir + '/' + pred_ids[t])
        x = image.img_to_array(img)[:, :, 1]
        preds.append(x)

    tt.close()

    tt = tqdm(range(len(pred_ids)), total=len(pred_ids),
              desc="Computing RLE ")

    pred_dict = {}
    for i, fn in enumerate(pred_ids):
        pred_dict_val = {fn[:-4]: rle(np.round(preds[i]))}
        pred_dict.update(pred_dict_val)
        tt.update()

    tt.close()

    print('Preparing submission...')
    sub = pd.DataFrame.from_dict(pred_dict, orient='index')
    sub.index.names = ['id']
    sub.columns = ['rle_mask']

    i = 1
    while os.path.exists(os.path.join(output_path, '%s-%s.csv' % (sub_prefix, i))):
        i += 1
    submission_path = get_max_unused_filename(output_path, sub_prefix, '.csv')
    sub.to_csv(submission_path)

    print('Submission saved to ' + submission_path)

    return submission_path

max_stars_repo_name: neurologists/ionic

max_issues_repo_name: neurologists/ionic

max_forks_repo_name: KaggleBuoys/ionic

reason: import numpy


text:
import numpy as np
import numpy.testing as npt
import nitime
import nitime.algorithms as tsa

def test_xcorr_zscored():
    """

    Test this function, which is not otherwise tested in the testing of the
    EventRelatedAnalyzer

    """

    cycles = 10
    l = 1024
    unit = 2 * np.pi / l
    t = np.arange(0, 2 * np.pi + unit, unit)
    signal = np.sin(cycles * t)
    events = np.zeros(t.shape)
    #Zero crossings:
    idx = np.where(np.abs(signal) < 0.03)[0]
    #An event occurs at the beginning of every cycle:
    events[idx[:-2:2]] = 1

    a = tsa.freq_domain_xcorr_zscored(signal, events, 1000, 1000)
    npt.assert_almost_equal(np.mean(a), 0, 1)
    npt.assert_almost_equal(np.std(a), 1, 1)

max_stars_repo_name: Xunius/nitime

max_issues_repo_name: Xunius/nitime

max_forks_repo_name: Xunius/nitime


text:
"""
===============================
Triangle plot with one data set
===============================

Plot triangle plot with one data set using some default settings.
"""
import coplot.plot_contours as plc
import matplotlib.pyplot as plt
import numpy as np

test1 = np.random.randn(50000, 3)

# If you want to use some default settings, you can use it like this:
fig = plc.Contours(test1).plot(labels=[r'$a$', r'$b$', r'$c$'], legend=True, subplots_adjust=False)
# and you can save the figure like this:
# pl.savefig('test_triangle_1.pdf', fig)
plt.show()

max_stars_repo_name: Guo-Jian-Wang/coplot

max_issues_repo_name: Guo-Jian-Wang/coplot

max_forks_repo_name: Guo-Jian-Wang/coplot


text:
from __future__ import division, print_function, absolute_import

from numpy import sqrt, inner, zeros, inf, finfo
from numpy.linalg import norm

from .utils import make_system

__all__ = ['minres']

def minres(A, b, x0=None, shift=0.0, tol=1e-5, maxiter=None,
           M=None, callback=None, show=False, check=False):
    """
    Use MINimum RESidual iteration to solve Ax=b

    MINRES minimizes norm(A*x - b) for a real symmetric matrix A.  Unlike
    the Conjugate Gradient method, A can be indefinite or singular.

    If shift != 0 then the method solves (A - shift*I)x = b

    Parameters
    ----------
    A : {sparse matrix, dense matrix, LinearOperator}
        The real symmetric N-by-N matrix of the linear system
        Alternatively, ``A`` can be a linear operator which can
        produce ``Ax`` using, e.g.,
        ``scipy.sparse.linalg.LinearOperator``.
    b : {array, matrix}
        Right hand side of the linear system. Has shape (N,) or (N,1).

    Returns
    -------
    x : {array, matrix}
        The converged solution.
    info : integer
        Provides convergence information:
            0  : successful exit
            >0 : convergence to tolerance not achieved, number of iterations
            <0 : illegal input or breakdown

    Other Parameters
    ----------------
    x0  : {array, matrix}
        Starting guess for the solution.
    tol : float
        Tolerance to achieve. The algorithm terminates when the relative
        residual is below `tol`.
    maxiter : integer
        Maximum number of iterations.  Iteration will stop after maxiter
        steps even if the specified tolerance has not been achieved.
    M : {sparse matrix, dense matrix, LinearOperator}
        Preconditioner for A.  The preconditioner should approximate the
        inverse of A.  Effective preconditioning dramatically improves the
        rate of convergence, which implies that fewer iterations are needed
        to reach a given error tolerance.
    callback : function
        User-supplied function to call after each iteration.  It is called
        as callback(xk), where xk is the current solution vector.

    References
    ----------
    Solution of sparse indefinite systems of linear equations,
        C. C. Paige and M. A. Saunders (1975),
        SIAM J. Numer. Anal. 12(4), pp. 617-629.
        https://web.stanford.edu/group/SOL/software/minres/

    This file is a translation of the following MATLAB implementation:
        https://web.stanford.edu/group/SOL/software/minres/minres-matlab.zip

    """
    A, M, x, b, postprocess = make_system(A, M, x0, b)

    matvec = A.matvec
    psolve = M.matvec

    first = 'Enter minres.   '
    last = 'Exit  minres.   '

    n = A.shape[0]

    if maxiter is None:
        maxiter = 5 * n

    msg = [' beta2 = 0.  If M = I, b and x are eigenvectors    ',   # -1
            ' beta1 = 0.  The exact solution is  x = 0          ',   # 0
            ' A solution to Ax = b was found, given rtol        ',   # 1
            ' A least-squares solution was found, given rtol    ',   # 2
            ' Reasonable accuracy achieved, given eps           ',   # 3
            ' x has converged to an eigenvector                 ',   # 4
            ' acond has exceeded 0.1/eps                        ',   # 5
            ' The iteration limit was reached                   ',   # 6
            ' A  does not define a symmetric matrix             ',   # 7
            ' M  does not define a symmetric matrix             ',   # 8
            ' M  does not define a pos-def preconditioner       ']   # 9

    if show:
        print(first + 'Solution of symmetric Ax = b')
        print(first + 'n      =  %3g     shift  =  %23.14e' % (n,shift))
        print(first + 'itnlim =  %3g     rtol   =  %11.2e' % (maxiter,tol))
        print()

    istop = 0
    itn = 0
    Anorm = 0
    Acond = 0
    rnorm = 0
    ynorm = 0

    xtype = x.dtype

    eps = finfo(xtype).eps

    x = zeros(n, dtype=xtype)

    # Set up y and v for the first Lanczos vector v1.
    # y  =  beta1 P' v1,  where  P = C**(-1).
    # v is really P' v1.

    y = b
    r1 = b

    y = psolve(b)

    beta1 = inner(b,y)

    if beta1 < 0:
        raise ValueError('indefinite preconditioner')
    elif beta1 == 0:
        return (postprocess(x), 0)

    beta1 = sqrt(beta1)

    if check:
        # are these too strict?

        # see if A is symmetric
        w = matvec(y)
        r2 = matvec(w)
        s = inner(w,w)
        t = inner(y,r2)
        z = abs(s - t)
        epsa = (s + eps) * eps**(1.0/3.0)
        if z > epsa:
            raise ValueError('non-symmetric matrix')

        # see if M is symmetric
        r2 = psolve(y)
        s = inner(y,y)
        t = inner(r1,r2)
        z = abs(s - t)
        epsa = (s + eps) * eps**(1.0/3.0)
        if z > epsa:
            raise ValueError('non-symmetric preconditioner')

    # Initialize other quantities
    oldb = 0
    beta = beta1
    dbar = 0
    epsln = 0
    qrnorm = beta1
    phibar = beta1
    rhs1 = beta1
    rhs2 = 0
    tnorm2 = 0
    gmax = 0
    gmin = finfo(xtype).max
    cs = -1
    sn = 0
    w = zeros(n, dtype=xtype)
    w2 = zeros(n, dtype=xtype)
    r2 = r1

    if show:
        print()
        print()
        print('   Itn     x(1)     Compatible    LS       norm(A)  cond(A) gbar/|A|')

    while itn < maxiter:
        itn += 1

        s = 1.0/beta
        v = s*y

        y = matvec(v)
        y = y - shift * v

        if itn >= 2:
            y = y - (beta/oldb)*r1

        alfa = inner(v,y)
        y = y - (alfa/beta)*r2
        r1 = r2
        r2 = y
        y = psolve(r2)
        oldb = beta
        beta = inner(r2,y)
        if beta < 0:
            raise ValueError('non-symmetric matrix')
        beta = sqrt(beta)
        tnorm2 += alfa**2 + oldb**2 + beta**2

        if itn == 1:
            if beta/beta1 <= 10*eps:
                istop = -1  # Terminate later

        # Apply previous rotation Qk-1 to get
        #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
        #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].

        oldeps = epsln
        delta = cs * dbar + sn * alfa   # delta1 = 0         deltak
        gbar = sn * dbar - cs * alfa   # gbar 1 = alfa1     gbar k
        epsln = sn * beta     # epsln2 = 0         epslnk+1
        dbar = - cs * beta   # dbar 2 = beta2     dbar k+1
        root = norm([gbar, dbar])
        Arnorm = phibar * root

        # Compute the next plane rotation Qk

        gamma = norm([gbar, beta])       # gammak
        gamma = max(gamma, eps)
        cs = gbar / gamma             # ck
        sn = beta / gamma             # sk
        phi = cs * phibar              # phik
        phibar = sn * phibar              # phibark+1

        # Update  x.

        denom = 1.0/gamma
        w1 = w2
        w2 = w
        w = (v - oldeps*w1 - delta*w2) * denom
        x = x + phi*w

        # Go round again.

        gmax = max(gmax, gamma)
        gmin = min(gmin, gamma)
        z = rhs1 / gamma
        rhs1 = rhs2 - delta*z
        rhs2 = - epsln*z

        # Estimate various norms and test for convergence.

        Anorm = sqrt(tnorm2)
        ynorm = norm(x)
        epsa = Anorm * eps
        epsx = Anorm * ynorm * eps
        epsr = Anorm * ynorm * tol
        diag = gbar

        if diag == 0:
            diag = epsa

        qrnorm = phibar
        rnorm = qrnorm
        if ynorm == 0 or Anorm == 0:
            test1 = inf
        else:
            test1 = rnorm / (Anorm*ynorm)    # ||r||  / (||A|| ||x||)
        if Anorm == 0:
            test2 = inf
        else:
            test2 = root / Anorm            # ||Ar|| / (||A|| ||r||)

        # Estimate  cond(A).
        # In this version we look at the diagonals of  R  in the
        # factorization of the lower Hessenberg matrix,  Q * H = R,
        # where H is the tridiagonal matrix from Lanczos with one
        # extra row, beta(k+1) e_k^T.

        Acond = gmax/gmin

        # See if any of the stopping criteria are satisfied.
        # In rare cases, istop is already -1 from above (Abar = const*I).

        if istop == 0:
            t1 = 1 + test1      # These tests work if tol < eps
            t2 = 1 + test2
            if t2 <= 1:
                istop = 2
            if t1 <= 1:
                istop = 1

            if itn >= maxiter:
                istop = 6
            if Acond >= 0.1/eps:
                istop = 4
            if epsx >= beta1:
                istop = 3
            # if rnorm <= epsx   : istop = 2
            # if rnorm <= epsr   : istop = 1
            if test2 <= tol:
                istop = 2
            if test1 <= tol:
                istop = 1

        # See if it is time to print something.

        prnt = False
        if n <= 40:
            prnt = True
        if itn <= 10:
            prnt = True
        if itn >= maxiter-10:
            prnt = True
        if itn % 10 == 0:
            prnt = True
        if qrnorm <= 10*epsx:
            prnt = True
        if qrnorm <= 10*epsr:
            prnt = True
        if Acond <= 1e-2/eps:
            prnt = True
        if istop != 0:
            prnt = True

        if show and prnt:
            str1 = '%6g %12.5e %10.3e' % (itn, x[0], test1)
            str2 = ' %10.3e' % (test2,)
            str3 = ' %8.1e %8.1e %8.1e' % (Anorm, Acond, gbar/Anorm)

            print(str1 + str2 + str3)

            if itn % 10 == 0:
                print()

        if callback is not None:
            callback(x)

        if istop != 0:
            break  # TODO check this

    if show:
        print()
        print(last + ' istop   =  %3g               itn   =%5g' % (istop,itn))
        print(last + ' Anorm   =  %12.4e      Acond =  %12.4e' % (Anorm,Acond))
        print(last + ' rnorm   =  %12.4e      ynorm =  %12.4e' % (rnorm,ynorm))
        print(last + ' Arnorm  =  %12.4e' % (Arnorm,))
        print(last + msg[istop+1])

    if istop == 6:
        info = maxiter
    else:
        info = 0

    return (postprocess(x),info)

if __name__ == '__main__':
    from numpy import zeros, arange
    from scipy.linalg import norm
    from scipy.sparse import spdiags

    n = 10

    residuals = []

    def cb(x):
        residuals.append(norm(b - A*x))

    # A = poisson((10,),format='csr')
    A = spdiags([arange(1,n+1,dtype=float)], [0], n, n, format='csr')
    M = spdiags([1.0/arange(1,n+1,dtype=float)], [0], n, n, format='csr')
    A.psolve = M.matvec
    b = zeros(A.shape[0])
    x = minres(A,b,tol=1e-12,maxiter=None,callback=cb)
    # x = cg(A,b,x0=b,tol=1e-12,maxiter=None,callback=cb)[0]

max_stars_repo_name: tuxcell/scipy

max_issues_repo_name: tuxcell/scipy

max_forks_repo_name: tuxcell/scipy

reason: from numpy,from scipy


text:
""" Python module containing a class for basic interface of an ML model, as well as functions
shared by multiple model types.
"""

import os
import logging
import numpy as np
from tqdm import tqdm

import tensorflow as tf

from .layer_functions import add_linear_layer

# -------------------------------------------------------------------------------------------------
# Workaround TF logger problem in TF 1.14
try:
    # Capirca uses Google's abseil-py library, which uses a Google-specific
    # wrapper for logging. That wrapper will write a warning to sys.stderr if
    # the Google command-line flags library has not been initialized.
    #
    # https://github.com/abseil/abseil-py/blob/pypi-v0.7.1/absl/logging/__init__.py#L819-L825
    #
    # This is not right behavior for Python code that is invoked outside of a
    # Google-authored main program. Use knowledge of abseil-py to disable that
    # warning; ignore and continue if something goes wrong.
    import absl.logging

    # https://github.com/abseil/abseil-py/issues/99
    logging.root.removeHandler(absl.logging._absl_handler)  #pylint: disable=protected-access
    # https://github.com/abseil/abseil-py/issues/102
    absl.logging._warn_preinit_stderr = False               #pylint: disable=protected-access
except Exception:                                           #pylint: disable=broad-except
    pass
# Workaround TF logger problem in TF 1.14
# -------------------------------------------------------------------------------------------------

class BaseModel():
    """ This class defines a common model interface. The models defined within this project
    should explicitly implement all depicted interface functions, to abstract the details
    away from the main scripts.

    :param model_settings: a dictionary containing the model settings. All dictionary key/value
        pairs will be set as class members, for further readibility (and IDE auto-complete :D)
    """

    # To add in the future:
    # - Add TensorBoard functionality
    # - Instead of "train_batch()" with an outer control loop, create a "train_generator()" that
    #   uses generators
    # - Similarly, crate a "predict_generator()"
    # - Add an MC dropout flag to the predict functions, to enable uncertainty prediction

    def __init__(self, model_settings):

        # Instatiates basic model settings
        self.model_folder = model_settings["model_folder"]
        self.input_name = model_settings["input_name"]
        self.input_shape = model_settings["input_shape"]
        self.input_type = model_settings["input_type"]
        self.output_name = model_settings["output_name"]
        self.output_type = model_settings["output_type"]
        self.output_shape = model_settings["output_shape"]
        self.validation_metric = model_settings["validation_metric"]
        self.batch_size = model_settings["batch_size"]
        self.early_stopping = model_settings["early_stopping"]
        self.max_epochs = model_settings["max_epochs"]
        self.learning_rate = model_settings["learning_rate"]
        self.learning_rate_decay = model_settings["learning_rate_decay"]
        self.optimizer_type = model_settings["optimizer_type"]

        # Instatiates optional model settings
        self.batch_size_inference = model_settings["batch_size_inference"] \
            if "batch_size_inference" in model_settings else self.batch_size
        self.val_eval_period = model_settings["val_eval_period"] \
            if "val_eval_period" in model_settings else 1
        self.dropout = model_settings["dropout"] \
            if "dropout" in model_settings else 0.
        self.fc_layers = model_settings["fc_layers"] \
            if "fc_layers" in model_settings else 0
        self.fc_neurons = model_settings["fc_neurons"] \
            if "fc_neurons" in model_settings else 0

        # Instantiates other common variables that might be used later
        self.current_validation_score = None
        self.current_learning_rate = None
        self.current_epoch = None
        self.epochs_not_improving = None
        self.best_validation_score = None
        self.saver = None
        self.session = None
        self.train_step = None
        self.learning_rate_var = None
        self.dropout_var = None
        self.is_training = None
        self.model_input = None
        self.model_output = None
        self.model_target = None

        # Runs basic initializations
        self._set_gpu(model_settings)

    # ---------------------------------------------------------------------------------------------
    # Model interface functions
    def set_graph(self):
        """Prototype: setup(self)

        Given the settings, sets up the model graph for training
        """
        raise NotImplementedError(
            "The model sub-class did not implement a 'set_graph()' function."
        )

    def train_epoch(self, X, Y):
        """Prototype: train_batch(self, X, Y)

        Trains on a single epoch, given the data (as numpy variables).
        """
        raise NotImplementedError(
            "The model sub-class did not implement a 'train_epoch()' function."
        )

    def epoch_end(self, X=None, Y=None):
        """Prototype: epoch_end(self, X, Y)

        Executes end of epoch logic (e.g. decay learning rate, try early stopping, ...) and
        returns whether the model should keep training or not
        """
        raise NotImplementedError(
            "The model sub-class did not implement a 'epoch_end()' function."
        )

    def predict(self, X, validation):
        """Prototype: predict(self, X)

        Returns the preditions on the given data
        """
        raise NotImplementedError(
            "The model sub-class did not implement a 'predict()' function."
        )

    def save(self, model_name):
        """Prototype: save(self, model_name)

        Stores all model data inside the specified folder, given the model name
        """
        raise NotImplementedError(
            "The model sub-class did not implement a 'save()' function."
        )

    def load(self, model_name):
        """Prototype: load(self, model_name)

        Loads all model data from the specified folder, given the model name
        """
        raise NotImplementedError(
            "The model sub-class did not implement a 'load()' function."
        )

    def close(self):
        """Prototype: close(self)

        Cleans up the session and any left over data
        """
        raise NotImplementedError(
            "The model sub-class did not implement a 'close()' function."
        )

    # ---------------------------------------------------------------------------------------------
    # Non-interface functions: model training
    def _train_epoch(self, X, Y):
        """Default training routine - trains the model for an epoch

        :param X: numpy array with the features
        :param Y: numpy array with the labels
        """
        #local import to avoid issues in nvidia jetson <-> sklearn
        from sklearn.utils import shuffle

        assert X.shape[0] == Y.shape[0], "X and Y have a different number of samples!"
        self.current_epoch += 1
        max_batches = int(np.ceil(X.shape[0] / self.batch_size))
        X, Y = shuffle(X, Y)
        train_string = "Training on epoch: {:4} || LR: {:3.2E} ||"

        for batch_idx in tqdm(
            range(max_batches),
            desc=train_string.format(self.current_epoch, self.current_learning_rate)
        ):
            start_batch = batch_idx * self.batch_size
            end_batch = min((batch_idx + 1) * self.batch_size, X.shape[0])
            self.train_step.run(
                feed_dict={
                    self.model_input: X[start_batch:end_batch, ...],
                    self.model_target: Y[start_batch:end_batch, ...],
                    self.dropout_var: self.dropout,
                    self.learning_rate_var: self.current_learning_rate,
                    self.is_training: True
                },
                session=self.session)

    def _epoch_end(self, y_true=None, y_pred=None):
        """ Default end of epoch routine - Performs end of epoch operations, such as decaying the
        learning rate. Some operations, such as the early stopping, require a validation set.

        :param y_true: ground truth, defaults to None
        :param y_pred: model predictions, defaults to None
        :returns: boolean indicating whether the model should keep training, validation score
        """
        #local import to avoid issues in nvidia jetson <-> sklearn
        from .metrics import score_predictions

        keep_training = True
        val_score = None

        # Decays LR
        self.current_learning_rate *= self.learning_rate_decay
        if y_true is not None and y_pred is not None:
            len_pred = y_pred.shape[0]
            val_score = score_predictions(y_true[:len_pred, :], y_pred, self.validation_metric)
            if self.early_stopping:
                # Evaluates early stopping, if requested
                keep_training = self._eval_early_stopping(val_score)

        if self.current_epoch >= self.max_epochs:
            keep_training = False
        return keep_training, val_score

    def _predict(self, X, validation):
        """ Returns the predictions on the given data

        :param X: numpy array with the features
        :param validation: boolean indicating whether these predictions have validation purposes
        :return: an numpy array with the predictions
        """
        if validation and self.current_epoch % self.val_eval_period != 0:
            return None
        max_batches = int(np.ceil(X.shape[0] / self.batch_size))
        predictions = []
        for batch_idx in range(max_batches):
            start_batch = batch_idx * self.batch_size_inference
            end_batch = min((batch_idx + 1) * self.batch_size_inference, X.shape[0])
            batch_predictions = self.model_output.eval(
                feed_dict={
                    self.model_input: X[start_batch:end_batch, ...],
                    self.dropout_var: 0.0,
                    self.is_training: False
                },
                session=self.session
            )
            batch_size = end_batch - start_batch
            predictions.extend([batch_predictions[idx, ...] for idx in range(batch_size)])
        return np.asarray(predictions)

    def _eval_early_stopping(self, validation_score):
        """ Evaluates the early stopping mechanism

        :param validation_score: the validation score
        :return: boolean depicting whether the model should keep training
        """
        keep_training = True
        if self.current_validation_score is None:
            self.current_validation_score = validation_score
            self.epochs_not_improving = 0
        elif validation_score < self.current_validation_score:
            self.current_validation_score = validation_score
            self.epochs_not_improving = 0
        else:
            self.epochs_not_improving += self.val_eval_period
        if self.epochs_not_improving >= self.early_stopping:
            logging.info("Early stopping - the model does not improve for %s epochs",
                self.epochs_not_improving)
            keep_training = False
        return keep_training

    # ---------------------------------------------------------------------------------------------
    # Non-interface functions: model input/output
    def _set_graph_io(self):
        """ Auxiliary function to "set_graph()". Sets the graph input and trainable target,
        as well as some other basic variables common to all model types
        """
        # The current learning rate
        self.learning_rate_var = tf.compat.v1.placeholder(tf.float32, shape=[])
        # Dropout probability
        self.dropout_var = tf.compat.v1.placeholder(tf.float32, name='dropout')
        # Boolean that is true when the model is training
        self.is_training = tf.compat.v1.placeholder(tf.bool, name='is_training')

        if self.input_type == "float":
            self.model_input = tf.compat.v1.placeholder(
                tf.float32,
                shape=[None] + self.input_shape,
                name=self.input_name,
            )
        elif self.input_type == "bool":
            original_input = tf.compat.v1.placeholder(
                tf.bool,
                shape=[None] + self.input_shape,
                name=self.input_name,
            )
            self.model_input = tf.cast(original_input, tf.float32)

        if self.output_type == "regression":
            self.model_target = tf.compat.v1.placeholder(
                tf.float32,
                shape=[None] + self.output_shape
            )
        elif self.output_type == "classification":
            self.model_target = tf.compat.v1.placeholder(
                tf.int64,
                shape=[None]
            )
        else:
            raise ValueError("Unknown 'output_type' ({}). Only 'classification' and 'regression' "
                "are accepted.".format(self.output_type))

    def _add_classification_output(self, input_data):
        """ Adds the last layer for classification models. Returns the trainable step.

        :param input_data: TF tensor with this layer's input
        :returns: the trainable step
        """
        # Defines the logits and the softmax output
        assert len(self.output_shape) == 1, "This function is not ready for outputs with higher "\
            "dimensionality"
        logits = add_linear_layer(input_data, self.output_shape[0])
        self.model_output = tf.nn.softmax(logits, name=self.output_name)

        # Defines the loss function [mean(cross_entropy(target_value - softmax))]
        cross_entropy = tf.reduce_mean(tf.nn.sparse_softmax_cross_entropy_with_logits(
            labels=self.model_target,
            logits=logits
        ))

        # Defines the optimizer and the train step
        optimizer = self._set_optimizer()
        train_step = optimizer.minimize(cross_entropy)
        return train_step

    def _add_regression_output(self, input_data):
        """ Adds the last layer for regression models. Returns the trainable step.

        :param input_data: TF tensor with this layer's input
        :returns: the trainable step
        """
        # Defines the regression output (used to learn), and its clipped version (used to predict)
        assert len(self.output_shape) == 1, "This function is not ready for outputs with higher "\
            "dimensionality"
        regression = add_linear_layer(input_data, self.output_shape[0], bias=0.5)
        self.model_output = tf.clip_by_value(regression, 0.0, 1.0, name=self.output_name)

        # Defines the loss function [MSE = mean(square(target_value - regression))]
        mse = tf.reduce_mean(tf.square(self.model_target - regression))

        # Defines the optimizer and the train step
        optimizer = self._set_optimizer()
        train_step = optimizer.minimize(mse)
        return train_step

    def _set_optimizer(self):
        """ Returns a TF optimizer, given the model settings
        """
        if self.optimizer_type == "ADAM":
            return tf.compat.v1.train.AdamOptimizer(learning_rate=self.learning_rate_var)
        elif self.optimizer_type == "NADAM":
            return tf.contrib.opt.NadamOptimizer(learning_rate=self.learning_rate_var)
        else:
            raise ValueError("{} is not a supported optimizer type. Supported optimizer types: "
                "ADAM, NADAM.".format(self.optimizer_type))

    # ---------------------------------------------------------------------------------------------
    # Non-interface functions: misc
    @staticmethod
    def _set_gpu(model_settings):
        """ If the option 'target_gpu' is defined in `model_settings`, sets that GPU
        """
        if "target_gpu" in model_settings:
            target_gpu = model_settings["target_gpu"]
            logging.info("[Using GPU #%s]", target_gpu)
            os.environ["CUDA_VISIBLE_DEVICES"] = str(target_gpu)

    def _prepare_model_for_training(self):
        """ Self-documenting :D
        """
        # Sets the saver and the session
        self.saver = tf.compat.v1.train.Saver()
        self._set_session()

        # Initializes TF variables
        self.session.run(tf.compat.v1.global_variables_initializer())
        trainable_parameters = int(np.sum(
            [np.product([var_dim.value for var_dim in var.get_shape()])
            for var in tf.compat.v1.trainable_variables()]
        ))
        logging.info("Model initialized with %s trainable parameters!", trainable_parameters)

        # Initializes other train-related variables
        self.current_learning_rate = self.learning_rate
        self.current_epoch = 0

    def _save(self, model_name):
        """ Default function to save a model

        :param model_name: the name of the model
        """
        if not os.path.exists(self.model_folder):
            os.makedirs(self.model_folder)
        save_path = os.path.join(self.model_folder, model_name)
        self.saver.save(self.session, save_path)

    def _load(self, model_name):
        """ Default function to load a model

        :param model_name: the name of the model
        """
        self._set_session()
        loader_path = os.path.join(self.model_folder, model_name)
        loader = tf.compat.v1.train.import_meta_graph(loader_path + '.meta')
        loader.restore(self.session, loader_path)

        # Redefines key model graph-related variables
        # (e.g. placeholders and operations that were saved by name)
        graph = tf.compat.v1.get_default_graph()
        self.model_input = graph.get_tensor_by_name(self.input_name + ":0")
        self.model_output = graph.get_tensor_by_name(self.output_name + ":0")
        self.dropout_var = graph.get_tensor_by_name("dropout:0")
        self.is_training = graph.get_tensor_by_name("is_training:0")

    def _set_session(self):
        """ Default function to set the session
        """
        self.session = tf.compat.v1.Session(
            config=tf.compat.v1.ConfigProto(gpu_options={"allow_growth": True})
        )

    def _close_session(self):
        """ Default function to close the session
        """
        self.session.close()


text:
import unittest

import numpy
import six
import warnings

import chainer
from chainer.backends import cuda
from chainer import functions
from chainer import gradient_check
from chainer import testing
from chainer.testing import attr
from chainer.testing import backend

def _as_noncontiguous_array(array):
    # TODO(niboshi): cupy + cudnn test fails in F.fixed_batch_normalization.
    # Fix it and use testing.array._as_noncontiguous_array.
    def as_noncontiguous_array(arr):
        if arr is None:
            return None
        if isinstance(arr, (numpy.ndarray, cuda.ndarray)):
            xp = chainer.backend.get_array_module(arr)
            return xp.asfortranarray(arr)
        return testing.array._as_noncontiguous_array(arr)

    if isinstance(array, (list, tuple)):
        return type(array)([as_noncontiguous_array(arr) for arr in array])
    return as_noncontiguous_array(array)

def _batch_normalization(
        inputs, running_mean=None, running_var=None, decay=None):
    x, gamma, beta, mean, var, eps, expander = inputs
    mean_expanded = mean[expander]
    std = numpy.sqrt(var + eps)[expander]
    y_expect = (gamma[expander] * (x - mean_expanded) / std + beta[expander])

    if running_mean is not None or running_var is not None:
        m = x.size // gamma.size
        adjust = m / max(m - 1., 1.)  # unbiased estimation
        if running_mean is not None:
            running_mean *= decay
            running_mean += (1 - decay) * mean
        if running_var is not None:
            running_var *= decay
            running_var += (1 - decay) * adjust * var

    return y_expect

@testing.parameterize(*(testing.product_dict(
    testing.product({
        'param_shape': [(3,), (3, 4), (3, 2, 3)],
        'ndim': [0, 1, 2],
    }) + [
        {'input_shape': (5, 4, 3, 2), 'axis': (0, 2, 3)},
        {'input_shape': (5, 4), 'axis': 0},
        {'input_shape': (5, 4, 3), 'axis': (0, 1)},
    ],
    testing.product({
        'dtype': [numpy.float32],
        'eps': [2e-5, 5e-1],
        'c_contiguous': [True, False],
        'running_statistics': [True, False],
    }),
) + testing.product({
    'param_shape': [(3,)],
    'ndim': [1],
    'eps': [2e-5, 5e-1],
    'dtype': [numpy.float16, numpy.float32, numpy.float64],
    'c_contiguous': [True, False],
    'running_statistics': [True, False],
})))
@backend.inject_backend_tests(
    ['test_forward', 'test_backward', 'test_double_backward'],
    # CPU tests
    testing.product({
        'use_cuda': [False],
        'use_ideep': ['never', 'always'],
    })
    # GPU tests
    + testing.product({
        'use_cuda': [True],
        'use_cudnn': ['never', 'always'],
        'cudnn_fast_batch_normalization': [True, False],
    })
    # ChainerX tests
    + [
        {'use_chainerx': True, 'chainerx_device': 'native:0'},
        {'use_chainerx': True, 'chainerx_device': 'cuda:0'},
    ])
class TestBatchNormalization(unittest.TestCase):

    def setUp(self):
        dtype = self.dtype

        if not hasattr(self, 'axis'):
            param_shape = self.param_shape
            ndim = self.ndim
            shape = (5,) + param_shape + (2,) * ndim
        else:
            aggr_axes = self.axis
            if isinstance(self.axis, int):
                aggr_axes = self.axis,
            param_shape = tuple(
                s
                for i, s in enumerate(self.input_shape)
                if i not in aggr_axes
            )
            shape = self.input_shape

        gamma = numpy.random.uniform(.5, 1, param_shape).astype(dtype)
        beta = numpy.random.uniform(-1, 1, param_shape).astype(dtype)
        x = numpy.random.uniform(-1, 1, shape).astype(dtype)
        gy = numpy.random.uniform(-1, 1, shape).astype(dtype)
        ggx = numpy.random.uniform(-1, 1, shape).astype(dtype)
        gggamma = numpy.random.uniform(-1, 1, param_shape).astype(dtype)
        ggbeta = numpy.random.uniform(-1, 1, param_shape).astype(dtype)

        if self.running_statistics:
            self.running_mean = numpy.random.uniform(
                -1, 1, param_shape).astype(dtype)
            self.running_var = numpy.random.uniform(
                -1, 1, param_shape).astype(dtype)
        else:
            self.running_mean = None
            self.running_var = None

        if not hasattr(self, 'axis'):
            head_ndim = gamma.ndim + 1
            aggr_axes = (0,) + tuple(six.moves.range(head_ndim, x.ndim))

            self.expander = (None, Ellipsis) + (None,) * ndim
        else:
            self.expander = tuple(
                None if i in aggr_axes else slice(None)
                for i in range(x.ndim)
            )

        mean = x.mean(axis=aggr_axes)
        var = x.var(axis=aggr_axes)

        self.decay = 0.9
        self.mean = mean
        self.var = var

        self.inputs = [x, gamma, beta]
        self.grad_outputs = [gy]
        self.grad_grad_inputs = [ggx, gggamma, ggbeta]

        self.bn_options = {
            'decay': self.decay,
            'eps': self.eps,
        }
        if hasattr(self, 'axis'):
            self.bn_options['axis'] = self.axis
        self.check_forward_options = {'atol': 1e-4, 'rtol': 1e-3}
        self.check_backward_options = {'dtype': numpy.float64}
        self.check_double_backward_options = {
            'dtype': numpy.float64, 'atol': 1e-3, 'rtol': 1e-2}
        if self.dtype == numpy.float16:
            self.check_forward_options = {'atol': 1e-2, 'rtol': 1e-2}
            self.check_backward_options = {
                'dtype': numpy.float64, 'atol': 1e-2, 'rtol': 1e-2}
            self.check_double_backward_options = {
                'dtype': numpy.float64, 'atol': 1e-2, 'rtol': 1e-2}

    def forward_cpu(self, inputs, running_mean, running_var):
        y_expect = _batch_normalization(
            inputs + [self.mean, self.var, self.eps, self.expander],
            running_mean, running_var, self.decay)
        return y_expect,

    def check_forward(self, inputs, backend_config):
        if self.running_statistics:
            running_mean_expected = self.running_mean.copy()
            running_var_expected = self.running_var.copy()
        else:
            running_mean_expected = None
            running_var_expected = None

        y_expected, = self.forward_cpu(
            inputs, running_mean_expected, running_var_expected)

        inputs = backend_config.get_array(inputs)
        running_mean = backend_config.get_array(self.running_mean)
        running_var = backend_config.get_array(self.running_var)

        if not self.c_contiguous:
            with backend_config:
                inputs = _as_noncontiguous_array(inputs)
                running_mean = _as_noncontiguous_array(running_mean)
                running_var = _as_noncontiguous_array(running_var)

        with backend_config:
            y = functions.batch_normalization(
                *inputs, running_mean=running_mean,
                running_var=running_var, **self.bn_options)
        assert y.data.dtype == self.dtype

        testing.assert_allclose(
            y_expected, y.data, **self.check_forward_options)
        if self.running_statistics:
            testing.assert_allclose(
                running_mean_expected, running_mean,
                **self.check_forward_options)
            testing.assert_allclose(
                running_var_expected, running_var,
                **self.check_forward_options)

    def test_forward(self, backend_config):
        self.check_forward(self.inputs, backend_config)

    def check_backward(self, inputs, grad_outputs, backend_config):
        inputs = backend_config.get_array(inputs)
        grad_outputs = backend_config.get_array(grad_outputs)
        if not self.c_contiguous:
            with backend_config:
                inputs = _as_noncontiguous_array(inputs)
                grad_outputs = _as_noncontiguous_array(grad_outputs)

        def f(*inputs):
            y = functions.batch_normalization(
                *inputs, **self.bn_options)
            return y,

        with backend_config:
            gradient_check.check_backward(
                f, inputs, grad_outputs,
                **self.check_backward_options)

    def test_backward(self, backend_config):
        self.check_backward(self.inputs, self.grad_outputs, backend_config)

    def check_double_backward(
            self, inputs, grad_outputs, grad_grad_inputs, backend_config):
        inputs = backend_config.get_array(inputs)
        grad_outputs = backend_config.get_array(grad_outputs)
        grad_grad_inputs = backend_config.get_array(grad_grad_inputs)
        if not self.c_contiguous:
            with backend_config:
                inputs = _as_noncontiguous_array(inputs)
                grad_outputs = _as_noncontiguous_array(grad_outputs)
                grad_grad_inputs = _as_noncontiguous_array(grad_grad_inputs)

        def f(*inputs):
            return functions.batch_normalization(
                *inputs, **self.bn_options)

        with backend_config:
            gradient_check.check_double_backward(
                f, inputs, grad_outputs, grad_grad_inputs,
                **self.check_double_backward_options)

    def test_double_backward(self, backend_config):
        self.check_double_backward(
            self.inputs, self.grad_outputs, self.grad_grad_inputs,
            backend_config)

@testing.parameterize(*(testing.product({
    'param_shape': [(3,), (3, 4), (3, 2, 3)],
    'ndim': [0, 1, 2],
    'eps': [2e-5, 5e-1],
    'dtype': [numpy.float32],
    'c_contiguous': [True, False],
}) + testing.product({
    'param_shape': [(3,)],
    'ndim': [1],
    'eps': [2e-5, 5e-1],
    'dtype': [numpy.float16, numpy.float32, numpy.float64],
    'c_contiguous': [True, False],
})))
@backend.inject_backend_tests(
    None,
    # CPU tests
    [{'use_cuda': False}]
    # GPU tests
    + testing.product({
        'use_cuda': [True],
        'use_cudnn': ['never', 'always'],
        'cudnn_fast_batch_normalization': [True, False],
    })
    # ChainerX tests
    + [
        {'use_chainerx': True, 'chainerx_device': 'native:0'},
        {'use_chainerx': True, 'chainerx_device': 'cuda:0'},
    ])
class TestFixedBatchNormalization(unittest.TestCase):

    def setUp(self):
        param_shape = self.param_shape
        dtype = self.dtype
        ndim = self.ndim

        gamma = numpy.random.uniform(.5, 1, param_shape).astype(dtype)
        beta = numpy.random.uniform(-1, 1, param_shape).astype(dtype)
        shape = (5,) + param_shape + (2,) * ndim
        x = numpy.random.uniform(-1, 1, shape).astype(dtype)
        mean = numpy.random.uniform(-1, 1, param_shape).astype(dtype)
        var = numpy.random.uniform(0.5, 1, param_shape).astype(dtype)

        gy = numpy.random.uniform(-1, 1, shape).astype(dtype)
        ggx = numpy.random.uniform(-1, 1, shape).astype(dtype)
        gggamma = numpy.random.uniform(-1, 1, param_shape).astype(dtype)
        ggbeta = numpy.random.uniform(-1, 1, param_shape).astype(dtype)
        ggmean = numpy.random.uniform(-1, 1, param_shape).astype(dtype)
        ggvar = numpy.random.uniform(-1, 1, param_shape).astype(dtype)

        self.decay = 0.0
        self.expander = (None, Ellipsis) + (None,) * ndim

        self.inputs = [x, gamma, beta, mean, var]
        self.grad_outputs = [gy]
        self.grad_grad_inputs = [ggx, gggamma, ggbeta, ggmean, ggvar]

        self.check_forward_options = {'atol': 1e-4, 'rtol': 1e-3}
        self.check_backward_options = {'dtype': numpy.float64}
        self.check_double_backward_options = {'dtype': numpy.float64}
        if self.dtype == numpy.float16:
            self.check_forward_options = {'atol': 1e-2, 'rtol': 1e-2}
            self.check_backward_options = {
                'dtype': numpy.float64, 'atol': 1e-2, 'rtol': 1e-2}
            self.check_double_backward_options = {
                'dtype': numpy.float64, 'atol': 1e-2, 'rtol': 1e-2}

    def forward_cpu(self, inputs):
        y_expect = _batch_normalization(inputs + [self.eps, self.expander])
        return y_expect,

    def check_forward(self, inputs, enable_backprop, backend_config):
        y_expected, = self.forward_cpu(inputs)

        inputs = backend_config.get_array(inputs)
        if not self.c_contiguous:
            with backend_config:
                inputs = _as_noncontiguous_array(inputs)

        with chainer.using_config('enable_backprop', enable_backprop):
            with backend_config:
                y = functions.fixed_batch_normalization(*inputs, eps=self.eps)
        assert y.data.dtype == self.dtype

        testing.assert_allclose(
            y_expected, y.data, **self.check_forward_options)

    def test_forward(self, backend_config):
        self.check_forward(self.inputs, False, backend_config)

    def test_forward_with_enable_backprop(self, backend_config):
        self.check_forward(self.inputs, True, backend_config)

    def check_backward(self, inputs, grad_outputs, backend_config):
        inputs = backend_config.get_array(inputs)
        grad_outputs = backend_config.get_array(grad_outputs)
        if not self.c_contiguous:
            with backend_config:
                inputs = _as_noncontiguous_array(inputs)
                grad_outputs = _as_noncontiguous_array(grad_outputs)

        def f(*inputs):
            y = functions.fixed_batch_normalization(*inputs, eps=self.eps)
            return y,

        with backend_config:
            gradient_check.check_backward(
                f, inputs, grad_outputs,
                **self.check_backward_options)

    def test_backward(self, backend_config):
        self.check_backward(self.inputs, self.grad_outputs, backend_config)

    def check_double_backward(
            self, inputs, grad_outputs, grad_grad_inputs, backend_config):
        inputs = backend_config.get_array(inputs)
        grad_outputs = backend_config.get_array(grad_outputs)
        grad_grad_inputs = backend_config.get_array(grad_grad_inputs)
        if not self.c_contiguous:
            with backend_config:
                inputs = _as_noncontiguous_array(inputs)
                grad_outputs = _as_noncontiguous_array(grad_outputs)
                grad_grad_inputs = _as_noncontiguous_array(grad_grad_inputs)

        def f(*inputs):
            return functions.fixed_batch_normalization(*inputs, eps=self.eps)

        with backend_config:
            gradient_check.check_double_backward(
                f, inputs, grad_outputs, grad_grad_inputs,
                **self.check_double_backward_options)

    def test_double_backward(self, backend_config):
        self.check_double_backward(
            self.inputs, self.grad_outputs, self.grad_grad_inputs,
            backend_config)

@testing.parameterize(*testing.product({
    'use_cudnn': ['always', 'auto', 'never'],
    'eps': [2e-5, 5e-1],
    # TODO(bkvogel): Check float16 support again in next cuDNN version.
    'dtype': [numpy.float32, numpy.float64],
}))
@attr.cudnn
class TestBatchNormalizationCudnnCall(unittest.TestCase):

    def setUp(self):
        ndim = 0
        param_shape = (3,)
        self.gamma = cuda.cupy.random.uniform(.5, 1,
                                              param_shape).astype(self.dtype)
        self.beta = cuda.cupy.random.uniform(-1, 1,
                                             param_shape).astype(self.dtype)
        shape = (7,) + param_shape + (2,) * ndim
        self.x = cuda.cupy.random.uniform(-1, 1, shape).astype(self.dtype)
        self.gy = cuda.cupy.random.uniform(-1, 1, shape).astype(self.dtype)
        self.args = [self.x, self.gamma, self.beta]
        head_ndim = self.gamma.ndim + 1
        self.aggr_axes = (0,) + tuple(six.moves.range(head_ndim, self.x.ndim))
        self.mean = self.x.mean(axis=self.aggr_axes)
        self.var = self.x.var(axis=self.aggr_axes) + self.eps
        with chainer.using_config('use_cudnn', self.use_cudnn):
            self.expect = chainer.should_use_cudnn('>=auto', 5000)

    def forward(self):
        return functions.batch_normalization(
            *[chainer.Variable(i) for i in self.args], eps=self.eps,
            running_mean=self.mean, running_var=self.var)

    def test_call_cudnn_forward(self):
        with chainer.using_config('use_cudnn', self.use_cudnn):
            with testing.patch(
                    'cupy.cudnn.batch_normalization_forward_training'
            ) as func:
                self.forward()
                self.assertEqual(func.called, self.expect)

    def test_call_cudnn_backward(self):
        with chainer.using_config('use_cudnn', self.use_cudnn):
            y = self.forward()
            y.grad = self.gy
            with testing.patch(
                    'cupy.cudnn.batch_normalization_backward'
            ) as func:
                y.backward()
                self.assertEqual(func.called, self.expect)

@attr.cudnn
class TestBatchNormalizationCudnnEps(unittest.TestCase):
    def setUp(self):
        ndim = 0
        param_shape = (3,)
        dtype = numpy.float32
        gamma = cuda.cupy.random.uniform(.5, 1, param_shape).astype(dtype)
        beta = cuda.cupy.random.uniform(-1, 1, param_shape).astype(dtype)
        shape = (7,) + param_shape + (2,) * ndim
        x = cuda.cupy.random.uniform(-1, 1, shape).astype(dtype)
        self.args = [x, gamma, beta]

    def test_valid(self):
        functions.batch_normalization(*self.args, eps=1e-5)

    def test_invalid(self):
        with self.assertRaises(RuntimeError):
            functions.batch_normalization(*self.args, eps=2e-6)

@attr.cudnn
class TestFixedBatchNormalizationCudnnEps(unittest.TestCase):
    def setUp(self):
        ndim = 0
        param_shape = (3,)
        dtype = numpy.float32
        gamma = cuda.cupy.random.uniform(.5, 1, param_shape).astype(dtype)
        beta = cuda.cupy.random.uniform(-1, 1, param_shape).astype(dtype)
        mean = cuda.cupy.random.uniform(-1, 1, param_shape).astype(dtype)
        var = cuda.cupy.random.uniform(-1, 1, param_shape).astype(dtype)
        shape = (7,) + param_shape + (2,) * ndim
        x = cuda.cupy.random.uniform(-1, 1, shape).astype(dtype)
        self.args = [x, gamma, beta, mean, var]

    def test_valid(self):
        functions.fixed_batch_normalization(*self.args, eps=1e-5)

    def test_invalid(self):
        with self.assertRaises(RuntimeError):
            functions.fixed_batch_normalization(*self.args, eps=2e-6)

class TestBatchNormalizationWarning(unittest.TestCase):
    def setUp(self):
        pass

    def create_batch(self, param_shape, x_shape):
        dtype = numpy.float32
        gamma = numpy.random.uniform(.5, 1, param_shape).astype(dtype)
        beta = numpy.random.uniform(-1, 1, param_shape).astype(dtype)
        x = numpy.random.uniform(-1, 1, x_shape).astype(dtype)
        args = [x, gamma, beta]
        return args

    def test_invalid_batch(self):
        args = self.create_batch((3,), (1, 3))
        with testing.assert_warns(UserWarning):
            functions.batch_normalization(*args)

    def test_invalid_batch_no_batch_axis(self):
        args = self.create_batch((1, 3,), (1, 3, 1))
        with testing.assert_warns(UserWarning):
            functions.batch_normalization(*args, axis=2)

    def test_valid_batch(self):
        args = self.create_batch((3,), (1, 3, 2, 2))
        with warnings.catch_warnings(record=True) as w:
            functions.batch_normalization(*args)
            assert len(w) == 0

    def test_valid_batch_no_batch_axis(self):
        args = self.create_batch((1, 3,), (1, 3, 2))
        with warnings.catch_warnings(record=True) as w:
            functions.batch_normalization(*args, axis=2)
            assert len(w) == 0

testing.run_module(__name__, __file__)

max_stars_repo_name: denjiry/chainer

max_issues_repo_name: nolfwin/chainer

max_forks_repo_name: nolfwin/chainer


text:
"""
This module contains all the code necessary for the VAE.

The general use is to instantiate a VAE with the appropriate
hyper parameters, then to train it with a dataset. You can
then load the resulting weights into a VAE later.

Much of this code was taken from here: https://blog.keras.io/building-autoencoders-in-keras.html
"""
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import keras
from keras.layers import Lambda, Input, Dense, Conv2D, UpSampling2D, MaxPooling2D, Flatten, Reshape
from keras.models import Model
from keras.datasets import mnist
from keras.losses import mse, binary_crossentropy
from keras import backend as K

import numpy as np
import argparse
import os

def _sampling(args):
    """
    Reparameterization trick by sampling from an isotropic unit Gaussian.

    This function combines `args` (which must be a mean and a log of variance)
    into a single distribution. We then sample from that distribution to
    produce a latent vector.

    :param args:    Two tensors of shape (batch, nembedding_dims) - the first of which
                    is the means and the second of which is the log of the variances.
    :returns:       z: The sampled latent vector.
    """
    z_mean, z_log_var = args
    batch = K.shape(z_mean)[0]
    dim = K.int_shape(z_mean)[1]

    epsilon = K.random_normal(shape=(batch, dim))
    stdev = K.exp(0.5 * z_log_var)

    #z_mean = K.print_tensor(z_mean, message="Z-Mean")
    #z_log_var = K.print_tensor(z_log_var, message="Z Log Var")
    #epsilon = K.print_tensor(epsilon, message="E")
    #stdev = K.print_tensor(stdev, message="STDEV")

    sample = z_mean + stdev * epsilon
    #sample = K.print_tensor(sample, message="Sample")

    return sample

class VariationalAutoEncoder:
    """
    A class to represent an encoder and decoder such that the encoder
    is constrained to learning an embedding of the input vectors
    that conforms to a normal distribution.
    """
    def __init__(self, input_shape, latent_dim, optimizer, loss, *,
                        kl_loss_prop=0.5, recon_loss_prop=0.5, std_loss_prop=0.0, encoder=None, decoder=None, inputlayer=None,
                        decoderinputlayer=None, save_intermediate_models=False, tbdir=None):
        """
        :param input_shape:                 (tuple) The shape of the input data.
        :param latent_dim:                  (int) The dimensionality of the latent vector space.
        :param optimizer:                   String representation of the optimizer.
        :param loss:                        String representation of the loss function.
        :param kl_loss_prop:                Value between 0 and 1.0 that shows how much of the whole VAE loss function to assign to KL loss
        :param recon_loss_prop:             Value between 0 and 1.0 that shows how much of the whole VAE loss function to
                                            assign to reconstructive loss
        :param std_loss_prop:               Value between 0 and 1.0 that shows how much of the whole VAE loss function
                                            to assign to the variance portion
        :param encoder:                     The encoder layer. This must be the result of a sequence of Keras functional calls.
                                            You will also need to pass in the encoder's inputlayer.
                                            If None provided, we use a reasonable default for the MNIST dataset.
                                            This is probably not what you want.
        :param decoder:                     The decoder layer. This must be the result of a sequence of Keras functional calls.
                                            You will also need to pass in decoderinputlayer. The input shape must be (latent_dim,)
                                            and the output shape must be `input_shape`.
                                            If None provided, we use a reasonable default for the MNIST dataset.
                                            This is probably not what you want.
        :param inputlayer:                  Necessary if encoder is not None. This must be a Keras Inputs layer.
        :param decoderinputlayer:           Necessary if decoder is not None. This must be a Keras Inputs layer.
        :param save_intermediate_models:    If `True`, we will save the most recent model after every epoch in a 'models' directory.
        :param tbdir:                       If truthy, must be a directory. This directory will be where we save the tensorboard information.
                                            If None or empty, we will not use tensorboard. To use tensorboard, simply run it from the command
                                            line and pass in --logdir=full/path/to/tbdir
        """
        if encoder is not None and inputlayer is None:
            raise ValueError("If `encoder` is not None, you must also pass in `inputlayer`.")
        if decoder is not None and decoderinputlayer is None:
            raise ValueError("If `decoder` is not None, you must also pass in `decoderinputlayer`.")

        self._latent_dim = latent_dim
        self._encoder, self._inputs, z_mean, z_log_var = self._build_encoder(input_shape, latent_dim, encoder=encoder, inputlayer=inputlayer)
        self._decoder = self._build_decoder(input_shape, latent_dim, decoder=decoder, decoderinputs=decoderinputlayer)
        self._outputs = self._decoder(self._encoder(self._inputs)[2])
        self._vae = Model(self._inputs, self._outputs, name='vae')
        flattened_input_shape = (np.product(np.array(input_shape)),)

        def _vae_loss(y_true, y_pred):
            """
            VAE loss that is broken out as a separate function. It is required to be broken
            out as a separate function for reasons not well understood. See:
            https://github.com/keras-team/keras/issues/10137
            """
            reconstruction_loss = self._build_loss(loss, flattened_input_shape)
            kl_loss = -0.5 * K.sum(1 + z_log_var - K.square(z_mean) - K.exp(z_log_var), axis=-1)
            stddev_loss = K.mean(K.square(z_log_var))
            vae_loss = K.mean((reconstruction_loss * recon_loss_prop) + (kl_loss * kl_loss_prop) + (stddev_loss * std_loss_prop))

            return vae_loss

        self._vae.compile(optimizer=optimizer, loss=_vae_loss)
        self._vae.summary()

        # Do callback initialization stuff
        if not os.path.isdir("models"):
            os.makedirs("models")
        if save_intermediate_models:
            saver = keras.callbacks.ModelCheckpoint("models/vae-weights.ep{epoch:02d}-loss{loss:.4f}.h5", period=1)
            self._callbacks = [saver]
        else:
            self._callbacks = []

        if tbdir:
            tb = keras.callbacks.TensorBoard(log_dir=tbdir,
                                                histogram_freq=0,   # every this many epochs, compute weights hist
                                                batch_size=32,      # Batch size for determining histogram
                                                write_graph=True,   # Should we write the whole network graph?
                                                write_grads=True,   # Gradient histograms
                                                write_images=False, # Not sure
                                                embeddings_freq=0,  # Every this many epochs, pass in embeddings_data to visualize embeddings
                                                update_freq=100)    # Every this many samples, we update tensor board
            self._callbacks.append(tb)

    def sample(self):
        """
        Returns a vector sampled from the latent VAE space. If you want to decode it using predict(),
        you must wrap it in a list:

        ```python
        sample = vae.sample()
        decoded_sample = vae.predict([sample])
        ```
        """
        mu = 0.0
        sigma = 1.0
        size = self._latent_dim
        return np.random.normal(mu, sigma, size)

    def sample_from_gaussian(self, mu: tuple, sigma: tuple):
        """
        Same as `sample` but instead of sampling from the normal distribution, samples from the given
        Gaussian.

        mu and sigma must be tuples of floats, of shape self._latent_dim.
        """
        return np.random.normal(mu, sigma)

    def predict(self, *args, **kwargs):
        """
        Do inference with the VAE. See the Keras predict() method's documentation: https://keras.io/models/model/#predict
        This method runs through the decoder specifically. The encoder is unused.
        To run something through the entire model, use `encode_decode()`.
        """
        return self._decoder.predict(*args, **kwargs)

    def encode_decode(self, x, **kwargs):
        """
        Do inference with the VAE. See the Keras predict() method's documentation: https://keras.io/models/model/#predict
        This method runs through the whole encoder and decoder. To just use the decoder, use `predict()`.
        """
        if len(x.shape) == 2:
            # Likely was passed a numpy array (rows, cols) - but we want a batch of images: (batch, rows, cols, channels)
            s = x.shape
            x = np.reshape(x, (1, s[0], s[1], 1))
        elif len(x.shape) == 3:
            # Who knows what we were given? Reject it and tell the user the error of their ways
            raise ValueError("Need either a numpy array of shape (rows, cols) or a batch of images (batch, rows, cols, channels). But got something with shape: {}".format(x.shape))
        return self._vae.predict(x, **kwargs)

    def fit(self, *args, **kwargs):
        """
        Train the VAE on the given data. See the Keras fit() method's documentation: https://keras.io/models/model/#fit
        """
        if 'callbacks' not in kwargs:
            kwargs['callbacks'] = self._callbacks
        else:
            kwargs['callbacks'].extend(self._callbacks)

        return self._vae.fit(*args, **kwargs)

    def fit_generator(self, datagen, batch_size, save_models=True, **kwargs):
        """
        Train the VAE on the given Sequence (data generator).
        """
        if 'callbacks' not in kwargs:
            kwargs['callbacks'] = self._callbacks
        else:
            kwargs['callbacks'].extend(self._callbacks)

        return self._vae.fit_generator(datagen, **kwargs)

    def load_weights(self, weightfpath):
        """
        Loads the given fpath and loads the weights found in that file.
        """
        self._vae.load_weights(weightfpath)

    def save_weights(self, fpath):
        """
        Saves the weights of this model at the given file path.
        """
        self._vae.save_weights(fpath)

    def _build_loss(self, loss, flattened_input_shape):
        """
        Builds and returns the loss function.
        """
        newshape = (-1, *flattened_input_shape)
        if loss.lower().strip() == "mse":
            reconstruction_loss = mse(K.reshape(self._inputs, newshape), K.reshape(self._outputs, newshape))
        else:
            raise ValueError("{} is not an allowed loss function.".format(loss))
        reconstruction_loss *= flattened_input_shape
        return reconstruction_loss

    def _build_encoder(self, input_shape, latent_dim, encoder=None, inputlayer=None):
        """
        Builds the encoder portion of the model and returns it.

        The default encoder is good for MNIST.
        """
        #                                                                       # MNIST dimensionality
        if encoder is None:
            inputs = Input(shape=input_shape, name="encoder_input")             # (-1, 28, 28, 1)
            x = Conv2D(16, (3, 3), activation='relu', padding='same')(inputs)   # (-1, 28, 28, 16)
            x = MaxPooling2D((2, 2), padding='same')(x)                         # (-1, 14, 14, 16)
            x = Conv2D(8, (3, 3), activation='relu', padding='same')(x)         # (-1, 14, 14, 8)
            x = MaxPooling2D((2, 2), padding='same')(x)                         # (-1, 7, 7, 8)
            x = Conv2D(8, (3, 3), activation='relu', padding='same')(x)         # (-1, 7, 7, 8)
            x = MaxPooling2D((2, 2), padding='same')(x)                         # (-1, 4, 4, 8)
            x = Flatten()(x)                                                    # (-1, 128)
            x = Dense(32, activation='relu')(x)                                 # (-1, 32)
        else:
            x = encoder
            inputs = inputlayer

        z_mean = Dense(latent_dim, name='z_mean')(x)
        z_log_var = Dense(latent_dim, name='z_log_var')(x)

        # use reparameterization trick to push the sampling out as input
        z = Lambda(_sampling, output_shape=(latent_dim,), name='z')([z_mean, z_log_var])  # (-1, 2)

        complete_encoder = Model(inputs, [z_mean, z_log_var, z], name='encoder')
        complete_encoder.summary()

        return complete_encoder, inputs, z_mean, z_log_var

    def _build_decoder(self, input_shape, latent_dim, decoder=None, decoderinputs=None):
        """
        Builds the decoder portion of the model and returns it.
        """
        if decoder is None:
            intermediate_dim = (4, 4, 8)  # hard-coded to MNIST
            latent_inputs = Input(shape=(latent_dim,), name='z_sampling')        # (-1, 2)
            x = Dense(32, activation='relu')(latent_inputs)                      # (-1, 32)
            x = Dense(np.product(intermediate_dim), activation='relu')(x)        # (-1, 128)
            x = Reshape(target_shape=intermediate_dim)(x)                        # (-1, 4, 4, 8)
            x = Conv2D(8, (3, 3), activation='relu', padding='same')(x)          # (-1, 4, 4, 8)
            x = UpSampling2D((2, 2))(x)                                          # (-1, 8, 8, 8)
            x = Conv2D(8, (3, 3), activation='relu', padding='same')(x)          # (-1, 8, 8, 8)
            x = UpSampling2D((2, 2))(x)                                          # (-1, 16, 16, 8)
            x = Conv2D(16, (3, 3), activation='relu')(x)                         # (-1, 14, 14, 16)
            x = UpSampling2D((2, 2))(x)                                          # (-1, 28, 28, 16)
            outputs = Conv2D(1, (3, 3), activation='sigmoid', padding='same')(x) # (-1, 28, 28, 1)
        else:
            latent_inputs = decoderinputs
            outputs = decoder

        # instantiate decoder model
        completed_decoder = Model(latent_inputs, outputs, name='decoder')
        completed_decoder.summary()

        return completed_decoder

if __name__ == '__main__':
    # If this module is run as a script, do a smoke test
    # Further tests are present in the tests directory of this project

    import matplotlib.pyplot as plt

    # Put together argparse
    parser = argparse.ArgumentParser()
    parser.add_argument("-a", "--latentdim", default=2, type=int, help="Dimensionality of the latent vector space")
    parser.add_argument("-b", "--batchsize", default=128, type=int, help="Number of input vectors to batch as a single training instance")
    parser.add_argument("-e", "--nepochs", default=50, type=int, help="Number of passes through the whole dataset for training")
    parser.add_argument("-o", "--optimizer", default="adam", type=str, help="Optimizer to use for training the VAE")
    parser.add_argument("-l", "--loss", default="mse", choices=["mse"], help="The loss function to use along with KL divergence")
    parser.add_argument("-s", "--save", default="vae_mnist.h5", type=str, help="File path to save the weights at. If empty, we will not save the weights.")
    parser.add_argument("-w", "--weights", type=str, help="Load the given file as weights and work in inference mode, rather than training first")
    args = parser.parse_args()

    # Validate the inputs
    if args.weights and not os.path.isfile(args.weights):
        print("Could not find weight file", args.weights)
        exit(1)

    # Get the MNIST data, along with some parameters about it
    (x_train, y_train), (x_test, y_test) = mnist.load_data()
    image_size = x_train.shape[1]
    original_dim = (28, 28, 1)

    # Reshape into the appropriate shapes and types
    x_train = np.reshape(x_train, [-1, *original_dim])
    x_test = np.reshape(x_test, [-1, *original_dim])
    x_train = x_train.astype('float32') / 255.0
    x_test = x_test.astype('float32') / 255.0

    # Create the VAE
    vae = VariationalAutoEncoder(input_shape=original_dim,
                                 latent_dim=args.latentdim,
                                 optimizer=args.optimizer,
                                 loss=args.loss)

    if args.weights:
        vae.load_weights(args.weights)
    else:
        # train the autoencoder
        vae.fit(x_train, x_train, epochs=args.nepochs, batch_size=args.batchsize, validation_data=(x_test, x_test))
        if args.save:
            vae.save_weights(args.save)

    if args.latentdim == 2:
        # Do some inference and plot it for human consumption

        # display a 30x30 2D manifold of digits
        n = 30
        digit_size = 28
        figure = np.zeros((digit_size * n, digit_size * n))

        # linearly spaced coordinates corresponding to the 2D plot
        # of digit classes in the latent space
        grid_x = np.linspace(-4, 4, n)
        grid_y = np.linspace(-4, 4, n)[::-1]

        for i, yi in enumerate(grid_y):
            for j, xi in enumerate(grid_x):
                z_sample = np.array([[xi, yi]])
                x_decoded = vae.predict(z_sample)
                digit = x_decoded[0].reshape(digit_size, digit_size)
                figure[i * digit_size: (i + 1) * digit_size,
                    j * digit_size: (j + 1) * digit_size] = digit

        plt.figure(figsize=(10, 10))
        start_range = digit_size // 2
        end_range = n * digit_size + start_range + 1
        pixel_range = np.arange(start_range, end_range, digit_size)
        sample_range_x = np.round(grid_x, 1)
        sample_range_y = np.round(grid_y, 1)
        plt.xticks(pixel_range, sample_range_x)
        plt.yticks(pixel_range, sample_range_y)
        plt.xlabel("z[0]")
        plt.ylabel("z[1]")
        plt.imshow(figure, cmap='Greys_r')
        plt.show()
    else:
        print("To see decoded samples, run with --latentdim 2")

max_stars_repo_name: MaxStrange/ArtieInfant

max_issues_repo_name: MaxStrange/ArtieInfant

max_forks_repo_name: MaxStrange/ArtieInfant


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

__author__ = "Christian Heider Nielsen"
__doc__ = r"""

           Created on 17-12-2020
           """

from typing import Sequence

import numpy

from apppath import ensure_existence
from draugr.numpy_utilities.datasets.splitting import Split
from draugr.numpy_utilities.signal_utilities.signal_statistics import root_mean_square
from draugr.numpy_utilities.signal_utilities.truncation import min_length_truncate_batch

LOG_MAGNITUDE_MULTIPLIER = 20
LOG_POWER_MULTIPLIER = 10

__all__ = ["mix_ratio"]

def mix_ratio(
    s1: Sequence,
    s2: Sequence,
    db_ratio: float = 0,
    *,
    log_multiplier: int = LOG_MAGNITUDE_MULTIPLIER,
) -> Sequence:
    """
    # Function to mix clean speech and noise at various SNR levels

    # Normalizing to -25 dB FS

    :param s1:
    :param s2:
    :param db_ratio:
    :param log_multiplier:
    :return:"""
    s1, s2 = min_length_truncate_batch((s1, s2))
    s1_max = numpy.abs(s1).max()
    mix = (s1 / s1_max) + (
        (s2 / numpy.abs(s2).max()) * (root_mean_square(s1) / root_mean_square(s2))
    ) / (10 ** (db_ratio / log_multiplier))
    return (mix / numpy.abs(mix).max()) * s1_max

if __name__ == "__main__":

    def asad():
        """ """
        from neodroidaudition.data.recognition.libri_speech import LibriSpeech
        from neodroidaudition.noise_generation.gaussian_noise import white_noise

        from pathlib import Path

        libri_speech = LibriSpeech(
            path=Path.home() / "Data" / "Audio" / "Speech" / "LibriSpeech"
        )
        files, sr = zip(*[(v[0].numpy(), v[1]) for _, v in zip(range(1), libri_speech)])
        assert all([sr[0] == s for s in sr[1:]])

        mixed = mix_ratio(files[0], white_noise(files[0].shape[-1]))
        print(mixed)

    def asadsa():
        """ """
        from draugr.torch_utilities import to_tensor
        from neodroidaudition.data.recognition.libri_speech import LibriSpeech
        from neodroidaudition.noise_generation.gaussian_noise import white_noise
        import torchaudio

        from pathlib import Path

        libri_speech = LibriSpeech(
            path=Path.home() / "Data" / "Audio" / "Speech" / "LibriSpeech",
            split=Split.Testing,
        )
        files, sr = zip(
            *[(v[0].numpy(), v[1]) for _, v in zip(range(20), libri_speech)]
        )
        assert all([sr[0] == s for s in sr[1:]])

        mix = files[0]
        for file in files[1:]:
            mix = mix_ratio(mix, file, 0)

        torchaudio.save(
            str(ensure_existence(Path.cwd() / "exclude") / f"mixed_even_babble.wav"),
            to_tensor(mix),
            int(sr[0]),
        )

        for ratio in range(-20, 20 + 1, 5):
            torchaudio.save(
                str(ensure_existence(Path.cwd() / "exclude") / f"mixed_{ratio}.wav"),
                to_tensor(mix_ratio(files[0], files[-1], ratio)),
                int(sr[0]),
            )

    def asadsa2():
        """ """
        from draugr.torch_utilities import to_tensor
        from neodroidaudition.data.recognition.libri_speech import LibriSpeech
        from neodroidaudition.noise_generation.gaussian_noise import white_noise
        import torchaudio
        from pathlib import Path

        libri_speech = LibriSpeech(
            path=Path.home() / "Data" / "Audio" / "Speech" / "LibriSpeech"
        )
        files, sr = zip(*[(v[0].numpy(), v[1]) for _, v in zip(range(1), libri_speech)])
        assert all([sr[0] == s for s in sr[1:]])

        normed = files[0]
        mixed = mix_ratio(normed, normed, 0)
        mixed2 = mix_ratio(mixed, mixed, 0)
        print(normed, mixed)
        print(mixed2, mixed)
        print(root_mean_square(normed))
        print(root_mean_square(mixed))
        print(root_mean_square(mixed2))
        assert numpy.allclose(normed, mixed)
        assert numpy.allclose(mixed2, mixed)
        torchaudio.save(
            str(ensure_existence(Path.cwd() / "exclude") / "mixed_same.wav"),
            to_tensor(mixed),
            int(sr[0]),
        )

    # asad()
    asadsa()
    # asadsa2()

max_stars_repo_name: cnHeider/draugr

max_issues_repo_name: cnHeider/draugr

max_forks_repo_name: cnHeider/draugr