text:
# Add const here when dev is over
NEARLIM_H_16 = 1.2
UPSAMPLIM_H_16 = 0.4

NEARLIM_RADIUS_16 = 3.5
UPSAMP_RADIUS_16 = 3.0

EXRAD_H = 0.05 # FMM exclusion radius / panel length
EDGELIM_H = 0.05 # Panel edge near criterion

# Bernstein radius computed in first quadrant
bernstein_radius_q1(z) = abs(z + sqrt(z^2 - 1))
q1(z) = abs(real(z)) + 1im*abs(imag(z))
bernstein_radius(z) = bernstein_radius_q1(q1(z))

rotate_and_scale(za, zb, z) = (za+zb)/(za-zb) + z*2/(zb-za)
complex2real(z) = [real(z), imag(z)]
tovec(x) = vec(x)
fromvec(x) = reshape(x, 2, :)

# Evaluate DLP with correction
function doublelayer(grid::DiscreteCurve,
                     density::Array{Float64},
                     target::Array{Float64},
                     alpha::Float64;
                     specquad=true)
    u = [0.0, 0.0]
    n = grid.panelorder
    coeffs = CurveDiscretization.map_panels(grid)
    for panel_idx = 1:grid.numpanels
        (up, rmin2, h) = dlp_panel_direct(grid, density, target,
                                          alpha, panel_idx)
        rmin = sqrt(rmin2)
        if specquad && rmin < 2*NEARLIM_H_16 *h*16/n
            # Nearby, investigate further:
            # Compute target point Bernstein radius
            z = target[1] + 1im*target[2]
            troot, _, converged = CurveDiscretization.invert_map(grid, coeffs, panel_idx, z,
                                                                 fall_back_to_initial=false)

            rho = bernstein_radius(troot)
            # Empirical thresholds for n=16

            #if specquad && rmin < NEARLIM_H_16 *h*16/n
            #upsample = rmin < UPSAMPLIM_H_16 *h*16/n
            if rho < NEARLIM_RADIUS_16^(16/n) # TODO: implement actual estimates here
                #@show troot, converged
                # Compute correction for this panel
                ucorr = dlp_panel_specquad(grid, density, target, alpha,
                                           panel_idx, troot)
                up = ucorr
            end # End near eval
        end
        u += up
    end # End panel loop
    return u
end

@inline function dlp_panel_direct(grid::DiscreteCurve,
                                  density::Array{Float64},
                                  target::Array{Float64},
                                  alpha::Float64,
                                  panel_idx::Int64)
    up = [0.0, 0.0]
    rmin2 = Inf
    h = 0.0
    for idx=grid.panelorder*(panel_idx-1) + (1:grid.panelorder)
        rvec1 = grid.points[1,idx]-target[1]
        rvec2 = grid.points[2,idx]-target[2]
        rnorm2 = rvec1*rvec1 + rvec2*rvec2
        rmin2 = min(rmin2, rnorm2)
        h += grid.dS[idx]
        ui = stresslet([rvec1, rvec2],
                       density[:,idx],
                       grid.normals[:,idx],
                       alpha)
        up += ui*grid.dS[idx]
    end
    return (up, rmin2, h)
end

function dlp_panel_specquad_weights(grid,za, zb, zt, zj,
                                    nuj, kappa,
                                    zpwj, zpj,
                                    troot,
                                    alpha::Float64,
                                    panel_idx::Int64,
                                    ifgrad::Bool,
                                    omega,
                                    DD,
                                    omega1, omega2,
                                    glpoints1, glweights1,
                                    glpoints2, glweights2,
                                    bcweights
                                    )
    n = length(glpoints1)
    n_sub = length(glpoints2)
    panel_length = sum(abs.(zpwj))
    dt = real(sum(zpwj ./ zpj))
    # Subdivide panel
    dt_max = LARGEALPHA_LIMIT_32*2/(alpha*panel_length)
    intervals = subdivide_interval_with_bisection(dt_max, troot, n_sub)
    # Compute weights on intervals
    S = zeros(2, 2*n)
    if ifgrad
        DUDXspec, DUDYspec = complex(zeros(1, 2*n)), complex(zeros(1, 2*n))
    end
    for idx = 1:length(intervals)-1
        ta = intervals[idx]
        tb = intervals[idx+1]
        troot_sub = (troot-ta)*2/(tb-ta) - 1.0
        # New source points, in original frame
        tsrc = ta + (1+glpoints2)*(tb-ta)/2
        # Interpolation to new source points
        Psub = bclag_interp_matrix(glpoints1, tsrc, bcweights)
        dt_sub = tb-ta
        if ta==-1 && tb==1
            za_sub, zb_sub = za, zb
        else
            za_sub = bclag_interp_eval(glpoints1, zj, ta, bcweights)
            zb_sub = bclag_interp_eval(glpoints1, zj, tb, bcweights)
        end
        nuj_sub = Psub*nuj
        zj_sub = Psub*zj
        zpj_sub = Psub*zpj
        kappa_sub = Psub*kappa
        zpwj_sub = zpj_sub .* glweights2 * dt/2 * (dt_sub/2)
        # nua and nub are not used, do this to make sure we notice if they are
        nua = nothing
        nub = nothing
        # Everything in place, do computation
        near_sub = bernstein_radius(troot_sub) < NEARLIM_RADIUS_16^(16/n_sub)
        if near_sub
            S_sub, DUDXspec_sub, DUDYspec_sub =
                complex_specquad_weights(zj_sub, zpwj_sub, nuj_sub, kappa_sub, za_sub, zb_sub, nua, nub,
                                         omega, omega1, omega2, zt, alpha,
                                         DD, n_sub, ifgrad)
        else
            S_sub = complex_direct_weights(zj_sub, zpwj_sub, nuj_sub, zt, alpha)
            if ifgrad
                DUDXspec_sub, DUDYspec_sub =
                    gradient_kernel_direct_weights(zj_sub, zpwj_sub, nuj_sub, zt, alpha)
            end
        end
        S += S_sub*kron(Psub, eye(2))
        if ifgrad
            DUDXspec += DUDXspec_sub*kron(Psub, eye(2))
            DUDYspec += DUDYspec_sub*kron(Psub, eye(2))
        end
    end
    # Join weights into one matrix and return
    if ifgrad
        W = vcat(S,
                 real(DUDXspec),
                 real(DUDYspec),
                 imag(DUDXspec),
                 imag(DUDYspec))
        return W
    else
        return S
    end
end

function complex_direct_weights(zj, zpwj, nuj, zt, alpha)
    # Build weights from panel to point
    n = length(zj)
    S = zeros(2, 2*n)
    for i=1:n
        r1 = real(zj[i]-zt)
        r2 = imag(zj[i]-zt)
        n1 = real(nuj[i])
        n2 = imag(nuj[i])
        dS = abs(zpwj[i])
        K11, K12, K21, K22 = dblkernel_noarrays(r1, r2, n1*dS, n2*dS, alpha)
        S[1, 1 + 2*(i-1)] = K11
        S[1, 2 + 2*(i-1)] = K12
        S[2, 1 + 2*(i-1)] = K21
        S[2, 2 + 2*(i-1)] = K22
    end
    return S
end

# Compute weights, using specquad if needed
function complex_specquad_weights(zj, zpwj, nuj, kappa, za, zb, nua, nub,
                                  omega, omega1, omega2, zt, alpha,
                                  DD, n, ifgrad)
    # Compute near eval coefficients
    wL, wC, wQ, wT = modifiedweights(za, zb, zt, zj)
    (wcorrL, wcorrC, wQ1, wQ2) =
        near_eval_weights(za, zb, zt, zj, nuj, zpwj, wL, wC, wQ)
    # Build weights from panel to point
    ucorr = [0.0, 0.0]
    S = zeros(2, 2*n)
    Kcorr = zeros(2,2)
    for i=1:n
        r1 = real(zj[i]-zt)
        r2 = imag(zj[i]-zt)
        n1 = real(nuj[i])
        n2 = imag(nuj[i])
        rnorm = sqrt(r1*r1 + r2*r2)
        rdotn = r1*n1 + r2*n2
        dS = abs(zpwj[i])
        # Direct quadrature:
        K11, K12, K21, K22 = dblkernel_noarrays(r1, r2, n1*dS, n2*dS, alpha)
        Kcorr[1,1] = K11
        Kcorr[1,2] = K12
        Kcorr[2,1] = K21
        Kcorr[2,2] = K22
        # Special quadrature correction:
        # Components of split stresslet
        uS11, uS21, uL11, uL21, uC11, uC21, uQ = stresslet_split_noarrays(r1, r2, 1.0, 0.0, n1, n2, alpha)
        uS12, uS22, uL12, uL22, uC12, uC22, uQ = stresslet_split_noarrays(r1, r2, 0.0, 1.0, n1, n2, alpha)

        # # Compute value as direct sum + corrections
        # # This could in theory avoid some cancellation,
        # # since uS is not explicitly used, and uL only
        # # multiplies a quadrature correction.

        # Actually, this avoids cancellation for targets close to nodes
        Kcorr[1,1] = uS11*dS
        Kcorr[1,2] = uS12*dS
        Kcorr[2,1] = uS21*dS
        Kcorr[2,2] = uS22*dS

        # Log
        wLi = wcorrL[i] #- log(rnorm)*dS
        Kcorr[1,1] += wLi*uL11
        Kcorr[2,1] += wLi*uL21
        Kcorr[1,2] += wLi*uL12
        Kcorr[2,2] += wLi*uL22
        # Cauchy
        wCi = wcorrC[i] #- rdotn/rnorm^2*dS
        Kcorr[1,1] += wCi*uC11
        Kcorr[2,1] += wCi*uC21
        Kcorr[1,2] += wCi*uC12
        Kcorr[2,2] += wCi*uC22
        # "Stresslet"
        A = wQ1[i]+wQ2[i]
        B = wQ1[i]-wQ2[i]
        Kcorr[1,1] += uQ*( real(A) )#- rdotn/rnorm^4*dS*r1*r1)
        Kcorr[1,2] += uQ*(-imag(B) )#- rdotn/rnorm^4*dS*r1*r2)
        Kcorr[2,1] += uQ*( imag(A) )#- rdotn/rnorm^4*dS*r1*r2)
        Kcorr[2,2] += uQ*( real(B) )#- rdotn/rnorm^4*dS*r2*r2)
        block_2x2_copy!(S, 1, i, Kcorr)
    end
    ## Gradient
    if ifgrad
        # Gradient spec weights
        partial_integration = true
        DUDX, DUDY =
            gradient_kernel_weights(zj, zpwj, nuj, kappa, za, zb, nua, nub,
                                    omega, omega1, omega2, zt, alpha,
                                    wcorrL, wC, wQ, wT, DD, partial_integration)
        return S, DUDX, DUDY
    else
        return S, [], []
    end
end

# On the fly computation of special quadrature
function dlp_panel_specquad(grid::DiscreteCurve,
                            density::Array{Float64},
                            target::Array{Float64},
                            alpha::Float64,
                            panel_idx::Int64,
                            troot)
    # For the correction weights we need complex arrays of
    # points on this panel
    n = grid.panelorder
    idx = n*(panel_idx-1) + (1:n)
    edges = grid.edges[:, panel_idx]
    za = edges[1] + 1im*edges[2]
    zb = edges[3] + 1im*edges[4]
    zt = target[1] + 1im*target[2]
    zj = grid.points[1,idx] + 1im*grid.points[2,idx]
    nuj = grid.normals[1,idx] + 1im*grid.normals[2,idx]
    weights = grid.weights[idx]
    zpwj = grid.dS[idx] .*nuj./1im
    omega = density[1,idx] + 1im*density[2,idx]
    # Resampling stuff
    Nsec = 32 # Fixed order
    glpoints1, glweights1 = gausslegendre(n)
    glpoints2, glweights2 = gausslegendre(Nsec)
    bcweights = bclag_interp_weights(glpoints1)
    # Panel subdivision
    # Subdivide panel if too long to resolve split kernel,
    # and target point in near quad region
    panel_length = sum(abs.(zpwj))
    dt_max = LARGEALPHA_LIMIT_32*2/(alpha*panel_length)
    intervals = subdivide_interval_with_bisection(dt_max, troot, Nsec)
    #intervals = subdivide_interval(dt_max, troot)
    # Go through subpanels and add contribution
    u = [0.0, 0.0]
    for idx = 1:length(intervals)-1
        ta = intervals[idx]
        tb = intervals[idx+1]
        troot_sec = (troot-ta)*2/(tb-ta) - 1.0
        # New source points, in original frame
        tsrc = ta + (1+glpoints2)*(tb-ta)/2
        # Interpolation to new source points
        Psec = bclag_interp_matrix(glpoints1, tsrc, bcweights)
        zj_sec = Psec*zj
        nuj_sec = Psec*nuj
        omega_sec = Psec*omega
        za_sec = bclag_interp_eval(glpoints1, zj, ta, bcweights)
        zb_sec = bclag_interp_eval(glpoints1, zj, tb, bcweights)
        zpwj_sec = glweights2*(tb-ta)/2 .* (Psec*(zpwj./glweights1))
        u += dlp_specquad_complex(za_sec, zb_sec, zt, zj_sec, nuj_sec, zpwj_sec, omega_sec, alpha, troot_sec)
    end
    return u
