text:
# homogeneous velocity and density model
nz = 51; nx = 75;
vel = 3000 * ones(nz, nx);  # m/s
vel[26:end,:] .= 3500;      # m/s
rho = 2000 * ones(nz, nx);  # kg/m^3

# vertical and horizontal grid size
dz = 10; dx = 10;

# time step size and maximum modelling length
dt = 0.001; tmax = 0.2;  # use second as unit

# receiver location
irx = collect(1 :2 :nx);
irz = 2 * ones(Int64, length(irx));

# source location
isx1 = floor(Int64, nx/2); isz1 = 2;
isx2 = collect(5:15:nx-5); ns=length(isx2); isz2 = 2*ones(ns);

# eight different test type
test_type = [(isz1, isx1, true , Float32, 1.0e-5 ),
             (isz1, isx1, true , Float64, 1.0e-10),
             (isz1, isx1, false, Float32, 1.0e-5 ),
             (isz1, isx1, false, Float64, 1.0e-10),
             (isz2, isx2, true , Float32, 1.0e-5 ),
             (isz2, isx2, true , Float64, 1.0e-10),
             (isz2, isx2, false, Float32, 1.0e-5 ),
             (isz2, isx2, false, Float64, 1.0e-10),
            ]

@testset "source side wavefield reconstruction forward and backward" begin

    for (isz, isx, free_surface, data_format, tol) in test_type

        # generate finite difference stencil
        params = TdParams(rho, vel, free_surface, dz, dx, dt, tmax;
                 data_format=data_format, fd_flag="taylor", order=3, npml=20, apml=900.);

        # create source term
        ns = length(isz)
        if ns == 1
           source_term = Source(isz, isx, params; ot=0.0, fdom=15.0, amp=100000);

        elseif ns > 1
           source_term = get_multi_sources(isz, isx, params; ot=0.3*rand(ns), fdom=15.0, amp=100000);
        end

        # save the boundary value and last wavefield
        path_bnd  = joinpath(homedir(), "Desktop/boundary_value.rsf");
        path_lwfd = joinpath(homedir(), "Desktop/last_wfd.rsf");

        # forward modeling of simultaneous sources
        multi_step_forward(source_term, params; path_bnd=path_bnd, path_lwfd=path_lwfd);

        # reconstruct source side  wavefield using boundary infomation
        source_side1 = sourceside_reconstruct_forward(path_bnd, source_term, params);

        # reconstruct source side wavefield using boundary infomation and last wavefield
        source_side2 = sourceside_reconstruct_backward(path_bnd, path_lwfd, source_term, params);

        # compute the difference
        @test sqrt(sum((source_side1 .- source_side2) .^ 2)) / sqrt(sum(source_side1 .^ 2)) < tol
    end
end

# ==============================================================================
#                  forward and adjoint of born approximation
# ==============================================================================
@testset "forward and adjoint born approimation for one shot" begin

    for (isz, isx, free_surface, data_format, tol) in test_type

        # generate finite difference stencil
        params = TdParams(rho, vel, free_surface, dz, dx, dt, tmax;
                 data_format=data_format, fd_flag="taylor", order=3, npml=20, apml=900.);

        # create source term
        ns = length(isz)
        if ns == 1
           source_term = Source(isz, isx, params; ot=0.0, fdom=15.0, amp=100000);

        elseif ns > 1
           source_term = get_multi_sources(isz, isx, params; ot=0.3*rand(ns), fdom=15.0, amp=100000);
        end

        # save the boundary value and last wavefield
        path_bnd  = joinpath(homedir(), "Desktop/boundary_value.rsf");
        path_lwfd = joinpath(homedir(), "Desktop/last_wfd.rsf");

        # forward modeling of simultaneous sources
        multi_step_forward(source_term, params; path_bnd=path_bnd, path_lwfd=path_lwfd);

        # forward born approximation
        m1   = randn(params.data_format, params.nz*params.nx)
        rec1 = born_approximation_forward(irz, irx, m1, path_bnd, source_term, params);

        # adjoint born approximation
        w   = ricker(15.0, params.dt)
        hl  = floor(Int64, length(w)/2)
        idx_o= [1]

        # generate data trace by trace
        rec2 = Recordings(irz, irx, params);
        for i = 1 : rec2.nr
            tmp = conv(randn(params.nt)*1000, w)
            copyto!(rec2.p, idx_o[1], tmp, hl+1, params.nt)
            idx_o[1] = idx_o[1] + params.nt
        end

        # adjoint
        m2 = born_approximation_adjoint(rec2, path_bnd, path_lwfd, source_term, params);

        # compute the difference
        @test (dot(m1,m2) - dot(rec1.p,rec2.p)) / (norm(m1) * norm(m2)) < tol
    end
end

# ==============================================================================
#         forward and adjoint of born approximation for multiple shots
# ==============================================================================
# test_type = [(true , Float32, 1.0e-6 ),
#              (true , Float64, 1.0e-12),
#              (false, Float32, 1.0e-6 ),
#              (false, Float64, 1.0e-12)
#             ]
#
# @testset "forward and adjoint born approimation for multiple shots" begin
#
#     for (free_surface, data_format, tol) in test_type
#
#         # generate finite difference stencil
#         params = TdParams(rho, vel, free_surface, dz, dx, dt, tmax;
#                  data_format=data_format, fd_flag="taylor", order=3, npml=20, apml=900.);
#
#         # create multiple source term
#         source_term = get_multi_sources(isz2, isx2, params; ot=0.0*rand(ns), fdom=15.0, amp=100000);
#
#         # pre-compute source side wavefield and their boundary
#         dir_sourceside = joinpath(homedir(), "Desktop/source_side");
#         get_wavefield_bound(dir_sourceside, source_term, params);
#
#         # forward born approximation
#         # generate input for forward operator
#         m1      = randn(params.data_format, params.nz*params.nx)
#         path_m1 = joinpath(homedir(), "Desktop/m1.rsf");
#         hdr1    = RegularSampleHeader(m1);
#         write_RSdata(path_m1, hdr1, m1);
#
#         dir_born1 = joinpath(homedir(), "Desktop/forward_born");
#         born_approximation_forward(dir_born1, path_m1, irz, irx, dir_sourceside, params);
#
#
#
#         # adjoint born approximation
#         w   = ricker(15.0, params.dt)
#         hl  = floor(Int64, length(w)/2)
#         idx_o= [1]
#
#         # generate data trace by trace
#         rec2 = Recordings(irz, irx, params);
#         for i = 1 : rec2.nr
#             tmp = conv(randn(params.nt)*1000, w)
#             copyto!(rec2.p, idx_o[1], tmp, hl+1, params.nt)
#             idx_o[1] = idx_o[1] + params.nt
#         end
#
#         # adjoint
#         m2 = born_approximation_adjoint(rec2, path_bnd, path_lwfd, source_term, params);
#
#         # compute the difference
#         @test (dot(m1,m2) - dot(rec1.p,rec2.p)) / (norm(m1) * norm(m2)) < tol
#     end
# end

# # ==============================================================================
# #            check the analytical gradient against numerical gradient
# # ==============================================================================
# using SeisPlot, SeisAcoustic, LinearAlgebra, DSP
#
# # size of model
# nz = 101; nx = 301;
#
# # 2-layers velocity model
# vel = 3000 * ones(nz, nx);  # m/s
# vel[51:end,:] .= 3500;
#
# # constant density model
# rho = 2000 * ones(nz, nx);  # kg/m^3
#
# # top boundary condition
# free_surface = false;
#
# # vertical and horizontal grid size
# dz = 10; dx = 10;
#
# # time step size and maximum modelling length
# dt = 0.001; tmax = 2.0;
#
# # precision
# data_format=Float64;
# order=3
#
# # organize these parameters into a structure
# params = TdParams(rho, vel, free_surface, dz, dx, dt, tmax;
#                   data_format=data_format, order=order);
#
# # initialize a source
# isz = 2; isx = 150;
# src = Source(isz, isx, params; amp=100000, fdom=15);
#
# # isx = collect(5:10:295); ns=length(isx);
# # isz = 2   * ones(ns);
# # ot  = 0.5 * rand(ns);
# # pol = rand(Bool, ns);
# # amp = 100000 * ones(ns);
# # for i = 1 : ns
# #     if !pol[i]
# #        amp[i] = - amp[i]
# #     end
# # end
# # srcs = get_multi_sources(isz, isx, params; amp=amp, ot=ot, fdom=5);
#
# # generate observed data
# irx = collect(1:2:params.nx);
# irz = 2 * ones(length(irx));
# dobs= Recordings(irz, irx, params);
#
# # generate recordings and boundary value
# @time multi_step_forward!(dobs, src, params);
# SeisPlotTX(dobs.p, cmap="gray", pclip=98,
#            xlabel="Traces", ylabel="Time (s)", dy=0.001)
#
# # ==============================================================================
# #                      the gradient based on adjoint method
# # initial velocity model
# vel0 = 3000 * ones(nz, nx);
# params0 = TdParams(rho, vel0, free_surface, dz, dx, dt, tmax;
#                    data_format=data_format, order=order);
#
# # generate synthetic data and save the boundary value
# path_bnd = joinpath(homedir(), "Desktop/bnd.rsb");
# path_wfd = joinpath(homedir(), "Desktop/wfd.rsb");
# dsyn = Recordings(irz, irx, params0);
#
# # the source-side wavefield is generated simultanes source
# isx = collect(5:10:295); ns=length(isx);
# isz = 2   * ones(ns);
# ot  = 0.5 * rand(ns);
# pol = rand(Bool, ns);
# amp = 100000 * ones(ns);
# for i = 1 : ns
#     if !pol[i]
#        amp[i] = - amp[i]
#     end
# end
# srcs = get_multi_sources(isz, isx, params; amp=amp, ot=ot, fdom=5);
#
# # born approximation forward
# multi_step_forward!(dsyn, srcs, params0, path_bnd=path_bnd, path_wfd=path_wfd);
# SeisPlotTX(dsyn.p, cmap="gray", pclip=98,
#            xlabel="Traces", ylabel="Time (s)", dy=0.001)
#
# # compute the residue
# dres = get_residue(dsyn, dobs);
# SeisPlotTX(dres.p, cmap="gray", pclip=98,
#            xlabel="Traces", ylabel="Time (s)", dy=0.001)
#
# # set the scatter
# # m = zeros(data_format, params.nz, params.nx)
# # m[50,:] .= 100000.;
# # m = vec(m);
#
# # born approximation forward
# m1 = 10000. * randn(params.nz * params.nx);
# rec1 = Recordings(irz, irx, params0);
# @time born_approximation_forward!(rec1, m1, path_bnd, srcs, params0);
# SeisPlotTX(rec1.p, cmap="gray", pclip=98,
#            xlabel="Traces", ylabel="Time (s)", dy=0.001)
#
# # make band-limited random recordings
# w    = ricker(15.0, params.dt)
# hl   = floor(Int64, length(w)/2)
# rec2 = Recordings(irz, irx, params);
# idx_o= [1]
# for i = 1 : rec2.nr
#     tmp = conv(randn(params0.nt)*1000, w)
#     copyto!(rec2.p, idx_o[1], tmp, hl+1, params0.nt)
#     idx_o[1] = idx_o[1] + params0.nt
# end
# SeisPlotTX(rec2.p, cmap="gray", pclip=98,
#            xlabel="Traces", ylabel="Time (s)", dy=0.001)
#
# # born approximation adjoint
# m2 = born_approximation_adjoint(rec2, path_bnd, path_wfd, srcs, params0);
#
# # dot product test
# tmp1 = dot(m1, m2)
#
# p1 = rec1.p[2:end,:];
# p2 = rec2.p[2:end,:];
# tmp2 = dot(vec(p1), vec(p2))
#
# (tmp1-tmp2) / tmp1
#
#
# # input and output for born approximation
# path_m1 = joinpath(dir_work, "cgls_model_space/image1.rsf");
# path_m2 = joinpath(dir_work, "cgls_model_space/image2.rsf");
# dir_d1  = joinpath(dir_work, "cgls_data_space/born_forward");
# dir_d2  = dir_obs;
#
# image  = 1000.0 * randn(fidiff.data_format, fidiff.nz*fidiff.nx);
# hdr    = RegularSampleHeader(image, title="random image");
# write_RSdata(path_m1, hdr, image);
#
# # apply operation
# born_params = (irz=irz, irx=irx, dir_sourceside=dir_sourceside, fidiff=fidiff, normalization_flag=true, mute_index=10);
# born_approximation(dir_d1, path_m1, 1; born_params...);
# born_approximation(path_m2, dir_d2, 2; born_params...);
#
# # test
# (hdr, m1) = read_RSdata(path_m1);
# (hdr, m2) = read_RSdata(path_m2);
# tmp1 = dot(m1, m2);
#
# tmp2 = [0.0]
# for i = 1 : length(src)
#     file_name = file_name = join(["recordings_" "$i" ".bin"])
#     rec1 = read_recordings(joinpath(dir_d1, file_name))
#     rec2 = read_recordings(joinpath(dir_d2, file_name))
#     tmp2[1] += dot(rec1.p, rec2.p)
# end
# (tmp2[1]-tmp1) / tmp2[1]
#
#
# # plot
# # rec = read_recordings(joinpath(dir_d1, "recordings_1.bin"));
# # SeisPlotTX(rec.p, dy=dt, cmap="gray", wbox=6, hbox=6);
# # (hdr, m2) = read_RSdata(path_m2);
# # SeisPlotTX(reshape(m2, fidiff.nz, fidiff.nx), cmap="gray", wbox=10, hbox=5);
# # m3 = laplace_filter(reshape(m2, fidiff.nz, fidiff.nx));
# # SeisPlotTX(m3, cmap="gray", wbox=10, hbox=5, pclip=96);
# #
# # tmp = zeros(fidiff.data_format, fidiff.nz, fidiff.nx)
# # for i2 = 1 : fidiff.nx
# #     for i1 = 2 : fidiff.nz
# #         tmp[i1,i2] = (rho[i1,i2]-rho[i1-1,i2]) / (rho[i1,i2]+rho[i1-1,i2])
# #     end
# # end
# # SeisPlotTX(tmp, cmap="gray", wbox=10, hbox=5);
#
# # examine the result
# # i = 50
# # p1 = dres.p[:,i] / maximum(abs.(dres.p[:,i]))
# # p2 = d_born.p[:,i] / maximum(abs.(d_born.p[:,i]))
# # figure(); plot(p1); plot(p2);


