text:
#==
Main server, DEPS
==#
module Jockey
using Crayons
using Genie, Logging, LoggingExtras
using Sockets
using Base.Threads: @spawn
using TOML
include("ServerController/Session.jl")
include("ServerController/Settings.jl")
"""
Developer API
## main()
Loads a virtual environment, starts Genie, and boots the server. Not meant to be
called directly, instead it should be
```
main()
```
"""
function main()
  Core.eval(Main, :(const UserApp = $(@__MODULE__)))

  Genie.genie(; context = @__MODULE__)

  Core.eval(Main, :(const Genie = UserApp.Genie))
  Core.eval(Main, :(using Genie))
end
#==
Polling
==#
"""
Developer API
## JockeyPoll(s::string)
Other calls:
- ?(JockeyPoll(s::String, status::Symbol))
JockeyPoll is Jockey's internal communication to the REPL via logging and
terminal output. This is of course vital to working server-side, logs are
important. Providing a string will provide appropriate Jockey return.
```
JockeyPoll("Have a nice day")
[Jockey] Have a nice day!
```
"""
function JockeyPoll(s::String)
  print(Crayon(foreground = :black, bold = true), "[")
  print(Crayon(foreground = :cyan, bold = true), "Jockey")
  print(Crayon(foreground = :black, bold = true), "]  ")
  print(Crayon(foreground = :cyan, bold = false), string(s, "\n"))
end
"""
Developer API
## JockeyPoll(s::string, status::Symbol)
Other calls:
- ?(JockeyPoll(s::String))
JockeyPoll is Jockey's internal communication to the REPL via logging and
terminal output. This is of course vital to working server-side, logs are
important. Providing a string will provide appropriate Jockey return.
The status symbol can be any of the following
- :critical
- :warn
- :request
```
JockeyPoll("Have a nice day")
[Jockey] Have a nice day!
In the future, these will also be **logged.**
```
"""
function JockeyPoll(s::String, status::Symbol)
  # warn, critical,  Request
  status_colours = Dict(:warn => [:red, :magenta],
   :critical => [:black, :red],
    :request => [:cyan, :green])
    status = status_colours[status]
  print(Crayon(foreground = :black, bold = true), "[")
  print(Crayon(foreground = status[1], bold = true), "Jockey")
  print(Crayon(foreground = :black, bold = true), "]  ")
  print(Crayon(foreground = status[2], bold = false), string(s, "\n"))
end
#==
Main run function (runjy)
==#
"""
Usage Functions
### runjy()
The runjy() server begins a Genie server hosting a new Jockey notebook server.
Ultimately, this is the foundation where one will configure server settings and
decide what sort of system they want to run jockey with. The core server settings
are controllable through the key-word arguments to this function.
##### Key-word arguments
- host = 127.0.0.1
IP address one wishes to host this connection on.
- port = 8000
Port which the server should run on.
- mode = :user
Available mode options are :user, :admin, :headless. Admin mode will grant access
 to a user account system, and Jockey authentication systems. Headless mode will
 begin a worker process, which can be further managed using the management
   interface.
- pwd = 123
Password for use at administrator login screen. **You can ignore this if not in
  admin mode!**
- admin = "admin"
"""
function runjy(;port = 8000, host = "127.0.0.1",
  mode = :user,
  pwd = "123",
  admin = "admin")
  HI = HostInfo(port, host, admin, pwd)
  global settings = JSettings(mode, HI)
  setting_dict = Dict("port" => string(port), "host" => host,
  "admin" => admin, "mode" => string(mode))
  open("config/env/settings.toml", "w") do io
           TOML.print(io, setting_dict)
         end

  JockeyPoll(string("New Jockey Server starting on port ", port))
  run(`./bin/server`)
end
#==
Session ASCII
==#
print(Crayon(foreground = :blue, bold = true), """
        .###
       ##########
    ###*    #######            #######
            #####################  ###
          ######################%  ####
          ,#####################
         #       ###/    ##   ###
         #(      #    ###       #*
        (#      ##
""")
print(Crayon(foreground = :cyan, bold = true), """
    &&&                        &&&
             *&&        #&&.   &&&          &&&
    &&&    &&%  &&&  &&&   ,   &&&*&&,   &&&   &&   &&   &&
    &&&   &&&   &&&  &&&       &&&&&&,   &&&&&&&&    && &&
    &&&     &&&&&&    &&&&&&&  &&&  %&&   &&&&&&&     &&&
&&&&&                                              &&&&

""")
println(Crayon(foreground = :light_gray, bold = false), "v. 0.0.2\n\n")
JockeyPoll("Use runjy() to start a new Jockey server, or use ?(runjy) for more information.")
export runjy
end


text:
# # Visual

# The function `visual` can be used to give data-independent visual information about the plot
# (plotting function or attributes).

# The available plotting functions are documented
# [here](http://makie.juliaplots.org/dev/examples/plotting_functions/). Refer to
# plotting functions using upper CamelCase for `visual`'s first argument (e.g.
# `visual(Scatter), visual(BarPlot)`). See the documentation of each plotting
# function to discover the available attributes. These attributes can be passed
# as additional keyword arguments to `visual`, or as part of the `mapping`
# you define.

# In fact `visual` can be used for any plotting function that is defined using
# the `@recipe` macro from Makie. This means that if you have a custom recipe
# defined in another package you can use it from AlgebraOfGraphics just like any
# of the plotting functions defined in Makie.

## Examples

using AlgebraOfGraphics, CairoMakie
set_aog_theme!()

df = (x=randn(1000), y=randn(1000))
plt = data(df) * mapping(:x, :y) * AlgebraOfGraphics.density(npoints=50)
draw(plt * visual(Heatmap)) # plot as heatmap (the default)

#

draw(plt * visual(colormap=:viridis)) # set a different colormap

#

draw(plt * visual(Contour)) # plot as contour

#

draw(plt * visual(Contour, linewidth=2)) # plot as contour with thicker lines

ext: jl


text:
export getHarmonicQ_m, getMorseQ_m, fitHarmonicParams, fitMorseParams, getHarmonicCapture, getMorseCapture

function getHarmonicQ_m(ħω,E0)
    """
    This function gets the value of ΔQ corresponding to an activationless
    Marcus regime (Q_m) for two harmonic potential energy surfaces

    Parameters
    -----------
    ħω : energy of vibration of the initial state, in eV
    E0 : vertical shift between PES, in eV

    Returns
    -------
    Q_m : The value of ΔQ corresponding to an activationless Marcus regime, in amu^0.5 Å^-1
    """
    amu = 931.4940954E6   # eV / c^2
    ħc = 0.19732697E-6    # eV m

    a = amu/2 * (ħω/ħc/1E10)^2
    return (E0/a)^0.5
end

function getMorseQ_m(a, b, E0)
    """
    This function gets the value of Q corresponding to an activationless
    Marcus regime (Q_m) for two Morse potentials

    Parameters
    -----------
    a : first parameter of the Morse potential of the initial state
    b : second parameter of the Morse potential of the final state
    ΔE : vertical shift between PES, in eV

    Returns
    -------
    Q_m : The value of ΔQ corresponding to an activationless Marcus regime, in amu^0.5 Å^-1
    """

    Q_m = (1/b)*log(1-(E0/a)^0.5)
    return Q_m
end

function fitHarmonicParams(ħω_i, ħω_f, ΔQ, ΔE)
    """
    This function fits the parameters to the harmonic potential energy surfaces and
    solves the potentials using a 1D S.E. solver

    Parameters
    ----------
    ħω_i: energy of vibration of the initial state
    ħω_f: energy of vibration of the final state

    ΔQ : horizontal shift between the PES', in amu^0.5 Å^-1
    ΔE : vertical shift between the PES', in eV

    Returns
    -------
    poti, potf: potential of initial and final states
    """
        # set nev
        nev_excited = 180
        nf_min = (1/(ħω_f))*((nev_excited+0.5)ħω_i+2ΔE)-0.5
        nev_relaxed = floor(nf_min+1)

        Q = range(-20-ΔQ, stop=20+ΔQ, length=5000) # ensure high enough to get good harmonic fit

        # potential 1
        potf = potential(); potf.name = "Relaxed state"
        potf.Q0 = ΔQ; potf.E0 = 0
        potf.nev = nev_relaxed # must overlap with initial excited state
        potf.func = x -> harmonic(x, ħω_f; E₀ = potf.E0, Q₀ = potf.Q0)
        potf.params["ħω"]= ħω_f
	    potf.Q = Q
        potf.E = potf.func.(Q)
        solve_pot!(potf)

        # potential 2
        poti = potential(); poti.name = "Excited state"
        poti.Q0 = 0; poti.E0 = ΔE;
        poti.nev = nev_excited # must be converged for partition function
        poti.func = x -> harmonic(x, ħω_i; E₀ = poti.E0, Q₀ = poti.Q0)
        poti.Q = Q
        poti.E = poti.func.(Q)
        poti.params["ħω"]= ħω_i
	    solve_pot!(poti)
        println("parameters fit to potentials")
        return poti, potf
end

function fitMorseParams(a_i, a_f, b_i, b_f, Q0, E0)
    """
    This function fits the parameters to the Morse potential energy surfaces and
    solves the potentials using a 1D S.E. solver

    Parameters
    ----------
    a_i: first parameter of the Morse potential of the initial state
    a_f: first parameter of the Morse potential of the final state
    b_i: second parameter of the Morse potential of the initial state
    b_f: second parameter of the Morse potential of the final state
    Q0 : horizontal shift between the PES', in amu^0.5 Å^-1
    E0 : vertical shift between the PES', in eV

    Returns
    -------
    poti, potf: potential of initial and final states
    """
    # set number of eigenvalues (nev) to solve for in potential
    nev_excited = 180 # must be converged for partition function
    nev_relaxed = 800 # ensure max(ε_f) - max(ε_i) > E0

    Q = range(-20-Q0, stop=20+Q0, length=5000) # ensure high enough to get good harmonic fit

    # potential 1
    poti = potential(); poti.name = "Excited state"
    poti.Q0 = 0; poti.E0 = E0;
    poti.nev = nev_excited
    poti.func = x -> morse(x, [a_i, b_i]; E₀ = poti.E0, Q₀ = poti.Q0)
    poti.Q = Q
    poti.E = poti.func.(Q)
    poti.params["a"]= a_i
    poti.params["b"]= b_i
    solve_pot!(poti)

    # potential 2
    potf = potential(); potf.name = "Relaxed state"
    potf.Q0 = Q0; potf.E0 = 0
    potf.nev = nev_relaxed # must overlap with initial excited state
    potf.func = x -> morse(x, [a_f, b_f]; E₀ = potf.E0, Q₀ = potf.Q0)
    potf.params["a"]= a_f
    potf.params["b"]= b_f
    potf.Q = Q
    potf.E = potf.func.(Q)
    solve_pot!(potf)

    println("parameters fit to potentials")
    return poti, potf
end

function getHarmonicCapture(poti, potf)
    """
    Get the capture coefficient in a charge trapping process from an initial to
    final state described by harmonic potential energy surfaces.

    Parameters
    ----------
    poti: initial potential
    potf: final potential

    Returns
    -------
    capt_coeff: capture coefficient in units of cm^3 /s
    barrier_height: classical barrier to a capture process in eV
    """
    Volume = 1e-21 # cm^3
    σ = 0.8*potf.params["ħω"]
    g = 1 # degeneracy
    temperature = [300]

    ħω= potf.params["ħω"]
    E0 = poti.E0
    println("E0: $E0")
    Q_m = getHarmonicQ_m(ħω, E0)
    Q0 = potf.Q0
    W = 0.068/(Q0-Q_m) # e-ph coupling

    cc = conf_coord(poti, potf) # i, f
    println("conf coordinate fitted")
    cc.W = W
    cc.g = g
    calc_overlap!(cc; σ = σ) # assigns overlap matrix and δ matrix to cc
    println("overlap calculated")
    calc_capt_coeff!(cc, Volume, temperature)
    println("capture coefficient calculated")

    cross = try
                find_crossing(potf,poti)[2]
            catch
                50.0 # if it can't converge, set barrier very high
            end

    if cross==50.0
        barrier_height = 50.0
    else
        barrier_height = cross-poti.E0
    end
    capt_coeff = cc.capt_coeff[1]

    return capt_coeff, barrier_height
end

function getMorseCapture(poti, potf)
    """
    Get the capture coefficient in a charge trapping process from an initial to
    final state described by harmonic potential energy surfaces.

    Parameters
    ----------
    poti: initial potential
    potf: final potential

    Returns
    -------
    capt_coeff: capture coefficient in units of cm^3 /s
    barrier_height: classical barrier to a capture process in eV
    """
    Volume = 1e-21 # cm^3
    σ = 0.01
    g = 1 # degeneracy
    temperature = [300]

    Q0 = potf.Q0
    E0 = poti.E0
    a = potf.params["a"]
    b = potf.params["b"]
    Q_m = getMorseQ_m(a, b, E0)

    W = 0.068/(Q0-Q_m) # e-ph coupling

    cc = conf_coord(poti, potf) # i, f
    println("conf coordinate fitted")
    cc.W = W
    cc.g = g
    calc_overlap!(cc; σ = σ) # assigns overlap matrix and δ matrix to cc
    println("overlap calculated")
    calc_capt_coeff!(cc, Volume, temperature)
    println("capture coefficient calculated")

    cross = try
                find_crossing(potf,poti)[2]
            catch
                50.0 # if it can't converge, set barrier very high
            end

    if cross==50.0
        barrier_height = 50.0
    else
        barrier_height = cross-poti.E0
    end
    capt_coeff = cc.capt_coeff[1]

    return capt_coeff, barrier_height
end


text:
@testset "Subpopulation plots" begin

