text:
# -*- coding: utf-8 -*-
"""
Created on Sun Jan 14 23:43:12 2018

Note that:
    This code is modified from the minresQLP algorithms:
        http://www.stanford.edu/group/SOL/software.html
    Edited to suit Newton-MR methods.

Contact:
    yang.liu2(AT)uq.edu.au

REFERENCES:
    Sou-Cheng T. Choi, Christopgher C. Paige, and Michael A. Saunders,
    MINRES-QLP: A Krylov subspace method for indefinite or singular symmetric
    systems, SIAM Journal of Scientific Computing, submitted on March 7, 2010.

    Sou-Cheng T. Choi's PhD Dissertation, Stanford University, 2006:
        http://www.stanford.edu/group/SOL/dissertations.html

    Yang Liu, Fred Roosta
    Stability Analysis of Newton-MR Under Hessian Perturbations,
        https://arxiv.org/abs/1909.06224

    Fred Roosta, Yang Liu, Peng Xu, Michael W. Mahoney
    Newton-MR: Newton's Method Without Smoothness or Convexity
        https://arxiv.org/abs/1810.00303

--------------------------------------------------------------------------
MinresQLP: Aim to obtain the min-length solution of symmetric
   (possibly singular) Ax=b or min||Ax-b||.

   X = MinresQLP(A,B) solves the system of linear equations A*X=B
   or the least-squares problem min norm(B-A*X) if A is singular.
   The N-by-N matrix A must be symmetric or Hermitian, but need not be
   positive definite or nonsingular.  It may be double or single.
   The rhs vector B must have length N.  It may be real or complex,
   double or single,

   X = MinresQLP(AFUN,B) accepts a function handle AFUN instead of
   the matrix A.  Y = AFUN(X) returns the matrix-vector product Y=A*X.
   In all of the following syntaxes, A can be replaced by AFUN.

   X = MinresQLP(A,B,RTOL) specifies a stopping tolerance.
   If RTOL=[] or is absent, a default value is used.
   (Similarly for all later input parameters.)
   Default RTOL=1e-6.

   X = MinresQLP(A,B,RTOL,MAXIT)
   specifies the maximum number of iterations.  Default MAXIT=N.

   X = MinresQLP(A,B,RTOL,MAXIT,M)
   uses a matrix M as preconditioner.  M must be positive definite
   and symmetric or Hermitian.  It may be a function handle MFUN
   such that Y=MFUN(X) returns Y=M\X.
   If M=[], a preconditioner is not applied.

   X = MinresQLP(A,B,RTOL,MAXIT,M,SHIFT)
   solves (A - SHIFT*I)X = B, or the corresponding least-squares problem
   if (A - SHIFT*I) is singular, where SHIFT is a real or complex scalar.
   Default SHIFT=0.

   X = MinresQLP(A,B,RTOL,MAXIT,M,SHIFT,MAXXNORM,ACONDLIM,TRANCOND)
   specifies three parameters associated with singular or
   ill-conditioned systems (A - SHIFT*I)*X = B.

   MAXXNORM is an upper bound on NORM(X).
   Default MAXXNORM=1e7.

   ACONDLIM is an upper bound on ACOND, an estimate of COND(A).
   Default ACONDLIM=1e15.

   TRANCOND is a real scalar >= 1.
   If TRANCOND>1,        a switch is made from MINRES iterations to
                         MINRES-QLP iterationsd when ACOND >= TRANCOND.
   If TRANCOND=1,        all iterations will be MINRES-QLP iterations.
   If TRANCOND=ACONDLIM, all iterations will be conventional MINRES
                         iterations (which are slightly cheaper).
   Default TRANCOND=1e7.

   FLAG:
   -1 (beta2=0)  B and X are eigenvectors of (A - SHIFT*I).
    0 (beta1=0)  B = 0.  The exact solution is X = 0.
    1 X solves the compatible/incompatible (possibly singular) system
    (A - SHIFT*I)X = B to the desired inexactness:
         <B, AX> < - ( 1 - rtol) ||B||^2
    2 X converged to an eigenvector of (A - SHIFT*I).
    3 XNORM exceeded MAXXNORM.
    4 ACOND exceeded ACONDLIM.
    5 MAXIT iterations were performed before one of the previous
      conditions was satisfied.
    9 The system appears to be exactly singular.  XNORM does not
      yet exceed MAXXNORM, but would if further iterations were
      performed.

    ITER:  the number of iterations performed.
    QLPITER:  the number of MINRES-QLP iterations.

    RELRES: Relative residuals for (A - SHIFT*I)X = B
    ANORM:  an estimate of the 2-norm of A-SHIFT*I.
    ACOND:  an estimate of COND(A-SHIFT*I,2).
    XNORM:  a recurred estimate of NORM(X).
"""

import numpy as np
import torch
from numpy.linalg import norm
from myCG import myCG

def MinresQLP(A, b, rtol=None, maxit=None, M=None, shift=None,maxxnorm=None,
              Acondlim=None, TranCond=None):
    if rtol is None:
        rtol = 1e-4
    if maxit is None:
        maxit = 100
    if shift is None:
        shift = 0
    if maxxnorm is None:
        maxxnorm = 1e7
    if Acondlim is None:
        Acondlim = 1e15
    if TranCond is None:
        TranCond = 1e7
    n = len(b)
    device = b.device
    # rtol = rtol*shift/2
    x0 = torch.zeros(n)
    x = x0.clone()
    Ab = Ax(A, b)
    r2 = b
    r3 = r2
    beta1 = torch.norm(r2)

    #function handle with M x r_hat = r
    if M is None:
        #test with M = lambda u: u
        noprecon = True
        pass
    else:
        noprecon = False
        r3 = Precond(M, r2)
        beta1 = torch.dot(r3, r2) #theta
        if beta1 <0:
            print('Error: "M" is indefinite!')
        else:
            beta1 = np.sqrt(beta1)

    ## Initialize
    flag0 = -2
    flag = -2
    iters = 0
    QLPiter = 0
    beta = 0
    tau = 0
    taul = 0
    phi = beta1
    betan = beta1
    gmin = 0
    cs = -1
    sn = 0
    cr1 = -1
    sr1 = 0
    cr2 = -1
    sr2 = 0
    dltan = 0
    eplnn = 0
    gama = 0
    gamal = 0
    gamal2 = 0
    eta = 0
    etal = 0
    etal2 = 0
    vepln = 0
    veplnl = 0
    veplnl2 = 0
    ul3 = 0
    ul2 = 0
    ul = 0
    u = 0
    rnorm = betan
    xnorm = 0
    xl2norm = 0
    Anorm = 0
    Acond = 1
#    relres = rnorm / (beta1 + 1e-50)
    relres = 1
#    x = np.zeros((n,1))
    w = torch.zeros(n)
    wl = torch.zeros(n)
    rk = b #new

    #b = 0 --> x = 0 skip the main loop
    while flag == flag0 and iters < maxit:
        #lanczos
        iters += 1
        betal = beta
        beta = betan
        v = r3/beta
        r3 = Ax(A, v)
        if shift == 0:
            pass
        else:
            r3 = r3 - shift*v

        alfa = torch.dot(r3, v)

        if iters > 1:
            r3 = r3 - r1*beta/betal

        r3 = r3 - r2*alfa/beta
        r1 = r2
        r2 = r3

        if noprecon:
            betan = norm(r3)
            if iters == 1:
                if betan == 0:
                    if alfa == 0:
                        flag = 0
                        print('WARNNING: flag = 0')
                        break
                    else:
                        flag = -1
                        print('WARNNING: flag = -1')
                        # Probbaly lost all the info, x=0 is true solution
                        x = b/alfa
                        break
        else:
            r3 = Precond(M, r2)
            betan = torch.dot(r3, r2)
            if betan > 0:
                betan = np.sqrt(betan)
            else:
                print('Error: "M" is indefinite or singular!')

        pnorm = np.sqrt(betal ** 2 + alfa ** 2 + betan ** 2)
#        vk = np.copy(v) #new
        vnew = r3/betan #new

        ## construct Lanczos V & T matrices s.t., AV = VT
        # if iters == 1:
        #     V = v
        #     alpha_diag = alfa
        #     beta_diag = betan
        #     Vn = np.append(V, vnew, axis=1)
        #     T = alfa
        # else:
        #     alpha_diag = np.append(alpha_diag, alfa)
        #     beta_diag = np.append(beta_diag, betan)
        #     V = Vn
        #     Vn = np.append(V, vnew, axis=1)
        #     T = np.diag(beta_diag[:-1],-1) + np.diag(
        #             alpha_diag,0) + np.diag(beta_diag[:-1],1)
        # t1 = np.zeros((1,iters))
        # t1[0,-1] = betan
        # T = np.append(T, t1, axis=0)

        # if iters == 11:
        #     print(' ')
        #previous left rotation Q_{k-1}
        dbar = dltan
        dlta = cs*dbar + sn*alfa
        epln = eplnn
        gbar = sn*dbar - cs*alfa
        eplnn = sn*betan
        dltan = -cs*betan
        dlta_QLP = dlta
        #current left plane rotation Q_k
        gamal3 = gamal2
        gamal2 = gamal
        gamal = gama
        oldcs = cs #new
        cs, sn, gama = SymGivens(gbar, betan)
        gama_tmp = gama
        taul2 = taul
        taul = tau
        tau = cs*phi
        oldphi = phi #new
        phi = sn*phi
        #previous right plane rotation P_{k-2,k}
        if iters > 2:
            veplnl2 = veplnl
            etal2 = etal
            etal = eta
            dlta_tmp = sr2*vepln - cr2*dlta
            veplnl = cr2*vepln + sr2*dlta
            dlta = dlta_tmp
            eta = sr2*gama
            gama = -cr2 *gama
        #current right plane rotation P{k-1,k}
        if iters > 1:
            cr1, sr1, gamal = SymGivens(gamal, dlta)
            vepln = sr1*gama
            gama = -cr1*gama

        ul4 = ul3
        ul3 = ul2
        if iters > 2:
            ul2 = (taul2 - etal2*ul4 - veplnl2*ul3)/gamal2
        if iters > 1:
            ul = (taul - etal*ul3 - veplnl *ul2)/gamal
        xnorm_tmp = np.sqrt(xl2norm**2 + ul2**2 + ul**2)
        if abs(gama) > np.finfo(np.double).tiny and xnorm_tmp < maxxnorm:
            u = (tau - eta*ul2 - vepln*ul)/gama
            if np.sqrt(xnorm_tmp**2 + u**2) > maxxnorm:
                u = 0
                flag = 3
                print('WARNNING: flag = 3')
        else:
            u = 0
            flag = 6
        xl2norm = np.sqrt(xl2norm**2 + ul2**2)
        xnorm = np.sqrt(xl2norm**2 + ul**2 + u**2)
        #update w&x
        #Minres
        if (Acond < TranCond) and flag != flag0 and QLPiter == 0:
            wl2 = wl
            wl = w
            w = (v - epln*wl2 - dlta_QLP*wl)/gama_tmp
            if xnorm < maxxnorm:
                x += tau*w
            else:
                flag = 3
                print('WARNNING: flag = 3')
        #Minres-QLP
        else:
            QLPiter += 1
            if QLPiter == 1:
                xl2 = torch.zeros(n)
                if (iters > 1):  # construct w_{k-3}, w_{k-2}, w_{k-1}
                    if iters > 3:
                        wl2 = gamal3*wl2 + veplnl2*wl + etal*w
                    if iters > 2:
                        wl = gamal_QLP*wl + vepln_QLP*w
                    w = gama_QLP*w
                    xl2 = x - wl*ul_QLP - w*u_QLP

            if iters == 1:
                wl2 = wl
                wl = v*sr1
                w = -v*cr1
            elif iters == 2:
                wl2 = wl
                wl = w*cr1 + v*sr1
                w = w*sr1 - v*cr1
            else:
                wl2 = wl
                wl = w
                w = wl2*sr2 - v*cr2
                wl2 = wl2*cr2 +v*sr2
                v = wl*cr1 + w*sr1
                w = wl*sr1 - w*cr1
                wl = v
            xl2 = xl2 + wl2*ul2
            oldx = x #new
            x = xl2 + wl*ul + w*u
            # new
            oldrk = rk
            rk = sn**2 * oldrk - phi * cs * vnew