end

function dlp_specquad_complex(za, zb, zt, zj, nuj, zpwj, omega, alpha, troot)
    N = length(zj)
    rho = bernstein_radius(troot)
    near = rho < NEARLIM_RADIUS_16^(16/N)
    if near
        # Compute near eval coefficients
        (wcorrL, wcorrC, wQ1, wQ2) =
            near_eval_weights(za, zb, zt, zj, nuj, zpwj)
    end
    ucorr = [0.0, 0.0]
    for i=1:N
        r = complex2real(zj[i]-zt)
        f = complex2real(omega[i])
        n = complex2real(nuj[i])
        dS = abs(zpwj[i])
        u = stresslet(r, f, n, alpha)
        ucorr += u*dS
        (uS, uL, uC, uQ) = stresslet_split(r, f, n, alpha)
        rnorm = norm(r)
        rdotf = dot(r,f)
        rdotn = dot(r,n)
        if near
            # Compute value as direct sum + corrections
            # This could in theory avoid som cancellation,
            # since uS is not explicitly used, and uL only
            # multiplies a quadrature correction.
            # Log
            ucorr += (wcorrL[i] - log(rnorm)*dS)*uL
            # Cauchy
            ucorr += (wcorrC[i] - rdotn/rnorm^2*dS)*uC
            # "Stresslet"
            tmp = omega[i]*wQ1[i] + conj(omega[i])*wQ2[i]
            ucorr += (complex2real(tmp) - rdotf*r*rdotn/rnorm^4*dS)*uQ
        end
    end
    return ucorr
end

# FMM evaluation of DLP, with near corrections
function doublelayer_fast(grid::DiscreteCurve,
                          density::Array{Float64},
                          targets::Array{Float64},
                          alpha::Float64;
                          ifgrad::Bool=false,
                          specquad::Bool=true,
                          weights=[],
                          precomp=())
    # First compute double layer at all points using FMM
    ndS = grid.normals .* grid.dS'
    if isempty(precomp)
        fmmpars, tree, sorted_pts = doublelayer_precomp(grid, targets, alpha)
    else
        fmmpars, tree, sorted_pts = precomp
    end
    FMM = fmm_stresslet_targ(fmmpars, tree, sorted_pts, density, ndS, alpha, ifgrad=ifgrad)
    if ifgrad
        u, u1grad, u2grad = FMM
    else
        u, u1grad, u2grad  = FMM, Float64[], Float64[]
    end

    if specquad
        if length(weights)==0
            weights = doublelayer_near_weights(grid, targets, alpha, ifgrad=ifgrad)
        end
        doublelayer_near_apply!(grid, density, u, u1grad, u2grad, weights, alpha, ifgrad)
    end
    if ifgrad
        return u, u1grad, u2grad
    else
        return u
    end
end

function doublelayer_precomp(grid, targets, alpha)
    exrad = Array{Float64}(undef, grid.numpoints)
    panel_lengths = Array{Float64}(undef, grid.numpanels)
    for i=1:grid.numpanels
        h = 0.0
        for j=1:grid.panelorder
            h += grid.dS[(i-1)*grid.panelorder + j]
        end
        for j=1:grid.panelorder
            exrad[(i-1)*grid.panelorder + j] = EXRAD_H*h
        end
    end
    return fmm_stresslet_prep(grid.points, targets, alpha, exrad=exrad)
end

function doublelayer_near_weights(grid::DiscreteCurve,
                                  targets::Array{Float64},
                                  alpha::Float64;
                                  ifgrad::Bool=false)
    # For each point on curve, keep a list of tuples: (target_idx, weights)
    WeightElement = Tuple{Int64,Array{Float64,2}}
    weights = Array{Array{WeightElement}}(grid.numpanels)
    for i=1:grid.numpanels
        weights[i] = Tuple[]
    end
    # Find which points are "near"
    # Limits experimentally determined for n=16
    nearlim_h = NEARLIM_H_16^(16/grid.panelorder)
    _, hmax = minmax_panel(grid)
    nnb, dist = nearest_nb_R(grid, targets, nearlim_h*hmax)
    numeval = size(targets, 2)
    # Compute Legendre coeffs of panels
    coeffs = CurveDiscretization.map_panels(grid)
    # Precompute upsampling quantities
    n = grid.panelorder
    n2 = 2n
    glpoints1, glweights1 = gausslegendre(n)
    glpoints2, glweights2 = gausslegendre(n2)
    # Barycentric interpolation weights
    bcweights1 = bclag_interp_weights(glpoints1)
    bcweights2 = bclag_interp_weights(glpoints2)
    # Interpolation matrices for endpoints values
    B1n2 = vec(bclag_interp_matrix(glpoints2, -1)).'
    B2n2 = vec(bclag_interp_matrix(glpoints2,  1)).'
    # Real to complex matrix
    OMEGAn2 = kron(eye(n2), [1.0 1.0im])
    # On-panel differentiation matrix
    DDn2 = legendre_diff_matrix(n2)
    # Matrices to send
    OMEGA1n2 = B1n2*OMEGAn2
    OMEGA2n2 = B2n2*OMEGAn2
    # Main loop over points
    for i=1:numeval
        if nnb[i] != 0
            near_point_idx = nnb[i];
            near_point = grid.points[:, near_point_idx]
            rmin = dist[i]
            target = targets[:,i]
            zt = target[1] + 1im*target[2]
            near_panel_idx = Int64(ceil(near_point_idx/grid.panelorder))
            # Check three closest panels
            neighbors = (near_panel_idx,
                         grid.prevpanel[near_panel_idx],
                         grid.nextpanel[near_panel_idx])
            for nb_idx = 1:3
                # Load nodes on panel
                panel_idx = neighbors[nb_idx]
                idx = n*(panel_idx-1) + (1:n)
                zj = grid.points[1,idx] + 1im*grid.points[2,idx]
                panel_length = sum(grid.dS[idx])
                # Check root (preimage) of target for this panel
                troot, _, converged = CurveDiscretization.invert_map(grid, coeffs, panel_idx, zt,
                                                                     fall_back_to_initial=false)
                rho = bernstein_radius(troot)
                rmin2 = minimum(abs2.(zt - zj))
                # Do nothing if ouside of near quadrature limit AND not excluded by FMM
                if rho > NEARLIM_RADIUS_16^(16/grid.panelorder) && rmin >= (EXRAD_H*panel_length)^2
                    continue
                end
                # Load panel data
                za = grid.edges[1, panel_idx] + 1im*grid.edges[2, panel_idx]
                zb = grid.edges[3, panel_idx] + 1im*grid.edges[4, panel_idx]
                idx = n*(panel_idx-1) + (1:n)
                nuj = grid.normals[1,idx] + 1im*grid.normals[2,idx]
                kappa = grid.curvature[idx]
                zpwj = grid.dS[idx] .* nuj./1im
                zpj = zpwj ./ grid.weights[idx]
                # Get direct weights, for subtraction
                if ifgrad
                    DUDXdir, DUDYdir = gradient_kernel_direct_weights(zj, zpwj, nuj, zt, alpha)
                end
                D = complex_direct_weights(zj, zpwj, nuj, zt, alpha)
                # Corner case: Points too close to discretization node can cause cancellation errors
                # FMM avoids points within EXRAD_H*panel_length, so no subtraction needed
                for j=1:n
                    if abs2(zj[j]-zt) < ( EXRAD_H*panel_length )^2
                        D[1:2, (1:2) + 2*(j-1)] = 0.0
                        if ifgrad
                            DUDXdir[1, (1:2) + 2*(j-1)] = 0.0
                            DUDYdir[1, (1:2) + 2*(j-1)] = 0.0
                        end
                    end
                end
                if ifgrad
                    Wsub = vcat(-D,
                                real(-DUDXdir),
                                real(-DUDYdir),
                                imag(-DUDXdir),
                                imag(-DUDYdir))
                else
                    Wsub = -D
                end
                # Check if too close to an edge
                prev_panel_length = sum(grid.dS[n*(grid.prevpanel[panel_idx]-1) + (1:n)])
                next_panel_length = sum(grid.dS[n*(grid.nextpanel[panel_idx]-1) + (1:n)])
                left_edge_rad  = EDGELIM_H * min(panel_length, prev_panel_length)
                right_edge_rad = EDGELIM_H * min(panel_length, next_panel_length)
                if abs(zt-za) > left_edge_rad && abs(zt-zb) > right_edge_rad
                    # Not close to an edge, compute special quadrature as usual
                    Wspec = dlp_panel_specquad_weights(grid,za, zb, zt, zj,
                                                       nuj, kappa,
                                                       zpwj, zpj,
                                                       troot, alpha, panel_idx, ifgrad,
                                                       OMEGAn2, DDn2,
                                                       OMEGA1n2, OMEGA2n2,
                                                       glpoints1, glweights1,
                                                       glpoints2, glweights2,
                                                       bcweights1
                                                       )
                    push!(weights[panel_idx], (i, Wspec+Wsub))
                else
                    if abs(zt-za) <= left_edge_rad
                        # Close to left edge, just subtract direct interaction
                        # Neighbor will deal with special quadrature
                        push!(weights[panel_idx], (i, Wsub))
                        continue
                    end
                    # Close to right edge, merge with right panel
                    next_panel_idx = grid.nextpanel[panel_idx]
                    zb = grid.edges[3, next_panel_idx] + 1im*grid.edges[4, next_panel_idx]
                    # index list spanning both panel
                    idx_wide = vcat(n*(panel_idx-1) + (1:n), n*(next_panel_idx-1) + (1:n))
                    # panels may be of different length in parameter
                    dt1 = grid.t_edges[2, panel_idx] - grid.t_edges[1, panel_idx]
                    dt2 = grid.t_edges[2, next_panel_idx] - grid.t_edges[1, next_panel_idx]
                    dt = dt1+dt2
                    ta, tb, tc = -1, dt1*2/dt-1, 1
                    g1 = ta + (1+glpoints1)*(tb-ta)/2
                    g2 = tb + (1+glpoints1)*(tc-tb)/2
                    # GL nodes of both panels, scaled into [-1,1]
                    glpoints_wide = vcat(g1, g2)
                    # Now we'll interpolate to 2*n nodes,
                    # find which fall into the first panel
                    pts1 = count(glpoints2 .< tb)
                    # Create interpolation from the two panels to nodes on new
                    P1 = bclag_interp_matrix(g1, glpoints2[1:pts1])
                    P2 = bclag_interp_matrix(g2, glpoints2[pts1+1:n2])
                    Pwide = zeros(n2,2n)
                    Pwide[1:pts1, 1:n] = P1
                    Pwide[pts1+1:n2, n+(1:n)] = P2
                    # Load boundary data for both panels
                    zj_wide = grid.points[1,idx_wide] + 1im*grid.points[2,idx_wide]
                    nuj_wide = grid.normals[1,idx_wide] + 1im*grid.normals[2,idx_wide]
                    kappa_wide = grid.curvature[idx_wide]
                    zpwj_wide = grid.dS[idx_wide] .* nuj_wide./1im
                    zpj_wide = zpwj_wide ./ grid.weights[idx_wide]
                    # Interpolate to merged panel
                    zj = Pwide*zj_wide
                    nuj = Pwide*nuj_wide
                    kappa = Pwide*kappa_wide
                    zpj = Pwide*zpj_wide
                    zpwj = zpj*dt/2.*glweights2
                    # Find root of zt for merged panel (safer than scaling)
                    coeff = legendre_matrix(n2)*rotate_and_scale(za, zb, zj)
                    troot_guess = 1im*imag(troot)/2
                    troot_new, _, _ = newton_legendre(coeff, rotate_and_scale(za, zb, zt),
                                                      troot_guess, 30, 1e-13)
                    # Get specquad weights for merged panel
                    Wwide = dlp_panel_specquad_weights(grid,za, zb, zt, zj,
                                                       nuj, kappa,
                                                       zpwj, zpj,
                                                       troot_new, alpha, panel_idx, ifgrad,
                                                       OMEGAn2, DDn2,
                                                       OMEGA1n2, OMEGA2n2,
                                                       glpoints2, glweights2,
                                                       glpoints2, glweights2,
                                                       bcweights2
                                                       )
                    # Multiply weights with interpolation and split to the right source panels
                    # Also subtract direct interaction from current panel
                    Wwide = Wwide * kron(Pwide, eye(2))
                    Wwide1 = Wwide[:,1:2*n]
                    Wwide2 = Wwide[:,2*n+1:4*n]
                    push!(weights[panel_idx], (i, Wwide1+Wsub))
                    push!(weights[next_panel_idx], (i, Wwide2))
                end
            end
        end
    end
    return weights
end

