text:
function test_scope(str, offset)
    server = LanguageServer.LanguageServerInstance(false,false,false)
    x = CSTParser.parse(str,true)
    s = LanguageServer.TopLevelScope(LanguageServer.ScopePosition("none", offset), LanguageServer.ScopePosition("none", 0), false, Dict(), LanguageServer.EXPR[], Symbol[], false, true, Dict{String,Set{String}}("toplevel" => Set{String}()), Dict{String,Set{String}}("toplevel" => Set{String}()), [])
    LanguageServer.toplevel(x, s, server)
    s.current.offset = 0
    y = LanguageServer._scope(x, s, server)
    if !isempty(s.stack) && last(s.stack) isa CSTParser.EXPR{CSTParser.MacroName}
        y = last(s.stack)
    end
    ns = isempty(s.namespace) ? "toplevel" : join(s.namespace, ".")
    yname = string(Expr(y))
    yname in keys(s.symbols) || yname in s.imported_names[ns]
end

function test_undefvar(str, offset = 0)
    server = LanguageServer.LanguageServerInstance(false,false,false)
    x = CSTParser.parse(str,true)
    s = LanguageServer.TopLevelScope(LanguageServer.ScopePosition("none", typemax(Int)), LanguageServer.ScopePosition("none", 0), false, Dict(), LanguageServer.EXPR[], Symbol[], false, true, Dict{String,Set{String}}("toplevel" => Set{String}()), Dict{String,Set{String}}("toplevel" => Set{String}()), [])
    LanguageServer.toplevel(x, s, server)
    L = LanguageServer.LintState([], [], [])
    s.current.offset = 0
    LanguageServer.lint(x, s, L, server, true)
    isempty(L.diagnostics)
end

@testset "scoping" begin
for f in [test_scope,test_undefvar]
    @testset "$(string(f))" begin
        @testset "for" begin
            @test f("""
            for iter in 1:10
                iter
            end
            """, 24)
            @test f("""
            for iter1 in 1:10, iter2 = 1:10, iter3 ∈ 1:30
                iter1
                iter2
                iter3
            end
            """, 53)
            @test f("""
            for iter1 in 1:10, iter2 = 1:10, iter3 ∈ 1:30
                iter1
                iter2
                iter3
            end
            """, 61)
            @test f("""
            for iter1 in 1:10, iter2 = 1:10, iter3 ∈ 1:30
                iter1
                iter2
                iter3
            end
            """, 69)
        end

        @testset "generators" begin
            @test f("""iter for iter in 1""", 2)
            @test f("""iter for iter in 1""", 11)
        end

        @testset "assignment" begin
            @test f("""
            var1 = 12345
            var1
            """, 14)
            @test f("""
            var1, var2 = 1,2
            var1
            var2
            """, 19)
            @test f("""
            var1, var2 = 1,2
            var1
            var2
            """, 24)
            @test f("""
            var1 = var2 = 12345
            var2
            """, 22)
        end

        @testset "anon func" begin
            @test f("arg -> arg * arg", 10)
            @test f("(arg1, arg2) -> arg1 * arg2", 18)
        end

        @testset "let" begin
            @test f("""
            let var1 = 1
                var1
            end
            """, 19)
            @test f("""
            let var1 = 1, var2 = 1
                var2
            end
            """, 30)
        end

        @testset "do" begin
            @test f("""
            sin([1]) do var1
                var1
            end
            """, 22)
            @test f("""
            sin([1]) do var1, var2
                var2
            end
            """, 28)
        end

        @testset "try" begin
            @test f("""
            try
                sin(1)
            catch err
                sin(err)
            end
            """, 34)
        end

        @testset "functions" begin
            @test f("""
            function func(x)
                func(x)
            end
            func(1)
            """, 24)
            @test f("""
            function func(x)
                func(x)
            end
            func(1)
            """, 35)
            @test f("""
            function func(arg = 1)
                arg
            end
            func(1)
            """, 29)
            @test f("""
            function func(arg::T) where T
                arg
            end
            func(1)
            """, 36)
        end

        @testset "macros" begin
            @test f("""
            macro mac(x)
                x
            end
            @mac 1
            """, 25)
        end

        @testset "datatypes" begin
            @test f("""
            struct name end
            name
            """, 19)
            @test f("""
            mutable struct name end
            name
            """, 26)
            @test f("""
            abstract type name end
            name
            """, 25)
            @test f("""
            primitive type name 8 end
            name
            """, 28)
        end

        @testset "modules" begin
            @test f("""
            module name end
            name
            """, 19)
            @test f("""
            baremodule name end
            name
            """, 23)
        end

        @testset "misc" begin
            @test f("""
            function(file)
                newstr = file
            end
            """, 29) # julia-vscode issue 241

            @test f("""
            using CSTParser: parse_operator
            parse_operator
            """, 34) # julia-vscode issue 222

            @test !f("""
            module A
            using CSTParser: parse_operator
            ParseState
            end
            """, 43) # julia-vscode issue 222
        end
    end
end
end


text:
# This script process the simulation data using non-linear least squares fitting
# for both the original Reference Region Model (NRRM) and the Constrained approach (CNRRM)

# Run this script by launching julia then entering: include("path/to/a02_DEMO_NonLinear.jl")

# Get current file path
curFilePath = Base.source_path() # Path to current file
workDir = dirname(curFilePath) # Our working directory (directory of current file)
## If the above doesn't work, then user has to manually enter the location of the CRRM directory
## i.e. wordDir = "/path/to/CLRRM/Directory"

# Figure out how many cores we have to work with
import Hwloc
topology = Hwloc.topology_load()
counts = Hwloc.histmap(topology)
nCores = counts[:Core]
if nprocs()<nCores
  addprocs(nCores-nprocs()+1)
end

## Load packages
using DCEMRI  # Needed for the levenberg-marquardt fitting
using MAT      # Needed for loading/saving .mat files

# Some required directories/files
warmupFile = workDir * "/data/QIBA-ToftsV6/QIBAv6-Mini4Jl-Warmup.mat"
auxCode = workDir * "/jlfiles/auxCode.jl"
refRegCode = workDir * "/jlfiles/refRegionFunction.jl"

inFile = workDir * "/data/demoData4Julia.mat" # File to load data from
outFile = workDir * "/data/demoDataFromJulia.mat" # File to save results to

# Load pre-requisite code
include(auxCode)

## Warmup run to make sure all is well
println("==== Warmup Run ====")
results = fitdata(datafile=warmupFile, models=[2])
rm("output.mat") # Delete output since it isn't needed

# If warmup succeeds, then proceed by loading Reference Region Models
include(refRegCode)

println("==== Processing data ====")

# Load .mat data
matData = matread(inFile)
# Assign data from .mat file to variables in Julia, for simplicity
Ct = matData["Ct"]
Crr = vec(matData["Crr"])
t = vec(matData["t"])

# Do the fitting
# Non-Linear Reference Region Model
# (Unconstrained and Constrained are both run if estKepRR is set to zero)
(pkParamsCN, residCN, estimatedKepRR, pkParamsN, residN, runtimeN, runtimeCN) = fitCNRRM(Ct, Crr, t, estKepRR=0.0, doTime=true)

# Output results to a .mat file
results = Dict()
results["estParamsNRRM"] = pkParamsN
results["estParamsCNRRM"] = pkParamsCN
results["runtimeNRRM"] = runtimeN
results["runtimeCNRRM"] = runtimeCN
results["estKepRR_N"] = estimatedKepRR
matwrite(outFile, results)

println("==== Finished ====")
println("==================")

ext: jl

max_stars_repo_name: MPUmri/CRRM

max_issues_repo_name: MPUmri/CRRM

max_forks_repo_name: MPUmri/CRRM


text:
#= %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Example: Find the surface origin of water for some interior box %
 %                                                                   %
 % This is equivalent to solving a sensitivity problem:              %
 % The total volume is V = v^T c , where v is the volume of a
 % given interior box,
 % and c is the fraction of volume from a given source which         %
 % satisfies the equation A c = δ.                                   %
 % Next, dV/d(δ) = A^(-T) v, and dV/d(δ) is exactly the volume       %
 % originating from each source.
 % Very similar mathematically to determining how the ocean is filled.                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
=#
using Revise
using TMI, PyPlot, PyCall

TMIversion = "modern_90x45x33_GH10_GH12"
A, Alu, γ, TMIfile, L, B = config_from_nc(TMIversion)

# % choose an interior location X (Xlon[lon], Xlat [lat], Xdepth [m depth]).
#% -7.38, 115.26E

xlon = 125.26; # deg E.
xlat = -6.38;  # deg N.
xdepth = 3000;  # meters.
loc = (xlon,xlat,xdepth)

origin = surfaceorigin(loc, Alu, γ)

# view the surface
cntrs = -8:1:0
tlabel="Surface origin of interior location [log10(fraction)]"
planviewplot(origin,cntrs,γ,titlelabel=tlabel)


text:
"""
    Fast adjust potential array data
    ================================

    *_r
        indicates that function arguments include (x, y, z), i.e., real coordinates.

    *_g
        indicates that function arguments include (xg, yg, zg), i.e., grid coordinates.

"""
struct FastAdjust2D
    pa::Array{Float64,3}
    el::Array{Bool,2}
    num::Int64
    nx::Int64
    ny::Int64
    dx::Float64
    dy::Float64
    x0::Float64
    y0::Float64
    xn::Float64
    yn::Float64
end;

struct FastAdjust3D
    pa::Array{Float64,4}
    el::Array{Bool,3}
    num::Int64
    nx::Int64
    ny::Int64
    nz::Int64
    dx::Float64
    dy::Float64
    dz::Float64
    x0::Float64
    y0::Float64
    z0::Float64
    xn::Float64
    yn::Float64
    zn::Float64
end;

FastAdjust = Union{FastAdjust2D, FastAdjust3D}

function h5read_pa(fil, r0::Vector{Float64}=Vector([0.0, 0.0, 0.0]))
    """ read a HDF5 file of potential array data
    """
    pa = h5read(fil, "pa")
    el = h5read(fil, "electrode")
    attrs = h5readattr(fil, "/")
    if ndims(pa) == 3
        num, ny, nx = size(pa)
        dx, dy = attrs["dx"], attrs["dy"]
        xn = r0[1] + (nx - 1) * dx
        yn = r0[2] + (ny - 1) * dy
        fa = FastAdjust2D(pa, el, num, nx, ny, dx, dy, r0[1], r0[2], xn, yn)
    elseif ndims(pa) == 4
        num, nz, ny, nx = size(pa)
        dx, dy, dz = attrs["dx"], attrs["dy"], attrs["dz"]
        xn = r0[1] + (nx - 1) * dx
        yn = r0[2] + (ny - 1) * dy
        zn = r0[3] + (nz - 1) * dz
        fa = FastAdjust3D(pa, el, num, nx, ny, nz, dx, dy, dz,
                        r0[1], r0[2], r0[3], xn, yn, zn)
    else
        ErrorException("file not compatible with FastAdjust2D or FastAdjust3D")
    end
    return fa
end

function potential(fa::FastAdjust, voltages::Vector{Float64})
    """ electric potential for a set of voltages
    """
    return dropdims(sum(fa.pa .* voltages, dims=1), dims=1)
end

function grid_r(fa::FastAdjust2D, r::Vector{Float64})
    """ convert real coordinates into grid coordinates
    """
    xg = 1.0 + (r[1] - fa.x0) / fa.dx
    yg = 1.0 + (r[2] - fa.y0) / fa.dy
    zg = NaN
    return xg, yg, zg
end

function grid_r(fa::FastAdjust3D, r::Vector{Float64})
    """ convert real coordinates into grid coordinates
    """
    xg = 1.0 + (r[1] - fa.x0) / fa.dx
    yg = 1.0 + (r[2] - fa.y0) / fa.dy
    zg = 1.0 + (r[3] - fa.z0) / fa.dz
    return xg, yg, zg
end

function electrode_g(fa::FastAdjust2D, xg::Float64, yg::Float64, zg::Float64)
    """ are any of the grid point adjacent to (xg, yg, zg) an electrode?
    """
    xn = Int64(floor(xg))
    yn = Int64(floor(yg))
    if (1 <= xn < fa.nx) && (1 <= yn < fa.ny)
        return any(fa.el[yn : yn + 1, xn : xn + 1])
    else
        # everywhere outside the pa is an electrode
        return true
    end
end

function electrode_g(fa::FastAdjust3D, xg::Float64, yg::Float64, zg::Float64)
    """ are any of the grid point adjacent to (xg, yg, zg) an electrode?
    """
    xn = Int64(floor(xg))
    yn = Int64(floor(yg))
    zn = Int64(floor(zg))
    if (1 <= xn < fa.nx) && (1 <= yn < fa.ny) && (1 <= zn < fa.nz)
        return any(fa.el[zn : zn + 1, yn : yn + 1, xn : xn + 1])
    else
        # everywhere outside the pa is an electrode
        return true
    end
end

