text:
using Test
using PETSc, MPI
MPI.Initialized() || MPI.Init()

@testset "DMDACreate1D" begin
    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)
    for petsclib in PETSc.petsclibs
        PETSc.initialize(petsclib)
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt = PETSc.inttype(petsclib)
        # Loop over all boundary types and try to use them
        for boundary_type in instances(PETSc.DMBoundaryType)
            @testset "$boundary_type" begin
                dof_per_node = 4
                stencil_width = 5

                # We test both setting and not setting the point distribution
                points_per_proc = [PetscInt(10 + i) for i in 1:mpisize]
                proc_global_offsets =
                    [PetscInt(0), accumulate(+, points_per_proc)...]

                global_size = proc_global_offsets[end]

                # left and right ghost points
                gl =
                    boundary_type == PETSc.DM_BOUNDARY_NONE && mpirank == 0 ?
                    0 : stencil_width
                gr =
                    boundary_type == PETSc.DM_BOUNDARY_NONE &&
                    mpirank == mpisize - 1 ? 0 : stencil_width

                # Set the points
                da = PETSc.DMDACreate1d(
                    petsclib,
                    comm,
                    boundary_type,
                    global_size,
                    dof_per_node,
                    stencil_width,
                    points_per_proc,
                )

                @test PETSc.gettype(da) == "da"
                @test PETSc.getdimension(da) == 1

                da_info = PETSc.getinfo(da)
                corners = PETSc.getcorners(da)

                ghost_corners = PETSc.getghostcorners(da)
                @test da_info.dim == 1
                @test da_info.global_size == [global_size, 1, 1]
                @test da_info.procs_per_dim == [mpisize, 1, 1]
                @test da_info.boundary_type == [
                    boundary_type,
                    PETSc.DM_BOUNDARY_NONE,
                    PETSc.DM_BOUNDARY_NONE,
                ]
                @test da_info.stencil_type == PETSc.DMDA_STENCIL_BOX
                @test da_info.stencil_width == stencil_width
                @test corners.lower ==
                      [proc_global_offsets[mpirank + 1] + 1, 1, 1]
                @test corners.upper == [proc_global_offsets[mpirank + 2], 1, 1]
                @test corners.size == [points_per_proc[mpirank + 1], 1, 1]
                @test ghost_corners.lower ==
                      [proc_global_offsets[mpirank + 1] + 1 - gl, 1, 1]
                @test ghost_corners.upper ==
                      [proc_global_offsets[mpirank + 2] + gr, 1, 1]
                @test ghost_corners.size ==
                      [points_per_proc[mpirank + 1] + gl + gr, 1, 1]

                # Do not set the points and test option parsing
                da_refine = 2
                da = PETSc.DMDACreate1d(
                    petsclib,
                    comm,
                    boundary_type,
                    global_size,
                    dof_per_node,
                    stencil_width,
                    nothing;
                    da_refine = da_refine,
                )
                @test PETSc.gettype(da) == "da"
                @test PETSc.getdimension(da) == 1

                da_info = PETSc.getinfo(da)

                @test da_info.dim == 1
                if boundary_type == PETSc.DM_BOUNDARY_PERIODIC
                    @test da_info.global_size ==
                          [2^da_refine * global_size, 1, 1]
                else
                    @test da_info.global_size ==
                          [2^da_refine * (global_size - 1) + 1, 1, 1]
                end
                @test da_info.procs_per_dim == [mpisize, 1, 1]
                @test da_info.boundary_type == [
                    boundary_type,
                    PETSc.DM_BOUNDARY_NONE,
                    PETSc.DM_BOUNDARY_NONE,
                ]
                @test da_info.stencil_type == PETSc.DMDA_STENCIL_BOX
                @test da_info.stencil_width == stencil_width
                # In this case we cannot check the numbers locally

                # TODO: Need a better test?
                ksp = PETSc.KSP(da)
                @test PETSc.gettype(ksp) == "gmres"

            end
        end
        PETSc.finalize(petsclib)
    end
end

@testset "DMDACreate2D" begin
    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)
    global_size_x = 100
    global_size_y = 45
    for petsclib in PETSc.petsclibs
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt = PETSc.inttype(petsclib)
        PETSc.initialize(petsclib)
        # Loop over all boundary types and stencil types
        for stencil_type in instances(PETSc.DMDAStencilType),
            boundary_type_y in instances(PETSc.DMBoundaryType),
            boundary_type_x in instances(PETSc.DMBoundaryType)

            # skip unsupported stencils
            stencil_type == PETSc.DMDA_STENCIL_BOX &&
                (
                    boundary_type_x == PETSc.DM_BOUNDARY_MIRROR ||
                    boundary_type_y == PETSc.DM_BOUNDARY_MIRROR
                ) &&
                continue

            @testset "$boundary_type_x, $boundary_type_y, $stencil_type" begin
                dof_per_node = 4
                stencil_width = 5

                # Set the points
                da = PETSc.DMDACreate2d(
                    petsclib,
                    comm,
                    boundary_type_x,
                    boundary_type_y,
                    stencil_type,
                    global_size_x,
                    global_size_y,
                    PETSc.PETSC_DECIDE,
                    PETSc.PETSC_DECIDE,
                    dof_per_node,
                    stencil_width,
                    nothing,
                    nothing,
                )
                @test PETSc.gettype(da) == "da"
                @test PETSc.getdimension(da) == 2

                da_info = PETSc.getinfo(da)

                @test da_info.global_size == [global_size_x, global_size_y, 1]
                @test da_info.dim == 2
                @test prod(da_info.procs_per_dim) == mpisize
                @test da_info.boundary_type ==
                      [boundary_type_x, boundary_type_y, PETSc.DM_BOUNDARY_NONE]
                @test da_info.stencil_type == stencil_type
                @test da_info.stencil_width == stencil_width

                # test refinement
                da_refine = 2
                da = PETSc.DMDACreate2d(
                    petsclib,
                    comm,
                    boundary_type_x,
                    boundary_type_y,
                    stencil_type,
                    global_size_x,
                    global_size_y,
                    PETSc.PETSC_DECIDE,
                    PETSc.PETSC_DECIDE,
                    dof_per_node,
                    stencil_width,
                    nothing,
                    nothing;
                    da_refine = da_refine,
                )
                @test PETSc.gettype(da) == "da"
                @test PETSc.getdimension(da) == 2

                da_info = PETSc.getinfo(da)

                # Compute refined global size
                ref_global_size_x =
                    boundary_type_x == PETSc.DM_BOUNDARY_PERIODIC ?
                    2^da_refine * global_size_x :
                    2^da_refine * (global_size_x - 1) + 1
                ref_global_size_y =
                    boundary_type_y == PETSc.DM_BOUNDARY_PERIODIC ?
                    2^da_refine * global_size_y :
                    2^da_refine * (global_size_y - 1) + 1

                @test da_info.global_size ==
                      [ref_global_size_x, ref_global_size_y, 1]
                @test prod(da_info.procs_per_dim) == mpisize
                @test da_info.boundary_type ==
                      [boundary_type_x, boundary_type_y, PETSc.DM_BOUNDARY_NONE]
                @test da_info.stencil_type == stencil_type
                @test da_info.stencil_width == stencil_width

                # TODO: Test with specific distribution of processors and sizes

                # TODO: Need a better test?
                ksp = PETSc.KSP(da)
                @test PETSc.gettype(ksp) == "gmres"
            end
        end
        PETSc.finalize(petsclib)
    end
end

@testset "DMDACreate3D" begin
    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)
    global_size_x = 12
    global_size_y = 13
    global_size_z = 14
    for petsclib in PETSc.petsclibs
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt = PETSc.inttype(petsclib)
        PETSc.initialize(petsclib)
        # Loop over all boundary types and stencil types
        for stencil_type in instances(PETSc.DMDAStencilType),
            boundary_type_z in instances(PETSc.DMBoundaryType),
            boundary_type_y in instances(PETSc.DMBoundaryType),
            boundary_type_x in instances(PETSc.DMBoundaryType)

            stencil_type == PETSc.DMDA_STENCIL_BOX &&
                (
                    boundary_type_x == PETSc.DM_BOUNDARY_MIRROR ||
                    boundary_type_y == PETSc.DM_BOUNDARY_MIRROR ||
                    boundary_type_z == PETSc.DM_BOUNDARY_MIRROR
                ) &&
                continue

            @testset "$boundary_type_x, $boundary_type_y, $boundary_type_z, $stencil_type" begin
                dof_per_node = 4
                stencil_width = 2

                # Set the points
                da = PETSc.DMDACreate3d(
                    petsclib,
                    comm,
                    boundary_type_x,
                    boundary_type_y,
                    boundary_type_z,
                    stencil_type,
                    global_size_x,
                    global_size_y,
                    global_size_z,
                    PETSc.PETSC_DECIDE,
                    PETSc.PETSC_DECIDE,
                    PETSc.PETSC_DECIDE,
                    dof_per_node,
                    stencil_width,
                    nothing,
                    nothing,
                    nothing,
                )
                @test PETSc.gettype(da) == "da"
                @test PETSc.getdimension(da) == 3

                da_info = PETSc.getinfo(da)

                @test da_info.global_size ==
                      [global_size_x, global_size_y, global_size_z]
                @test da_info.dim == 3
                @test prod(da_info.procs_per_dim) == mpisize
                @test da_info.boundary_type ==
                      [boundary_type_x, boundary_type_y, boundary_type_z]
                @test da_info.stencil_type == stencil_type
                @test da_info.stencil_width == stencil_width

                # test refinement
                da_refine = 2
                da = PETSc.DMDACreate3d(
                    petsclib,
                    comm,
                    boundary_type_x,
                    boundary_type_y,
                    boundary_type_z,
                    stencil_type,
                    global_size_x,
                    global_size_y,
                    global_size_z,
                    PETSc.PETSC_DECIDE,
                    PETSc.PETSC_DECIDE,
                    PETSc.PETSC_DECIDE,
                    dof_per_node,
                    stencil_width,
                    nothing,
                    nothing,
                    nothing;
                    da_refine = da_refine,
                )
                @test PETSc.gettype(da) == "da"
                @test PETSc.getdimension(da) == 3

                da_info = PETSc.getinfo(da)

                # Compute refined global size
                ref_global_size_x =
                    boundary_type_x == PETSc.DM_BOUNDARY_PERIODIC ?
                    2^da_refine * global_size_x :
                    2^da_refine * (global_size_x - 1) + 1
                ref_global_size_y =
                    boundary_type_y == PETSc.DM_BOUNDARY_PERIODIC ?
                    2^da_refine * global_size_y :
                    2^da_refine * (global_size_y - 1) + 1
                ref_global_size_z =
                    boundary_type_z == PETSc.DM_BOUNDARY_PERIODIC ?
                    2^da_refine * global_size_z :
                    2^da_refine * (global_size_z - 1) + 1

                @test da_info.global_size ==
                      [ref_global_size_x, ref_global_size_y, ref_global_size_z]
                @test prod(da_info.procs_per_dim) == mpisize
                @test da_info.boundary_type ==
                      [boundary_type_x, boundary_type_y, boundary_type_z]
                @test da_info.stencil_type == stencil_type
                @test da_info.stencil_width == stencil_width

                # TODO: Test with specific distribution of processors and sizes

                # TODO: Need a better test?
                ksp = PETSc.KSP(da)
                @test PETSc.gettype(ksp) == "gmres"
            end
        end
        PETSc.finalize(petsclib)
    end