#            print('pTr', -rk.T.dot(x)<0, -rk.T.dot(x))
#            print('pTg', -b.T.dot(x)<0, -b.T.dot(x)) # <0, >0#
#            print(' ')
#            print('r', norm(b - Ax(A,x) - rk))
#            oldArk = oldphi * (gbar * vk + dltan * vnew )
#            print('Ar', norm(Ax(A, b - Ax(A,oldx)) - oldArk))
#    #        print(norm(oldArk), norm(A(b - matvec(oldx))), Arnorm)
#            oldrAr = oldrk.T.dot(oldArk)
#            oldrk2 = b - Ax(A,oldx)
#            oldrAr2 = oldrk2.T.dot(Ax(A,oldrk2))
#            nc = - oldphi**2 * oldcs * gbar
#            if nc/norm(oldrk)**2 < -shift/2:
#                print('nc',iters, nc)
#                x = oldx
#            print('rAr', oldrAr, oldrAr2, nc)
#            print('iters', iters)

        #next right plane rotation P{k-1,k+1}
        gamal_tmp = gamal
        cr2, sr2, gamal = SymGivens(gamal, eplnn)
        #transfering from Minres to Minres-QLP
        gamal_QLP = gamal_tmp
        vepln_QLP = vepln
        gama_QLP = gama
        ul_QLP = ul
        u_QLP = u
        ## Estimate various norms
        abs_gama = abs(gama)
        Anorm = max([Anorm, pnorm, gamal, abs_gama])
        if iters == 1:
            gmin = gama
            gminl = gmin
        elif iters > 1:
            gminl2 = gminl
            gminl = gmin
            gmin = min([gminl2, gamal, abs_gama])
        Acondl = Acond
        Acond = Anorm / gmin
        rnorml = rnorm
        relresl = relres
        if flag != 9:
            rnorm = phi

        # inexactness of Newton-MR
        # relres = 1 - torch.dot(Ab, x)/beta1**2 # <g, Hp> < - ( 1 - rtol) ||g||^2
        relres = rnorm / beta1 # relative residual ||Hp+g|| < rtol ||g||

        ## See if any of the stopping criteria are satisfied.
        epsx = Anorm * xnorm * np.finfo(float).eps
        if (flag == flag0) or (flag == 6):
            if iters >= maxit:
                flag = 5 #exit before maxit
            if Acond >= Acondlim:
                flag = 4 #Huge Acond
                print('WARNNING: Acondlim exceeded!')
            if xnorm >= maxxnorm:
                flag = 3 #xnorm exceeded
                print('WARNNING: maxxnorm exceeded!')
            if epsx >= beta1:
                flag = 2 #x = eigenvector
            if relres <= np.sqrt(rtol):
                flag = 1 #Trustful Solution
        if flag == 3 or flag == 4:
            print('WARNNING: possibly singular!')
            #possibly singular
            iters = iters - 1
            Acond = Acondl
            rnorm = rnorml
            relres = relresl

#    print('betan', betan)
    return x,relres,iters

def modified_MinresQLP(A, b, rtol=None, maxit=None, M=None, shift=None,
                       maxxnorm=None, Acondlim=None, TranCond=None):
    """
    Modifed MINRES-QLP, modified to be able to return Negative curvature
    direction & stopping condition is changge to satisfy Negative curvature
    Line search
    """
    if rtol is None:
        rtol = 1e-4
    if maxit is None:
        maxit = 100
    if shift is None:
        shift = 0
    if maxxnorm is None:
        maxxnorm = 1e7
    if Acondlim is None:
        Acondlim = 1e15
    if TranCond is None:
        TranCond = 1e7
    n = len(b)
    x0 = np.zeros((n,1))
    x = np.copy(x0)
    b = b.reshape(n,1)
    Ab = Ax(A, b)
    r2 = b
    r3 = r2
    beta1 = norm(r2)

    #function handle with M x r_hat = r
    if M is None:
        #test with M = lambda u: u
        noprecon = True
        pass
    else:
        noprecon = False
        r3 = Precond(M, r2)
        beta1 = np.matmul(r3.T,r2) #theta
        if beta1 <0:
            print('Error: "M" is indefinite!')
        else:
            beta1 = np.sqrt(beta1)

    ## Initialize
    flag0 = -2
    flag = -2
    iters = 0
    QLPiter = 0
    beta = 0
    tau = 0
    taul = 0
    phi = beta1
    betan = beta1
    gmin = 0
    cs = -1
    sn = 0
    cr1 = -1
    sr1 = 0
    cr2 = -1
    sr2 = 0
    dltan = 0
    eplnn = 0
    gama = 0
    gamal = 0
    gamal2 = 0
    eta = 0
    etal = 0
    etal2 = 0
    vepln = 0
    veplnl = 0
    veplnl2 = 0
    ul3 = 0
    ul2 = 0
    ul = 0
    u = 0
    rnorm = betan
    xnorm = 0
    xl2norm = 0
    Anorm = 0
    Acond = 1
#    relres = rnorm / (beta1 + 1e-50)
    relres = 1
#    x = np.zeros((n,1))
    w = np.zeros((n,1))
    wl = np.zeros((n,1))
    rk = b
    #b = 0 --> x = 0 skip the main loop
    d_type = 'SOL'
    while flag == flag0 and iters < maxit:
        #lanczos
        iters += 1
        betal = beta
        beta = betan
        v = r3/beta
        r3 = Ax(A, v)
        if shift == 0:
            pass
        else:
            r3 = r3 - shift*v

        alfa = np.real(np.matmul(r3.T,v))

        if iters > 1:
            r3 = r3 - r1*beta/betal

        r3 = r3 - r2*alfa/beta
        r1 = r2
        r2 = r3

        if noprecon:
            betan = norm(r3)
            if iters == 1:
                if betan == 0:
                    if alfa == 0:
                        flag = 0
                        print('WARNNING: flag = 0')
                        break
                    else:
                        flag = -1
                        print('WARNNING: flag = -1')
                        # Probbaly lost all the info, x=0 is true solution
                        x = b/alfa
                        break
        else:
            r3 = Precond(M, r2)
            betan = np.matmul(r2.T, r3)
            if betan > 0:
                betan = np.sqrt(betan)
            else:
                print('Error: "M" is indefinite or singular!')

        vk = v
        vnew = r3/betan
        pnorm = np.sqrt(betal ** 2 + alfa ** 2 + betan ** 2)
        #previous left rotation Q_{k-1}
        dbar = dltan
        dlta = cs*dbar + sn*alfa
        epln = eplnn
        gbar = sn*dbar - cs*alfa
        eplnn = sn*betan
        dltan = -cs*betan
        dlta_QLP = dlta
        #current left plane rotation Q_k
        gamal3 = gamal2
        gamal2 = gamal
        gamal = gama
        oldcs = cs
        cs, sn, gama = SymGivens(gbar, betan)
        gama_tmp = gama
        taul2 = taul
        taul = tau
        tau = cs*phi
        oldphi = phi
        phi = sn*oldphi
        #previous right plane rotation P_{k-2,k}
        if iters > 2:
            veplnl2 = veplnl
            etal2 = etal
            etal = eta
            dlta_tmp = sr2*vepln - cr2*dlta
            veplnl = cr2*vepln + sr2*dlta
            dlta = dlta_tmp
            eta = sr2*gama
            gama = -cr2 *gama
        #current right plane rotation P{k-1,k}
        if iters > 1:
            cr1, sr1, gamal = SymGivens(gamal, dlta)
            vepln = sr1*gama
            gama = -cr1*gama

        ul4 = ul3
        ul3 = ul2
        if iters > 2:
            ul2 = (taul2 - etal2*ul4 - veplnl2*ul3)/gamal2
        if iters > 1:
            ul = (taul - etal*ul3 - veplnl *ul2)/gamal
        xnorm_tmp = np.sqrt(xl2norm**2 + ul2**2 + ul**2)
        if abs(gama) > np.finfo(np.double).tiny and xnorm_tmp < maxxnorm:
            u = (tau - eta*ul2 - vepln*ul)/gama
            if np.sqrt(xnorm_tmp**2 + u**2) > maxxnorm:
                u = 0
                flag = 3
                print('WARNNING: flag = 3')
        else:
            u = 0
            flag = 6
        xl2norm = np.sqrt(xl2norm**2 + ul2**2)
        xnorm = np.sqrt(xl2norm**2 + ul**2 + u**2)
        #update w&x
        #Minres
        if (Acond < TranCond) and flag != flag0 and QLPiter == 0:
            wl2 = wl
            wl = w
            w = (v - epln*wl2 - dlta_QLP*wl)/gama_tmp
            if xnorm < maxxnorm:
                x += tau*w
            else:
                flag = 3
                print('WARNNING: flag = 3')
        #Minres-QLP
        else:
            QLPiter += 1
            if QLPiter == 1:
                xl2 = np.zeros((n,1))
                if (iters > 1):  # construct w_{k-3}, w_{k-2}, w_{k-1}
                    if iters > 3:
                        wl2 = gamal3*wl2 + veplnl2*wl + etal*w
                    if iters > 2:
                        wl = gamal_QLP*wl + vepln_QLP*w
                    w = gama_QLP*w
                    xl2 = x - wl*ul_QLP - w*u_QLP

            if iters == 1:
                wl2 = wl
                wl = v*sr1
                w = -v*cr1
            elif iters == 2:
                wl2 = wl
                wl = w*cr1 + v*sr1
                w = w*sr1 - v*cr1
            else:
                wl2 = wl
                wl = w
                w = wl2*sr2 - v*cr2
                wl2 = wl2*cr2 +v*sr2
                v = wl*cr1 + w*sr1
                w = wl*sr1 - w*cr1
                wl = v
            xl2 = xl2 + wl2*ul2
            oldx = x #new
            x = xl2 + wl*ul + w*u
#            print(f_nc(A, x, shift))
            # new
            oldrk = rk
            rk = sn**2 * oldrk - phi * cs * vnew