text:
using Pkg

if !haskey(Pkg.installed(), "GaussDCA")
    Pkg.add(PackageSpec(url="https://github.com/carlobaldassi/GaussDCA.jl"))
end

msa = map(Residue, rand(1:21,100,20))
dca = gaussdca(msa, min_separation=2)

@test  isnan(dca[1,1])
@test  isnan(dca[1,2])
@test !isnan(dca[1,3])

ext: jl


text:
module Snowflakes

#core opt
using JuMP
using MathOptInterface

#data extraction
using XLSX
using DataFrames
using Query
using Distances

#support
using Random
using Combinatorics
using UUIDs
using Statistics

include("Structures.jl")
include("Base.jl")
include("Column.jl")
include("Settings.jl")

export b
export extract!
export stats
export initStab
export root
export origin

export set_default_optimizer!
export get_default_optimizer
export reset_default_optimizer

export master
export buildMaster
export getDuals
export sub
export buildSub
export getCols
export updateStab!
export checkStab
export colGen

end


text:
#Two actions, zero observation, one timestep long, action-dependent +1/-1 reward: The first
#env to exercise the policy! If my agent can't learn to pick the better action, there's
#something wrong with either my advantage calculations, my policy loss or my policy update.
#That's three things, but it's easy to work out by hand the expected values for each one and
# check that the values produced by your actual code line up with them.
import AlphaZero.GI
using Flux: onehot
struct GameSpec <: GI.AbstractGameSpec end

mutable struct GameEnv <: GI.AbstractGameEnv
  terminated :: Bool
  numero :: Int
end

const INITSTATE = (; terminated=false, numero=0)

GI.init(::GameSpec, state=INITSTATE) = GameEnv(state.terminated, state.numero)

GI.spec(::GameEnv) = GameSpec()

GI.two_players(::GameSpec) = false

const ACTIONS = collect(1:2)

GI.actions(::GameSpec) = ACTIONS

GI.actions_mask(g::GameEnv) = trues(2)

GI.current_state(g::GameEnv) = (; terminated=g.terminated, numero=g.numero)

GI.white_playing(g::GameEnv) = true

GI.game_terminated(g::GameEnv) = g.terminated

function GI.white_reward(g::GameEnv)
  if g.terminated
    return g.numero == 1 ? 1. : -1.
  else
    return 0.
  end
end

function GI.play!(g::GameEnv, action)
  g.numero += action
  g.terminated = true
end

GI.vectorize_state(::GameSpec, state) = collect(onehot(state.numero, 0:2))


text:
n = parse(Int64, readline())
A = map((x) -> parse(Int64,x) , split(readline(), " "))

for i in n:-1:1
    print(A[i]," ")
end
println()

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 07.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
"""
    open(vis::Visualizer; host::IPAddr = ip"127.0.0.1", start_browser = true,
                          default_port = 8700, max_retries = 500)

Start a server for the visualizer so that it can be accessed from a browser
using the given host URL. This method will try to search for an open port,
starting at `default_port` and then trying `max_retries` additional port
numbers. If `start_browser` is true, then a new browser window will be opened.

    open(vis::Visualizer, port::Integer; host::IPAddr = ip"127.0.0.1", start_browser = true)

Start a server for the visualizer so that it can be accessed from a browser
using the given host URL and port. If `start_browser` is true, then a new
browser window will be opened.

"""
Base.open(vis::Visualizer, args...; kw...) = open(vis.core, args...; kw...)

function Base.open(core::CoreVisualizer, port::Integer;
                   host::IPAddr = ip"127.0.0.1", start_browser = true)
    @async WebIO.webio_serve(Mux.page("/", req -> core.scope), host, port)
    url = "http://$host:$port"
    @info("Serving MeshCat visualizer at $url")
    start_browser && open_url(url)
end

function open_url(url)
    try
        if Sys.iswindows()
            run(`cmd.exe /C "start $url"`)
        elseif Sys.isapple()
            run(`open $url`)
        elseif Sys.islinux()
            run(`xdg-open $url`)
        end
    catch e
        println("Could not open browser automatically: $e")
        println("Please open the following URL in your browser:")
        println(url)
    end
end

function Base.open(core::CoreVisualizer;
                   host::IPAddr=ip"127.0.0.1", default_port=8700, max_retries=500, kw...)
    for port in default_port:(default_port + max_retries)
        server = try
            listen(host, port)
        catch e
            if e isa Base.IOError
                continue
            end
        end
        close(server)
        # It is *possible* that a race condition could occur here, in which
        # some other process grabs the given port in between the close() above
        # and the open() below. But it's unlikely and would not be terribly
        # damaging (the user would just have to call open() again).
        return open(core, port; host=host, kw...)
    end
end


text:
# Fortran 90

function f90setup()
    quote
        const f90_filename = Ref{AbstractString}("code.f90")
        function set_f90_filename!(filename::AbstractString)
            return f90_filename[] = filename
        end

        const f90_chunks = Dict{AbstractString,Vector{String}}()
        function f90_add_code!(f90_code::AbstractString)
            return push!(get!(f90_chunks, f90_filename[], String[]), f90_code)
        end
        function f90_get_code()
            iobuffer = IOBuffer()
            allcode = Dict{AbstractString,String}()
            for (filename, chunks) in f90_chunks
                for chunk in chunks
                    println(iobuffer, chunk)
                end
                allcode[filename] = String(take!(iobuffer))
            end
            return allcode
        end
        function f90_write_code!()
            println("Generating Fortran code:")
            allcode = f90_get_code()
            for (filename, content) in allcode
                println("Generating $filename...")
                write(filename, content)
            end
            println("Done.")
            return nothing
        end
    end
end
export f90setup

function f90newfile(filename::AbstractString, f90_stmts::AbstractString)
    f90_code = """
               $f90_stmts
               """
    quote
        set_f90_filename!($filename)
        f90_add_code!($f90_code)
    end
end
export f90newfile

function f90code(f90_stmts::AbstractString)
    f90_code = """
               $f90_stmts
               """
    quote
        f90_add_code!($f90_code)
    end
end
export f90code

function f90function(name::FnName, result::FnResult, arguments::AbstractVector{FnArg}, f90_stmts::AbstractString)
    julia_code = quote
        function $(name.julia_name)($([:($(arg.julia_name)::$(arg.initial_julia_type)) for arg in arguments]...))
            res = ccall(($(QuoteNode(Symbol(name.cxx_name, "_"))), $(name.cxx_library)), $(result.julia_type),
                        ($([arg.julia_type for arg in arguments if !arg.skip]...),),
                        $([arg.convert_from_initial(arg.julia_name) for arg in arguments if !arg.skip]...))
            return $(result.convert_to_final(:res))::$(result.final_julia_type)
        end
    end
    julia_code = clean_code(julia_code)
    simple_julia_code = simplify_code(julia_code)

    f90_code = """
        ! $(join(split(string(simple_julia_code), "\n"), "\n! "))
        $(result.cxx_type == "void" ? "subroutine" : "$(result.cxx_type) function") $(name.cxx_name)($(join(["$(arg.cxx_name)" for arg in arguments if !arg.skip], ", ")))
            implicit none
            $(join(["$(arg.cxx_type) :: $(arg.cxx_name)" for arg in arguments if !arg.skip], "\n    "))
            $f90_stmts
        end $(result.cxx_type == "void" ? "subroutine" : "function") $(name.cxx_name)
        """
    quote
        f90_add_code!($f90_code)
        $julia_code
    end
end
export f90function


text:
doc"""
  `Unums.Gflags` is an internal type strtucture thats carries calculation
  override information about the gnum.  Certain special values (infinity, NaN,
  zero, etc.) can trigger overrides of the default calculation and proceed to
  a faster calculation.  This needs to be in a separate, to allow the gnum
  itself to be immutable, and code-optimized.
"""
type Gflags
  lower::UInt16
  upper::UInt16
end

doc"""
  `Unums.Gnum{ESS,FSS}` is an internal type structure that's not exposed to
  the global scope.  It represents the g-layer from the "End of Error".  Some
  of the status flags are stored in the flags parameter of the 'lower' term.

  this will be mostly used by the internal `@glayer` directive which will take
  an expression and convert it to a series of functions, which, at compile-time,
  will allocate a set of 'global' gnum registers, then pass them through a chain
  of mutator functions which will keep the gnum value in the registers, then
  output the whole thing as a single Utype value as needed (Unum or Ubound)

  The Gnum value also keeps a 'buffer' variable to store calculated unums for
  certain operations.  This variable is DISTINCT from the 'scratchpad' in that
  it does NOT carry extra precision, and should only be used when operations
  need to use an intermediate value before or after a calculation.  For example,
  mmr - (val) requires calculating a lower bound starting with bigexact.  Or
  multiplication of ubounds requires storing unum variables for comparison.
"""
abstract Gnum{ESS,FSS} <: Real

immutable GnumSmall{ESS,FSS} <: Gnum{ESS,FSS}
  lower::UnumSmall{ESS,FSS}
  upper::UnumSmall{ESS,FSS}
  #the buffer a preallocated location to store unum values.
  buffer::UnumSmall{ESS,FSS}
  #flags is a section that holds flagged values about the lower and upper unums
  flags::Gflags
end

immutable GnumLarge{ESS,FSS} <: Gnum{ESS,FSS}
  lower::UnumLarge{ESS,FSS}
  upper::UnumLarge{ESS,FSS}
  #the buffer a preallocated location to store unum values.
  buffer::UnumLarge{ESS,FSS}
  #flags is a section that holds flagged values about the lower and upper unums]
  flags::Gflags
end

#these g-layer values go into the scratchpad to indicate properties of the gnum.
GNUM_ONESIDED_MASK = 0x8000
#throws a bit saying this number is NaN.
GNUM_NAN_MASK  = 0x4000
GNUM_SFLAGS_MASK = 0xC000

#g-layer flags that apply to both "lower" and "higher" slots.
#temporary flag saying to ignore calculations on this side of the gnum.
GNUM_IGNORE_SIDE_MASK = 0x2000

#################################
#precedence of bits:
# NAN > SBIT > IGNORE
# A nan is automatically single, even if sbit isn't thrown.
# An SBIT automatically ignores the second value, even if ignore isn't thrown.

#informational flags that report on the identity of the number.
GNUM_FLAG_MASK = 0x0F00
GNUM_INF_MASK  = 0x0800
GNUM_MMR_MASK  = 0x0400
GNUM_SSS_MASK  = 0x0200
GNUM_ZERO_MASK = 0x0100

#make lower and upper side designations more elegant.
const LOWER_UNUM = Val{:lower}
const UPPER_UNUM = Val{:upper}
const SCRATCHPAD = Val{:scratchpad}
const BUFFER = Val{:buffer}
#other values used for instruction directives in mults and divs
const DISCARD_PRIMARY = Val{:discard_primary}
const DISCARD_SECONDARY = Val{:discard_secondary}

#create a global scratchpad array.
const GLOBAL_SCRATCHPAD_SIZE = __cell_length(11) + (__cell_length(11) >> 1)
const GLOBAL_SCRATCHPAD = zeros(UInt64, GLOBAL_SCRATCHPAD_SIZE)

#scratches an operation.
macro scratch_this_operation!(s)
  esc(quote
    nan!($s)
    ignore_both_sides!($s)
    return
  end)
end