end

@testset "creatematrix" begin
    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)
    for petsclib in PETSc.petsclibs
        PETSc.initialize(petsclib)
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt = PETSc.inttype(petsclib)
        boundary_type = PETSc.DM_BOUNDARY_NONE
        dof_per_node = 1
        stencil_width = 1
        number_points = 10
        points_per_proc = [PetscInt(10) for i in 1:mpisize]
        global_size = sum(points_per_proc)
        # Set the points
        da = PETSc.DMDACreate1d(
            petsclib,
            comm,
            boundary_type,
            global_size,
            dof_per_node,
            stencil_width,
            points_per_proc,
        )
        mat = PETSc.creatematrix(da)

        # Build the 1-D Laplacian FD matrix
        Sten = PETSc.MatStencil{PetscInt}
        col = Vector{Sten}(undef, 2)
        row = Vector{Sten}(undef, 2)
        val = Vector{PetscScalar}(undef, 4)
        corners = PETSc.getcorners(da)

        for i in corners.lower[1]:min(corners.upper[1], global_size - 1)
            row[1] = Sten(i = i)
            row[2] = Sten(i = i + 1)
            col[1] = Sten(i = i)
            col[2] = Sten(i = i + 1)
            val .= [-1, 1, 1, -1]
            PETSc.MatSetValuesStencil!(mat, row, col, val, PETSc.ADD_VALUES)
        end

        PETSc.assemblybegin(mat)
        PETSc.assemblyend(mat)

        for i in corners.lower[1]:corners.upper[1]
            if i == 1
                @test mat[i, i:(i + 1)] == [-1, 1]
            elseif i == global_size
                @test mat[i, (i - 1):i] == [1, -1]
            else
                @test mat[i, (i - 1):(i + 1)] == [1, -2, 1]
            end
        end
        PETSc.finalize(petsclib)
    end
end

@testset "DM Vectors and Coordinates" begin
    comm = MPI.COMM_WORLD
    mpirank = MPI.Comm_rank(comm)
    mpisize = MPI.Comm_size(comm)
    for petsclib in PETSc.petsclibs
        PETSc.initialize(petsclib)
        PetscScalar = PETSc.scalartype(petsclib)
        PetscInt = PETSc.inttype(petsclib)
        boundary_type = PETSc.DM_BOUNDARY_NONE
        dof_per_node = 1
        stencil_width = 1
        number_points = 10
        points_per_proc = [PetscInt(10) for i in 1:mpisize]
        global_size = sum(points_per_proc)
        # Set the points
        da = PETSc.DMDACreate1d(
            petsclib,
            comm,
            boundary_type,
            global_size,
            dof_per_node,
            stencil_width,
            points_per_proc,
        )

        corners = PETSc.getcorners(da)

        # Create the local and global vectors
        local_vec = PETSc.createlocalvector(da)
        global_vec = PETSc.createglobalvector(da)

        # Fill everything with some data
        fill!(local_vec, mpirank)
        fill!(global_vec, mpisize)

        # Add the local values to the global values
        PETSc.update!(global_vec, local_vec, PETSc.ADD_VALUES)

        # end points added with neighbor due to ghost of size 1
        bot_val = mpisize + mpirank + (mpirank == 0 ? 0 : mpirank - 1)
        top_val = mpisize + mpirank + (mpirank == mpisize - 1 ? 0 : mpirank + 1)
        @test global_vec[corners.lower[1]] == bot_val
        @test global_vec[corners.upper[1]] == top_val

        # Center is just self plus the global
        for i in (corners.lower[1] + 1):(corners.upper[1] - 1)
            @test global_vec[i] == mpisize + mpirank
        end

        # reset the local values with the global values
        PETSc.update!(local_vec, global_vec, PETSc.INSERT_VALUES)

        # My first value and my ghost should be the bot/top values
        @test local_vec[1] == bot_val
        @test local_vec[2] == bot_val
        @test local_vec[end - 1] == top_val
        @test local_vec[end] == top_val

        # interior is just self plus the global
        for i in 3:(length(local_vec) - 2)
            @test local_vec[i] == mpisize + mpirank
        end

        # Test DM Coordinates
        coord_da = PETSc.getcoordinateDM(da)
        # Crank it up to 11!
        xmin, xmax = 0, 11
        PETSc.setuniformcoordinates!(coord_da, (xmin,), (xmax,))
        coord_vec = PETSc.getcoordinateslocal(coord_da)
        Δx = (xmax - xmin) / (global_size - 1)

        # Figure out the values we should have in the coordinate vector
        ghost_lower = corners.lower[1] - (mpirank == 0 ? 0 : 1)
        ghost_upper = corners.upper[1] + (mpirank == mpisize - 1 ? 0 : 1)
        for (loc, glo) in enumerate(ghost_lower:ghost_upper)
            @test coord_vec[loc] ≈ (glo - 1) * Δx
        end
        PETSc.finalize(petsclib)
    end
end

nothing


text:
module FDTD2D
export run

"""Calculate parameters needed for simulation and return them as a dictionary."""
function calculate_params(init_params)
  params = deepcopy(init_params)

  params["box_size"] = Dict{String,Float64}()
  params["box_size"]["x"] = params["x_bounds"][2] - params["x_bounds"][1]
  params["box_size"]["y"] = params["y_bounds"][2] - params["y_bounds"][1]

  params["space_step"] = Dict{String,Float64}()
  params["space_step"]["x"] = params["box_size"]["x"]/params["matrix_size"]["x"]
  params["space_step"]["y"] = params["box_size"]["y"]/params["matrix_size"]["y"]

  params["time_step"] = 0.5*sqrt(params["space_step"]["x"]^2 + params["space_step"]["y"]^2)

  params["half_cfl"] = Dict{String,Float64}()
  params["half_cfl"]["x"] = 0.5*params["time_step"]/params["space_step"]["x"]
  params["half_cfl"]["y"] = 0.5*params["time_step"]/params["space_step"]["y"]

  params["y1"] = range(params["y_bounds"][1], length=params["matrix_size"]["y"], step=params["space_step"]["y"]) |> collect
  params["y2"] = range(params["y_bounds"][1]+0.5*params["space_step"]["y"], length=params["matrix_size"]["y"]-1, step=params["space_step"]["y"]) |> collect

  return params
end

"""Initialize data"""
function init_data(params)
  data = Dict{String, Array{Float64,2}}()
  for name in ("ex","ey","ezx","ezy","ez")
    data[name] = zeros(Float64, params["matrix_size"]["x"], params["matrix_size"]["y"])
  end
  for name in ("hx","hy","hzx","hzy","hz")
    data[name] = zeros(Float64, params["matrix_size"]["x"]-1, params["matrix_size"]["y"]-1)
  end
  return data
end

"""Generate fields at left boundary (x_min = p["x_bounds"][0]) at time moment 'time'"""
function generate_fields_x_min!(d, time, p)
  d["ey"][2,:] -= p["laser_pulse_y_shape"].(time, p["x_bounds"][1] + p["space_step"]["x"], p["y1"])
  d["ezx"][2,:] -= p["laser_pulse_z_shape"].(time, p["x_bounds"][1] + p["space_step"]["x"], p["y1"])

  d["hzx"][2,:] -= p["laser_pulse_y_shape"].(time + 0.5*p["time_step"], p["x_bounds"][1] + 1.5*p["space_step"]["x"], p["y2"])
  d["hy"][2,:] -= p["laser_pulse_z_shape"].(time + 0.5*p["time_step"], p["x_bounds"][1] + 1.5*p["space_step"]["x"], p["y2"])

  d["ez"][2,:] = d["ezx"][2,:] + d["ezy"][2,:]
  d["hz"][2,:] = d["hzx"][2,:] + d["hzy"][2,:]
end

"""Update E_x"""
function update_ex!(ex, hz, v)
  for j = 2:size(ex,2)-1
    for i = 2:size(ex,1)-1
      ex[i,j] += v*(hz[i,j] + hz[i-1,j] - hz[i,j-1] - hz[i-1,j-1])
    end
  end
end

"""Update E_y"""
function update_ey!(ey, hz, v)
  for j = 2:size(ey,2)-1
    for i = 2:size(ey,1)-1
      ey[i,j] -= v*(hz[i,j] - hz[i-1,j] + hz[i,j-1] - hz[i-1,j-1])
    end
  end
end

"""Update H_x"""
function update_hx!(hx, ez, v)
  for j = 1:size(hx,2)
    for i = 1:size(hx,1)
      hx[i,j] -= v*(ez[i+1,j+1] + ez[i,j+1] - ez[i+1,j] - ez[i,j])
    end
  end
end

"""Update H_y"""
function update_hy!(hy, ez, v)
  for j = 1:size(hy,2)
    for i = 1:size(hy,1)
      hy[i,j] += v*(ez[i+1,j+1] - ez[i,j+1] + ez[i+1,j] - ez[i,j])
    end
  end
end

"""Make step"""
function make_step!(d, p)
  update_ex!(d["ex"], d["hz"], p["half_cfl"]["y"])
  update_ey!(d["ey"], d["hz"], p["half_cfl"]["x"])
  update_ey!(d["ezx"], d["hy"], -p["half_cfl"]["x"])
  update_ex!(d["ezy"], d["hx"], -p["half_cfl"]["y"])
  d["ez"] = d["ezx"] + d["ezy"]

  update_hx!(d["hx"], d["ez"], p["half_cfl"]["y"])
  update_hy!(d["hy"], d["ez"], p["half_cfl"]["x"])
  update_hy!(d["hzx"], d["ey"], -p["half_cfl"]["x"])
  update_hx!(d["hzy"], d["ex"], -p["half_cfl"]["y"])
  d["hz"] = d["hzx"] + d["hzy"]
end

function output(d, k, p)
  return true
end