subpop1 = Subpopulation( "Subpop 1" )
subpop2 = Subpopulation( "Subpop 2" )
subpop3 = Subpopulation( "Subpop 3" )

setSubpopulationSourceNode!.( [ subpop1, subpop2, subpop3 ], [ "", "Master",
    "Branch B" ] )

addSubpopulationCondition!( subpop1, MPcondition( "tenure", >=, 60 ) )
addSubpopulationCondition!( subpop2, MPcondition( "started as", ==,
    "A junior" ) )
addSubpopulationCondition!( subpop3, MPcondition( "isCareer", ==, "yes" ) )

subpopPlot( mpSim, 12, subpop1, subpop2, subpop3 )
subpopAgePlot( mpSim, 12, 12, :age, subpop1, subpop2, subpop3 )
subpopAgePlot( mpSim, 12, 12, :tenure, subpop1, subpop2, subpop3 )
subpopAgePlot( mpSim, 12, 12, :timeInNode, subpop1, subpop2, subpop3 )
#=
subpopPlot( mpSim, 12, subpop1, subpop2, subpop3, showPlot=false,
    filename="plot/subpop plot test" )
subpopAgePlot( mpSim, 12, 12, :age, subpop1, subpop2, subpop3, showPlot=false,
    filename="plot/subpop age plot test 1" )
subpopAgePlot( mpSim, 12, 12, :tenure, subpop1, subpop2, subpop3,
    showPlot=false, filename="plot/subpop age plot test 2" )
subpopAgePlot( mpSim, 12, 12, :timeInNode, subpop1, subpop2, subpop3,
    showPlot=false, filename="plot/subpop age plot test 3" )
=#
end  # @testset "Subpopulation plots"


text:
# planet_parameters = names(CLIMA.Parameters.Planet, all=true)
planet_parameters = ["molmass_dryair","R_d","kappa_d","cp_d","cv_d","ρ_cloud_liq","ρ_cloud_ice","molmass_water","molmass_ratio","R_v","cp_v","cp_l","cp_i","cv_v","cv_l","cv_i","T_freeze","T_min","T_max","T_icenuc","T_triple","T_0","LH_v0","LH_s0","LH_f0","e_int_v0","e_int_i0","press_triple","ρ_ocean","cp_ocean","planet_radius","day","Omega","grav","year_anom","orbit_semimaj","TSI","MSLP"]

function explicit_import(filename, _usingmodule)

  parameters_found = []
  contents = open(f->read(f, String), filename)
  for pp in planet_parameters
    if occursin(pp, contents)
      push!(parameters_found, pp)
    end
  end

  i = findfirst("$(_usingmodule)\n", contents)
  if i ≠ nothing && !isempty(parameters_found)
    contents = replace(contents, "$(_usingmodule)\n" => "$(_usingmodule): "*join(parameters_found, ", ")*"\n")
  end
  if isempty(parameters_found)
    contents = replace(contents, "$(_usingmodule)\n" => "")
  end
  open(filename,"w") do io
    print(io, contents)
  end
end

code_dir = joinpath(@__DIR__, "..")
all_files = [joinpath(root, f) for (root, dirs, files) in Base.Filesystem.walkdir(code_dir) for f in files]
_usingmodule = "using CLIMA.PlanetParameters"

filter!(x->endswith(x, ".jl"), all_files)
filter!(x->!endswith(x, "runtests.jl"), all_files)
filter!(x->!occursin("$(basename(@__FILE__))", x), all_files)
filter!(x->occursin("$(_usingmodule)", open(f->read(f, String), x)), all_files)
filter!(x->!occursin("$(_usingmodule):", open(f->read(f, String), x)), all_files)

@show code_dir
for f in all_files
  file = joinpath(last(split(f, "..")))
  @show file
end

for driver in all_files
  explicit_import(driver, _usingmodule)
end


text:
# Testing discrete univariate distributions

module TestTruncate

using Distributions
using ForwardDiff: Dual
import JSON
using Test

function verify_and_test_drive(jsonfile, selected, n_tsamples::Int,lower::Int,upper::Int)
    R = JSON.parsefile(jsonfile)
    for dct in R
        ex = dct["expr"]
        dsym = Symbol(dct["dtype"])
        if dsym in [:Skellam, :NormalInverseGaussian]
            continue
        end

        dname = string(dsym)

        dsymt = Symbol("Truncated($(dct["dtype"]),$lower,$upper")
        dnamet = string(dsym)

        # test whether it is included in the selected list
        # or all are selected (when selected is empty)
        if !isempty(selected) && !(dname in selected)
            continue
        end

        # perform testing
        dtype = eval(dsym)
        dtypet = Truncated
        d0 = eval(Meta.parse(ex))
        if minimum(d0) > lower || maximum(d0) < upper
            continue
        end

        println("    testing Truncated($(ex),$lower,$upper)")
        d = Truncated(eval(Meta.parse(ex)),lower,upper)
        if dtype != TruncatedNormal
            @assert isa(dtype, Type) && dtype <: UnivariateDistribution
            @test isa(d, dtypet)
        end

        # verification and testing
        verify_and_test(d, dct, n_tsamples)
    end
end

_parse_x(d::DiscreteUnivariateDistribution, x) = round(Int, x)
_parse_x(d::ContinuousUnivariateDistribution, x) = Float64(x)

_json_value(x::Number) = x

_json_value(x::AbstractString) =
    x == "inf" ? Inf :
    x == "-inf" ? -Inf :
    x == "nan" ? NaN :
    error("Invalid numerical value: $x")

function verify_and_test(d::UnivariateDistribution, dct::Dict, n_tsamples::Int)
    # verify stats
    @test minimum(d) ≈ max(_json_value(dct["minimum"]),d.lower)
    @test maximum(d) ≈ min(_json_value(dct["maximum"]),d.upper)
    @test extrema(d) == (minimum(d), maximum(d))

    # verify logpdf and cdf at certain points
    pts = dct["points"]
    for pt in pts
        x = _parse_x(d, pt["x"])
        lp = d.lower <= x <= d.upper ? Float64(pt["logpdf"]) - d.logtp : -Inf
        cf = x <= d.lower ? 0.0 : x >= d.upper ? 1.0 : (Float64(pt["cdf"]) - d.lcdf)/d.tp
        if !isa(d, Distributions.Truncated{Distributions.StudentizedRange{Float64},Distributions.Continuous})
            @test isapprox(logpdf(d, x), lp, atol=sqrt(eps()))
        end
        @test isapprox(cdf(d, x)   , cf, atol=sqrt(eps()))
        # NOTE: some distributions use pdf() in StatsFuns.jl which have no generic support yet
        if !(typeof(d) in [Distributions.Truncated{Distributions.NoncentralChisq{Float64},Distributions.Continuous},
                           Distributions.Truncated{Distributions.NoncentralF{Float64},Distributions.Continuous},
                           Distributions.Truncated{Distributions.NoncentralT{Float64},Distributions.Continuous},
                           Distributions.Truncated{Distributions.StudentizedRange{Float64},Distributions.Continuous}])
            @test isapprox(logpdf(d, Dual(float(x))), lp, atol=sqrt(eps()))
        end
        # NOTE: this test is disabled as StatsFuns.jl doesn't have generic support for cdf()
        # @test isapprox(cdf(d, Dual(x))   , cf, atol=sqrt(eps()))
    end

    # test if truncated quantile function can be evaluated consistently for different types at certain points
    @test isapprox(quantile(d, 0), quantile(d, Float32(0)))
    @test isapprox(quantile(d, 1), quantile(d, Float32(1.0)))
    @test isapprox(quantile(d, Float64(0.3)), quantile(d, Float32(0.3)))
    @test isapprox(quantile(d, Float64(0.7)), quantile(d, Float32(0.7)))

    try
        m = mgf(d,0.0)
        @test m == 1.0
        m = mgf(d,Dual(0.0))
        @test m == 1.0
    catch e
        isa(e, MethodError) || throw(e)
    end
    try
        c = cf(d,0.0)
        @test c == 1.0
        c = cf(d,Dual(0.0))
        @test c == 1.0
        # test some extra values: should all be well-defined
        for t in (0.1,-0.1,1.0,-1.0)
            @test !isnan(cf(d,t))
            @test !isnan(cf(d,Dual(t)))
        end
    catch e
        isa(e, MethodError) || throw(e)
    end

    # generic testing
    if isa(d, Cosine)
        n_tsamples = floor(Int, n_tsamples / 10)
    end
end

## main

for c in ["discrete",
          "continuous"]

    title = string(uppercase(c[1]), c[2:end])
    println("    [$title]")
    println("    ------------")
    jsonfile = joinpath(dirname(@__FILE__), "ref", "$(c)_test.ref.json")
    verify_and_test_drive(jsonfile, ARGS, 10^6,3,5)
    println()
end

end


text:
using DelimitedFiles
using Test
include("../Analysis/AlphaShapes/AlphaShapes.jl")
X = rand(1000,2).*100;
writedlm("./delaunay-test.txt",X,',')

function Density(X)
    tess,inds = AlphaShapes.GetDelaunayTriangulation(X,true);
	d = zeros(size(X,1))
	for i in 1:size(X,1)
    	d[i] = AlphaShapes.WeightedDTFELocalDensity(i,tess,inds)
	end
	return d
end

D_julia = Density(X)
run(`g++ -std=c++11 -O3 delaunator-test.cpp -o test`)
run(`./test delaunay-test.txt`)
D_cpp = readdlm("out.txt",'\n')[:,1]
D_julia.-D_cpp#

@show @test sum(abs.(D_julia .- D_cpp))/length(D_julia) < 1e-6

run(`rm test out.txt delaunay-test.txt`)

max_stars_repo_name: harveydevereux/CUDA-Whirligigs

max_issues_repo_name: harveydevereux/CUDA-Whirligigs

max_forks_repo_name: harveydevereux/CUDA-Whirligigs


text:
#####
##### Nothing happens when your boundary condition is nothing
#####

fill_west_halo!(c, ::Nothing, arch, barrier, grid, args...) = NoneEvent()
fill_east_halo!(c, ::Nothing, arch, barrier, grid, args...) = NoneEvent()

fill_north_halo!(c, ::Nothing, arch, barrier, grid, args...) = NoneEvent()
fill_south_halo!(c, ::Nothing, arch, barrier, grid, args...) = NoneEvent()

fill_top_halo!(c, ::Nothing, arch, barrier, grid, args...)    = NoneEvent()
fill_bottom_halo!(c, ::Nothing, arch, barrier, grid, args...) = NoneEvent()


text:
using ArgParse

function args_pq(s)
    @add_arg_table s begin
        "--M"
            help = "M"
            arg_type = Int
            default = -1
        "--Ks"
            help = "Ks"
            arg_type = Int
            default = -1
        "--D"
            help = "D"
            arg_type = Int
            default = -1
        "--pq_codebook", "-b"
            help = "path"
            arg_type = String
            default = "-1"
        "--pq_codes","-c"
            help = "path"
            arg_type = String
            default = "-1"

        "--metric","-m"
            help = "D"
            arg_type = String
            default = "-1"
        end
        return s
end

function args_data(s)
    @add_arg_table s begin

        "--queries","-q"
        help = "path"
        arg_type = String
        default = "-1"

        "--datan"
        help = "data name"
        arg_type = String
        default = "-1"

        "--tr100"
        help = "path"
        arg_type = String
        default = "-1"
    end
    return s
end

function args_vq(s)
    @add_arg_table s begin
        "--vq_codebook"
        help = "path"
        arg_type = String
        default = "-1"

        "--vq_codes"
        help = "path"
        arg_type = String
        default = "-1"
    end
end

function args_recall(s)
    @add_arg_table s begin
        "--num_leaves_to_search"
        help = "num_leaves_to_search"
        arg_type = Int
        default = -1

        "--topk"
        help = "num_leaves_to_search"
        arg_type = Int
        default = -1

    end
end