macro init_sflags()
  esc(:(sflags::UInt16))
end

macro preserve_sflags(s, expr)
  esc(quote
    sflags = $s.scratchpad.flags & GNUM_SFLAGS_MASK
    $expr
    $s.scratchpad.flags |= sflags
  end)
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
module Omniscape
using Circuitscape
using LinearAlgebra.BLAS
using Base.Threads
using StatsBase
using Statistics

include("config.jl")
include("consts.jl")
include("functions.jl")
include("run_omniscape.jl")
include("errors_warnings.jl")
export run_omniscape

end


text:
# In the various algorithms we store the tolerance as a generic Real in order to
# construct these objects using keyword arguments in a type-stable manner. At
# the beginning of the corresponding algorithm, the actual tolerance will be
# converted to the concrete numeric type appropriate for the problem at hand

# Orthogonalization and orthonormalization
"""
    abstract type Orthogonalizer

Supertype of a hierarchy for representing different orthogonalization strategies or
algorithms.

See also: [`ClassicalGramSchmidt`](@ref), [`ModifiedGramSchmidt`](@ref),
[`ClassicalGramSchmidt2`](@ref), [`ModifiedGramSchmidt2`](@ref),
[`ClassicalGramSchmidtIR`](@ref), [`ModifiedGramSchmidtIR`](@ref).
"""
abstract type Orthogonalizer end
abstract type Reorthogonalizer <: Orthogonalizer end

# Simple
"""
    ClassicalGramSchmidt()

Represents the classical Gram Schmidt algorithm for orthogonalizing different vectors,
typically not an optimal choice.
"""
struct ClassicalGramSchmidt <: Orthogonalizer
end

"""
    ModifiedGramSchmidt()

Represents the modified Gram Schmidt algorithm for orthogonalizing different vectors,
typically a reasonable choice for linear systems but not for eigenvalue solvers with a
large Krylov dimension.
"""
struct ModifiedGramSchmidt <: Orthogonalizer
end

# A single reorthogonalization always
"""
    ClassicalGramSchmidt2()

Represents the classical Gram Schmidt algorithm with a second reorthogonalization step
always taking place.
"""
struct ClassicalGramSchmidt2 <: Reorthogonalizer
end

"""
    ModifiedGramSchmidt2()

Represents the modified Gram Schmidt algorithm with a second reorthogonalization step
always taking place.
"""
struct ModifiedGramSchmidt2 <: Reorthogonalizer
end

# Iterative reorthogonalization
"""
    ClassicalGramSchmidtIR(η::Real = 1/sqrt(2))

Represents the classical Gram Schmidt algorithm with iterative (i.e. zero or more)
reorthogonalization until the norm of the vector after an orthogonalization step has not
decreased by a factor smaller than `η` with respect to the norm before the step. The
default value corresponds to the Daniel-Gragg-Kaufman-Stewart condition.
"""
struct ClassicalGramSchmidtIR{S<:Real} <: Reorthogonalizer
    η::S
end
ClassicalGramSchmidtIR() = ClassicalGramSchmidtIR(1/sqrt(2)) # Daniel-Gragg-Kaufman-Stewart

"""
    ModifiedGramSchmidtIR(η::Real = 1/sqrt(2))

Represents the modified Gram Schmidt algorithm with iterative (i.e. zero or more)
reorthogonalization until the norm of the vector after an orthogonalization step has not
decreased by a factor smaller than `η` with respect to the norm before the step. The
default value corresponds to the Daniel-Gragg-Kaufman-Stewart condition.
"""
struct ModifiedGramSchmidtIR{S<:Real} <: Reorthogonalizer
    η::S
end
ModifiedGramSchmidtIR() = ModifiedGramSchmidtIR(1/sqrt(2)) # Daniel-Gragg-Kaufman-Stewart

# Solving eigenvalue problems
abstract type KrylovAlgorithm end

# General purpose; good for linear systems, eigensystems and matrix functions
"""
    Lanczos(; krylovdim = KrylovDefaults.krylovdim, maxiter = KrylovDefaults.maxiter,
        tol = KrylovDefaults.tol, orth = KrylovDefaults.orth, eager = false, verbosity = 0)

Represents the Lanczos algorithm for building the Krylov subspace; assumes the linear
operator is real symmetric or complex Hermitian. Can be used in `eigsolve` and
`exponentiate`. The corresponding algorithms will build a Krylov subspace of size at most
`krylovdim`, which will be repeated at most `maxiter` times and will stop when the norm of
the residual of the Lanczos factorization is smaller than `tol`. The orthogonalizer `orth`
will be used to orthogonalize the different Krylov vectors. Eager mode, as selected by
`eager = true`, means that the algorithm that uses this Lanczos process (e.g. `eigsolve`)
can try to finish its computation before the total Krylov subspace of dimension `krylovdim`
is constructed. Default verbosity level `verbosity` is zero, meaning that no output will be
printed.

Use `Arnoldi` for non-symmetric or non-Hermitian linear operators.

See also: `factorize`, `eigsolve`, `exponentiate`, `Arnoldi`, `Orthogonalizer`
"""
struct Lanczos{O<:Orthogonalizer, S<:Real} <: KrylovAlgorithm
    orth::O
    krylovdim::Int
    maxiter::Int
    tol::S
    eager::Bool
    verbosity::Int
end
Lanczos(; krylovdim::Int = KrylovDefaults.krylovdim,
            maxiter::Int = KrylovDefaults.maxiter,
            tol::Real = KrylovDefaults.tol,
            orth::Orthogonalizer = KrylovDefaults.orth,
            eager::Bool = false,
            verbosity::Int = 0) =
    Lanczos(orth, krylovdim, maxiter, tol, eager, verbosity)

"""
    GKL(; krylovdim = KrylovDefaults.krylovdim, maxiter = KrylovDefaults.maxiter,
        tol = KrylovDefaults.tol, orth = KrylovDefaults.orth, verbosity = 0)

Represents the Golub-Kahan-Lanczos bidiagonalization algorithm for sequentially building a
Krylov-like factorization of a general matrix or linear operator with a bidiagonal reduced
matrix. Can be used in `svdsolve`. The corresponding algorithm builds a Krylov subspace of
size at most `krylovdim`, which will be repeated at most `maxiter` times and will stop when
the norm of the residual of the Arnoldi factorization is smaller than `tol`. The
orthogonalizer `orth` will be used to orthogonalize the different Krylov vectors. Default
verbosity level `verbosity` is zero, meaning that no output will be printed.

See also: `svdsolve`, `Orthogonalizer`
"""
struct GKL{O<:Orthogonalizer, S<:Real} <: KrylovAlgorithm
    orth::O
    krylovdim::Int
    maxiter::Int
    tol::S
    eager::Bool
    verbosity::Int
end
GKL(; krylovdim::Int = KrylovDefaults.krylovdim,
        maxiter::Int = KrylovDefaults.maxiter,
        tol::Real = KrylovDefaults.tol,
        orth::Orthogonalizer = KrylovDefaults.orth,
        eager::Bool = false,
        verbosity::Int = 0) =
    GKL(orth, krylovdim, maxiter, tol, eager, verbosity)

"""
    Arnoldi(; krylovdim = KrylovDefaults.krylovdim, maxiter = KrylovDefaults.maxiter,
        tol = KrylovDefaults.tol, orth = KrylovDefaults.orth, verbosity = 0)

Represents the Arnoldi algorithm for building the Krylov subspace for a general matrix or
linear operator. Can be used in `eigsolve` and `exponentiate`. The corresponding algorithms
will build a Krylov subspace of size at most `krylovdim`, which will be repeated at most
`maxiter` times and will stop when the norm of the residual of the Arnoldi factorization is
smaller than `tol`. The orthogonalizer `orth` will be used to orthogonalize the different
Krylov vectors. Eager mode, as selected by `eager = true`, means that the algorithm that
uses this Arnoldi process (e.g. `eigsolve`) can try to finish its computation before the
total Krylov subspace of dimension `krylovdim` is constructed. Default verbosity level
`verbosity` is zero, meaning that no output will be printed.

Use `Lanczos` for real symmetric or complex Hermitian linear operators.

See also: [`eigsolve`](@ref), [`exponentiate`](@ref), [`Lanczos`](@ref),
[`Orthogonalizer`](@ref)
"""
struct Arnoldi{O<:Orthogonalizer, S<:Real} <: KrylovAlgorithm
    orth::O
    krylovdim::Int
    maxiter::Int
    tol::S
    eager::Bool
    verbosity::Int
end
Arnoldi(; krylovdim::Int = KrylovDefaults.krylovdim,
            maxiter::Int = KrylovDefaults.maxiter,
            tol::Real = KrylovDefaults.tol,
            orth::Orthogonalizer = KrylovDefaults.orth,
            eager::Bool = false,
            verbosity::Int = 0) =
    Arnoldi(orth, krylovdim, maxiter, tol, eager, verbosity)

"""
    GolubYe(; krylovdim = KrylovDefaults.krylovdim, maxiter = KrylovDefaults.maxiter,
        tol = KrylovDefaults.tol, orth = KrylovDefaults.orth, verbosity = 0)

Represents the Golub-Ye algorithm for solving hermitian (symmetric) generalized eigenvalue
problems `A x = λ B x` with positive definite `B`, without the need for inverting `B`.
Builds a Krylov subspace of size `krylovdim` starting from an estimate `x` by acting with
`(A - ρ(x) B)`, where `ρ(x) = dot(x, A*x)/dot(x, B*x)`, and employing the Lanczos
algorithm. This process is repeated at most `maxiter` times. In every iteration `k>1`, the
subspace will also be expanded to size `krylovdim+1` by adding ``x_k - x_{k-1}``, which is
known as the LOPCG correction and was suggested by Money and Ye. With `krylovdim = 2`, this
algorithm becomes equivalent to `LOPCG`.
"""
struct GolubYe{O<:Orthogonalizer, S<:Real} <: KrylovAlgorithm
    orth::O
    krylovdim::Int
    maxiter::Int
    tol::S
    verbosity::Int
end
GolubYe(; krylovdim::Int = KrylovDefaults.krylovdim,
            maxiter::Int = KrylovDefaults.maxiter,
            tol::Real = KrylovDefaults.tol,
            orth::Orthogonalizer = KrylovDefaults.orth,
            verbosity::Int = 0) =
    GolubYe(orth, krylovdim, maxiter, tol, verbosity)

# Solving linear systems specifically
abstract type LinearSolver <: KrylovAlgorithm end

"""
    CG(; maxiter = KrylovDefaults.maxiter, tol = KrylovDefaults.tol)

Construct an instance of the conjugate gradient algorithm with specified parameters, which
can be passed to `linsolve` in order to iteratively solve a linear system with a positive
definite (and thus symmetric or hermitian) coefficient matrix or operator. The `CG` method
will search for the optimal `x` in a Krylov subspace of maximal size `maxiter`, or stop when
`norm(A*x - b) < tol`. Default verbosity level `verbosity` is zero, meaning that no output
will be printed.

See also: [`linsolve`](@ref), [`MINRES`](@ref), [`GMRES`](@ref), [`BiCG`](@ref),
[`BiCGStab`](@ref)
"""
struct CG{S<:Real} <: LinearSolver
    maxiter::Int
    tol::S
    verbosity::Int
end
CG(; maxiter::Integer = KrylovDefaults.maxiter,
        tol::Real = KrylovDefaults.tol,
        verbosity::Int = 0) =
    CG(maxiter, tol, verbosity)

"""
    GMRES(; krylovdim = KrylovDefaults.krylovdim, maxiter = KrylovDefaults.maxiter,
        tol = KrylovDefaults.tol, orth::Orthogonalizer = KrylovDefaults.orth)

Construct an instance of the GMRES algorithm with specified parameters, which can be passed
to `linsolve` in order to iteratively solve a linear system. The `GMRES` method will search
for the optimal `x` in a Krylov subspace of maximal size `krylovdim`, and repeat this
process for at most `maxiter` times, or stop when `norm(A*x - b) < tol`. In building the
Krylov subspace, `GMRES` will use the orthogonalizer `orth`. Default verbosity level
`verbosity` is zero, meaning that no output will be printed.

Note that in the traditional nomenclature of `GMRES`, the parameter `krylovdim` is referred
to as the restart parameter, and `maxiter` is the number of outer iterations, i.e. restart
cycles. The total iteration count, i.e. the number of expansion steps, is roughly
`krylovdim` times the number of iterations.

See also: [`linsolve`](@ref), [`BiCG`](@ref), [`BiCGStab`](@ref), [`CG`](@ref),
[`MINRES`](@ref)
"""
struct GMRES{O<:Orthogonalizer,S<:Real} <: LinearSolver
    orth::O
    maxiter::Int
    krylovdim::Int
    tol::S
    verbosity::Int