#            print('r', norm(b - Ax(A,x) - rk))
#            oldArk = oldphi * (gbar * vk + dltan * vnew )
#            print('Ar', norm(Ax(A, b - Ax(A,oldx)) - oldArk))
    #        print(norm(oldArk), norm(A(b - matvec(oldx))), Arnorm)
#            oldrAr = oldrk.T.dot(oldArk)
#            oldrk2 = b - Ax(A,oldx)
#            oldrAr2 = oldrk2.T.dot(Ax(A,oldrk2))
            nc = - oldphi**2 * oldcs * gbar
#            print('rAr', oldrAr, nc, norm(oldrk)**2, nc/norm(oldrk)**2, -shift/2)
#            print('iters', iters)
            if nc/norm(oldrk)**2 < -shift/2:
#                x = oldx
                d_type = 'NC'
                print('perp', f_nc(A, oldrk, shift), iters)
#                if iters > 1:
#                    p_perp = oldrk.T.dot(oldx2)/norm(oldx2)**2*oldx2
#                    p_nc = oldrk - p_perp
#    #                    tr = oldrk + tr_perp
#                    print('perp', f_nc(A, oldrk, shift), f_nc(A, p_nc, shift), f_nc(A, p_perp, shift), iters)
#                if T_nc == 1:
#                    tr = oldrk
#                    R_nc = oldrk
#                else:
##                    R_nc = np.append(R_nc, oldrk, axis=1)
##                    r_nc = oldrk + r_nc
#    #                print(norm(x0))
#    #                print('rAr', oldrAr2/norm(oldrk2)**2, oldrAr2)
#    #                x0 = r9
##                    oldx_perp = oldx - oldx.T.dot(oldrk)/norm(oldrk)**2*oldrk
##                    tr = oldrk + tr_perp
##                    print('perp', f_nc(A, oldrk), f_nc(A, p_nc))
##                    rAr_new, lamda_new = f_nc(A, tr)
##                    print('combine', rAr_new, lamda_new)
#                    print(' ')
##                return x, r_nc, nc, iters, d_type
                return oldx, oldrk, nc, iters, d_type
##                print(norm(tr))
#            else:
#                T_pc = T_pc + 1
#                if T_pc == 1:
#                    R_pc = oldrk
#                else:
#                    R_pc = np.append(R_pc, oldrk, axis=1)

        #next right plane rotation P{k-1,k+1}
        gamal_tmp = gamal
        cr2, sr2, gamal = SymGivens(gamal, eplnn)
        #transfering from Minres to Minres-QLP
        gamal_QLP = gamal_tmp
        vepln_QLP = vepln
        gama_QLP = gama
        ul_QLP = ul
        u_QLP = u
        ## Estimate various norms
        abs_gama = abs(gama)
        Anorm = max([Anorm, pnorm, gamal, abs_gama])
        if iters == 1:
            gmin = gama
            gminl = gmin
        elif iters > 1:
            gminl2 = gminl
            gminl = gmin
            gmin = min([gminl2, gamal, abs_gama])
        Acondl = Acond
        Acond = Anorm / gmin
        rnorml = rnorm
        relresl = relres
        if flag != 9:
            rnorm = phi

        # inexactness of Newton-MR
#        relres = 1 - Ab.T.dot(x)/beta1**2 # <g, Hp> < - ( 1 - rtol) ||g||^2
        relres = rnorm / beta1 # relative residual ||Hp+g|| < rtol ||g||
        print('rtol',relres - rtol)
#        relres = -rnorm/xnorm/shift*2
#        print(relres)

        ## See if any of the stopping criteria are satisfied.
        epsx = Anorm * xnorm * np.finfo(float).eps
        if (flag == flag0) or (flag == 6):
            if iters >= maxit:
                flag = 5 #exit before maxit
            if Acond >= Acondlim:
                flag = 4 #Huge Acond
                print('WARNNING: Acondlim exceeded!')
            if xnorm >= maxxnorm:
                flag = 3 #xnorm exceeded
                print('WARNNING: maxxnorm exceeded!')
            if epsx >= beta1:
                flag = 2 #x = eigenvector
            if relres <= rtol:
                flag = 1 #Trustful Solution
        if flag == 3 or flag == 4:
            print('WARNNING: possibly singular!')
            #possibly singular
            iters = iters - 1
            Acond = Acondl
            rnorm = rnorml
            relres = relresl

#    for i in range(T_pc):
#        rrr = R_pc[:,i].reshape(-1, 1)
#        x_pc = x_pc + x.T.dot(rrr)/norm(rrr)**2*rrr
#        print('T_pc', i)
#
#    for i in range(T_nc):
#        rrr = R_nc[:,i].reshape(-1, 1)
#        x_nc = x_nc + x.T.dot(rrr)/norm(rrr)**2*rrr
#        print('T_nc', i)

#    return x,relres,iters
    return x, x, alfa, iters, d_type
#    return x, x_pc, x_nc, alfa, iters, d_type

def Ax(A, x):
    if callable(A):
        Ax = A(x)
    else:
        Ax =torch.mv(A, x)
    return Ax

def Precond(M, r):
    if callable(M):
        h = myCG(M, r, 10e-6, 1000)[0]
#        h = cg(M, r)
    else:
        h = torch.mv(M, r)
    return h

def f_nc(A, v, shift=0):
    nc = v.T.dot(Ax(A, v)) - shift*v.T.dot(v)
    lamda = nc/(v.T.dot(v))
    return nc, lamda

def SymGivens(a, b, device="cpu"):
    """This is used by the function minresQLP.

    Arguments:
        a (float): A number.
        b (float): A number.
        device (torch.device): The device tensors will be allocated to.
    """
    if not torch.is_tensor(a):
        a = torch.tensor(float(a), device=device)
    if not torch.is_tensor(b):
        b = torch.tensor(float(b), device=device)
    if b == 0:
        if a == 0:
            c = 1
        else:
            c = torch.sign(a)
        s = 0
        r = torch.abs(a)
    elif a == 0:
        c = 0
        s = torch.sign(b)
        r = torch.abs(b)
    elif torch.abs(b) > torch.abs(a):
        t = a / b
        s = torch.sign(b) / torch.sqrt(1 + t ** 2)
        c = s * t
        r = b / s
    else:
        t = b / a
        c = torch.sign(a) / torch.sqrt(1 + t ** 2)
        s = c * t
        r = a / c
    return c, s, r

def main():
    """Run an example of minresQLP."""
    from time import time
    # torch.manual_seed(1)
    # device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")
    device = torch.device("cpu")
    n = 5
    A = torch.randn(n, n, device=device)
    b = torch.randn(n, device=device)
    # print(b)
    # x = torch.dot(b,b)
    # x = torch.mm(b.T, b)
    x, relres, iters= MinresQLP(A + A.T, b, 1E-2, 100)
    print(x)

if __name__ == '__main__':
    main()

max_stars_repo_name: syangliu/CUTEst-test

max_issues_repo_name: syangliu/CUTEst-test

max_forks_repo_name: syangliu/CUTEst-test

reason: import numpy,from numpy


text:
import os
import numpy as np
from tqdm import tqdm
import tensorflow as tf

epochs = 25
batch_size = 4
data_zh = "data/zh.tsv"

# 数据处理&读取数据
class DealData(object):

    # 初始化
    def __init__(self, filename=data_zh):

        with open(filename, 'r', encoding='utf-8') as fout:
            data = fout.readlines()[:100]

        self.__input_data = []
        self.__label_data = []
        for i in tqdm(range(len(data))):
            key, pny, hanzi = data[i].split('\t')
            self.__input_data.append(pny.split(' '))
            self.__label_data.append(hanzi.strip('\n').split(' '))

        print('文件读取正常----->>>1')

    # 获取词汇
    def get_vocab(self, data):
        print('词汇获取正常----->>>2')
        vocab = ['<PAD>']
        for line in tqdm(data):
            for char in line:
                if char not in vocab:
                    vocab.append(char)
        return vocab

    # 获取输入值与标签值
    def pyn_han(self):
        pny2id = self.get_vocab(self.__input_data)  # 获取输入值
        han2id = self.get_vocab(self.__label_data)  # 获取标签值
        print('PNY:', pny2id[:10])
        print('HAN:', han2id[:10])
        return pny2id, han2id

    # 获取输入数量与标签数量
    def input_lable_num(self):

        print('输入与输出获取正常----->>>3')
        pny2id, han2id = self.pyn_han()
        input_num = [[pny2id.index(pny) for pny in line] for line in tqdm(self.__input_data)]
        label_num = [[han2id.index(han) for han in line] for line in tqdm(self.__label_data)]
        return input_num, label_num

    # 文件分批读取
    def get_batch(self, input_data, label_data, batch_size):
        print('文件分批读取正常----->>>4')
        batch_num = len(input_data) // batch_size
        for k in range(batch_num):
            begin = k * batch_size
            end = begin + batch_size
            input_batch = input_data[begin:end]
            label_batch = label_data[begin:end]
            max_len = max([len(line) for line in input_batch])
            input_batch = np.array([line + [0] * (max_len - len(line)) for line in input_batch])
            label_batch = np.array([line + [0] * (max_len - len(line)) for line in label_batch])
            yield input_batch, label_batch

    # # input与labelbatch
    # def input_label_batch(self):
    #
    #     input_num, label_num = self.input_lable_num()
    #     batch = self.get_batch(input_num, label_num, 4)
    #     input_batch, label_batch = next(batch)
    #     # print(input_batch)
    #     # print(label_batch)

    # 参数设定
    def create_hparams(self):
        print('参数设定正常----->>>5')
        params = tf.contrib.training.HParams(
            num_heads=8,
            num_blocks=6,
            # vocab
            input_vocab_size=50,
            label_vocab_size=50,
            # embedding size
            max_length=100,
            hidden_units=512,
            dropout_rate=0.2,
            lr=0.0003,
            is_training=True)
        return params