"""Run the simulations"""
function run(init_params_dict)
  params_dict = calculate_params(init_params_dict)
  data = init_data(params_dict)
  for k = 1:params_dict["time_steps"]
    time = k*params_dict["time_step"]
    generate_fields_x_min!(data, time, params_dict)
    make_step!(data, params_dict)
    output(data, k, params_dict)
  end
  return data
end

end

ext: jl


text:
## --------------   Access methods

Base.show(io :: IO, e :: AbstractEntrySwitches) =
    print(io, typeof(e));
Base.show(io :: IO, e :: AbstractEntryDecision) =
    print(io, typeof(e),
        ":  preference scale ",  round(entry_pref_scale(e), digits = 2),
        "  ",  n_colleges(e), " colleges  ",  n_locations(e), " locations.");

min_entry_prob(switches :: AbstractEntrySwitches) = switches.minEntryProb;
max_entry_prob(switches :: AbstractEntrySwitches) = switches.maxEntryProb;
fix_entry_probs!(switches :: AbstractEntrySwitches) =
    switches.fixEntryProbs = true;
# Colleges that can only be attended locally
local_only_colleges(switches :: AbstractEntrySwitches) =
    switches.localOnlyIdxV;

"""
	$(SIGNATURES)

Fix preference scale parameter for entry decision.
"""
function fix_entry_pref_scale!(switches :: AbstractEntrySwitches{F1}, pScale) where F1
    switches.entryPrefScale = pScale;
    switches.calEntryPrefScale = false;
end

@forward AbstractEntryDecision.switches (
    min_entry_prob, max_entry_prob,
    capacities, capacities_c, capacity, limited_capacity,
    local_only_colleges, n_colleges, n_locations, n_types, value_local,
    fix_entry_probs!,
    type_mass_jl, type_mass_j
);

@forward AbstractEntryResults.switches (
    n_locations, n_types, n_colleges,
    capacities, capacities_c,
    type_mass_j
)

entry_probs_fixed(e :: AbstractEntryDecision) = e.switches.fixEntryProbs;

# Do not forward this
entry_pref_scale(e :: AbstractEntryDecision) = e.entryPrefScale;

"""
	$(SIGNATURES)

Mark which colleges can be attended by locals only.
"""
function set_local_only_colleges!(switches :: AbstractEntrySwitches, icV)
    if !isempty(icV)
        @assert all(icV .>= 1)  &&  all(icV .<= n_colleges(switches));
        switches.localOnlyIdxV = deepcopy(icV);
    end
    return nothing
end

"""
	$(SIGNATURES)

Number of locations.
"""
n_locations(switches :: AbstractEntrySwitches) = switches.nLocations;
# n_locations(e :: AbstractEntryResults) = n_locations(e.switches);
# n_locations(switches :: EntryDecisionSwitches{F1}) where F1 =
#     size(switches.typeMass_jlM, 2);

value_local(switches :: AbstractEntrySwitches{F1}) where F1 = switches.valueLocal;
# value_local(a :: EntryDecision{F1}) where F1 = a.valueLocal;

"""
	$(SIGNATURES)

Total mass of students (not conditional on entry).
"""
total_mass(switches :: AbstractEntrySwitches{F1}) where F1 =
    n_types(switches) * n_locations(switches);

"""
	$(SIGNATURES)

Number of student types.
"""
n_types(switches :: AbstractEntrySwitches{F1}) where F1 = switches.nTypes;

"""
	$(SIGNATURES)

Number of colleges.
"""
n_colleges(switches :: AbstractEntrySwitches{F1}) where F1 = switches.nColleges;

"""
	$(SIGNATURES)

Mass of each type. Only plays a role when colleges have capacities. Set to 1 otherwise. By location, if available.
"""
type_mass_jl(switches :: AbstractEntrySwitches{F1}, j :: Integer) where F1 =
    switches.typeMass_jlM[j, :];

type_mass_jl(switches :: AbstractEntrySwitches{F1}, j, l) where F1 =
    switches.typeMass_jlM[j, l];

# By location, if any
type_mass_jl(e :: AbstractEntryResults{F1}) where F1 = type_mass_jl(e.switches);
type_mass_jl(switches :: AbstractEntrySwitches{F1}) where F1 =
    switches.typeMass_jlM;

# type_mass_j(e :: AbstractEntryResults{F1}) where F1 = type_mass_j(e.switches);
type_mass_j(switches :: AbstractEntrySwitches{F1}) where F1 =
    vec(sum(switches.typeMass_jlM, dims = 2));

# type_mass_j(e :: AbstractEntryResults{F1}, j :: Integer) where F1 =
#     type_mass_j(e.switches, j);
type_mass_j(switches :: AbstractEntrySwitches{F1}, j :: Integer) where F1 =
    sum(switches.typeMass_jlM[j,:]);

"""
	$(SIGNATURES)

College capacities. Set to an arbitrary large number for entry mechanisms where capacities do not matter.
"""
capacities(switches :: AbstractEntrySwitches{F1}) where F1 =
    switches.capacity_clM;
    # capacities(e :: AbstractEntryResults{F1}) where F1 = capacities(e.switches);

"""
	$(SIGNATURES)

Set college capacities.
"""
function set_capacities!(switches :: AbstractEntrySwitches{F1},
    capacity_clM :: Matrix{F1}) where F1

    @assert size(capacity_clM) == size(switches.capacity_clM)
    switches.capacity_clM = capacity_clM;
    return nothing
end

"""
	$(SIGNATURES)

Increase one college capacities by a factor.
"""
increase_capacity!(switches :: AbstractEntrySwitches{F1},
    cFactor :: F1, cIdx :: Integer) where F1 =
    switches.capacity_clM[cIdx,:] .*= cFactor;

"""
	$(SIGNATURES)

Increase all college capacities by a common factor.
"""
increase_capacities!(switches :: AbstractEntrySwitches{F1}, cFactor) where F1 =
    switches.capacity_clM .*= cFactor;

"""
	$(SIGNATURES)

Capacities; summed across locations
"""
capacities_c(switches :: AbstractEntrySwitches{F1}) where F1 =
    vec(sum(switches.capacity_clM, dims = 2));
    # capacities_c(e :: AbstractEntryResults{F1}) where F1 = capacities_c(e.switches);

"""
	$(SIGNATURES)

Capacity of one college by location.
"""
capacity(switches :: AbstractEntrySwitches{F1}, iCollege :: Integer) where F1 =
    switches.capacity_clM[iCollege, :];

capacity(switches :: AbstractEntrySwitches{F1}, iCollege :: Integer,
    l :: Integer) where F1 = switches.capacity_clM[iCollege, l];

limited_capacity(switches:: AbstractEntrySwitches{F1}) where F1 =
    any(capacities(switches) .< CapacityInf);

"""
    $(SIGNATURES)

Change the number of types
Useful for solving with a subset of types, keeping in mind that
capacity constraints likely won't bind with fewer types.
"""
function subset_types!(switches :: EntryDecisionSwitches{F1}, typeV :: AbstractVector{I1}) where {F1, I1 <: Integer}

    @assert maximum(typeV) <= n_types(switches)
    switches.nTypes = length(typeV);
    switches.typeMass_jlM = switches.typeMass_jlM[typeV, :];
    @assert validate_es(switches)
    return nothing
end

"""
	$(SIGNATURES)

Adjust the `EntryDecision` so that only a subset of types is kept.
Useful for solving with a subset of types, keeping in mind that
capacity constraints likely won't bind with fewer types.
`deepcopy(entryS)` first to keep the original `EntryDecision` unchanged.
"""
subset_types!(entryS :: EntryDecision{F1}, typeV :: AbstractVector{I1}) where {F1, I1 <: Integer} =
    subset_types!(entryS.switches, typeV);

# --------------


text:
using PyCall
using DFTK

# Setup ASE quantum espresso calculator
ENV["ASE_ESPRESSO_COMMAND"] = "/usr/bin/pw.x -in PREFIX.pwi > PREFIX.pwo"
ENV["ESPRESSO_PSEUDO"] = joinpath(@__DIR__, "qe_pseudos")

const PSEUDOMAP = Dict(
   "Si" => "Si.pbe-hgh.UPF",
   "O"  => "O.pbe-hgh.UPF",
   "Al" => "Al.pbe-hgh.UPF",
   "Cu" => "Cu.pbe-d-hgh.UPF",
   "Li" => "Li.pz-s-hgh.UPF",
   "Ga" => "Ga.pbe-d-hgh.UPF",
   "As" => "As.pbe-hgh.UPF",
   "H" => "H.pbe-hgh.UPF",
)

"""Ecut is in Hartree, tol is energy tolerance in Hartree, kspacing is in inverse Bohrs"""
function run_qe(lattice, atoms; Ecut=20, tol=1e-11, mixing_mode="local-TF",
                fileprefix="espresso", kspacing=0.3, smearing="gaussian",
                temperature=0, mixing_beta=0.7, mixing_ndim=10,
                maxiter=500, kwargs...)
    espresso = pyimport("ase.calculators.espresso")

    kgrid = kgrid_size_from_minimal_spacing(lattice, kspacing)
    @info "Autoselected kgrid_size", kgrid
    system = ase_atoms(lattice, atoms)
    calc = espresso.Espresso(
        label=fileprefix,
        input_dft="pbe",
        pseudopotentials=PSEUDOMAP,
        kpts=kgrid,
        ecutwfc=Ecut / DFTK.units.Ry,  # QE uses Rydbergs
        tstress=true,                  # Compute stresses
        tprnfor=true,                  # Compute forces
        mixing_mode=mixing_mode,
        mixing_beta=mixing_beta,
        conv_thr=tol / DFTK.units.Ry,
        occupations=ifelse(temperature > 0, "smearing", "fixed"),
        smearing=smearing,
        degauss=temperature / DFTK.units.Ry,
        electron_maxstep=maxiter,
        mixing_ndim=mixing_ndim,
        # verbosity="high",
        kwargs...
    )
    system.calc = calc

    system.get_potential_energy() * DFTK.units.eV
end

function parse_pwo(fileprefix)
    fn = isfile(fileprefix) ? fileprefix : fileprefix * ".pwo"
    iterations = []
    current_iter = nothing

    for line in readlines(fn)
        miter = match(r"^ *iteration # *([0-9]+)", line)
        if miter !== nothing  # New iteration starts
            if !isnothing(current_iter)
                push!(iterations, (; current_iter...))
            end
            current_iter = Dict{Symbol, Any}(:n_iter => parse(Int, miter[1]))
        end

        menergy = match(r"^ *total energy *= *([0-9.eE+-]+) *Ry", line)
        if menergy !== nothing
            current_iter[:total_energy] = DFTK.units.Ry * parse(Float64, menergy[1])
        end

        maccuracy = match(r"^ *estimated scf accuracy *< *([0-9.eE+-]+) *Ry", line)
        if maccuracy !== nothing
            current_iter[:accuracy] = DFTK.units.Ry * parse(Float64, maccuracy[1])
        end
    end

    iterations