end
GMRES(; krylovdim::Integer = KrylovDefaults.krylovdim,
        maxiter::Integer = KrylovDefaults.maxiter,
        tol::Real = KrylovDefaults.tol,
        orth::Orthogonalizer = KrylovDefaults.orth,
        verbosity::Int = 0) =
    GMRES(orth, maxiter, krylovdim, tol, verbosity)

# TODO
"""
    MINRES(; maxiter = KrylovDefaults.maxiter, tol = KrylovDefaults.tol)

Construct an instance of the conjugate gradient algorithm with specified parameters, which
can be passed to `linsolve` in order to iteratively solve a linear system with a real
symmetric or complex hermitian coefficient matrix or operator. The `MINRES` method will
search for the optimal `x` in a Krylov subspace of maximal size `maxiter`, or stop when
`norm(A*x - b) < tol`. In building the Krylov subspace, `MINRES` will use the
orthogonalizer `orth`. Default verbosity level `verbosity` is zero, meaning that no output
will be printed.

!!! warning "Not implemented yet"

See also: [`linsolve`](@ref), [`CG`](@ref), [`GMRES`](@ref), [`BiCG`](@ref), [`BiCGStab`](@ref)
"""
struct MINRES{S<:Real} <: LinearSolver
    maxiter::Int
    tol::S
    verbosity::Int
end
MINRES(; maxiter::Integer = KrylovDefaults.maxiter,
            tol::Real = KrylovDefaults.tol,
            verbosity::Int = 0) =
    MINRES(maxiter, tol, verbosity)

"""
    BiCG(; maxiter = KrylovDefaults.maxiter, tol = KrylovDefaults.tol)

Construct an instance of the Biconjugate gradient algorithm with specified parameters, which
can be passed to `linsolve` in order to iteratively solve a linear system general linear
map, of which the adjoint can also be applied. The `BiCG` method will search for the
optimal `x` in a Krylov subspace of maximal size `maxiter`, or stop when `norm(A*x - b) <
tol`. Default verbosity level `verbosity` is zero, meaning that no output will be printed.

!!! warning "Not implemented yet"

See also: [`linsolve`](@ref), [`BiCGStab`](@ref), [`GMRES`](@ref), [`CG`](@ref),
[`MINRES`](@ref)
"""
struct BiCG{S<:Real} <: LinearSolver
    maxiter::Int
    tol::S
    verbosity::Int
end
BiCG(; maxiter::Integer = KrylovDefaults.maxiter,
        tol::Real = KrylovDefaults.tol,
        verbosity::Int = 0) =
    BiCG(maxiter, tol, verbosity)

"""
    BiCGStab(; maxiter = KrylovDefaults.maxiter, tol = KrylovDefaults.tol)

Construct an instance of the Biconjugate gradient algorithm with specified parameters, which
can be passed to `linsolve` in order to iteratively solve a linear system general linear
map. The `BiCGStab` method will search for the optimal `x` in a Krylov subspace of maximal
size `maxiter`, or stop when `norm(A*x - b) < tol`. Default verbosity level `verbosity` is
zero, meaning that no output will be printed.

!!! warning "Not implemented yet"

See also: [`linsolve`](@ref), [`BiCG`](@ref), [`GMRES`](@ref), [`CG`](@ref),
[`MINRES`](@ref)
"""
struct BiCGStab{S<:Real} <: LinearSolver
    maxiter::Int
    tol::S
    verbosity::Int
end
BiCGStab(; maxiter::Integer = KrylovDefaults.maxiter,
            tol::Real = KrylovDefaults.tol,
            verbosity::Int = 0) =
    BiCGStab(maxiter, tol, verbosity)

# Solving eigenvalue systems specifically
abstract type EigenSolver <: KrylovAlgorithm end

struct JacobiDavidson <: EigenSolver
end

# Default values
"""
    module KrylovDefaults
        const orth = KrylovKit.ModifiedGramSchmidtIR()
        const krylovdim = 30
        const maxiter = 100
        const tol = 1e-12
    end

A module listing the default values for the typical parameters in Krylov based algorithms:
*   `orth`: the orthogonalization routine used to orthogonalize the Krylov basis in the
    `Lanczos` or `Arnoldi` iteration
*   `krylovdim`: the maximal dimension of the Krylov subspace that will be constructed
*   `maxiter`: the maximal number of outer iterations, i.e. the maximum number of times the
    Krylov subspace may be rebuilt
*   `tol`: the tolerance to which the problem must be solved, based on a suitable error
    measure, e.g. the norm of some residual.
!!! warning
    The default value of `tol` is a `Float64` value, if you solve problems in `Float32` or
    `ComplexF32` arithmetic, you should always specify a new `tol` as the default value
    will not be attainable.
"""
module KrylovDefaults
    using ..KrylovKit
    const orth = KrylovKit.ModifiedGramSchmidt2() # conservative choice
    const krylovdim = 30
    const maxiter = 100
    const tol = 1e-12
end


text:
module LibCEED

using StaticArrays, UnsafeArrays, Requires, Preferences
using Base: RefValue

# import low-level C interface
include("C.jl")
import .C

export @interior_qf,
    @witharray,
    @witharray_read,
    Abscissa,
    AbscissaAndWeights,
    Basis,
    BasisCollocated,
    COPY_VALUES,
    Ceed,
    CeedDim,
    CeedInt,
    CeedScalar,
    CeedVector,
    CeedVectorActive,
    CeedVectorNone,
    Context,
    CopyMode,
    EVAL_CURL,
    EVAL_DIV,
    EVAL_GRAD,
    EVAL_INTERP,
    EVAL_NONE,
    EVAL_WEIGHT,
    ElemRestriction,
    ElemRestrictionNone,
    EvalMode,
    GAUSS,
    GAUSS_LOBATTO,
    HEX,
    LINE,
    MEM_DEVICE,
    MEM_HOST,
    MemType,
    NORM_1,
    NORM_2,
    NORM_MAX,
    NOTRANSPOSE,
    NormType,
    OWN_POINTER,
    Operator,
    PRISM,
    PYRAMID,
    QFunction,
    QFunctionNone,
    QUAD,
    QuadMode,
    RequestImmediate,
    RequestOrdered,
    STRIDES_BACKEND,
    TET,
    TRANSPOSE,
    TRIANGLE,
    Topology,
    TransposeMode,
    USE_POINTER,
    UserQFunction,
    add_input!,
    add_output!,
    apply!,
    apply_add!,
    apply,
    assemble,
    assemble_add_diagonal!,
    assemble_diagonal!,
    axpy!,
    ceedversion,
    create_composite_operator,
    create_elem_restriction,
    create_elem_restriction_strided,
    create_evector,
    create_h1_basis,
    create_identity_qfunction,
    create_interior_qfunction,
    create_lvector,
    create_tensor_h1_basis,
    create_tensor_h1_lagrange_basis,
    create_vectors,
    det,
    extract_array,
    extract_context,
    gauss_quadrature,
    get_libceed_path,
    get_preferred_memtype,
    getcompstride,
    getnumelements,
    getelementsize,
    getlvectorsize,
    getmultiplicity!,
    getmultiplicity,
    getdimension,
    getgrad,
    getgrad1d,
    getinterp,
    getinterp1d,
    getnumcomponents,
    getnumnodes,
    getnumnodes1d,
    getnumqpts,
    getnumqpts1d,
    getqref,
    getqweights,
    getresource,
    gettopology,
    getvoigt!,
    getvoigt,
    iscuda,
    isdeterministic,
    isrelease,
    lobatto_quadrature,
    norm,
    pointwisemult!,
    reciprocal!,
    scale!,
    set_context!,
    set_cufunction!,
    set_data!,
    set_field!,
    set_libceed_path!,
    setarray!,
    setvalue!,
    setvoigt!,
    setvoigt,
    syncarray!,
    takearray!,
    use_prebuilt_libceed!,
    witharray,
    witharray_read

include("Globals.jl")
include("Ceed.jl")
include("CeedVector.jl")
include("Basis.jl")
include("ElemRestriction.jl")
include("Quadrature.jl")
include("Context.jl")
include("UserQFunction.jl")
include("QFunction.jl")
include("Request.jl")
include("Operator.jl")
include("Misc.jl")

const minimum_libceed_version = v"0.8.0"

function __init__()
    libceed_version = ceedversion()
    if minimum_libceed_version > libceed_version
        @warn("""
              Incompatible libCEED version.
              LibCEED.jl requires libCEED version at least $minimum_libceed_version.
              The version of the libCEED library is $libceed_version."
              """)
    end
    @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" include("Cuda.jl")
    set_globals()
end

"""
    ceedversion()

Returns a `VersionNumber` corresponding to the version of the libCEED library currently used.
"""
function ceedversion()
    major = Ref{Cint}()
    minor = Ref{Cint}()
    patch = Ref{Cint}()
    release = Ref{Bool}()
    C.CeedGetVersion(major, minor, patch, release)
    return VersionNumber(major[], minor[], patch[])
end

"""
    isrelease()

Returns true if the libCEED library is a release build, false otherwise.
"""
function isrelease()
    major = Ref{Cint}()
    minor = Ref{Cint}()
    patch = Ref{Cint}()
    release = Ref{Bool}()
    C.CeedGetVersion(major, minor, patch, release)
    return release[]
end

"""
    set_libceed_path!(path::AbstractString)
    set_libceed_path!(:prebuilt)
    set_libceed_path!(:default)

Sets the path of the libCEED dynamic library. `path` should be the absolute path to the library
file.

`set_libceed_path!(:prebuilt)` indicates to LibCEED.jl to use the prebuilt version of the libCEED
library (bundled with libCEED\\_jll).

`set_libceed_path!(:default)` indicates to LibCEED.jl to use the default library. This usually has
the same effect as `set_libceed_path!(:prebuilt)`, unless a different path has been specified in the
depot-wide preferences or using
[Overrides.toml](https://pkgdocs.julialang.org/dev/artifacts/#Overriding-artifact-locations).

This function sets the library path as a _preference_ associated with the currently active
environment. Changes will take effect after restarting the Julia session. See the [Preferences.jl
documentation](https://github.com/JuliaPackaging/Preferences.jl) for more information.
"""
function set_libceed_path!(path::AbstractString)
    set_preferences!(C.libCEED_jll, "libceed_path" => path; force=true)
    @info("""
          Setting the libCEED library path to $path.
          Restart the Julia session for changes to take effect.
          """)
end

function set_libceed_path!(sym::Symbol)
    if sym == :prebuilt
        set_preferences!(C.libCEED_jll, "libceed_path" => nothing; force=true)
        @info("""
              Using the prebuilt libCEED binary.
              Restart the Julia session for changes to take effect.
              """)
    elseif sym == :default
        delete_preferences!(Preferences.get_uuid(C.libCEED_jll), "libceed_path"; force=true)
        @info("""
              Deleting preference for libCEED library path.
              Restart the Julia session for changes to take effect.
              """)
    else
        error("set_libceed_path(::Symbol) must be called with :prebuilt or :default.")
    end
end

"""
    use_prebuilt_libceed!()

Indicates that the prebuilt version of the libCEED library (bundled with libCEED\\_jll) should be
used.

Equivalent to `set_libceed_path!(:prebuilt)`.
"""
use_prebuilt_libceed!() = set_libceed_path!(:prebuilt)

"""
    get_libceed_path()

Returns the path to the currently used libCEED library. A different libCEED library can be used by
calling [`set_libceed_path!`](@ref) or by using a depot-wide Overrides.toml file.
"""
get_libceed_path() = C.libCEED_jll.libceed_path

end # module

max_stars_repo_name: jedbrown/libCEED

max_issues_repo_name: jedbrown/libCEED

max_forks_repo_name: jedbrown/libCEED


text:
using PFKernels
using PFKernels.PowerSystem
using PFKernels.ForwardDiff
using LinearAlgebra

AT = PFKernels.AT

function mynorm(x)
    t1s{N} = ForwardDiff.Dual{Nothing,Float64, N} where N
    n = length(x)
    FT = t1s{n}
    V = Vector
    F = zeros(n)
    adF = V{FT}(undef, n)
    adx = V{FT}(undef, n)
    # ForwardDiff.seed!(adx, x, one(ForwardDiff.Partials{n,Float64}))
    ForwardDiff.seed!(adx, x, ForwardDiff.construct_seeds(ForwardDiff.Partials{n,Float64}))
    dF = AT(adF)
    dx = AT(adx)
    PFKernels.kernel!(PFKernels.backend, dF, dx)
    nrm2 = t1s{n}(0.0)
    for v in dF
      nrm2 += v
    end
    @show typeof(nrm2)
    @show nrm2
    return ForwardDiff.value(nrm2), ForwardDiff.partials(nrm2)
end

# Inject GPU kernels
PFKernels.include(joinpath(dirname(@__FILE__), "kernels.jl"))

x = rand(10)

nrm2 = x -> norm(x)^2
res = nrm2(x)
dres = ForwardDiff.gradient(nrm2, x)
@show x
F, dF = mynorm(x)
@assert(F ≈ res)
@assert(dF ≈ dres)