# 模型各分层的实现
class BuildLayer(object):

    # layer norm层
    def normalize(self,
                  inputs,
                  epsilon = 1e-8,
                  scope="ln",
                  reuse=None):
        print('layer norm层正常----->>>6')
        with tf.variable_scope(scope, reuse=reuse):
            inputs_shape = inputs.get_shape()
            params_shape = inputs_shape[-1:]

            mean, variance = tf.nn.moments(inputs, [-1], keep_dims=True)
            beta= tf.Variable(tf.zeros(params_shape))
            gamma = tf.Variable(tf.ones(params_shape))
            normalized = (inputs - mean) / ( (variance + epsilon) ** (.5) )
            outputs = gamma * normalized + beta

        return outputs

    # embedding层
    def embedding(self,
                  inputs,
                  vocab_size,
                  num_units,
                  zero_pad=True,
                  scale=True,
                  scope="embedding",
                  reuse=None):
        print('embedding层正常----->>>7')
        with tf.variable_scope(scope, reuse=reuse):
            lookup_table = tf.get_variable('lookup_table',
                                           dtype=tf.float32,
                                           shape=[vocab_size, num_units],
                                           initializer=tf.contrib.layers.xavier_initializer())
            if zero_pad:
                lookup_table = tf.concat((tf.zeros(shape=[1, num_units]),
                                          lookup_table[1:, :]), 0)
            outputs = tf.nn.embedding_lookup(lookup_table, inputs)

            if scale:
                outputs = outputs * (num_units ** 0.5)

        return outputs

    # multihead层
    def multihead_attention(self,
                            emb,
                            queries,
                            keys,
                            num_units=None,
                            num_heads=8,
                            dropout_rate=0,
                            is_training=True,
                            causality=False,
                            scope="multihead_attention",
                            reuse=None):
        print('multihead层正常----->>>8')
        with tf.variable_scope(scope, reuse=reuse):
            # Set the fall back option for num_units
            if num_units is None:
                num_units = queries.get_shape().as_list[-1]

            # Linear projections
            Q = tf.layers.dense(queries, num_units, activation=tf.nn.relu)  # (N, T_q, C)
            K = tf.layers.dense(keys, num_units, activation=tf.nn.relu)  # (N, T_k, C)
            V = tf.layers.dense(keys, num_units, activation=tf.nn.relu)  # (N, T_k, C)

            # Split and concat
            Q_ = tf.concat(tf.split(Q, num_heads, axis=2), axis=0)  # (h*N, T_q, C/h)
            K_ = tf.concat(tf.split(K, num_heads, axis=2), axis=0)  # (h*N, T_k, C/h)
            V_ = tf.concat(tf.split(V, num_heads, axis=2), axis=0)  # (h*N, T_k, C/h)

            # Multiplication
            outputs = tf.matmul(Q_, tf.transpose(K_, [0, 2, 1]))  # (h*N, T_q, T_k)

            # Scale
            outputs = outputs / (K_.get_shape().as_list()[-1] ** 0.5)

            # Key Masking
            key_masks = tf.sign(tf.abs(tf.reduce_sum(emb, axis=-1)))  # (N, T_k)
            key_masks = tf.tile(key_masks, [num_heads, 1])  # (h*N, T_k)
            key_masks = tf.tile(tf.expand_dims(key_masks, 1), [1, tf.shape(queries)[1], 1])  # (h*N, T_q, T_k)

            paddings = tf.ones_like(outputs) * (-2 ** 32 + 1)
            outputs = tf.where(tf.equal(key_masks, 0), paddings, outputs)  # (h*N, T_q, T_k)

            # Causality = Future blinding
            if causality:
                diag_vals = tf.ones_like(outputs[0, :, :])  # (T_q, T_k)
                tril = tf.contrib.linalg.LinearOperatorTriL(diag_vals).to_dense()  # (T_q, T_k)
                masks = tf.tile(tf.expand_dims(tril, 0), [tf.shape(outputs)[0], 1, 1])  # (h*N, T_q, T_k)

                paddings = tf.ones_like(masks) * (-2 ** 32 + 1)
                outputs = tf.where(tf.equal(masks, 0), paddings, outputs)  # (h*N, T_q, T_k)

            # Activation
            outputs = tf.nn.softmax(outputs)  # (h*N, T_q, T_k)

            # Query Masking
            query_masks = tf.sign(tf.abs(tf.reduce_sum(emb, axis=-1)))  # (N, T_q)
            query_masks = tf.tile(query_masks, [num_heads, 1])  # (h*N, T_q)
            query_masks = tf.tile(tf.expand_dims(query_masks, -1), [1, 1, tf.shape(keys)[1]])  # (h*N, T_q, T_k)
            outputs *= query_masks  # broadcasting. (N, T_q, C)

            # Dropouts
            outputs = tf.layers.dropout(outputs, rate=dropout_rate, training=tf.convert_to_tensor(is_training))

            # Weighted sum
            outputs = tf.matmul(outputs, V_)  # ( h*N, T_q, C/h)

            # Restore shape
            outputs = tf.concat(tf.split(outputs, num_heads, axis=0), axis=2)  # (N, T_q, C)

            # Residual connection
            outputs += queries

            # Normalize
            outputs = self.normalize(outputs)  # (N, T_q, C)

        return outputs

    # feedforward
    def feedforward(self,
                    inputs,
                    num_units=[2048, 512],
                    scope="multihead_attention",
                    reuse=None):
        print('feedforward层正常----->>>9')
        with tf.variable_scope(scope, reuse=reuse):
            # Inner layer
            params = {"inputs": inputs, "filters": num_units[0], "kernel_size": 1,
                      "activation": tf.nn.relu, "use_bias": True}
            outputs = tf.layers.conv1d(**params)
            # Readout layer
            params = {"inputs": outputs, "filters": num_units[1], "kernel_size": 1,
                      "activation": None, "use_bias": True}
            outputs = tf.layers.conv1d(**params)
            # Residual connection
            outputs += inputs
            # Normalize
            outputs = self.normalize(outputs)

        return outputs

    # label_smoothing
    def label_smoothing(self, inputs, epsilon=0.1):
        print('label_smoothing正常----->>>9')
        K = inputs.get_shape().as_list()[-1]  # number of channels
        return ((1 - epsilon) * inputs) + (epsilon / K)

dealdat = DealData()
arg = dealdat.create_hparams()
pny2id, han2id = dealdat.pyn_han()
input_num, label_num = dealdat.input_lable_num()
arg.input_vocab_size = len(pny2id)
arg.label_vocab_size = len(han2id)

# 搭建模型
class Graph(object):
    def __init__(self, is_training=True):

        print('Graph初始化正常----->>>10')
        # self.dealdat = DealData()
        # self.arg = self.dealdat.create_hparams()
        # self.__pny2id, self.__han2id = self.dealdat.pyn_han()
        # self.__input_num, self.__label_num = self.dealdat.input_lable_num()
        # self.arg.input_vocab_size = len(self.__pny2id)
        # self.arg.label_vocab_size = len(self.__han2id)

        tf.reset_default_graph()
        self.is_training = arg.is_training
        self.hidden_units = arg.hidden_units
        self.input_vocab_size = arg.input_vocab_size
        self.label_vocab_size = arg.label_vocab_size
        self.num_heads = arg.num_heads
        self.num_blocks = arg.num_blocks
        self.max_length = arg.max_length
        self.lr = arg.lr
        self.dropout_rate = arg.dropout_rate

        self.laybuild = BuildLayer()
        # input
        self.x = tf.placeholder(tf.int32, shape=(None, None))
        self.y = tf.placeholder(tf.int32, shape=(None, None))
        # embedding
        self.emb = self.laybuild.embedding(self.x,
                                           vocab_size=self.input_vocab_size,
                                           num_units=self.hidden_units,
                                           scale=True,
                                           scope="enc_embed")
        self.enc = self.emb + self.laybuild.embedding(tf.tile(tf.expand_dims(tf.range(tf.shape(self.x)[1]), 0),
                                                              [tf.shape(self.x)[0], 1]),
                                                      vocab_size=self.max_length,
                                                      num_units=self.hidden_units,
                                                      zero_pad=False,
                                                      scale=False,
                                                      scope="enc_pe")
        # Dropout
        self.enc = tf.layers.dropout(self.enc,
                                     rate=self.dropout_rate,
                                     training=tf.convert_to_tensor(self.is_training))

        # Blocks
        for i in range(self.num_blocks):
            with tf.variable_scope("num_blocks_{}".format(i)):
                # Multihead Attention
                self.enc = self.laybuild.multihead_attention(emb=self.emb,
                                                             queries=self.enc,
                                                             keys=self.enc,
                                                             num_units=self.hidden_units,
                                                             num_heads=self.num_heads,
                                                             dropout_rate=self.dropout_rate,
                                                             is_training=self.is_training,
                                                             causality=False)

        # Feed Forward
        self.outputs = self.laybuild.feedforward(self.enc, num_units=[4 * self.hidden_units, self.hidden_units])

        # Final linear projection
        self.logits = tf.layers.dense(self.outputs, self.label_vocab_size)
        self.preds = tf.to_int32(tf.argmax(self.logits, axis=-1))
        self.istarget = tf.to_float(tf.not_equal(self.y, 0))
        self.acc = tf.reduce_sum(tf.to_float(tf.equal(self.preds, self.y)) * self.istarget) / (
            tf.reduce_sum(self.istarget))
        tf.summary.scalar('acc', self.acc)

        if is_training:
            # Loss
            self.y_smoothed = self.laybuild.label_smoothing(tf.one_hot(self.y, depth=self.label_vocab_size))
            self.loss = tf.nn.softmax_cross_entropy_with_logits(logits=self.logits, labels=self.y_smoothed)
            self.mean_loss = tf.reduce_sum(self.loss * self.istarget) / (tf.reduce_sum(self.istarget))

            # Training Scheme
            self.global_step = tf.Variable(0, name='global_step', trainable=False)
            self.optimizer = tf.train.AdamOptimizer(learning_rate=self.lr, beta1=0.9, beta2=0.98, epsilon=1e-8)
            self.train_op = self.optimizer.minimize(self.mean_loss, global_step=self.global_step)

            # Summary
            tf.summary.scalar('mean_loss', self.mean_loss)
            self.merged = tf.summary.merge_all()

# 模型训练
def train_model():

    print('模型训练正常----->>>11')
    g = Graph(arg)
    saver = tf.train.Saver()
    with tf.Session() as sess:
        merged = tf.summary.merge_all()
        sess.run(tf.global_variables_initializer())
        if os.path.exists('logs/model.meta'):
            saver.restore(sess, 'logs/model')
        writer = tf.summary.FileWriter('tensorboard/lm', tf.get_default_graph())
        for k in range(epochs):
            total_loss = 0
            batch_num = len(input_num) // batch_size
            batch = dealdat.get_batch(input_num, label_num, batch_size)
            for i in range(batch_num):
                input_batch, label_batch = next(batch)
                feed = {g.x: input_batch, g.y: label_batch}
                cost, _ = sess.run([g.mean_loss, g.train_op], feed_dict=feed)
                total_loss += cost
                if (k * batch_num + i) % 10 == 0:
                    rs=sess.run(merged, feed_dict=feed)
                    writer.add_summary(rs, k * batch_num + i)
            if (k+1) % 5 == 0:
                print('epochs', k+1, ': average loss = ', total_loss/batch_num)
        saver.save(sess, 'logs/model')
        writer.close()

# 模型推断
def model_predict():
    print('模型推断正常----->>>12')
    arg.is_training = False
    g = Graph(arg)
    saver = tf.train.Saver()

    with tf.Session() as sess:
        saver.restore(sess, 'logs/model')
        while True:
            line = input('输入测试拼音: ')
            if line == 'exit': break
            line = line.strip('\n').split(' ')
            print('--->', type(line))
            x = np.array([pny2id.index(pny) for pny in line])
            x = x.reshape(1, -1)
            preds = sess.run(g.preds, {g.x: x})
            got = ''.join(han2id[idx] for idx in preds[0])
            print(got)

if __name__ == '__main__':

    # train_model()
    model_predict()

ext: py

max_stars_repo_name: AustinKuture/MachineLearning

max_issues_repo_name: AustinKuture/MachineLearning

max_forks_repo_name: AustinKuture/MachineLearning

reason: import numpy


text:
from unittest import TestCase