text:
@testset "Hilbert symbols" begin
  v = [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, 1, 1,
       -1, -1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1, 1, 1, 1, -1, 1, -1, 1, -1, 1, 1,
       -1, -1, -1, -1, 1, 1, -1, 1, 1, -1, -1, 1, 1, 1, -1, 1, -1, -1, 1, 1 ]

  for i in 1:8
    for j in 1:8
      @test hilbert_symbol(i, j, 2) == v[(i - 1) * 8 + j]
      @test hilbert_symbol(fmpz(i), fmpz(j), 2) == v[(i - 1) * 8 + j]
      @test hilbert_symbol(fmpz(i), fmpz(j), fmpz(2)) == v[(i - 1) * 8 + j]
      @test hilbert_symbol(fmpq(i), fmpq(j), 2) == v[(i - 1) * 8 + j]
      @test hilbert_symbol(fmpq(i), fmpq(j), fmpz(2)) == v[(i - 1) * 8 + j]
    end
  end

  for p in PrimesSet(3, 100)
    for a in 1:100
      for b in 1:00
        h = hilbert_symbol(a, b, p)
        a = fmpz(a)
        b = fmpz(b)
        r = (-1)^(valuation(a, p) * valuation(b, p)) * a^(valuation(b, p)) * b^(valuation(a, p))
        @test h == jacobi_symbol(r, fmpz(p))
      end
    end
  end

  Qx, x = PolynomialRing(FlintQQ, "x")
  K, b = NumberField(x^3-3*x-1, "a")
  OK = maximal_order(K)
  for P in prime_ideals_up_to(OK, 200)
    @test hilbert_symbol(b, -3, P) == 1
  end

  # Test where Magma div(x, y) differs from julia div(x, y) (internally)
  K, a = CyclotomicRealSubfield(8, "a") # x^2 - 2
  z = 9278908160780559301//4*a+6561375391013480455//2
  w = K(-2)
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  @test hilbert_symbol(z, w, p) == 1

  k, a = rationals_as_number_field()
  kt, t = k["t"]
  LS, LSb = number_field(t^2 - 2)
  LNS, LNSb = number_field([t^2 - 2])
  KS, KSb = number_field(x^2 - 2)
  KNS, KNSb = number_field([x^2 - 2])
  p = prime_decomposition(maximal_order(k), 2)[1][1]
  pKS = prime_decomposition(maximal_order(KS), 2)[1][1]
  pKNS = prime_decomposition(maximal_order(KNS), 2)[1][1]
  pLS = prime_decomposition(maximal_order(LS), p)[1][1]
  pLNS = prime_decomposition(maximal_order(LNS), p)[1][1]

  q = prime_decomposition(maximal_order(k), 3)[1][1]
  qKS = prime_decomposition(maximal_order(KS), 3)[1][1]
  qKNS = prime_decomposition(maximal_order(KNS), 3)[1][1]
  qLS = prime_decomposition(maximal_order(LS), q)[1][1]
  qLNS = prime_decomposition(maximal_order(LNS), q)[1][1]

  for i in 1:10
    z = rand(KS, -2:2)
    v = quadratic_defect(z, pKS)
    #@test v == quadratic_defect(KNS(coordinates(z)), pKNS)
    #@test v == quadratic_defect(LS(k.(coordinates(z))), pLS)
    #@test v == quadratic_defect(LNS(k.(coordinates(z))), pLNS)
    v = quadratic_defect(z, qKS)
    #@test v == quadratic_defect(KNS(coordinates(z)), qKNS)
    #@test v == quadratic_defect(LS(k.(coordinates(z))), qLS)
    #@test v == quadratic_defect(LNS(k.(coordinates(z))), qLNS)
  end

  for i in 1:8
    for j in 1:8
      v = hilbert_symbol(i, j, 2)
      @test v == hilbert_symbol(k(i), k(j), p)
      #@test v == hilbert_symbol(KNS(i), KNS(i), pKNS)
      #@test v == hilbert_symbol(LS(i), LS(i), pLS)
      #@test v == hilbert_symbol(LNS(i), LNS(i), pLNS)
    end
  end
end


text:
using Slack

"""
    sendtoslack(data::String, user_endpoint::String)
Takes in a string of data and the endpoint url and returns the HTTP response.
The response will be "ok" if the request is sent successfully.
"""
function sendtoslack(data::String, user_endpoint::String)
    HTTP.setuseragent!("HTTP.jl")

    params = Dict("text"=>"$(data)")
    base_url = "https://hooks.slack.com"

    endpoint = user_endpoint
    url = base_url * endpoint

    r = HTTP.post(url,
                 ["Content-Type" => "application/json"],
                 JSON.json(params))

    return (String(r.body))
end

"""
    sendattachmenttoslack(data::String, user_endpoint::String)
Takes in a string of data and the endpoint url and returns the HTTP response.
The response will be "ok" if the request is sent successfully.
"""
function sendattachmenttoslack(data, user_endpoint::String)
    HTTP.setuseragent!("HTTP.jl")

    params = data
    base_url = "https://hooks.slack.com"

    endpoint = user_endpoint
    url = base_url * endpoint

    r = HTTP.post(url,
                 ["Content-Type" => "application/json"],
                 JSON.json(params))

    return (String(r.body))
end


text:
module EDutils

using Iterators
using LinearMaps
using TensorOperations

"""
Applies operator `op` to state `s` at sites `n` to `n+r-1`,
where `r` is the number of sites on which `op` acts nontrivially (the range).

Both `s` and `op` must be tensors with one index per site:

A state on `Ns` sites has `ndim(s) == Ns` .

An operator taking `Ns_op` sites to `Ns_op_out` sites has `ndim(op) == Ns_op_out + Ns_op`.
The `Ns_op` indices acting on the state are assumed to be the last `Ns_op` indices.

Periodic boundary conditions are used to act on sites past the end of the chain.
"""
function apply_op_at!(s, res, op, n::Int; α=1.0, β=0.0, Ns_op_out=0, s_sz=size(s))
    Ns = length(s_sz)
    op_sz = size(op)

    if Ns_op_out == 0
        Ns_op = ndims(op) ÷ 2
        Ns_op_out = Ns_op
    else
        Ns_op = ndims(op) - Ns_op_out
    end

    if n == 1
        #Make this more obviously just matrix multiplication where possible. Why does tensorcontract!() fail to do this by itself?
        s = reshape(s, (prod(s_sz[1:Ns_op]), prod(s_sz[Ns_op+1:end])))
        op = reshape(op, (prod(op_sz[1:Ns_op_out]), size(s, 1)))
        res_r = reshape(res, (size(op,1), size(s, 2)))
        tensorcontract!(α, op, [:o1, :o2], 'N', s, [:o2, :sR], 'N', β, res_r, [:o1, :sR])
    else
        #Must avoid crossing the max. tuple-length limit for type-inference.
        #This is the general case. Some of these dimensions will be of size 1.

        s = reshape(s, (prod(s_sz[1:n + Ns_op - 1 - Ns]), #wrapped operator index
                        prod(s_sz[max(n + Ns_op - Ns, 1):n - 1]), #rest left
                        prod(s_sz[n:min(n + Ns_op - 1, Ns)]), #operator index
                        prod(s_sz[n + Ns_op:Ns]))) #rest right

        op_wrap_ind = min(Ns_op_out, Ns - n + 1)
        op = reshape(op, (prod(op_sz[1:op_wrap_ind]), #out
                          prod(op_sz[op_wrap_ind+1:Ns_op_out]), #out, wrapped
                          size(s, 3), size(s, 1))) #in

        res_r = reshape(res, (size(op, 2), size(s, 2), size(op, 1), size(s, 4))) #this should not alter the storage

        opinds = [:op_o1, :op_o2, :op1, :op2]
        inds = [:op2, :sL, :op1, :sR]
        inds_out = [:op_o2, :sL, :op_o1, :sR]

        tensorcontract!(α, op, opinds, 'N', s, inds, 'N', β, res_r, inds_out)
    end

    res
end

"""
Translates a state by `n` sites (site 1 becomes site N, site 2 becomes site 1, and so on).
"""
function translate_pbc!(s, st, n::Int, a::Number, b::Number, s_sz)
    Ns = length(s_sz)
    st_sz = size(st)

    #ensure type-inference works by reducing to a minimal number of indicies.
    if n >= 0
        s = reshape(s, (prod(s_sz[1:n]), prod(s_sz[n+1:end])))
    else
        s = reshape(s, (prod(s_sz[1:Ns+n]), prod(s_sz[Ns+n+1:end])))
    end

    st = reshape(st, (size(s,2), size(s,1)))

    if a == 1.0 && b == 0.0
        transpose!(st, s)
    else
        tensoradd!(a, s, (:s1, :s2), b, st, (:s2, :s1)) #at least in Julia 0.5, tensoradd!() is slower, particularly for n < 0
    end

    st = reshape(st, st_sz)

    st, ntuple(j->s_sz[mod1(j+n,Ns)], Ns)
end

function translate_pbc!(s, st, n::Int; a=1.0, b=0.0, s_sz=size(s), ret_size::Bool=false)
    st, s_sz = translate_pbc!(s, st, n, a, b, s_sz)
    ret_size ? (st, s_sz) : st
end

function translate_pbc(s, n::Int; s_sz=size(s), ret_size::Bool=false)
    res = zeros(s)
    translate_pbc!(s, res, n, s_sz=s_sz, ret_size=ret_size)
end

function apply_global_prod_op!(s::AbstractArray, ops::Vector{AbstractArray}; tmp = zeros(s), s_sz=size(s))
    for n in 1:length(s_sz)
        apply_op_at!(s, tmp, ops[n], 1, s_sz=s_sz)
        s, s_sz = translate_pbc!(tmp, s, 1, 1.0, 0.0, s_sz)
    end
    s
end

function apply_global_prod_op!(s::AbstractArray, lop::AbstractArray; tmp = zeros(s), s_sz=size(s))
    apply_global_prod_op!(s, AbstractArray[lop for j in 1:length(s_sz)], tmp=tmp, s_sz=s_sz)
end

#NOTE: Modifies s!
function apply_local_op!{T}(s::Vector{T}, s_sz::Tuple{Vararg{Int}}, H::Vector{Array}, s_tmp::Vector{T}, res_tmp1::Vector{T}, res_tmp2::Vector{T})
    fill!(res_tmp1, 0.0)
    for n in 1:length(s_sz)
        apply_op_at!(s, res_tmp1, H[n], 1, β=1.0, s_sz=s_sz)

        translate_pbc!(res_tmp1, res_tmp2, 1, 1.0, 0.0, s_sz) #translating by +1 is usually faster than -1
        tmp = res_tmp1
        res_tmp1 = res_tmp2
        res_tmp2 = tmp

        if n < length(s_sz)
            translate_pbc!(s, s_tmp, 1, 1.0, 0.0, s_sz)
            tmp = s
            s = s_tmp
            s_tmp = tmp
        end
    end
    res_tmp1
end

function apply_local_op{T}(s::Vector{T}, s_sz::Tuple{Vararg{Int}}, H::Vector{Array})
    s_vlen = prod(s_sz)
    res1 = zeros(T, s_vlen)
    res2 = zeros(T, s_vlen)
    s1 = copy(s)
    s2 = zeros(T, s_vlen)
    apply_local_op!(s1, s_sz, H, s2, res1, res2)
end

function apply_local_op(s::Vector, s_sz::Tuple{Vararg{Int}}, Hloc::Array)
    Ns = length(s_sz)
    H = Array[Hloc for n in 1:Ns]
    apply_local_op(s, s_sz, H)
end

#The part of the Hamiltonian spectrum obtained is determined by the spectrum and the choice of "which" parameter.
#To get the low-energy states, ensure the spectrum of H is negative.
#Resolving momenta requires H to have real spectrum and be positive or negative definite.
function eigs_localHam_circle(H::Vector{Array}; momenta::Bool=false, mom_dir=:L, symm_func=nothing, nev::Int=41, ncv::Int=max(20, 2*nev+1), which=:LM)
    Ns = length(H)
    Nh = ndims(H[1]) ÷ 2
    D = size(H[1], 1)

    s_sz = ntuple(x->D, Ns)
    s_vlen = prod(s_sz)
    @show Ns, Nh, s_sz

    #We work with vectors because Julia freaks when the number of array dimensions gets too high (type-inference fails).
    res1 = zeros(Complex128, s_vlen)
    res2 = zeros(Complex128, s_vlen)
    s1 = zeros(Complex128, s_vlen)
    s2 = zeros(Complex128, s_vlen)
    s3 = zeros(Complex128, s_vlen)
    function f!(dest::AbstractVector{Complex128}, s_in::AbstractVector{Complex128})
        copy!(s1, s_in)
        #s_in may be a subarray... make sure we have an Array, otherwise TensorOps complains.
        #Also creating a reshapedsubarray can be slooooow.

        if momenta
            s2 = translate_pbc!(s1, s2, mom_dir == :L ? 1 : -1, s_sz=s_sz)
            tmp = s1
            s1 = s2
            s2 = tmp
        end

        if symm_func != nothing
            s1 = symm_func(s1, s2, s3, s_sz)
        end

        Hsr = apply_local_op!(s1, s_sz, H, s2, res1, res2)
        copy!(dest, Hsr)
        dest
    end

    #FIXME: Calling convention changed since 0.2.0
    Hop = LinearMap(f!, s_vlen, Complex128, ishermitian=false, ismutating=true)

    ev, eV, nconv, niter, nmult, resid = eigs(Hop, nev=nev, ritzvec=true, which=which, ncv=ncv)

    if momenta
        Hev = -abs(ev)
        Tev = ev ./ Hev
        return Hev, Tev, eV
    else
        return ev, eV
    end
end

function eigs_localHam_circle(Hloc::Array, Ns::Int; momenta::Bool=false, mom_dir=:L, symm_func=nothing, nev::Int=41, ncv::Int=max(20, 2*nev+1))

    Nh = ndims(Hloc) ÷ 2
    D = size(Hloc, 1)
    Hlocr = reshape(Hloc, (D^Nh, D^Nh))
    em = maximum(real(eigvals(Hlocr)))

    Hlocr = Hlocr - I * em #offset so that we can find eigenvalues with largest magnitude
    Hloc = reshape(Hlocr, size(Hloc))

    hs = Array[Hloc for n in 1:Ns]

    res = eigs_localHam_circle(hs, momenta=momenta, mom_dir=mom_dir, symm_func=symm_func, nev=nev)

    res = (res[1] + em * Ns, res[2:end]...)

    res
end

function get_Pmoms(ps::Vector{Float64})
    (s, tmp1, tmp2, s_sz)->begin
        copy!(tmp2, s)
        fill!(s, 0.0)

        N = length(s_sz)
        for n in 1:N
            tmp1 = translate_pbc!(tmp2, tmp1, 1, s_sz=s_sz)
            tmp1, tmp2 = (tmp2, tmp1)

            LinAlg.BLAS.axpy!(sum(cis(-n * p * 2pi / N) for p in ps), tmp2, s)
        end

        scale!(s, 1/N)
        s
    end
end

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

const pauliX = [0.0 1.0;
                1.0 0.0]
const pauliY = 1.0im*[0.0 -1.0;
                      1.0  0.0]
const pauliZ = [1.0  0.0;
                0.0 -1.0]

function ham_ising_comp(λ::Number=1.0, hz::Number=1.0, hx::Number=0.0)
    -(hz*pauliZ + hx*pauliX), -λ * kron(pauliX, pauliX)
end

function ham_ising(λ::Number=1.0, hz::Number=1.0, hx::Number=0.0)
    h1, h2 = ham_ising_comp(λ, hz, hx)

    h2 + 1/2 * kron(h1, eye(2)) + 1/2 * kron(eye(2), h1)