text:
"""
    IGRound

Abstract base type for dispatched InstaRound rounds.
"""
abstract type IGRound end

units = [
    "K", "M", "B", "t", "q", "Q", "s", "S", "o",
    "n", "d", "U", "D", "T", "Qt", "Qd", "Sd", "St",
    "O", "N", "v", "c"
]

unit_names = [
    "Thousand",
    "Million",
    "Billion",
    "Trillion",
    "Quadrillion",
    "Quintillion",
    "Sextillion",
    "Septillion",
    "Octillion",
    "Nonillion",
    "Decillion",
    "Undecillion",
    "Duodecillion",
    "Tredecillion",
    "Quattuordecillion",
    "Quindecillion",
    "Sexdecillion",
    "Septendecillion",
    "Octodecillion",
    "Novemdecillion",
    "Vigintillion",
    "Unvigintillion",
]


text:
using DifferentialEquations
using Plots
include("../src/ToyModelsMD.jl")
using .ToyModelsMD

mb = MullerBrown()
langevin = Langevin(1.0,1.0)
system = ThermostattedSystem(mb,langevin)

#fd(du,u,p,t) = fderivative(du,u,p,t;sys=system)
#gd(du,u,p,t) = gderivative(du,u,p,t;therm=langevin)

fd(du,u,p,t) = fderivative(du,u,system)
gd(du,u,p,t) = gderivative(du,u,langevin)

u0 = [0.63,0.03]
prob = SDEProblem(fd,gd,u0,(0.0,0.0015))
sol = solve(prob,EM(),dt=0.0001)
println(sol.u)
#plt = plot(sol)
#display(plt)


text:
module Contact

using LinearAlgebra
using RigidBodyDynamics.Spatial
using StaticArrays

# base types
export SoftContactModel,
    SoftContactState,
    SoftContactStateDeriv,
    ContactPoint,
    ContactEnvironment,
    HalfSpace3D

# interface functions
export num_states,
    location,
    contact_model,
    contact_dynamics!,
    reset!,
    point_inside,
    separation,
    detect_contact

# specific models
export HuntCrossleyModel,
    hunt_crossley_hertz

export ViscoelasticCoulombModel,
    ViscoelasticCoulombState,
    ViscoelasticCoulombStateDeriv

# defaults
export DefaultContactPoint,
    DefaultSoftContactState,
    DefaultSoftContactStateDeriv

# Normal force model/state/state derivative interface:
# * `normal_force(model, state, z::Number, ż::Number)`: return the normal force given penetration `z` and penetration velocity `ż`
# * `num_states(model)`: return the number of state variables associated with the model
# * `state(model, statevec, frame)`: create a new state associated with the model
# * `state_derivative(model, statederivvec, frame)`: create a new state derivative associated with the model
# * `reset!(state)`: resets the contact state
# * `dynamics!(state_deriv, model, state, fnormal)`: sets state_deriv given the current state and the normal force

# Friction model/state/state derivative interface:
# * `friction_force(model, state, fnormal::Number, tangential_velocity::FreeVector3D)`: return the friction force given normal force `fnormal` and tangential velocity `tangential_velocity`
# * `num_states(model)`: return the number of state variables associated with the model
# * `state(model, statevec, frame)`: create a new state associated with the model
# * `state_derivative(model, statederivvec, frame)`: create a new state derivative associated with the model
# * `reset!(state)`: resets the contact state
# * `dynamics!(state_deriv, model, state, ftangential)`: sets state_deriv given the current state and the tangential force force

## SoftContactModel and related types.
struct SoftContactModel{N, F}
    normal::N
    friction::F
end

struct SoftContactState{N, F}
    normal::N
    friction::F
end

struct SoftContactStateDeriv{N, F}
    normal::N
    friction::F
end

normal_force_model(model::SoftContactModel) = model.normal
friction_model(model::SoftContactModel) = model.friction

normal_force_state(state::SoftContactState) = state.normal
friction_state(state::SoftContactState) = state.friction

normal_force_state_deriv(deriv::SoftContactStateDeriv) = deriv.normal
friction_state_deriv(deriv::SoftContactStateDeriv) = deriv.friction

num_states(model::SoftContactModel) = num_states(normal_force_model(model)) + num_states(friction_model(model))

for (fun, returntype) in [(:state, :SoftContactState), (:state_derivative, :SoftContactStateDeriv)]
    @eval function $returntype(model::SoftContactModel, vec::AbstractVector, frame::CartesianFrame3D)
        nnormal = num_states(normal_force_model(model))
        nfriction = num_states(friction_model(model))
        vecnormal = view(vec, 1 : nnormal - 1)
        vecfriction = view(vec, 1 + nnormal : nnormal + nfriction)
        $returntype($fun(normal_force_model(model), vecnormal, frame), $fun(friction_model(model), vecfriction, frame))
    end
end

reset!(state::SoftContactState) = (reset!(normal_force_state(state)); reset!(friction_state(state)))
zero!(deriv::SoftContactStateDeriv) = (zero!(friction_state_deriv(deriv)); zero!(friction_state_deriv(deriv)))

## ContactPoint
mutable struct ContactPoint{T, M <: SoftContactModel}
    location::Point3D{SVector{3, T}}
    model::M
end
location(point::ContactPoint) = point.location
contact_model(point::ContactPoint) = point.model

function contact_dynamics!(state_deriv::SoftContactStateDeriv, state::SoftContactState,
        model::SoftContactModel, penetration::Number, velocity::FreeVector3D, normal::FreeVector3D)
    @boundscheck penetration >= 0 || error("penetration must be nonnegative")

    z = penetration
    ż = -dot(velocity, normal) # penetration velocity
    fnormal = max(normal_force(normal_force_model(model), normal_force_state(state), z, ż), 0)
    dynamics!(normal_force_state_deriv(state_deriv), normal_force_model(model), normal_force_state(state), fnormal)

    tangential_velocity = velocity + ż * normal
    ftangential = friction_force(friction_model(model), friction_state(state), fnormal, tangential_velocity)
    dynamics!(friction_state_deriv(state_deriv), friction_model(model), friction_state(state), ftangential)

    fnormal * normal + ftangential
end

## Models with no state
reset!(::Nothing) = nothing
zero!(::Nothing) = nothing

## Normal contact models
struct HuntCrossleyModel{T}
    # (2) in Marhefka, Orin, "A Compliant Contact Model with Nonlinear Damping for Simulation of Robotic Systems"
    k::T
    λ::T
    n::T
end

function hunt_crossley_hertz(; k = 50e3, α = 0.2)
    λ = 3/2 * α * k # (12) in Marhefka, Orin
    HuntCrossleyModel(k, λ, 3/2)
end

num_states(::HuntCrossleyModel) = 0
state(::HuntCrossleyModel, ::AbstractVector, ::CartesianFrame3D) = nothing
state_derivative(::HuntCrossleyModel, ::AbstractVector, ::CartesianFrame3D) = nothing

function normal_force(model::HuntCrossleyModel, ::Nothing, z, ż)
    zn = z^model.n
    f = model.λ * zn * ż + model.k * zn # (2) in Marhefka, Orin (note: z is penetration, returning repelling force)
end

dynamics!(ẋ::Nothing, model::HuntCrossleyModel, state::Nothing, fnormal::Number) = nothing

# Friction models
struct ViscoelasticCoulombModel{T}
    # See section 11.8 of Featherstone, "Rigid Body Dynamics Algorithms", 2008
    μ::T
    k::T
    b::T
end

mutable struct ViscoelasticCoulombState{V}
    tangential_displacement::FreeVector3D{V} # Use a 3-vector; technically only need one state for each tangential direction, but this is easier to work with.
end

mutable struct ViscoelasticCoulombStateDeriv{V}
    deriv::FreeVector3D{V}
end

num_states(::ViscoelasticCoulombModel) = 3
function state(::ViscoelasticCoulombModel, statevec::AbstractVector, frame::CartesianFrame3D)
    ViscoelasticCoulombState(FreeVector3D(frame, statevec))
end
function state_derivative(::ViscoelasticCoulombModel, statederivvec::AbstractVector, frame::CartesianFrame3D)
    ViscoelasticCoulombStateDeriv(FreeVector3D(frame, statederivvec))
end

reset!(state::ViscoelasticCoulombState) = fill!(state.tangential_displacement.v, 0)
zero!(deriv::ViscoelasticCoulombStateDeriv) = fill!(deriv.deriv.v, 0)

function friction_force(model::ViscoelasticCoulombModel, state::ViscoelasticCoulombState,
        fnormal, tangential_velocity::FreeVector3D)
    T = typeof(fnormal)
    μ = model.μ
    k = model.k
    b = model.b
    x = convert(FreeVector3D{SVector{3, T}}, state.tangential_displacement)
    v = tangential_velocity

    # compute friction force that would be needed to avoid slip
    fstick = -k * x - b * v

    # limit friction force to lie within friction cone
    fstick_norm² = dot(fstick, fstick)
    fstick_max_norm² = (μ * fnormal)^2
    ftangential = if fstick_norm² > fstick_max_norm²
        fstick * sqrt(fstick_max_norm² / fstick_norm²)
    else
        fstick
    end
end

function dynamics!(ẋ::ViscoelasticCoulombStateDeriv, model::ViscoelasticCoulombModel, state::ViscoelasticCoulombState, ftangential::FreeVector3D)
    k = model.k
    b = model.b
    x = state.tangential_displacement
    @framecheck ẋ.deriv.frame x.frame
    ẋ.deriv.v .= (-k .* x.v .- ftangential.v) ./ b
end

## VectorSegment: type of a view of a vector
const VectorSegment{T} = SubArray{T,1,Array{T, 1},Tuple{UnitRange{Int64}},true} # TODO: a bit too specific

const DefaultContactPoint{T} = ContactPoint{T,SoftContactModel{HuntCrossleyModel{T},ViscoelasticCoulombModel{T}}}
const DefaultSoftContactState{T} = SoftContactState{Nothing, ViscoelasticCoulombState{VectorSegment{T}}}
const DefaultSoftContactStateDeriv{T} = SoftContactStateDeriv{Nothing, ViscoelasticCoulombStateDeriv{VectorSegment{T}}}

# Contact detection
# TODO: should probably move this somewhere else

mutable struct HalfSpace3D{T}
    point::Point3D{SVector{3, T}}
    outward_normal::FreeVector3D{SVector{3, T}}

    function HalfSpace3D(point::Point3D{SVector{3, T}}, outward_normal::FreeVector3D{SVector{3, T}}) where {T}
        @framecheck point.frame outward_normal.frame
        new{T}(point, normalize(outward_normal))
    end
end

frame(halfspace::HalfSpace3D) = halfspace.point.frame

function HalfSpace3D(point::Point3D, outward_normal::FreeVector3D)
    T = promote_eltype(point, outward_normal)
    HalfSpace3D(convert(Point3D{SVector{3, T}}, point), convert(FreeVector3D{SVector{3, T}}, outward_normal))
end

Base.eltype(::Type{HalfSpace3D{T}}) where {T} = T
separation(halfspace::HalfSpace3D, p::Point3D) = dot(p - halfspace.point, halfspace.outward_normal)
point_inside(halfspace::HalfSpace3D, p::Point3D) = separation(halfspace, p) <= 0
detect_contact(halfspace::HalfSpace3D, p::Point3D) = separation(halfspace, p), halfspace.outward_normal

# ContactEnvironment
mutable struct ContactEnvironment{T}
    halfspaces::Vector{HalfSpace3D{T}}
    ContactEnvironment{T}() where {T} = new{T}(HalfSpace3D{T}[])
end

Base.push!(environment::ContactEnvironment, halfspace::HalfSpace3D) = push!(environment.halfspaces, halfspace)
Base.length(environment::ContactEnvironment) = length(environment.halfspaces)

end


text:
using MPI
using Test
using LinearAlgebra
using Random
using KernelAbstractions, StaticArrays

using ClimateMachine
using ClimateMachine.SystemSolvers
using ClimateMachine.SystemSolvers:
    band_lu_kernel!,
    band_forward_kernel!,
    band_back_kernel!,
    DGColumnBandedMatrix

import ClimateMachine.MPIStateArrays: array_device

ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()

function band_to_full(B, p, q)
    _, n = size(B)

    A = similar(B, n, n) # assume square
    fill!(A, 0)
    for j in 1:n, i in max(1, j - q):min(j + p, n)
        A[i, j] = B[q + i - j + 1, j]
    end
    A
end