function potential_g(fa::FastAdjust2D, voltages::Vector{Float64}, xg::Float64, yg::Float64, zg::Float64)
    """ electric potential at grid coord
    """
    try
        # integer coords (floor rounding)
        xn = Int64(floor(xg))
        yn = Int64(floor(yg))
        # get potential
        wx = xg - xn
        wy = yg - yn
        pa = fa.pa[:, yn: yn + 1, xn : xn + 1]
        phi = dropdims(sum(pa .* voltages, dims=1), dims=1)
        result = (1.0 - wx) * (1.0 - wy) * phi[1, 1] +
                 wx * (1.0 - wy) * phi[1, 2] +
                 (1.0 - wx) * wy * phi[2, 1] +
                 wx * wy * phi[2, 2]
        return result
    catch
        return NaN
    end
end

function potential_g(fa::FastAdjust3D, voltages::Vector{Float64}, xg::Float64, yg::Float64, zg::Float64)
    """ electric potential at grid coord
    """
    try
        # integer coords (floor rounding)
        xn = Int64(floor(xg))
        yn = Int64(floor(yg))
        zn = Int64(floor(zg))
        pa = fa.pa[:, zn: zn + 1, yn: yn + 1, xn : xn + 1]
        phi = dropdims(sum(pa .* voltages, dims=1), dims=1)
        ## interpolate along x
        wx = (xg - xn)
        c00 = phi[1, 1, 1] * (1 - wx) .+ phi[1, 1, 2] * wx
        c01 = phi[2, 1, 1] * (1 - wx) .+ phi[2, 1, 2] * wx
        c10 = phi[1, 2, 1] * (1 - wx) .+ phi[1, 2, 2] * wx
        c11 = phi[2, 2, 1] * (1 - wx) .+ phi[2, 2, 2] * wx
        ## interpolate along y
        wy = (yg - yn)
        c0 = c00 * (1 - wy) .+ c10 * wy
        c1 = c01 * (1 - wy) .+ c11 * wy
        ## interpolate along z
        wz = (zg - zn)
        return c0 * (1 - wz) .+ c1 * wz
    catch
        return NaN
    end
end

function field_g(fa::FastAdjust2D, voltages::Vector{Float64}, xg::Float64, yg::Float64, zg::Float64)
    """ electric field at grid coord
    """
    ex = (potential_g(fa, voltages, xg + 0.5, yg, zg) -
          potential_g(fa, voltages, xg - 0.5, yg, zg)) / fa.dx
    ey = (potential_g(fa, voltages, xg, yg + 0.5, zg) -
          potential_g(fa, voltages, xg, yg - 0.5, zg)) / fa.dy
    ez = 0.0
    return ex, ey, ez
end

function field_g(fa::FastAdjust3D, voltages::Vector{Float64}, xg::Float64, yg::Float64, zg::Float64)
    """ electric field at grid coord
    """
    ex = (potential_g(fa, voltages, xg + 0.5, yg, zg) -
          potential_g(fa, voltages, xg - 0.5, yg, zg)) / fa.dx
    ey = (potential_g(fa, voltages, xg, yg + 0.5, zg) -
          potential_g(fa, voltages, xg, yg - 0.5, zg)) / fa.dy
    ez = (potential_g(fa, voltages, xg, yg, zg + 0.5) -
          potential_g(fa, voltages, xg, yg, zg - 0.5)) / fa.dz
    return ex, ey, ez
end

function amp_field_g(fa::FastAdjust, voltages::Vector{Float64}, xg::Float64, yg::Float64, zg::Float64)
    """ amplitude of the field at grid coord
    """
    ex, ey, ez = field_g(fa, voltages, xg, yg, zg)
    return sqrt(ex^2.0 + ey^2.0 + ez^2.0)
end

function grad_field_g(fa::FastAdjust2D, voltages::Vector{Float64}, xg::Float64, yg::Float64, zg::Float64)
    """ gradient in the amplitude of the electric field at grid coord
    """
    dfx = (amp_field_g(fa, voltages, xg + 0.5, yg, zg) -
           amp_field_g(fa, voltages, xg - 0.5, yg, zg)) / fa.dx
    dfy = (amp_field_g(fa, voltages, xg, yg + 0.5, zg) -
           amp_field_g(fa, voltages, xg, yg - 0.5, zg)) / fa.dy
    dfz = 0.0
    return dfx, dfy, dfz
end

function grad_field_g(fa::FastAdjust3D, voltages::Vector{Float64}, xg::Float64, yg::Float64, zg::Float64)
    """ gradient in the amplitude of the electric field at grid coord
    """
    dfx = (amp_field_g(fa, voltages, xg + 0.5, yg, zg) -
           amp_field_g(fa, voltages, xg - 0.5, yg, zg)) / fa.dx
    dfy = (amp_field_g(fa, voltages, xg, yg + 0.5, zg) -
           amp_field_g(fa, voltages, xg, yg - 0.5, zg)) / fa.dy
    dfz = (amp_field_g(fa, voltages, xg, yg, zg + 0.5) -
           amp_field_g(fa, voltages, xg, yg, zg - 0.5)) / fa.dz
    return dfx, dfy, dfz
end


text:
include("../diagonal_sbp.jl")
include("legacy_functions.jl")

if length(ARGS) != 0
    level = parse(Int,ARGS[1])
    Iterations = parse(Int,ARGS[2])
else
    level = 10
    Iterations = 10000
end

using LinearAlgebra
using SparseArrays
using Plots
using IterativeSolvers
using BenchmarkTools
using MAT

function e(i,n)
    # A = Matrix{Float64}(I,n,n)
    # return A[:,i]
    out = spzeros(n)
    out[i] = 1.0
    return out
end

function eyes(n)
    # return Matrix{Float64}(I,n,n)
    out = spzeros(n,n)
    for i in 1:n
        out[i,i] = 1.0
    end
    return out
end

function u(x,y)
    return sin.(π*x .+ π*y)
end

function Diag(A)
    # Self defined function that is similar to Matlab Diag
    return Diagonal(A[:])
end

h_list_x = [1/2^1, 1/2^2, 1/2^3, 1/2^4, 1/2^5, 1/2^6, 1/2^7, 1/2^8, 1/2^9, 1/2^10, 1/2^11, 1/2^12, 1/2^13, 1/2^14]
h_list_y = [1/2^1, 1/2^2, 1/2^3, 1/2^4, 1/2^5, 1/2^6, 1/2^7, 1/2^8, 1/2^9, 1/2^10, 1/2^11, 1/2^12, 1/2^13, 1/2^14]

function Operators_2d(i, j, h_list_x, h_list_y; p=2)
    hx = h_list_x[i];
    hy = h_list_y[j];

    x = range(0,step=hx,1);
    y = range(0,step=hy,1);
    m_list = 1 ./h_list_x;
    n_list = 1 ./h_list_y;

    # Matrix Size
    N_x = Integer(m_list[i]);
    N_y = Integer(n_list[j]);

    (D1x, HIx, H1x, r1x) = diagonal_sbp_D1(p,N_x,xc=(0,1));
    (D2x, S0x, SNx, HI2x, H2x, r2x) = diagonal_sbp_D2(p,N_x,xc=(0,1));

    (D1y, HIy, H1y, r1y) = diagonal_sbp_D1(p,N_y,xc=(0,1));
    (D2y, S0y, SNy, HI2y, H2y, r2y) = diagonal_sbp_D2(p,N_y,xc=(0,1));

    # BSx = sparse(SNx - S0x);
    # BSy = sparse(SNy - S0y);
    BSx = SNx - S0x
    BSy = SNy - S0y

    # Forming 2d Operators
    # e_1x = sparse(e(1,N_x+1));
    # e_Nx = sparse(e(N_x+1,N_x+1));
    # e_1y = sparse(e(1,N_y+1));
    # e_Ny = sparse(e(N_y+1,N_y+1));
    e_1x = e(1,N_x+1);
    e_Nx = e(N_x+1,N_x+1);
    e_1y = e(1,N_x+1);
    e_Ny = e(N_y+1,N_y+1);

    # I_Nx = sparse(eyes(N_x+1));
    # I_Ny = sparse(eyes(N_y+1));
    I_Nx = eyes(N_x+1);
    I_Ny = eyes(N_y+1);

    e_E = kron(e_Nx,I_Ny);
    e_W = kron(e_1x,I_Ny);
    e_S = kron(I_Nx,e_1y);
    e_N = kron(I_Nx,e_Ny);

    E_E = kron(sparse(Diag(e_Nx)),I_Ny);   # E_E = e_E * e_E'
    E_W = kron(sparse(Diag(e_1x)),I_Ny);
    E_S = kron(I_Nx,sparse(Diag(e_1y)));
    E_N = sparse(kron(I_Nx,sparse(Diag(e_Ny))));

    D1_x = kron(D1x,I_Ny);
    D1_y = kron(I_Nx,D1y);

    D2_x = kron(D2x,I_Ny);
    D2_y = kron(I_Nx,D2y);
    D2 = D2_x + D2_y

    HI_x = kron(HIx,I_Ny);
    HI_y = kron(I_Nx,HIy);

    H_x = kron(H1x,I_Ny);
    H_y = kron(I_Nx,H1y);

    BS_x = kron(BSx,I_Ny);
    BS_y = kron(I_Nx,BSy);

    HI_tilde = kron(HIx,HIx);
    H_tilde = kron(H1x,H1y);

    return (D1_x, D1_y, D2_x, D2_y, D2, HI_x, HI_y, BS_x, BS_y, HI_tilde, H_tilde, I_Nx, I_Ny, e_E, e_W, e_S, e_N, E_E, E_W, E_S, E_N)
end