import numpy as np
from dicom_parser.image import Image
from dicom_parser.utils.siemens.mosaic import Mosaic
from tests.fixtures import (
    TEST_RSFMRI_IMAGE_PATH,
    TEST_RSFMRI_IMAGE_VOLUME,
    TEST_RSFMRI_SERIES_PIXEL_ARRAY,
)

class MosaicTestCase(TestCase):
    @classmethod
    def setUpClass(cls):
        cls.image = Image(TEST_RSFMRI_IMAGE_PATH)
        cls.mosaic = Mosaic(cls.image._data, cls.image.header)

    def test_init_read_series_header_info(self):
        csa_header = self.mosaic.series_header_info
        self.assertIsInstance(csa_header, dict)

    def test_init_read_volume_shape(self):
        value = self.mosaic.volume_shape
        expected = 96, 96, 64
        self.assertTupleEqual(value, expected)

    def test_init_read_mosaic_dimensions(self):
        value = self.mosaic.mosaic_dimensions
        expected = 8, 8
        self.assertTupleEqual(value, expected)

    def test_ascending_attribute_with_ascending(self):
        value = self.mosaic.ascending
        self.assertTrue(value)

    # TODO: Find descending mosaic to test with
    # def test_ascending_attribute_with_descending(self):
    #     value = self.mosaic.ascending
    #     self.assertFalse(value)

    def test_folded_data(self):
        volume = self.mosaic.fold()
        expected = np.load(TEST_RSFMRI_IMAGE_VOLUME)
        self.assertTrue(np.array_equal(volume, expected))

    def test_folded_data_is_same_as_nifti(self):
        volume = self.mosaic.fold()
        nii_data = np.load(TEST_RSFMRI_SERIES_PIXEL_ARRAY)
        nii_volume = nii_data[:, :, :, 0]
        self.assertTrue(np.array_equal(volume, nii_volume))


text:
# -*- coding: utf-8 -*-
# Run this app with `python3 app.py` and
# make sure you activated the venc
# visit http://127.0.0.1:8050/ in your web browser.

import dash
import dash_core_components as dcc
import dash_html_components as html
import plotly.express as px
import pandas as pd
import numpy as np
from datetime import datetime as dt
import statsmodels.api as sm
import plotly.express as px

from urllib.request import urlopen
import json
import plotly

external_stylesheets = ['https://codepen.io/chriddyp/pen/bWLwgP.css']

app = dash.Dash(__name__, external_stylesheets=external_stylesheets)

# assume you have a "long-form" data frame
# see https://plotly.com/python/px-arguments/ for more options

# DATA manipulation
df_clean = pd.read_csv("data/clean_data.csv")
df_clean = df_clean[['date', 'Name', 'state', 'FIPS', 'positive', 'AQI', 'date_str']]
df_clean['date'] = pd.to_datetime(df_clean['date'], format='%Y/%m/%d')
df_clean.sort_values(by=["date", "state"])
df_clean["rank"] =  df_clean[['positive', 'AQI']].apply(lambda x: x[0]/(9.041993e+06) + x[1]/(153.285714) , axis=1)

top_10_states = df_clean.loc[df_clean["date"] == dt.strptime('2020-09-21', '%Y-%m-%d')].sort_values(by = 'positive').tail(10)['state']
df_top10_select = df_clean.loc[(df_clean["state"].isin(top_10_states)) & (df_clean["date"] >= dt.strptime('2020-03-15', '%Y-%m-%d')),].copy()

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

############### US RANKS
rank = px.choropleth(df_clean,
                           color='rank',
                           locationmode = 'USA-states',
                           locations='state',
                           color_continuous_scale="thermal_r",
                           hover_name = 'Name',
                           hover_data = ['AQI', 'positive'],
                           scope = 'usa',template = 'plotly_dark',
                           animation_frame = 'date_str',
                           animation_group = 'date_str',
                           width = 1600,
                           height = 800,)
rank.update_layout(margin={"r":0,"t":50,"l":0,"b":0})

############## US MAP positvive
positive = px.choropleth(df_clean,
                           color='positive',
                           locationmode = 'USA-states',
                           locations='state',
                           color_continuous_scale="thermal_r",
                           hover_name = 'Name',
                           scope = 'usa',
                           animation_frame = 'date_str',
                           animation_group = 'date_str',
                           template = 'plotly_dark',
                           labels={'tests_per_1000':'Tests Per 1000 People'},
                           #title = "US State-By-State Corona Cases"
                            )
positive.update_layout(margin={"r":0,"t":50,"l":0,"b":0})

aqi = px.choropleth(df_clean,
                           color='AQI',
                           locationmode = 'USA-states',
                           locations='state',
                           color_continuous_scale="thermal_r",
                           hover_name = 'Name',
                           hover_data = ['positive'],
                           scope = 'usa',
                           animation_frame = 'date_str',
                           animation_group = 'date_str', template = 'plotly_dark', )
aqi.update_layout(margin={"r":0,"t":45,"l":0,"b":0})

bubbles = px.scatter(df_clean, x="date_str", y="AQI",
	               size="positive", color="state",
                 hover_name="Name", size_max=60,
                 template = 'plotly_dark',
                 title = "Correlation of AQI vs Positives in Bubbles")

ols = px.scatter(df_top10_select,
                 x="rank", y="AQI",
                 color = "state", trendline = "ols",template = 'plotly_dark',
                 hover_name = "Name", title = "OLS timeless RANK vs AQI")

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

app.layout = html.Div(
    id="root",
    children=[
        html.Div(
            id="header",
            children=[
                html.H4(id = "NAME1", children="Data Analysis of Corona vs AQI"),
                html.P(
                    id="description",
                    children="Projects in Data Science Final Project by Andriy, Daniella, and Nurgazy.",
                ),
            ],
        ), ###### RANK MAP ############3333
        html.Div(
            className="app-container",
            children=[
                html.Div(

                    children=[
                        html.Div(
                            className="heatmap-container",
                            children=[
                                html.P(
                                    "US States Ranks over Time",
                                    className ="heatmap-title",
                                ),
                                dcc.Graph(

                                    className='county-choropleth',
                                    figure=rank
                                )
                            ],
                        )
                    ]
                )
            ]
        ),
        html.Div(
            className="float-container",
            children=[
                html.Div(
                    children=[
                        html.Div(
                            className="heatmap-container",
                            children=[
                                html.P(
                                    "US States AQI over time",
                                    className ="heatmap-title",
                                ),
                                dcc.Graph(

                                    className='county-choropleth',
                                    figure=aqi
                                )
                            ],
                        ),
                        html.Div(
                            className="heatmap-container",

                            children=[
                                html.P(
                                    "US States Positive Corona Cases over time",
                                    className ="heatmap-title",
                                ),
                                dcc.Graph(

                                    className='county-choropleth',
                                    figure=positive
                                )
                            ],
                        ),
                        html.Div(
                            className="heatmap-container",

                            children=[
                                html.P(
                                    "Correlation of AQI and Positive Tests",
                                    className ="heatmap-title",
                                ),
                                dcc.Graph(
                                    className='county-choropleth',
                                    figure=bubbles
                                )
                            ],
                        ),
                        html.Div(
                            className="heatmap-container",

                            children=[
                                html.P(
                                    "Regression fit on the data",
                                    className ="heatmap-title",
                                ),
                                dcc.Graph(
                                    className='county-choropleth',
                                    figure=ols
                                )
                            ],
                        )
                    ]
                )
            ]
        )
    ]
)

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

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

reason: import numpy,import statsmodels


text:
import numpy
from collections import Counter

import png

def otsu(histogram, total_pixels):
    mu_total = sum(g * histogram.get(g) for g in range(256))
    max_t = 0
    max_value = 0
    p1 = 0.0
    mu1 = 0.0
    for t in range(256 - 1):
        p1 += histogram.get(t)
        p2 = total_pixels - p1
        if p1 == 0 or p2 == 0:
            continue
        mu1 += histogram.get(t) * t
        mu2 = mu_total - mu1

        tmp = p1 * p2 * (mu1 / p1 - mu2 / p2) ** 2
        if tmp > max_value:
            max_value = tmp
            max_t = t
    return max_t

def apply_threshold(png_file, threshold):
    with open(png_file, 'r') as f:
        bitmap = png.Reader(f).asDirect()
        new_bitmap = []
        for row in bitmap[2]:
            shape = len(row)
            new_row = numpy.empty(shape)
            for i in range(shape):
                new_row[i] = 0 if row[i] <= threshold else 255
            new_bitmap.append(new_row)
    return new_bitmap

def binarize(png_file, output_file):
    with open(png_file, 'r') as f:
        reader = png.Reader(f)
        data = reader.asDirect()
        bitmap = data[2]
        counter = Counter()
        for row in bitmap:
            counter.update(row)
        total_pixels = data[0] * data[1]
        threshold = otsu(counter, total_pixels)
    new_bitmap = apply_threshold(png_file, threshold)
    with open(output_file, 'wb') as f:
        w = png.Writer(data[0], data[1], greyscale=True, bitdepth=8, gamma=1.0, background=(1,))
        w.write(f, new_bitmap)

binarize("../resources/graysample.png", "../resources/graysample2.png")

max_stars_repo_name: hlibbabii/uni-dca

max_issues_repo_name: hlibbabii/uni-dca

max_forks_repo_name: hlibbabii/uni-dca


text:
################################################################################
# Peach - Computational Intelligence for Python
# Jose Alexandre Nalon
#
# This file: nn/nn.py
# Basic topologies of neural networks
################################################################################

# Doc string, reStructuredText formatted:
__doc__ = """
Basic topologies of neural networks.

This sub-package implements various neural network topologies, see the complete
list below. These topologies are implemented using the ``Layer`` class of the
``base`` sub-package. Please, consult the documentation of that module for more
information on layers of neurons. The neural nets implemented here don't derive
from the ``Layer`` class, instead, they have instance variables to take control
of them. Thus, there is no base class for networks. While subclassing the
classes of this module is usually safe, it is recomended that a new kind of
net is developed from the ground up.
"""

################################################################################
from numpy import sum, abs, reshape, sqrt, argmin
import random

from base import *
from af import *
from lrules import *