function run_columnwiselu_test(FT, N)

    Nq = N .+ 1
    Nq1 = Nq[1]
    Nq2 = Nq[2]
    Nqv = Nq[3]
    nstate = 3
    nvertelem = 5
    nhorzelem = 4
    eband = 2

    m = n = Nqv * nstate * nvertelem
    p = q = Nqv * nstate * eband - 1

    Random.seed!(1234)
    AB = rand(FT, Nq1, Nq2, p + q + 1, n, nhorzelem)

    AB[:, :, q + 1, :, :] .+= 10 # Make A's diagonally dominate

    Random.seed!(5678)
    b = rand(FT, Nq1, Nq2, Nqv, nstate, nvertelem, nhorzelem)
    x = similar(b)

    perm = (4, 3, 5, 1, 2, 6)
    bp = reshape(PermutedDimsArray(b, perm), n, Nq1, Nq2, nhorzelem)
    xp = reshape(PermutedDimsArray(x, perm), n, Nq1, Nq2, nhorzelem)

    d_F = ArrayType(AB)
    d_F = DGColumnBandedMatrix{
        3,
        N,
        nstate,
        nhorzelem,
        nvertelem,
        eband,
        false,
        typeof(d_F),
    }(
        d_F,
    )

    groupsize = (Nq1, Nq2)
    ndrange = (Nq1, Nq2, nhorzelem)

    event = Event(array_device(d_F.data))
    event = band_lu_kernel!(array_device(d_F.data), groupsize, ndrange)(
        d_F,
        dependencies = (event,),
    )
    wait(array_device(d_F.data), event)

    F = Array(d_F.data)

    for h in 1:nhorzelem, j in 1:Nq2, i in 1:Nq1
        B = AB[i, j, :, :, h]
        G = band_to_full(B, p, q)
        GLU = lu!(G, Val(false))

        H = band_to_full(F[i, j, :, :, h], p, q)

        @assert H ≈ G

        xp[:, i, j, h] .= GLU \ bp[:, i, j, h]
    end

    b = reshape(b, Nq1 * Nq2 * Nqv, nstate, nvertelem * nhorzelem)
    x = reshape(x, Nq1 * Nq2 * Nqv, nstate, nvertelem * nhorzelem)

    d_x = ArrayType(b)

    event = Event(array_device(d_x))
    event = band_forward_kernel!(array_device(d_x), groupsize, ndrange)(
        d_x,
        d_F,
        dependencies = (event,),
    )

    event = band_back_kernel!(array_device(d_x), groupsize, ndrange)(
        d_x,
        d_F,
        dependencies = (event,),
    )
    wait(array_device(d_x), event)

    result = x ≈ Array(d_x)
    return result
end

@testset "Columnwise LU test" begin
    for FT in (Float64, Float32)
        for N in ((1, 1, 1), (1, 1, 2), (2, 2, 1))
            result = run_columnwiselu_test(FT, N)
            @test result
        end
    end
end

max_forks_repo_name: jm-c/CLIMA


text:
export Wait

"""
    Wait{N, T} <: TrivialGate{N}
    Wait{N}(t)

Wait the experimental signals for time `t` (empty run).
"""
struct Wait{N, T} <: TrivialGate{N}
    t::T
    Wait{N}(t::T) where {N,T} = new{N, T}(t)
end
Yao.print_block(io::IO, d::Wait) = print(io, "Wait → $(d.t)")

export EchoBlock
struct EchoBlock{N,OT<:IO} <: PrimitiveBlock{N}
    sym::Symbol
    io::OT
end

EchoBlock(nbits::Int, sym::Symbol) = EchoBlock{nbits, typeof(stdout)}(sym, stdout)
EchoBlock(nbits::Int) = EchoBlock(nbits, :ECHO)
EchoBlock() = nbits->EchoBlock(nbits)

_apply!(reg::AbstractRegister, ec::EchoBlock) = (println(ec.io, "apply!(::$(typeof(reg)), $(ec.sym))"); reg)
Yao.mat(::Type{T}, ec::EchoBlock{N}) where {T,N} = (println(ec.io, "mat(::Type{$T}, $(ec.sym))"); IMatrix{1<<N}())
Yao.ishermitian(ec::EchoBlock{N}) where N = (println(ec.io, "ishermitian($(ec.sym))"); true)
Yao.isunitary(ec::EchoBlock{N}) where N = (println(ec.io, "isunitary($(ec.sym))"); true)
Yao.isreflexive(ec::EchoBlock{N}) where N = (println(ec.io, "isreflexive($(ec.sym))"); true)
Yao.getiparams(ec::EchoBlock{N}) where N = (println(ec.io, "getiparams($(ec.sym))");())
Yao.print_block(io::IO, ec::EchoBlock) = print(io, "EchoBlock($(ec.sym))")

max_stars_repo_name: QuantumBFS/YaoExtensions

max_issues_repo_name: QuantumBFS/YaoExtensions

max_forks_repo_name: QuantumBFS/YaoExtensions


text:
# Exercise 13-5
function cumulsum(t)
    res = Array{Float64, 1}()
    for ele in t
        if ele == t[begin]
            push!(res, ele)
        else
            push!(res, ele + res[end])
        end
    end
    return res
end

function histogram(s)
    d = Dict()
    for c in s
        d[c] = get(d, c, 0) + 1
    end
    d
end

function choosefromhist(hist)
    chars, counts = Array{Char, 1}(), Array{Float64, 1}()
    for (k, v) in hist
        push!(chars, k)
        push!(counts, v)
    end
    prob = cumulsum(counts / sum(counts))
    rv = rand()
    for (i, p) in enumerate(prob)
        if rv < p
            return chars[i]
        end
    end
end

hist = histogram("abbbcccccddddddddeeeeeeeeeeeee")
# println(choosefromhist(hist))

# Exercise 13-9
function ex13_9(hist)
    charfreq = collect(hist)
    sort!(charfreq, by = x -> x.second)
    logf, logr = Array{Float64, 1}(), Array{Float64, 1}()
    for (r, pair) in enumerate(charfreq)
        char, f = pair.first, pair.second
        println("$char with f:$f, log f: $(log(f)), log r: $(log(r))")
        push!(logf, log(f))
        push!(logr, log(r))
    end
    return logf, logr
end

# logf, logr = ex13_9(hist)
# using Plots
# plot(logr, logf)


text:
module TestSKPreprocessing

using Test
using AutoMLPipeline.SKPreprocessors
using AutoMLPipeline.Utils
using Statistics
using DataFrames

const IRIS = getiris()
const X = IRIS[:,1:4] |> DataFrame
const Y = IRIS[:,5] |> Vector

# "KernelCenterer","MissingIndicator","KBinsDiscretizer","OneHotEncoder",
const preprocessors = [
     "DictionaryLearning", "FactorAnalysis", "FastICA", "IncrementalPCA",
     "KernelPCA", "LatentDirichletAllocation", "MiniBatchDictionaryLearning",
     "MiniBatchSparsePCA", "NMF", "PCA",
     "TruncatedSVD",
     "VarianceThreshold",
     "SimpleImputer",
     "Binarizer", "FunctionTransformer",
     "MultiLabelBinarizer", "MaxAbsScaler", "MinMaxScaler", "Normalizer",
     "OrdinalEncoder", "PolynomialFeatures", "PowerTransformer",
     "QuantileTransformer", "RobustScaler", "StandardScaler"
 ]

function fit_test(preproc::String,in::DataFrame,out::Vector)
	_preproc=SKPreprocessor(Dict(:preprocessor=>preproc))
	fit!(_preproc,in,out)
	@test _preproc.model != Dict()
	return _preproc
end

function transform_test(preproc::String,in::DataFrame,out::Vector)
	_preproc=SKPreprocessor(Dict(:preprocessor=>preproc))
	fit!(_preproc,in,out)
	res = transform!(_preproc,in)
	@test size(res)[1] == size(out)[1]
end

@testset "scikit preprocessors fit test" begin
    for cl in preprocessors
	println(cl)
	fit_test(cl,X,Y)
    end
end

@testset "scikit preprocessors transform test" begin
    for cl in preprocessors
	println(cl)
	transform_test(cl,X,Y)
    end
end

function skptest()
    features = X
    labels = Y

    pca = SKPreprocessor(Dict(:preprocessor=>"PCA",:impl_args=>Dict(:n_components=>3)))
    fit!(pca,features)
    @test transform!(pca,features) |> x->size(x,2) == 3

    svd = SKPreprocessor(Dict(:preprocessor=>"TruncatedSVD",:impl_args=>Dict(:n_components=>2)))
    fit!(svd,features)
    @test transform!(svd,features) |> x->size(x,2) == 2

    ica = SKPreprocessor(Dict(:preprocessor=>"FastICA",:impl_args=>Dict(:n_components=>2)))
    fit!(ica,features)
    @test transform!(ica,features) |> x->size(x,2) == 2

    stdsc = SKPreprocessor(Dict(:preprocessor=>"StandardScaler",:impl_args=>Dict()))
    fit!(stdsc,features)
    @test abs(mean(transform!(stdsc,features) |> Matrix)) < 0.00001

    minmax = SKPreprocessor(Dict(:preprocessor=>"MinMaxScaler",:impl_args=>Dict()))
    fit!(minmax,features)
    @test mean(transform!(minmax,features) |> Matrix) ≈ 0.4486931104833648
end
@testset "scikit preprocessor fit/transform test with real data" begin
    skptest()
end

end


text:
import Random

include("City.jl")

# abstract type MapModification end
#
# struct AddCity <: MapModification
#     city::City
# end
#
# struct DeleteCity <: MapModification
#     city::City
# end

mutable struct Map

    cities::Vector{City}

    transitions::Dict{City, Union{Vector{CityIndex}, Nothing}}

end

struct RandomCreation end

function Map()
    return Map(Vector{City}(), Dict{City, Union{Vector{CityIndex}, Nothing}}())
end

#Is it better for Space complexity ? Time Complexity ? to not recrate space memory every time

# function Map(nbrVille)
#     return Map(Vector{City}(undef, nbrVille), Vector{Vector{Float64}}( undef, nbrVille ))
# end

# TODO : Find how to use the macro correctly

# TODO :
# - Precise  which city are connected to each other ! To make a perfect map
# - Give ptentially a set distance to each other city, like a Dict{City, distance}
#   Then if there isn't that : the distance are calculate by x^2 + y^2
#   else : the fixed one are used

function city!(map::Map, x::Float64, y::Float64, transition::Union{Vector{Int64}, Nothing} = nothing,  pheroInit::Float64=2.)

    nbrVille::Int64 = length(map.cities)
    city::City = City(x, y, nbrVille+=1)

    if transition != nothing
        transition = [CityIndex(index) for index in transition]
    end
    #I'm there right now !

    # ways::Dict{City, Way} = Dict(map.cities[ind] => Way(distance(city, map.cities[ind]), pheroInit) for ind = 1:nbrVille-1)

    push!(map.cities, city)
    push!(map.transitions, city => transition)

    # push!(map.listModification, AddCity(city))

    return city
end

#TODO : Find how tu use the macro correctly

function Map(nbrVille::Int64, borneX::Int64 = 50, borneY::Int64 = 50)
    m = Map()
    for ind = 1:nbrVille
        cityX = borneX*rand()
        cityY = borneY*rand()
        city!(m, cityX, cityY)
    end
    return m
end

function createWays(map::Map, pheroInit::Float64)
    listWays::Dict{CityIndex, Dict{CityIndex, Way}} = Dict{CityIndex, Dict{CityIndex, Way}}()
    for city in map.cities
        nextCitiesIndex = map.transitions[city]
        if nextCitiesIndex != nothing
            push!(listWays, city.index => Dict(cityIndex => Way(distance(city, map.cities[cityIndex.value]), pheroInit) for cityIndex in (map(nextCityIndex-> nextCityIndex != city.index, nextCitiesIndex))))
        else
            listCities = map.cities[:]
            deleteat!(listCities, city.index.value)
            ways = Dict(nextCity.index => Way(distance(city, nextCity), pheroInit) for nextCity in listCities)
            push!(listWays, city.index => ways)
        end
    end
    return listWays
end