end

function weylops(p::Int)
    om = cis(2π / p)
    U = diagm(Complex128[om^j for j in 0:p-1])
    V = diagm(ones(p - 1), 1)
    V[end, 1] = 1
    U, V, om
end

function weylops_basis(p::Int)
    U, V, om = weylops(p)
    ops = Matrix[(U^j for j in 1:p-1)..., (V^j for j in 1:p-1)..., (U^j * V^k for j in 1:p-1 for k in 1:p-1)...]
    lbls = [("U^$j" for j in 1:p-1)..., ("V^$j" for j in 1:p-1)..., ("U^$(j)V^$(k)" for j in 1:p-1 for k in 1:p-1)...]
    ops, lbls
end

function ham_vpotts_comp(p::Int, J::Float64=1.0, h::Float64=1.0, hU::Number=0.0)
    U, V, om = weylops(p)

    UUh = kron(U, U')

    -0.5 * h * (V + V') - hU * 1/3 * (U + U' + I), -0.5 * J * (UUh + UUh')
end

function ham_vpotts(p::Int, J::Float64=1.0, h::Float64=1.0, hU::Number=0.0)
    h1, h2 = ham_vpotts_comp(p, J, h, hU)

    h2 + 1/2 * kron(h1, eye(p)) + 1/2 * kron(eye(p), h1)
end

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

function ed(N::Int, h::Array, nev::Int; symm_func=nothing, ncv::Int=max(20, 2*nev+1))
    @time en_pbc, mom_pbc, eV_pbc = eigs_localHam_circle(h, N, momenta=true, symm_func=symm_func, nev=nev, ncv=ncv)
    return en_pbc, mom_pbc, eV_pbc
end

function ed!(ev_data, Ns::Vector{Int}, h::Array, nev::Int; symm_func=nothing)
    for N in Ns
        println("Starting $N")
        if !haskey(ev_data, N)
            ev_data[N] = ed(N, h, nev, symm_func=symm_func)
        end

        eV = ev_data[N][3]
        eVisU = vecnorm(eV' * eV - I)
        println("Is eV isometry? $(eVisU)")
    end
    ev_data
end

end

max_stars_repo_name: amilsted/spinchain-virasoro-demo

max_issues_repo_name: amilsted/spinchain-virasoro-demo

max_forks_repo_name: amilsted/spinchain-virasoro-demo


text:
const Shift = 1
const ReduceL = 2
const ReduceR = 3

type State
    tokens::Vector{Token}
    step::Int
    top::Int
    right::Int
    left::State
    lc::State
    rc::State

    State(tokens) = new(tokens, 1, 1, 2)
end

Base.done(s::State) = isdefined(s,:left) && s.right == length(s.tokens)+1

function Base.next(s::State, acts::Vector{Int})
    scores = s.scorefun(s, acts)
    states = State[]
    for i = 1:length(acts)
        push!(states, expand(s, acts[i], scores[i]))
    end
    states
end

function Base.next(s::State, act::Int, score::Float64)
    top, lcs, rcs, left, right = begin
        if act == shift
            s.right, Int[], Int[], s, s.right+1
        elseif act == reducel
            _lcs = length(s.lcs) == 0 ? [s.left.top] : [s.left.top,s.lcs[1]]
            s.top, _lcs, s.rcs, s.left.left, s.right
        elseif act == reducer
            _rcs = length(s.left.rcs) == 0 ? [s.top] : [s.top,s.left.rcs[1]]
            s.left.top, s.left.lcs, _rcs, s.left.left, s.right
        else
            throw("Invalid action: $(act)")
        end
    end
    score = s.score + score
    State(top, lcs, rcs, left, right, s.step+1, s, act, s.tokens, s.scorefun, score, nothing)
end

function expand_gold(s::State)
    isfinal(s) && return State[]
    if s.left == nothing
        act = shift
    else
        tokens = s.tokens
        s0, s1 = tokens[s.top], tokens[s.left.top]
        if s1.headid == s.top
            act = reducel
        elseif s0.headid == s.left.top
            reducable = all(i -> tokens[i].headid != s.top, s.right:length(tokens))
            act = reducable ? reducer : shift
        elseif s.right <= length(tokens)
            act = shift
        else
            throw("Invalid")
        end
    end
    expand(s, [act])
end

function expand_pred(s::State)
    isfinal(s) && return State[]
    acts = Int[]
    s.left == nothing || push!(acts, reducel, reducer)
    s.right <= length(s.tokens) && push!(acts, shift)
    expand(s, acts)
end

function toheads(s::State)
    heads = fill(0, length(s.tokens))
    while s != nothing
        length(s.lcs) > 0 && (heads[s.lcs[1]] = s.top)
        length(s.rcs) > 0 && (heads[s.rcs[1]] = s.top)
        s = s.prev
    end
    heads
end

max_stars_repo_name: hshindo/JukaiNLP


text:
struct IBP{T<:Real}
    α::T
    # TODO: IBP{T} -> IBP below and Line 11 could be removed
    function IBP{T}(α) where {T<:Real}
        @assert α > 0 "α is not positve"
        return new{T}(α)
    end
end
IBP(α::T) where {T<:Real} = IBP{T}(α)

"""
    rand(ibp::IBP, n::Int, k_max::Int)

Sample from IBP using the stick-breaking construction.
The returned sample is of size `n`-by-`k_max`,
which is not a Julia convension by a math convension.
"""
function rand(ibp::IBP, n::Int, k_max::Int)
    ν = rand(Distributions.Beta(ibp.α, 1), k_max)
    p = break_stick_ibp(ν)
    Z = hcat(rand.(Distributions.Bernoulli.(p), n)...)
    return Z
end

"""
    rand(ibp::IBP, n::Int)

Sample from IBP using the buffet metaphor.
The returned sample is of size `n`-by-`Kmax`,
which is not a Julia convension by a math convension.
"""
function rand(ibp::IBP, n::Int)
    k_init = rand(Distributions.Poisson(ibp.α))
    Z = zeros(Int, n, k_init)
    Z[1,:] .= 1
    for i = 2:n
        for k = 1:size(Z, 2)
            mk = sum(Z[:,k])
            Z[i,k] = rand(Distributions.Bernoulli(mk / i))
        end
        k_new = rand(Distributions.Poisson(ibp.α / i))
        if k_new > 0
            Z = hcat(Z, zeros(Int, n, k_new))
            Z[i,end-k_new+1:end] .= 1
        end
    end
    return Z
end


text:
################################################################################
#
#          NfOrd/Elem.jl : Elements of orders of number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015, 2016: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

export deepcopy, call, parent, elem_in_nf, elem_in_basis, discriminant, hash,
       ==, zero, one, iszero, isone, show, -, +, *, divexact, ^, mod, powermod,
       representation_mat, trace, norm, rand, rand!, add!, mul!, minkowski_map,
       conjugates_arb, conjugates_arb_log, t2

################################################################################
#
#  Deepcopy
#
################################################################################

doc"""
***
    deepcopy(x::NfOrdElem) -> NfOrdElem

> Returns a copy of $x$.
"""
function deepcopy(x::NfOrdElem)
  z = parent(x)()
  z.elem_in_nf = deepcopy(x.elem_in_nf)
  if x.has_coord
    z.has_coord = true
    z.elem_in_basis = fmpz[ deepcopy(y) for y in x.elem_in_basis ]
  end
  return z
end

################################################################################
#
#  Parent object overloading
#
################################################################################

#if VERSION < v"0.5.0-"
#  doc"""
#  #      call(O::NfOrd, a::nf_elem, check::Bool = true) -> NfOrdElem
#
#  > Given an element $a$ of the ambient number field of $\mathcal O$, this
#  > function coerces the element into $\mathcal O$. It will be checked that $a$
#  > is contained in $\mathcal O$ if and only if `check` is `true`.
#  """
#  Base.call{T <: NfOrd}(O::T, a::nf_elem, check::Bool = true) = nothing
#
#  for T in (:NfOrdGen, :NfMaxOrd)
#    @eval begin
#      function Base.call(O::$T, a::nf_elem, check::Bool = true)
#        if check
#          (x,y) = _check_elem_in_order(a,O)
#          !x && error("Number field element not in the order")
#          return NfOrdElem{$T}(O, deepcopy(a), fmpz[ deepcopy(x) for x in y])
#        else
#          return NfOrdElem{$T}(O, deepcopy(a))
#        end
#      end
#    end
#  end
#
#  doc"""
#  #      call(O::NfOrd, a::Union{fmpz, Integer}) -> NfOrdElem
#
#  > Given an element $a$ of type `fmpz` or `Integer`, this
#  > function coerces the element into $\mathcal O$. It will be checked that $a$
#  > is contained in $\mathcal O$ if and only if `check` is `true`.
#  """
#  Base.call{T <: NfOrd}(O::T, a::Union{fmpz, Integer}) = nothing
#
#  for T in (:NfOrdGen, :NfMaxOrd)
#    @eval begin
#      function Base.call(O::$T, a::Union{fmpz, Integer})
#        return NfOrdElem{$T}(O, nf(O)(a))
#      end
#    end
#  end
#
#  doc"""
#  #      call(O::NfOrd, arr::Array{fmpz, 1})
#
#  > Returns the element of $\mathcal O$ with coefficient vector `arr`.
#  """
#  Base.call{T <: NfOrd}(O::T, a::Array{fmpz, 1}) = nothing
#
#  for T in (:NfOrdGen, :NfMaxOrd)
#    @eval begin
#      function Base.call(O::$T, arr::Array{fmpz, 1})
#        return NfOrdElem{$T}(O, fmpz[ deepcopy(x) for x in arr ])
#      end
#    end
#  end
#
#
#  doc"""
#  #      call{T <: Integer}(O::NfOrd, arr::Array{T, 1})
#
#  > Returns the element of $\mathcal O$ with coefficient vector `arr`.
#  """
#  Base.call{T <: NfOrd, S <: Integer}(O::T, arr::Array{S, 1}) = nothing
#
#  for T in (:NfOrdGen, :NfMaxOrd)
#    @eval begin
#      function Base.call{S <: Integer}(O::$T, arr::Array{S, 1})
#        return NfOrdElem{$T}(O, deepcopy(arr))
#      end
#    end
#  end
#
#  doc"""
#  #      call(O::NfOrd, a::nf_elem, arr::Array{fmpz, 1}) -> NfOrdElem
#
#  > This function constructs the element of $\mathcal O$ with coefficient vector
#  > `arr`. It is assumed that the corresponding element of the ambient number
#  > field is $a$.
#  """
#  Base.call{T <: NfOrd}(O::T, a::nf_elem, arr::Array{fmpz, 1}) = nothing
#
#  for T in (:NfOrdGen, :NfMaxOrd)
#    @eval begin
#      function Base.call(O::$T, a::nf_elem, arr::Array{fmpz, 1})
#        return NfOrdElem{$T}(O, deepcopy(a), fmpz[ deepcopy(x) for x in arr])
#      end
#    end
#  end
#
#  doc"""
#  #      call(O::NfOrd) -> NfOrdElem
#
#  > This function constructs a new element of $\mathcal O$ which is set to $0$.
#  """
#  Base.call{T <: NfOrd}(O::T)
#
#  for T in (:NfOrdGen, :NfMaxOrd)
#    @eval begin
#      Base.call(O::$T) = NfOrdElem{$T}(O)
#    end
#  end
#end
#if VERSION > v"0.5.0-"
  doc"""
  ***
        call(O::NfOrd, a::nf_elem, check::Bool = true) -> NfOrdElem

  > Given an element $a$ of the ambient number field of $\mathcal O$, this
  > function coerces the element into $\mathcal O$. It will be checked that $a$
  > is contained in $\mathcal O$ if and only if `check` is `true`.
  """
  Base.call(O::NfOrdGen, a::nf_elem, check::Bool = true) = begin
    if check
      (x,y) = _check_elem_in_order(a,O)
      !x && error("Number field element not in the order")
      return NfOrdElem{NfOrdGen}(O, deepcopy(a), fmpz[ deepcopy(x) for x in y])
    else
      return NfOrdElem{NfOrdGen}(O, deepcopy(a))
    end
  end

  doc"""
  ***
        call(O::NfOrd, a::nf_elem, check::Bool = true) -> NfOrdElem

  > Given an element $a$ of the ambient number field of $\mathcal O$, this
  > function coerces the element into $\mathcal O$. It will be checked that $a$
  > is contained in $\mathcal O$ if and only if `check` is `true`.
  """
  Base.call(O::NfMaxOrd, a::nf_elem, check::Bool = true) = begin
    if check
      (x,y) = _check_elem_in_order(a,O)
      !x && error("Number field element not in the order")
      return NfOrdElem{NfMaxOrd}(O, deepcopy(a), fmpz[ deepcopy(x) for x in y])
    else
      return NfOrdElem{NfMaxOrd}(O, deepcopy(a))
    end
  end

  doc"""
  ***
        call(O::NfOrd, a::Union{fmpz, Integer}) -> NfOrdElem

  > Given an element $a$ of type `fmpz` or `Integer`, this
  > function coerces the element into $\mathcal O$. It will be checked that $a$
  > is contained in $\mathcal O$ if and only if `check` is `true`.
  """
  Base.call(O::NfOrdGen, a::Union{fmpz, Integer}) = begin
    return NfOrdElem{NfOrdGen}(O, nf(O)(a))
  end

  doc"""
  ***
        call(O::NfOrd, a::Union{fmpz, Integer}) -> NfOrdElem

  > Given an element $a$ of type `fmpz` or `Integer`, this
  > function coerces the element into $\mathcal O$. It will be checked that $a$
  > is contained in $\mathcal O$ if and only if `check` is `true`.
  """
  Base.call(O::NfMaxOrd, a::Union{fmpz, Integer}) = begin
    return NfOrdElem{NfMaxOrd}(O, nf(O)(a))
  end

  doc"""
  ***
        call(O::NfOrd, arr::Array{fmpz, 1})

  > Returns the element of $\mathcal O$ with coefficient vector `arr`.
  """
  Base.call(O::NfOrdGen, arr::Array{fmpz, 1}) = begin
    return NfOrdElem{NfOrdGen}(O, fmpz[ deepcopy(x) for x in arr ])
  end

  doc"""
  ***
        call(O::NfOrd, arr::Array{fmpz, 1})

  > Returns the element of $\mathcal O$ with coefficient vector `arr`.
  """
  Base.call(O::NfMaxOrd, arr::Array{fmpz, 1}) = begin
    return NfOrdElem{NfMaxOrd}(O, fmpz[ deepcopy(x) for x in arr ])
  end

  doc"""
  ***
        call{T <: Integer}(O::NfOrd, arr::Array{T, 1})

  > Returns the element of $\mathcal O$ with coefficient vector `arr`.
  """
  Base.call{S <: Integer}(O::NfOrdGen, arr::Array{S, 1}) = begin
    return NfOrdElem{NfOrdGen}(O, deepcopy(arr))
  end

  doc"""
  ***
        call{T <: Integer}(O::NfOrd, arr::Array{T, 1})

  > Returns the element of $\mathcal O$ with coefficient vector `arr`.
  """
  Base.call{S <: Integer}(O::NfMaxOrd, arr::Array{S, 1}) = begin
    return NfOrdElem{NfMaxOrd}(O, deepcopy(arr))
  end

  doc"""
  ***
        call(O::NfOrd, a::nf_elem, arr::Array{fmpz, 1}) -> NfOrdElem

  > This function constructs the element of $\mathcal O$ with coefficient vector
  > `arr`. It is assumed that the corresponding element of the ambient number
  > field is $a$.
  """
  Base.call(O::NfOrdGen, a::nf_elem, arr::Array{fmpz, 1}) = begin
    return NfOrdElem{NfOrdGen}(O, deepcopy(a), fmpz[ deepcopy(x) for x in arr])
  end

  doc"""
  ***
        call(O::NfOrd, a::nf_elem, arr::Array{fmpz, 1}) -> NfOrdElem

  > This function constructs the element of $\mathcal O$ with coefficient vector
  > `arr`. It is assumed that the corresponding element of the ambient number
  > field is $a$.
  """
  Base.call(O::NfMaxOrd, a::nf_elem, arr::Array{fmpz, 1}) = begin
    return NfOrdElem{NfMaxOrd}(O, deepcopy(a), fmpz[ deepcopy(x) for x in arr])
  end

  doc"""
  ***
        call(O::NfOrd) -> NfOrdElem

  > This function constructs a new element of $\mathcal O$ which is set to $0$.
  """
  Base.call(O::NfOrdGen) = NfOrdElem{NfOrdGen}(O)

  doc"""
  ***
        call(O::NfOrd) -> NfOrdElem

  > This function constructs a new element of $\mathcal O$ which is set to $0$.
  """
  Base.call(O::NfMaxOrd) = NfOrdElem{NfMaxOrd}(O)
#end

################################################################################
#
#  Parent
#
################################################################################

doc"""
***
    parent(a::NfOrdElem) -> NfOrd

> Returns the order of which $a$ is an element.
"""
parent(a::NfOrdElem) = a.parent

################################################################################
#
#  Element in number field
#
################################################################################

doc"""
***
    elem_in_nf(a::NfOrdElem) -> nf_elem

> Returns the element $a$ considered as an element of the ambient number field.
"""
function elem_in_nf(a::NfOrdElem)
  if isdefined(a, :elem_in_nf)
    return deepcopy(a.elem_in_nf)
  end
  error("Not a valid order element")
end

################################################################################
#
#  Coordinates
#
################################################################################

doc"""
***
    elem_in_basis(a::NfOrdElem) -> Array{fmpz, 1}

> Returns the coefficient vector of $a$.
"""
function elem_in_basis(a::NfOrdElem)
  if a.has_coord
    return fmpz[ deepcopy(x) for x in a.elem_in_basis]
  else
    (x,y) = _check_elem_in_order(a.elem_in_nf, parent(a))
    !x && error("Not a valid order element")
    a.elem_in_basis = fmpz[ deepcopy(x) for x in y]
    a.has_coord = true
    @hassert :NfOrd 2 a == dot(a.elem_in_basis, basis(parent(a)))
    return fmpz[ deepcopy(x) for x in a.elem_in_basis]
  end
end

################################################################################
#
#  Discriminant
#
################################################################################

doc"""
***
    discriminant(B::Array{NfOrdElem, 1}) -> fmpz

> Returns the discriminant of the family $B$.
"""
function discriminant{T}(B::Array{NfOrdElem{T}, 1})
  length(B) == 0 && error("Number of elements must be non-zero")
  length(B) != degree(parent(B[1])) &&
        error("Number of elements must be $(degree(parent(B[1])))")
  O = parent(B[1])
  A = MatrixSpace(FlintZZ, degree(O), degree(O))()
  for i in 1:degree(O)
    for j in 1:degree(O)
      A[i,j] = FlintZZ(trace(B[i]*B[j]))
    end
  end
  return det(A)
end

################################################################################
#
#  Hashing
#
################################################################################

hash(x::NfOrdElem) = hash(x.elem_in_nf)

################################################################################
#
#  Equality
#
################################################################################

doc"""
***
    ==(x::NfOrdElem, y::NfOrdElem) -> Bool

> Returns whether $x$ and $y$ are equal.
"""
 ==(x::NfOrdElem, y::NfOrdElem) = parent(x) == parent(y) &&
                                            x.elem_in_nf == y.elem_in_nf

################################################################################
#
#  Special elements
#
################################################################################

doc"""
***
    zero(O::NfOrd) -> NfOrdElem

> Returns the zero element of $\mathcal O$.
"""
zero(O::NfOrd) = O(fmpz(0))

doc"""
***
    one(O::NfOrd) -> NfOrdElem

> Returns the zero element of $\mathcal O$.
"""
one(O::NfOrd) = O(fmpz(1))

doc"""
***
    zero(a::NfOrdElem) -> NfOrdElem

> Returns the zero element of the parent of $a$.
"""
zero(a::NfOrdElem) = parent(a)(0)

doc"""
***
    one(O::NfOrd) -> NfOrdElem

> Returns the one element of the parent of $a$.
"""
one(a::NfOrdElem) = parent(a)(1)

################################################################################
#
#  isone/iszero
#
################################################################################

doc"""
***
    isone(a::NfOrd) -> Bool

> Tests if $a$ is one.
"""
isone(a::NfOrdElem) = isone(a.elem_in_nf)

doc"""
***
    iszero(a::NfOrd) -> Bool

> Tests if $a$ is one.
"""
iszero(a::NfOrdElem) = iszero(a.elem_in_nf)

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, a::NfOrdElem)
  print(io, a.elem_in_nf)
end

################################################################################
#
#  Unary operations
#
################################################################################

doc"""
***
    -(x::NfOrdElem) -> NfOrdElem

> Returns the additive inverse of $x$.
"""
function -(x::NfOrdElem)
  z = parent(x)()
  z.elem_in_nf = - x.elem_in_nf
  if x.has_coord
    z.elem_in_basis = map(y -> -y, x.elem_in_basis)
    z.has_coord = true
  end
  return z
end

###############################################################################
#
#  Binary operations
#
###############################################################################

doc"""
***
    *(x::NfOrdElem, y::NfOrdElem) -> NfOrdElem

> Returns $x \cdot y$.
"""
function *(x::NfOrdElem, y::NfOrdElem)
  parent(x) != parent(y) && error("Wrong parents")
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf*y.elem_in_nf
  return z
end

doc"""
***
    +(x::NfOrdElem, y::NfOrdElem) -> NfOrdElem

> Returns $x + y$.
"""
function +(x::NfOrdElem, y::NfOrdElem)
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf + y.elem_in_nf
  if x.has_coord && y.has_coord
    z.elem_in_basis =
      fmpz[ x.elem_in_basis[i] + y.elem_in_basis[i] for i = 1:degree(parent(x))]
    z.has_coord = true
  end
  return z
end

doc"""
***
    -(x::NfOrdElem, y::NfOrdElem) -> NfOrdElem

> Returns $x - y$.
"""
function -(x::NfOrdElem, y::NfOrdElem)
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf - y.elem_in_nf
  if x.has_coord && y.has_coord
    z.elem_in_basis =
      fmpz[ x.elem_in_basis[i] - y.elem_in_basis[i] for i = 1:degree(parent(x))]
    z.has_coord = true
  end
  return z
end

################################################################################
#
#  Ad hoc operations
#
################################################################################

function *(x::NfOrdElem, y::Integer)
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf * y
  if x.has_coord
    z.elem_in_basis = map(z -> y*z, x.elem_in_basis)
    z.has_coord = true
  end
  return z
end

*(x::Integer, y::NfOrdElem) = y * x

function *(x::NfOrdElem, y::fmpz)
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf * y
  if x.has_coord
    z.elem_in_basis = map(z -> y*z, x.elem_in_basis)
    z.has_coord = true
  end
  return z
end

*(x::fmpz, y::NfOrdElem) = y * x

Base.dot(x::fmpz, y::NfOrdElem) = y*x

function +(x::NfOrdElem, y::Integer)
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf + y
  return z
end

+(x::Integer, y::NfOrdElem) = y + x

function +(x::NfOrdElem, y::fmpz)
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf + y
  return z
end

+(x::fmpz, y::NfOrdElem) = y + x

function -(x::NfOrdElem, y::fmpz)
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf - y
  return z
end

function -(x::NfOrdElem, y::Int)
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf - y
  return z
end

-(x::fmpz, y::NfOrdElem) = -(y - x)

-(x::Integer, y::NfOrdElem) = -(y - x)

################################################################################
#
#  Adhoc exact division
#
################################################################################

function divexact(x::NfOrdElem, y::Int)
  z = parent(x)()
  z.elem_in_nf = divexact(x.elem_in_nf, y)
  return z
end

function divexact(x::NfOrdElem, y::fmpz)
  z = parent(x)()
  z.elem_in_nf = divexact(x.elem_in_nf, y)
  return z
end

################################################################################
#
#  Exponentiation
#
################################################################################

doc"""
***
    ^(x::NfOrdElem, y::Union{fmpz, Int})

> Returns $x^y$.
"""
function ^(x::NfOrdElem, y::Union{fmpz, Int})
  z = parent(x)()
  z.elem_in_nf = x.elem_in_nf^y
  return z
end

################################################################################
#
#  Modular reduction
#
################################################################################

doc"""
***
    mod(a::NfOrdElem, m::Union{fmpz, Int}) -> NfOrdElem

> Reduces the coefficient vector of $a$ modulo $m$ and returns the corresponding
> element. The coefficient vector of the result will have entries $x$ with
> $0 \leq x \leq m$.
"""
function mod(a::NfOrdElem, m::Union{fmpz, Int})
  ar = elem_in_basis(a)
  for i in 1:degree(parent(a))
    ar[i] = mod(ar[i],m)
  end
  return parent(a)(ar)
end

################################################################################
#
#  Modular exponentiation
#
################################################################################

doc"""
***
    powermod(a::NfOrdElem, i::fmpz, m::Union{fmpz, Int}) -> NfOrdElem

> Returns the element $a^i$ modulo $m$.
"""

function powermod(a::NfOrdElem, i::fmpz, p::fmpz)
  if i == 0 then
    return one(parent(a))
  end
  if i == 1
    b = mod(a,p)
    return b
  end
  if mod(i,2) == 0
    j = div(i, 2)
    b = powermod(a, j, p)
    b = b^2
    b = mod(b,p)
    return b
  end
  b = mod(a*powermod(a,i-1,p),p)
  return b
end

doc"""
***
    powermod(a::NfOrdElem, i::Integer, m::Integer) -> NfOrdElem

> Returns the element $a^i$ modulo $m$.
"""
powermod(a::NfOrdElem, i::Integer, m::Integer)  = powermod(a, fmpz(i), fmpz(m))

doc"""
***
    powermod(a::NfOrdElem, i::fmpz, m::Integer) -> NfOrdElem

> Returns the element $a^i$ modulo $m$.
"""
powermod(a::NfOrdElem, i::fmpz, m::Integer)  = powermod(a, i, fmpz(m))

doc"""
***
    powermod(a::NfOrdElem, i::Integer, m::fmpz) -> NfOrdElem

> Returns the element $a^i$ modulo $m$.
"""
powermod(a::NfOrdElem, i::Integer, m::fmpz)  = powermod(a, fmpz(i), m)

################################################################################
#
#  Representation matrices
#
################################################################################

doc"""
***
    representation_mat(a::NfOrdElem) -> fmpz_mat

> Returns the representation matrix of the element $a$.
"""
function representation_mat(a::NfOrdElem)
  O = parent(a)
  A = representation_mat(a, nf(parent(a)))
  A = basis_mat(O)*A*basis_mat_inv(O)
  !(A.den == 1) && error("Element not in order")
  return A.num
end

doc"""
***
    representation_mat(a::NfOrdElem, K::AnticNumberField) -> FakeFmpqMat

> Returns the representation matrix of the element $a$ considered as an element
> of the ambient number field $K$. It is assumed that $K$ is the ambient number
> field of the order of $a$.
"""
function representation_mat(a::NfOrdElem, K::AnticNumberField)
  nf(parent(a)) != K && error("Element not in this field")
  d = den(a.elem_in_nf)
  b = d*a.elem_in_nf
  A = representation_mat(b)
  z = FakeFmpqMat(A, d)
  return z
end

#  nf(parent(a)) != K && error("Element not in this field")
#  return representation_mat(a.elem_in_nf)
#end

################################################################################
#
#  Trace
#
################################################################################

doc"""
***
    trace(a::NfOrdElem) -> fmpz

> Returns the trace of $a$.
"""
function trace(a::NfOrdElem)
  return FlintZZ(trace(elem_in_nf(a)))
end

################################################################################
#
#  Norm
#
################################################################################

doc"""
***
    norm(a::NfOrdElem) -> fmpz

> Returns the norm of $a$.
"""
function norm(a::NfOrdElem)
  return FlintZZ(norm(elem_in_nf(a)))
end

################################################################################
#
#  Random element generation
#
################################################################################

function rand!{T <: Integer}(z::NfOrdElem, O::NfOrd, R::UnitRange{T})
  y = O()
  ar = rand(R, degree(O))
  B = basis(O)
  mul!(z, ar[1], B[1])
  for i in 2:degree(O)
    mul!(y, ar[i], B[i])
    add!(z, z, y)
  end
  return z
end

doc"""
***
    rand(O::NfOrd, R::UnitRange{Integer}) -> NfOrdElem

> Computes a coefficient vector with entries uniformly distributed in `R` and returns
> the corresponding element of the order.
"""
function rand{T <: Integer}(O::NfOrd, R::UnitRange{T})
  z = O()
  rand!(z, O, R)
  return z
end

function rand!(z::NfOrdElem, O::NfOrd, n::Union{Integer, fmpz})
  return rand!(z, O, -n:n)
end

doc"""
***
    rand(O::NfOrd, n::Union{Integer, fmpz}) -> NfOrdElem

> Computes a coefficient vector with entries uniformly distributed in
> $\{-n,\dotsc,-1,0,1,\dotsc,n\}$ and returns the corresponding element of the
> order $\mathcal O$.
"""
function rand(O::NfOrd, n::Integer)
  return rand(O, -n:n)
end

function rand(O::NfOrd, n::fmpz)
  z = O()
  rand!(z, O, BigInt(n))
  return z
end

function rand!(z::NfOrdElem, O::NfOrd, n::fmpz)
  return rand!(z, O, BigInt(n))
end

################################################################################
#
#  Unsafe operations
#
################################################################################

function add!(z::NfOrdElem, x::NfOrdElem, y::NfOrdElem)
  z.elem_in_nf = x.elem_in_nf + y.elem_in_nf
  if x.has_coord && y.has_coord
    for i in 1:degree(parent(x))
      z.elem_in_basis[i] = x.elem_in_basis[i] + y.elem_in_basis[i]
    end
    z.has_coord = true
  else
    z.has_coord = false
  end
  nothing
end

function mul!(z::NfOrdElem, x::NfOrdElem, y::NfOrdElem)
  z.elem_in_nf = x.elem_in_nf * y.elem_in_nf
  z.has_coord = false
  nothing
end

function mul!(z::NfOrdElem, x::fmpz, y::NfOrdElem)
  z.elem_in_nf = x * y.elem_in_nf
  if y.has_coord
    for i in 1:degree(parent(z))
      z.elem_in_basis[i] = x*y.elem_in_basis[i]
    end
    z.has_coord = true
  else
    z.has_coord = false
  end
  nothing
end

mul!(z::NfOrdElem, x::Integer, y::NfOrdElem) =  mul!(z, ZZ(x), y)

mul!(z::NfOrdElem, x::NfOrdElem, y::Integer) = mul!(z, y, x)

function add!(z::NfOrdElem, x::fmpz, y::NfOrdElem)
  z.elem_in_nf = y.elem_in_nf + x
  z.has_coord = false
  nothing
end

add!(z::NfOrdElem, x::NfOrdElem, y::fmpz) = add!(z, y, x)

function add!(z::NfOrdElem, x::Integer, y::NfOrdElem)
  z.elem_in_nf = x + y.elem_in_nf
  z.has_coord = false
  nothing
end

add!(z::NfOrdElem, x::NfOrdElem, y::Integer) = add!(z, y, x)

mul!(z::NfOrdElem, x::NfOrdElem, y::fmpz) = mul!(z, y, x)

################################################################################
#
#  Base cases for dot product of vectors
#
################################################################################

dot(x::NfOrdElem, y::Int) = x*y

dot(x::Int, y::NfOrdElem) = dot(y, x)

dot(x::NfOrdElem, y::fmpz) = x*y

################################################################################
#
#  Conversion
#
################################################################################

Base.call(K::AnticNumberField, x::NfOrdElem) = elem_in_nf(x)

################################################################################
#
#  Minkowski embedding
#
################################################################################

doc"""
***
    minkowski_map(a::NfOrdElem, abs_tol::Int) -> Array{arb, 1}

> Returns the image of $a$ under the Minkowski embedding.
> Every entry of the array returned is of type `arb` with radius less then
> `2^-abs_tol`.
"""
function minkowski_map(a::NfOrdElem, abs_tol::Int = 32)
  # Use a.elem_in_nf instead of elem_in_nf(a) to avoid copying the data.
  # The function minkowski_map does not alter the input!
  return minkowski_map(a.elem_in_nf, abs_tol)
end

################################################################################
#
#  Conjugates
#
################################################################################

doc"""
***
    conjugates_arb(x::NfOrdElem, abs_tol::Int) -> Array{acb, 1}

> Compute the the conjugates of `x` as elements of type `acb`.
> Recall that we order the complex conjugates
> $\sigma_{r+1}(x),...,\sigma_{r+2s}(x)$ such that
> $\sigma_{i}(x) = \overline{\sigma_{i + s}(x)}$ for $r + 2 \leq i \leq r + s$.
>
> Every entry `y` of the array returned satisfies `radius(real(y)) < 2^-abs_tol`,
> `radius(imag(y)) < 2^-abs_tol` respectively.
"""
function conjugates_arb(x::NfOrdElem, abs_tol::Int = 32)
  # Use a.elem_in_nf instead of elem_in_nf(a) to avoid copying the data.
  # The function minkowski_map does not alter the input!
  return conjugates_arb(x.elem_in_nf, abs_tol)
end

doc"""
***
    conjugates_arb_log(x::NfOrdElem, abs_tol::Int) -> Array{arb, 1}

> Returns the elements
> $(\log(\lvert \sigma_1(x) \rvert),\dotsc,\log(\lvert\sigma_r(x) \rvert),
> \dotsc,2\log(\lvert \sigma_{r+1}(x) \rvert),\dotsc,
> 2\log(\lvert \sigma_{r+s}(x)\rvert))$ as elements of type `arb` radius
> less then `2^-abs_tol`.
"""
function conjugates_arb_log(x::NfOrdElem, abs_tol::Int = 32)
  return conjugates_arb_log(x.elem_in_nf, abs_tol)
end

################################################################################
#
#  T2
#
################################################################################

doc"""
***
    t2(x::NfOrdElem, abs_tol::Int = 32) -> arb

> Return the $T_2$-norm of $x$. The radius of the result will be less than
> `2^-abs_tol`.
"""
function t2(x::NfOrdElem, abs_tol::Int = 32)
  return t2(x.elem_in_nf, abs_tol)
end

################################################################################
#
#  Promotion
#
################################################################################

Base.promote_rule{T <: Integer, S <: NfOrd}(::Type{NfOrdElem{S}}, ::Type{T}) =
        NfOrdElem{S}


text:
import Base: getindex, setindex!

function getindex(cluster::MixedCluster, i::Integer)
    cluster.engines[i]
end

function setindex!(cl::MixedCluster, egn::Engine, i::Integer)
    cl.engines[i] = egn
end

function getindex(cluster::UniformCluster, _)
    cluster.engine
end


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")

using Plots

##This code plots the simscale results
function get_data(filename)
    f1 = open(filename)
    tempdata = readlines(f1)
    print(tempdata[1],'\n')
    #use this for "Dataset/FINAL_DATA_.txt" for other files may not
    #tempdata[1] = tempdata[1][4:end]  #cutting of beggining bc some weird chars
    close(f1)
    data = []
    regvec = zeros(length(tempdata),
                   length(split(tempdata[1], ",")))
    #parsing datavector to float array... not regularized
    for i in 1:length(tempdata)
        regvec[i,:] =  parse.(Float64,split(tempdata[i], ","))
    end
    return regvec
end

#getting simscale results
fname = "Dataset/FINAL_DATA_reorganized.txt"
data_tmp = get_data(fname)

#series color is gradient based on altitude
alt=zeros(1,length(data_tmp[:,1]));mach=zeros(1,length(data_tmp[:,1]));
fault=zeros(1,length(data_tmp[:,1]));pressure=zeros(1,length(data_tmp[:,1]));

machline=[];faultline=[];pressureline=[];
specified_mach = 0.3;specified_alt=1525.0
for i in 1:length(data_tmp[:,1])
    alt[1,i] = data_tmp[i,1]
    mach[1,i] = data_tmp[i,2]
    fault[1,i] = data_tmp[i,3]
    pressure[1,i] = data_tmp[i,4]/1000.0
    if mach[1,i] == specified_mach && alt[1,i] == specified_alt
        push!(machline, mach[1,i])
        push!(pressureline, pressure[1,i])
        push!(faultline, fault[1,i])
    end
end

plotly()
scatter(fault,mach,pressure,
        title="SIMSCALE Results",
        xlabel="Fault",
        ylabel="Mach",
        zlabel="Pressure (kPa)",
        c=:alpine,
        colorbar_title="Alt",
        marker_z=alt,
        clims=(1525.0,12000.0),
        markersize=2,
        label=nothing)
plot!(faultline,machline,pressureline,
        label="Constant Mach/Alt",
        linecolor=:red)

#second plot
fault=[];mach=[];alt=[];pressure=[];

machline=[];faultline=[];pressureline=[];
specified_mach = 0.3;specified_alt=1525.0
for i in 1:length(data_tmp[:,1])
    if data_tmp[i,1]==specified_alt
        push!(alt, data_tmp[i,1])
        push!(mach, data_tmp[i,2])
        push!(fault, data_tmp[i,3])
        push!(pressure, data_tmp[i,4]/1000.0)
    end
    if data_tmp[i,2] == specified_mach && data_tmp[i,1] == specified_alt
        push!(machline, data_tmp[i,2])
        push!(pressureline, data_tmp[i,4]/1000.0)
        push!(faultline, data_tmp[i,3])
    end
end
gr()
scatter(fault,mach,pressure,
    title="Altitude Slice (1525m)",
    xlabel="Fault",
    ylabel="Mach",
    zlabel="Pressure (kPa)",
    label=nothing)
plot!(faultline,machline,pressureline,
        label="Constant Mach/Alt",
        linecolor=:red)

max_stars_repo_name: NoahLopezGit/PitotTubeNeuralNetwork

max_issues_repo_name: NoahLopezGit/PitotTubeNeuralNetwork

max_forks_repo_name: NoahLopezGit/PitotTubeNeuralNetwork


text:
# imported from LinSolvers.jl

export LinSolverCreator, BackslashLinSolverCreator;
export FactorizeLinSolverCreator, DefaultLinSolverCreator;
export GMRESLinSolverCreator;

export create_linsolver;

abstract type LinSolverCreator ; end

"""
    struct BackslashLinSolverCreator <: LinSolverCreator

Creator to for the `BackslashLinSolver`, i.e., usage of backslash to make linear solves.
Specify objects of this type if you want the solver to use backslash.

See also: [`BackslashLinSolver`](@ref), [`create_linsolver`](@ref)
"""
struct BackslashLinSolverCreator <: LinSolverCreator
end

"""
    create_linsolver(creator::LinSovlerCreator,nep,λ)

Creates a `LinSolver` instance for the `nep` corresponding
which is evaluated in `λ`. The type of the output is
decided by dispatch and the type of the `LinSolverCreator`.

See also: `LinSolver`, `FactorizeLinSolverCreator`,
`BackslashLinSolverCreator`, `DefaultLinSolverCreator`,
`GMRESLinSolverCreator`.
"""
function create_linsolver(creator::BackslashLinSolverCreator,nep,λ)
    return BackslashLinSolver(nep,λ);
end

"""
    FactorizeLinSolverCreator(;umfpack_refinements,max_factorizations,nep,precomp_values)

`FactorizeLinSolverCreator`-objects can instantiate `FactorizeLinSolver`
objects via the `create_linsolver` function.

The `FactorizeLinSolver` is based on `factorize`-calls.
The time point of the call to `factorize` can be controlled by parameters
to `FactorizeLinSolverCreator`:

* By default, the `factorize` call is carried out by the instantiation of the `FactorizeLinSolver`, i.e., when the NEP-solver calls `create_linsolver`.

* You can also precompute the factorization, at the time point when you instantiate `FactorizeLinSolverCreator`. If you set `precomp_values::Vector{Number}` to a non-empty vector, and set `nep` kwarg, the factorization (of all λ-values in the `precomp_values`) will be computed  when the `FactorizeLinSolverCreator` is instantiated. If the NEP-solver calls a `create_linsolver` with a λ-value from that vector, the factorization will be used (otherwise it will be computed).

Further recycling is possible. If the variable `max_factorizations` is set
to a positive value, the object will store that many factorizations
for possible reuse. Every `lin_solve`-call then computes a factorization,
unless a `lin_solve`-call for that `λ` has been computed earlier.
This procedure can at most store `max_factorization` (which can be set `Inf`).

See also: [`FactorizeLinSolver`](@ref), [`create_linsolver`](@ref)

"""
struct FactorizeLinSolverCreator{T_values,T_factor} <: LinSolverCreator
    umfpack_refinements::Int;
    recycled_factorizations::Dict{T_values,T_factor};
    max_factorizations::Int
    function FactorizeLinSolverCreator(;umfpack_refinements::Int=2,
                                       max_factorizations=0,
                                       nep=nothing,
                                       precomp_values=[]
                                       )

        if (precomp_values isa Number) # Helper: Allow values to be a scalar. Put in a vector
            precomp_values=[precomp_values];
        end

        if (size(precomp_values,1)>0 && nep==nothing)
            error("When you want to precompute factorizations you need to supply the keyword argument `nep`");
        end

        # Compute all the factorizations
        precomp_factorizations=map(s-> factorize(compute_Mder(nep,s)), precomp_values);

        # Put them in a dict
        T_from=eltype(precomp_values);
        local T_to,T_from
        if (size(precomp_values,1)>0)
            T_to=eltype(precomp_factorizations);
        else
            T_to=Any;
            T_from=Number;
        end

        dict=Dict{T_from,T_to}();
        for i=1:size(precomp_values,1)
            dict[precomp_values[i]]=precomp_factorizations[i];
        end

        return new{T_from,T_to}(umfpack_refinements,dict,max_factorizations)

    end
end

# For the moment, Factorize is the default behaviour
"""
    DefaultLinSolverCreator

This is the default linear solver if no other is specified (for most methods).
It is a `FactorizeLinSolverCreator`.

See also: [`LinSolver`](@ref), [`create_linsolver`](@ref),
[`lin_solve`](@ref), [`FactorizeLinSolverCreator`](@ref), [`FactorizeLinSolver`](@ref)

"""
DefaultLinSolverCreator = FactorizeLinSolverCreator

function create_linsolver(creator::FactorizeLinSolverCreator,nep,λ)
    # Let's see if we find it in recycled_factorizations
    if (λ in keys(creator.recycled_factorizations))
        Afact=creator.recycled_factorizations[λ];
        return FactorizeLinSolver(Afact,creator.umfpack_refinements);
    else
        solver=FactorizeLinSolver(nep,λ,creator.umfpack_refinements);
        if  (length(keys(creator.recycled_factorizations))
             < creator.max_factorizations )
            # Let's save the factorization
            creator.recycled_factorizations[λ]=solver.Afact;
        end
        return solver;
    end

end

struct GMRESLinSolverCreator{T} <: LinSolverCreator where {T}
    kwargs::T
end

"""
    GMRESLinSolverCreator(;kwargs...)

This is the creator for the GMRES-method. Instantiate this
object if you want to use GMRES as your linear system solver.
The `kwargs` are stored and used as keyword arguments in the
call to gmres. See list
of keyword in the [IterativeSolvers.jl manual](https://juliamath.github.io/IterativeSolvers.jl/dev/linear_systems/gmres/).

"""
function GMRESLinSolverCreator(;kwargs...)
    return GMRESLinSolverCreator{typeof(kwargs)}(kwargs)
end

function create_linsolver(creator::GMRESLinSolverCreator,nep, λ)
    return GMRESLinSolver{typeof(λ)}(nep, λ, creator.kwargs)
end


text:
export MNISTData

immutable MNISTData{B}
	numClasses::Int
	numTrain::Int
	numTest::Int
	trainImages::AFArray{B}
	testImages::AFArray{B}
	trainLabels::AFArray{B}
	testLabels::AFArray{B}
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
module data

export readReportsFromJson, ReportDataset, readDataset, load_preselected_candidates, readDataset, PreselectedCandidates, load_training_pairs, load_negative_candidates

using DataStructures
using preprocessing
using util
using JSON
using Dates

function readReportsFromJson(filePath, preOpt::PreprocessingOption; stopId::Int64=0)
    l = Dict{Int64, Report}()
    vocab = Vocab(Dict{String, UInt32}(), Set())
    lastCreationts = -1
    reportJsons = JSON.parsefile(filePath)

    for reportJson in reportJsons
        reportId = reportJson["bug_id"]
        dup_id = reportJson["dup_id"]

        if isnothing(dup_id)
            dup_id = reportId
        end

        # println("Report: $(reportId)")

        curCreationts = reportJson["creation_ts"]
        creation_day_ts = calculateDayTimestamp(curCreationts)

        stacks = preprocessStacktrace(reportId, reportJson["stacktrace"], vocab, preOpt)
        l[reportId] = Report(reportId, dup_id, curCreationts, creation_day_ts, stacks)

        if stopId > 0
            if lastCreationts > curCreationts
                error("You cannot remove reports created after a report if the json is not sorted")
            end

            if stopId == reportId
                break
            end
        end

        lastCreationts = curCreationts
    end

    return (l,vocab)
end

struct ReportDataset
    info::String
    bugIds::Vector{Int64}
    duplicateIds::Vector{Int64}
end

function readDataset(file)
    open(file, "r") do io
        info = strip(readline(io))
        bugIds = [parse(Int64, id) for id in split(strip(readline(io)))]
        duplicateIds = [parse(Int64, id) for id in split(strip(readline(io)))]

        return ReportDataset(info, bugIds, duplicateIds)
    end
end

function load_training_pairs(path)
    pairs = MutableLinkedList{Tuple{Int64,Int64,Int64}}()

    open(path, "r") do io
        for line in eachline(io)
            split_line = split(strip(line), ",")
            push!(pairs, (parse(Int64, split_line[1]), parse(Int64, split_line[2]), parse(Int64, split_line[3])))
        end
    end

    return pairs
end

function load_negative_candidates(path)
    negative_candidates = Vector{Vector{Int64}}()

    open(path, "r") do io
        for line in eachline(io)
            split_line = split(strip(line), " ")

            v = Vector{Int64}(undef, length(split_line))

            for i=1:length(split_line)
                v[i] = parse(Int64, split_line[i])
            end

            push!(negative_candidates, v)
        end
    end

    return negative_candidates
end

struct PreselectedCandidates
    query_id::Int64
    candidates::Vector{Int64}
end

function load_preselected_candidates(path)
    preselectedList = MutableLinkedList{PreselectedCandidates}()

    open(path, "r") do io
        for line in eachline(io)
            split_line = split(strip(line), " ")

            v = Vector{Int64}(undef, length(split_line) - 1)

            for i=2:length(split_line)
                v[i-1] = parse(Int64, split_line[i])
            end

            push!(preselectedList, PreselectedCandidates(parse(Int64, split_line[1]), v))
        end
    end

    return preselectedList
end

end

max_stars_repo_name: irving-muller/FaST

max_issues_repo_name: irving-muller/FaST

max_forks_repo_name: irving-muller/FaST


text:
module RocketHelpersUpdatesTest

using Test
using Rocket

import Rocket: GenericUpdatesStatus, UInt8UpdatesStatus, getustorage
import Rocket: cstatus, vstatus, ustatus, cstatus!, vstatus!, ustatus!
import Rocket: fill_cstatus!, fill_vstatus!, fill_ustatus!, all_cstatus, all_vstatus, all_ustatus

@testset "Helpers Updates" begin

    println("Testing: updates")

    @testset begin
        for i in 9:100
            T = typeof(ntuple(_ -> 1, i))
            u = getustorage(T)
            @test typeof(u) <: GenericUpdatesStatus
        end
    end

    @testset begin
        for i in 1:8
            T = typeof(ntuple(_ -> 1, i))
            u = getustorage(T)
            @test typeof(u) <: UInt8UpdatesStatus
        end
    end

    @testset begin
        for n in 1:12
            T = typeof(ntuple(_ -> 1, n))
            u = getustorage(T)

            for i in 1:n
                for _ in 1:2
                    cstatus!(u, i, true)
                    @test cstatus(u, i) === true
                    cstatus!(u, i, false)
                    @test cstatus(u, i) === false
                    vstatus!(u, i, true)
                    @test vstatus(u, i) === true
                    vstatus!(u, i, false)
                    @test vstatus(u, i) === false
                    ustatus!(u, i, true)
                    @test ustatus(u, i) === true
                    ustatus!(u, i, false)
                    @test ustatus(u, i) === false

                    fill_cstatus!(u, true)
                    @test all_cstatus(u) === true
                    for k in 1:n
                        @test cstatus(u, k) === true
                    end
                    cstatus!(u, i, false)
                    @test all_cstatus(u) === false
                    fill_cstatus!(u, false)
                    @test all_cstatus(u) === false
                    for k in 1:n
                        @test cstatus(u, k) === false
                    end

                    fill_vstatus!(u, true)
                    @test all_vstatus(u) === true
                    for k in 1:n
                        @test vstatus(u, k) === true
                    end
                    vstatus!(u, i, false)
                    @test all_vstatus(u) === false
                    fill_vstatus!(u, false)
                    @test all_vstatus(u) === false
                    for k in 1:n
                        @test vstatus(u, k) === false
                    end

                    fill_ustatus!(u, true)
                    @test all_ustatus(u) === true
                    for k in 1:n
                        @test ustatus(u, k) === true
                    end
                    ustatus!(u, i, false)
                    @test all_ustatus(u) === false
                    fill_ustatus!(u, false)
                    @test all_ustatus(u) === false
                    for k in 1:n
                        @test ustatus(u, k) === false
                    end

                end
            end
        end

    end

end

end


text:
module ChangePointMean

using Random
using Statistics
using Plots

export mcppv                # p-value of existence of change point
export mcptime              # index of change point (new segment)
export mcpoint              # Single change point in Mean
export mcplast              # Most recent change point (rightmost)
export mcpall               # All change point
export mcplot               # plotting routine

include("meancpd.jl")

end


text:
using XUnit: @testset, @testcase

@testset "server-metrics.jl" begin
  include("server-metrics.jl")
end
@testset "statsd-exporter.jl" begin
  include("statsd-exporter.jl")
end
@testset "prometheus-exporter.jl" begin
  include("prometheus-exporter.jl")
end


text:
# ---
# title: 50. Pow(x, n)
# id: problem50
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Binary Search
# link: <https://leetcode.com/problems/powx-n/description/>
# hidden: true
# ---
#
# Implement [pow( _x_ , _n_
# )](http://www.cplusplus.com/reference/valarray/pow/), which calculates  _x_
# raised to the power _n_ (i.e. xn).
#
#
#
# **Example 1:**
#
#
#
#     Input: x = 2.00000, n = 10
#     Output: 1024.00000
#
#
# **Example 2:**
#
#
#
#     Input: x = 2.10000, n = 3
#     Output: 9.26100
#
#
# **Example 3:**
#
#
#
#     Input: x = 2.00000, n = -2
#     Output: 0.25000
#     Explanation: 2-2 = 1/22 = 1/4 = 0.25
#
#
#
#
# **Constraints:**
#
#   * `-100.0 < x < 100.0`
#   * `-231 <= n <= 231-1`
#   * `-104 <= xn <= 104`
#
#
## @lc code=start
using LeetCode

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


text:
using Parameters
include("macros.jl")

const Optional{F} = Union{F,Nothing}

struct ObjectGraphNode
    type::Symbol
    children::Vector{ObjectGraphNode}
    parents::Vector{ObjectGraphNode}
    attributes::Set{Symbol}
end

const ObjectGraph = Dict{Symbol, ObjectGraphNode}

struct Object
    type::Symbol
    name::String
    id::Symbol
    attrs::Dict{Symbol, Any}
    world
end

struct Rule
    name::Symbol
    rule::Function
    any_actor::Bool
end

struct AliasOf
    attr::Symbol
    other::Symbol
end

struct Rulebook
    name::Symbol
    default_outcome::Optional{Any}
    first_rules::Vector{Rule}
    rules::Vector{Rule}
    last_rules::Vector{Rule}
    variables::Dict{Symbol,Any}
end

struct Activity
    name::Symbol
    before_rules::Rulebook
    carry_out_rules::Rulebook
    after_rules::Rulebook
    activity_variables::Dict{Symbol, Any}
end

mutable struct Action
    name::Symbol
    understand_as::Vector{String}
    applies_to::Int64
    set_action_variables::Rulebook
    before_rules::Rulebook
    check_rules::Rulebook
    carry_out_rules::Rulebook
    report_rules::Rulebook
    action_variables::Dict{Symbol, Any}
end

@with_kw mutable struct World
    message_indent::Int64 = 0
    print_counter::Int64 = 0
    object_hierarchy::ObjectGraph = Dict()
    attribute_store::Dict{Symbol, Union{Dict{Symbol, Any}, Set{Any}}} = Dict()
    objects::Dict{Symbol, Object} = Dict()
    objects_by_type::Dict{Symbol, Vector{Object}} = Dict()
    rulebooks::Dict{Symbol, Rulebook} = Dict()
    title::String = ""
    active_rulebook::Vector{Rulebook} = []
    player::Optional{Object}
    first_room::Optional{Object}
    actions::Dict{Symbol, Action} = Dict()
    activities::Dict{Symbol, Activity} = Dict()
    room_descriptions::Symbol = :sometimes_abbreviated
    darkness_witnessed::Bool = false
end

struct DynamicString
    str::String
    subs::Dict{String,Function}
    world::World
    function DynamicString(w::World, str2::String, s::Dict)
        str = str2
        world = w
        subs = Dict()
        for k in keys(s)
            subs[string("{", k, "}")] = s[k]
        end
        new(str, subs)
    end
end

const Text = Union{String, DynamicString}


text:
# purpose: read ASCII values from a stream
function parse_digits(io::IO, c::UInt8, N::UInt8)
  v = zero(UInt32)
  p = 0x00989680
  i = 0x00
  m = 0x00000001

  while is_u8_digit(c)
    eof(io) && break
    if i < N
      p = div(p, 0x0a)
      i += 0x01
      v += (c-0x30)*p
      m = p
    end
    c = read(io, UInt8)
  end
  skip(io, -1)
  return v, m
end

function stream_int(io::IO, L::Int64)
  p = 10^(L-1)
  n = zero(Int64)
  @inbounds for i = 1:L
    c = read(io, UInt8)
    if is_u8_digit(c)
      n += p*(c-0x30)
    end
    p = div(p, 10)
  end
  return n
end

function stream_float(io::IO, c::UInt8)
  read_state = 0x00
  zz = 0x00000000
  p0 = 0x00989680
  v = zz
  m = 0x00000001
  d = zz
  x = zz
  q = zz
  msgn = true
  esgn = true
  read_exp = false
  while true
    c = read(io, UInt8)
    if is_u8_digit(c)
      if read_exp
        x,q = parse_digits(io, c, 0x04)
      else
        v,m = parse_digits(io, c, 0x08)
      end

    # '-'
    elseif c == 0x2d
      if read_exp
        esgn = false
      else
        msgn = false
      end

    # '+'
    elseif c == 0x2b
      continue

    # '.'
    elseif c == 0x2e
      c = read(io, UInt8)
      d, u = parse_digits(io, c, 0x08)

    # 'E', 'F', 'e', 'f'
    elseif c in (0x45, 0x46, 0x65, 0x66)
      read_exp = true
      sgn = false

    else
      break
    end
  end
  v1 = Float32(div(v, m)) + Float32(d) / Float32(p0)
  if q > zz
    if esgn
      v3 = 10.0f0^Float32(div(x,q))
    else
      v3 = 10.0f0^(-1.0f0*Float32(div(x,q)))
    end
  else
    v3 = one(Float32)
  end
  if msgn
    return v1 * v3
  else
    return -v1 * v3
  end
end

function stream_time(io::IO, T::Array{Y,1}) where Y<:Integer
  fill!(T, zero(Y))
  i = mark(io)
  k = 1
  while true
    c = read(io, UInt8)

    if c in (0x2d, 0x2e, 0x3a, 0x54)
      L = position(io)-i-1
      reset(io)
      T[k] = stream_int(io, L)
      skip(io, 1)
      i = mark(io)
      k += 1

    # exit on any non-digit character except (',', '.', ':', 'T')
    elseif c < 0x30 || c > 0x39
      L = position(io)-i-1
      if L > 0
        reset(io)
        T[k] = stream_int(io, L)
        skip(io, 1)
      end
      break
    end
  end
  y = T[1]
  if T[2] == zero(Y)
    return y2μs(y[1])
  else
    j = md2j(y, T[2], T[3])
    return mktime(y, j, T[4], T[5], T[6], T[7])
  end
end


text:
const compSize = 100
const boundary = 10 # Nodes with PS
const compDX = compSize * Float64.(collect(1:(compSize - 2)))
const dRdRMaxRMaxR = 1.0 / compSize / compSize

" Setup the parameters for the full TAM receptor model. "
function compParamm(params::Vector{T})::comprates{T} where {T}
    @assert all(params .>= 0.0)

    TAM, hetRs = paramTAMrate(view(params, 5:12))

    out = comprates{T}(TAM, params[1], params[2], params[3], params[4], hetRs)

    return detailedBalance(out)
end

function TAMreact(du::Vector, u::Vector, r::comprates, t; reaction = true)
    # If we're dealing with the PDE form
    if length(du) > 100
        @views for ii = 1:27
            duu = view(du, ii:27:length(du))
            uu = u[ii:27:end]

            if boundLigC[ii] == 0
                @. duu[2:(end - 1)] =
                    (-4.0 * uu[2:(end - 1)] + (2.0 - 1.0 / compDX) * uu[1:(end - 2)] + (2.0 + 1.0 / compDX) * uu[3:end]) / 2 / dRdRMaxRMaxR
            else
                @. duu[2:(boundary - 1)] =
                    (
                        -4.0 * uu[2:(boundary - 1)] +
                        (2.0 - 1.0 / compDX[1:(boundary - 2)]) * uu[1:(boundary - 2)] +
                        (2.0 + 1.0 / compDX[1:(boundary - 2)]) * uu[3:boundary]
                    ) / 2 / dRdRMaxRMaxR
                duu[boundary] = -4.0 * (uu[boundary] - uu[boundary - 1]) / dRdRMaxRMaxR

                duu[boundary + 1] = 4.0 * (uu[boundary + 2] - uu[boundary + 1]) / dRdRMaxRMaxR
                @. duu[(boundary + 2):(end - 1)] =
                    (
                        -4.0 * uu[(boundary + 2):(end - 1)] +
                        (2.0 - 1.0 / compDX[(boundary + 1):end]) * uu[(boundary + 1):(end - 2)] +
                        (2.0 + 1.0 / compDX[(boundary + 1):end]) * uu[(boundary + 3):end]
                    ) / 2 / dRdRMaxRMaxR

                # Diffusion into compartment
                duu[boundary + 1] -= uu[boundary + 1] / dRdRMaxRMaxR
                duu[boundary] += uu[boundary + 1] / dRdRMaxRMaxR
            end

            duu[1] = 4.0 * (uu[2] - uu[1]) / dRdRMaxRMaxR
            duu[end] = -4.0 * (uu[end] - uu[end - 1]) / dRdRMaxRMaxR
        end

        rmul!(du, r.diff)
    else
        fill!(du, 0.0)
    end

    if reaction
        @views for ii = 0:(Int(length(u) / 27) - 1)
            if ii < boundary
                gass = r.gasCur * r.gasPart
            else
                gass = r.gasCur
            end

            idx = (27 * ii + 1):(27 * ii + 27)
            compartmentReact(view(u, idx), view(du, idx), gass, nothing, r)
        end
    end
end


text:
module Learn_Discrete_BayesNet

# package code goes here

include("L_DVBN.jl")

end # module


text:
# This file includes an example file to plot the interpolation errors.

using FractalTools
using Makie
using StaticArrays
using GeometryBasics
import GeometryBasics: Ngon
using Random
Random.seed!(0)

# Construct ngon
ngon = Ngon(SVector([Point(BigFloat(cos(θ)), BigFloat(sin(θ))) for θ in 0 : π / 2 : 2π - π / 2]...))
f(x, y) = x^2 + y^2
pnts3d = getdata(f, ngon, 100)

# Construct interpolant
itp = interpolate(pnts3d)

# Calculate interpolation error
pnts2d1 = project(pnts3d)
fvals1 = map(pnt -> Point(pnt..., f(pnt...)), pnts2d1)
ivals1 = map(pnt -> Point(pnt..., itp(pnt...)), pnts2d1)
evals1 = map(pnt -> Point(pnt..., relerf(pnt...)), pnts2d1)

@show maximum(getindex.(evals1, 3))

pnts3d2 = getdata(f, ngon, 400)
pnts2d2 = project(pnts3d2)
fvals2 = map(pnt -> Point(pnt..., f(pnt...)), pnts2d2)
ivals2 = map(pnt -> Point(pnt..., itp(pnt...)), pnts2d2)
evals2 = map(pnt -> Point(pnt..., relerf(pnt...)), pnts2d2)

@show maximum(getindex.(evals2, 3))

# Plots
fig = Figure()
ax11 = LScene(fig)
ax12 = LScene(fig)
ax13 = LScene(fig)
ax21 = LScene(fig)
ax22 = LScene(fig)
ax23 = LScene(fig)
fig[1, 1] = vgrid!(hgrid!(ax11, ax12, ax13),hgrid!(ax21, ax22, ax23))

msh21 = project(triangulate(pnts3d)[2])
scatter!(ax11, pnts3d, color=:black)
wireframe!(ax11, msh21)
trisurf!(ax11, fvals1, meshcolor3=last.(fvals1))
trisurf!(ax12, ivals1, meshcolor3=last.(ivals1))
wireframe!(ax12, msh21)
trisurf!(ax13, evals1, meshcolor3=last.(evals1))
wireframe!(ax13, msh21)

msh22 = project(triangulate(pnts3d2)[2])
scatter!(ax21, pnts3d2, color=:black)
trisurf!(ax21, fvals2, meshcolor3=last.(fvals2))
wireframe!(ax21, msh22)
trisurf!(ax22, ivals2, meshcolor3=last.(ivals2))
wireframe!(ax22, msh22)
trisurf!(ax23, evals2, meshcolor3=last.(evals2))
wireframe!(ax23, msh22)
display(fig)


text:
import .OrdinaryDiffEq

struct DiffEqIntegrator{T<:AbstractScalarOrVec{<:AbstractFloat}, DiffEqSolver} <: AbstractLeapfrog{T}
    ϵ::T
    solver::DiffEqSolver
end

function step(
    integrator::DiffEqIntegrator,
    h::Hamiltonian,
    z::P,
    n_steps::Int=1;
    fwd::Bool=n_steps > 0,  # simulate hamiltonian backward when n_steps < 0
    res::Union{Vector{P}, P}=z) where {P<:PhasePoint}

    @unpack θ, r = z

    # For DynamicalODEProblem `u` is `θ` and `v` is `r`
    # f1 is dr/dt RHS function
    # f2 is dθ/dt RHS function
    v0, u0 = r, θ

    f1(v, u, p, t) = -∂H∂θ(h, u).gradient
    f2(v, u, p, t) = ∂H∂r(h, v)

    ϵ = fwd ? step_size(integrator) : -step_size(integrator)
    tspan = (0.0, sign(n_steps))
    problem = OrdinaryDiffEq.DynamicalODEProblem(f1, f2, v0, u0, tspan)
    diffeq_integrator = OrdinaryDiffEq.init(problem, integrator.solver, save_everystep=false, save_start=false, save_end=false, dt=ϵ)

    for i in 1:abs(n_steps)
        OrdinaryDiffEq.step!(diffeq_integrator)
        solution = diffeq_integrator.u.x  # (r, θ) at the proposed step
        z = phasepoint(h, solution[2], solution[1])
        !isfinite(z) && break
        if res isa Vector
            res[i] = z
        else
            res = z
        end
    end
    return res
end


text:
using Interpolations
include("UnitConversions.jl")

"""
    Fluid(gravityAPI::T, gasGravity::T, solutionGOR::T)

Creates a struct to represent a crude oil sample.
All properties must be given at standard conditions (60F, 14.5 psia)
"""
mutable struct Fluid{T}
    gravityAPI::T
    gasGravity::T
    solutionGOR::T
end

"""
    oilGravity(fluid::Fluid)

Computer the oil gravity using the API of the fluid.
"""
oilGravity(fluid::Fluid) = 141.5/(131.5+fluid.gravityAPI)

"""
    oilGravity(density::Float64)

Compute the oil gravity if density in lb/cuft is given
"""
oilGravity(density) = density/62.4

"""
    oilDensity(f::Fluid, temperature::Float64)

Compute the oil density in lb/cuft given a fluid and a temperature in Fahrenheit
"""
function oilDensity(fluid::Fluid, temperature)
    oilGrav = oilGravity(fluid)
    rs = fluid.solutionGOR
    gasGrav = fluid.gasGravity
    num = 62.4*oilGrav+0.0136*rs*gasGrav
    den = 0.972+0.000147*(rs*((gasGrav/oilGrav)^0.5)+(1.25*(temperature)))^1.175
    return num/den
end

"""
    adiabaticOilModulus(f::Fluid, temperature::Float64, pressure::Float64)

Compute the adiabatic oil modulus in psi using the Arco formula.
Temperature must be given in Fahrenheit.
Pressure must be given in psia
"""
function adiabaticOilBulkModulus(fluid::Fluid, temperature, pressure)
    a,b,c,d,e,f = 1.286e6, 13.55, 4.122e4, 4.53e3, 10.59, 3.228
    temp = fahrenhetToRankine(temperature)
    ka = a+b*pressure-c*(temp^0.5)-d*(fluid.gravityAPI)-e*(fluid.gravityAPI^2)+f*temp*(fluid.gravityAPI)
end

"""
    pressureGradient(density::Float64)

Compute the pressure gradient given a density in lb/cuft
"""
pressureGradient(density::Float64) = density/144.0

"""
    adiabaticWaterBulkModulus(temperature::Float, pressure::Float)

Calculates the adiabatic bulk modulus for water.
Temperature must be in Fahrenheit
Pressure must be in psia.
"""
function adiabaticWaterBulkModulus(temperature::Float64, pressure::Float64)
    tempArray = [32, 50, 68, 86, 104, 122, 140, 158, 176, 194, 212]
    kosArray = [289, 308, 323, 333, 340, 345, 348, 348, 341, 342, 336]
    ipt = LinearInterpolation(tempArray, kosArray)
    ko = ipt(temperature)
    return 1000*ko+3.4*pressure
end

"""
    gasDensity(fluid::Fluid, pressure::Float64, temperature::Float64)

Calculates gas density in lb/cuft.
pressure must be in psia
temperature must be in Fahrenheit
"""
function gasDensity(fluid::Fluid, pressure::Float64, temperature::Float64)
        ma = 28.96*fluid.gasGravity
        temp = fahrenhetToRankine(temperature)
        return (pressure*ma)/(10.73*temp)
end

"""
    mixtureDensity(oilDensity, waterDensity, waterFraction)
"""
function mixtureDensity(oilDensity::Float64, waterDensity::Float64, waterFraction::Float64)
    return oilDensity*(1-waterFraction)+waterDensity*waterFraction
end

"""
    mixtureBulkModulus(oilDensity, waterDensity, waterFraction)
"""
function mixtureBulkModulus(oilBulkModulus::Float64, waterBulkModulus::Float64, waterFraction::Float64)
    return oilBulkModulus*(1-waterFraction)+waterBulkModulus*waterFraction
end

max_stars_repo_name: leytzher/Shot

max_issues_repo_name: leytzher/Shot

max_forks_repo_name: leytzher/Shot


text:
"""Computes the value of the function ackerman(n,m)"""
function ack(m :: Int64, n :: Int64) :: Int64
    if m == 0
        return n + 1
    elseif n == 0
        return ack(m - 1, 1)
    else
        return ack(m - 1, ack(m, n - 1))
    end
end