end


text:
test = "target area: x=20..30, y=-10..-5"
puzzle = "target area: x=253..280, y=-73..-46"

function parse_input(input::AbstractString)
    re = r"^target area: x=(-?[0-9]+)..(-?[0-9]+), y=(-?[0-9]+)..(-?[0-9]+)$"
    m = match(re, input)
    @assert m !== nothing
    coords = parse.(Int, m.captures)
    return ((coords[3], coords[4]), (coords[1], coords[2]))
end

function q1((ymin, ymax), (xmin, xmax))
    @assert ymax < 0
    vx = round(Int, - 1/2 + sqrt(1/4 + 2xmin), RoundUp)
    vy = -ymin - 1
    @assert vx * (vx - 1) / 2 < xmin <= vx * (vx + 1) / 2
    return vy * (vy + 1) ÷ 2
end

function q2((ymin, ymax), (xmin, xmax))
    @assert xmin > 0 && ymax < 0
    vxmin = round(Int, - 1/2 + sqrt(1/4 + 2xmin), RoundUp)
    vxmax = xmax
    @assert vxmin * (vxmin - 1) / 2 < xmin <= vxmin * (vxmin + 1) / 2
    vymin = ymin
    vymax = -ymin - 1

    return count(vymin:vymax .=> (vxmin:vxmax)') do v
        return shoot(v, (ymin, ymax), (xmin, xmax))
    end
end

function shoot((vy, vx), (ymin, ymax), (xmin, xmax))::Bool
    (x, y) = (0, 0)
    while x <= xmax && y >= ymin
        x >= xmin && y <= ymax && return true
        x += vx
        y += vy
        vx = max(vx - 1, 0)
        vy -= 1
    end
    return false
end

let v = parse_input(test)
    @assert q1(v...) == 45
    @assert q2(v...) == 112
end

v = parse_input(puzzle)
println("Q1: ", q1(v...))
println("Q2: ", q2(v...))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
using QuadraticTools
using Test

@testset "QuadraticTools.jl" begin
    # test the validity checker
    @test_throws DomainError QuadraticTools._validitycheck(0)
    @test_nowarn QuadraticTools._validitycheck(0.1)
    @test_nowarn QuadraticTools._validitycheck(0 + 5im)

    # test the delta calculation
    @test calcdelta(1, 1, -1) == 5
    @test calcdelta(1, 1, 1) == -3
    @test calcdelta(2, -7, -5) == 89

    # test conversions from the expanded form
    @test tofactored(Expanded(1, 5, 6)) == Factored(1, -3.0, -2.0)
    @test tovertex(Expanded(1, 5, 6)) == Vertex(1, -2.5, -0.25)
    @test calcvalue(Expanded(1, 5, 6), 10) == 156

    # test conversions from the factored form
    @test toexpanded(Factored(1, -3.0, -2.0)) == Expanded(1, 5.0, 6.0)
    @test tovertex(Factored(1, -3.0, -2.0)) == Vertex(1, -2.5, -0.25)
    @test calcvalue(Factored(1, -3.0, -2.0), 10) == 156.0
    @test tofactored(Vertex(-1, -2.5, -0.25)) == Factored(-1, -2.5 - 0.5im, -2.5 + 0.5im)

    # test conversions from the vertex form
    @test toexpanded(Vertex(1, -2.5, -0.25)) == Expanded(1, 5.0, 6.0)
    @test tofactored(Vertex(1, -2.5, -0.25)) == Factored(1, -3.0, -2.0)
    @test calcvalue(Vertex(1, -2.5, -0.25), 10) == 156.0
end


text:
export IdDict

"""
    IdDict{T}

A dictionary for converting key::T into integer id.

## 👉 Example
```julia
dict = IdDict{String}()
push!(dict, "abc") == 1
push!(dict, "def") == 2
push!(dict, "abc") == 1
dict["abc"] == 1

getkey(dict, id1) == "abc"

count(dict, id1) == 2
```
"""
type IdDict{T}
    key2id::Dict{T,Int}
    id2key::Vector{T}
    id2count::Vector{Int}

    IdDict() = new(Dict{T,Int}(), T[], Int[])
end
IdDict() = IdDict{Any}()

function IdDict{T}(data::Vector{T})
    d = IdDict{T}()
    for x in data
        push!(d, x)
    end
    d
end

function read(path)
    data = map(x -> chomp(x), open(readlines,path))
    IdDict(data)
end

Base.count(d::IdDict, id::Int) = d.id2count[id]

Base.getkey(d::IdDict, id::Int) = d.id2key[id]

Base.getindex(d::IdDict, key) = d.key2id[key]

Base.get(d::IdDict, key, default=0) = get(d.key2id, key, default)

Base.length(d::IdDict) = length(d.key2id)

function Base.push!(d::IdDict, key)
    if haskey(d.key2id, key)
        id = d.key2id[key]
        d.id2count[id] += 1
    else
        id = length(d.id2key) + 1
        d.key2id[key] = id
        push!(d.id2key, key)
        push!(d.id2count, 1)
    end
    id
end

Base.append!(d::IdDict, keys::Vector) = map(k -> push!(d,k), keys)


text:
include("spatial_types.jl")

struct Point2 <: AbstractPoint end
struct Point3 <: AbstractPoint end

function Point end


text:
#using Infiltrator

function fold(dict::Dict)
    io = IOBuffer()
    TOML.print(io, dict)
    return take!(io)
end

unfold(bytes::Vector{UInt8}) = TOML.parse(String(copy(bytes)))

struct Contract
    signatures::Vector
end

Contract(signatures::Vector{Vector{UInt8}}) = Dict[unfold(s) for s in signatures]

PeaceCypher.id(c::Contract, notary::Notary) = [id(s, notary) for s in c.signatures]

function PeaceCypher.verify(document, c::Contract, notary::Notary)
    for signature in c.signatures
        verify(document, signature, notary) || return false
    end
    return true
end

struct Mixer <: Agent
    mixer::ID
end

struct MetaData
    mix::ID
    gk::ID
    hash::Hash
end

function Base.Dict(m::MetaData)
    dict = Dict()
    dict["mix"] = string(m.mix, base=16)
    dict["gk"] = string(m.gk, base=16)
    dict["N"] = m.hash.N
    dict["hash"] = string(m.hash.hash, base=16)
    return dict
end

function MetaData(dict::Dict)

    mix = parse(BigInt, dict["mix"], base=16)
    gk = parse(BigInt, dict["gk"], base=16)
    N = dict["N"]
    hash = parse(BigInt, dict["hash"], base=16)

    return MetaData(ID(mix), ID(gk), Hash(N, hash))
end

struct SynchronicBraid <: Braid
    metadata::MetaData
    pseudonyms::Vector{ID}
end

function SynchronicBraid(ballot::Ballot)
    metadata = MetaData(unfold(ballot.metadata))
    pseudonyms = ID[]

    for i in 1:size(ballot.votes, 2)
        bytes = ballot.votes[:,i]
        str = String(copy(bytes))
        id = parse(BigInt, str, base=16)
        push!(pseudonyms, ID(id))
    end

    return SynchronicBraid(metadata, pseudonyms)
end

# struct SynchronicBraid{ID} <: Braid
#     braider::ID #Braider{ID}
#     mix::ID #Mixer{ID}
#     hash::Hash
#     pseudonyms::Vector{ID}
# end

pseudonyms(braid::SynchronicBraid) = braid.pseudonyms

function validate(state::State, braid::SynchronicBraid, signers::Vector)
#    @infiltrate
    # for id in braid.pseudonyms
    #     id in state.pseudonyms || return false
    # end

    length(unique(braid.pseudonyms)) == length(unique(signers)) == length(signers) || return false

    braid.metadata.gk in state.braiders || return false

    for signer in signers
        signer in state.pseudonyms || return false
    end

    return true
end

function Base.push!(tlog::BraidChain, row::Transaction{T}) where T <: Braid
    #signers = [id(s, tlog.crypto) for s in row.signatures]
    signers = id(row.signature, tlog.notary)

    @assert validate(tlog.state, row.document, signers)

    ### At this point we are also able to validate also that ledger is also validated by hashes

    @assert validate(tlog, row.document.metadata.hash)

    @assert verify(row.document, row.signature, tlog.notary)

    push!(tlog.state, row.document, signers)
    push!(tlog.transactions, row)
end


text:
using TidalFluxCalibrations, Base.Dates
using Base.Test

include("models.jl")


text:
import Base.parse, Base.range
abstract NWBData
abstract Acquisition <: NWBData
abstract Analysis <: NWBData
abstract Epochs <: NWBData
abstract General <: NWBData
abstract Processing <: NWBData
abstract Stimulus <: NWBData
abstract TimeSeries <: Acquisition
abstract AbstractElectricalSeries <: TimeSeries

type ElectricalSeries{T<:Unitful.Voltage,S<:Unitful.Time} <: AbstractElectricalSeries
	data::Array{T,2}
	help::String
	resolution::Float64
	electrode_idx::Array{Int64,1}
	name::String
  timestamps::Union{AbstractArray{S,1}, StepRange{S,S},Range{S}}
  rate::Unitful.Frequency{Float64}
  start_time::Unitful.Time{Float64}
end

function ElectricalSeries(data, help::String, resolution::Float64, electrode_idx::Array{Int64,1}, name::String, rate::Unitful.Frequency{Int64})
  timestamps = (0//1)u"s":(1//rate):(size(data,1)//rate)
  ElectricalSeries(data, help, resolution, electrode_idx, name, timestamps, float(rate),0.0u"s")
end

type SpikeEventSeries{T<:Unitful.Voltage, S<:Unitful.Time} <: AbstractElectricalSeries
	data::Array{T, 3}
	help::String
  name::String
	timestamps::Union{AbstractArray{S,1}, StepRange{S,S}}
end

function SpikeEventSeries{T<:Unitful.Voltage}(D::Array{T,2}, help, name, timestamps)
	wf = zeros(T, size(D,1), 1, size(D,2))
	wf[:,1,:] = D
	SpikeEventSeries(wf, help, name, timestamps)
end

type SpatialSeries{T<:Unitful.Length,S<:Unitful.Time} <: TimeSeries
	data::Array{T, 2}
	help::String
	name::String
	timestamps::Union{AbstractArray{S,1}, StepRange{S,S}}
end


text:
export solve

using LinearAlgebra: Symmetric

#
# Generic interface to solve a quadratic optimization problem of the form:
#
#  min  1/2 <Qx,x> + <q,x>
# l≤x≤u
#
solve(Q::Symmetric{<:Real},
      q::AbstractVector{<:Real},
      x_init::AbstractVector{<:Real},
      bc::BoundConstraints{<:Real,1},
      conf::Abstract_BC_QuadSolver_Conf) = @assert(false,"To implement")


text:
include("MaterialConstructor.jl")
include("FEMConstructor.jl")
include("DislocationConstructor.jl")
include("IntegrationConstructor.jl")

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
using ArgParse

include("./InsarTimeseries.jl")
using .InsarTimeseries

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--overwrite"
        action = :store_true
        help = "Erase previous processed stack files"
        "--zero"
        action = :store_true
        help = "Run just the igram zeroing step"
        "--igram-path"
        default = "."
        help = "Path containing .int files"
        "--geo-path"
        default = "../"
        help = "Path containing .geo files"
        "--ref-station"
        arg_type = String
        help =
            "Name of a gps reference station to shift the stack to." *
            " Will overwrite the current shifted dataset if specified"
        "--ref-row"
        arg_type = Int
        "--ref-col"
        arg_type = Int
        "--window"
        arg_type = Int
        range_tester = x -> (x >= 1)
        help = "Size of window to use to shift stack"
        default = 5
    end

    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    println("Parsed args:")
    for (arg, val) in parsed_args
        @show arg, val
    end

    overwrite = parsed_args["overwrite"]
    zero = parsed_args["zero"]
    igram_path = parsed_args["igram-path"]
    geo_path = parsed_args["geo-path"]
    ref_station = parsed_args["ref-station"]
    ref_row = parsed_args["ref-row"]
    ref_col = parsed_args["ref-col"]
    window = parsed_args["window"]

    if zero
        mask_filename = abspath(joinpath(igram_path, InsarTimeseries.MASK_FILENAME))
        @time InsarTimeseries.create_mask_stacks(
            igram_path,
            overwrite = overwrite,
            geo_path = geo_path,
            mask_filename = mask_filename,
        )
        zero_masked_areas(igram_path, input_exts = [".int", ".cc"], overwrite = overwrite)
    else
        @time InsarTimeseries.prepare_stacks(
            igram_path,
            overwrite = overwrite,
            geo_path = geo_path,
            ref_row = ref_row,
            ref_col = ref_col,
            ref_station = ref_station,
            window = window,
        )
    end

end

main()


text:
##################################################################
#  html to pdf conversion example
##################################################################

using Wkhtmltox

# Initialize conversion setup
pdf_init(0)

# Sets output parameters
ps = PdfSettings("out" => "c:/temp/example.pdf", # output file
                 "orientation" => "Landscape")   # orientation

ps["out"] = "c:/temp/example.pdf" # can be set also like this

ps["size.pageSize"]  # settings can be read like this

# Sets source parameters
os = PdfObject()

# 'page' is the source file
os["page"] = joinpath(dirname(@__FILE__),"../examples/example.html")

# Create the converter
conv = PdfConverter(ps, os) # pass the pdf settings and source settings

# push!(conv, os2) # additional sources can be added with 'push!()'

# Do conversion
run(conv)

# free objects
conv = nothing
pdf_deinit()


text:
module TestModularityBrim
using Test
using EcologicalNetworks
using LinearAlgebra

    N = web_of_life("M_PL_026")
    modules = lp(N)[2]
    m = brim(N, modules)

    @test Q(m...) ≈ 0.5 atol = 0.3
    @test Qr(m...) ≈ 0.5 atol = 0.3

    N = BipartiteProbabilisticNetwork(Matrix{Float64}(I, (50, 50)))
    modules = brim(lp(N)...)
    @test Q(modules...) > 0.95

    N = BipartiteNetwork(Matrix{Bool}(I, (50, 50)))
    modules = brim(lp(N)...)
    @test Q(modules...) > 0.95

end


text:
"""
    Schedule(ts, values, interpolations; unit = :epoch; kwargs...)

Describes how the values of a hyperparameter change over the
training.

## Examples

```julia
# one-cycle schedule over 10 epochs from 0.01 over 0.1 to 0.001
es = length(learner.data.training)
Schedule([0es, 2es, 10es], [0.01, 0.1, 0.001], [sineio(), sineio()])
```

`Schedule` is an alias for `Animations.Animation`, see the
[Animations.jl documentation](https://jkrumbiegel.github.io/Animations.jl/dev)
for more detailed information.
"""
const Schedule = Animation

"""
    Scheduler(schedules...)

Callback for hyperparameter scheduling. Takes pairs of [`HyperParameter`](#)
types and [`Schedule`](#)s.

See [the tutorial](/documents/docs/tutorials/hyperparameters.md) for more information.

## Example

```julia
es = length(learner.data.training)
lrschedule = Schedule([0es, 10es], [0.1, 0.001], Animations.sineio())
scheduler = Scheduler(
    LearningRate => lrschedule
)
```

See also [`Schedule`](#).
"""
mutable struct Scheduler <: Callback
    schedules::Dict{Type{<:HyperParameter},Animations.FiniteLengthAnimation}
    step::Int
    Scheduler(args...; kwargs...) = new(Dict(args...; kwargs...), 1)
end

Base.show(io::IO, scheduler::Scheduler) = print(
    io, "Scheduler(", join(keys(scheduler.schedules), ", "), ")")

function stateaccess(scheduler::Scheduler)
    # TODO: implement proper merging of permissions
    if length(keys(scheduler.schedules)) == 0
        hpstateaccess = (;)
    else
        hpstateaccess = merge(stateaccess.(keys(scheduler.schedules))...)
    end
    return (
        data = Read(), cbstate = (; hyperparams=Write(), history = Read()),
        hpstateaccess...
    )
end

function init!(::Scheduler, learner)
    if !haskey(learner.cbstate, :hyperparams)
        learner.cbstate.hyperparams = MVHistory()
    end
end

function on(::StepBegin, phase::AbstractTrainingPhase, scheduler::Scheduler, learner)
    step = scheduler.step
    for (H, animation) in scheduler.schedules
        value = Animations.at(animation, scheduler.step)
        sethyperparameter!(learner, H, value)
        push!(
            learner.cbstate.hyperparams,
            Symbol(H),
            learner.cbstate.history[phase].steps,
            value)
    end
    scheduler.step += 1
end

"""
    onecycle(nsteps, max_val, [start_val, end_val; pct_start])

Creates a one-cycle [`Schedule`](#) over `nsteps` steps from `start_val`
over `max_val` to `end_val`.

## Examples

```julia

epochlength = length(traindataiter)
cb = Scheduler(LearningRate => onecycle(10epochlength, 0.01))
learner = Learner(<args>..., cb)
```
"""
function onecycle(
        nsteps, max_val;
        pct_start=0.25,
        div=25, divfinal=1e5,
        start_val=max_val / div, end_val=max_val / divfinal)
    return Animation(
        [0, nsteps * pct_start, nsteps],
        [start_val, max_val, end_val],
        [Animations.sineio(), Animations.sineio()]
    )
end


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ c7b8a7d7-eda7-4cae-96c8-6c4e52d6e73d
using StructArrays, Query

# ╔═╡ 523d7e80-5c67-47a7-8ed7-69fa06b6b0f3
begin
	using XLSX, DataFrames
	using Statistics: mean
end

# ╔═╡ 9ae4f227-26cc-47c6-b95d-8f7811386171
#split a dataframe by case and dataset, then combine with summary stats
function case_dataset_pivot(summary, method)
	f = filter(row->row.method==method, summary)
	groups = groupby(f, [:case, :dataset])
	combine(groups, nrow => :n_problems,
		:total_elapsed_time => mean => :aver_run_time,
		:highest_rtol => mean => :aver_tol,
		:lowest_gap => mean => :aver_gap)
end

# ╔═╡ db530169-b2d3-48bf-a69a-49b4c070513b
begin
	# table helper functions
	tbl(s) = "
" * join(s) * "
"
	td(s) = "$(join(s))	"
	tr(s) = "
$(join(s))
"

	init_table(df) = fill("", (6*9, 6))
end

# ╔═╡ c9f0573c-ce4c-4114-8bee-5c52a3016246
# the combined dataset/case pivot table is a tall rectangle where every entry is a row
# we use this function to write each record to a specific location of a 2D summary
# table
""" write the passed dataframe to the record object, at the location specifed by df[:case] and df[:dataset]"""
function store!(record, df)
	str(i) = "$i"
	j = df[:case]
	i = ((df[:dataset] - 1) * 6) + 1

	record[i, j] = str(df[:case])
	record[i + 1, j] = str(df[:dataset])
	record[i + 2, j] = str(df[:n_problems])
	record[i + 3, j] = str(df[:aver_run_time])
	record[i + 4, j] = str(df[:aver_tol])
	record[i + 5, j] = str(df[:aver_gap])
end

# ╔═╡ e4f32a9c-34ff-4990-bd8b-2508eecd6e03
""" turn a combined GroupedDataframe into a flat table """
function flat_tbl(g)
	record = init_table(g)
	map(val -> store!(record, val), eachrow(g))
	record
end

# ╔═╡ 001d0435-b172-47c9-86e4-316a6a1104e0
show_tbl(g) = HTML(tbl(mapreduce(r -> tr(td.(r)), *, eachrow(flat_tbl(g)))))

# ╔═╡ 96a6702a-d181-45c7-b715-de9341eb31a6

# ╔═╡ 46b258b5-2e68-4739-a77b-a583afbe3dd6
"""turn an array into a head=>tail pair, used to add titles to the XLSX sheet when
making it a DataFrame"""
function col_to_pair(col)
	h, t = Iterators.peel(col)
	h => collect(t)
end

# ╔═╡ 8e2669a4-5382-45bc-b91e-2fe27b2a6910
#read the excel file we are summarizing
begin
	f = XLSX.readxlsx("../decision_tree_D.xlsx")
	sh = f["Sheet1"]
	dfm = DataFrame(map(col_to_pair, eachcol(DataFrame(sh[:]))))
	dfm
end

# ╔═╡ dde4a7e7-2298-440c-ae72-3546aaf8ef56
vasko_sum_table = case_dataset_pivot(dfm, "base");

# ╔═╡ 853f8eab-e829-423c-9a11-8153ff14090d
show_tbl(vasko_sum_table)

# ╔═╡ 8668715e-a505-4a59-9934-94d9ba824d95

# ╔═╡ 2ed608a9-b0f5-4c1d-85de-a0000f5009b2

# ╔═╡ 76515738-b71a-4d65-8583-3924d160813e

# ╔═╡ Cell order:
# ╠═853f8eab-e829-423c-9a11-8153ff14090d
# ╠═dde4a7e7-2298-440c-ae72-3546aaf8ef56
# ╠═9ae4f227-26cc-47c6-b95d-8f7811386171
# ╠═db530169-b2d3-48bf-a69a-49b4c070513b
# ╠═001d0435-b172-47c9-86e4-316a6a1104e0
# ╠═e4f32a9c-34ff-4990-bd8b-2508eecd6e03
# ╠═c9f0573c-ce4c-4114-8bee-5c52a3016246
# ╠═8e2669a4-5382-45bc-b91e-2fe27b2a6910
# ╠═96a6702a-d181-45c7-b715-de9341eb31a6
# ╠═46b258b5-2e68-4739-a77b-a583afbe3dd6
# ╠═c7b8a7d7-eda7-4cae-96c8-6c4e52d6e73d
# ╟─523d7e80-5c67-47a7-8ed7-69fa06b6b0f3
# ╠═8668715e-a505-4a59-9934-94d9ba824d95
# ╠═2ed608a9-b0f5-4c1d-85de-a0000f5009b2
# ╠═76515738-b71a-4d65-8583-3924d160813e

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
using Pkg
Pkg.activate("examples")
Pkg.instantiate()

import GPLinearODEMaker; GLOM = GPLinearODEMaker
using SymEngine

# process for se_kernel_base
@vars δ λ
GLOM.kernel_coder(GLOM.se_kernel_base(λ, δ), "se")

#process for matern52_kernel_base
@vars δ λ
GLOM.kernel_coder(GLOM.matern52_kernel_base(λ, δ), "m52")

@vars δ λ
GLOM.kernel_coder(GLOM.pp_kernel_base(λ, δ), "pp"; cutoff_var="λ")

@vars δ λ1 λ2 sratio
GLOM.kernel_coder(GLOM.matern52_kernel_base(λ1, δ) + sratio * sratio * GLOM.matern52_kernel_base(λ2, δ), "m52_m52")

@vars δ λ1 λ2 sratio
GLOM.kernel_coder(GLOM.se_kernel_base(λ1, δ) + sratio * sratio * GLOM.se_kernel_base(λ2, δ), "se_se")

# process for periodic_kernel_base
# @vars δ se_P λ
# GLOM.kernel_coder(periodic_kernel_base([se_P, λ], δ), "periodic")
# @vars δ λ
# GLOM.kernel_coder(se_kernel_base(λ, δ), "se"; periodic_var="δ")

# process for quasi_periodic_kernel_base
# @vars δ se_λ qp_P p_λ
# GLOM.kernel_coder(quasi_periodic_kernel_base([se_λ, qp_P, p_λ], δ), "quasi_periodic")
@vars δ δp se_λ p_amp
GLOM.kernel_coder(GLOM.se_kernel_base(se_λ, δ) * GLOM.se_kernel_base(1 / p_amp, δp), "qp"; periodic_var="δp")

#process for rq_kernel_base
@vars δ α μ
GLOM.kernel_coder(GLOM.rq_kernel_base([α, μ], δ), "rq")

@vars δ α μ
GLOM.kernel_coder(GLOM.rm52_kernel_base([α, μ], δ), "rm52")

# @vars δ P λ α
# GLOM.kernel_coder(periodic_rq_kernel_base([P, λ, α], δ), "periodic_rq_kernel")
@vars δ λ α
GLOM.kernel_coder(GLOM.rq_kernel_base([α, λ], δ), "rq_per"; periodic_var="δ")

# kernel_coder can use π_sym to know not to numerically evaluate π
@vars δ λ π_sym
custom_cos(λ, δ) = cos(2 * π_sym * δ / λ)
GLOM.kernel_coder(custom_cos(λ, δ), "cos")

@vars σ
GLOM.kernel_coder(σ^2, "scale")


text:
function hyperbolic_cross_weights(y::AbstractArray{T,1},nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2}) where {T <: AbstractFloat, S <: Integer}

  interpolation_matrix = zeros(size(nodes,1),size(nodes,1))

  n = 2*maximum(multi_index)+1

  unique_multi_index = sort(unique(multi_index))

  base_polynomials        = Array{Array{T,2},1}(undef,length(unique_multi_index))
  unique_base_polynomials = Array{Array{T,2},1}(undef,length(unique_multi_index))

  @inbounds for k in axes(nodes,1)

    # Construct the base polynomials

    for i in eachindex(unique_multi_index)
      base_polynomials[i] = chebyshev_polynomial(n-1,nodes[k,:])
    end

    # Compute the unique polynomial terms from the base polynomials

    for i in eachindex(unique_multi_index)
      if unique_multi_index[i] == 0
        unique_base_polynomials[i] = base_polynomials[i][:,1:1]
      else
        unique_base_polynomials[i] = base_polynomials[i][:,2*unique_multi_index[i]:2*unique_multi_index[i]+1]
      end
    end

    # Construct a row of the interplation matrix

    l = 1
    @inbounds @views for j in axes(multi_index,1)
      new_polynomials = unique_base_polynomials[multi_index[j,1]+1][1,:]
      for i = 2:size(multi_index,2)
        new_polynomials = kron(new_polynomials,unique_base_polynomials[multi_index[j,i]+1][i,:])
      end
      m = length(new_polynomials)
      interpolation_matrix[k,l:l+m-1] = new_polynomials
      l += m
    end

  end

  weights = interpolation_matrix\y

  return weights

end

function hyperbolic_cross_weights(y::AbstractArray{T,1},nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  # Normalize nodes to the [-1.0 1.0] interval

  d = size(multi_index,2)
  grid = similar(nodes)
  for i = 1:d
    grid[:,i] = normalize_node(nodes[:,i],domain[:,i])
  end

  weights = hyperbolic_cross_weights(y,grid,multi_index)

  return weights

end

function hyperbolic_cross_weights_threaded(y::AbstractArray{T,1},nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2}) where {T<:AbstractFloat,S<:Integer}

  interpolation_matrix = zeros(size(nodes,1),size(nodes,1))

  n = 2*maximum(multi_index)+1

  unique_multi_index = sort(unique(multi_index))

  @inbounds @sync @qthreads for k in axes(nodes,1)

    base_polynomials        = Array{Array{T,2},1}(undef,length(unique_multi_index))
    unique_base_polynomials = Array{Array{T,2},1}(undef,length(unique_multi_index))

    # Construct the base polynomials

    for i in eachindex(unique_multi_index)
      base_polynomials[i] = chebyshev_polynomial(n-1,nodes[k,:])
    end

    # Compute the unique polynomial terms from the base polynomials

    for i in eachindex(unique_multi_index)
      if unique_multi_index[i] == 0
        unique_base_polynomials[i] = base_polynomials[i][:,1:1]
      else
        unique_base_polynomials[i] = base_polynomials[i][:,2*unique_multi_index[i]:2*unique_multi_index[i]+1]
      end
    end

    # Construct a row of the interplation matrix

    l = 1
    @inbounds @views for j in axes(multi_index,1)
      new_polynomials = unique_base_polynomials[multi_index[j,1]+1][1,:]
      for i = 2:size(multi_index,2)
        new_polynomials = kron(new_polynomials,unique_base_polynomials[multi_index[j,i]+1][i,:])
      end
      m = length(new_polynomials)
      interpolation_matrix[k,l:l+m-1] = new_polynomials
      l += m
    end

  end

  weights = interpolation_matrix\y

  return weights

end

function hyperbolic_cross_weights_threaded(y::AbstractArray{T,1},nodes::Union{Array{T,1},Array{T,2}},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {T<:AbstractFloat,S<:Integer}

  # Normalize nodes to the [-1.0 1.0] interval

  d = size(multi_index,2)
  grid = similar(nodes)
  for i = 1:d
    grid[:,i] = normalize_node(nodes[:,i],domain[:,i])
  end

  weights = hyperbolic_cross_weights_threaded(y,grid,multi_index)

  return weights

end

function hyperbolic_cross_weights(y::AbstractArray{T,1},inverse_interpolation_matrix::Array{T,2}) where {T<:AbstractFloat}

  weights = inverse_interpolation_matrix*y

  return weights

end


text:
# ---
# title: 887. Super Egg Drop
# id: problem887
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Binary Search, Dynamic Programming
# link: <https://leetcode.com/problems/super-egg-drop/description/>
# hidden: true
# ---
#
# You are given `K` eggs, and you have access to a building with `N` floors from
# `1` to `N`.
#
# Each egg is identical in function, and if an egg breaks, you cannot drop it
# again.
#
# You know that there exists a floor `F` with `0 <= F <= N` such that any egg
# dropped at a floor higher than `F` will break, and any egg dropped at or below
# floor `F` will not break.
#
# Each _move_ , you may take an egg (if you have an unbroken one) and drop it
# from any floor `X` (with `1 <= X <= N`).
#
# Your goal is to know  **with certainty**  what the value of `F` is.
#
# What is the minimum number of moves that you need to know with certainty what
# `F` is, regardless of the initial value of `F`?
#
#
#
# **Example 1:**
#
#
#
#     Input: K = 1, N = 2
#     Output: 2
#     Explanation:
#     Drop the egg from floor 1.  If it breaks, we know with certainty that F = 0.
#     Otherwise, drop the egg from floor 2.  If it breaks, we know with certainty that F = 1.
#     If it didn't break, then we know with certainty F = 2.
#     Hence, we needed 2 moves in the worst case to know what F is with certainty.
#
#
# **Example 2:**
#
#
#
#     Input: K = 2, N = 6
#     Output: 3
#
#
# **Example 3:**
#
#
#
#     Input: K = 3, N = 14
#     Output: 4
#
#
#
#
# **Note:**
#
#   1. `1 <= K <= 100`
#   2. `1 <= N <= 10000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Copyright (c) 2018 Maarten Pronk, MIT license
# @references
# 1. `https://github.com/evetion/GeoArrays.jl/blob/master/src/geoarray.jl`
# Base.IndexStyle(::Type{T}) where {T<:AbstractGeoArray} = IndexLinear()
# Base.convert(::Type{Array{T, 3}}, A::AbstractGeoArray{T}) where {T} = convert(Array{T,3}, ga.A)
# using Stars

for f in (:isnan,)
    m = :Base
    @eval $m.$f(ga::AbstractGeoArray) = GeoArray(ga, vals=$m.$f.(ga.A))
end

Base.parent(ga::AbstractGeoArray) = ga.A
Base.iterate(ga::AbstractGeoArray) = iterate(ga.A)
Base.length(ga::AbstractGeoArray) = length(ga.A)
Base.size(ga::AbstractGeoArray) = size(ga.A)
Base.eltype(::Type{AbstractGeoArray{T}}) where {T} = T
Base.map(f, ga::AbstractGeoArray) = GeoArray(ga, vals=map(f, ga.A))

Base.show(io::IO, ::MIME"text/plain", ga::AbstractGeoArray) = show(io, ga)
function Base.show(io::IO, ga::AbstractGeoArray)
    info = gdalinfo(ga)

    crs = GeoFormatTypes.val(ga.crs)
    wkt = length(crs) == 0 ? "undefined CRS" : "CRS $crs"
    print(io, "size      : $(join(size(ga), "x")) $(typeof(ga.A)) \n")
    # print(io, "dimensions: $(info["dim"]) (nrow, ncol, nlyr)\n")
    print(io, "resolution: $(info["cellsize"]) (x, y)\n")
    print(io, "extent    : $(st_bbox(ga)) [xmin, ymin, xmax, ymax] \n")
    print(io, "affine    : $(ga.f) \n")
    print(io, "crs       : $(wkt) \n")
    print(io, "names     : $(ga.names) \n")
    print(io, "time      : $(ga.time) \n")
    print(io, "----------------\n")
    r_summary(ga.A[:, :, 1])
end

## OPERATIONS ------------------------------------------------------------------
"""
Check whether two `GeoArrays`s `a` and `b` are geographically equal,
although not necessarily in content.
"""
Base.:(==)(a::AbstractGeoArray, b::AbstractGeoArray) = size(a) == size(b) && a.f == b.f && a.crs == b.crs
Base.:(==)(a::AbstractGeoArray, b::Real) = GeoArray(a, vals=a.A .== b)
Base.:!(a::AbstractGeoArray) = GeoArray(a, vals=.!a.A)

Base_ops = ((:Base, :+), (:Base, :-), (:Base, :*), (:Base, :/))

Base_lgl = ((:Base, :>), (:Base, :<), (:Base, :>=), (:Base, :<=),
    (:Base, :!=),
    (:Base, :&), (:Base, :|))

for (m, f) in Base_ops
    # _f = Symbol(m, ".:", f)
    @eval begin
        $m.$f(a::AbstractGeoArray, b::AbstractGeoArray) = begin
            a == b || throw(
                DimensionMismatch("Can't operate on non-geographic-equal `AbstractGeoArray`s"))
            GeoArray(a, vals=$m.$f.(a.A, b.A))
        end

        $m.$f(a::AbstractGeoArray, b::Real) = GeoArray(a, vals=$m.$f.(a.A, b))
        $m.$f(a::Real, b::AbstractGeoArray) = GeoArray(a, vals=$m.$f.(a, b.A))
    end
end

# 对于lgl型运算，增加NaN的处理
for (m, f) in Base_lgl
    # _f = Symbol(m, ".:", f)
    @eval begin
        $m.$f(a::AbstractGeoArray, b::AbstractGeoArray) = begin
            a == b || throw(
                DimensionMismatch("Can't operate on non-geographic-equal `AbstractGeoArray`s"))
            GeoArray(a, vals=$m.$f.(a.A, b.A))
        end

        $m.$f(a::AbstractGeoArray, b::Real) = begin
            lgl = .&($m.$f.(a.A, b), .!isnan.(a.A))
            GeoArray(a, vals=lgl)
        end
        $m.$f(a::Real, b::AbstractGeoArray) = begin
            lgl = .&($m.$f.(a, b.A), .!isnan.(b.A))
            GeoArray(a, vals=lgl)
        end
    end
end

# export -, +, *, /, isless, equals, isequal, ==

## Math OPERATIONS -------------------------------------------------------------
for (m, f) in ((:Base, :sum), (:Base, :maximum), (:Base, :minimum), (:Statistics, :mean))
    @eval begin
        $m.$f(ga::AbstractGeoArray; dims=3, kw...) = begin
            arr = $m.$f(ga.A; dims=dims, kw...)
            GeoArray(ga, vals=arr)
        end
    end
end

export sum, mean, maximum, minimum


text:
if basename(pwd()) == "aoc"
    cd("2016/24")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
# Load packages.
using Conv
using FocusedBlindDecon
using Gadfly

# We consider an illustrative synthetic experiment with the following parameters.
ntg=30 # number of time samples in `g`
nr=20 # number of receivers
nt=40*ntg # time samples in records `d`
nts=nt # samples in `s`
#+

# The aim is to reconstruct the true `g`, i.e., `gobs`, that we are going to design below. This design
# is of particular interest e.g., in seismic inversion and room acoustics as they reveal
# the arrival of energy, propagated from an impulsive source, at the receivers.
# Here, the arrivals curve linearly and hyperbolically, depending on `c` and have onsets
# depending on `bfrac`. Their amplitudes are determined by `afrac`.
gobs=zeros(ntg,nr) # allocate
FBD.toy_direct_green!(gobs, c=4.0, bfrac=0.20, afrac=1.0); # add arrival 1
FBD.toy_reflec_green!(gobs, c=1.5, bfrac=0.35, afrac=-0.6); # add arrival 2
plotg=(x,args...)->spy(x, Guide.xlabel("channel"), Guide.ylabel("time"),args...) # define a plot recipe
p1=plotg(gobs, Guide.title("True g"))

# The source signature `s` for the experiment is arbitrary: we simply use a Gaussian random signal.
sobs=randn(nts)
plot(y=sobs,x=1:nts,Geom.line, Guide.title("arbitrary source"), Guide.xlabel("time"))

# The next task is to generate synthetic observed records `dobs`:
# first lets construct a linear operator `S`; then applying `S` on `g` will result in measurements `d`.
# This task can be skipped if measured `dobs` are already available.
S=Conv.S(sobs, gsize=[ntg,nr], dsize=[nt,nr], slags=[nts-1,0]);
dobs=reshape(S*vec(gobs), nt, nr);
#+

# We first need to allocate a parameter variable `pa`, where the inputs `gobs` and `sobs` are optional.
pa=P_fbd(ntg, nt, nr, nts, dobs=dobs, gobs=gobs, sobs=sobs)
#+

# The we perform LSBD i.e., least-squares fitting, without regularization, of `dobs` to jointly
# optimize the arrays `g` and `s`.
FBD.lsbd!(pa)
#+

# We extract `g` from `pa` and plot to notice that it doesn't match `gobs`.
p2=plotg(pa[:g], Guide.title("LSBD g"))

# Instead, we perform FBD that uses the focusing functionals to regularize `lsbd!`.
FBD.fbd!(pa)
#+

# Notice that the extract impulse responses are closer to `gobs`, except for a scaling factor and an overall translation in time.
p3=plotg(pa[:g], Guide.title("FBD g"))


text:
# currently only creates a symlink
function mklink(src::AbstractString, dest::AbstractString; soft=true, overwrite=true)
    if ispath(src)
        if ispath(dest) && soft && overwrite
            rm(dest, recursive=true)
        end

        if !ispath(dest)
            @unix_only begin
                run(`ln -s $(src) $(dest)`)
            end
            @windows_only begin
                if isfile(src)
                    run(`mklink $(dest) $(src)`)
                else
                    run(`mklink /D $(dest) $(src)`)
                end
            end
        elseif !soft
            error("$(dest) already exists.")
        end
    else
        error("$(src) is not a valid path")
    end
end

function copy(src::AbstractString, dest::AbstractString; soft=true, overwrite=true)
    if ispath(src)
        if ispath(dest) && soft && overwrite
            rm(dest, recursive=true)
        end

        if !ispath(dest)
            # Shell out to copy directories cause this isn't supported
            # in v0.3 and I don't feel like copying all of that code into this
            # project. This if should be deleted when 0.3 is deprecated
            cp(src, dest)
        elseif !soft
            error("$(dest) already exists.")
        end
    else
        error("$(src) is not a valid path")
    end
end

@doc doc"""
    We overload download for our tests in order to make sure we're just download. The
    julia builds once.
""" ->
function Base.download(src::AbstractString, dest::AbstractString, overwrite)
    if !ispath(dest) || overwrite
        download(src, dest)
    end
    return dest
end

function get_playground_dir(config::Config, dir::AbstractString, name::AbstractString)
    if dir == "" && name == ""
        return abspath(joinpath(pwd(), config.default_playground_path))
    elseif dir == "" && name != ""
        return abspath(joinpath(config.dir.store, name))
    elseif dir != ""
        return abspath(dir)
    end
end

function get_playground_name(config::Config, dir::AbstractString)
    root_path = abspath(dir)
    name = ""

    for p in readdir(config.dir.store)
        file_path = joinpath(config.dir.store, p)
        if islink(file_path)
            if abspath(readlink(file_path)) == root_path
                name = p
                break
            end
        end
    end

    return name
end

function julia_url(version::VersionNumber, os::Symbol=OS_NAME, arch::Integer=WORD_SIZE)
    # Cannibalized from https://github.com/travis-ci/travis-build/blob/master/lib/travis/build/script/julia.rb

    if os === :Linux && arch == 64
        os_arch = "linux/x64"
        ext = "linux-x86_64.tar.gz"
        nightly_ext = "linux64.tar.gz"
    elseif os === :Linux && arch == 32
        os_arch = "linux/x32"
        ext = "linux-i686.tar.gz"
        nightly_ext = "linux32.tar.gz"
    elseif os === :Darwin && arch == 64
        os_arch = "osx/x64"
        ext = "osx10.7+.dmg"
        nightly_ext = "osx.dmg"
    elseif os === :Windows && arch == 64
        os_arch = "winnt/x64"
        ext = "win64.exe"
        nightly_ext = ext
    elseif os === :Windows && arch == 32
        os_arch = "winnt/x86"
        ext = "win32.exe"
        nightly_ext = ext
    else
        error("Julia does not support $arch-bit $os")
    end

    # Note: We could probably get specific revisions if we really wanted to.
    # eg. https://status.julialang.org/download/osx10.7+ redirects to
    # https://s3.amazonaws.com/julianightlies/bin/osx/x64/0.5/julia-0.5.0-2bb94d6f99-osx.dmg

    major_minor = "$(version.major).$(version.minor)"
    future_release = Base.nextpatch(NIGHTLY)  # Note: Nightly expected to be a prerelease
    if version >= future_release
        throw(ArgumentError("The version $version exceeds the latest known nightly build $NIGHTLY"))
    elseif version >= NIGHTLY
        url = "s3.amazonaws.com/julianightlies/bin/$os_arch/julia-latest-$nightly_ext"
    elseif version.patch == 0 && version == Base.upperbound(version)
        url = "s3.amazonaws.com/julialang/bin/$os_arch/$major_minor/julia-$major_minor-latest-$ext"
    else
        url = "s3.amazonaws.com/julialang/bin/$os_arch/$major_minor/julia-$version-$ext"
    end

    return "https://$url"
end

function julia_url(version::AbstractString, os::Symbol=OS_NAME, arch::Integer=WORD_SIZE)
    if version == "nightly"
        ver = NIGHTLY
    elseif version == "release"
        latest_release = VersionNumber(NIGHTLY.major, NIGHTLY.minor - 1, 0, (), ("",))
        ver = latest_release
    else
        ver = VersionNumber(version)
    end

    return julia_url(ver, os, arch)
end


text:
#=
Author: Miles Cochran-Branson
Date: 2/15/22
=#

module LUComptonScatter

#dependencies
import LinearAlgebra
import DataFrames
import Findpeaks
import Smoothers
import Interpolations

#functions for peak finding and analysis
include("analysisFuncs.jl")
export expected_energy
export read_Txt
export get_xy_data
export select_part
export weighted_mean
export find_peak_ends
export peak_parameters
export get_peak_value
export gauss
export fit_to_gauss

#functions for data fitting
include("chiSquredFitting.jl")
export χ²
export compton
export scan
export scan_box

end # module


text:
module TestMultivariateStats

using Test
using MLJBase
import MultivariateStats, Random
using MLJModels.MultivariateStats_
import Dates

using LinearAlgebra
import Random.seed!
seed!(1234)

@testset "Ridge" begin
    ## SYNTHETIC DATA TEST

    # Define some linear, noise-free, synthetic data:
    intercept = -42.0
    coefficients = Float64[1, 3, 7]
    n = 1000
    A = randn(n, 3)
    Xtable = table(A)
    y = A*coefficients

    # Train model on all data with no regularization and no
    # standardization of target:
    ridge = RidgeRegressor(lambda=0.0)

    fitresult, report, cache = fit(ridge, 0, Xtable, y)

    # Training error:
    yhat = predict(ridge, fitresult, Xtable)
    @test norm(yhat - y)/sqrt(n) < 1e-12

    # Get the true intercept?
    fr = fitted_params(ridge, fitresult)
    @test abs(fr.intercept) < 1e-10
    @test norm(fr.coefficients - coefficients) < 1e-10

    info_dict(ridge)
end

X, y = @load_crabs

@testset "PCA" begin
    X_array = matrix(X)
    pratio = 0.9999

    # MultivariateStats PCA
    pca_ms = MultivariateStats.fit(MultivariateStats.PCA,
                                   permutedims(X_array),
                                   pratio=pratio)
    Xtr_ms = permutedims(MultivariateStats.transform(pca_ms,
                                                     permutedims(X_array)))

    # MLJ PCA
    pca_mlj = PCA(pratio=pratio)
    fitresult, _, _ = fit(pca_mlj, 1, X)
    Xtr_mlj = matrix(transform(pca_mlj, fitresult, X))

    @test Xtr_mlj ≈ Xtr_ms
end

@testset "KernelPCA" begin
    X_array = matrix(X)

    # MultivariateStats KernelPCA
    kpca_ms = MultivariateStats.fit(MultivariateStats.KernelPCA
                                  , permutedims(X_array))
    Xtr_ms = permutedims(MultivariateStats.transform(kpca_ms,
                                                     permutedims(X_array)))

    # MLJ KernelPCA
    kpca_mlj = KernelPCA()
    fitresult, _, _ = fit(kpca_mlj, 1, X)
    Xtr_mlj = matrix(transform(kpca_mlj, fitresult, X))

    @test Xtr_mlj ≈ Xtr_ms
end

@testset "ICA" begin
    X_array = matrix(X)
    k = 5
    tolerance = 5.0

    # MultivariateStats ICA
    seed!(1234) # winit gets randomly initialised
    ica_ms = MultivariateStats.fit(MultivariateStats.ICA,
                                   permutedims(X_array),
                                   k;
                                   tol=tolerance)
    Xtr_ms = permutedims(MultivariateStats.transform(ica_ms,
                                                     permutedims(X_array)))

    # MLJ ICA
    seed!(1234) # winit gets randomly initialised
    ica_mlj = ICA(;k=k, tol=tolerance)
    fitresult, _, _ = fit(ica_mlj, 1, X)
    Xtr_mlj = matrix(transform(ica_mlj, fitresult, X))

    @test Xtr_mlj ≈ Xtr_ms
end

@testset "MulticlassLDA" begin
    Xfull, y = @load_smarket
    X = selectcols(Xfull, [:Lag1,:Lag2])

    train  = selectcols(Xfull, :Year) .< Dates.Date(2005)
    test   = .!train
    Xtrain = selectrows(X, train)
    ytrain = selectrows(y, train)
    Xtest  = selectrows(X, test)
    ytest  = selectrows(y, test)

    LDA_model = LDA()
    fitresult, = fit(LDA_model, 1, Xtrain, ytrain)
    class_means, projection_matrix = fitted_params(LDA_model, fitresult)

    preds = predict(LDA_model, fitresult, Xtest)

    mce = cross_entropy(preds, ytest) |> mean

    @test 0.685 ≤ mce ≤ 0.695

    @test round.(class_means', sigdigits = 3) == [0.0428 0.0339; -0.0395 -0.0313]

    d = info_dict(LDA)
    @test d[:input_scitype] == Table(Continuous)
    @test d[:target_scitype] == AbstractVector{<:Finite}
    @test d[:name] == "LDA"
end

@testset "MLDA-2" begin
    Random.seed!(1125)
    X1 = -2 .+ randn(100, 2)
    X2 = randn(100, 2)
    X3 = 2 .+ randn(100, 2)
    y1 = ones(100)
    y2 = 2ones(100)
    y3 = 3ones(100)
    X = vcat(X1, X2, X3)
    y = vcat(y1, y2, y3)
    p = Random.randperm(300)
    X = X[p, :]
    y = y[p]
    X = table(X)
    y = categorical(y)
    train, test = partition(eachindex(y), 0.7)
    Xtrain = selectrows(X, train)
    ytrain = selectrows(y, train)
    Xtest = selectrows(X, test)
    ytest = selectrows(y, test)

    lda_model = LDA()
    fitresult, = fit(lda_model, 1, Xtrain, ytrain)
    preds = predict_mode(lda_model, fitresult, Xtest)
    mcr = misclassification_rate(preds, ytest)
    @test mcr ≤ 0.15
end

@testset "BayesianMulticlassLDA" begin
    Xfull, y = @load_smarket
    X = selectcols(Xfull, [:Lag1,:Lag2])

    train  = selectcols(Xfull, :Year) .< Dates.Date(2005)
    test   = .!train
    Xtrain = selectrows(X, train)
    ytrain = selectrows(y, train)
    Xtest  = selectrows(X, test)
    ytest  = selectrows(y, test)

    BLDA_model = BayesianLDA()
    fitresult, = fit(BLDA_model, 1, Xtrain, ytrain)
    class_means, projection_matrix, priors = fitted_params(BLDA_model, fitresult)

    preds = predict(BLDA_model, fitresult, Xtest)

    mce = cross_entropy(preds, ytest) |> mean

    @test 0.685 ≤ mce ≤ 0.695

    @test round.(class_means', sigdigits = 3) ==
        [0.0428 0.0339; -0.0395 -0.0313]

    d = info_dict(BayesianLDA)
    @test d[:input_scitype] == Table(Continuous)
    @test d[:target_scitype] == AbstractVector{<:Finite}
    @test d[:name] == "BayesianLDA"
end

@testset "BayesianSubspaceLDA" begin
    X, y = @load_iris

    LDA_model = BayesianSubspaceLDA()

    fitresult, _, report = fit(LDA_model, 1, X, y)
    class_means,projection_matrix,prior_probabilities =
        fitted_params(LDA_model, fitresult)
    preds=predict(LDA_model, fitresult, X)
    predicted_class = predict_mode(LDA_model, fitresult, X)
    mcr = misclassification_rate(predicted_class, y)
    mce = cross_entropy(preds, y) |> mean

    @test mean(abs.(class_means' - [5.006 3.428 1.462 0.246;
                          5.936 2.770 4.260 1.326;
                          6.588 2.974 5.552 2.026])) < 0.01

    @test mean(abs.(projection_matrix ≈ [0.0675721   0.00271023;
  					   0.127666    0.177718;
 					  -0.180211   -0.0767255;
 				      -0.235382    0.231435])) < 0.05

    @test round.(prior_probabilities, sigdigits=7) ==
        [0.3333333, 0.3333333, 0.3333333]
    @test round.(mcr, sigdigits=1) == 0.02

    @test round.(report.explained_variance_ratio, digits=4) == [0.9915, 0.0085]
    @test 0.04 ≤ mce ≤ 0.045

    d = info_dict(BayesianSubspaceLDA)
    @test d[:input_scitype] == Table(Continuous)
    @test d[:target_scitype] == AbstractVector{<:Finite}
    @test d[:name] == "BayesianSubspaceLDA"
end

@testset "SubspaceLDA" begin
    Random.seed!(1125)
    X1 = -2 .+ randn(100, 2)
    X2 = randn(100, 2)
    X3 = 2 .+ randn(100, 2)
    y1 = ones(100)
    y2 = 2ones(100)
    y3 = 3ones(100)
    X = vcat(X1, X2, X3)
    y = vcat(y1, y2, y3)
    p = Random.randperm(300)
    X = X[p, :]
    y = y[p]
    X = table(X)
    y = categorical(y)
    train, test = partition(eachindex(y), 0.7)
    Xtrain = selectrows(X, train)
    ytrain = selectrows(y, train)
    Xtest = selectrows(X, test)
    ytest = selectrows(y, test)

    lda_model = SubspaceLDA()
    fitresult, = fit(lda_model, 1, Xtrain, ytrain)
    preds = predict_mode(lda_model, fitresult, Xtest)
    mcr = misclassification_rate(preds, ytest)
    @test mcr ≤ 0.15

    # MultivariateStats Linear Discriminant Analysis transform
     proj    = fitresult[1].projw
     projLDA = fitresult[1].projLDA
     proj =  proj * projLDA
     XWt = matrix(X) * proj
     tlda_ms = table(XWt, prototype=X)

    # MLJ Linear Discriminant Analysis transform
     tlda_mlj = transform(lda_model, fitresult, X)

     @test tlda_mlj == tlda_ms

     d = info_dict(SubspaceLDA)
     @test d[:input_scitype] == Table(Continuous)
     @test d[:target_scitype] == AbstractVector{<:Finite}
     @test d[:name] == "SubspaceLDA"
end

end
true