text:
#
# Function that reads atom information from PDB or mmCIF files
#
function read_atom(
    record::String;
    mmCIF::Bool = false,
    mmCIF_fields::Indexes_mmCIF_fields = Indexes_mmCIF_fields(),
)

    atom = Atom()

    if length(record) < 6
        return nothing
    end

    if record[1:4] == "ATOM" || record[1:6] == "HETATM"

        if !mmCIF

            atom.name = strip(record[13:16])

            atom.resname = strip(record[17:21])
            resname = alternate_conformation(atom)
            if resname == nothing
                return nothing
            else
                atom.resname = resname
            end

            atom.chain = strip(record[22:22])
            if atom.chain == " "
                atom.chain = "0"
            end

            atom.index = 1
            try
                atom.index_pdb = parse_int(record[7:11])
            catch
                atom.index_pdb = 0
            end
            atom.resnum = parse_int(record[23:26])
            atom.x = parse(Float64, record[31:38])
            atom.y = parse(Float64, record[39:46])
            atom.z = parse(Float64, record[47:54])
            try
                atom.beta = parse(Float64, record[61:66])
            catch
                atom.beta = 0.0
            end
            try
                atom.occup = parse(Float64, record[56:60])
            catch
                atom.occup = 0.0
            end
            atom.model = 1
            try
                atom.segname = strip(record[73:76])
                if length(atom.segname) == 0
                    atom.segname = "-"
                end
            catch
                atom.segname = "-"
            end

        else # if mmCIF

            mmcif_data = split(record)
            atom.index = 1
            try
                atom.index_pdb = parse(Int, mmcif_data[mmCIF_fields.index])
            catch
                atom.index_pdb = 0
            end
            atom.name = mmcif_data[mmCIF_fields.name]
            atom.resname = mmcif_data[mmCIF_fields.resname]
            resname = alternate_conformation(atom)
            if resname == nothing
                return nothing
            else
                atom.resname = resname
            end
            atom.chain = mmcif_data[mmCIF_fields.chain]
            try
                atom.resnum = parse(Int, mmcif_data[mmCIF_fields.resnum])
            catch
                atom.resnum = 0
            end
            try
                atom.segname = mmcif_data[mmCIF_fields.segname]
            catch
                atom.segname = ""
            end
            atom.x = parse(Float64, mmcif_data[mmCIF_fields.x])
            atom.y = parse(Float64, mmcif_data[mmCIF_fields.y])
            atom.z = parse(Float64, mmcif_data[mmCIF_fields.z])
            atom.beta = parse(Float64, mmcif_data[mmCIF_fields.beta])
            atom.occup = parse(Float64, mmcif_data[mmCIF_fields.occup])
            atom.model = 1

        end

        return atom

    else

        return nothing

    end

end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
module DrivenJumpProblems

export DrivenJumpProblem

using DiffEqJump
using CommonSolve
using PWS: Trajectory

"""
    IdentityMap()

A simple helper type that always returns `i` when indexed as `IdentityMap()[i]`.
"""
struct IdentityMap end
@inline Base.getindex(x::IdentityMap, i::Integer) = i

mutable struct TrajectoryCallback{T, IndexMap}
    traj::T
    index::Int
    index_map::IndexMap
end

TrajectoryCallback(traj, index_map = IdentityMap()) = TrajectoryCallback(traj, 1, index_map)

function (tc::TrajectoryCallback)(integrator::DiffEqBase.DEIntegrator) # affect!
    traj = tc.traj
    tc.index = min(tc.index + 1, length(traj.t))
    cond_u = traj.u[tc.index]
    for i in eachindex(cond_u)
        integrator.u[tc.index_map[i]] = cond_u[i]
    end
    # it is important to call this to properly update reaction rates
    DiffEqJump.reset_aggregated_jumps!(integrator, nothing, integrator.cb, update_jump_params=false)
    nothing
end

function (tc::TrajectoryCallback)(u, t::Real, i::DiffEqBase.DEIntegrator)::Bool # condition
    @inbounds tcb = tc.traj.t[tc.index]
    while tc.index < length(tc.traj.t) && t > tcb
        tc.index += 1
        @inbounds tcb = tc.traj.t[tc.index]
    end
    t == tcb
end

"""
    DrivenJumpProblem(jump_problem, driving_trajectory; index_map = IdentityMap(), save_jumps=false)

Create a `DrivenJumpProblem` from a base `jump_problem` and a `driving_trajectory`.

The optional argument `index_map` specifies how the components from `driving_trajectory`
map onto the components of the `jump_problem`. The default `IdentityMap` maps the N components
of the `driving_trajectory` onto the first N components of the jump problem, leaving the
remaining components of the jump problem unaltered by the driving trajectory.
"""
struct DrivenJumpProblem{Prob,Cb}
    prob::Prob
    callback::Cb

    function DrivenJumpProblem(jump_problem::JP, driving_trajectory; index_map = IdentityMap(), save_jumps=false) where {JP}
        tcb = TrajectoryCallback(Trajectory(driving_trajectory), index_map)
        callback = DiscreteCallback(tcb, tcb, save_positions=(false, save_jumps))
        new{JP, typeof(callback)}(jump_problem, callback)
    end
end

function CommonSolve.init(prob::DrivenJumpProblem; kwargs...)
    prob.callback.condition.index = 1
    tstops = prob.callback.condition.traj.t
    from = searchsortedfirst(tstops, prob.prob.prob.tspan[1])
    to = searchsortedlast(tstops, prob.prob.prob.tspan[2])
    DiffEqBase.init(prob.prob, SSAStepper(), callback=prob.callback, tstops=tstops[from:to], save_start=false; kwargs...)
end

function CommonSolve.solve(prob::DrivenJumpProblem; kwargs...)
    integrator = init(prob; kwargs...)
    solve!(integrator)
    integrator.sol
end

Base.summary(io::IO, prob::DrivenJumpProblem) = string(DiffEqBase.parameterless_type(prob)," with problem ",DiffEqBase.parameterless_type(prob.prob))
function Base.show(io::IO, mime::MIME"text/plain", A::DrivenJumpProblem)
    println(io,summary(A))
end

# We define a new jump aggregator to deal with a jump process that is driven
# by an input signal. The external driving signal is supplied as a trajectory of
# jump times t_i and values u_i. The reaction rates are determined by mass
# action kinetics.

function time_to_next_jump(remaining_input_trajectory)
    t_jump = zero(Float64)
    r = randexp()

    while true
        u_s, t_s = current_values(remaining_input_trajectory)
        u_s_next, t_s_next = iterate(remaining_input_trajectory)
        sumrate = eval_sumrate(u_s, u_x)
        next_sumrate = eval_sumrate(u_s_next, u_x)

        Δt = t_s_next - t_s
        if Δt * (sumrate + next_sumrate) / 2 > r
            # the next jump does not occur before t_s_next
            r -= Δt * (sumrate + next_sumrate) / 2
            continue
        else
            # now we have to find the precise jump time between t_s and t_s_next
            # which is the solution to the equation r = ∫_[t, t+Δt] sumrate(t) dt
            # when solved for Δt
            sumrate_slope = (next_sumrate - sumrate) / Δt
            t_jump = t_s + solve_quadratic(0.5 * sumrate_slope, sumrate, -r)[2]
            break
        end
    end

    t_jump
end

# ================== UTILITY FUNCTIONS ==================

# computes ab-cd in a numerically stable way
@inline function diff_of_products(a, b, c, d)
    w = d * c
    e = fma(-d, c, w)
    f = fma(a, b, -w)
    f + e
end

# computes real roots of quadratic equation ax^2 + bx + c = 0, assuming discriminant is positive
@inline function solve_quadratic(a, b, c)
    discriminant = sqrt(diff_of_products(b, b, 4 * a, c))
    q = -0.5 * (b + copysign(discriminant, b))
    q / a, c / q
end

end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct DeterministicInternal <: ForecastInternal
        label::String
        resolution::Dates.Period
        initial_time::Dates.DateTime
        time_series_uuid::UUIDs.UUID
        horizon::Int
        internal::InfrastructureSystemsInternal
    end

A deterministic forecast for a particular data field in a Component.