################################################################################
class FeedForward(list):
    '''
    Classic completely connected neural network.

    A feedforward neural network is implemented as a list of layers, each layer
    being a ``Layer`` object (please consult the documentation on the ``base``
    module for more information on layers). The layers are completely connected,
    which means that every neuron in one layers is connected to every other
    neuron in the following layer.

    There is a number of learning methods that are already implemented, but in
    general, any learning class derived from ``FFLearning`` can be used. No
    other kind of learning can be used. Please, consult the documentation on the
    ``lrules`` (*learning rules*) module.
    '''
    def __init__(self, layers, phi=Linear, lrule=BackPropagation, bias=False):
        '''
        Initializes a feedforward neural network.

        A feedforward network is implemented as a list of layers, completely
        connected.

        :Parameters:
          layers
            A list of integers containing the shape of the network. The first
            element of the list is the number of inputs of the network (or, as
            somebody prefer, the number of input neurons); the number of outputs
            is the number of neurons in the last layer. Thus, at least two
            numbers should be given.
          phi
            The activation functions to be used with each layer of the network.
            Please consult the ``Layer`` documentation in the ``base`` module
            for more information. This parameter can be a single function or a
            list of functions. If only one function is given, then the same
            function is used in every layer. If a list of functions is given,
            then the layers use the functions in the sequence given. Note that
            heterogeneous networks can be created that way. Defaults to
            ``Linear``.
          lrule
            The learning rule used. Only ``FFLearning`` objects (instances of
            the class or of the subclasses) are allowed. Defaults to
            ``BackPropagation``. Check the ``lrules`` documentation for more
            information.
          bias
            If ``True``, then the neurons are biased.
        '''
        list.__init__(self, [ ])
        layers = list(layers)
        for n, m in zip(layers[:-1], layers[1:]):
            self.append(Layer((m, n), bias=bias))
        self.phi = phi
        self.__n = len(self)
        self.__lrule = lrule
        if isinstance(lrule, FFLearning):
            self.__lrule = lrule
        else:
            try:
                issubclass(lrule, FFLearning)
                self.__lrule = lrule()
            except TypeError:
                raise ValueError, 'uncompatible learning rule'

    def __getnlayers(self):
        return self.__n
    nlayers = property(__getnlayers, None)
    '''Number of layers of the neural network. Not writable.'''

    def __getbias(self):
        r = [ ]
        for l in self:
            r.append(l.bias)
        return tuple(r)
    bias = property(__getbias, None)
    '''A tuple containing the bias of each layer. Not writable.'''

    def __gety(self):
        return self[-1].y
    y = property(__gety, None)
    '''A list of activation values for each neuron in the last layer of the
    network, ie., the answer of the network. This property is available only
    after the network is fed some input.'''

    def __getphi(self):
        r = [ ]
        for l in self:
            r.append(l.phi)
        return tuple(r)
    def __setphi(self, phis):
        try:
            phis = tuple(phis)
            for w, v in zip(self, phis):
                w.phi = v
        except TypeError:
            for w in self:
                w.phi = phis
    phi = property(__getphi, __setphi)
    '''Activation functions for every layer in the network. It is a list of
    ``Activation`` objects, but can be set with only one function. In this case,
    the same function is used for every layer.'''

    def __call__(self, x):
        '''
        The feedforward method of the network.

        The ``__call__`` interface should be called if the answer of the neuron
        network to a given input vector ``x`` is desired. *This method has
        collateral effects*, so beware. After the calling of this method, the
        ``y`` property is set with the activation potential and the answer of
        the neurons, respectivelly.

        :Parameters:
          x
            The input vector to the network.

        :Returns:
          The vector containing the answer of every neuron in the last layer, in
          the respective order.
        '''
        for w in self:
            x = w(x)
        return self[-1].y

    def learn(self, x, d):
        '''
        Applies one example of the training set to the network.

        Using this method, one iteration of the learning procedure is made with
        the neurons of this network. This method presents one example (not
        necessarilly of a training set) and applies the learning rule over the
        network. The learning rule is defined in the initialization of the
        network, and some are implemented on the ``lrules`` method. New methods
        can be created, consult the ``lrules`` documentation but, for
        ``FeedForward`` instances, only ``FFLearning`` learning is allowed.

        Also, notice that *this method only applies the learning method!* The
        network should be fed with the same input vector before trying to learn
        anything first. Consult the ``feed`` and ``train`` methods below for
        more ways to train a network.

        :Parameters:
          x
            Input vector of the example. It should be a column vector of the
            correct dimension, that is, the number of input neurons.
          d
            The desired answer of the network for this particular input vector.
            Notice that the desired answer should have the same dimension of the
            last layer of the network. This means that a desired answer should
            be given for every output of the network.

        :Returns:
          The error obtained by the network.
        '''
        self.__lrule(self, x, d)
        return sum(abs(d - self.y))

    def feed(self, x, d):
        '''
        Feed the network and applies one example of the training set to the
        network.

        Using this method, one iteration of the learning procedure is made with
        the neurons of this network. This method presents one example (not
        necessarilly of a training set) and applies the learning rule over the
        network. The learning rule is defined in the initialization of the
        network, and some are implemented on the ``lrules`` method. New methods
        can be created, consult the ``lrules`` documentation but, for
        ``FeedForward`` instances, only ``FFLearning`` learning is allowed.

        Also, notice that *this method feeds the network* before applying the
        learning rule. Feeding the network has collateral effects, and some
        properties change when this happens. Namely, the ``y`` property is set.
        Please consult the ``__call__`` interface.

        :Parameters:
          x
            Input vector of the example. It should be a column vector of the
            correct dimension, that is, the number of input neurons.
          d
            The desired answer of the network for this particular input vector.
            Notice that the desired answer should have the same dimension of the
            last layer of the network. This means that a desired answer should
            be given for every output of the network.

        :Returns:
          The error obtained by the network.
        '''
        self(x)
        return self.learn(x, d)

    def train(self, train_set, imax=2000, emax=1e-5, randomize=False):
        '''
        Presents a training set to the network.

        This method automatizes the training of the network. Given a training
        set, the examples are shown to the network (possibly in a randomized
        way). A maximum number of iterations or a maximum admitted error should
        be given as a stop condition.

        :Parameters:
          train_set
            The training set is a list of examples. It can have any size and can
            contain repeated examples. In fact, the definition of the training
            set is open. Each element of the training set, however, should be a
            two-tuple ``(x, d)``, where ``x`` is the input vector, and ``d`` is
            the desired response of the network for this particular input. See
            the ``learn`` and ``feed`` for more information.
          imax
            The maximum number of iterations. Examples from the training set
            will be presented to the network while this limit is not reached.
            Defaults to 2000.
          emax
            The maximum admitted error. Examples from the training set will be
            presented to the network until the error obtained is lower than this
            limit. Defaults to 1e-5.
          randomize
            If this is ``True``, then the examples are shown in a randomized
            order. If ``False``, then the examples are shown in the same order
            that they appear in the ``train_set`` list. Defaults to ``False``.
        '''
        i = 0
        error = 1
        s = len(train_set)
        while i<imax and error>emax:
            if randomize:
                x, d = random.choice(train_set)
            else:
                x, d = train_set[i%s]
            error = self.feed(x, d)
            i = i+1
        return error

################################################################################
class SOM(Layer):
    '''
    A Self-Organizing Map (SOM).

    A self-organizing map is a type of neural network that is trained via
    unsupervised learning. In particular, the self-organizing map finds the
    neuron closest to an input vector -- this neuron is the winning neuron, and
    it is the answer of the network. Thus, the SOM is usually used for
    classification and pattern recognition.

    The SOM is a single-layer network, so this class subclasses the ``Layer``
    class. But some of the properties of a ``Layer`` object are not available or
    make no sense in this context.
    '''
    def __init__(self, shape, lrule=Competitive):
        '''
        Initializes a self-organizing map.

        A self-organizing map is implemented as a layer of neurons. There is no
        connection among the neurons. The answer to a given input is the neuron
        closer to the given input. ``phi`` (the activation function) ``v`` (the
        activation potential) and ``bias`` are not used.

        :Parameters:
          shape
            Stablishes the size of the SOM. It must be a two-tuple of the
            format ``(m, n)``, where ``m`` is the number of neurons in the
            layer, and ``n`` is the number of inputs of each neuron. The neurons
            in the layer all have the same number of inputs.
          lrule
            The learning rule used. Only ``SOMLearning`` objects (instances of
            the class or of the subclasses) are allowed. Defaults to
            ``Competitive``. Check the ``lrules`` documentation for more
            information.
        '''
        Layer.__init__(self, shape, phi=None, bias=False)
        self.__lrule = lrule
        self.__y = None
        self.__phi = None
        if isinstance(lrule, SOMLearning):
            self.__lrule = lrule
        else:
            try:
                issubclass(lrule, SOMLearning)
                self.__lrule = lrule()
            except TypeError:
                raise ValueError, 'uncompatible learning rule'

    def __gety(self):
        if self.__y is None:
            raise ValueError, "activation unavailable"
        else:
            return self.__y
    y = property(__gety, None)
    '''The winning neuron for a given input, the answer of the network. This
    property is available only after the network is fed some input.'''

    def __call__(self, x):
        '''
        The response of the network to a given input.

        The ``__call__`` interface should be called if the answer of the neuron
        network to a given input vector ``x`` is desired. *This method has
        collateral effects*, so beware. After the calling of this method, the
        ``y`` property is set with the activation potential and the answer of
        the neurons, respectivelly.

        :Parameters:
          x
            The input vector to the network.

        :Returns:
          The winning neuron.
        '''
        x = reshape(x, (1, self.inputs))
        dist = sqrt(sum((x - self.weights)**2, axis=1))
        self.__y = argmin(dist)
        return self.y

    def learn(self, x):
        '''
        Applies one example of the training set to the network.

        Using this method, one iteration of the learning procedure is made with
        the neurons of this network. This method presents one example (not
        necessarilly of a training set) and applies the learning rule over the
        network. The learning rule is defined in the initialization of the
        network, and some are implemented on the ``lrules`` method. New methods
        can be created, consult the ``lrules`` documentation but, for
        ``SOM`` instances, only ``SOMLearning`` learning is allowed.

        Also, notice that *this method only applies the learning method!* The
        network should be fed with the same input vector before trying to learn
        anything first. Consult the ``feed`` and ``train`` methods below for
        more ways to train a network.

        :Parameters:
          x
            Input vector of the example. It should be a column vector of the
            correct dimension, that is, the number of input neurons.

        :Returns:
          The error obtained by the network.
        '''
        self.__lrule(self, x)
        return sum(abs(x - self.y))

    def feed(self, x):
        '''
        Feed the network and applies one example of the training set to the
        network.

        Using this method, one iteration of the learning procedure is made with
        the neurons of this network. This method presents one example (not
        necessarilly of a training set) and applies the learning rule over the
        network. The learning rule is defined in the initialization of the
        network, and some are implemented on the ``lrules`` method. New methods
        can be created, consult the ``lrules`` documentation but, for
        ``SOM`` instances, only ``SOMLearning`` learning is allowed.

        Also, notice that *this method feeds the network* before applying the
        learning rule. Feeding the network has collateral effects, and some
        properties change when this happens. Namely, the ``y`` property is set.
        Please consult the ``__call__`` interface.

        :Parameters:
          x
            Input vector of the example. It should be a column vector of the
            correct dimension, that is, the number of input neurons.

        :Returns:
          The error obtained by the network.
        '''
        self(x)
        return self.learn(x)

    def train(self, train_set, imax=2000, emax=1e-5, randomize=False):
        '''
        Presents a training set to the network.

        This method automatizes the training of the network. Given a training
        set, the examples are shown to the network (possibly in a randomized
        way). A maximum number of iterations or a maximum admitted error should
        be given as a stop condition.

        :Parameters:
          train_set
            The training set is a list of examples. It can have any size and can
            contain repeated examples. In fact, the definition of the training
            set is open. Each element of the training set, however, should be a
            input vector of the correct dimensions, See the ``learn`` and
            ``feed`` for more information.
          imax
            The maximum number of iterations. Examples from the training set
            will be presented to the network while this limit is not reached.
            Defaults to 2000.
          emax
            The maximum admitted error. Examples from the training set will be
            presented to the network until the error obtained is lower than this
            limit. Defaults to 1e-5.
          randomize
            If this is ``True``, then the examples are shown in a randomized
            order. If ``False``, then the examples are shown in the same order
            that they appear in the ``train_set`` list. Defaults to ``False``.
        '''
        i = 0
        error = 1
        s = len(train_set)
        while i<imax and error>emax:
            if randomize:
                x = random.choice(train_set)
            else:
                x = train_set[i%s]
            error = self.feed(x)
            i = i+1
        return error