# Apply precomputed near eval weights
function doublelayer_near_apply!(grid::DiscreteCurve,
                                 density::Array{Float64},
                                 u, u1grad, u2grad,
                                 weights,
                                 alpha::Float64,
                                 ifgrad)
    n = grid.panelorder
    for panel_idx = 1:grid.numpanels
        for i=1:length(weights[panel_idx])
            point_idx = weights[panel_idx][i][1]
            Wspec = weights[panel_idx][i][2]
            uspec1, uspec2 = 0.0, 0.0
            idx0 = 2*n*(panel_idx-1)
            if ifgrad
                uspec1, uspec2, du1dx, du1dy, du2dx, du2dy = matvec_6xN(Wspec, density, idx0, 2*n)
                u1grad[1,point_idx] += du1dx
                u1grad[2,point_idx] += du1dy
                u2grad[1,point_idx] += du2dx
                u2grad[2,point_idx] += du2dy
            else
                uspec1, uspec2 = matvec_2xN(Wspec, density, idx0, 2*n)
            end
            u[1,point_idx] += uspec1
            u[2,point_idx] += uspec2
        end
    end
end
# Helper functions for the above
function matvec_2xN(A, x, idx, n)
    f1,f2 = 0.0, 0.0
    @simd for i=1:n
        @inbounds f1 += A[1,i]*x[idx+i]
        @inbounds f2 += A[2,i]*x[idx+i]
    end
    return f1,f2
end
function matvec_6xN(A, x, idx, n)
    f1,f2,f3,f4,f5,f6 = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
    @simd for i=1:n
        @inbounds f1 += A[1,i]*x[idx+i]
        @inbounds f2 += A[2,i]*x[idx+i]
        @inbounds f3 += A[3,i]*x[idx+i]
        @inbounds f4 += A[4,i]*x[idx+i]
        @inbounds f5 += A[5,i]*x[idx+i]
        @inbounds f6 += A[6,i]*x[idx+i]
    end
    return f1,f2,f3,f4,f5,f6
end

# DLP on curve
# !This does not include correction weights for log kernel
function doublelayer_self(grid::DiscreteCurve,
                          density::Array{Float64},
                          alpha::Float64;
                          zerolimit::Bool=true)
    npoints = grid.numpoints
    u = zeros(2,npoints)
    for i=1:npoints
        ui = 0
        zi = grid.points[:,i]
        for j = 1:i-1
            rij = grid.points[:, j] - zi
            ui += stresslet(rij, density[:,j],
                            grid.normals[:,j], alpha) * grid.dS[j]
        end
        for j = i+1:npoints
            rij = grid.points[:, j] - zi
            ui += stresslet(rij, density[:,j],
                            grid.normals[:,j], alpha) * grid.dS[j]
        end
        u[:,i] = ui
    end
    if zerolimit
        for i=1:npoints
            # Self
            ni = grid.normals[:,i]
            ti = [ni[2],-ni[1]]
            T2C0 = -1/pi
            u[:,i] += -T2C0*grid.curvature[i]*0.5*grid.dS[i] *
                dot(ti, density[:,i])*ti
        end
    end
    return u
end