# Arguments
- `label::String`: label of component parameter forecasted
- `resolution::Dates.Period`
- `initial_time::Dates.DateTime`: forecast availability time
- `time_series_uuid::UUIDs.UUID`: reference to time series data; timestamp - scalingfactor
- `horizon::Int`: length of this forecast
- `internal::InfrastructureSystemsInternal`
"""
mutable struct DeterministicInternal <: ForecastInternal
    "label of component parameter forecasted"
    label::String
    resolution::Dates.Period
    "forecast availability time"
    initial_time::Dates.DateTime
    "reference to time series data; timestamp - scalingfactor"
    time_series_uuid::UUIDs.UUID
    "length of this forecast"
    horizon::Int
    internal::InfrastructureSystemsInternal
end

function DeterministicInternal(label, resolution, initial_time, time_series_uuid, horizon, )
    DeterministicInternal(label, resolution, initial_time, time_series_uuid, horizon, InfrastructureSystemsInternal(), )
end

function DeterministicInternal(; label, resolution, initial_time, time_series_uuid, horizon, )
    DeterministicInternal(label, resolution, initial_time, time_series_uuid, horizon, )
end

"""Get DeterministicInternal label."""
get_label(value::DeterministicInternal) = value.label
"""Get DeterministicInternal resolution."""
get_resolution(value::DeterministicInternal) = value.resolution
"""Get DeterministicInternal initial_time."""
get_initial_time(value::DeterministicInternal) = value.initial_time
"""Get DeterministicInternal time_series_uuid."""
get_time_series_uuid(value::DeterministicInternal) = value.time_series_uuid
"""Get DeterministicInternal horizon."""
get_horizon(value::DeterministicInternal) = value.horizon
"""Get DeterministicInternal internal."""
get_internal(value::DeterministicInternal) = value.internal


text:
__precompile__(true)

module YAML

import Base: isempty, length, show

import Base: iterate

using Base64: base64decode
using Dates
using Printf

include("scanner.jl")
include("parser.jl")
include("composer.jl")
include("constructor.jl")
include("writer.jl") # write Julia dictionaries to YAML files

# const _constructor = Union{Nothing, Dict}

function load(ts::TokenStream, constructor::Constructor = SafeConstructor())
    events = EventStream(ts)
    node = compose(events)
    return construct_document(constructor, node)
end

function load(input::IO, constructor::Constructor = SafeConstructor())
    load(TokenStream(input), constructor)
end

load(input, constructors::Dict, multi_constructors::Dict = Dict()) = load(input, SafeConstructor(constructors, multi_constructors))

mutable struct YAMLDocIterator
    input::IO
    ts::TokenStream
    constructor::Constructor
    next_doc

    function YAMLDocIterator(input::IO, constructor::Constructor = SafeConstructor())
        it = new(input, TokenStream(input), constructor, nothing)
        it.next_doc = eof(it.input) ? nothing : load(it.ts, it.constructor)
        return it
    end
end

YAMLDocIterator(input::IO, constructors::Dict, multi_constructors::Dict = Dict()) = YAMLDocIterator(input, SafeConstructor(constructors, multi_constructors))

# Old iteration protocol:

start(it::YAMLDocIterator) = nothing

function next(it::YAMLDocIterator, state)
    doc = it.next_doc
    if eof(it.input)
        it.next_doc = nothing
    else
        reset!(it.ts)
        it.next_doc = load(it.ts, it.constructor)
    end
    return doc, nothing
end

done(it::YAMLDocIterator, state) = it.next_doc === nothing

# 0.7 iteration protocol:

iterate(it::YAMLDocIterator) = next(it, start(it))
iterate(it::YAMLDocIterator, s) = done(it, s) ? nothing : next(it, s)

function load_all(input::IO, args...)
    YAMLDocIterator(input, args...)
end

function load(input::AbstractString, constructor::Constructor = SafeConstructor())
    load(IOBuffer(input), constructor)
end

function load_all(input::AbstractString, args...)
    load_all(IOBuffer(input), args...)
end

function load_file(filename::AbstractString, args...)
    open(filename, "r") do input
        load(input, args...)
    end
end

function load_all_file(filename::AbstractString, args...)
    open(filename, "r") do input
        load_all(input, args...)
    end
end

end  # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "pugixml"
version = v"1.11.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://github.com/zeux/pugixml/releases/download/v1.11/pugixml-1.11.tar.gz", "26913d3e63b9c07431401cf826df17ed832a20d19333d043991e611d23beaa2c")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/pugixml*
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_SHARED_LIBS=ON \
    -S ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libpugixml", :libpugixml)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is 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.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using TensorFlow
import TensorFlow.API: transpose_, pack

type OpsBlock
    inputs::Tuple{Vararg{Tensor}}
    oplist::Vector{Function}
    opargs::Vector{Any}
    outs::Vector{Tensor}
    output::Tensor
end

function OpsBlock(inputs::Tuple{Vararg{Tensor}}, oplist::Vector{Function}; opargs::Vector{Any}=Any[])
    outs = map(f -> transpose_(f(inputs..., opargs ...)), oplist) #FIXME: remove the double transpose, but this is due to pack packing to the wrong dimension
    output = transpose_(pack(Tensor(outs)))

    OpsBlock(inputs, oplist,  opargs, outs, output)
end

function num_ops(ops::OpsBlock)
    length(ops.oplist)
end

function out(ops::OpsBlock)
    ops.output
end


text:
# Practical Machine learning
# Bayesian learning - Naive Bayes example
# Chapter 9
module NaiveBayes

export NBModel,
       MultinomialNB,
       GaussianNB,
       fit,
       predict,
       predict_proba

include("nbtypes.jl")
include("core.jl")

end


text:
const HTML_TEMPLATES = LittleDict{Symbol,BlockTemplate}(e.opening => e for e in [
   BlockTemplate(:COMMENT, :COMMENT_OPEN, :COMMENT_CLOSE),
   BlockTemplate(:SCRIPT,  :SCRIPT_OPEN,  :SCRIPT_CLOSE ),
   BlockTemplate(:DBB,     :DBB_OPEN,     :DBB_CLOSE)
   ])


text:
module TopOptProblems

using JuAFEM, StaticArrays, LinearAlgebra, SparseArrays
using InpParser, Setfield, CuArrays, CUDAnative, GPUUtils
using CUDAdrv: CUDAdrv
#using Makie
#using GeometryTypes

import JuAFEM: assemble!
import GPUUtils: whichdevice

abstract type AbstractTopOptProblem end

const dev = CUDAdrv.device()
const ctx = CUDAdrv.CuContext(dev)

#include("gpu_utils.jl")
include("utils.jl")
include("grids.jl")
include("metadata.jl")
include("problem_types.jl")
include("matrices_and_vectors.jl")
include("penalties.jl")
include("assemble.jl")
#include("makie.jl")

export PointLoadCantilever, HalfMBB, LBeam, TieBeam, InpStiffness, StiffnessTopOptProblem, AbstractTopOptProblem, GlobalFEAInfo, ElementFEAInfo, YoungsModulus, assemble, assemble_f!, RaggedArray, ElementMatrix, rawmatrix, bcmatrix

end # module


text:
# # Postprocessing
#
# ![](heat_square_fluxes.png)

# ## Introduction
#
# After running a simulation, we usually want to visualize the results in different ways.
# The `L2Projector` and the `PointEvalHandler` build a pipeline for doing so. With the `L2Projector`,
# integration point quantities can be projected to the nodes. The `PointEvalHandler` enables evaluation of
# the finite element approximated function in any coordinate in the domain. Thus with the combination of both functionalities,
# both nodal quantities and integration point quantities can be evaluated in any coordinate, allowing for example
# cut-planes through 3D structures or cut-lines through 2D-structures.
#
# This example continues from the Heat equation example, where the temperature field was
# determined on a square domain. In this example, we first compute the heat flux in each
# integration point (based on the solved temperature field) and then we do an L2-projection
# of the fluxes to the nodes of the mesh. By doing this, we can more easily visualize
# integration points quantities. Finally, we visualize the temperature field and the heat fluxes along a cut-line.
#
# The L2-projection is defined as follows: Find projection ``q(\boldsymbol{x}) \in L_2(\Omega)`` such that
# ```math
# \int v q \ \mathrm{d}\Omega = \int v d \ \mathrm{d}\Omega \quad \forall v \in L_2(\Omega),
# ```
# where ``d`` is the quadrature data to project. Since the flux is a vector the projection function
# will be solved with multiple right hand sides, e.g. with ``d = q_x`` and ``d = q_y`` for this 2D problem.
#
# Ferrite has functionality for doing much of this automatically, as displayed in the code below.
# In particular [`L2Projector`](@ref) for assembling the left hand side, and
# [`project`](@ref) for assembling the right hand sides and solving for the projection.

# ## Implementation
#
# Start by simply running the Heat equation example to solve the problem
include("heat_equation.jl");

# Next we define a function that computes the heat flux for each integration point in the domain.
# Fourier's law is adopted, where the conductivity tensor is assumed to be isotropic with unit
# conductivity ``\lambda = 1 ⇒ q = - \nabla u``, where ``u`` is the temperature.
function compute_heat_fluxes(cellvalues::CellScalarValues{dim,T}, dh::DofHandler, a) where {dim,T}

    n = getnbasefunctions(cellvalues)
    cell_dofs = zeros(Int, n)
    nqp = getnquadpoints(cellvalues)

    ## Allocate storage for the fluxes to store
    q = [Vec{2,T}[] for _ in 1:getncells(dh.grid)]

    for (cell_num, cell) in enumerate(CellIterator(dh))
        q_cell = q[cell_num]
        celldofs!(cell_dofs, dh, cell_num)
        aᵉ = a[cell_dofs]
        reinit!(cellvalues, cell)

        for q_point in 1:nqp
            q_qp = - function_gradient(cellvalues, q_point, aᵉ)
            push!(q_cell, q_qp)
        end
    end
    return q
end
#md nothing # hide

# Now call the function to get all the fluxes.
q_gp = compute_heat_fluxes(cellvalues, dh, u);

# Next, create an [`L2Projector`](@ref) using the same interpolation as was used to approximate the
# temperature field. On instantiation, the projector assembles the coefficient matrix `M` and
# computes the Cholesky factorization of it. By doing so, the projector can be reused without
# having to invert `M` every time.
projector = L2Projector(ip, grid);

# Project the integration point values to the nodal values
q_projected = project(projector, q_gp, qr; project_to_nodes=false); # TODO: this should be default.

# ## Exporting to VTK
# To visualize the heat flux, we export the projected field `q_projected`
# to a VTK-file, which can be viewed in e.g. [ParaView](https://www.paraview.org/).
vtk_grid("heat_equation_flux", grid) do vtk
    vtk_point_data(vtk, projector, q_projected, "q")
end;

# ## Point Evaluation
# ![](heat_square_pointevaluation.png)

# Consider a cut-line through the domain, like the black line in the figure above.
# We will evaluate the temperature and the heat flux distribution along a horizontal line.
points = [Vec((x, 0.75)) for x in range(-1.0, 1.0, length=101)];

# First, we need to generate a `PointEvalHandler`. This will find and store the cells
# containing the input points.
ph = PointEvalHandler(grid, points);

# After the L2-Projection, the heat fluxes `q_projected` are stored in the DoF-ordering
# determined by the projector's internal DoFHandler, so to evalute the flux `q` at our points
# we give the `PointEvalHandler`, the `L2Projector` and the values `q_projected`.
q_points = get_point_values(ph, projector, q_projected);

# We can also extract the field values, here the temperature, right away from the result
# vector of the simulation, that is stored in `u`. These values are stored in the order of
# our initial DofHandler so the input is not the `PointEvalHandler`, the original `DofHandler`,
# the dof-vector `u`, and (optionally for single-field problems) the name of the field.
# from the `L2Projection`, the values are stored in the order of the degrees of freedom.
u_points = Ferrite.get_point_values(ph, dh, u, :u);

# Now, we can plot the temperature and flux values with the help of any plotting library, e.g. Plots.jl.
# To do this, we need to import the package:
import Plots

# Firstly, we are going to plot the temperature values along the given line.
Plots.plot(getindex.(points,1), u_points, xlabel="x (coordinate)", ylabel="u (temperature)", label=nothing)

# Secondly, the horizontal heat flux (i.e. the first component of the heat flux vector) is plotted.
Plots.plot(getindex.(points,1), getindex.(q_points,1), xlabel="x (coordinate)", ylabel="q_x (flux in x-direction)", label=nothing)

#md # ## [Plain Program](@id postprocessing-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [postprocessing.jl](postprocessing.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
# packages -
using LinearAlgebra

# my codes -
include("Expa.jl")
include("Stoichiometric.jl")
include("Parser.jl")


text:
# (C) Copyright IBM 2017, 2018.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.
#
# Modifications copyright Spencer King, 2020

using Qiskit

Qiskit.init()

# making first circuit: bell state
qc1 = Qiskit.QuantumCircuit(2, 2)
qc1.h(0)
qc1.cx(0, 1)
qc1.measure([0,1], [0,1])

# making another circuit: superpositions
qc2 = Qiskit.QuantumCircuit(2, 2)
qc2.h([0,1])
qc2.measure([0,1], [0,1])

# setting up the backend
print("(BasicAER Backends)")
print(Qiskit.get_basic_backends())

# running the job
job_sim = Qiskit.execute([qc1, qc2], Qiskit.get_basic_backend("qasm_simulator"))
sim_result = job_sim.result()

# Show the results
print(sim_result.get_counts(qc1))
print(sim_result.get_counts(qc2))


text:
using Janitor
using Test

@testset "Janitor.jl" begin
    # Write your tests here.
end


text:
"""
Snape's Riddle - Author JK Rowling
Danger lies before you, while safety lies behind,
Two of us will help you, whichever you would find,
One among us seven will let you move ahead,
Another will transport the drinker back instead,
Two among our number hold only nettle wine,
Three of us are killers, waiting hidden in line.
Choose, unless you wish to stay here for evermore,
To help you in your choice, we give you these clues four:
First, however slyly the poison tries to hide
You will always find some on nettle wine’s left side;
Second, different are those who stand at either end,
But if you would move onwards, neither is your friend;
Third, as you see clearly, all are different size,
Neither dwarf nor giant holds death in their insides;
Fourth, the second left and the second on the right
Are twins once you taste them, though different at first sight.
"""
using AbstractLogic

# Codes:
# 1 is netted wine - NW
# 2 move ahead     - MA
# 3 move back      - MB
# 4 poison         - PO

# Seven Potions a,b,c,d,e,f,g each have 1 of 4 possible values
# a is left of b which is left of c etc.
hpset = logicalparse("a, b, c, d, e, f, g  ∈  NW, MA, MB, PO")

# Two among our number hold only nettle wine, #code 1
hpset = logicalparse("{{i}} == 'NW' {{2}}", logicset=hpset); showfeasible(hpset)
# Equivalent to logicalparse(["a|b|c|d|e|f|g == 1 {2}"], logicset=hpset)

# One among us seven will let you move ahead, #code 2
hpset = logicalparse("{{i}} == 'MA' {{1}}", logicset=hpset); showfeasible(hpset)
#hpset = logicalparse(["a|b|c|d|e|f|g == 2 {1}"], logicset=hpset)

#Another will transport the drinker back instead, #code 3
hpset = logicalparse("{{i}} == 'MB' {{1}}", logicset=hpset); showfeasible(hpset)

# Three of us are killers #code 4
hpset = logicalparse("{{i}} == 'PO' {{3}}", logicset=hpset); showfeasible(hpset)
# This last specification does not reduce the feasible set

# First, however slyly the poison tries to hide
# You will always find some on nettle wine’s left side;
# So if b is poison then a is nettle wine
# The i-1! indicates that outofbounds matches are considered "false"
hpset = logicalparse("{{i}} == 'NW' ==> {{i-1!}} == 'PO'", logicset=hpset); showfeasible(hpset)
# 60 possible outcomes

# Second, different are those who stand at either end,
hpset = logicalparse("a != g", logicset=hpset)

# But if you would move onwards, neither is your friend;
hpset = logicalparse("a,g != 'MA'", logicset=hpset); showfeasible(hpset)
# 30 possible outcomes

# Third, as you see clearly, all are different size,
# Neither dwarf nor giant holds death in their insides;
# c dwarf, f are giant
hpset = logicalparse("c,f != 'PO'", logicset=hpset); showfeasible(hpset)

# Fourth, the second left and the second on the right
# Are twins once you taste them, though different at first sight.
hpset = logicalparse("b == f", logicset=hpset)

search("='MA'", hpset) # c is a match, for the third bottle
search("='MB'", hpset) # g is a match, for the last bottle


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiBatchV1beta1JobTemplateSpec <: SwaggerModel
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiBatchV1JobSpec } # spec name: spec

    function IoK8sApiBatchV1beta1JobTemplateSpec(;metadata=nothing, spec=nothing)
        o = new()
        validate_property(IoK8sApiBatchV1beta1JobTemplateSpec, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiBatchV1beta1JobTemplateSpec, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        o
    end
end # type IoK8sApiBatchV1beta1JobTemplateSpec

const _property_map_IoK8sApiBatchV1beta1JobTemplateSpec = Dict{Symbol,Symbol}(Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"))
const _property_types_IoK8sApiBatchV1beta1JobTemplateSpec = Dict{Symbol,String}(Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiBatchV1JobSpec")
Base.propertynames(::Type{ IoK8sApiBatchV1beta1JobTemplateSpec }) = collect(keys(_property_map_IoK8sApiBatchV1beta1JobTemplateSpec))
Swagger.property_type(::Type{ IoK8sApiBatchV1beta1JobTemplateSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiBatchV1beta1JobTemplateSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiBatchV1beta1JobTemplateSpec }, property_name::Symbol) =  _property_map_IoK8sApiBatchV1beta1JobTemplateSpec[property_name]

function check_required(o::IoK8sApiBatchV1beta1JobTemplateSpec)
    true
end

function validate_property(::Type{ IoK8sApiBatchV1beta1JobTemplateSpec }, name::Symbol, val)
end


text:
function main()
    d = setup_data()

end

function train(d::Dataset)
    opt = SGD()
    for epoch = 1:nepochs
        println("epoch: $(epoch)")
        opt.rate = 0.0075 / epoch
        loss = fit(train_x, train_y, model, crossentropy, opt)
        println("loss: $(loss)")

        test_z = map(x -> predict(model,x), test_x)
        acc = accuracy(test_y, test_z)
        println("test acc.: $(acc)")
        println("")
    end
end

include("intdict.jl")
include("io.jl")
include("model.jl")

max_stars_repo_name: hshindo/Merlin-Examples

max_issues_repo_name: hshindo/Merlin-Examples

max_forks_repo_name: hshindo/Merlin-Examples