################################################################################
# Test
if __name__ == "__main__":
    pass

reason: from numpy


text:
import numpy; from numpy import array
from algopy import UTPM, zeros

def F(x):
    y = zeros(3, dtype=x)
    y[0] = x[0]*x[1]
    y[1] = x[1]*x[2]
    y[2] = x[2]*x[0]
    return y

x0 = array([1,3,5],dtype=float)
x1 = array([1,0,0],dtype=float)

D = 2; P = 1
x = UTPM(numpy.zeros((D,P,3)))
x.data[0,0,:] = x0
x.data[1,0,:] = x1

# normal function evaluation
y0 = F(x0)

# UTP function evaluation
y = F(x)

print('y0 = ', y0)
print('y  = ', y)
print('y.shape =', y.shape)
print('y.data.shape =', y.data.shape)
print('dF/dx(x0) * x1 =', y.data[1,0])

import numpy; from numpy import log, exp, sin, cos, abs
import algopy; from algopy import UTPM, dot, inv, zeros

def f(x):
    A = zeros((2,2),dtype=x)
    A[0,0] = numpy.log(x[0]*x[1])
    A[0,1] = numpy.log(x[1]) + exp(x[0])
    A[1,0] = sin(x[0])**2 + abs(cos(x[0]))**3.1
    A[1,1] = x[0]**cos(x[1])
    return log( dot(x.T,  dot( inv(A), x)))

x = numpy.array([3.,7.])
x = UTPM.init_jacobian(x)

y = f(x)

print('normal function evaluation f(x) = ',y.data[0,0])
print('Jacobian df/dx = ', UTPM.extract_jacobian(y))


text:
import os
import numpy as np

datafolder="/home/jliu447/lossycompression/NYX/"

ebs=[i*1e-4 for i in range(1,10)]+[i*1e-3 for i in range(1,10)]+[i*1e-2 for i in range(1,11)]
fields=["baryon_density","temperature","dark_matter_density"]
cr=np.zeros((len(ebs)+1,len(fields)+1),dtype=np.float32)
psnr=np.zeros((len(ebs)+1,len(fields)+1),dtype=np.float32)
for i,eb in enumerate(ebs):
    cr[i+1][0]=eb
    psnr[i+1][0]=eb
    for j,field in enumerate(fields):

        filename="%s_3.dat.log10" % str(field)
        filepath=os.path.join(datafolder,filename)
        comm="sz -z -f -i %s -M REL -R %f -3 512 512 512" % (filepath,eb)
        os.system(comm)
        szpath=filepath+".sz"
        comm="sz -x -f -i %s -s %s -3 512 512 512 -a>temp.txt" % (filepath,szpath)
        os.system(comm)
        with open("temp.txt","r") as f:
            lines=f.read().splitlines()
            #print(lines)

            p=eval(lines[4].split(',')[0].split('=')[1])

            r=eval(lines[7].split('=')[1])

            cr[i+1][j+1]=r
            psnr[i+1][j+1]=p

        comm="rm -f %s" % szpath
        os.system(comm)
        comm="rm -f %s" % szpath+".out"
        os.system(comm)
        os.system("rm -f temp.txt")

np.savetxt("sz_nyx3log10_cr.txt",cr,delimiter='\t')
np.savetxt("sz_nyx3log10_psnr.txt",psnr,delimiter='\t')


text:
import numpy as np
import tensorflow as tf

from collections import namedtuple
from util_conv import *

def map_str_to_func(str_activ):
    if str_activ == 'elu':
        activ = tf.nn.elu
    elif str_activ == 'relu':
        activ = tf.nn.relu
    elif str_activ == 'tanh':
        activ = tf.tanh
    elif str_activ == 'sigmoid':
        activ = tf.sigmoid
    else:
        raise ValueError("Unsupported activation function ({})".format(str_activ))
    return activ

def get_activ_func(activ):
    if activ is None:
        return None
    elif type(activ) is str:
        return map_str_to_func(activ)
    elif type(activ) is list:
        return [map_str_to_func(x) for x in activ]
    else:
        raise ValueError("Definition of activation units should be string or list of string")

class decoder_interface(object):
    def get_shape(self):
        raise NotImplementedError # expected to return (height, width, layers)

class decoder_pic(decoder_interface):
    HParams = namedtuple('HParams',
                         'latent_size, output_shape, activ_unit, last_activ_unit, do_batch_norm')

    def __init__(self, hps):
        self.hps = hps
        self.activ = get_activ_func(self.hps.activ_unit)
        self.last_active = get_activ_func(self.hps.last_activ_unit)

    def get_shape(self):
        return self.hps.output_shape
    def get_depth(self):
        return convert_shape_to_depth(self.hps.output_shape[0:2])

    def build_decoder(self, z, is_training, concat_dict=None):
        latent_size = self.hps.latent_size
        output_shape = self.hps.output_shape
        activ = self.activ
        last_active = self.last_active
        do_BN = self.hps.do_batch_norm

        if len(output_shape) != 3:
            raise ValueError('output_shape should be (height, width, channel)')
        if type(last_active) is list:
            if len(last_active) != output_shape[2]:
                raise ValueError("# of elements of last_active should be same with that of output_shape")

        return self.decoder(z, latent_size, output_shape, activ=activ, last_active=last_active,
                            is_training=is_training, concat_dict=concat_dict, do_BN=do_BN)

    # concat_dict: python dictionary of (height, data), it will be concatenated (along the 3rd dim.) with a layer having same height
    # assume height == width
    def decoder_unit(self, layer_name, x, out_filters, concat_dict=None, is_training=True, activ=tf.nn.elu, do_BN=False):
        height, width = x.get_shape().as_list()[1:3] # height and width of input
        assert height == width
        size = height

        # concatenate the provided data if available
        if (concat_dict is not None) and size in concat_dict.keys():
            data_concat = concat_dict[size]
            x = tf.concat([x, data_concat], 3) # concatenate with the resized partial input
            do_concat = True
        # increase the resolution
        x = conv_transpose_unit('conv_t_' + layer_name, x, 3, out_filters, 2,
                                is_training=is_training, activ_func=activ, do_batch_norm=do_BN)
        #additional convolution
        x = conv_unit('conv_shuf_' + layer_name, x, 3, out_filters, 1,
                      is_training=is_training, activ_func=activ, do_batch_norm=do_BN) # convolution while keeping the resolution
        return x

    def decoder_unit_interp(self, layer_name, x, out_filters, concat_dict=None, is_training=True, activ=tf.nn.elu, do_BN=False):
        height, width = x.get_shape().as_list()[1:3] # height and width of input
        assert height == width
        size = height
        size_up = size * 2

        x = tf.image.resize_images(x, [size_up, size_up])
        # concatenate the provided data if available
        if (concat_dict is not None) and size_up in concat_dict.keys():
            print(x.shape,concat_dict[size_up].shape)
            data_concat = concat_dict[size_up]
            x = tf.concat([x, data_concat], 3)
        # increase the resolution
        x = conv_unit('conv_shuf_{}_1'.format(layer_name), x, 3, out_filters, 1,
                                is_training=is_training, activ_func=activ, do_batch_norm=do_BN)
        x = conv_unit('conv_shuf_{}_2'.format(layer_name), x, 3, out_filters, 1,
                      is_training=is_training, activ_func=activ, do_batch_norm=do_BN)
        return x
    def decoder(self, z, latent_size, output_shape, activ=tf.nn.elu, last_active=None,  is_training=True, concat_dict=None, do_BN=False):
        depth = convert_shape_to_depth(output_shape[0:2])
        last_layers = output_shape[2]
        #filters = {2:128, 4:128, 8:64, 16:64, 32:64, 64:64, 128:32}
        #filters = {2:512, 4:256, 8:128, 16:64, 32:64, 64:64, 128:32}
        filters = {2:1024, 4:512, 8:256, 16:128, 32:64, 64:64, 128:32}
        #size = 4
        #x = fully_connected('linear_decoder', z, 2048)
        #x = tf.reshape(x, [-1, size, size, 128]) # convert the code to 1x1 image
        size = 1
        x = tf.reshape(z, [-1, size, size, latent_size]) # convert the code to 1x1 image

        # transposed convolutions
        for i in range(depth):
            size = 2<<(i) # target resolution(just for naming the operation)
            out_filters = filters[size]
            #if i < 3:
            #    x = self.decoder_unit(str(size), x, out_filters, concat_dict, is_training, activ, do_BN)
            #else:
            #    x = self.decoder_unit_interp(str(size), x, out_filters, concat_dict, is_training, activ, do_BN)

            x = self.decoder_unit(str(size), x, out_filters, concat_dict, is_training, activ, do_BN)
            #x = self.decoder_unit_interp(str(size), x, out_filters, concat_dict, is_training, activ, do_BN)

        # final deconvolution
        #x = self.decoder_unit_interp(str(1<<depth), x, filters[2],  concat_dict, is_training, activ, do_BN)

        # final 1x1 convolution for final channels
        if last_layers == 1:
            x = conv_unit('final_conv', x, 1, last_layers, 1,
                          is_training=is_training, activ_func=last_active, do_batch_norm=False)
        else:
            x = conv_unit('final_conv', x, 1, last_layers, 1,
                          is_training=is_training, activ_func=None, do_batch_norm=False)
            xs = list(tf.split(x, last_layers, axis=3))
            for i, _ in enumerate(xs):
                xs[i] = last_active[i](xs[i])
            x = tf.concat(xs,axis=3)

        return x


text:
import glob
import json
import logging
import numpy as np
import os
from pathlib import Path
import random
import re
import tree  # pip install dm_tree
from typing import List, Optional, TYPE_CHECKING, Union
from urllib.parse import urlparse
import zipfile

try:
    from smart_open import smart_open
except ImportError:
    smart_open = None

from ray.rllib.offline.input_reader import InputReader
from ray.rllib.offline.io_context import IOContext
from ray.rllib.policy.policy import Policy
from ray.rllib.policy.sample_batch import (
    DEFAULT_POLICY_ID,
    MultiAgentBatch,
    SampleBatch,
)
from ray.rllib.utils.annotations import override, PublicAPI
from ray.rllib.utils.compression import unpack_if_needed
from ray.rllib.utils.spaces.space_utils import clip_action, normalize_action
from ray.rllib.utils.typing import Any, FileType, SampleBatchType

if TYPE_CHECKING:
    from ray.rllib.evaluation import RolloutWorker

logger = logging.getLogger(__name__)

WINDOWS_DRIVES = [chr(i) for i in range(ord("c"), ord("z") + 1)]