"Return dense system matrix"
function system_matrix(grid::DiscreteCurve,
                       alpha::Float64)
    nvecdS = vec(grid.dS' .* grid.normals)
    nvec = vec(grid.normals)
    L = sum(grid.dS)
    Dmat = doublelayer_matrix(grid, alpha)
    ## Add add rank correction
    # NNT = nvec*nvecdS' / L
    # Dmat += NNT
    @. Dmat += nvec*nvecdS' / L
    # Add identity
    #Dmat += eye(2*grid.numpoints)/2
    for i=1:2*grid.numpoints
        Dmat[i,i] += 1/2
    end
    return  Dmat
end

"Return system matvec function handle"
function system_matvec(grid::DiscreteCurve,
                       alpha::Float64)
    ndS = grid.normals .* grid.dS'
    L = sum(grid.dS)
    fmmpars, tree, sorted_pts = fmm_stresslet_prep(grid.points, grid.points, alpha)
    corrmat = correction_block_matrix(grid, alpha)
    return x -> x/2 + corrmat*x +
        tovec( fmm_stresslet_self(fmmpars, tree, sorted_pts, fromvec(x), ndS, alpha) ) +
        tovec(grid.normals)*dot(x, tovec(ndS))/L
end

"Return block matrix that corrects a naive direct sum"
function correction_block_matrix(grid::DiscreteCurve,
                                 alpha::Float64)
    return specquad_block_matrix(grid, alpha; correction=true)
end

"Return block matrix containing elements of full system matrix for
neighboring panels"
function system_block_matrix(grid::DiscreteCurve,
                                 alpha::Float64)
    return specquad_block_matrix(grid, alpha;
                                 correction=false,
                                 system_matrix = true)
end

function specquad_block_matrix(grid::DiscreteCurve,
                               alpha::Float64;
                               correction::Bool=false,
                               system_matrix::Bool=false)
    # Compute sparse block matrix that contains all interactions requiring
    # special quadrature
    #
    # correction=true:
    # Subtract naive point-to-point sum (excluding self interactions) from
    # blocks, such that results can be used to correct an FMM
    #
    # system_matrix=true:
    # Include dbl layer jump and null space correction

    npoints = grid.numpoints
    panelorder = grid.panelorder
    L = sum(grid.dS)
    # Allocate arrays for constructing sparse matrix
    # only make room for nearest neighbor interactions,
    # this will have to grow for close geometries
    blocksize = 2*panelorder
    nzblocks = 3*grid.numpanels
    indptr = Array{Int64}(grid.numpanels+1)
    indices = Array{Int64}(nzblocks)
    data = zeros(nzblocks, blocksize, blocksize)
    B1 = Array{Float64}(blocksize, blocksize) # Panel-to-panel block
    B2 = Array{Float64}(blocksize, blocksize) # Panel-to-panel block
    idx = 1;
    for trg_panel = 1:grid.numpanels
        indptr[trg_panel] = idx -1 # -1 for Python indexing
        neighbors = (grid.prevpanel[trg_panel],
                     trg_panel,
                     grid.nextpanel[trg_panel])
        for src_panel = neighbors
            if correction
                # Compute block that represents naive sum
                p2p_block!(B1, grid, src_panel, trg_panel, alpha, specquad=false, upsampling=false)
                if src_panel==trg_panel
                    for i=1:panelorder
                        B1[2*(i-1)+(1:2), 2*(i-1)+(1:2)] = 0.0
                    end
                end
            end
            # Compute correct block
            p2p_block!(B2, grid, src_panel, trg_panel, alpha)
            if system_matrix
                # Add null space correction
                for i=1:panelorder
                    trg_idx = i + panelorder*(trg_panel-1)
                    ni1, ni2 = grid.normals[1:2, trg_idx]
                    for j=1:panelorder
                        src_idx = j + panelorder*(src_panel-1)
                        nj1, nj2 = grid.normals[1:2, src_idx]
                        dSoverL = grid.dS[src_idx] / L
                        B2[2*(i-1)+1, 2*(j-1)+1] += ni1*nj1*dSoverL
                        B2[2*(i-1)+2, 2*(j-1)+1] += ni2*nj1*dSoverL
                        B2[2*(i-1)+1, 2*(j-1)+2] += ni1*nj2*dSoverL
                        B2[2*(i-1)+2, 2*(j-1)+2] += ni2*nj2*dSoverL
                    end
                end
                if src_panel==trg_panel
                    # Add jump
                    for i=1:2*panelorder
                        B2[i,i] += 1.0/2
                    end
                end
            end
            # Insert into structure
            indices[idx] = src_panel -1
            data[idx,:,:] = B2
            if correction
                @. data[idx,:,:] -= B1
            end
            idx += 1
        end
        indptr[trg_panel+1] = idx -1
    end
    # Assemble BSR matrix
    s = pyimport("scipy.sparse")
    A = BSR_Matrix(s[:bsr_matrix]((data, indices, indptr),
                                  blocksize=(blocksize, blocksize)))
    return A
end

function doublelayer_matrix(grid::DiscreteCurve,
                            alpha::Float64)
    npoints = grid.numpoints
    panelorder = grid.panelorder
    A = Array{Float64}(2*npoints, 2*npoints) # Output matrix
    B = Array{Float64}(2*panelorder, 2*panelorder) # Panel-to-panel block
    for src_panel = 1:grid.numpanels
        for trg_panel = 1:grid.numpanels
            p2p_block!(B, grid, src_panel, trg_panel, alpha)
            block_nxn_copy!(A, trg_panel, src_panel, B, 2*panelorder)
        end
    end
    return A
end

@inline function p2p_block!(Bin, grid, src_panel, trg_panel, alpha;
                            upsampling=false,
                            specquad=true)
    panelorder = numpoints = grid.panelorder
    self_panel = (trg_panel==src_panel)
    prev_panel = (trg_panel==grid.prevpanel[src_panel])
    next_panel = (trg_panel==grid.nextpanel[src_panel])
    if self_panel || prev_panel || next_panel
        # The upsample flag corresponds to Scheme B in Helsing & Holst,
        # which computes panel interactions using upsampled grids.
        # This improves accuracy for large alpha when subdivision are not used
        # Subdivisions make upsampling obsolete, so turned off by default
        upsample = upsampling
        near = specquad
        if self_panel
            d = 0.0
        elseif prev_panel
            d = 1.0
        else
            d = -1.0
        end
    else
        upsample = false
        near = false
    end
    # Setup buffer block
    if upsample
        B = Array{Float64}(2*2*panelorder, 2*2*panelorder)
    else
        B = Bin;
    end
    # Load sources
    src_idx = (1:panelorder) + panelorder*(src_panel-1)
    zsrc = grid.points[1:2, src_idx]
    nsrc = grid.normals[1:2, src_idx]
    wsrc = grid.dS[src_idx]
    # Load targets
    trg_idx = (1:panelorder) + panelorder*(trg_panel-1)
    ztrg = grid.points[1:2, trg_idx]
    if upsample
        P = upsampling_matrix(panelorder, 2*panelorder)
        Q = upsampling_matrix(2*panelorder, panelorder)
        zsrc = (P*zsrc')'
        nsrc = (P*nsrc')'
        ztrg = (P*ztrg')'
        numpoints *= 2
        # Take GL weights out of interpolation
        _, glweights1 = gausslegendre(panelorder)
        _, glweights2 = gausslegendre(2*panelorder)
        wsrc = glweights2 .* (P*(wsrc ./ glweights1))
    end
    # Compute interaction block between src and trg panels
    if near
        ## Near, need some kind of special quadrature
        # Transformations
        trg_dt = grid.t_edges[2, trg_panel] - grid.t_edges[1, trg_panel]
        src_dt = grid.t_edges[2, src_panel] - grid.t_edges[1, src_panel]
        scale = trg_dt / src_dt
        trans = -d*(1+scale)
        # If panel is short enough w.r.t. alpha, use direct computation
        # If panel is too long, compute quadrature using subdivisions of src panel
        panel_max = LARGEALPHA_LIMIT_16/alpha
        panel_length = sum(wsrc)
        if panel_length <= panel_max
            ## Direct
            dblkernel_block!(B, ztrg, zsrc, nsrc, wsrc, alpha)
            add_logcorr_block!(B, ztrg, zsrc, nsrc, wsrc, scale, trans, alpha)
            if self_panel
                # Correct 2x2 blocks on diagonal, is now NaN
                ksrc = grid.curvature[src_idx] # Curvature
                if upsample
                    ksrc = P*ksrc
                end
                for i=1:numpoints
                    ni = nsrc[:,i]
                    ti = [ni[2],-ni[1]]
                    T2C0 = -1/pi
                    K = (-T2C0*ksrc[i]*0.5*wsrc[i])*(ti*ti')
                    block_2x2_copy!(B, i, i, K)
                end
            end
        else
            ## Using subdivisions
            fill!(B, 0.0) # Need to initialize, since we add to it in loop
            # Seems subdivisions are fine with using 16 points
            numsec = 16
            gcpoints, gcweights = gausslegendre(numsec)
            glpoints, glweights = gausslegendre(numpoints)
            bcweights = bclag_interp_weights(glpoints)
            for i=1:numpoints # Target point loop
                # Location of target in parametrization of source panel
                ttrg = trans + scale*glpoints[i]
                # Allowed panel size, in parametrization
                hsub = LARGEALPHA_LIMIT_16*2/(alpha*panel_length)
                # Don't let min subpanel size be too close to 2.0
                #hsub = min(1.9, hsub)
                # Get subdivision of source interval
                #intervals = subdivide_interval(hsub, ttrg)
                intervals = subdivide_interval_with_bisection(hsub, ttrg, numsec)

                R = zeros(2, 2*numsec) # Output row
                for idx = 1:length(intervals)-1
                    fill!(R, 0.0)
                    ta = intervals[idx]
                    tb = intervals[idx+1]
                    # New source points, in original frame
                    tsrc = ta + (1+gcpoints)*(tb-ta)/2
                    # Interpolation to new source points
                    Psec = bclag_interp_matrix(glpoints, tsrc, bcweights)
                    zsrc_sec = (Psec*zsrc')'
                    nsrc_sec = (Psec*nsrc')'
                    wsrc_sec = gcweights*(tb-ta)/2 .* (Psec*(wsrc ./ glweights))
                    # Direct double layer interaction
                    dblkernel_block!(R, ztrg[1:2,i], zsrc_sec, nsrc_sec, wsrc_sec, alpha)
                    # Target point, in new frame
                    ttrg_sec = (ttrg-ta)*2/(tb-ta) - 1.0
                    # tdist = Distance in parametrization of src panel
                    tdist = abs.(ttrg_sec - gcpoints)
                    # Very empirical limit for when neighbors need special quad due to log kernel:
                    # 0.2 for 32, 0.4 for 16
                    # Too small gives nearly sing error
                    # Too large makes special quad inaccurate, though not very sensitive
                    if minimum(tdist)<0.2*32/numsec
                        # Compute Helsing log-correction and apply to row
                        WfrakL = WfrakLcomp(ttrg_sec, gcpoints)
                        KL = Array{Float64}(2,2)
                        for j=1:numsec
                            nb_wcorr = WfrakL[1, j]/gcweights[j] - log(tdist[j])
                            zi = ztrg[1:2, i]
                            zj = zsrc_sec[1:2, j]
                            nj = nsrc_sec[1:2, j]
                            wj = wsrc_sec[j]
                            rij = zj.-zi
                            dblkernel_log!(rij, nj, alpha, KL)
                            block_2x2_add!(R, 1, j, KL.*(nb_wcorr*wj))
                        end
                    end
                    # Add row*interpolation to output block
                    B[2*(i-1)+(1:2), :] += R*kron(Psec, eye(2))
                end
            end

        end
    else
        # Not near, direct interaction
        dblkernel_block!(B, ztrg, zsrc, nsrc, wsrc, alpha)
    end
    if upsample
        Bin .= kron(Q, eye(2)) * B * kron(P, eye(2))
    end
end

@inline function dblkernel_block!(B, ztrg, zsrc, nsrc, wsrc, alpha)
    K = Array{Float64}(2,2)
    rji = Array{Float64}(2)
    nj = Array{Float64}(2)
    numtrg = size(ztrg, 2)
    numsrc = size(zsrc, 2)
    for j=1:numsrc
        @. nj = nsrc[:,j]*wsrc[j]
        for i=1:numtrg
            @. rji = zsrc[:,j]-ztrg[:,i]
            dblkernel!(rji, nj, alpha, K)
            block_2x2_copy!(B, i, j, K)
        end
    end
end

@inline function block_nxn_copy!(A, i, j, K, n)
    for k=1:n
        for l=1:n
            A[n*(i-1)+k, n*(j-1)+l] = K[k,l]
        end
    end
end

@inline function block_2x2_copy!(A, i, j, K)
    for k=1:2
        for l=1:2
            A[2*(i-1)+k, 2*(j-1)+l] = K[k,l]
        end
    end
end

@inline function block_2x2_add!(A, i, j, K)
    for k=1:2
        for l=1:2
            A[2*(i-1)+k, 2*(j-1)+l] += K[k,l]
        end
    end
end


text:
#
# update velocities using berendsen bath
#

function update_velocities!(v :: Array{Float64} ,kavg :: Float64,
                            f :: Array{Float64}, flast :: Array{Float64},
                            input :: MDInput)

  n = input.n
  dt = input.dt
  vx = 0.
  vy = 0.
  lambda = sqrt( 1 + (dt/input.tau)*(input.kavg_target/kavg-1) )
  for i in 1:n
    v[i,1] = v[i,1]*lambda
    v[i,2] = v[i,2]*lambda
    vx = vx + v[i,1]
    vy = vy + v[i,2]
  end
  # Remove drift
  vx = vx / n
  vy = vy / n
  for i in 1:n
    v[i,1] = v[i,1] - vx
    v[i,2] = v[i,2] - vy
  end

  # Updating velocities
  for i in 1:n
    v[i,1] = v[i,1] + 0.5*(f[i,1]+flast[i,1])*dt
    v[i,2] = v[i,2] + 0.5*(f[i,2]+flast[i,2])*dt
  end

end

ext: jl


text:
# ---------------------------------------------------------------------------- #
#
#   meshF3D.jl
#
#   Type for 3D frame meshes
#   Inherits from the abstract "meshF" type
#
#   intrico
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    MeshF3D

Mesh3D type:
Type for 3D meshes consisting of tetrahedrons.
"""
struct MeshF3D <: MeshF

    dim::Int64                  # Dimension of the problem

    n::Int64                    # Number of nodes

    p::Matrix{Float64}          # Nodal locations
    t::Matrix{Int64}            # Triangle - node connectivity
    t2f::Matrix{Int64}          # Triangle - face connectivity
    f::Matrix{Int64}            # face - node/triangle connectivity
    fb::Matrix{Int64}           # Boundary face info

    e::Matrix{Int64}            # Edge connectivity
    e2f::Vector{Vector{Int64}}  # Edge to face connectivity
    f2e::Matrix{Int64}          # Face to edge connectivity
    nodes::Array{Float64,3}     # Nodes on which solution is evaluated -- needed to
                                #  compute Jacobian on face

    n2e::Vector{Vector{Int64}}  # Node to edge connectivity

    boundtags::Vector{String}   # boundary tags

end

"""
    MeshF3D( mesh::divido.Mesh3D )

Constructor that generates frame from `mesh`.
"""
function MeshF( mesh::divido.Mesh3D )

    e, e2f, f2e  = genEdgesF3D( mesh.f )

    n2e = genNEconnec3D( e, mesh.n )

    MeshF3D( 3, mesh.n, mesh.p, mesh.t, mesh.t2f, mesh.f, mesh.fb, e, e2f, f2e,
             mesh.nodes, n2e, mesh.boundtags )

end

"""
    MeshF3D( name::String )

Constructor for frame that reads in a mesh from `name`.
"""
function MeshF3D( name::String )

    if name[end-3:end] == ".su2"
        (p_, t_, bel_, tags_ ) = divido.readSU2_3D( name )
    elseif name[end-4:end] == ".mesh" # FEFLOA
        (p_, t_, bel_, tags_ ) = divido.readFEFLOA_3D( name )
    else
        error("Unknown mesh type")
    end

    porder_ = divido.P1()

    (f_, t2f_, nodes_, ploc_, tloc_, trorder_, fb_) = divido.genmesh3D( porder_, p_, t_, bel_ )

    n_ = size( p_, 1 )

    e_, e2f_, f2e_ = genEdgesF3D( mesh.f )

    n2e = genNEconnec3D( e_, n_ )

    MeshF3D( 3, n_, p_, t_, t2f_, f_, fb_, e_, e2f_, f2e_, nodes_, n2e, tags_ )

end

"""
genEdgesF3D( f::Matrix{Int64} )

Generates edge connectivity given face connectivity (`f`).
"""
    function genEdgesF3D( f::Matrix{Int64} )

    nf = size(f,1)

    edges = vcat( f[:,[1,2]], f[:,[2,3]], f[:,[3,1]] ) # This holds all edges (but multiple copies)
    faces  = vcat( 1:nf, 1:nf, 1:nf )

    boundsf =  f[:,5]
    boundsf[ boundsf .> 0 ] = 0
    boundsf = -boundsf
    bounds  =  vcat( boundsf, boundsf, boundsf )
    ne = size(edges,1)

    boundsUni = fill( 0, ne, 2 )

    # sort in ascending order
    edges = sort( edges, 2 )

    # This index links to the first unique index in the array edges
    ix  = groupslices( edges, 1 )

    # Find unique vector index
    # Need to find an index that links from edges to the unique edges (without gaps)
    jx = Array{Int64}( size(ix) )
    jx[1] = 0
    lind = 0 # maximum unique index
    mind = 0 # maximum index in edges that had unique index
    ne = 0  # number of unique elements

    indUni = fill( false, size(ix,1) ) # index of unique elements

    for ii = 1:length(ix)
        temp    = ix[ii] - lind - 1
        temp2   = ix[ii] - mind - 1
        jx[ii]  = temp
        mind    = max(mind,ix[ii])

        if temp2 >= 0 # Found a unique index
            lind = lind + 1
            ne  = ne + 1
            indUni[ii] = true
            # add boundary information
            boundsUni[ii,1] = bounds[ii]
        end

        if temp2 < 0 # Need to link to the unique index
            jx[ii] = jx[ix[ii]]
            # add boundary information
            if bounds[ii] > 0
                if boundsUni[ ix[ii] ,1] > 0
                    boundsUni[ ix[ii] ,2] = bounds[ii]
                else
                    boundsUni[ ix[ii] ,1] = bounds[ii]
                end
            end
        end

    end

    e = fill( 0::Int64, ne, 4 ) # Last two are for boundary edges
              # We need two because some faces share two
              # boundaries

    e[:,1:2] = edges[indUni,:]

    # add boundary information
    e[:,3:4]   = sort( boundsUni[indUni,:], 2 )

    # make e2f, f2e
    e2f = Vector{Vector{Int64}}( ne )
    for jj in 1:ne
        e2f[jj] = fill( 0, 0 )
    end
    f2e = Matrix{Int64}( nf, 3 )
    for jj = 1:3*nf

        ind  = faces[jj]
        indf = Int64(ceil(jj / nf))

        inde = ix[jj] - jx[jj]

        append!( e2f[inde], ind )

        f2e[ ind, indf ] = inde

    end

    return e, e2f, f2e

end

"""
    genNEconnec3D( e::Matrix{Int64}, n::Int64 )

Generates node to edge connectivity given edge connectivity (`e`).
"""
function genNEconnec3D( e::Matrix{Int64}, n::Int64 )

    n2e = Vector{Vector{Int64}}(n)
    for ii in 1:n
        n2e[ii] = fill( 0, 0 )
    end

    for jj in 1:size(e,1)

        ii1 = e[jj,1]
        ii2 = e[jj,2]

        if isempty( n2e[ii1] .== jj ) || !any( n2e[ii1] .== jj )
            append!(n2e[ii1], jj)
        end
        if isempty( n2e[ii2] .== jj ) || !any( n2e[ii2] .== jj )
            append!(n2e[ii2], jj)
        end

    end

    for kk in 1:n
        n2e[kk] = sort( n2e[kk] )
    end

    return n2e

end


text:
import Glucose
using Glucose: optim_all_vars_exo
using SuParameters
using RunTools
using Omega: normal, uniform, ciid
using DiffEqFlux, OrdinaryDiffEq, Flux, Optim, Plots, AdvancedHMC, MCMCChains
using JLD, StatsPlots

function runparams()
  φ = Params()
  φ.simname = "train"
  φ.runname = ciid(randrunname)
  φ.tags = ["odeoptim", "firsttry"]
  φ.logdir = ciid(ω -> logdir(runname = φ.runname(ω), tags = φ.tags))
  φ.runfile = @__FILE__
  φ.gitinfo = current_commit(@__FILE__)
  φ
end

"Optimization Parameters"
function optparams()
  Params((binsize = uniform([10]),
          numwarmup = uniform([500]),
          numsamples = uniform([500]),
          initstepsize = uniform([0.45]),
          odesolver = uniform([Tsit5]) # Rosenbrock23, Trapezoid
  ))
end

function train(params)
  samples, losses = bayes_all_vars_exo(params.binsize,
                                       params.numwarmup,
                                       params.numsamples,
                                       params.initstepsize,
                                       log_dir = params.logdir)
end

function allparams()
  φ = Params()
  merge(φ, runparams(), optparams()) # netparams()
end

"Run with julia -L bayesparams.jl -E 'hyper(;)' -- --queue"
function hyper(; params = Params(), n = 2)
  params_ = allparams()
  paramsamples = rand(params_, n)
  display.(paramsamples)
  control(train, paramsamples)
end

function testrun()
  params_ = allparams()
  train(rand(params_))
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
immutable Stop <: Exception end

function run_task(f)
   try f()
   catch e
       if isa(e, Stop)
           println("Your program was gracefully terminated.")
       else
           rethrow()
       end
   end
end

function f()
   throw(Stop())
end

function g()
   f()
end

run_task() do
   g()
end
println()

run_task() do
   error("This is bad.")
end


text:
samesize_float = Dict(Int8=>Float16, UInt8=>Float16, Int16=>Float16, UInt16=>Float16,
                      Int32=>Float32, UInt32=>Float32, Int64=>Float64, UInt64=>Float64)

# GlobalPool(x, aggr, batch, size=nothing) # aggr=sum, mean, max
#
# TopKPool()

# struct MaxPool
#
# end
#
# function MaxPool(adj::AbstractMatrix, ch::Pair{<:Integer,<:Integer}, σ = identity)
#     MaxPool()
# end
#
# (m::MaxPool)(X::AbstractMatrix) = maxpool(m.cluster, X)
#
#
#
# struct MeanPool
#
# end
#
# (m::MeanPool)(X::AbstractMatrix) = meanpool(m.cluster, X)

function sumpool(cluster::Array{Int}, X::Array{T}) where {T<:Real}
    dims = _pooling_dim_check(cluster, X)
    c = length(Set(cluster))
    Y = zeros(T, dims..., c)
    scatter_add!(Y, X, cluster)
    Y
end

function subpool(cluster::Array{Int}, X::Array{T}) where {T<:Real}
    dims = _pooling_dim_check(cluster, X)
    c = length(Set(cluster))
    Y = zeros(T, dims..., c)
    scatter_sub!(Y, X, cluster)
    Y
end

function prodpool(cluster::Array{Int}, X::Array{T}) where {T<:Real}
    dims = _pooling_dim_check(cluster, X)
    c = length(Set(cluster))
    Y = ones(T, dims..., c)
    scatter_mul!(Y, X, cluster)
    Y
end

function divpool(cluster::Array{Int}, X::Array{T}) where {T<:Real}
    dims = _pooling_dim_check(cluster, X)
    c = length(Set(cluster))
    FT = (T <: Integer) ? samesize_float[T] : T
    Y = ones(FT, dims..., c)
    scatter_div!(Y, FT.(X), cluster)
    Y
end

function maxpool(cluster::Array{Int}, X::Array{T}) where {T<:Real}
    dims = _pooling_dim_check(cluster, X)
    c = length(Set(cluster))
    Y = fill(typemin(T), dims..., c)
    scatter_max!(Y, X, cluster)
    Y
end

function minpool(cluster::Array{Int}, X::Array{T}) where {T<:Real}
    dims = _pooling_dim_check(cluster, X)
    c = length(Set(cluster))
    Y = fill(typemax(T), dims..., c)
    scatter_min!(Y, X, cluster)
    Y
end

function meanpool(cluster::Array{Int}, X::Array{T}) where {T<:Real}
    dims = _pooling_dim_check(cluster, X)
    c = length(Set(cluster))
    FT = (T <: Integer) ? samesize_float[T] : T
    Y = zeros(FT, dims..., c)
    scatter_mean!(Y, FT.(X), cluster)
    Y
end

function _pooling_dim_check(cluster::Array{Int}, X::Array{T}) where {T<:Real}
    dim_c = size(cluster)
    d = length(dim_c)
    dim_X = size(X)
    n = length(dim_X) - d
    @assert n > 0 "X must have more dimensions than cluster."
    @assert dim_c == dim_X[n+1:end] "X must have the same latter dimension with cluster."
    dim_X[1:n]
end

pool(op::Symbol, cluster::Array, X::Array) = pool(Val(op), cluster, X)
pool(::Val{:add}, cluster::Array, X::Array) = sumpool(cluster, X)
pool(::Val{:sub}, cluster::Array, X::Array) = subpool(cluster, X)
pool(::Val{:mul}, cluster::Array, X::Array) = prodpool(cluster, X)
pool(::Val{:div}, cluster::Array, X::Array) = divpool(cluster, X)
pool(::Val{:max}, cluster::Array, X::Array) = maxpool(cluster, X)
pool(::Val{:min}, cluster::Array, X::Array) = minpool(cluster, X)
pool(::Val{:mean}, cluster::Array, X::Array) = meanpool(cluster, X)


text:
#=
This file is for updating the NCDataset database that stores
input values to the moist thermodynamic state constructors
which have caused convergence issues. Updating this database
allows us to optimize the convergence rate of the moist
thermodynamics constructor for a variety of realistic input
values.
=#

using NCDatasets
using DelimitedFiles

"""
    compress_targz(file)

Platform-independent file compression
"""
function compress_targz(file)
  if Sys.iswindows()
    error("Needs to be implemented")
  else
    run(`tar -zcvf $file $(readdir())`)
  end
end

const archive = true
const clean = true
const append = false
const create = true

# Error check
@assert !(append && create) # cannot create and append
@assert !(append && clean)  # cannot clean and append

folder = joinpath(@__DIR__,"MTConstructorData")
output_file = joinpath(@__DIR__,"MTConstructorDataZipped.tar.gz")
mkpath(folder)

constructors = Dict("PhaseEquil"=>(:e_int, :ρ, :q_tot),)
get_nc(k) = joinpath(folder,"test_data_$(k).nc")
get_data_to_append(k) = joinpath(@__DIR__,"test_data_$(k).csv")

if clean
  for k in keys(constructors)
    rm(get_nc(k); force=true)
  end
end
FT = Float64
if create || append
  for (k,v) in constructors
    if isfile(get_data_to_append(k))
      if append
        ds = Dataset(get_nc(k),"a")
      elseif create
        ds = Dataset(get_nc(k),"c")
      end
      data_all = readdlm(get_data_to_append(k), ',')
      for (i,_v) in enumerate(v)
        vardata = Array{FT}(data_all[2:end,i])
        s = string(_v)
        defVar(ds, s, vardata, ("datapoint",))
      end
      close(ds)
    end
  end
end

if archive
  cd(folder) do
    compress_targz(output_file)
  end
end


text:
function read_input(input_filename)
    lines = readlines(input_filename)
    grid = zeros(Int, length(lines), length(lines[1]))
    for iy in 1:length(lines)
        for ix in 1:length(lines[iy])
            grid[iy, ix] = parse(Int, lines[iy][ix])
        end
    end
    return grid
end

function is_local_minimum(x, y, grid)
    if x == 1
        local_min_x = (grid[x,y] < grid[x+1,y])
    elseif x == size(grid)[1]
        local_min_x = (grid[x,y] < grid[x-1,y])
    else
        local_min_x = (grid[x,y] < grid[x+1,y]) && (grid[x,y] < grid[x-1,y])
    end
    if y == 1
        local_min_y = (grid[x,y] < grid[x,y+1])
    elseif y == size(grid)[2]
        local_min_y = (grid[x,y] < grid[x,y-1])
    else
        local_min_y = (grid[x,y] < grid[x,y+1]) && (grid[x,y] < grid[x,y-1])
    end
    return (local_min_x && local_min_y)
end

function find_basin_coords(x, y, grid, basin_coords)
    # basin coords is a set
    # recursive function walks in all directions and adds coords to list
    # if : value at new grid square is not 9, AND new grid square isn't already in set,
    # AND new grid square is not off the ends of the grid.
    if grid[x,y] != 9
        push!(basin_coords, (x, y))
    end
    if ((x-1) in 1:size(grid)[1]) && (grid[x-1,y] != 9) && (!).(in((x-1,y),basin_coords))

        basin_coords = find_basin_coords(x-1, y, grid, basin_coords)
    end
    if ((x+1) in 1:size(grid)[1]) && (grid[x+1,y] != 9)  && (!).(in((x+1,y),basin_coords))
        basin_coords = find_basin_coords(x+1, y, grid, basin_coords)
    end
    if ((y-1) in 1:size(grid)[2]) && (grid[x,y-1] != 9)  && (!).(in((x,y-1),basin_coords))
        basin_coords = find_basin_coords(x, y-1, grid, basin_coords)
    end
    if ((y+1) in 1:size(grid)[2]) && (grid[x,y+1] != 9) && (!).(in((x,y+1),basin_coords))
        basin_coords = find_basin_coords(x, y+1, grid, basin_coords)
    end
    return basin_coords
end

function find_all_local_minima(grid)
    minima_coords = []
    for ix in 1:size(grid)[1]
        for iy in 1:size(grid)[2]
            if is_local_minimum(ix, iy, grid)
                push!(minima_coords,(ix,iy))
            end
        end
    end
    return minima_coords
end

function find_all_basin_sizes(grid)
    minima = find_all_local_minima(grid)
    basin_sizes = []
    for minimum in minima
        basin_coordinates = Set([])
        push!(basin_sizes, length(find_basin_coords(minimum[1], minimum[2],grid, basin_coordinates)))
    end
    return sort(basin_sizes, rev=true)

end

grid = read_input("input.txt")
basin_sizes = find_all_basin_sizes(grid)
answer = prod(basin_sizes[1:3])
println("answer is ",answer)


text:
#!/usr/bin/env julia

"Generates the index.md file from the README.md file"
function gen_index()
  whole_file = open(readstring, "README.md")

  header = """
  ---
  layout: default
  ---

  """

  open(f -> write(f, header * whole_file), "index.md", "w")
end

gen_index()
println("Generated index.md from README.md")

max_stars_repo_name: Saintis/netrunner-code-viewer

max_issues_repo_name: Saintis/netrunner-code-viewer

max_forks_repo_name: Saintis/netrunner-code-viewer


text:
S4s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1])
    Symmetry.Symel("C_2^1", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("S_4^1", [0.0 -1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_4^3", [0.0 1.0 0.0; -1.0 0.0 -0.0; 0.0 0.0 -1.0])]

S6s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1])
    Symmetry.Symel("i", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("C_3^1", [-0.4999999999999998 -0.8660254037844387 0.0; 0.8660254037844387 -0.4999999999999998 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("C_3^2", [-0.5000000000000003 0.8660254037844384 0.0; -0.8660254037844384 -0.5000000000000003 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("S_6^1", [0.5000000000000001 -0.8660254037844386 0.0; 0.8660254037844386 0.5000000000000001 0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_6^5", [0.49999999999999944 0.8660254037844392 0.0; -0.8660254037844392 0.49999999999999944 0.0; 0.0 0.0 -1.0])]

S8s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1])
    Symmetry.Symel("C_4^1", [0.0 -1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("C_2^1", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("C_4^3", [0.0 1.0 0.0; -1.0 0.0 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("S_8^1", [0.7071067811865476 -0.7071067811865475 0.0; 0.7071067811865475 0.7071067811865476 0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_8^3", [-0.7071067811865474 -0.7071067811865477 -0.0; 0.7071067811865477 -0.7071067811865474 0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_8^5", [-0.7071067811865479 0.7071067811865471 0.0; -0.7071067811865471 -0.7071067811865479 -0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_8^7", [0.707106781186547 0.707106781186548 0.0; -0.707106781186548 0.707106781186547 0.0; 0.0 0.0 -1.0])]

S10s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1])
    Symmetry.Symel("i", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("C_5^1", [0.30901699437494745 -0.9510565162951535 0.0; 0.9510565162951535 0.30901699437494745 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("C_5^2", [-0.8090169943749473 -0.5877852522924731 0.0; 0.5877852522924731 -0.8090169943749473 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("C_5^3", [-0.8090169943749475 0.587785252292473 0.0; -0.587785252292473 -0.8090169943749475 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("C_5^4", [0.30901699437494734 0.9510565162951535 0.0; -0.9510565162951535 0.30901699437494734 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("S_10^1", [0.8090169943749475 -0.5877852522924731 0.0; 0.5877852522924731 0.8090169943749475 0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_10^3", [-0.3090169943749474 -0.9510565162951536 -0.0; 0.9510565162951536 -0.3090169943749474 0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_10^7", [-0.30901699437494756 0.9510565162951538 0.0; -0.9510565162951538 -0.30901699437494756 -0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_10^9", [0.8090169943749477 0.5877852522924735 0.0; -0.5877852522924735 0.8090169943749477 0.0; 0.0 0.0 -1.0])]

"""S12s = [Symmetry.Symel("E", [1 0 0; 0 1 0; 0 0 1])
    Symmetry.Symel("C_6^1", [0.5000000000000001 -0.8660254037844386 0.0; 0.8660254037844386 0.5000000000000001 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("C_3^1", [-0.4999999999999998 -0.8660254037844388 0.0; 0.8660254037844388 -0.4999999999999998 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("C_2^1", [-1.0 0.0 0.0; 0.0 -1.0 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("C_3^2", [-0.5000000000000006 0.8660254037844385 0.0; -0.8660254037844385 -0.5000000000000006 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("C_6^5", [0.49999999999999944 0.8660254037844392 0.0; -0.8660254037844392 0.49999999999999944 0.0; 0.0 0.0 1.0])
    Symmetry.Symel("S_12^1", [0.8660254037844387 -0.49999999999999994 0.0; 0.49999999999999994 0.8660254037844387 0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_4^1", [0.0 -1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_12^5", [-0.8660254037844385 -0.5000000000000006 -0.0; 0.5000000000000006 -0.8660254037844385 0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_12^7", [-0.866025403784439 0.49999999999999944 0.0; -0.49999999999999944 -0.866025403784439 -0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_4^3", [0.0 1.0 0.0; -1.0 0.0 -0.0; 0.0 0.0 -1.0])
    Symmetry.Symel("S_12^11", [0.8660254037844383 0.500000000000001 0.0; -0.500000000000001 0.8660254037844383 0.0; 0.0 0.0 -1.0])]
"""
S4cn = ["A", "B", "E"]
S4ct = C4ct
S6cn = ["Ag", "Eg", "Au", "Eu"]
S6ct = C3hct
S8cn = ["A", "B", "E1", "E2", "E3"]
S8ct = [1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0;
        1.0 -1.0  1.0 -1.0  1.0 -1.0  1.0 -1.0;
        2.0  rt2  0.0 -rt2 -2.0 -rt2  0.0  rt2;
        2.0  0.0 -2.0  0.0  2.0  0.0 -2.0  0.0;
        2.0 -rt2  0.0  rt2 -2.0  rt2  0.0 -rt2]
S10cn = ["Ag", "E1g", "E2g", "Au", "E1u", "E2u"]
S10ct = C5hct


text:
module PENBBControl

	using Sockets
        using HTTP, JSON
        #using PENScintAnalysis

	include("HV_control.jl")
	include("motor_base.jl")
	include("motor_commands.jl")
	#include("PENBBControl.jl")

end


text:
# Examples presented in class - Lecture 1
# written by Fabricio - 07/09/2018

using JuMP, Cbc
# JuMP is for implementing math. programming models;
# Cbc is for solving them.

# Example 1 - resource allocation

# Problem data
i = 1:2 # i=1: Seattle; i=2: San Diego
j = 1:3 # j=1: New York; j=2: Chicago; j=3: Miami

C = [350 600] # Capacities of the factories
D = [325 300 275] # Demand of clients
T = [2.5 1.7 1.8
     3.5 1.9 1.4] # Transportation costs

# Model implementation
# Creates a model and informs the solver to be used
m = Model(solver = CbcSolver())

# Decision variable for the total transported
@variable(m, x[i,j] >= 0)

# Capacity constraint
@constraint(m, cap[i = 1:2], sum(x[i,j] for j = 1:3) <= C[i])
# Demand constraint
@constraint(m, dem[j = 1:3], sum(x[i,j] for i = 1:2) >= D[j])

# Total distribution cost that we want to minimise
@objective(m, Min, sum(T[i,j]*x[i,j] for i=1:2, j=1:3))

println(m) # Prints the mathematical model for debugging
solve(m) # Solve the model

# Prints the optimal solution
println("\nDistribution plan: \n", getvalue(x))

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

# Example 2 - portfolio optimisation

using Ipopt, Plots
# Ipopt is for solving nonlinear problems - different technology
# Plots is for plotting (!)

# Read the daily prices data from a .csv file.
data = readcsv("prices.csv")
#First row has the names of the stocks
stock_names = data[1, 1:end-2]
#Last two columns has data and US$ rate
prices_data = data[2:end, 1:end-2]

#Having a peek at the data

cs = plot(prices_data, label = stock_names)
plotlyjs()

gui(cs)
#Returns are calculated as (p(t+1) - p(t))/p(t-1)
returns_data = diff(prices_data) ./ prices_data[1:end-1,:]

#Number of days and stocks in data
T, n =  size(returns_data)
#Calculates expected return and covariance
μ, Σ =  mean(returns_data, 1), cov(returns_data)
# Input form the model: minimum average return required.
r_min = 0.2

port = Model(solver=IpoptSolver())

#allocation variables
@variable(port, 0 <= x[1:n] <= 1);

#notice the division by T to correct the average
@constraint(port, sum(μ[j]*x[j] for j=1:n) >= r_min/T);

@objective(port, Min , sum(x[i]*Σ[i,j]*x[j] for i=1:n,j=1:n));

solve(port);

alloc = getvalue(x)
# Should be the same as r_min
ret = μ*alloc*T
# Looking at risk as the st. deviation of the returns.
risk = sqrt(alloc'*Σ*alloc*T)

# Organising the data for plotting. In an array, they are
# considered different series (with different colours)
stock_names = convert(Array{String}, stock_names)
alloc = convert(Array{Float64}, alloc)
# plots a pie chart
pie(stock_names, alloc, legend = false)

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

# Example 3 - Robust Knapsack Problem

using Distributions, ECOS, LaTeXStrings
# Distributions includes probability distributions,
# ECOS is a solver that we need because of the stucture of the problem
# LatexStrings is to write LaTeX in our plots

N = 18

#Input data
value = [50, 40, 70, 55, 80, 35, 65, 50, 60, 85, 20, 45, 55, 25, 80, 45, 45, 65]
weight_average = [3.5, 4, 5.5, 5, 6, 4.5, 6, 4, 5.5, 7, 5, 4.5, 7, 3.5, 5.5, 3.5, 4, 6.5]
weight_stdev = weight_average*0.3

#Randomly generating initial data using normal distribution
weight_data = Array{Float64}(100,N)
for j= 1:N
   weight_data[:,j] = rand(Normal(weight_average[j], weight_stdev[j]), 100)
end

#Input data: weight limit (capacity)
capacity = 20

# Protection elipsoid: adds 50% of weight as a protection level
P = 0.5*weight_average

# Declaring the model and solving it as a function.
function solve_robust_model(N, value, weight_average, P, Γ)
    m = Model(solver = ECOSSolver()) #creates the model, select the solver

    @variable(m, 0 <= x[1:N] <= 1 ) # creates the binary variables x, one for each item

    @constraint(m, sum(weight_average[j]*x[j] for j = 1:N) + Γ*norm(P'*x) <= capacity) # declare the knapsack constraint

    @objective(m, Max, sum(value[j]*x[j] for j = 1:N)) # declare the objective function

    solve(m)

    return getvalue(x)
end

# This function simulates the item selection againts feasibility.
function feasibility_estimate(solution, repetitions)
    feasible_count = 0
    actual_weight = zeros(18)
    for n= 1:repetitions
        for j= 1:N
           # generate random weights according to distribution
           actual_weight[j] = rand(Normal(weight_average[j], weight_stdev[j]))
        end
        #if total weight more than capacity => problem infeasible.
        if actual_weight'*solution <= capacity
            feasible_count += 1
        end
    end
    return feasible_count/repetitions
end

# generate a range of Gammas to try
Γ_range = linspace(0,1,10)

# storing the results of each run
feas = []
total_value = []
Γ_used = []

# for each Gamma, solve te model and store the results
for Γ in Γ_range
    println(Γ)
    x = solve_robust_model(N, value, weight_average, P, Γ)
    push!(Γ_used, Γ)
    push!(feas, feasibility_estimate(x,5000))
    push!(total_value, value'x)
end

# potting the results from simulation
p1 = plot(Γ_used, feas, xlabel = L"\Gamma", ylabel = "feas. prob.", legend=false)
p2 = plot(Γ_used, total_value, xlabel = L"\Gamma", ylabel = "total value", legend=false, color=:orange)
plot(p1,p2, layout = (2,1))

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

# Example 4 - Classification

# Generating random multivariate data.
data11 = rand(Normal(8, 4), 100)
data12 = rand(Normal(12, 2), 100)
data1 = [data11 data12]

data01 = rand(Normal(2, 6), 100)
data02 = rand(Normal(5, 2), 100)
data0 = [data01 data02]

scatter(data1[:,1], data1[:,2], label="Pos. obs.")
scatter!(data0[:,1], data0[:,2], label="Neg. obs.")

function classify(data1, data0, δ)
    m = Model(solver = IpoptSolver())

    @variable(m, u[1:100] >= 0)
    @variable(m, v[1:100] >= 0)
    @variable(m, a[1:2])
    @variable(m, b)

    @constraint(m, error1[i=1:100], sum(a[j]*data1[i,j] for j = 1:2) - b - u[i] <= -1)
    @constraint(m, error0[i=1:100], sum(a[j]*data0[i,j] for j = 1:2) - b + v[i] >= 1)

# This parameter trades off accuracy and robustness for the classifier
    @objective(m, Min, sum(u[i] for i=1:100) + sum(v[i] for i=1:100)
        + δ*(sum(a[j]^2 for j=1:2)))

    solve(m);

    return getvalue(a), getvalue(b), getvalue(u), getvalue(v)
end

a, b, u1, v1 = classify(data1, data0, 0.01);

# Calculating the points on plane defined by a and b
x1 = linspace(-10,20,100);
x2 = (b - a[1]*x1)/ a[2];
# boundaries of the slab
x2u = (b + 1 - a[1]*x1)/ a[2];
x2d = (b - 1 - a[1]*x1)/ a[2];
plot!(x1, [x2u x2 x2d],
    xlabel = "x_1",
    ylabel = "x_2",
    label= [ "", "Classifier", ""] ,
    color=[:green :green :green],
    linestyle=[ :dot :solid :dot],
    )

# another classifier, with δ=50
a, b, u2, v2 = classify(data1, data0, 100);
x2 = (b - a[1]*x1)/ a[2];
x2u = (b + 1 - a[1]*x1)/ a[2];
x2d = (b - 1 - a[1]*x1)/ a[2];
plot!(x1, [x2u x2 x2d],
    xlabel = "x_1",
    ylabel = "x_2",
    legend = false,
    color=[:purple :purple :purple],
    linestyle=[:dot :solid :dot],
    )


text:
using RbfInterpolationTools
using Test

X = hcat(rand(100), rand(100))

@testset "Basic Propieties" begin
    z_obs = [1 for (x, y) in zip(rand(100), rand(100))]
    inter = RbfInterpolationTools.Interpolator(X = X, y = z_obs, basis_name = "InvMQ", ϵ = 0.01)
    new_domain = (hcat(rand(1000), rand(1000)))
    @test isapprox(RbfInterpolationTools.interpolatePoint(inter, [rand(), rand()]), 1.0; atol = 0.1)
    inter = RbfInterpolationTools.Interpolator(X = X, y = z_obs .* 0, basis_name = "Cubic")
    app_z = RbfInterpolationTools.interpolateMesh(inter, new_domain)
    for i = 1:1000
        @test isapprox(app_z[i], 0.0; atol = 0.1)
    end
    z_obs = [x for x in 0:1/999:1]
    inter = RbfInterpolationTools.Interpolator(X = hcat(0:1/999:1), y = z_obs, basis_name = "ThinPlate")
    @test isapprox(RbfInterpolationTools.gradientPoint(inter, [0.5])[1], 1.0; atol = 0.1)

end


text:
using DataAssimilation
using DifferentialEquations
using DiffEqSensitivity
using LinearAlgebra
using Zygote
using Statistics
using Distributions
using Plots
using Random

# set random seed for reproducible results
Random.seed!(1)

## ----- setup true solution ------------------------------------

# define function for RHS of lorenz system
function lorenz!(du, u, p, t)
    x,y,z = u
    σ,ρ,β = p

    du[1] = σ*(y - x)
    du[2] = x*(ρ - z) - y
    du[3] = x*y - β*z
end

θ = [10.0, 28.0, (8/3)]  # true parameter values
dt = 0.01  # time step for simulation (this may be unnecessary)
tend = 10.0  # max time for total simulation
ts_sim = (0.0, tend)
u0True = [1.0, 1.0, 1.0]  # true initial condition

# integrate to obtain the true solution
prob_true = ODEProblem(lorenz!, u0True, ts_sim, θ, saveat=dt)
sol_true = solve(prob_true)

# set up observation uncertainties
σ_m = 0.15  # i.e. observation std
# set up observation error covariance matrix
R = σ_m^2 .* I(3)

# set up details for measurement points
dt_m = 0.2
tm_end = 2.0  # final observation time
ts_m = dt_m:dt_m:tm_end

# define observation operator
function h(u)
    w = u
    return w
end

# define observation operator jacobian. We should probably generate this symbolically
function Dh(u)
    n = size(u, 1)
    D = I(n)
    return D
end

# populate observation vector
w = hcat([h(sol_true(tᵢ))+rand(MvNormal(zeros(3), R)) for tᵢ ∈ ts_m]...)

# set up data assimilation
u0b = [2.0; 3.0; 4.0]  # guess at initial condition
σ_b = 0.1  # model error std
B = σ_b^2 .* I(3)  # initial model error covariance
Q = 0.0 * I(3)  # process noise (assumed zero)

# integrate to obtain the true solution
prob_no_analysis = ODEProblem(lorenz!, u0b, ts_sim, θ, saveat=dt)
sol_no_analysis = solve(prob_no_analysis)

ekf_prob = AssimilationProblem(u0b,
                               ts_sim,
                               lorenz!,
                               θ,
                               dt,
                               Q,
                               B,
                               R,
                               dt_m,
                               h,
                               Dh,
                               w
                               )
uₐ, P = EKFsolve(ekf_prob)

# var_names = ["x(t)", "y(t)", "z(t)"]
# plots = []
# for i ∈ 1:3
#     p = plot(sol_true, vars=(0,i), color=:black, alpha=0.75, linewidth=2, label="true system")
#     plot!(p, ts_m, w[i, :], seriestype=:scatter, color="light green", label="observation" )
#     plot!(p, sol_no_analysis, vars=(0, i), linestyle=:dashdot, color="purple", label="forecast")
#     plot!(p, 0:dt:tend, ua[i, :], color="royal blue", alpha=1, linewidth = 1, label="analysis")
#     xlabel!(p, "t")
#     ylabel!(p, var_names[i])
#     push!(plots, p)
# end

# final_plot = plot(plots..., layout=(3,1))

# savefig("EKF_using_zygote.svg")
# savefig("EKF_using_zygote.png")
# display(final_plot)


text:
using Downloads
using Pkg.Artifacts
using ZipFile
""" unzip(file,exdir)
    From https://discourse.julialang.org/t/how-to-extract-a-file-in-a-zip-archive-without-using-os-specific-tools/34585/5?u=hellemo
"""
function unzip(file, exdir="")
    fileFullPath = isabspath(file) ? file : joinpath(pwd(), file)
    basePath = dirname(fileFullPath)
    outPath = (exdir == "" ? basePath : (isabspath(exdir) ? exdir : joinpath(pwd(), exdir)))
    isdir(outPath) ? "" : mkdir(outPath)
    zarchive = ZipFile.Reader(fileFullPath)
    for f in zarchive.files
        fullFilePath = joinpath(outPath, f.name)
        if (endswith(f.name, "/") || endswith(f.name, "\\"))
            mkdir(fullFilePath)
        else
            write(fullFilePath, read(f))
        end
    end
    return close(zarchive)
end

artifact_toml = realpath(joinpath(@__DIR__, "..", "Artifacts.toml"))
chmod(artifact_toml, 0o644)
ucanaccess_hash = artifact_hash("ucanaccess", artifact_toml)

# Download and unzip to Artifact if not already existing
if ucanaccess_hash === nothing || !artifact_exists(ucanaccess_hash)
    ucanaccess_hash = create_artifact() do artifact_dir
        tmp_file = Downloads.download("https://sourceforge.net/projects/ucanaccess/files/UCanAccess-5.0.1.bin.zip/download")
        unzip(tmp_file, artifact_dir) # Unpack to artifact_dir
    end
    bind_artifact!(artifact_toml, "ucanaccess", ucanaccess_hash; force=true)
end

max_stars_repo_name: hellemo/Access

max_issues_repo_name: hellemo/Access

max_forks_repo_name: hellemo/Access


text:
# julia Mapping.jl --parallel no --lmax 5
##################################################
# Parameters of the grid in (a,j)--space
# This is the grid where the DRR & DNR coefficients are computed
# Evalue a between aminMeasure (mpc) and amaxMeasure (mpc)
##################################################
using HDF5 # To have access to .hf5 files
########################################
include("../sources/Main.jl") # Loading the main code
########################################
jminMeasure, jmaxMeasure = 0.001,1.0 # Range in j where the Djj are computed
aminMeasure, amaxMeasure = 0.1,1000.0 # Range in j where the Djj are computed
nbjMeasure = 30 # Number of j for which the Djj are computed
nbaMeasure = 30 # Number of a for which the Djj are computed
nbajGrid = nbjMeasure*nbaMeasure # Number of (a,j) for which the Djj are computed
tabjMeasure = exp.(range(log(jminMeasure),length=nbjMeasure,log(jmaxMeasure)))
tabaMeasure = exp.(range(log(aminMeasure),length=nbaMeasure,log(amaxMeasure)))

const tabajGrid  = zeros(Float64,2,nbajGrid) # Location (a,j) of the grid points where the diffusion coefficients are computed
const tabDRRjjGrid = zeros(Float64,nbajGrid) # Values of the DRR_jj coefficients on the (a,j)-grid
const tabDNRjjGrid = zeros(Float64,nbajGrid) # Values of the DNR_jj coefficients on the (a,j)-grid
const tabDjjGrid = zeros(Float64,nbajGrid) # Values of the total D_jj coefficients on the (a,j)-grid
const tabDNRJJRedGrid = zeros(Float64,nbajGrid) # Values of the DNR_jj coefficient$

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

function tabajGrid!()
    index = 1
    for ia=1:nbaMeasure
    aMeasure = tabaMeasure[ia]
        for ij=1:nbjMeasure
            jMeasure = tabjMeasure[ij]
            tabajGrid[1,index], tabajGrid[2,index] = aMeasure, jMeasure
            index += 1
        end
    end
end

function init_tabDRRGrid!()
    for iGrid=1:nbajGrid
        tabDRRjjGrid[iGrid] = 0.0
    end
end

function init_tabDNRGrid!()
    for iGrid=1:nbajGrid
        tabDNRjjGrid[iGrid] = 0.0

        tabDNRJJRedGrid[iGrid] = 0.0
    end
end

function init_tabDGrid!()
    for iGrid=1:nbajGrid
        tabDjjGrid[iGrid] = 0.0
    end
end

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

function tabDRRGrid!()
    init_tabDRRGrid!() # Making sure that the grid is initially set to 0
    if (PARALLEL == "yes") # Computation is made in parallel
        Threads.@threads for iGrid=1:nbajGrid # Loop over the elements of the (a,j)-grid
            a, j = tabajGrid[1,iGrid], tabajGrid[2,iGrid] # Current (a,j) location
            if j<=jlc(a) # Do not compute DRR within the loss cone
                tabDRRjjGrid[iGrid] = 0.0
            elseif j<1.0
                tabSMARes_parallel = zeros(Float64,2,2) # Container of tabSMARes for the current thread
                IntTable_K_parallel = IntTable_create!()

                DRR = DRR_jj(a,j,IntTable_K_parallel,tabSMARes_parallel)
                tabDRRjjGrid[iGrid] = DRR # Computing DRR_jj with the thread's containers
            else
                tabDRRjjGrid[iGrid] = 0.0
            end
        end
    else # Computation is not made in parallel
        for iGrid=1:nbajGrid # Loop over the elements of the (a,j)-grid
            a, j = tabajGrid[1,iGrid], tabajGrid[2,iGrid] # Current (a,j) location
            if j<=jlc(a) # Do not compute DRR within the loss cone
                tabDRRjjGrid[iGrid] = 0.0
            elseif j<1.0
                DRR = DRR_jj(a,j)
                tabDRRjjGrid[iGrid] = DRR # Computing DRR_jj with the global containers
            else
                tabDRRjjGrid[iGrid] = 0.0
            end
        end
    end
end

function tabDNRGrid!()
    init_tabDNRGrid!() # Making sure that the grid is initially set to 0
    if (PARALLEL == "yes") # Computation is made in parallel
        Threads.@threads for iGrid=1:nbajGrid # Loop over the elements of the (a,j)-grid
            a, j = tabajGrid[1,iGrid], tabajGrid[2,iGrid] # Current (a,j) location
#            if j<=jlc(a) # Do not compute DNR within the loss cone
 #               tabDNRjjGrid[iGrid] = 0.0

  #              tabDNRJJRedGrid[iGrid] = 0.0
   #         elseif j<1.0
            if j<1.0
                DNR = DNR_jj(a,j)
                tabDNRjjGrid[iGrid] = DNR # Computing DNR_jj with the thread's containers

                tabDNRJJRedGrid[iGrid] = DNR_JJ_red(a,j)
            else
                tabDNRjjGrid[iGrid] = 0.0
                tabDNRJJRedGrid[iGrid] = 0.0
            end
        end
    else # Computation is not made in parallel
        for iGrid=1:nbajGrid # Loop over the elements of the (a,j)-grid
            a, j = tabajGrid[1,iGrid], tabajGrid[2,iGrid] # Current (a,j) location
 #            if j<=jlc(a) # Do not compute DNR within the loss cone
 #               tabDNRjjGrid[iGrid] = 0.0

  #              tabDNRJJRedGrid[iGrid] = 0.0
   #         elseif j<1.0
            if j<1.0
                DNR = DNR_jj(a,j)
                tabDNRjjGrid[iGrid] = DNR # Computing DNR_jj with the global containers

                tabDNRJJRedGrid[iGrid] = DNR_JJ_red(a,j)
            else
                tabDNRjjGrid[iGrid] = 0.0
                tabDNRJJRedGrid[iGrid] = 0.0
            end
        end
    end
end

function tabDGrid!()
    init_tabDGrid!() # Making sure that the grid is initially set to 0
    if (PARALLEL == "yes") # Computation is made in parallel
        Threads.@threads for iGrid=1:nbajGrid # Loop over the elements of the (a,j)-grid
            a, j = tabajGrid[1,iGrid], tabajGrid[2,iGrid] # Current (a,j) location
  #          if j<=jlc(a) # Do not compute D_jj within the loss cone
   #             tabDjjGrid[iGrid] = 0.0
    #        elseif j<1.0
            if j<1.0
                DNR = tabDNRjjGrid[iGrid]
                DRR = tabDRRjjGrid[iGrid]
                tabDjjGrid[iGrid] = DNR + DRR # Computing D_jj with the thread's containers
            else
                tabDjjGrid[iGrid] = 0.0
            end
        end
    else # Computation is not made in parallel
        for iGrid=1:nbajGrid # Loop over the elements of the (a,j)-grid
            a, j = tabajGrid[1,iGrid], tabajGrid[2,iGrid] # Current (a,j) location
        #    if j<=jlc(a) # Do not compute D_jj within the loss cone
       #         tabDjjGrid[iGrid] = 0.0
       #     elseif j<1.0
            if j<1.0
                DNR = tabDNRjjGrid[iGrid]
                DRR = tabDRRjjGrid[iGrid]
                tabDjjGrid[iGrid] = DNR + DRR # Computing D_jj with the thread's containers
            else
                tabDjjGrid[iGrid] = 0.0
            end
        end
    end
end

########################################
namefile = "../data/Dump_Diffusion_Coefficients.hf5"
########################################
# Function that writes to .hf5 files
########################################
function writedump!(namefile)
    file = h5open(namefile,"w") # Opening the file
    write(file,"tabaj",tabajGrid) # Dumping the grid of (a,j)
    write(file,"tabj",tabjMeasure) # Dumping the grid of j
    write(file,"taba",tabaMeasure) # Dumping the grid of a
    write(file,"tabDRRjj",tabDRRjjGrid) # Dumping the DRRjj(a,j) for the current used value of lmax
    write(file,"tabDNRjj",tabDNRjjGrid) # Dumping the DNRjj(a,j) for the current used value of lmax
    write(file,"tabDjj",tabDjjGrid) # Dumping the total Djj(a,j) for the current used value of lmax
    write(file,"amin",aminMeasure) # Dumping aminMeasure
    write(file,"amax",amaxMeasure) # Dumping amaxMeasure
    write(file,"jmin",jminMeasure) # Dumping jminMeasure
    write(file,"jmax",jmaxMeasure) # Dumping jmaxMeasure
    write(file,"nbj",nbjMeasure) # Dumping nbjMeasure
    write(file,"nba",nbaMeasure) # Dumping nbaMeasure

    write(file,"tabDNRJJRed",tabDNRJJRedGrid)

    write(file,"G",G)
    write(file,"cvel",cvel)
    write(file,"mBH",mBH)

    close(file) # Closing the file
end

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

tabajGrid!()

println("Computing SRR diffusion coefficients...")
@time tabDRRGrid!()
println("Computing NR diffusion coefficients...")
@time tabDNRGrid!()
println("Computing total diffusion coefficients...")
@time tabDGrid!()

########################################
println("Saving data...")
writedump!(namefile) # Dumping the computed table


text:
using SupervisedLearning

X = rand(10, 10000000)
t = round(rand(10000000) * 4) + 1
ce = MultivalueClassEncoding(["V1","V2","V3","V4","V5"])
ce2 = OneOfKClassEncoding(["V1","V2","V3","V4","V5"])
t2 = labelencode(ce2, labeldecode(ce, t))

msg("EncodedInMemoryLabeledDataSource: Instantiate from array (10x10000000) with targets (1x10000000)")
ds = EncodedInMemoryLabeledDataSource(X, t, ce)
@spc_time ds = EncodedInMemoryLabeledDataSource(X, t, ce)

msg("EncodedInMemoryLabeledDataSource: Access all features")
@spc_time tmp = features(ds)

msg("EncodedInMemoryLabeledDataSource: Access features minibatch n=1000")
@spc_time tmp = features(ds, 10000, 1000)

msg("EncodedInMemoryLabeledDataSource: Access features minibatch n=10000")
@spc_time tmp = features(ds, 10000, 10000)

msg("EncodedInMemoryLabeledDataSource: Access all targets")
@spc_time tmp = targets(ds)

msg("EncodedInMemoryLabeledDataSource: Access targets minibatch n=1000")
@spc_time tmp = targets(ds, 10000, 1000)

msg("EncodedInMemoryLabeledDataSource: Access targets minibatch n=10000")
@spc_time tmp = targets(ds, 10000, 10000)

msg("EncodedInMemoryLabeledDataSource: splitTrainTest")
@spc_time trainDs, testDs = splitTrainTest!(ds, p_train = .7)

# ===================================================================

msg("EncodedInMemoryLabeledDataSource: Instantiate from array (10x10000000) with targets (5x10000000)")
ds = EncodedInMemoryLabeledDataSource(X, t2, ce2)
@spc_time ds = EncodedInMemoryLabeledDataSource(X, t2, ce2)

msg("EncodedInMemoryLabeledDataSource: Access all features")
@spc_time tmp = features(ds)

msg("EncodedInMemoryLabeledDataSource: Access minibatch n=1000")
@spc_time tmp = features(ds, 10000, 1000)

msg("EncodedInMemoryLabeledDataSource: Access minibatch n=10000")
@spc_time tmp = features(ds, 10000, 10000)

msg("EncodedInMemoryLabeledDataSource: Access all targets")
@spc_time tmp = targets(ds)

msg("EncodedInMemoryLabeledDataSource: Access targets minibatch n=1000")
@spc_time tmp = targets(ds, 10000, 1000)

msg("EncodedInMemoryLabeledDataSource: Access targets minibatch n=10000")
@spc_time tmp = targets(ds, 10000, 10000)

msg("EncodedInMemoryLabeledDataSource: splitTrainTest")
@spc_time trainDs, testDs = splitTrainTest!(ds, p_train = .7)


text:
using MPI
using ClimateMachine
using Logging
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGmethods
using ClimateMachine.DGmethods.NumericalFluxes
using ClimateMachine.MPIStateArrays
using ClimateMachine.LowStorageRungeKuttaMethod
using LinearAlgebra
using ClimateMachine.GenericCallbacks:
    EveryXWallTimeSeconds, EveryXSimulationSteps
using ClimateMachine.ODESolvers

MPI.Initialized() || MPI.Init()
mpicomm = MPI.COMM_WORLD

# set up domain
topl = StackedBrickTopology(
    mpicomm,
    (0:10, 0:10, 0:3);
    periodicity = (false, false, false),
    boundary = ((1, 1), (1, 2), (1, 2)),
)

@show MPI.Comm_rank(mpicomm) length(topl.realelems)


text:
function cartesian_iter(state, res, A, Asize)
    for i in CartesianIndices(Asize)
        res[i] = A[i]
    end
    return
end

function clmap!(ctx, f, out, b)
    i = linear_index(ctx) # get the kernel index it gets scheduled on
    out[i] = f(b[i])
    return
end

function ntuple_test(ctx, result, ::Val{N}) where N
    result[1] = ntuple(Val(N)) do i
        Float32(i) * 77f0
    end
    return
end

function ntuple_closure(ctx, result, ::Val{N}, testval) where N
    result[1] = ntuple(Val(N)) do i
        Float32(i) * testval
    end
    return
end

@testsuite "base" (AT, eltypes)->begin
    @testset "copyto!" begin
        x = fill(0f0, (10, 10))
        y = rand(Float32, (20, 10))
        a = AT(x)
        b = AT(y)
        r1 = CartesianIndices((1:7, 3:8))
        r2 = CartesianIndices((4:10, 3:8))
        copyto!(x, r1, y, r2)
        copyto!(a, r1, b, r2)
        @test x == Array(a)
        r2 = CartesianIndices((4:11, 3:8))
        @test_throws ArgumentError copyto!(a, r1, b, r2)

        r2 = CartesianIndices((4:10, 3:8))
        x2 = fill(0f0, (10, 10))
        copyto!(x2, r1, b, r2)
        @test x2 == x

        fill!(a, 0f0)
        copyto!(a, r1, y, r2)
        @test Array(a) == x

        x = fill(0f0, (10,))
        y = rand(Float32, (20,))
        a = AT(x)
        b = AT(y)
        r1 = CartesianIndices((1:7,))
        r2 = CartesianIndices((4:10,))
        copyto!(x, r1, y, r2)
        copyto!(a, r1, b, r2)
        @test x == Array(a)
        r2 = CartesianIndices((4:11,))
        @test_throws ArgumentError copyto!(a, r1, b, r2)

        x = fill(0f0, (10,))
        y = rand(Float32, (20,))
        a = AT(x)
        b = AT(y)
        r1 = 1:7
        r2 = 4:10
        copyto!(x, r1, y, r2)
        copyto!(a, r1, b, r2)
        @test x == Array(a)

        x = fill(0f0, (10,))
        y = fill(1f0, (10,))
        a = AT(x)
        b = AT(y)
        copyto!(a, b)
        @test Float32.(y) == Array(a)

        # wrapped gpu array to wrapped gpu array
        x = rand(Float32, 4, 4)
        a = AT(x)
        b = view(a, 2:3, 2:3)
        c = AT{eltype(b)}(undef, size(b))
        copyto!(c, b)
        Array(c) == Array(b)

        # wrapped gpu array to cpu array
        z = Array{eltype(b)}(undef, size(b))
        copyto!(z, b)
        @test z == Array(b)

        # cpu array to wrapped gpu array
        copyto!(b, z)

        # bug in copyto!
        ## needless N type parameter
        @test compare((x,y)->copyto!(y, selectdim(x, 2, 1)), AT, ones(Float32, 2, 2, 2), zeros(Float32, 2, 2))
        ## inability to copyto! smaller destination
        ## (this was broken on Julia <1.5)
        @test compare((x,y)->copyto!(y, selectdim(x, 2, 1)), AT, ones(Float32, 2, 2, 2), zeros(Float32, 3, 3))

        if (Float32 in eltypes && Float64 in eltypes)
            # mismatched types
            let src = rand(Float32, 4)
                dst = AT{Float64}(undef, size(src))
                copyto!(dst, src)
                @test Array(dst) == src
            end
            let dst = Array{Float64}(undef, 4)
                src = AT(rand(Float32, size(dst)))
                copyto!(dst, src)
                @test Array(src) == dst
            end
        end
    end

    @testset "cat" begin
        @test compare(hcat, AT, rand(Float32, 3), rand(Float32, 3))
        @test compare(hcat, AT, rand(Float32, ), rand(Float32, 1, 3))
        @test compare(hcat, AT, rand(Float32, 1, 3), rand(Float32))
        @test compare(hcat, AT, rand(Float32, 3), rand(Float32, 3, 3))
        @test compare(hcat, AT, rand(Float32, 3, 3), rand(Float32, 3))
        @test compare(hcat, AT, rand(Float32, 3, 3), rand(Float32, 3, 3))
        #@test compare(hcat, AT, rand(Float32, ), rand(Float32, 3, 3))
        #@test compare(hcat, AT, rand(Float32, 3, 3), rand(Float32))

        @test compare(vcat, AT, rand(Float32, 3), rand(Float32, 3))
        @test compare(vcat, AT, rand(Float32, 3, 3), rand(Float32, 3, 3))
        @test compare(vcat, AT, rand(Float32, ), rand(Float32, 3))
        @test compare(vcat, AT, rand(Float32, 3), rand(Float32))
        @test compare(vcat, AT, rand(Float32, ), rand(Float32, 3, 3))
        #@test compare(vcat, AT, rand(Float32, 3, 3), rand(Float32))

        @test compare((a,b) -> cat(a, b; dims=4), AT, rand(Float32, 3, 4), rand(Float32, 3, 4))
    end

    @testset "reshape" begin
        @test compare(reshape, AT, rand(Float32, 10), Ref((10,)))
        @test compare(reshape, AT, rand(Float32, 10), Ref((10,1)))
        @test compare(reshape, AT, rand(Float32, 10), Ref((1,10)))

        @test_throws Exception reshape(AT(rand(Float32, 10)), (10,2))
    end

    @testset "reinterpret" begin
        a = rand(ComplexF32, 22)
        A = AT(a)
        af0 = reinterpret(Float32, a)
        Af0 = reinterpret(Float32, A)
        @test Array(Af0) == af0
        a = rand(ComplexF32, 4, 4)
        A = AT(a)

        a = rand(ComplexF32, 10 * 10)
        A = AT(a)
        af0 = reshape(reinterpret(Float32, vec(a)), (20, 10))
        Af0 = reshape(reinterpret(Float32, vec(A)), (20, 10))
        @test Array(Af0) == af0
    end

    AT <: AbstractGPUArray && @testset "ntuple test" begin
        result = AT(Vector{NTuple{3, Float32}}(undef, 1))
        gpu_call(ntuple_test, result, Val(3))
        @test Array(result)[1] == (77, 2*77, 3*77)
        x = 88f0
        gpu_call(ntuple_closure, result, Val(3), x)
        @test Array(result)[1] == (x, 2*x, 3*x)
    end

    AT <: AbstractGPUArray && @testset "cartesian iteration" begin
        Ac = rand(Float32, 32, 32)
        A = AT(Ac)
        result = fill!(copy(A), 0.0f0)
        gpu_call(cartesian_iter, result, A, size(A))
        Array(result) == Ac
    end

    AT <: AbstractGPUArray && @testset "Custom kernel from Julia function" begin
        x = AT(rand(Float32, 100))
        y = AT(rand(Float32, 100))
        gpu_call(clmap!, -, x, y; target=x)
        jy = Array(y)
        @test map!(-, jy, jy) ≈ Array(x)
    end

    @testset "map" begin
        @test compare((a, b)-> map(+, a, b),    AT, rand(Float32, 10), rand(Float32, 10))
        @test compare((a, b)-> map!(-, a, b),   AT, rand(Float32, 10), rand(Float32, 10))
        @test compare((a, b, c, d)-> map!(*, a, b, c, d), AT, rand(Float32, 10), rand(Float32, 10), rand(Float32, 10), rand(Float32, 10))
    end

    @testset "repeat" begin
        @test compare(a-> repeat(a, 5, 6),  AT, rand(Float32, 10))
        @test compare(a-> repeat(a, 5),     AT, rand(Float32, 10))
        @test compare(a-> repeat(a, 5),     AT, rand(Float32, 5, 4))
        @test compare(a-> repeat(a, 4, 3),  AT, rand(Float32, 10, 15))
        @test compare(a-> repeat(a, 0),     AT, rand(Float32, 10))
        @test compare(a-> repeat(a, 0),     AT, rand(Float32, 5, 4))
        @test compare(a-> repeat(a, 4, 0),  AT, rand(Float32, 10, 15))
    end

    @testset "permutedims" begin
        @test compare(x->permutedims(x, [1, 2]), AT, rand(Float32, 4, 4))

        inds = rand(1:100, 150, 150)
        @test compare(x->permutedims(view(x, inds, :), (3, 2, 1)), AT, rand(Float32, 100, 100))
    end

    @testset "circshift" begin
        @test compare(x->circshift(x, (0,1)), AT, reshape(Vector(1:16), (4,4)))
    end

    @testset "copy" begin
        a = AT([1])
        b = copy(a)
        fill!(b, 0)
        @test Array(b) == [0]
        @test Array(a) == [1]
    end

    @testset "input output" begin
        # compact=false to avoid type aliases
        replstr(x, kv::Pair...) = sprint((io,x) -> show(IOContext(io, :compact => false, :limit => true, :displaysize => (24, 80), kv...), MIME("text/plain"), x), x)
        showstr(x, kv::Pair...) = sprint((io,x) -> show(IOContext(io, :limit => true, :displaysize => (24, 80), kv...), x), x)

        @testset "showing" begin
            # vectors and non-vector arrays showing
            # are handled differently in base/arrayshow.jl
            A = AT(Int64[1])
            B = AT(Int64[1 2;3 4])

            msg = replstr(A)
            @test occursin(Regex("^1-element $AT{Int64,\\s?1.*}:\n 1\$"), msg)

            # # result of e.g. `print` differs on 32bit and 64bit machines
            # due to different definition of `Int` type
            # print([1]) shows as [1] on 64bit but Int64[1] on 32bit
            msg = showstr(A)
            @test msg == "[1]" || msg == "Int64[1]"

            msg = replstr(B)
            @test occursin(Regex("^2×2 $AT{Int64,\\s?2.*}:\n 1  2\n 3  4\$"), msg)

            msg = showstr(B)
            @test msg == "[1 2; 3 4]" || msg == "Int64[1 2; 3 4]"

            # the printing of Adjoint depends on global state
            msg = replstr(A')
            @test occursin(Regex("^1×1 Adjoint{Int64,\\s?$AT{Int64,\\s?1.*}}:\n 1\$"), msg) ||
                occursin(Regex("^1×1 LinearAlgebra.Adjoint{Int64,\\s?$AT{Int64,\\s?1.*}}:\n 1\$"), msg) ||
                occursin(Regex("^1×1 adjoint\\(::$AT{Int64,\\s?1.*}\\) with eltype Int64:\n 1\$"), msg)
        end
    end
end


text:
function ADMFile()
    return ADMFile(())
end

function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_conventions(obj::NDMFile)
    return jcall(obj, "getConventions", IERSConventions, ())
end

function get_creation_date(obj::NDMFile)
    return jcall(obj, "getCreationDate", AbsoluteDate, ())
end

function get_data_context(obj::NDMFile)
    return jcall(obj, "getDataContext", DataContext, ())
end

function get_format_version(obj::NDMFile)
    return jcall(obj, "getFormatVersion", jdouble, ())
end

function get_header_comment(obj::NDMFile)
    return jcall(obj, "getHeaderComment", List, ())
end

function get_mission_reference_date(obj::NDMFile)
    return jcall(obj, "getMissionReferenceDate", AbsoluteDate, ())
end

function get_mu(obj::NDMFile)
    return jcall(obj, "getMu", jdouble, ())
end

function get_originator(obj::NDMFile)
    return jcall(obj, "getOriginator", JString, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function set_conventions(obj::NDMFile, arg0::IERSConventions)
    return jcall(obj, "setConventions", void, (IERSConventions,), arg0)
end

function set_creation_date(obj::NDMFile, arg0::AbsoluteDate)
    return jcall(obj, "setCreationDate", void, (AbsoluteDate,), arg0)
end

function set_data_context(obj::NDMFile, arg0::DataContext)
    return jcall(obj, "setDataContext", void, (DataContext,), arg0)
end

function set_format_version(obj::NDMFile, arg0::jdouble)
    return jcall(obj, "setFormatVersion", void, (jdouble,), arg0)
end

function set_header_comment(obj::NDMFile, arg0::List)
    return jcall(obj, "setHeaderComment", void, (List,), arg0)
end

function set_mission_reference_date(obj::NDMFile, arg0::AbsoluteDate)
    return jcall(obj, "setMissionReferenceDate", void, (AbsoluteDate,), arg0)
end

function set_mu(obj::NDMFile, arg0::jdouble)
    return jcall(obj, "setMu", void, (jdouble,), arg0)
end

function set_originator(obj::NDMFile, arg0::JString)
    return jcall(obj, "setOriginator", void, (JString,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end