function Assembling_matrix(level;p=2)
    i = j = level
    hx = h_list_x[i];
    hy = h_list_y[j];

    x = range(0,step=hx,1);
    y = range(0,step=hy,1);
    m_list = 1 ./h_list_x;
    n_list = 1 ./h_list_y;

    N_x = Integer(m_list[i]);
    N_y = Integer(n_list[j]);

    Nx = N_x + 1;
    Ny = N_y + 1;

    (D1_x, D1_y, D2_x, D2_y, D2, HI_x, HI_y, BS_x, BS_y, HI_tilde, H_tilde, I_Nx, I_Ny, e_E, e_W, e_S, e_N, E_E, E_W, E_S, E_N) = Operators_2d(i,j,h_list_x,h_list_y,p=p);
    analy_sol = u(x,y');

    # Penalty Parameters
    tau_E = -13/hx;
    tau_W = -13/hx;
    tau_N = -1;
    tau_S = -1;

    beta = 1;

    # Forming SAT terms

    ## Formulation 1
    SAT_W = tau_W*HI_x*E_W + beta*HI_x*BS_x'*E_W;
    SAT_E = tau_E*HI_x*E_E + beta*HI_x*BS_x'*E_E;

    # SAT_S = tau_S*HI_y*E_S*D1_y
    # SAT_N = tau_N*HI_y*E_N*D1_y

    SAT_S = tau_S*HI_y*E_S*BS_y;
    SAT_N = tau_N*HI_y*E_N*BS_y;

    SAT_W_r = tau_W*HI_x*E_W*e_W + beta*HI_x*BS_x'*E_W*e_W;
    SAT_E_r = tau_E*HI_x*E_E*e_E + beta*HI_x*BS_x'*E_E*e_E;
    SAT_S_r = tau_S*HI_y*E_S*e_S;
    SAT_N_r = tau_N*HI_y*E_N*e_N;

    (alpha1,alpha2,alpha3,alpha4,beta) = (tau_N,tau_S,tau_W,tau_E,beta);

    g_W = sin.(π*y);
    g_E = -sin.(π*y);
    g_S = -π*cos.(π*x);
    # g_N = -π*cos.(π*x)
    g_N = π*cos.(π*x .+ π);

    # Solving with CPU
    A = D2 + SAT_W + SAT_E + SAT_S + SAT_N;

    b = -2π^2*u(x,y')[:] + SAT_W_r*g_W + SAT_E_r*g_E + SAT_S_r*g_S + SAT_N_r*g_N;

    A = -H_tilde*A;
    b = -H_tilde*b;

    return (A,b,H_tilde,Nx,Ny)
end

function prolongation_matrix(N)
    # SBP preserving
    # N = 2^level + 1
    odata = spzeros(2*N-1,N)
    for i in 1:2*N-1
        if i % 2 == 1
            odata[i,div(i+1,2)] = 1
        else
            odata[i,div(i,2)] = 1/2
            odata[i,div(i,2)+1] = 1/2
        end
    end
    return odata
end

function restriction_matrix(N)
    # SBP preserving
    odata = spzeros(div(N+1,2),N)
    odata[1,1] = 1/2
    odata[1,2] = 1/2
    odata[end,end-1] = 1/2
    odata[end,end] = 1/2
    for i in 2:div(N+1,2)-1
        odata[i,2*i-2] = 1/4
        odata[i,2*i-1] = 1/2
        odata[i,2*i] = 1/4
    end
    return odata
end

function restriction_matrix_normal(N)
    # SBP preserving
    odata = spzeros(div(N+1,2),N)
    odata[1,1] = 1/2
    odata[1,2] = 1/4
    odata[end,end-1] = 1/4
    odata[end,end] = 1/2
    for i in 2:div(N+1,2)-1
        odata[i,2*i-2] = 1/4
        odata[i,2*i-1] = 1/2
        odata[i,2*i] = 1/4
    end
    return odata
end

function prolongation_2d(N)
    prolongation_1d = prolongation_matrix(N)
    prolongation_2d = kron(prolongation_1d,prolongation_1d)
    return prolongation_2d
end

function restriction_2d(N)
    restriction_1d = restriction_matrix_normal(N)
    # restriction_1d = restriction_matrix(N)
    restriction_2d = kron(restriction_1d,restriction_1d)
    return restriction_2d
end

function restriction_2d_GPU(N)
    return CUDA.CUSPARSE.CuSparseMatrixCSC(restriction_2d(N))
end

function prolongation_2d_GPU(N)
    return CUDA.CUSPARSE.CuSparseMatrixCSC(prolongation_2d(N))
end

function CG_CPU(A,b,x;maxiter=length(b),abstol=sqrt(eps(real(eltype(b)))),direct_sol=0,H_tilde=0)
    r = b - A * x;
    p = r;
    rsold = r' * r
    # Ap = spzeros(length(b))
    Ap = similar(b);

    num_iter_steps = 0
    norms = [sqrt(rsold)]
    errors = []
    if direct_sol != 0 && H_tilde != 0
        append!(errors,sqrt(direct_sol' * A * direct_sol))
    end
    # @show rsold
    for step = 1:maxiter
        num_iter_steps += 1
        mul!(Ap,A,p);
        alpha = rsold / (p' * Ap)
        x .= x .+ alpha * p;
        r .= r .- alpha * Ap;
        rsnew = r' * r
        append!(norms,sqrt(rsnew))
        if direct_sol != 0 && H_tilde != 0
            error = sqrt((x - direct_sol)' * A * (x - direct_sol))
            append!(errors,error)
        end
        if sqrt(rsnew) < abstol
              break
        end
        p = r + (rsnew / rsold) * p;
        rsold = rsnew
    end
    return (num_iter_steps,norms,errors)
end

function regularCG!(A,b,x,H1,exact;maxiter=length(b)^2,abstol=sqrt(eps(real(eltype(b)))))
    r = b - A * x;
    p = r;
    rsold = r'*r

    num_iter_steps = 0
    norms = [sqrt(r'*H1*r)]

    diff = x-exact
    err = sqrt(diff'*A*diff)

    E = [err]
    for step = 1:maxiter
        Ap = A*p;
        num_iter_steps += 1
        alpha = rsold/(p'*Ap)
        x .= x .+ alpha * p;

        diff = x-exact
        err = sqrt(diff'*A*diff)
        append!(E, err)

        r .= r .- alpha * Ap;
        rsnew = r' * r

        norm_r = sqrt(rsnew)
        append!(norms,norm_r)

        if norm_r < abstol
            break
        end

        beta = rsnew/rsold;
        p = r + beta * p;
        rsold = rsnew
    end

    return E, num_iter_steps, norms
end

function jacobi_brittany!(x,A,b;maxiter=3, ω = 2/3)

    Pinv = Diagonal(1 ./ diag(A))
    P = Diagonal(diag(A))
    Q = A-P

    for j in 1:maxiter
        x[:] = ω * Pinv*(b .- Q*x[:]) + (1 - ω)*x[:]
    end
end

function modified_richardson!(x,A,b;maxiter=3,ω=0.15)
    for _ in 1:maxiter
        x[:] .= x[:] + ω*(b .- A*x[:])
    end
end

function modified_richardson_GPU!(x_GPU,A_GPU,b_GPU;maxiter=3,ω=0.15)
    for _ in 1:maxiter
        x_GPU .= x_GPU .+ ω*(b_GPU .- A_GPU*x_GPU)
    end
end

function Two_level_multigrid(A,b,Nx,Ny,A_2h;nu=3,NUM_V_CYCLES=1,SBPp=2)
    v_values = Dict(1=>zeros(Nx*Ny))
    rhs_values = Dict(1 => b)
    N_values = Dict(1 => Nx)
    N_values[2] = div(Nx+1,2)

    x = zeros(length(b));
    v_values[1] = x

    for cycle_number in 1:NUM_V_CYCLES
        # jacobi_brittany!(v_values[1],A,b;maxiter=nu);
        modified_richardson!(v_values[1],A,b,maxiter=nu)
        r = b - A*v_values[1];
        f = restriction_2d(Nx) * r;
        v_values[2] = A_2h \ f

        # println("Pass first part")
        e_1 = prolongation_2d(N_values[2]) * v_values[2];
        v_values[1] = v_values[1] + e_1;
        # println("After coarse grid correction, norm(A*x-b): $(norm(A*v_values[1]-b))")
        # jacobi_brittany!(v_values[1],A,b;maxiter=nu);
        modified_richardson!(v_values[1],A,b,maxiter=nu)
    end
    return (v_values[1],norm(A * v_values[1] - b))
end

function Two_level_multigrid_GPU(A_GPU,b_GPU,Nx,Ny,A_2h;nu=3,NUM_V_CYCLES=1,SBPp=2)
    v_values = Dict(1=>CuArray(zeros(Nx*Ny)))
    rhs_values = Dict(1 => b_GPU)
    N_values = Dict(1 => Nx)
    N_values[2] = div(Nx+1,2)

    x = CuArray(zeros(length(b_GPU)));
    v_values[1] = x

    for cycle_number in 1:NUM_V_CYCLES
        # jacobi_brittany!(v_values[1],A,b;maxiter=nu);
        modified_richardson_GPU!(v_values[1],A_GPU,b_GPU,maxiter=nu)
        r = b_GPU - A_GPU*v_values[1];
        f = Array(restriction_2d_GPU(Nx) * r);
        v_values[2] = CuArray(A_2h \ f)

        # println("Pass first part")
        e_1 = prolongation_2d_GPU(N_values[2]) * v_values[2];
        v_values[1] = v_values[1] + e_1;
        # println("After coarse grid correction, norm(A*x-b): $(norm(A*v_values[1]-b))")
        # jacobi_brittany!(v_values[1],A,b;maxiter=nu);
        modified_richardson_GPU!(v_values[1],A_GPU,b_GPU,maxiter=nu)
    end
    return (v_values[1],norm(A_GPU * v_values[1] - b_GPU))
end

function precond_matrix(A, b, A_2h; m=3, solver="jacobi")
    #pre and post smoothing
    N = length(b)
    IN = sparse(Matrix(I, N, N))
    P = Diagonal(diag(A))
    Pinv = Diagonal(1 ./ diag(A))
    Q = P-A
    L = A - triu(A)
    U = A - tril(A)

    if solver == "jacobi"
       ω = 2/3
        H = ω*Pinv*Q + (1-ω)*IN
        R = ω*Pinv
        R0 = ω*Pinv
    elseif solver == "ssor"
        ω = 1.4  #this is just a guess. Need to compute ω_optimal (from jacobi method)
        B1 = (P + ω*U)\Matrix(-ω*L + (1-ω)*P)
        B2 = (P + ω*L)\Matrix(-ω*U + (1-ω)*P)
        H = B1*B2
        X = (P+ω*L)\Matrix(IN)

        R = ω*(2-ω)*(P+ω*U)\Matrix(P*X)
        R0 = ω*(2-ω)*(P+ω*U)\Matrix(P*X)
    elseif solver == "richardson"
        ω =ω_richardson
        H = IN - ω*A
        R = ω*IN
        # R = spzeros(N,N)
        R0 = ω*IN
    elseif solver == "richardson_chebyshev" #TODO: FIX ME FOR CHEB
        ω =ω_richardson
        H = IN - ω*A
        R = ω*IN
        R0 = ω*IN
    elseif solver == "chebyshev" #TODO: FIX ME FOR CHEB
        ω =ω_richardson
        H = IN - ω*A
        R = ω*IN
        R0 = ω*IN
    else
    end

    for i = 1:m-1
        R += H^i * R0
    end

    # (A_2h, b_2h, x_2h, H1_2h) = get_operators(SBPp, 2*h);
    A_2h = A_2h
    I_r = standard_restriction_matrix_2D(N)

    I_p = standard_prolongation_matrix_2D(length(b_2h))
    M = H^m * (R + I_p * (A_2h\Matrix(I_r*(IN - A * R)))) + R
    M_alternative = R + R + I_p * (A_2h \ Matrix(I_r * (IN - A*R)))

    M_v0 = I_p * (A_2h \ Matrix(I_r)) # No richardson iteration
    M_v1_no_post = ω*IN + I_p * (A_2h \ Matrix(I_r * (IN - ω*A)))

    M_v1 = (IN - ω*A)*(ω*IN + I_p * (A_2h \ Matrix(I_r * (IN - ω*A)))) + ω*IN # one pre and post richardson iteration

    M_v1_alternative = H * (R0 + I_p * (A_2h \ Matrix(I_r * H))) + R0 # alternative form

    M_test = H^m * R + R + H^m * (I_p * (A_2h \ Matrix(I_r * H^m))) # need to change, a generalized representation of M for richardson iteration

    # return (M, R, H, I_p, A_2h, I_r, IN)
    return (M_test,R,H,I_p,A_2h,I_r,IN)
end

function mg_preconditioned_CG(A,b,x;A_2h = A_2h_lu,maxiter=length(b),abstol=sqrt(eps(real(eltype(b)))),NUM_V_CYCLES=1,nu=3,use_galerkin=true,direct_sol=0,H_tilde=0,SBPp=2)
    Nx = Ny = Int(sqrt(length(b)))
    level = Int(log(2,Nx-1))
    # (A_2h,b_2h,H_tilde_2h,Nx_2h,Ny_2h) = Assembling_matrix(level-1,p=SBPp);
    # A_2h = lu(A_2h)
    r = b - A * x;
    # (M, R, H, I_p, A_2h, I_r, IN) = precond_matrix(A,b;m=nu,solver="jacobi",p=p)
    z = Two_level_multigrid(A,r,Nx,Ny,A_2h;nu=nu,NUM_V_CYCLES=1)[1]
    # z = M*r
    p = z;
    num_iter_steps = 0
    norms = [norm(r)]
    errors = []
    if direct_sol != 0 && H_tilde != 0
        append!(errors,sqrt(direct_sol' * A * direct_sol))
    end

    rzold = r'*z

    for step = 1:maxiter
    # for step = 1:5
        num_iter_steps += 1
        alpha = rzold / (p'*A*p)
        x .= x .+ alpha * p;
        r .= r .- alpha * A*p
        rs = r' * r
        append!(norms,sqrt(rs))
        if direct_sol != 0 && H_tilde != 0
            error = sqrt((x - direct_sol)' * A * (x - direct_sol))
            # @show error
            append!(errors,error)
        end
        if sqrt(rs) < abstol
            break
        end
        z = Two_level_multigrid(A,r,Nx,Ny,A_2h;nu=nu,NUM_V_CYCLES=1)[1]
        # z = M*r
        rznew = r'*z
        beta = rznew/(rzold);
        p = z + beta * p;
        rzold = rznew
    end
    # @show num_iter_steps
    return num_iter_steps, norms, errors
end

function mg_preconditioned_CG_GPU(A_GPU,b_GPU,x_GPU;A_2h = A_2h_lu,maxiter=length(b),abstol=sqrt(eps(real(eltype(b)))),NUM_V_CYCLES=1,nu=3,use_galerkin=true,direct_sol=0,H_tilde=0,SBPp=2)
    Nx = Ny = Int(sqrt(length(b_GPU)))
    level = Int(log(2,Nx-1))
    # (A_2h,b_2h,H_tilde_2h,Nx_2h,Ny_2h) = Assembling_matrix(level-1,p=SBPp);
    # A_2h = lu(A_2h)
    r_GPU = b_GPU - A_GPU * x_GPU;
    # (M, R, H, I_p, A_2h, I_r, IN) = precond_matrix(A,b;m=nu,solver="jacobi",p=p)
    z_GPU = Two_level_multigrid_GPU(A_GPU,r_GPU,Nx,Ny,A_2h;nu=nu,NUM_V_CYCLES=1)[1]
    # z = M*r
    p_GPU = z_GPU;
    num_iter_steps = 0
    norms = [norm(r_GPU)]
    errors = []
    if direct_sol != 0 && H_tilde != 0
        append!(errors,sqrt(direct_sol' * A * direct_sol))
    end

    rzold = r_GPU'*z_GPU

    for step = 1:maxiter
    # for step = 1:5
        num_iter_steps += 1
        alpha = rzold / (p_GPU'*A_GPU*p_GPU)
        x_GPU .= x_GPU .+ alpha * p_GPU;
        r_GPU .= r_GPU .- alpha * (A_GPU*p_GPU)
        rs = r_GPU' * r_GPU
        append!(norms,sqrt(rs))
        # if direct_sol != 0 && H_tilde != 0
        #     error = sqrt((x_GPU - direct_sol)' * A_GPU * (x_GPU - direct_sol))
        #     # @show error
        #     append!(errors,error)
        # end
        if sqrt(rs) < abstol
            break
        end
        z_GPU = Two_level_multigrid_GPU(A_GPU,r_GPU,Nx,Ny,A_2h;nu=nu,NUM_V_CYCLES=1)[1]
        # z = M*r
        rznew = r_GPU'*z_GPU
        beta = rznew/(rzold);
        p_GPU .= z_GPU .+ beta * p_GPU;
        rzold = rznew
    end
    # @show num_iter_steps
    return num_iter_steps, norms, errors
end

function test_preconditioned_CG(;level=level,nu=3,ω=2/3,SBPp=2)
    (A,b,H_tilde,Nx,Ny) = Assembling_matrix(level,p=SBPp);
    (A_2h,b_2h,H_tilde_2h,Nx_2h,Ny_2h) = Assembling_matrix(level-1,p=SBPp);
    A_2h_lu = lu(A_2h)
    direct_sol = A\b
    reltol = sqrt(eps(real(eltype(b))))
    x = zeros(Nx*Ny);
    abstol = norm(A*x-b) * reltol
    h = 1/(Nx-1)
    ω_richardson = 0.15
    p = SBPp
    (M, R, H, I_p, A_2h, I_r, IN) = precond_matrix(A,b,A_2h;m=nu,solver="richardson")

    cond_A_M = cond(M*A)
    x = zeros(Nx*Ny);
    iter_mg_cg, norm_mg_cg, error_mg_cg = mg_preconditioned_CG(A,b,x;A_2h = A_2h_lu, maxiter=length(b),abstol=abstol,NUM_V_CYCLES=1,nu=nu,use_galerkin=true,direct_sol=direct_sol,H_tilde=H_tilde,SBPp=SBPp)
    error_mg_cg_bound_coef = (sqrt(cond_A_M) - 1) / (sqrt(cond_A_M) + 1)
    error_mg_cg_bound = error_mg_cg[1] .* 2 .* error_mg_cg_bound_coef .^ (0:1:length(error_mg_cg)-1)
    scatter(log.(10,error_mg_cg),label="error_mg_cg", markercolor = "darkblue")
    xlabel!("Iterations")
    ylabel!("Log errors, base 10")
    plot!(log.(10,error_mg_cg_bound),label="error_mg_cg_bound",linecolor = "darkblue")

    cond_A = cond(Matrix(A))
    x0 = zeros(Nx*Ny)
    (E_cg, num_iter_steps_cg, norms_cg) = regularCG!(A,b,x0,H_tilde,direct_sol;maxiter=20000,abstol=abstol)

    scatter!(log.(10,E_cg),label="error_cg", markercolor = "darksalmon")
    error_cg_bound_coef = (sqrt(cond_A) - 1) / (sqrt(cond_A) + 1)
    error_cg_bound = E_cg[1] .* 2 .* error_cg_bound_coef .^ (0:1:length(E_cg)-1)
    plot!(log.(10,error_cg_bound),label="error_cg_bound", linecolor = "darksalmon")

    savefig("convergence_richardson.png")
    # my_solver = "jacobi"
    # x0 = zeros(Nx*Ny)
    # (E_mgcg, num_iter_steps_mgcg, norms_mgcg) = MGCG!(A,b,x0,H_tilde,direct_sol,my_solver;smooth_steps = nu,maxiter=20000,abstol=abstol)

    # plot(error_mg_cg,label="error_mg_cg")
    # plot!(error_mg_cg_bound,label="error_mg_cg_bound")
end


text:
# This file is a part of SimilaritySearch.jl

using SimilaritySearch
using Test

function testsorted(res, Vsorted)
    @info "========================", (typeof(res), length(res), maxlength(res)), length(Vsorted)
    @show res
    @show Vsorted
    @info "======================== starting ============="
    @test collect(res) == Vsorted
    @test minimum(res) == first(Vsorted)[end]
    @test maximum(res) == last(Vsorted)[end]
    @test argmin(res) == first(Vsorted)[1]
    @test argmax(res) == last(Vsorted)[1]

    @show res

    @test res.id == first.(Vsorted)
    @test res.dist == last.(Vsorted)

    pop!(Vsorted)
    pop!(res)
    @test res.id == first.(Vsorted)
    @test res.dist == last.(Vsorted)
    @test collect(res) == Vsorted

    popfirst!(Vsorted)
    popfirst!(res)

    @info "b   collect id:" => res.id
    @info "b collect dist:" => res.dist
    @info "b ========" => first.(Vsorted)
    @info "b ========" => last.(Vsorted)
    @show res.id
    @show res.dist
    @test res.id == first.(Vsorted)
    @test res.dist == last.(Vsorted)
    @test collect(res) == Vsorted

end

function testsorted2(res, st, Vsorted)
    @info "========================", (typeof(res), length(res, st), length(res), maxlength(res), st), length(Vsorted)
    @show res
    @show Vsorted
    @info "======================== starting ============="
    @test length(res, st) == length(res)
    @test minimum(res, st) == minimum(res)
    @test maximum(res, st) == maximum(res)
    @test argmin(res, st) == argmin(res)
    @test argmax(res, st) == argmax(res)

    @show res.id, res.dist, st, Vsorted
    @test collect(res) == Vsorted
    @test minimum(res, st) == first(Vsorted)[end]
    @test maximum(res, st) == last(Vsorted)[end]
    @test argmin(res, st) == first(Vsorted)[1]
    @test argmax(res, st) == last(Vsorted)[1]

    @show res
    @show st

    @info "collect:" => collect(idview(res, st))
    @info "========" => first.(Vsorted)

    @test collect(idview(res, st)) == first.(Vsorted)
    @test collect(distview(res, st)) == last.(Vsorted)

    pop!(Vsorted)
    _, st = pop!(res, st)
    @test collect(idview(res, st)) == first.(Vsorted)
    @test collect(distview(res, st)) == last.(Vsorted)
    @test length(res, st) == length(res)
    @test collect(res) == Vsorted

    popfirst!(Vsorted)
    _, st = popfirst!(res, st)

    @info "b   collect id:" => collect(idview(res, st))
    @info "b collect dist:" => collect(distview(res, st))
    @info "b ========" => first.(Vsorted)
    @info "b ========" => last.(Vsorted)
    @show res.id, res.dist, st
    @info collect(res)
    @info "c   collect id:" => collect(idview(res, st))
    @info "c collect dist:" => collect(distview(res, st))

    @test collect(idview(res, st)) == first.(Vsorted)
    @test collect(distview(res, st)) == last.(Vsorted)
    @test length(res, st) == length(res)
    @test collect(res) == Vsorted

end

function create_random_array(n, k)
    V = rand(Float32, n)
    Vsorted = sort!([i => v for (i, v) in enumerate(V)], by=x->x[2])[1:k]
    V, Vsorted
end

@testset "shifted vector-based result set" begin
    k = 10
    res = KnnResult(k)
    res2 = KnnResultShift(k)
    st = initialstate(res2)
    V, Vsorted = create_random_array(50, k)
    for i in eachindex(V)
        push!(res, i, V[i])
        st = push!(res2, st, i, V[i])
    end

    testsorted(res, copy(Vsorted))
    testsorted2(res2, st, Vsorted)
end

max_stars_repo_name: sadit/NNS

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/NNS


text:
using DNC
using Flux
using Zygote

using Flux.Data: DataLoader
include("repeatcopy.jl")

nbits = 4
minrepeats = 1
maxrepeats = 1
minlength = 1
maxlength = 1

X = nbits+2
Y = nbits+1
N, W, R = 16, 16, 2
controllerout = 64
clipvalue = 20

niter = 700
batchsize = 16
seqs = [RepeatCopy(;
            nbits=nbits,
            minrepeats=minrepeats,
            maxrepeats=maxrepeats,
            minlength=minlength,
            maxlength=maxlength)
        for i in 1:(niter*batchsize)]

batcheddata = RepeatCopyBatchLoader(seqs, batchsize=batchsize)

controller = Chain(Dense(DNC.inputsize(X, R, W), 64), Dense(64, 64), Dense(64, controllerout))

model = Dnc(controller, X, Y, controllerout, N, W, R, batchsize; clipvalue=clipvalue)

loss(rc::RepeatCopy; printoutput=false) = loss(model, rc; printoutput=printoutput)
loss(batch::Tuple; printoutput=false) = loss(model, batch...; printoutput=printoutput)

opt = RMSProp(1e-3)

evalcb = ThrottleIterations(100) do
    idx = rand(1:(length(seqs)-batchsize))
    loss(Base.iterate(batcheddata, idx)[1]; printoutput=true)
end

@time mytrain!(loss, Flux.params(model), batcheddata, opt; cb=evalcb)


text:
# Copyright 2021 The FirstOrderLp Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

mutable struct TwoSidedQpProblem
  "Lower bounds on the variables."
  variable_lower_bound::Vector{Float64}
  "Upper bounds on the variables."
  variable_upper_bound::Vector{Float64}
  "Lower bounds on the constraints."
  constraint_lower_bound::Vector{Float64}
  "Upper bounds on the constraints."
  constraint_upper_bound::Vector{Float64}
  "The constraint matrix."
  constraint_matrix::SparseMatrixCSC{Float64,Int64}
  "The constant term in the objective."
  objective_offset::Float64
  "The objective vector."
  objective_vector::Vector{Float64}
  "The objective matrix."
  objective_matrix::SparseMatrixCSC{Float64,Int64}
end

"""
Transforms a quadratic program into a `QuadraticProgrammingProblem` struct.
The `TwoSidedQpProblem` is destructively modified in place to avoid creating a
copy.

# Returns
A QuadraticProgrammingProblem struct containing the quadratic programming
problem.
"""
function transform_to_standard_form(qp::TwoSidedQpProblem)
  two_sided_rows_to_slacks(qp)

  is_equality_row = qp.constraint_lower_bound .== qp.constraint_upper_bound
  is_geq_row = .!is_equality_row .& isfinite.(qp.constraint_lower_bound)
  is_leq_row = .!is_equality_row .& isfinite.(qp.constraint_upper_bound)

  # Two-sided rows should be removed by two_sided_rows_to_slacks.
  @assert !any(is_geq_row .& is_leq_row)

  num_equalities = sum(is_equality_row)

  if num_equalities + sum(is_geq_row) + sum(is_leq_row) !=
     length(qp.constraint_lower_bound)
    error("Not all constraints have finite bounds on at least one side.")
  end

  # Flip the signs of the leq rows in place.
  # TODO: Skip this pass if there are no leq rows.
  for idx in 1:SparseArrays.nnz(qp.constraint_matrix)
    if is_leq_row[qp.constraint_matrix.rowval[idx]]
      qp.constraint_matrix.nzval[idx] *= -1
    end
  end

  new_row_to_old = [findall(is_equality_row); findall(.!is_equality_row)]
  if new_row_to_old != 1:size(qp.constraint_matrix, 1)
    row_permute_in_place(qp.constraint_matrix, invperm(new_row_to_old))
  end

  right_hand_side = copy(qp.constraint_lower_bound)
  right_hand_side[is_leq_row] .= .-qp.constraint_upper_bound[is_leq_row]
  permute!(right_hand_side, new_row_to_old)

  return QuadraticProgrammingProblem(
    qp.variable_lower_bound,
    qp.variable_upper_bound,
    qp.objective_matrix,
    qp.objective_vector,
    qp.objective_offset,
    qp.constraint_matrix,
    right_hand_side,
    num_equalities,
  )
end

"""
Transforms a `TwoSidedQpProblem` in-place to remove two-sided constraints, i.e.,
constraints with lower and upper bounds that are both finite and not equal.
For each constraint, a slack variable is added, and the constraint is changed
to an equality: `l <= a'x <= u` becomes `a'x - s = 0, l <= s <= u`.
"""
function two_sided_rows_to_slacks(qp::TwoSidedQpProblem)
  two_sided_rows = findall(
    isfinite.(qp.constraint_lower_bound) .&
    isfinite.(qp.constraint_upper_bound) .&
    (qp.constraint_lower_bound .!= qp.constraint_upper_bound),
  )
  if isempty(two_sided_rows)
    return
  end

  slack_matrix = sparse(
    two_sided_rows,  # row indices
    1:length(two_sided_rows),  # column indices
    fill(-1, length(two_sided_rows)),  # nonzeros
    length(qp.constraint_lower_bound),  # number of rows
    length(two_sided_rows),  # number of columns
  )
  qp.variable_lower_bound =
    [qp.variable_lower_bound; qp.constraint_lower_bound[two_sided_rows]]
  qp.variable_upper_bound =
    [qp.variable_upper_bound; qp.constraint_upper_bound[two_sided_rows]]
  qp.objective_vector = [qp.objective_vector; zeros(length(two_sided_rows))]
  qp.constraint_matrix = [qp.constraint_matrix slack_matrix]
  qp.constraint_lower_bound[two_sided_rows] .= 0
  qp.constraint_upper_bound[two_sided_rows] .= 0

  new_num_variables = length(qp.variable_lower_bound)
  row_indices, col_indices, nonzeros = SparseArrays.findnz(qp.objective_matrix)
  qp.objective_matrix = sparse(
    row_indices,
    col_indices,
    nonzeros,
    new_num_variables,
    new_num_variables,
  )
  return
end

"""
Reads an MPS or QPS file using the QPSReader package and transforms it into a
`QuadraticProgrammingProblem` struct.

# Arguments
- `filename::String`: the path of the file. The file extension is ignored,
  except that if the filename ends with ".gz", then it will be uncompressed
  using GZip. Accepted formats are documented in the README of the QPSReader
  package.
- `fixed_format::Bool`: If true, parse as a fixed-format file.

# Returns
A QuadraticProgrammingProblem struct.
"""
function qps_reader_to_standard_form(
  filename::String;
  fixed_format::Bool = false,
)
  if endswith(filename, ".gz")
    io = GZip.gzopen(filename)
  else
    io = open(filename)
  end

  format = fixed_format ? :fixed : :free

  mps = Logging.with_logger(Logging.NullLogger()) do
    QPSReader.readqps(io, mpsformat = format)
  end
  close(io)

  constraint_matrix =
    sparse(mps.arows, mps.acols, mps.avals, mps.ncon, mps.nvar)
  # The reader returns only the lower triangle of the objective matrix. We have
  # to symmetrize it.
  obj_row_index = Int[]
  obj_col_index = Int[]
  obj_value = Float64[]
  for (i, j, v) in zip(mps.qrows, mps.qcols, mps.qvals)
    push!(obj_row_index, i)
    push!(obj_col_index, j)
    push!(obj_value, v)
    if i != j
      push!(obj_row_index, j)
      push!(obj_col_index, i)
      push!(obj_value, v)
    end
  end
  objective_matrix =
    sparse(obj_row_index, obj_col_index, obj_value, mps.nvar, mps.nvar)
  @assert mps.objsense == :notset

  return transform_to_standard_form(
    TwoSidedQpProblem(
      mps.lvar,
      mps.uvar,
      mps.lcon,
      mps.ucon,
      constraint_matrix,
      mps.c0,
      mps.c,
      objective_matrix,
    ),
  )
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Libdl, LibGit2_jll

@testset "LibGit2_jll" begin
    major = Ref{Cint}(0)
    minor = Ref{Cint}(0)
    patch = Ref{Cint}(0)
    @test ccall((:git_libgit2_version, libgit2), Cint, (Ref{Cint}, Ref{Cint}, Ref{Cint}), major, minor, patch) == 0
    @test VersionNumber(major[], minor[], patch[]) == v"1.1.0"
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
module ShaleDrillingLikelihood_DataStructure_Test

using Revise
using ShaleDrillingLikelihood
using Test
using Random
using StatsBase
using InteractiveUtils
using StatsFuns

using ShaleDrillingLikelihood: SimulationDraws, _u, _v, SimulationDrawsMatrix, SimulationDrawsVector,
    ObservationRoyalty, DataRoyalty, _y, _x, _xbeta, num_choices, _num_x,
    DataProduce, _xsum, obs_ptr, group_ptr, _nu,
    groupstart, grouplength, grouprange, obsstart, obsrange, obslength,
    ObservationProduce, ObservationGroup,
    _i, _data, _num_obs, update_nu!, Observation, update!,
    SimulationDraw,
    _ψ1, _ψ2, _ρ

@testset "SimulationDraws" begin

    Random.seed!(1234)
    k = 3
    nobs = 10
    M = 5

    u = rand(M, nobs)
    v = rand(M, nobs)
    qm = zeros(M)
    x = SimulationDraws(u,v, similar(u), similar(v), qm, similar(qm), similar(qm), similar(qm), reshape(qm,1,:))

    @test size(x) == (M,nobs)
    @test isa(x, SimulationDrawsMatrix)
    vw = view(x, 1)
    @test isa(vw, SimulationDrawsVector)
    @test size(vw) == (M,)

    @test size(SimulationDraws(M, nobs, 1)) == (M,nobs,)

    prodsim(s::SimulationDraw) = _ψ1(s) * _ψ2(s)

    N = 100_000
    thetarho = 0.5
    xbar = mean(prodsim(SimulationDraw(randn(), randn(), thetarho)) for i in 1:N)
    xtrue = _ρ(thetarho)
    @test abs(xbar-xtrue) < 0.01
end

@testset "DataRoyalty" begin

    nobs = 10
    k = 3
    rates = [1/8, 3/16, 1/4]
    L = length(rates)

    X = rand(k,nobs)
    y = collect(1:nobs)
    @views sample!(collect(1:L), y[L+1:end])

    data = DataRoyalty(RoyaltyModel(), y, X, rates)

    @test length(data) == nobs == _num_obs(data)
    @test _num_x(data) == k
    @test_throws BoundsError data[nobs+1]
    @test_throws BoundsError data[0]
    @test iterate(data, nobs+1) == nothing
    @test iterate(data, 1) == iterate(data)
    let obsi = 0
        for (i,grp) in enumerate(data)
            obsi += 1
            @test Observation(grp,1) == Observation(data,i) # == first(grp) == last(grp)
        end
        @test obsi == length(data)
    end
end

@testset "DataProduce" begin

    @testset "DataProduce Deterministic" begin
        n = 10
        k = 2
        nwell = 3
        ngroups = 3
        nobs = n*nwell
        beta = rand(2)
        alphapsi = 0.5
        sigeta = 0.3
        sigu = 0.5
        theta = vcat(alphapsi, beta, sigeta, sigu)

        x = rand(k, nobs)
        y = x' * beta + rand(nobs)

        xsum = hcat(collect(sum(x[:, (j-1)*n .+ (1:n)]; dims=2) for j in 1:nwell)...)
        xpnu = similar(xsum)
        nusum = zeros(size(xsum,2))
        nusumsq = similar(nusum)

        nu = y - x'*beta
        obsptr = [j*n+1 for j in 0:nwell]
        groupptr = vcat(1,2,fill(nwell+1,ngroups+1-2))

        d = DataProduce(ProductionModel(),y,x,xsum,nu,xpnu,nusum,nusumsq,obsptr,groupptr)
        update!(d, theta)
        obs = Observation(d,2)
        @test isa(obs, ObservationProduce)

        @test length(d) == ngroups == length(groupptr)-1
        @test _num_x(d) == k
        @test length(d) == last(groupptr)-1 == length(obsptr)-1

        @test obsstart(d,2) == n+1
        @test obsstart(d,3) == 2*n + 1
        rng = obsstart(d,2) : obsstart(d,3)-1
        @test isa(rng, UnitRange)
        @test obsrange(d,2) == rng

        @test obslength(d,2) == length(rng)
        @test obsstart(d,2) == first(rng)

        @test length(obs) == n
        @test _y(obs) == y[rng]
        @test _x(obs) == x[:, rng]
        @test dropdims(sum(_x(obs); dims=2); dims=2) == _xsum(obs)

        grp = ObservationGroup(d,1)
        grouprange(grp)
        @test Observation(grp,1) == Observation(d,1)
        @test iterate(grp) == (Observation(d,1), 2,)
        @test iterate(grp, iterate(grp)[2]) == nothing

        grp = ObservationGroup(d,2)
        @test length(grp) == 2
        @test Observation(grp,1) == Observation(d,2)
        @test Observation(grp,2) == Observation(d,3)
        @test iterate(grp) == (Observation(d,2), 2,)
        @test iterate(grp, 2) == (Observation(d,3), 3,)
        @test iterate(grp, 3) == nothing

        grp = ObservationGroup(d,3)
        @test length(grp) == 0
        @test iterate(grp) == nothing

        ii, jj = (0,0)
        for g in d
            ii += 1
            for o in g
                jj+=1
            end
        end
        @test ii == length(d) == ngroups == length(groupptr)-1
        @test jj == _num_obs(d) == nwell
    end

    @testset "DataProduce random" begin
        numgroups = 11
        maxwells = 5
        mint = 4
        maxt = 25
        k = 3

        beta = rand(3)
        alphapsi = 0.5
        sigeta = 0.3
        sigu = 0.5

        # create group variables
        psi = rand(numgroups)
        grouplens = vcat(0, collect(0:maxwells)..., sample(0:maxwells, numgroups-maxwells-1))
        @test length(grouplens) == numgroups+1
        groups = 1 .+ cumsum(grouplens)

        @test all(groups[2:end] .- groups[1:end-1] .== grouplens[2:end])

        # create wells per group
        nwells = last(groups)-1
        us = rand(nwells)
        obslens = vcat(0, sample(mint:maxt, nwells))
        obsptr = 1 .+ cumsum(obslens)
        nobs = last(obsptr)-1

        @test all(obsptr[2:end] .- obsptr[1:end-1] .== obslens[2:end])

        # create observation epsilons
        e = rand(nobs)
        x = rand(k,nobs)
        xsum = zeros(k, nwells)
        xpnu = similar(xsum)
        nusum = zeros(size(xsum,2))
        nusumsq = similar(nusum)
        nu = zeros(nobs)
        y = x'*beta

        for i = 1:numgroups
            wellptrs = groups[i]:groups[i+1]-1
            for k in wellptrs
                obsrng = obsptr[k]:obsptr[k+1]-1
                # sum!(xsum[:,k:k], x[:, obsrng])
                # y[obsrng] .+= sigu .* us[k] .+ alphapsi .* psi[i]
            end
        end

        data = DataProduce(ProductionModel(),y,x,xsum,nu,xpnu,nusum,nusumsq,obsptr,groups)

        for g in data
            i = _i(g)
            for (k,o) in enumerate(g)
                j = getindex(grouprange(g), k)
                let y = _y(o)
                    y .+= sigu .* us[j] .+ alphapsi .* psi[i]
                end
                sum!(reshape(_xsum(o), :, 1), _x(o))
                _nu(o) .= _y(o) .- _x(o)'*beta
            end
        end

        for j = 1:nwells
            o = Observation(data,j)
            @test all(_xsum(o) .== sum(_x(o); dims=2))
        end

        @test _y(data) - _x(data)'*beta == _nu(data)

        @test !all(_nu(data) .== 0)

        dd = DataProduce(numgroups, maxwells, mint:maxt, vcat(rand(3), 0.5, 0.5, 0.3))

        # @code_warntype ShaleDrillingLikelihood.update_xsum!(dd)
        # @code_warntype ShaleDrillingLikelihood.update_xpnu!(dd)
        # @code_warntype update_nu!(dd, ProductionModel(), vcat(rand(3), 0.5, 0.5, 0.3))

    end
end

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule pandoc_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("pandoc")
JLLWrappers.@generate_main_file("pandoc", UUID("c5432543-76ad-5c9d-82bf-db097047a5e2"))
end  # module pandoc_jll


text:
function drive_forward(inputs::BehaviorInputs)
    run_continuous(inputs.robot.motors.right, 60)
    run_continuous(inputs.robot.motors.left, 60)
end

function turn_right(inputs::BehaviorInputs)
    run_continuous(inputs.robot.motors.right, -70)
    run_continuous(inputs.robot.motors.left, 30)
end

function halt(inputs::BehaviorInputs)
    stop(inputs.robot.motors.right, "brake")
    stop(inputs.robot.motors.left, "brake")
    stop(inputs.robot.head, "coast")
end

function look_right(inputs::BehaviorInputs)
    servo_absolute(inputs.robot.head, Degrees(3 * 80))
end

function look_left(inputs::BehaviorInputs)
    servo_absolute(inputs.robot.head, Degrees(-3 * 80))
end

const mapping_behaviors = Dict(:forward => drive_forward,
                               :turn_right => turn_right,
                               :look_right => look_right,
                               :look_left => look_left,
                               :halt => halt)

function mapping_state_machines(timeout=Second(30))
    driving_machine = StateMachine(:forward, :halt,
        [(:forward, (inputs) -> inputs.sensor_data.ultrasound < 0.25, :turn_right),
         (:turn_right, (inputs) -> inputs.sensor_data.ultrasound > 0.5, :forward),
         (:forward, (inputs) -> inputs.time > timeout, :halt),
         (:turn_right, (inputs) -> inputs.time > timeout, :halt)])

    head_machine = StateMachine(:look_right, :halt,
        [(:look_right, (inputs) -> inputs.sensor_data.head_angle < -pi/4, :look_left),
         (:look_left, (inputs) -> inputs.sensor_data.head_angle > pi/4, :look_right)])

    (driving_machine, head_machine)
end


text:
export search

abstract type AbstractDeductiveIndex end
abstract type AbstractDeductiveIndexTable end

"""
    Index(table::T, identifier_map::Dict{Int, E}) where {E, T}

Generic index container type which maps between index identifiers and associated values, as well as references the
root table of the index. In a linear index the `table` field contains the only index table, but in recursive indices
this table is merely the root of the index where the search starts.
"""
mutable struct Index{E, T <: AbstractDeductiveIndexTable} <: AbstractDeductiveIndex
    table::T
    identifier_map::Dict{Int, E}
    current_id::Int

    Index(table::T, identifier_map::Dict{Int, E}) where {E, T} = new{E, T}(table, identifier_map, 1)
end
"""
    add!(idx::Index{E, T}, el::E) where {E, T}

Adds an element to an index, making it searchable.
"""
function add!(idx::Index{E, T}, el::E) where {E, T}
    el_id = idx.current_id
    idx.identifier_map[el_id] = el
    idx.current_id += 1
    add!(idx.table, (el_id, el))
end

"""
    search(idx::Index{E, T}, pattern::E) where {E, T}

Searches an index for a given pattern, which is not necessarily an exact match. This function merely passes the search call to
the root index table, so details with how pattern matching is handled will be clarified by the `search` function documentation
for a specific index table type.
"""
search(idx::Index{E, T}, pattern::E) where {E, T} = Set{E}([idx.identifier_map[id] for id ∈ search(idx.table, pattern)])


text:
# the current julia logo and graphics

# not exported

const julia_blue    = (0.251, 0.388, 0.847)
const julia_purple  = (0.584, 0.345, 0.698)
const julia_green   = (0.22, 0.596, 0.149)
const julia_red     = (0.796, 0.235, 0.2)

# to be deprecated
const darker_blue    = (0.251, 0.388, 0.847)
const lighter_blue   = (0.4, 0.51, 0.878)
const darker_purple  = (0.584, 0.345, 0.698)
const lighter_purple = (0.667, 0.475, 0.757)
const darker_green   = (0.22, 0.596, 0.149)
const lighter_green  = (0.376, 0.678, 0.318)
const darker_red     = (0.796, 0.235, 0.2)
const lighter_red    = (0.835, 0.388, 0.361)
const purples        = (darker_purple, lighter_purple)
const greens         = (darker_green, lighter_green)
const reds           = (darker_red, lighter_red)

"""
    julialogo(;
        action=:fill,
        color=true,
        bodycolor=colorant"black")

Draw the Julia logo. The default action is to fill the logo and use the colors:

    julialogo()

If `color` is `false`, the `bodycolor` color is used for the logo..

The logo can be difficult to position well due to its asymmetric design.
The `0/0` point is at the top left corner, the total width is 315pt, the total height is 214pts. The optical center is _somewhere_ between 163-180pts in x, and 96-114pts in y. The gap to the left edge of "j"s descender is 16; the distance between the left edge of the "j" (not the descender) and the right edge of the "a" is at 268pts.

So, to place the logo by locating its center at a point, do this:

```
gsave()
# scale() first if required
translate(-165, -114) # anything between (x: -163 to -180, y: -96 to -114)
julialogo()
grestore()
```

To use the logo as a clipping mask:

    julialogo(action=:clip)

(In this case the `color` setting is automatically ignored.)

"""
function julialogo(;action=:fill,
        color=true,
        bodycolor=colorant"black")
    # save current color
    r, g, b, a = Luxor.get_current_redvalue(), Luxor.get_current_greenvalue(), Luxor.get_current_bluevalue(), Luxor.get_current_alpha()
    # "j"
    setcolor(bodycolor)
    move(72.872, 177.311)
    curve(72.872, 184.847, 72.024, 190.932, 70.329, 195.567)
    curve(68.633, 200.202, 66.222, 203.8, 63.094, 206.362)
    curve(59.967, 208.925, 56.217, 210.639, 51.846, 211.506)
    curve(47.476, 212.373, 42.615, 212.806, 37.264, 212.806)
    curve(30.029, 212.806, 24.49, 211.675, 20.647, 209.414)
    curve(16.803, 207.154, 14.881, 204.441, 14.881, 201.275)
    curve(14.881, 198.638, 15.955, 196.415, 18.103, 194.606)
    curve(20.251, 192.797, 23.134, 191.893, 26.751, 191.893)
    curve(29.464, 191.893, 31.631, 192.628, 33.251, 194.097)
    curve(34.872, 195.567, 36.209, 197.018, 37.264, 198.449)
    curve(38.47, 200.032, 39.487, 201.087, 40.316, 201.615)
    curve(41.145, 202.142, 41.899, 202.406, 42.577, 202.406)
    curve(44.009, 202.406, 45.102, 201.558, 45.855, 199.862)
    curve(46.609, 198.167, 46.985, 194.87, 46.985, 189.971)
    line(46.985, 97.051)
    line(72.872, 89.929)
    line(72.872, 177.311)
    newsubpath()

    "u"
    move(109.739, 92.416)
    line(109.739, 152.215)
    curve(109.739, 153.874, 110.06, 155.437, 110.7, 156.907)
    curve(111.341, 158.376, 112.226, 159.639, 113.357, 160.694)
    curve(114.487, 161.749, 115.806, 162.596, 117.313, 163.237)
    curve(118.821, 163.878, 120.441, 164.198, 122.174, 164.198)
    curve(124.134, 164.198, 126.36, 163.101, 129.069, 161.202)
    curve(133.363, 158.194, 135.965, 156.127, 135.965, 153.685)
    curve(135.965, 153.099, 135.965, 92.416, 135.965, 92.416)
    line(161.738, 92.416)
    line(161.738, 177.311)
    line(135.965, 177.311)
    line(135.965, 169.398)
    curve(132.574, 172.262, 128.956, 174.56, 125.113, 176.293)
    curve(121.27, 178.027, 117.539, 178.893, 113.922, 178.893)
    curve(109.702, 178.893, 105.783, 178.196, 102.165, 176.802)
    curve(98.548, 175.408, 95.383, 173.505, 92.67, 171.093)
    curve(89.957, 168.682, 87.828, 165.856, 86.283, 162.615)
    curve(84.739, 159.375, 83.966, 155.908, 83.966, 152.215)
    line(83.966, 92.416)
    line(109.739, 92.416)
    newsubpath()

    "l"
    move(197.881, 177.311)
    line(172.221, 177.311)
    line(172.221, 58.278)
    line(197.881, 51.156)
    line(197.881, 177.311)
    newsubpath()

    "i"
    move(208.603, 97.051)
    line(234.376, 89.929)
    line(234.376, 177.311)
    line(208.603, 177.311)
    line(208.603, 97.051)
    newsubpath()

    "a"
    move(288.225, 133.451)
    curve(285.738, 134.506, 283.232, 135.731, 280.707, 137.124)
    curve(278.183, 138.519, 275.884, 140.045, 273.812, 141.703)
    curve(271.74, 143.361, 270.062, 145.132, 268.782, 147.015)
    curve(267.501, 148.9, 266.86, 150.859, 266.86, 152.894)
    curve(266.86, 154.476, 267.067, 156.002, 267.481, 157.472)
    curve(267.896, 158.941, 268.48, 160.204, 269.234, 161.259)
    curve(269.988, 162.314, 270.816, 163.162, 271.721, 163.802)
    curve(272.625, 164.443, 273.605, 164.763, 274.66, 164.763)
    curve(276.77, 164.763, 278.899, 164.123, 281.047, 162.841)
    curve(283.194, 161.56, 285.587, 159.94, 288.225, 157.981)
    line(288.225, 133.451)
    newsubpath()
    move(314.111, 177.311)
    line(288.225, 177.311)
    line(288.225, 170.528)
    curve(286.792, 171.734, 285.399, 172.846, 284.042, 173.863)
    curve(282.686, 174.88, 281.16, 175.766, 279.464, 176.519)
    curve(277.768, 177.273, 275.866, 177.857, 273.755, 178.272)
    curve(271.645, 178.686, 269.158, 178.893, 266.295, 178.893)
    curve(262.375, 178.893, 258.853, 178.328, 255.725, 177.198)
    curve(252.597, 176.067, 249.941, 174.522, 247.756, 172.563)
    curve(245.571, 170.604, 243.893, 168.287, 242.725, 165.611)
    curve(241.558, 162.936, 240.973, 160.015, 240.973, 156.85)
    curve(240.973, 153.61, 241.595, 150.671, 242.838, 148.033)
    curve(244.082, 145.396, 245.778, 143.022, 247.925, 140.911)
    curve(250.073, 138.801, 252.579, 136.917, 255.443, 135.259)
    curve(258.306, 133.601, 261.377, 132.075, 264.655, 130.681)
    curve(267.934, 129.287, 271.344, 128.006, 274.886, 126.838)
    curve(278.427, 125.67, 281.932, 124.558, 285.399, 123.503)
    line(288.225, 122.825)
    line(288.225, 114.459)
    curve(288.225, 109.034, 287.188, 105.19, 285.116, 102.929)
    curve(283.044, 100.668, 280.274, 99.538, 276.807, 99.538)
    curve(272.738, 99.538, 269.912, 100.518, 268.329, 102.477)
    curve(266.747, 104.436, 265.955, 106.81, 265.955, 109.599)
    curve(265.955, 111.181, 265.786, 112.726, 265.447, 114.233)
    curve(265.108, 115.741, 264.523, 117.059, 263.695, 118.19)
    curve(262.866, 119.32, 261.679, 120.225, 260.134, 120.903)
    curve(258.589, 121.581, 256.649, 121.92, 254.312, 121.92)
    curve(250.695, 121.92, 247.756, 120.884, 245.495, 118.812)
    curve(243.234, 116.739, 242.104, 114.12, 242.104, 110.955)
    curve(242.104, 108.016, 243.102, 105.284, 245.099, 102.76)
    curve(247.097, 100.235, 249.79, 98.068, 253.182, 96.26)
    curve(256.573, 94.451, 260.492, 93.019, 264.938, 91.964)
    curve(269.384, 90.909, 274.094, 90.382, 279.068, 90.382)
    curve(285.173, 90.382, 290.429, 90.928, 294.838, 92.021)
    curve(299.246, 93.114, 302.883, 94.677, 305.746, 96.712)
    curve(308.609, 98.747, 310.72, 101.196, 312.076, 104.06)
    curve(313.433, 106.923, 314.111, 110.127, 314.111, 113.668)
    line(314.111, 177.311)
    (action == :clip) ? newsubpath() : do_action(action)

    color && setcolor(julia_red) #  red
    circle(Point(240.272, 68.091), Point(205.272, 68.091), :path)
    (action == :clip) ? newsubpath() : do_action(action)

    color && setcolor(julia_blue) # blue
    circle(Point(77.954, 68.091), Point(42.954, 68.091), :path)
    (action == :clip) ? newsubpath() : do_action(action)

    color && setcolor(julia_purple) # dark purple
    circle(Point(282.321, 68.091), Point(247.321, 68.091), :path)
    (action == :clip) ? newsubpath() : do_action(action)

    color && setcolor(julia_green) # green
    circle(Point(261.299, 31.672), Point(226.299, 31.672), :path)

    if action ≠ :clip
        do_action(action)
    else
        clip()
    end
    # restore saved color
    setcolor(r, g, b, a)
end

"""
    juliacircles(radius=100)

Draw the three Julia circles ("dots") in color centered at the origin.

The distance of the centers of each circle from the origin is `radius`.

The optional keyword argument `outercircleratio` (default 0.75) determines the radius of each circle relative to the main radius. So the default is to draw circles of radius 75 points around a larger circle of radius 100.

Return the three centerpoints.

The `innercircleratio` (default 0.65) no longer does anything useful (it used to draw the smaller circles) and will be deprecated.
"""
function juliacircles(radius=100; outercircleratio=0.75, innercircleratio=0.65)
    # clockwise, from bottom left
    color_sequence = [julia_red, julia_green, julia_purple]
    points = ngon(O, radius, 3, pi/6, vertices=true)
    @layer begin
        for (n, p) in enumerate(points)
            setcolor(color_sequence[n]...)
            circle(p, outercircleratio * radius, :fill)
        end
    end
    return points
end


text:
import Base: show

"""

# CmdStanModels

Base model specialized in:

### Method
```julia
*  SampleModel                  : StanSample.jl
*  OptimizeModel                : StanOptimize.jl
*  VariationalModel             : StanVariational.jl
*  DiagnoseModel                : StanDiagnose.jl
```
"""
abstract type CmdStanModels end

"""

# RandomSeed

Random number generator seed value

### Method
```julia
RandomSeed(;seed=-1)
```
### Optional arguments
```julia
* `seed::Int`           : Starting seed value
```
"""
mutable struct RandomSeed
  seed::Int64
end
RandomSeed(;seed::Number=-1) = RandomSeed(seed)

"""

# Init

Default bound for parameter initial value interval (if not found in init file)

### Method
```julia
Init(;bound=2)
```
### Optional arguments
```julia
* `bound::Number`           : Set interval to [-bound, bound]
```
"""
mutable struct Init
  bound::Int64
end
Init(;bound::Int64=2) = Init(bound)

"""

# Output

Default settings for output portion in cmdstan command

### Method
```julia
Output(;refresh=100)
```
### Optional arguments
```julia
* `refresh::Int64`                     : Output refresh rate
```
"""
mutable struct Output
  refresh::Int64
end
Output(;refresh=100) = Output(refresh)


text:
using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators, DomainSets, Distributions, NetCDF
using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())
u_exact = (x,t) -> exp.(-t) * cos.(x)
@parameters t x y
@variables u(..)
Dx = Differential(x)
Dy = Differential(y)
#Dxx = Differential(x)^2
#Dyy = Differential(y)^2
Dt = Differential(t)

#= ncinfo("/Users/Minwoo/Desktop/emis/emis_mole_all_20160701_cb6_bench.nc")
ncinfo("/Users/Minwoo/Desktop/emis/METDOT3D_160701.nc")
#101×81×35×25 Array{Float32, 4} /// 35 x 25???  35- layer =#

t_min= 0.
t_max = 24
x_min = 0.
x_max = 12000*101
#may be 101?
y_min = 0.
y_max = 12000*81

#Getting real dat from CMAQ
emisission_file = "/Users/Minwoo/Desktop/emis/METDOT3D_160701.nc"

nx_emis = ncgetatt(emisission_file, "Global", "NCOLS")
ny_emis = ncgetatt(emisission_file, "Global", "NROWS")
nt_emis = 24
dx_emis = ncgetatt(emisission_file, "Global", "XCELL")
dy_emis = ncgetatt(emisission_file, "Global", "YCELL")

#Getting Uhat data
r1 = ncread("/Users/Minwoo/Desktop/emis/METDOT3D_160701.nc","UHAT_JD")

function u_hat(x,y)
    i = Int(ceil((x - x_min)/dx_emis))
    j = Int(ceil((y - y_min)/dy_emis))
    #println(x, " ",y," ", i," ", j)
    u = r1[:, :, 1, 1]'
    return u[j,i]
end

@register u_hat(x, y)

#Getting Vhat data
r2 = ncread("/Users/Minwoo/Desktop/emis/METDOT3D_160701.nc","VHAT_JD")

function v_hat(x,y)
    i = Int(ceil((x - x_min)/dx_emis))
    j = Int(ceil((y - y_min)/dy_emis))
    #println(x, " ",y," ", i," ", j)
    v = r2[:, :, 1, 1]'
    return v[j,i]
end
@register v_hat(x, y)

#So2 Concetration
r3 = ncread("/Users/Minwoo/Desktop/emis/emis_mole_all_20160701_cb6_bench.nc","SO2")
emis = r3[:, :, 1, 1]'

function emissions(x,y)
    i = Int(ceil((x - x_min)/dx_emis))
    j = Int(ceil((y - y_min)/dy_emis))
    #println(x, " ",y," ", i," ", j)
    return emis[j,i]
end

#= function emissions(x,y,t)
    i = Int(ceil((x - x_min)/dx_emis))
    j = Int(ceil((y - y_min)/dy_emis))
    k = maximum([Int(ceil((t - t_min)/dt_emis)), 1])
    #println(x, " ",y," ",t," ", i," ", j," ",k)
    return emis[i,j,k]
end =#

@register emissions(x, y)

#@register emissions(x, y, t)

# 3D PDE

eq  = Dt(u(t,x,y)) ~ u_hat(x,y)*Dx(u(t,x,y)) + v_hat(x,y)*Dy(u(t,x,y)) + emissions(x, y)
bcs = #[u(t_min,x,y) ~ analytic_sol_func(t_min,x,y),
       [u(t_min,x,y) ~ 0,
       u(t,x_min,y) ~ 0,
       u(t,x_max,y) ~ 0,
       u(t,x,y_min) ~ 0,
       u(t,x,y_max) ~ 0]

# Space and time domains
domains = [t ∈ IntervalDomain(t_min,t_max),
           x ∈ IntervalDomain(x_min,x_max),
           y ∈ IntervalDomain(y_min,y_max)]
pdesys = PDESystem([eq],bcs,domains,[t,x,y],[u(t,x,y)])

# Method of lines discretization
dx = 24000.0; dy = 24000.0; dt = 1;
#dx = 12000
discretization = MOLFiniteDifference([x=>dx,y=>dy],t)
prob = ModelingToolkit.discretize(pdesys,discretization)
#sol = solve(prob,Tsit5())
#sol = solve(prob,saveat = dt, Tsit5())
sol = solve(prob,saveat = dt, progress = true, progress_steps = 1,Tsit5())

using Plots
_, xs,ys = [infimum(d.domain):dx:supremum(d.domain) for d in domains]

u_sol = reshape(sol.u[1], length(xs)-1,length(ys)-1)
p1 = plot(xs[2:length(xs)], ys[2:length(ys)], u_sol', linetype=:contourf,title = "solution")
u_sol = reshape(sol.u[2], length(xs)-1,length(ys)-1)
p2 = plot(xs[2:length(xs)], ys[2:length(ys)], u_sol', linetype=:contourf,title = "solution")
u_sol = reshape(sol.u[3], length(xs)-1,length(ys)-1)
p3 = plot(xs[2:length(xs)], ys[2:length(ys)], u_sol', linetype=:contourf,title = "solution")
u_sol = reshape(sol.u[4], length(xs)-1,length(ys)-1)
p4 = plot(xs[2:length(xs)], ys[2:length(ys)], u_sol', linetype=:contourf,title = "solution")
plottogether = plot(p1,p2,p3,p4)

julia> p1 = heatmap(xs[2:length(xs)-1], ys[2:length(ys)-1], u_sol', title = "solution")

max_stars_repo_name: EarthSciML/tutorials

max_issues_repo_name: EarthSciML/tutorials

max_forks_repo_name: EarthSciML/tutorials


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_m2(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_m2!(A_copy)
end

@inline function exp_ps_m2!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=4
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: A2 P1 C0 P0
    # Computing A2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing P1 = x*I+x*A+x*A2
    coeff1=0.5
    coeff2=0.16666666666666666
    coeff3=0.041666666666666664
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots2
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing C0 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating P1 in slot 3
    # Deallocating A2 in slot 2
    # Computing P0 = x*I+x*A+x*C0
    coeff1=1.0
    coeff2=1.0
    coeff3=1.0
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating C0 in slot 4
    return memslots1 # Returning P0
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using BinDeps

@BinDeps.setup

libmpfit = library_dependency("libmpfit", aliases=["mpfit"])

# build from source
provides(Sources, URI("https://www.physics.wisc.edu/~craigm/idl/down/cmpfit-1.3a.tar.gz"),
         libmpfit)

builddir = joinpath(BinDeps.depsdir(libmpfit), "src", "cmpfit-1.3a", "build")
libdir = BinDeps.libdir(libmpfit)

libfilename = ""
if Sys.isapple()
    libfilename = "libmpfit.dylib"
elseif Sys.isunix()
    libfilename = "libmpfit.so"
elseif Sys.iswindows()
    libfilename = "libmpfit.dll"
end

provides(SimpleBuild,
         (@build_steps begin
          FileRule(joinpath(libdir, libfilename),
                   @build_steps begin
                   GetSources(libmpfit)
                   CreateDirectory(builddir)
                   CreateDirectory(libdir)
                   @build_steps begin
                   ChangeDirectory(builddir)
                   `pwd`
                   `cp ../../../cmake.mpfit ../CMakeLists.txt`
                   `cmake ../`
                   `make`
                   `cp $libfilename $libdir`
                   end
                   end
                   )
          end), libmpfit)

BinDeps.@install Dict(:libmpfit => :libmpfit)


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dashtextareaautocomplete

"""
    dashtextareaautocomplete(;kwargs...)

A DashTextareaAutocomplete component.
DashTextareaAutocomplete

Simple `@webscopeio/react-textarea-autocomplete` wrapper for Dash
enabling auto-completion in multi-line `<textarea>` elements.
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `className` (String; optional): Class names of the <textarea> (from `react-textarea-autocomplete`).
- `containerClassName` (String; optional): Class names of the textarea container (from `react-textarea-autocomplete`).
- `containerStyle` (Dict; optional): Styles of the textarea container (from `react-textarea-autocomplete`).
- `dropdownStyle` (Dict; optional): Styles of the dropdown wrapper.
- `itemStyle` (Dict; optional): Styles of the items wrapper.
- `listStyle` (Dict; optional): Style of the list wrapper (from `react-textarea-autocomplete`).
- `loaderStyle` (Dict; optional): Style of the loader wrapper (from `react-textarea-autocomplete`).
- `minChar` (Real; optional): Number of characters that user should type for trigger a suggestion.
Defaults to 1. (from `react-textarea-autocomplete`)
- `placeholder` (String; optional): Provides a hint to the user of what can be entered in the <textarea> field.
- `style` (Dict; optional): Style of the <textarea>.(from `react-textarea-autocomplete`).
- `triggerChar` (String; optional): Character that triggers auto-completion machinery.
Defaults to `:`. (from `react-textarea-autocomplete`)
- `value` (String; optional): The value displayed in the <textarea>.
- `wordList` (Array; required): List of string available for auto-completion.
"""
function dashtextareaautocomplete(; kwargs...)
        available_props = Symbol[:id, :className, :containerClassName, :containerStyle, :dropdownStyle, :itemStyle, :listStyle, :loaderStyle, :minChar, :placeholder, :style, :triggerChar, :value, :wordList]
        wild_props = Symbol[]
        return Component("dashtextareaautocomplete", "DashTextareaAutocomplete", "dash_textarea_autocomplete", available_props, wild_props; kwargs...)
end


text:
"""
    report(res::OperationModelResults)

 function make_report(res::OperationModelResults; kwargs...)
 This function uses weave to either generate a LaTeX or HTML
 file based on the report_design.jmd (julia markdown) file
 that it reads. Out_path in the weave function dictates
 where the created file gets exported.

# Example
using Requirements
using PowerSimulations
using Weave

results = solve_op_model!(OpModel)
out_path = "/Users/lhanig/GitHub"

report(results, out_path)

# kwargs:
doctype = "md2html" to create an HTML, default is PDF via latex

jmd = "custom file path/report_design.jmd"
report(results, out_path; jmd = jmd)
jmd has a default of ".../pwd()/report_design/report_design.jmd"

"""

function report(res::OperationModelResults,out_path::String; kwargs...)

    doctype = get(kwargs, :doctype, "md2pdf")
    default_string = joinpath(pwd(), "src/utils/report_design/report_design.jmd")
    jmd = get(kwargs, :jmd, default_string)
    args = Dict("res" => res, "variables" => res.variables)
    Weave.weave(jmd, out_path=out_path, latex_cmd = "xelatex",
                doctype = doctype, args = args)

end


text:
maxpoolNd(x, poolSz) = maxpool(x, poolSz)
function maxpoolNd(x, poolSz::Tuple{<:Integer})
    shaped = reshape(x, (size(x, 1), 1, size(x)[2:end]...))
    subsampled = maxpool(shaped, (poolSz..., 1))
    ax = axes(subsampled)
    return subsampled[ax[1], 1, ax[3:end]...]
end

# extend maxpool to use rational types, in case pooling by 2x is too steep
import NNlib.maxpool

"""
    r = RationPool(resSize, k=2; nExtraDims=1, poolType = MeanPool)

slight extension of the Flux pooling methods to subsample at a rational rate
`resSize`. `k` is the window size over which to apply the pooling, and poolType
is the type of pooling, either MaxPool or MeanPool. nExtraDims counts the
number of dimensions uninvolved in the convolution; normally this is 2, the
last dimension for nExamples, and the penultimate for channels. You can expect
pooling to work on sizes up to 5 total dimensions.
"""
struct RationPool{A,B}
    m::A # should inheret from MaxPool
    resSize::B # a tuple of ints and rationals
end

function Base.show(io::IO, m::RationPool)
    print(io, "RationPool(windowSize=$(m.m.k), poolingRate=$(m.resSize))")
end

function RationPool(resSize::NTuple{N,Union{<:Integer,Rational{<:Integer}}}, k=2; nExtraDims=2, poolType=MeanPool) where {N}
    effResSize = (resSize..., ntuple(ii -> 1 // 1, min(nExtraDims - 2, 5))...)
    subBy = map(ki -> ((ki == 1) ? 1 : k), effResSize) # any non-trivial dim
    # should subsample at a rate of k
    m = poolType(subBy, pad=0, stride=1)
    # SamePad() means that sz = inputsize / stride
    RationPool{typeof(m),typeof(resSize)}(m, resSize)
end

RationPool(resSize::Union{<:Integer,Rational{<:Integer}}, k=3;
    nExtraDims=2) =
    RationPool((resSize,), k; nExtraDims=nExtraDims)

import Base: getindex
Base.getindex(X::RationPool, i::Union{AbstractArray,<:Integer}) = X.resSize[i]
import Base: ndims
ndims(r::MaxPool{N,M}) where {N,M} = N
ndims(r::MeanPool{N,M}) where {N,M} = N
nPoolDims(r::RationPool{A,<:NTuple{N,<:Any}}) where {A,N} = N
import Flux: outdims
outdims(r::RationPool) = 3
using Zygote: hook
function (r::RationPool)(x::AbstractArray{<:Any,N}) where {N}
    Nd = nPoolDims(r)
    Nneed = ndims(r.m) + 2
    if N <= Nneed
        # x needs to be slightly padded
        extraDims = ntuple(ii -> 1, Nneed - N)
        partial = reshape(x, (size(x)[1:Nd]..., extraDims..., size(x)[Nd+1:end]...))
    else
        # x is too big and needs to be temporarily reshaped
        extraDims = ntuple(ii -> 1, 2)
        partial = reshape(x, (size(x)[1:Nd]..., extraDims..., prod(size(x)[Nd+1:end])))
    end
    partial = r.m(partial)
    address = map(stopAtExactly_WithRate_FromSize_, size(partial)[1:nPoolDims(r)], r.resSize, size(x)[1:nPoolDims(r)])
    ax = axes(partial)
    if N <= Nneed
        endAxes = ax[Nd+Nneed-N+1:end] # grab the stuff after extraDims
        return partial[address..., extraDims..., endAxes...]
    else
        return reshape(partial[address..., extraDims..., ax[Nd+3:end]...], (length.(address)..., size(x)[Nd+1:end]...))
    end
end

"""
    poolSize(k, sizes)
if we're pooling at rates `k` in each dimension (e.g. `(3//2, 3//2)`), how many entries should we expect in the next layer
"""
function poolSize(r::RationPool, sizes)
    resSize = [x for x in r.resSize]
    [poolSingle(resSize[i], sizes[i]) for i in 1:length(resSize)]
end

function poolSize(kks, sizs)
    map(poolSingle, kks, sizs)
end

function poolSingle(kk, siz)
    return length(stopAtExactly_WithRate_(siz, kk))
end

function stopAtExactly_WithRate_(i, subBy)
    tmp = round.(Int, range(1, stop=i, length=round(Int, i / subBy)))
    return tmp
end
function stopAtExactly_WithRate_FromSize_(i, subBy, orig)
    tmp = round.(Int, range(1, stop=i, length=round(Int, orig / subBy)))
    return tmp
end


text:
using Distributions, LinearAlgebra, ForwardDiff, Seaborn, Pandas

# Evaluates target distribution π(q)
function π_q(q)
    q₁ = q[1]
    q₂ = q[2]
    in_bound = (-π/2 < q₁ < π/2) & (-π/2 < q₂ < π/2)
    return in_bound*(2/π * exp(-2*(q₁^2 + q₂^2)) + (sin(q₁*q₂)*sin(q₁)*cos(q₂))^2)
end

# Evaluates conditional probability π(p|q)
function π_pq(p,q)
    return pdf(MvNormal([0,0], 1.0*Matrix(I, 2, 2)), p)
end

# Evaluates Kinetic enery K(p,q)
function K(p,q)
    return -log(π_pq(p,q))
end

# Evaluates Potential energy V(q)
function V(q)
    return -log(π_q(q))
end

# Compute Hamiltonian at q,p
function H(p,q)
    return K(p,q) + V(q)
end

# Gets a trajectory of length L with n steps
function get_trajectory(q,p,L,n,dV_dq)

    # Step size
    ϵ = L/n

    q_trajectory = zeros(n,2)
    p_trajectory = zeros(n,2)

    q_trajectory[1,:] = q
    p_trajectory[1,:] = p

    for j=2:n
        ψ = p_trajectory[j-1,:] - ϵ/2 * dV_dq(q_trajectory[j-1,:])
        q_trajectory[j,:] = q_trajectory[j-1,:] + ϵ*ψ
        p_trajectory[j,:] = ψ - ϵ/2*dV_dq(q_trajectory[j,:])
    end

    return q_trajectory, p_trajectory
end

function hmc_alg(q₀, p₀)
    # Number of samples to return
    n = 25000

    q_samples = zeros(n,2)
    p_samples = zeros(n,2)
    q_samples[1,:] = q₀
    p_samples[1,:] = p₀

    trajectory_length = 1
    trajectory_partitions = 5

    # Gradient of V
    dV_dq(x) = ForwardDiff.gradient(V, x)

    for j=2:n

        q = q_samples[j-1,:]
        p = rand(MvNormal([0,0],1.0*Matrix(I, 2, 2)))
        q_proposal_trajectory, p_proposal_trajectory =
        get_trajectory(q, p, trajectory_length, trajectory_partitions, dV_dq)

        # Proposal
        q_L = q_proposal_trajectory[lastindex(q_proposal_trajectory)÷2,:]
        p_L = p_proposal_trajectory[lastindex(p_proposal_trajectory)÷2,:]

        # Compute acceptance probabilty
        α = min(1, exp(H(p,q) - H(-p_L,q_L)))

        if rand() < α
            q_samples[j,:] = q_L
            p_samples[j,:] = -p_L
        else
            q_samples[j,:] = q
            p_samples[j,:] = p
        end

    end

    return q_samples
end

samples = hmc_alg([0,0],[1,1])

using PyCall
using Plots
@pyimport seaborn as sns
pygui(true)
data = Pandas.DataFrame(Dict(:q₁=>samples[:,1], :q₂=>samples[:,2]))
rc("axes",facecolor = "#f7ece2")
rc("figure",facecolor = "#f7ece2")
rc("savefig",facecolor = "#f7ece2")
sns.jointplot(x = "q₁", y = "q₂", data = data, s = 1, alpha = 0.18)

using SpecialFunctions
q_values = -π/2:0.01:π/2

# Plot of q₁ marginal
function q₁_marginal(q)
    return √(2/π)*erf(π/√2)*exp(-2q^2) +
    sin(q)^2*(-π*q+π*q^3+sin(π*q))/(4*q^3-4*q)
end

q₁_values = [q₁_marginal(q) for q in q_values]

Plots.plot(q_values, q₁_values, legend = false, grid = false, xlabel = L"q_1", ylabel = L"\pi_{q_1}(q_1)", bg = RGB(247/255, 236/255, 226/255))

Plots.savefig("hmc_ex_2_marginal_1_plot.svg")

# Plot of q₂ marginal
function q₂_marginal(q)
    return √(2/π)*erf(π/√2)*exp(-2q^2) +
    1/4*cos(q)^2*(π + (1-2*q^2)*sin(π*q)/(q^3-q))
end

q₂_values = [q₂_marginal(q) for q in q_values]

Plots.plot(q_values, q₂_values, legend = false, grid = false, xlabel = L"q_2", ylabel = L"\pi_{q_2}(q_2)", bg = RGB(247/255, 236/255, 226/255))
Plots.savefig("hmc_ex_2_marginal_2_plot.svg")

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
include("../src/density.jl")

## Parameters for the spectral function plot
μ = 0.30            # Chemical potential
T = 0.0             # Temperature

nPts = 600;         # Number of points in the spectral function curve
x_pos = 100;  # Number of unit cells to be plotted in positive x  direction for neutral SV
# y_pos = 10;        #Number of unit cells to be plotted in positive y direction for both SVs

# Cut settings
xs = -x_pos:x_pos

#slice energy
ω = 0.05

# Coordinates of the UC's
XS = repeat(xs, 1, length(xs))
YS = permutedims(XS)

# Energies used in the spectral function
ω_min = -1;
ω_max = 10;
ωs = range(ω_min, ω_max, length = nPts)

# Plotting for different values of the potential
Delta_1 = LocalPotential(-1, Location(-0, 0))
Delta_2 = LocalPotential(-1.5, Location(0, 0))
Delta_3 = LocalPotential(-2.6, Location(1, 0))
#List of local potentials to simulate the charged and neutral SVs
POTENTIAL = [Delta_1]

s = AtomsSystem(μ, T, POTENTIAL)


text:
using Base.Threads
using StringDistances
using LightXML
using BenchmarkTools
using Random
using Plots; gr()

function normalize_str(c)
    c = replace(c, "\n" => "")
    c = replace(c, "\t" => "")
    c = replace(c, "  " => "")
    c = strip(c)
    c = uppercase(c)
    c
end

function distance_mat(c, d)
    n = length(c)
    n = length(c)
    m = Matrix(undef, n, n)

    @threads for i in 1:size(m)[1]
        @inbounds m[i, :] .= d.(c[i], c)
    end

    m
end

function centers_distance!(p, c, d, m)
    for k in 1:length(c)
        @inbounds m[:, k] .= d.(p[c[k]], p)
    end
end

function row_minimum(m)
    mindist = Array{Int64}(undef, size(m, 1))
    @threads for i in 1:size(m, 1)
        @inbounds mindist[i] = findfirst(m[i,:] .== minimum(m[i,:]))
    end
    mindist
end

function kmeans(strings, distfun, k, convtol = 1e-2, maxiter=25)
    centers = randperm(length(strings))[1:k]
    centers_mat = Matrix{Int64}(undef, length(strings), length(centers))

    converged = false
    iters = 0
    prior = 0
    distances = []
    while !converged | iters < maxiter
        iters += 1
        centers_distance!(strings, centers, d, centers_mat)
        assignments = row_minimum(centers_mat)

        cdist = 0
        for k in 1:length(centers)
            c_p = [strings[i] for i in 1:length(strings) if assignments[i] == k]
            if length(c_p) > 0
                dm = distance_mat(c_p, d)
                md = mean(dm, dims=2)
                new_center = findfirst(md[:] .== minimum(md))
                centers[k] = new_center
                cdist += sum(md)
            end
        end

        if iters > 1
            converged = abs(cdist/prior - 1) <= convtol
        end

        prior = cdist
        append!(distances, [cdist])
    end

    return centers, row_minimum(centers_mat), distances
end

files = readdir("blogs")
xposts = []

for f in files
    x = parse_file(joinpath("blogs", files[1]))
    r = root(x)
    append!(xposts, get_elements_by_tagname(r, "post"))
end

d = Levenshtein()
posts = normalize_str.(content.(xposts))
res = kmeans(posts, d, 3, 1e-5, 10)

max_stars_repo_name: Matthew-Amos/fuzzy-kmeans

max_issues_repo_name: Matthew-Amos/fuzzy-kmeans

max_forks_repo_name: Matthew-Amos/fuzzy-kmeans


text:
module OctreeBH
using StaticArrays

export Node, Data, AbstractData, GravTreeGather
export buildtree, get_scatter_ngb_tree, get_gather_ngb_tree
export gravity_treewalk!, gravity_bruteforce!
export nearest
export isLeaf
export get_max_tree_depth, set_max_depth_AMR!, balance_all_level!

include("particle_data.jl");
include("tree.jl");
include("ngb.jl");
include("gravity_treewalk.jl");
include("balance.jl")

end #module OctreeBH


text:
using DataStructures

mutable struct World
    universe::Array{Any,1} # Things to classify
    predicates::Array{Function,1} # Define classes
    reporting::Bool
    agenda::PriorityQueue{Any,Real,Base.Order.ReverseOrdering{Base.Order.ForwardOrdering}}
end

function explore_concepts(W::World)
    while !(isempty(W.agenda))
        current = dequeue!(W.agenda)
        if (W.reporting)
            println(W.agenda.xs)
            println("Calling: $(current[1])($(current[2]),W)")
        end
        #Call unqueued task
        current[1](current[2],W)
    end
end


text:
"""
# BBSPI

Julia implementation of Bit Banged
[SPI](https://en.wikipedia.org/wiki/Serial_Peripheral_Interface)

See https://github.com/notinaboat/PiADXL345.jl for usage example.

## TODO

 * Try putting the `delay` function in `struct SPISlave`.
   Does the delay loop still get inlined nicely?

"""
module BBSPI

function delay end

struct CPHA0 end
struct CPHA1 end

"""
    SPISlave(; cs=chip_select_output_pin,
              clk=clock_output_pin,
             mosi=master_output_pin,
             miso=master_input_pin)

Connect to SPISlave using GPIO pins.

Methods must be defined for `setindex!(::PinType, v)` and `getindex(::PinType)`.

e.g.

```
struct Pin
    pin::UInt8
end
Base.setindex!(p::Pin, v) = gpioWrite(p.pin, v)
Base.getindex(p::Pin) = gpioRead(p.pin)
```
"""
struct SPISlave{ClockPhase,CST,CLKT,MOSIT,MISOT}
    chip_select::CST
    clock::CLKT
    master_out::MOSIT
    master_in::MISOT

    function SPISlave{CPHA}(; cs::CST=nothing,
                             clk::CLKT=nothing,
                            mosi::MOSIT=nothing,
                            miso::MISOT=nothing,
                           ) where {CPHA,CST,CLKT,MOSIT,MISOT}

        new{CPHA,CST,CLKT,MOSIT,MISOT}(cs,clk,mosi,miso)
    end
end

SPISlave(; kwargs...) = SPISlave{CPHA1}(;kwargs...)

"""
    transfer(::SPISlave, tx_buffer, [rx_buffer])

Transfer bytes from `tx_buffer` to SPISlave.
Store response in `rx_buffer`.
"""
function transfer(s, tx, rx=UInt8[])::Nothing

    @assert isidle(s)

    rx_len = length(rx)
    tx_len = length(tx)

    s.chip_select[] = 1
    delay(s)

    for i in 1:max(tx_len, rx_len)
        byte = i <= tx_len ? tx[i] : UInt8(0)
        byte = transfer_byte(s, byte)
        if i <= rx_len
            rx[i] = byte
        end
    end

    s.master_out[] = 0
    s.chip_select[] = 0
    delay(s)

    @assert isidle(s)
    nothing
end

isidle(s) = s.chip_select[] == 0 && s.clock[] == 0

function write_bit(pin, tx)
    pin[] = (tx & UInt8(0x80) == 0) ? 0 : 1
    tx << 1
end

function read_bit(pin, rx)
    rx << 1 | pin[]
end

function transfer_bit(s::SPISlave{CPHA1}, tx, rx)

    s.clock[] = 1
    tx = write_bit(s.master_out, tx)
    delay(s)

    s.clock[] = 0
    rx = read_bit(s.master_in, rx)
    delay(s)

    tx, rx
end

function transfer_bit(s::SPISlave{CPHA0}, tx, rx)

    tx = write_bit(s.master_out, tx)
    delay(s)

    s.clock[] = 1
    rx = read_bit(s.master_in, rx)
    delay(s)

    s.clock[] = 0

    tx, rx
end

function transfer_byte(s, tx)

    rx = UInt8(0)
    for i in 1:8
        tx, rx = transfer_bit(s, tx, rx)
    end
    rx
end

end # module