def _adjust_obs_actions_for_policy(json_data: dict, policy: Policy) -> dict:
    """Handle nested action/observation spaces for policies.

    Translates nested lists/dicts from the json into proper
    np.ndarrays, according to the (nested) observation- and action-
    spaces of the given policy.

    Providing nested lists w/o this preprocessing step would
    confuse a SampleBatch constructor.
    """
    for k, v in policy.view_requirements.items():
        if k not in json_data:
            continue
        if policy.config.get("_disable_action_flattening") and (
            k == SampleBatch.ACTIONS or v.data_col == SampleBatch.ACTIONS
        ):
            json_data[k] = tree.map_structure_up_to(
                policy.action_space_struct,
                lambda comp: np.array(comp),
                json_data[k],
                check_types=False,
            )
        elif policy.config.get("_disable_preprocessor_api") and (
            k == SampleBatch.OBS or v.data_col == SampleBatch.OBS
        ):
            json_data[k] = tree.map_structure_up_to(
                policy.observation_space_struct,
                lambda comp: np.array(comp),
                json_data[k],
                check_types=False,
            )
    return json_data

def from_json_data(json_data: Any, worker: Optional["RolloutWorker"]):
    # Try to infer the SampleBatchType (SampleBatch or MultiAgentBatch).
    if "type" in json_data:
        data_type = json_data.pop("type")
    else:
        raise ValueError("JSON record missing 'type' field")

    if data_type == "SampleBatch":
        if worker is not None and len(worker.policy_map) != 1:
            raise ValueError(
                "Found single-agent SampleBatch in input file, but our "
                "PolicyMap contains more than 1 policy!"
            )
        for k, v in json_data.items():
            json_data[k] = unpack_if_needed(v)
        if worker is not None:
            policy = next(iter(worker.policy_map.values()))
            json_data = _adjust_obs_actions_for_policy(json_data, policy)
        return SampleBatch(json_data)
    elif data_type == "MultiAgentBatch":
        policy_batches = {}
        for policy_id, policy_batch in json_data["policy_batches"].items():
            inner = {}
            for k, v in policy_batch.items():
                inner[k] = unpack_if_needed(v)
            if worker is not None:
                policy = worker.policy_map[policy_id]
                inner = _adjust_obs_actions_for_policy(inner, policy)
            policy_batches[policy_id] = SampleBatch(inner)
        return MultiAgentBatch(policy_batches, json_data["count"])
    else:
        raise ValueError(
            "Type field must be one of ['SampleBatch', 'MultiAgentBatch']", data_type
        )

# TODO(jungong) : use DatasetReader to back JsonReader, so we reduce
#     codebase complexity without losing existing functionality.
@PublicAPI
class JsonReader(InputReader):
    """Reader object that loads experiences from JSON file chunks.

    The input files will be read from in random order.
    """

    @PublicAPI
    def __init__(
        self, inputs: Union[str, List[str]], ioctx: Optional[IOContext] = None
    ):
        """Initializes a JsonReader instance.

        Args:
            inputs: Either a glob expression for files, e.g. `/tmp/**/*.json`,
                or a list of single file paths or URIs, e.g.,
                ["s3://bucket/file.json", "s3://bucket/file2.json"].
            ioctx: Current IO context object or None.
        """
        logger.info(
            "You are using JSONReader. It is recommended to use "
            + "DatasetReader instead for better sharding support."
        )

        self.ioctx = ioctx or IOContext()
        self.default_policy = self.policy_map = None
        if self.ioctx.worker is not None:
            self.policy_map = self.ioctx.worker.policy_map
            self.default_policy = self.policy_map.get(DEFAULT_POLICY_ID)

        if isinstance(inputs, str):
            inputs = os.path.abspath(os.path.expanduser(inputs))
            if os.path.isdir(inputs):
                inputs = [os.path.join(inputs, "*.json"), os.path.join(inputs, "*.zip")]
                logger.warning(f"Treating input directory as glob patterns: {inputs}")
            else:
                inputs = [inputs]

            if any(urlparse(i).scheme not in [""] + WINDOWS_DRIVES for i in inputs):
                raise ValueError(
                    "Don't know how to glob over `{}`, ".format(inputs)
                    + "please specify a list of files to read instead."
                )
            else:
                self.files = []
                for i in inputs:
                    self.files.extend(glob.glob(i))
        elif isinstance(inputs, (list, tuple)):
            self.files = list(inputs)
        else:
            raise ValueError(
                "type of inputs must be list or str, not {}".format(inputs)
            )
        if self.files:
            logger.info("Found {} input files.".format(len(self.files)))
        else:
            raise ValueError("No files found matching {}".format(inputs))
        self.cur_file = None

    @override(InputReader)
    def next(self) -> SampleBatchType:
        batch = self._try_parse(self._next_line())
        tries = 0
        while not batch and tries < 100:
            tries += 1
            logger.debug("Skipping empty line in {}".format(self.cur_file))
            batch = self._try_parse(self._next_line())
        if not batch:
            raise ValueError(
                "Failed to read valid experience batch from file: {}".format(
                    self.cur_file
                )
            )

        return self._postprocess_if_needed(batch)

    def read_all_files(self) -> SampleBatchType:
        """Reads through all files and yields one SampleBatchType per line.

        When reaching the end of the last file, will start from the beginning
        again.

        Yields:
            One SampleBatch or MultiAgentBatch per line in all input files.
        """
        for path in self.files:
            file = self._try_open_file(path)
            while True:
                line = file.readline()
                if not line:
                    break
                batch = self._try_parse(line)
                if batch is None:
                    break
                yield batch

    def _postprocess_if_needed(self, batch: SampleBatchType) -> SampleBatchType:
        if not self.ioctx.config.get("postprocess_inputs"):
            return batch

        if isinstance(batch, SampleBatch):
            out = []
            for sub_batch in batch.split_by_episode():
                out.append(self.default_policy.postprocess_trajectory(sub_batch))
            return SampleBatch.concat_samples(out)
        else:
            # TODO(ekl) this is trickier since the alignments between agent
            #  trajectories in the episode are not available any more.
            raise NotImplementedError(
                "Postprocessing of multi-agent data not implemented yet."
            )

    def _try_open_file(self, path):
        if urlparse(path).scheme not in [""] + WINDOWS_DRIVES:
            if smart_open is None:
                raise ValueError(
                    "You must install the `smart_open` module to read "
                    "from URIs like {}".format(path)
                )
            ctx = smart_open
        else:
            # Allow shortcut for home directory ("~/" -> env[HOME]).
            if path.startswith("~/"):
                path = os.path.join(os.environ.get("HOME", ""), path[2:])

            # If path doesn't exist, try to interpret is as relative to the
            # rllib directory (located ../../ from this very module).
            path_orig = path
            if not os.path.exists(path):
                path = os.path.join(Path(__file__).parent.parent, path)
            if not os.path.exists(path):
                raise FileNotFoundError(f"Offline file {path_orig} not found!")

            # Unzip files, if necessary and re-point to extracted json file.
            if re.search("\\.zip$", path):
                with zipfile.ZipFile(path, "r") as zip_ref:
                    zip_ref.extractall(Path(path).parent)
                path = re.sub("\\.zip$", ".json", path)
                assert os.path.exists(path)
            ctx = open
        file = ctx(path, "r")
        return file

    def _try_parse(self, line: str) -> Optional[SampleBatchType]:
        line = line.strip()
        if not line:
            return None
        try:
            batch = self._from_json(line)
        except Exception:
            logger.exception(
                "Ignoring corrupt json record in {}: {}".format(self.cur_file, line)
            )
            return None

        # Clip actions (from any values into env's bounds), if necessary.
        cfg = self.ioctx.config
        if cfg.get("clip_actions") and self.ioctx.worker is not None:
            if isinstance(batch, SampleBatch):
                batch[SampleBatch.ACTIONS] = clip_action(
                    batch[SampleBatch.ACTIONS], self.default_policy.action_space_struct
                )
            else:
                for pid, b in batch.policy_batches.items():
                    b[SampleBatch.ACTIONS] = clip_action(
                        b[SampleBatch.ACTIONS],
                        self.ioctx.worker.policy_map[pid].action_space_struct,
                    )
        # Re-normalize actions (from env's bounds to zero-centered), if
        # necessary.
        if (
            cfg.get("actions_in_input_normalized") is False
            and self.ioctx.worker is not None
        ):

            # If we have a complex action space and actions were flattened
            # and we have to normalize -> Error.
            error_msg = (
                "Normalization of offline actions that are flattened is not "
                "supported! Make sure that you record actions into offline "
                "file with the `_disable_action_flattening=True` flag OR "
                "as already normalized (between -1.0 and 1.0) values. "
                "Also, when reading already normalized action values from "
                "offline files, make sure to set "
                "`actions_in_input_normalized=True` so that RLlib will not "
                "perform normalization on top."
            )

            if isinstance(batch, SampleBatch):
                pol = self.default_policy
                if isinstance(
                    pol.action_space_struct, (tuple, dict)
                ) and not pol.config.get("_disable_action_flattening"):
                    raise ValueError(error_msg)
                batch[SampleBatch.ACTIONS] = normalize_action(
                    batch[SampleBatch.ACTIONS], pol.action_space_struct
                )
            else:
                for pid, b in batch.policy_batches.items():
                    pol = self.policy_map[pid]
                    if isinstance(
                        pol.action_space_struct, (tuple, dict)
                    ) and not pol.config.get("_disable_action_flattening"):
                        raise ValueError(error_msg)
                    b[SampleBatch.ACTIONS] = normalize_action(
                        b[SampleBatch.ACTIONS],
                        self.ioctx.worker.policy_map[pid].action_space_struct,
                    )
        return batch

    def _next_line(self) -> str:
        if not self.cur_file:
            self.cur_file = self._next_file()
        line = self.cur_file.readline()
        tries = 0
        while not line and tries < 100:
            tries += 1
            if hasattr(self.cur_file, "close"):  # legacy smart_open impls
                self.cur_file.close()
            self.cur_file = self._next_file()
            line = self.cur_file.readline()
            if not line:
                logger.debug("Ignoring empty file {}".format(self.cur_file))
        if not line:
            raise ValueError(
                "Failed to read next line from files: {}".format(self.files)
            )
        return line

    def _next_file(self) -> FileType:
        # If this is the first time, we open a file, make sure all workers
        # start with a different one if possible.
        if self.cur_file is None and self.ioctx.worker is not None:
            idx = self.ioctx.worker.worker_index
            total = self.ioctx.worker.num_workers or 1
            path = self.files[round((len(self.files) - 1) * (idx / total))]
        # After the first file, pick all others randomly.
        else:
            path = random.choice(self.files)
        return self._try_open_file(path)

    def _from_json(self, data: str) -> SampleBatchType:
        if isinstance(data, bytes):  # smart_open S3 doesn't respect "r"
            data = data.decode("utf-8")
        json_data = json.loads(data)
        return from_json_data(json_data, self.ioctx.worker)

max_stars_repo_name: mgelbart/ray

max_issues_repo_name: mgelbart/ray

max_forks_repo_name: mgelbart/ray