text:
using Stats
using Base.Test

## moments

wv = weights(ones(5) * 2.0)

@test_approx_eq skewness(1:5) 0.0
@test_approx_eq skewness([1, 2, 3, 4, 5]) 0.0
@test_approx_eq skewness([1, 2, 2, 2, 5]) 1.1731251294063556
@test_approx_eq skewness([1, 4, 4, 4, 5]) -1.1731251294063556

@test_approx_eq skewness([1, 2, 2, 2, 5], wv) 1.1731251294063556

@test_approx_eq kurtosis(1:5) -1.3
@test_approx_eq kurtosis([1, 2, 3, 4, 5]) -1.3
@test_approx_eq kurtosis([1, 2, 3, 3, 2]) -1.1530612244897953

@test_approx_eq kurtosis([1, 2, 3, 4, 5], wv) -1.3

## variability

@test_approx_eq variation([1:5]) 0.527046276694730

@test_approx_eq sem([1:5]) 0.707106781186548

@test_approx_eq mad([1:5], 3) 1.4826
@test_approx_eq mad([1:5]) 1.4826
@test_approx_eq mad(1:5) 1.4826

# min/max related

@test minmax(1:5) == (1,5)
@test minmax([3, 2, 7, 1, 4]) == (1,7)
@test minmax([3.0, NaN, 1.0, NaN, 2.0]) == (1.0, 3.0)

@test_approx_eq midrange([1, 2, 3, 4, 5]) 3.0
@test_approx_eq midrange([1, 2, 2, 2, 5]) 3.0
@test_approx_eq midrange([1, 4, 4, 4, 5]) 3.0
@test_approx_eq midrange([1:5]) 3.0

@test_approx_eq range([1, 2, 3, 4, 5]) 4.0
@test_approx_eq range([1, 2, 2, 2, 5]) 4.0
@test_approx_eq range([1, 4, 4, 4, 5]) 4.0
@test_approx_eq range([1:15]) 14.0

# quantile & friends

@test_approx_eq quantile(1:5) [1:5]
@test_approx_eq nquantile(1:5, 2) [1, 3, 5]
@test_approx_eq nquantile(1:5, 4) [1:5]

@test_approx_eq percentile([1:5], 25) 2.0
@test_approx_eq percentile([1:5], [25, 50, 75]) [2.0, 3.0, 4.0]
@test_approx_eq iqr(1:5) 2.0

# mode & modes

@test mode([1, 2, 3, 3, 2, 2, 1], 1:3) == 2
@test modes([1, 2, 3, 3, 2, 2, 1], 1:3) == [2]
@test modes([1, 3, 2, 3, 3, 2, 2, 1], 1:3) == [2, 3]

@test mode([1, 2, 3, 3, 2, 2, 1]) == 2
@test modes([1, 2, 3, 3, 2, 2, 1]) == [2]
@test sort(modes([1, 3, 2, 3, 3, 2, 2, 1])) == [2, 3]

# summarystats

s = summarystats(1:5)
@test isa(s, Stats.SummaryStats)
@test s.min == 1.0
@test s.max == 5.0
@test_approx_eq s.mean 3.0
@test_approx_eq s.median 3.0
@test_approx_eq s.q25 2.0
@test_approx_eq s.q75 4.0


text:
module Outputs

export SubroutineError

struct SubroutineError
    name::String
    cerr::String
    msg::String
end

include("PWscf.jl")
include("PHonon.jl")

end

ext: jl


text:
export HierarchicalVector, partition, nlevels

partition(x)=x
nlevels(A)=0

##
# Represent a binary hierarchical vector
##

abstract type AbstractHierarchicalArray{SV,T,HS,N} <: AbstractArray{T,N} end

const AbstractHierarchicalVector{S,T,HS} = AbstractHierarchicalArray{S,T,HS,1}

mutable struct HierarchicalVector{S,T,HS} <: AbstractHierarchicalVector{S,T,HS}
    data::HS
    HierarchicalVector{S,T,HS}(data::HS) where {S,T,HS} = new{S,T,HS}(data)
end

HierarchicalVector(data::Tuple{S1,S2}) where {S1,S2} = HierarchicalVector{promote_type(S1,S2),promote_type(eltype(S1),eltype(S2)),Tuple{S1,S2}}(data)
HierarchicalVector(data::Tuple{HierarchicalVector{S1,T1,HS1},HierarchicalVector{S2,T2,HS2}}) where {S1,S2,T1,T2,HS1,HS2} = HierarchicalVector{promote_type(S1,S2),promote_type(eltype(S1),eltype(S2),T1,T2),Tuple{HierarchicalVector{S1,T1,HS1},HierarchicalVector{S2,T2,HS2}}}(data)
HierarchicalVector(data::Tuple{S1,HierarchicalVector{S,T,HS}}) where {S,S1,T,HS} = HierarchicalVector{promote_type(S,S1),promote_type(eltype(S),eltype(S1),T),Tuple{S1,HierarchicalVector{S,T,HS}}}(data)
HierarchicalVector(data::Tuple{HierarchicalVector{S,T,HS},S1}) where {S,S1,T,HS} = HierarchicalVector{promote_type(S,S1),promote_type(eltype(S),eltype(S1),T),Tuple{HierarchicalVector{S,T,HS},S1}}(data)

HierarchicalVector(data::Vector) = HierarchicalVector(data,round(Int,log2(length(data))))

function HierarchicalVector(data::Vector,n::Int)
    @assert length(data) == 2^n
    if n == 1
        return HierarchicalVector(tuple(data...))
    elseif n ≥ 2
        return HierarchicalVector((HierarchicalVector(data[1:2^(n-1)],n-1),HierarchicalVector(data[1+2^(n-1):end],n-1)))
    end
end

Base.similar(H::HierarchicalVector) = HierarchicalVector(map(similar,data(H)))
Base.similar(H::HierarchicalVector{SS,V,T}, ::Type{S}) where {SS,V,T,S} =
    HierarchicalVector(map(A->similar(A,S),data(H)))
Base.similar(H::HierarchicalVector{SS,V,T}, S) where {SS,V,T} =
    HierarchicalVector(map(A->similar(A,S),data(H)))

data(H::HierarchicalVector) = H.data

nlevels(H::HierarchicalVector) = 1+mapreduce(nlevels,max,data(H))

partition(H::HierarchicalVector) = data(H)

function partition(H::Vector{HierarchicalVector})
    n = length(H)
    H11,H12 = partition(H[1])
    H1,H2 = fill(H11,n),fill(H12,n)
    for i=1:n
        H1[i],H2[i] = partition(H[i])
    end
    H1,H2
end

collectdata(H::HierarchicalVector{S,T,NTuple{2,S}}) where {S,T} = collect(data(H))
collectdata(H::HierarchicalVector{S,T,Tuple{S,HierarchicalVector{S,T,HS}}}) where {S,T,HS} = vcat(H.data[1],collectdata(H.data[2]))
collectdata(H::HierarchicalVector{S,T,Tuple{HierarchicalVector{S,T,HS},S}}) where {S,T,HS} = vcat(collectdata(H.data[1]),H.data[2])
function collectdata(H::HierarchicalVector{S}) where S
    ret = S[]
    append!(ret,mapreduce(collectdata,vcat,data(H)))
    ret
end

convert(::Type{HierarchicalVector{S,T,HS}},M::HierarchicalVector) where {S,T,HS} = HierarchicalVector(convert(Vector{S},collectdata(M)))
Base.promote_rule(::Type{HierarchicalVector{S,T,HS}},::Type{HierarchicalVector{SS,TT,HSS}}) where {S,T,HS,SS,TT,HSS}=HierarchicalVector{promote_type(S,SS),promote_type(T,TT),promote_type(HS,HSS)}

function Base.size(H::HierarchicalVector{S}) where S<:AbstractVector
    H1,H2 = data(H)
    (size(H1)[1]+size(H2)[1],)
end

function Base.size(H::HierarchicalVector{S}) where S<:Number
    H1,H2 = data(H)
    if typeof(H1) <: S && typeof(H2) <: S
        return (2,)
    elseif typeof(H1) <: S
        return (1+size(H2)[1],)
    elseif typeof(H2) <: S
        return (size(H1)[1]+1,)
    else
        return (size(H1)[1]+size(H2)[1],)
    end
end

function Base.getindex(H::HierarchicalVector{S},i::Int) where S<:Union{Number,AbstractVector}
    H1,H2 = data(H)
    m1,m2 = length(H1),length(H2)
    if 1 ≤ i ≤ m1
        return H1[i]
    elseif m1 < i ≤ m1+m2
        return H2[i-m1]
    else
        throw(BoundsError())
    end
end

function Base.setindex!(H::HierarchicalVector{S},v,i::Int) where S<:Union{Number,AbstractVector}
    H1,H2 = data(H)
    m1,m2 = length(H1),length(H2)
    if 1 ≤ i ≤ m1
        H1[i] = v
    elseif m1 < i ≤ m1+m2
        H2[i-m1] = v
    else
        throw(BoundsError())
    end
end

Base.getindex(H::HierarchicalVector{S},ir::AbstractRange) where {S<:Union{Number,AbstractVector}} =
    eltype(H)[H[i] for i=ir]
Base.full(H::HierarchicalVector{S}) where {S<:Union{Number,AbstractVector}}=H[1:size(H,1)]

# algebra

for op in (:+,:-)
    @eval begin
        $op(a::Bool,H::HierarchicalVector{S,Bool}) where {S} = error("Not callable")
        $op(H::HierarchicalVector{S,Bool},a::Bool) where {S} = error("Not callable")

        $op(H::HierarchicalVector) = HierarchicalVector(map($op,data(H)))

        $op(H::HierarchicalVector,J::HierarchicalVector) = HierarchicalVector(map($op,data(H),data(J)))
        $op(H::HierarchicalVector,J::Vector) = $op(full(H),J)
        $op(J::Vector,H::HierarchicalVector) = $op(J,full(H))
        broadcasted(::DefaultArrayStyle{1}, ::typeof($op), H::HierarchicalVector, a::Number) =
            HierarchicalVector(($op.(H.data[1],a), $op.(H.data[2],a)))
        broadcasted(::DefaultArrayStyle{1}, ::typeof($op), a::Number, H::HierarchicalVector) =
            HierarchicalVector(($op.(a,H.data[1]), $op.(a,H.data[2])))
    end
end

*(H::HierarchicalVector,a::Number) = HierarchicalVector((H.data[1]*a,H.data[2]*a))
*(a::Number,H::HierarchicalVector) = H*a

for op in (:(LinearAlgebra.dot),:dotu)
    @eval begin
        $op(H::HierarchicalVector,J::HierarchicalVector) = $op(H.data[1],J.data[1])+$op(H.data[2],J.data[2])
        $op(H::HierarchicalVector,J::Vector) = $op(full(H),J)
        $op(J::Vector,H::HierarchicalVector) = $op(J,full(H))
    end
end

Base.cumsum(H::HierarchicalVector) = HierarchicalVector((cumsum(H.data[1]),sum(H.data[1])+cumsum(H.data[2])))
Base.conj!(H::HierarchicalVector) = (map(conj!,data(H));H)
Base.copy!(H::HierarchicalVector,J::HierarchicalVector) = (map(copy!,data(H),data(J));H)

for op in (:(Base.zero),:(Base.ones),:(Base.abs),:(Base.abs2),:(Base.conj),:(Base.copy))
    @eval begin
        $op(H::HierarchicalVector) = HierarchicalVector(map($op,data(H)))
    end
end

mul!(b::Vector,A::Matrix,H::HierarchicalVector) = mul!(b,A,full(H))
mul!(b::Vector,A::LowRankMatrix,H::HierarchicalVector) = mul!(b,A,full(H))

#=
        function $op{S,T,HS,SS,TT,HSS}(H::Array{HierarchicalVector{S,T,HS}}, J::Array{HierarchicalVector{SS,TT,HSS}})
            ret = similar(H)
            for i in eachindex(H,J)
                @inbounds ret[i] = $op(H[i], J[i])
            end
            ret
        end
        function $op{S,T,HS,V}(H::Array{HierarchicalVector{S,T,HS}}, J::Array{V})
            ret = similar(H)
            for i in eachindex(H,J)
                @inbounds ret[i] = $op(H[i], J[i])
            end
            ret
        end
        $op{S,T,HS,V}(J::Array{V}, H::Array{HierarchicalVector{S,T,HS}}) = $op(H,J)
=#


text:
using Luna
import Test: @test, @testset

radius = 125e-6 # HCF core radius
flength = 3 # HCF length

gas = :Ar
pressure = 80e-3 # gas pressure in bar

λ0 = 800e-9 # central wavelength of the pump pulse
τfwhm = 10e-15 # FWHM duration of the pump pulse
energy = 60e-6 # energy in the pump pulse

τ0 = Tools.τfw_to_τ0(τfwhm, :gauss)

out = prop_capillary(radius, flength, gas, pressure; λ0, τfwhm, energy,
                     modes=:HE11, trange=1000e-15, λlims=(150e-9, 4e-6), kerr=false, plasma=false)

t, Et = Processing.getEt(out)
Ptend = abs2.(Et[:, end]) # pulse profile at end

fwend = Maths.fwhm(t, Ptend; method=:spline) # FWHM duration at end
τ0end = Tools.τfw_to_τ0(fwend, :gauss) # Gaussian τ0 duration at end

mode = Capillary.MarcatiliMode(radius, gas, pressure) # propagation mode
β2 = Modes.dispersion(mode, 2, PhysData.wlfreq(λ0)) # dispersion of the mode
GDD = β2*flength # total GDD accumulated

τ0calc = sqrt((GDD/τ0)^2 + τ0^2) # analytical formula for Gaussian duration after GDD

@test isapprox(τ0calc, τ0end; rtol=1e-2)

energy_calc = Capillary.transmission(radius, λ0, flength)*energy # calculated energy at end
@test isapprox(energy_calc, Processing.energy(out)[end]; rtol=1e-2)


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

# source path in tasks
path = Base.source_path()::String # this variable is leaked to the source script
@test endswith(path, joinpath("test","test_sourcepath.jl"))
@test let ct = current_task()
    yieldto(@task yieldto(ct, Base.source_path()))
end == path
@test let ct = current_task()
    yieldto(@task schedule(ct, Base.source_path()))
end == path
@test let ct = current_task(), t = @task Base.source_path()
    schedule(ct)
    yieldto(t)
    fetch(t)
end == path
@test isabspath(@__FILE__)

max_stars_repo_name: DragonXDoom/julia

max_issues_repo_name: DragonXDoom/julia

max_forks_repo_name: DragonXDoom/julia


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

mutable struct RateLimiter <: SwaggerModel
    bandwidth::Union{ Nothing, TokenBucket } # bandwidth
    ops::Union{ Nothing, TokenBucket } # ops

    function RateLimiter(;bandwidth=nothing, ops=nothing)
        o = new()
        set_field!(o, :bandwidth, bandwidth)
        set_field!(o, :ops, ops)
        o
    end
end # type RateLimiter

const _name_map_RateLimiter = Dict{String,Symbol}(["bandwidth"=>:bandwidth, "ops"=>:ops])
const _field_map_RateLimiter = Dict{Symbol,String}([:bandwidth=>"bandwidth", :ops=>"ops"])
Swagger.name_map(::Type{ RateLimiter }) = _name_map_RateLimiter
Swagger.field_map(::Type{ RateLimiter }) = _field_map_RateLimiter

function check_required(o::RateLimiter)
    true
end

function validate_field(o::RateLimiter, name::Symbol, val)
end


text:
module SolverBenchmarking

export summarize_results, run_solver_benchmark, BenchmarkParams

using JSON, DataFrames, CSV
using MINLPLib, JuMP, EAGO

"""
Holds parameters for solver benchmarking.
"""
Base.@kwdef mutable struct BenchmarkParams
    time::Float64 = 200.0
    threads_per::Int = 1
    rerun::Bool = true
    has_obj_bnd::Bool = false
end

"""
Custom error type to handle exceptions
and log information that occur when solving.
"""
struct InstanceError <: Exception
    path::String
    solver_name::String
    model_name::String
end

"""
Defines storage routine when an instance error occurs.
"""
function store_error(d::InstanceError)
    data = Dict{String,Any}()
    data["SolverName"] = d.solver_name
    data["InstanceName"] = d.model_name
    data["TerminationStatus"] = string(MOI.OTHER_ERROR)
    data["PrimalStatus"] = string(MOI.OTHER_RESULT_STATUS)
    data["ObjBound"] = 0.0
    data["ObjValue"] = 0.0
    data["SolveTime"] = 200.0
    data["CompletedSolveTime"] = 200.0
    data["Error"] = true

    file_path = joinpath(d.path, "$(d.solver_name)_$(d.model_name).json")
    open(file_path,"w") do f
        JSON.print(f, data)
    end
    return
end

#=
function capture_stdout(f)
    stdout_orig = stdout
    (rd, wr) = redirect_stdout()
    f()
    close(wr)
    redirect_stdout(stdout_orig)
    read(rd, String)
end
=#

"""
Fetches the model at "lib\\model_name", attaches the optimizers factory "opt_factory", sets benchmark parameters,
solves the optimization problem, then tidies up the memory space by finalizing the background then garbage collecting.
"""
function fetch_and_optimize(lib::String, model_name::String, opt_factory::Any, params::BenchmarkParams)
    m = fetch_model(lib, model_name)
    set_optimizer(m, opt_factory)

    set_time_limit_sec(m, params.time)
    set_silent(m)

    t = time()
    try
        JuMP.optimize!(m)
    catch
    end
    t = time() - t

    finalize(backend(m))
    GC.gc()
    return m, t
end

"""
If the rerun parameter is set to `true` or no JSON file ("p\\\$(solver_name)_\$(model_name).json") exists
then fetch the model, optimize it and store results to a JSON file at the path ("p\\\$(solver_name)_\$(model_name).json").
"""
function benchmark_problem(opt_factory::Any, solver_name::String, lib::String, libout::String, model_name::String, p::String, params::BenchmarkParams)
    file_path = joinpath(p, libout, "$(solver_name)_$(model_name).json")
    if !isfile(file_path) || params.rerun
        m, t = fetch_and_optimize(lib, model_name, opt_factory, params)

        obj_val = Inf
        s_time = params.time
        t_status = MOI.OPTIMIZE_NOT_CALLED
        p_status = MOI.NO_SOLUTION
        try
            obj_val = objective_value(m)
            s_time = solve_time(m)
            t_status = termination_status(m)
            p_status = primal_status(m)
        catch e
            obj_val = Inf
            s_time = params.time
            t_status = MOI.OPTIMIZE_NOT_CALLED
            p_status = MOI.NO_SOLUTION
        end

        obj_bnd = -Inf
        try
            obj_bnd = objective_bound(m)
        catch e
            obj_bnd = -Inf
        end

        data = Dict{String,Any}()
        data["SolverName"] = solver_name
        data["InstanceName"] = model_name
        data["TerminationStatus"] = string(t_status)
        data["PrimalStatus"] = string(p_status)
        data["ObjBound"] = params.has_obj_bnd ? obj_bnd : -Inf
        data["ObjValue"] = obj_val
        data["SolveTime"] = s_time
        data["SolveTimeRaw"] = t
        data["CompletedSolveTime"] = (t_status == MOI.OPTIMAL) ? s_time : 100.0
        data["Error"] = false

        open(file_path,"w") do f
            JSON.print(f, data)
        end
    end
    return
end

"""
Main routine for running solver benchmark.
- If no directory `p` exists, make it.
- If no directory `p\\lib` exists, make it.
- Write the parameter file to the file path.
For each solver factory in "s", and each problem in "lib".
"""
function run_solver_benchmark(p::String, s::Dict{String,Any}, lib::String, libout::String; params = BenchmarkParams())

    # store parameter file
    !isdir(p) && mkdir(p)
    !isdir(joinpath(p, lib)) && mkdir(joinpath(p, lib))
    file_path = joinpath(p, libout, "parameters.json")
    data = Dict{String,Any}()
    data["time"] = params.time
    data["threads_per"] = params.threads_per
    data["rerun"] = params.rerun
    open(file_path,"w") do f
        JSON.print(f, data)
    end

    # benchmark solvers
    n = MINLPLib.fetch_names(lib)
    for (k,v) in s
        println("Running all benchmarks problems for solver = $(k)")
        @show length(n)
        for ni in n
            println("Running problem = $ni in $lib")
            benchmark_problem(v, k, lib, libout, ni, p, params)
        end
    end
    println("Benchmark Suite Complete")
    return
end

"""
"""
function summarize_results(lib::String, p::String)
    folder_path = joinpath(p, lib)
    df = DataFrame(SolverName = String[], InstanceName = String[], TerminationStatus = String[], PrimalStatus = String[],
                   ObjBound = Float64[], ObjValue = Float64[], SolveTime = Float64[], SolveTimeRaw = Float64[],
                   CompletedSolveTime = Float64[], Error = Bool[])
    for file in readdir(folder_path)
        if !endswith(file, ".json") || file === "parameters.json"
            continue
        end
        d = joinpath(folder_path, file) |> open |> JSON.parse
        obj_bnd = isnothing(d["ObjBound"]) ? -Inf : d["ObjBound"]
        obj_val = isnothing(d["ObjValue"]) ? Inf : d["ObjValue"]
        push!(df, (d["SolverName"], d["InstanceName"], d["TerminationStatus"], d["PrimalStatus"], obj_bnd, obj_val, d["SolveTime"], d["SolveTimeRaw"], d["CompletedSolveTime"], d["Error"]))
    end

    param_dict = joinpath(p, lib, "parameters.json") |> open |> JSON.parse
    df.RelativeGap = abs.(df.ObjValue - df.ObjBound)./min.(abs.(df.ObjValue), abs.(df.ObjBound))
    df.AbsoluteGap = df.ObjValue - df.ObjBound
    df.SolvedInTime = df.SolveTime .< param_dict["time"]
    CSV.write(joinpath(folder_path, "result_summary.csv"), df)
    return
end

end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
Sets of Prinz conductances from 2003 and 2004 papers
"""

# parameters from Prinz 2003 https://journals.physiology.org/doi/pdf/10.1152/jn.00641.2003, figure 3.
burster_channels = [Prinz.Na(400.0 * Prinz_conv), Prinz.CaS(8.0 * Prinz_conv), Prinz.CaT(0.0 * Prinz_conv), Prinz.H(0.04 * Prinz_conv),
    Prinz.Ka(50.0 * Prinz_conv), Prinz.KCa(20.0 * Prinz_conv), Prinz.Kdr(50.0 * Prinz_conv), Prinz.leak(0.0)]

# parameters from Xolotl, which in turn come from Prinz 2004 "Similar network activity from disparate circuit parameters"
AB2_ch = [Prinz.Na(100 * Prinz_conv), Prinz.CaS(6 * Prinz_conv), Prinz.CaT(2.5 * Prinz_conv), Prinz.H(0.01 * Prinz_conv),
    Prinz.Ka(50 * Prinz_conv), Prinz.KCa(5 * Prinz_conv), Prinz.Kdr(100 * Prinz_conv), Prinz.leak(0.0)]
LP4_ch = [Prinz.Na(100 * Prinz_conv), Prinz.CaS(4 * Prinz_conv), Prinz.CaT(0.0), Prinz.H(0.05 * Prinz_conv), Prinz.Ka(20 * Prinz_conv),
    Prinz.KCa(0.0), Prinz.Kdr(2.5 * Prinz_conv), Prinz.leak(0.03 * Prinz_conv)]
PY1_ch = [Prinz.Na(100 * Prinz_conv), Prinz.CaS(2 * Prinz_conv), Prinz.CaT(2.4 * Prinz_conv), Prinz.H(0.05 * Prinz_conv), Prinz.Ka(50 * Prinz_conv),
    Prinz.KCa(0.0), Prinz.Kdr(125 * Prinz_conv), Prinz.leak(0.01 * Prinz_conv)]

# parameters from wiktor phillip's Conductor demos
AB_channels = [Prinz.Na(100.0 * Prinz_conv), Prinz.CaS(3.0 * Prinz_conv), Prinz.CaT(1.3 * Prinz_conv), Prinz.H(0.5 * Prinz_conv),
    Prinz.Ka(5.0 * Prinz_conv), Prinz.KCa(10.0 * Prinz_conv), Prinz.Kdr(20.0 * Prinz_conv), Prinz.leak(0.01 * Prinz_conv)]

"""
Sets of Liu conductances from various sources.
"""

# parameters from Drion et. al. 2015 "Dynamic Input Conductances Shape Neuronal Spiking" https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4586923/
Burster1 = [Liu.Na(700.0 * Liu_conv), Liu.CaS(4.0 * Liu_conv), Liu.CaT(2.0 * Liu_conv), Liu.Ka(50.0 * Liu_conv), Liu.KCa(40.0 * Liu_conv),
    Liu.Kdr(70.0 * Liu_conv), Liu.H(0.03 * Liu_conv), Liu.leak(0.01 * Liu_conv)]

Burster2 = [Liu.Na(700.0 * Liu_conv), Liu.CaS(2.25 * Liu_conv), Liu.CaT(6.25 * Liu_conv), Liu.Ka(85.0 * Liu_conv), Liu.KCa(50.0 * Liu_conv),
    Liu.Kdr(90.0 * Liu_conv), Liu.H(0.03 * Liu_conv), Liu.leak(0.01 * Liu_conv)]

# parameters from O'Leary (2014) 10.1016/j.neuron.2014.04.002
# Burster. These don't have to multiplied by Liu_conv
AB_channels = [Liu.Na(1831.0), Liu.CaS(27.0), Liu.CaT(23.0), Liu.H(10.1), Liu.Ka(246.0), Liu.KCa(980.0), Liu.Kdr(610.0), Liu.leak(0.99)]

"""
Hodgkin-Huxley model conductances
tonic spiking when an applied current Iapp>2 is used
"""
HH = [Liu.Na(1000), Liu.leak(0.1), Liu.Kdr(200)]


text:
"""
    Articulations{T}

A state type for the depth-first search that finds the articulation points in a graph.
"""
mutable struct Articulations{T<:Integer}
    low::Vector{T}
    depth::Vector{T}
    articulation_points::BitArray
    id::T
end

function Articulations(g::AbstractGraph)
    n = nv(g)
    T = typeof(n)
    return Articulations(zeros(T, n), zeros(T, n), falses(n), zero(T))
end

"""
    visit!(state, g, u, v)

Perform a depth first search storing the depth (in `depth`) and low-points
(in `low`) of each vertex.
Call this function repeatedly to complete the DFS (see [`articulation`](@ref) for usage).
"""
function visit!(state::Articulations, g::AbstractGraph, u::Integer, v::Integer)
    children = 0
    state.id += 1
    state.depth[v] = state.id
    state.low[v] = state.depth[v]

    for w in out_neighbors(g, v)
        if state.depth[w] == 0
            children += 1
            visit!(state, g, v, w)

            state.low[v] = min(state.low[v], state.low[w])
            if state.low[w] >= state.depth[v] && u != v
                state.articulation_points[v] = true
            end

        elseif w != u
            state.low[v] = min(state.low[v], state.depth[w])
        end
    end
    if u == v && children > 1
        state.articulation_points[v] = true
    end
end

"""
    articulation(g)

Compute the [articulation points](https://en.wikipedia.org/wiki/Biconnected_component)
of a connected graph `g` and return an array containing all cut vertices.
"""
function articulation(g::AbstractGraph)
    state = Articulations(g)
    for u in vertices(g)
        if state.depth[u] == 0
            visit!(state, g, u, u)
        end
    end
    return find(state.articulation_points)
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
module potentials

export Potential, BoxPotential

type Potential
    A::Function
    B::Function
    E::Function
    basis_function::Function
end

fcos(s, n::Int, i::Int, j::Int) = cospi((i+j)/2)*sinc((i+j+2*n*s)/2)
fsin(s, n::Int, i::Int, j::Int) = sinpi((i+j+n)/2)*sinc((i+j+n*s)/2)

function A(s, n::Int, i::Int, j::Int)
    @assert n>=0
    @assert i>0
    @assert j>0
    # i += 1
    # j += 1
    # n += 1
    k = (i==j ? 0.5 : 0.)
    k + (-1)^n/4 * (fcos(s,n,i,j) + fcos(s,-n,i,j) - fcos(s,n,i,-j) - fcos(s,-n,i,-j))
end

function B(s, n::Int, i::Int, j::Int)
    @assert n>=0
    @assert i>0
    @assert j>0
    # i += 1
    # j += 1
    # n += 1
    return 0.5 * (-fsin(s,n,i,j) + fsin(s,n,-i,j) + fsin(s,n,i,-j) + fsin(s,-n,i,+j))
end

E(E0, i::Int) = E0*i^2

function basis_function(index::Int, x)
    k = pi*index/2
    if (x<-1 || x>1)
        return zero(x)
    else
        return 1./sqrt(2)*sin(k*(x+1))
    end
end

basis_function(index::Int, x::Vector) = map(xi->basis_function(index, xi), x)

const BoxPotential = Potential(A, B, E, basis_function)

end


text:
# indicator of a PSD

"""
  IndPSD()

Returns the function `g = ind{A : A ⪰ 0}`, i.e. the indicator of the positive semidefinite cone.
The argument to the function can be either a Symmetric or Hermitian object.
From Julia 0.5, the argument can also be an AbstractVector{Float64} holding a symmetric matrix in (lower triangular) packed storage.
"""

immutable IndPSD <: IndicatorConvexCone end

function (f::IndPSD){T <: RealOrComplex}(X::HermOrSym{T})
  F = eigfact(X);
  for i in eachindex(F.values)
    #Do we allow for some tolerance here?
    if F.values[i] <= -1e-14
      return +Inf
    end
  end
  return 0.0
end

function prox!{T <: RealOrComplex}(f::IndPSD, X::HermOrSym{T}, Y::HermOrSym{T}, gamma::Real=1.0)
  n = size(X, 1);
  F = eigfact(X);
  for i in eachindex(F.values)
    F.values[i] = max.(0.0, F.values[i]);
  end
  for i = 1:n
    for j = 1:n
      Y.data[i,j] = 0.0
      for k = 1:n
        Y.data[i,j] += F.vectors[i,k]*F.values[k]*F.vectors[j,k]
      end
    end
  end
  return 0.0
end

fun_name(f::IndPSD) = "indicator of positive semidefinite cone"
fun_dom(f::IndPSD) = "Symmetric, Hermitian, AbstractArray{Float64}"
fun_expr(f::IndPSD) = "x ↦ 0 if A ⪰ 0, +∞ otherwise"
fun_params(f::IndPSD) = "none"

function prox_naive{T <: RealOrComplex}(f::IndPSD, X::HermOrSym{T}, gamma::Real=1.0)
  F = eigfact(X);
  return F.vectors * diagm(max.(0.0, F.values)) * F.vectors', 0.0;
end

### Below: with AbstractVector argument

function (f::IndPSD){T <: Float64}(x::AbstractVector{T})
  Z = dspev!('N', 'L', copy(x))
  for i in 1:length(Z)
    #Do we allow for some tolerance here?
    if Z[i] <= -1e-14
      return +Inf
    end
  end
  return 0.0
end

function prox!(f::IndPSD, x::AbstractVector{Float64}, y::AbstractVector{Float64}, gamma::Real=1.0)
  y[:] = x              # Copy x since dspev! corrupts input
  (W, Z) = dspevV!('L', y)
  W = max.(W, 0.0)         # NonNeg eigenvalues
  M = Z.*W'             # Equivalent to Z*diagm(W) without constructing W matrix
  M = M*Z'              # Now let M = Z*diagm(W)*Z'
  n = length(W)
  k = 1
  for j in 1:n, i in j:n  # Store lower diagonal of M in y
    y[k] = M[i,j]
    k = k+1
  end
  return 0.0
end

function prox_naive{T<:Float64}(f::IndPSD, x::AbstractVector{T}, gamma::Real=1.0)
  n = Int(sqrt(1/4+2*length(x))-1/2)  # Formula for size of matrix
  X = Array{T,2}(n,n)
  k = 1
  for j = 1:n, i = j:n                # Store y in M
    X[i,j] = x[k]                   # Lower half
    if i != j
      X[j,i] = x[k]               # Strictly upper half
    end
    k = k+1
  end
  X, v = prox_naive(f, Symmetric(X), gamma)
  y = similar(x)
  k = 1
  for j = 1:n, i = j:n                # Store Lower half of X in y
    y[k] = X[i,j]
    k = k+1
  end
  return y, 0.0
end


text:
__precompile__(true)
module ParallelPOMCPOW

importall POMDPs
using Parameters
using POMCPOW
using POMDPToolbox
using StaticArrays

export
    ParallelPOMCPOWSolver

struct ParallelPOMCPOWSolver
    powsolver::POMCPOWSolver
    batch_size::Int
end

struct ParallelPOMCPOWPlanner{P<:POMCPOWPlanner}
    powplanner::P
    batch_size::Int
end

function solve(s::ParallelPOMCPOWSolver, m::POMDP)
    powplanner = solve(s.powsolver, m)
    return ParallelPOMCPOWPlanner(powplanner, s.batch_size)
end

include("belief.jl")
include("message.jl")
include("worker.jl")

action(p::ParallelPOMCPOWPlanner, b) = first(action_info(p, b))

function POMDPToolbox.action_info(p::ParallelPOMCPOWPlanner, b)
    pp = p.powplanner
    ps = pp.solver
    @assert ps.eps == 0.01 # not supported; should be default
    @assert !ps.enable_action_pw # action progressive widening not supported
    m = pp.problem
    info = Dict{Symbol, Any}()

    acts = SVector(collect(iterator(actions(m, b)))...)
    na = length(acts)

    report_channel = RemoteChannel(()->Channel{WorkReport{na}}(10*nprocs()))
    if nprocs() > 1
        workers = 2:nprocs()
    else
        workers = [1]
    end
    # should these channels be on 1 or on the remote? I think on 1 since it will be idle the most?
    update_channels = Dict(w => RemoteChannel(()->Channel{GlobalUpdate{na}}(10*nprocs()), 1) for w in workers)

    futures = Dict(w => remotecall(run_worker,
                                   w, p, b, acts,
                                   update_channels[w],
                                   report_channel) for w in workers)

    qs = @SVector(zeros(Float64, na)) # XXX should be init_q
    adjusted_ns = @SVector(zeros(Int, na))
    worker_ns = Dict(w=>@SVector(zeros(Int, na)) for w in workers) # number of simulations for that action
    worker_n_obs = Dict(w=>@SVector(ones(Int, na)) for w in workers)
    worker_qs = Dict(w=>@SVector(zeros(Float64, na)) for w in workers)

    max_ns_per_obs = MVector{na, Int}()
    min_ns_per_obs = MVector{na, Int}()
    w_max_ns_per_obs = MVector{na, Int}()

    new_branches = Dict{Int, Nullable{Int}}()

    start_ns = time_ns()
    end_ns = start_ns + 1e9*ps.max_time
    for (w,c) in update_channels
        put!(c, GlobalUpdate(adjusted_ns, qs, Nullable{Int}(), false))
    end
    while true
        # process all waiting reports
        wait(report_channel)
        while isready(report_channel)
            r = take!(report_channel)
            if r.error
                error("error on worker")
            end
            worker_ns[r.id] = r.n_sims
            worker_n_obs[r.id] = r.n_obs
        end

        if time_ns() > end_ns
            break
        end

        # calculate global qs: each observation branch gets equal weight
        q_sums = zeros(Float64, na)
        n_obs = zeros(Int, na)
        fill!(max_ns_per_obs, typemin(Int))
        fill!(min_ns_per_obs, typemax(Int))
        for w in workers
            wno = worker_n_obs[w]
            n_obs += wno
            qs += worker_qs[w].*wno
            ns_per_obs = div.(worker_ns[w], wno)
            newmax = ns_per_obs .> max_ns_per_obs
            w_max_ns_per_obs[newmax] = w
            max_ns_per_obs[newmax] = ns_per_obs[newmax]
            min_ns_per_obs[:] = min.(ns_per_obs, min_ns_per_obs)
        end
        qs = convert(SVector{na, Float64}, q_sums./n_obs)

        # calculate adjusted global ns
        # find minimum ns per obs and act as if all obs get the same
        adjusted_ns = convert(SVector{na, Int}, min_ns_per_obs.*n_obs)

        if sum(adjusted_ns) >= ps.tree_queries
            break
        end

        # add new branches
        empty!(new_branches)
        for ai in 1:na
            if n_obs[ai] <= ps.k_observation*adjusted_ns[ai]^ps.alpha_observation
                # find worker with the maximum sims per obs for that action
                w = w_max_ns_per_obs[ai]
                new_branches[w] = ai
            end
        end

        # send instructions
        for (w, c) in update_channels
            put!(c, GlobalUpdate(adjusted_ns, qs, get(new_branches, w, Nullable{Int}()), false))
        end
    end

    info[:search_time_us] = (time_ns() - start_ns)/1000
    # terminate all worker tasks
    for (w, c) in update_channels
        put!(c, GlobalUpdate(adjusted_ns, qs, Nullable{Int}(), true))
    end
    info[:tree_queries] = sum(sum(ns for (w, ns) in worker_ns))
    info[:adjusted_queries] = sum(adjusted_ns)

    for f in futures
        fetch(f)
    end

    # choose the best
    return acts[indmax(qs)], info
end

function normalized_ns

end

end # module


text:
using Base.Test

using MPI

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

root = size-1
val = rank == root ? sum(0:size-1) : nothing
@test MPI.Reduce(rank, MPI.SUM, root, comm) == val

val = rank == root ? size-1 : nothing
@test MPI.Reduce(rank, MPI.MAX, root, comm) == val

val = rank == root ? 0 : nothing
@test MPI.Reduce(rank, MPI.MIN, root, comm) == val

val = rank == root ? size : nothing
@test MPI.Reduce(1, +, root, comm) == val

mesg = collect(1.0:5.0)
sum_mesg = MPI.Reduce(mesg, MPI.SUM, root, comm)
sum_mesg = rank == root ? sum_mesg : size*mesg
@test isapprox(norm(sum_mesg-size*mesg), 0.0)

MPI.Finalize()


text:
using SPN
using MNIST
using MultivariateStats
using Gadfly

x, y = traindata()
@show size(x)
@show size(y)

pcaModel = fit(PCA, x; maxoutdim=9)
xTransformed = transform(pcaModel, x)
# swap rows and columns:
x = xTransformed'

# Normalize the dataset such that each
# variable (column) has mean 0 and variance 1.

function normalizeDataset!(x::AbstractArray)
    x = x .- mean(x,1)
    x = x
    return x ./ std(x,1)
end

x = normalizeDataset!(x)
spn = structureLearnPoon(x,3,3,nsum=5,nleaf=10,baseres=1)
x = x[1:100,:]
llhvals = parameterLearnEM!(spn, x)

plot(y=llhvals, Geom.point, Geom.line)


text:
function add(obj::FieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "add", FieldMatrix, (FieldMatrix,), arg0)
end

function add_to_entry(obj::FieldMatrix, arg0::jint, arg1::jint, arg2::FieldElement)
    return jcall(obj, "addToEntry", void, (jint, jint, FieldElement), arg0, arg1, arg2)
end

function copy(obj::FieldMatrix)
    return jcall(obj, "copy", FieldMatrix, ())
end

function copy_sub_matrix(obj::FieldMatrix, arg0::Vector{jint}, arg1::Vector{jint}, arg2::Vector{Vector{FieldElement}})
    return jcall(obj, "copySubMatrix", void, (Vector{jint}, Vector{jint}, Vector{Vector{FieldElement}}), arg0, arg1, arg2)
end

function copy_sub_matrix(obj::FieldMatrix, arg0::jint, arg1::jint, arg2::jint, arg3::jint, arg4::Vector{Vector{FieldElement}})
    return jcall(obj, "copySubMatrix", void, (jint, jint, jint, jint, Vector{Vector{FieldElement}}), arg0, arg1, arg2, arg3, arg4)
end

function create_matrix(obj::FieldMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "createMatrix", FieldMatrix, (jint, jint), arg0, arg1)
end

function get_column(obj::FieldMatrix, arg0::jint)
    return jcall(obj, "getColumn", Vector{FieldElement}, (jint,), arg0)
end

function get_column_matrix(obj::FieldMatrix, arg0::jint)
    return jcall(obj, "getColumnMatrix", FieldMatrix, (jint,), arg0)
end

function get_column_vector(obj::FieldMatrix, arg0::jint)
    return jcall(obj, "getColumnVector", FieldVector, (jint,), arg0)
end

function get_data(obj::FieldMatrix)
    return jcall(obj, "getData", Vector{Vector{FieldElement}}, ())
end

function get_entry(obj::FieldMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "getEntry", FieldElement, (jint, jint), arg0, arg1)
end

function get_field(obj::FieldMatrix)
    return jcall(obj, "getField", Field, ())
end

function get_row(obj::FieldMatrix, arg0::jint)
    return jcall(obj, "getRow", Vector{FieldElement}, (jint,), arg0)
end

function get_row_matrix(obj::FieldMatrix, arg0::jint)
    return jcall(obj, "getRowMatrix", FieldMatrix, (jint,), arg0)
end

function get_row_vector(obj::FieldMatrix, arg0::jint)
    return jcall(obj, "getRowVector", FieldVector, (jint,), arg0)
end

function get_sub_matrix(obj::FieldMatrix, arg0::Vector{jint}, arg1::Vector{jint})
    return jcall(obj, "getSubMatrix", FieldMatrix, (Vector{jint}, Vector{jint}), arg0, arg1)
end

function get_sub_matrix(obj::FieldMatrix, arg0::jint, arg1::jint, arg2::jint, arg3::jint)
    return jcall(obj, "getSubMatrix", FieldMatrix, (jint, jint, jint, jint), arg0, arg1, arg2, arg3)
end

function get_trace(obj::FieldMatrix)
    return jcall(obj, "getTrace", FieldElement, ())
end

function multiply(obj::FieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "multiply", FieldMatrix, (FieldMatrix,), arg0)
end

function multiply_entry(obj::FieldMatrix, arg0::jint, arg1::jint, arg2::FieldElement)
    return jcall(obj, "multiplyEntry", void, (jint, jint, FieldElement), arg0, arg1, arg2)
end

function operate(obj::FieldMatrix, arg0::FieldVector)
    return jcall(obj, "operate", FieldVector, (FieldVector,), arg0)
end

function operate(obj::FieldMatrix, arg0::Vector{FieldElement})
    return jcall(obj, "operate", Vector{FieldElement}, (Vector{FieldElement},), arg0)
end

function power(obj::FieldMatrix, arg0::jint)
    return jcall(obj, "power", FieldMatrix, (jint,), arg0)
end

function pre_multiply(obj::FieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "preMultiply", FieldMatrix, (FieldMatrix,), arg0)
end

function pre_multiply(obj::FieldMatrix, arg0::FieldVector)
    return jcall(obj, "preMultiply", FieldVector, (FieldVector,), arg0)
end

function pre_multiply(obj::FieldMatrix, arg0::Vector{FieldElement})
    return jcall(obj, "preMultiply", Vector{FieldElement}, (Vector{FieldElement},), arg0)
end

function scalar_add(obj::FieldMatrix, arg0::FieldElement)
    return jcall(obj, "scalarAdd", FieldMatrix, (FieldElement,), arg0)
end

function scalar_multiply(obj::FieldMatrix, arg0::FieldElement)
    return jcall(obj, "scalarMultiply", FieldMatrix, (FieldElement,), arg0)
end

function set_column(obj::FieldMatrix, arg0::jint, arg1::Vector{FieldElement})
    return jcall(obj, "setColumn", void, (jint, Vector{FieldElement}), arg0, arg1)
end

function set_column_matrix(obj::FieldMatrix, arg0::jint, arg1::FieldMatrix)
    return jcall(obj, "setColumnMatrix", void, (jint, FieldMatrix), arg0, arg1)
end

function set_column_vector(obj::FieldMatrix, arg0::jint, arg1::FieldVector)
    return jcall(obj, "setColumnVector", void, (jint, FieldVector), arg0, arg1)
end

function set_entry(obj::FieldMatrix, arg0::jint, arg1::jint, arg2::FieldElement)
    return jcall(obj, "setEntry", void, (jint, jint, FieldElement), arg0, arg1, arg2)
end

function set_row(obj::FieldMatrix, arg0::jint, arg1::Vector{FieldElement})
    return jcall(obj, "setRow", void, (jint, Vector{FieldElement}), arg0, arg1)
end

function set_row_matrix(obj::FieldMatrix, arg0::jint, arg1::FieldMatrix)
    return jcall(obj, "setRowMatrix", void, (jint, FieldMatrix), arg0, arg1)
end

function set_row_vector(obj::FieldMatrix, arg0::jint, arg1::FieldVector)
    return jcall(obj, "setRowVector", void, (jint, FieldVector), arg0, arg1)
end

function set_sub_matrix(obj::FieldMatrix, arg0::Vector{Vector{FieldElement}}, arg1::jint, arg2::jint)
    return jcall(obj, "setSubMatrix", void, (Vector{Vector{FieldElement}}, jint, jint), arg0, arg1, arg2)
end

function subtract(obj::FieldMatrix, arg0::FieldMatrix)
    return jcall(obj, "subtract", FieldMatrix, (FieldMatrix,), arg0)
end

function transpose(obj::FieldMatrix)
    return jcall(obj, "transpose", FieldMatrix, ())
end

function walk_in_column_order(obj::FieldMatrix, arg0::FieldMatrixChangingVisitor)
    return jcall(obj, "walkInColumnOrder", FieldElement, (FieldMatrixChangingVisitor,), arg0)
end

function walk_in_column_order(obj::FieldMatrix, arg0::FieldMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInColumnOrder", FieldElement, (FieldMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_column_order(obj::FieldMatrix, arg0::FieldMatrixPreservingVisitor)
    return jcall(obj, "walkInColumnOrder", FieldElement, (FieldMatrixPreservingVisitor,), arg0)
end

function walk_in_column_order(obj::FieldMatrix, arg0::FieldMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInColumnOrder", FieldElement, (FieldMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_optimized_order(obj::FieldMatrix, arg0::FieldMatrixChangingVisitor)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldMatrixChangingVisitor,), arg0)
end

function walk_in_optimized_order(obj::FieldMatrix, arg0::FieldMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_optimized_order(obj::FieldMatrix, arg0::FieldMatrixPreservingVisitor)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldMatrixPreservingVisitor,), arg0)
end

function walk_in_optimized_order(obj::FieldMatrix, arg0::FieldMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInOptimizedOrder", FieldElement, (FieldMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_row_order(obj::FieldMatrix, arg0::FieldMatrixChangingVisitor)
    return jcall(obj, "walkInRowOrder", FieldElement, (FieldMatrixChangingVisitor,), arg0)
end

function walk_in_row_order(obj::FieldMatrix, arg0::FieldMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInRowOrder", FieldElement, (FieldMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_row_order(obj::FieldMatrix, arg0::FieldMatrixPreservingVisitor)
    return jcall(obj, "walkInRowOrder", FieldElement, (FieldMatrixPreservingVisitor,), arg0)
end

function walk_in_row_order(obj::FieldMatrix, arg0::FieldMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInRowOrder", FieldElement, (FieldMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end


text:
abstract type AbstractGBGraph{T} <: Graphs.AbstractGraph{UInt64} end

Base.eltype(::AbstractGBGraph) = UInt64

mutable struct GBWeightedGraph{T} <: AbstractGBGraph{T}
    A::GBMatrix{T} # The adjacency matrix of the graph
    nvertices::Int64
    function GBWeightedGraph{T}(adjmatrix::GBMatrix{T}, nvertices=size(adjmatrix, 1)) where {T}
        drows, dcols = size(adjmatrix)
        isequal(drows, dcols) || error("Adjacency matrix must be square.")
        issymmetric(adjmatrix) || error("Adjacency matrix must be symmetric.")
        gbset(adjmatrix, :format, :byrow) # we want graphs stored by row for the most part.
        return new{T}(adjmatrix, nvertices)
    end
end

weighttype(::GBWeightedGraph{T}) where T = T

Graphs.ne(g::GBWeightedGraph) = nnz(g.A)
Graphs.nv(g::GBWeightedGraph) = size(g.A, 1)

vertices(g::GBWeightedGraph) = one(UInt64):UInt64(nv(g))

Graphs.edgetype(::GBWeightedGraph{T}) where T = GBWeightedEdge{T}
Graphs.edges(g::GBWeightedGraph) = (GBWeightedEdge(x...) for x in zip(findnz(triu(g.A))...))

Graphs.weights(g::GBWeightedGraph) = g.A

Graphs.inneighbors(g::GBWeightedGraph, x...) = outneighbors(g, x...)

function Graphs.add_edge!(g::GBWeightedGraph, e::GBWeightedEdge)
    s, d, weight = Tuple(e)
    @inbounds g.A[s, d] = weight
    @inbounds g.A[d, s] = weight
    return true
end

function Graphs.rem_edge!(g::GBWeightedGraph, e::GBWeightedEdge)
    deleteat!(g.A, e.dst, e.src)
    deleteat!(g.A, e.src, e.dst)
    return true
end

Base.:(==)(g::GBWeightedGraph, h::GBWeightedGraph) = g.A == h.A

Graphs.is_directed(::Type{<:GBWeightedGraph}) = false

Graphs.has_edge(g::GBWeightedGraph, e::GBWeightedEdge) = g.A[dst(e), src(e)] !== nothing

# handles single-argument edge constructors such as pairs and tuples
Graphs.has_edge(g::GBWeightedGraph{T}, x) where T = has_edge(g, edgetype(g)(x))
Graphs.add_edge!(g::GBWeightedGraph{T}, x) where T = add_edge!(g, edgetype(g)(x))

# handles two-argument edge constructors like src,dst
Graphs.has_edge(g::GBWeightedGraph, x, y) = has_edge(g, edgetype(g)(x, y, 0))
Graphs.add_edge!(g::GBWeightedGraph, x, y) = add_edge!(g, edgetype(g)(x, y, 1))
Graphs.add_edge!(g::GBWeightedGraph, x, y, z) = add_edge!(g, edgetype(g)(x, y, z))
Graphs.has_vertex(g::GBWeightedGraph, v::Integer) = v in vertices(g)

function Graphs.rem_edge!(g::GBWeightedGraph{T}, u::Integer, v::Integer) where {T}
    rem_edge!(g, edgetype(g)(T(u), T(v), one(T)))
end

# IMPLEMENT WITH PACK/UNPACK.
function Graphs.outneighbors(g::GBWeightedGraph, v::Integer)
    row = findnz(g.A[v, :])[2]


text:
"""
    PearsonResiduals(model, data)
Calculates Pearson residuals between model predicted values and the actual response variable values.
Usage:
```julia-repl
PearsonResiduals(model, data)
```
Arguments:
- `model` : The `GLM` model.
- `data` : The `DataFrame` containing data that was used as input to the model.
"""

function PearsonResiduals(model::StatsModels.TableRegressionModel, data::DataFrame)

    # Generate predictions

    f_hat = predict(model)

    # Parse response vector

    y_name = string(formula(model).lhs)
    y = Array(data[!, names(data, y_name)])

    # Calculate residuals as per https://www.datascienceblog.net/post/machine-learning/interpreting_generalized_linear_models/

    r = (y .- f_hat) ./ .√(f_hat)
    r = sum(r .^ 2)
    return r
end

"""
    EstimateDispersionParameter(residuals, model)
Estimates the dispersion parameter ϕ by standardising the sum of squared Pearson residuals against the residual degrees of freedom.
Usage:
```julia-repl
EstimateDispersionParameter(residuals, model)
```
Arguments:
- `residuals` : The sum of squared Pearson residuals.
- `model` : The `GLM` model.
"""

function EstimateDispersionParameter(residuals::Float64, model::StatsModels.TableRegressionModel)

    # Calculate dispersion/scale parameter estimate by dividing Pearson residuals by residual degrees of freedom

    ϕ = residuals / (dof_residual(model) - 1) # This aligns calculation with R's df.residuals function
    println("\nDispersion parameter (ϕ) for model taken to be " * string(round(ϕ, digits = 5)))
    println("Standard errors are multiplied by " * string(round(sqrt(ϕ), digits = 5)) * " to adjust for dispersion parameter (ϕ)")
    return ϕ
end

"""
    coefarray(model, ϕ; level)
Calculates relevant statistics for inference based of estimates and dispersion-adjusted standard errors and returns results in a concatenated array.
Usage:
```julia-repl
coefarray(model, ϕ; level)
```
Arguments:
- `model` : The `GLM` model.
- `ϕ` : The estimated dispersion parameter.
- `level` : The desired degree of confidence.
"""

function coefarray(model::StatsModels.TableRegressionModel, ϕ::Real; level::Real=0.95)

    # NOTE: Function modified from https://docs.juliahub.com/AxisIndices/AHOcZ/0.6.3/coeftable/

    cc = coef(model)
    se = stderror(model) * ϕ
    tt = cc ./ se
    p = ccdf.(Ref(FDist(1, dof_residual(model))), abs2.(tt))
    ci = se * quantile(TDist(dof_residual(model)), (1 - level) / 2)
    ct = hcat(coefnames(model), cc, se, tt, p, cc + ci, cc - ci)
    return ct
end


text:
module operators_sparse

import Base: ==, *, /, +, -
import ..operators

using ..bases, ..states, ..operators, ..sparsematrix

export SparseOperator

"""
Sparse array implementation of Operator.

The matrix is stored as the julia built-in type SparseMatrixCSC
in the data field.
"""
type SparseOperator <: Operator
    basis_l::Basis
    basis_r::Basis
    data::SparseMatrixCSC{Complex128}
    function SparseOperator(b1::Basis, b2::Basis, data)
        if length(b1) != size(data, 1) || length(b2) != size(data, 2)
            throw(DimensionMismatch())
        end
        new(b1, b2, data)
    end
end

SparseOperator(b::Basis, data::SparseMatrixCSC{Complex128}) = SparseOperator(b, b, data)
SparseOperator(b::Basis, data::Matrix{Complex128}) = SparseOperator(b, sparse(data))
SparseOperator(op::DenseOperator) = SparseOperator(op.basis_l, op.basis_r, sparse(op.data))

SparseOperator(b1::Basis, b2::Basis) = SparseOperator(b1, b2, spzeros(Complex128, length(b1), length(b2)))
SparseOperator(b::Basis) = SparseOperator(b, b)

Base.full(a::SparseOperator) = DenseOperator(a.basis_l, a.basis_r, full(a.data))

Base.sparse(a::Operator) = throw(ArgumentError("Direct conversion from $(typeof(a)) not implemented. Use sparse(full(op)) instead."))
Base.sparse(a::SparseOperator) = deepcopy(a)
Base.sparse(a::DenseOperator) = SparseOperator(a.basis_l, a.basis_r, sparse(a.data))

==(x::SparseOperator, y::SparseOperator) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && (x.data == y.data)

*(a::SparseOperator, b::SparseOperator) = (check_multiplicable(a.basis_r, b.basis_l); SparseOperator(a.basis_l, b.basis_r, a.data*b.data))
*(a::SparseOperator, b::Number) = SparseOperator(a.basis_l, a.basis_r, complex(b)*a.data)
*(a::Number, b::SparseOperator) = SparseOperator(b.basis_l, b.basis_r, complex(a)*b.data)

/(a::SparseOperator, b::Number) = SparseOperator(a.basis_l, a.basis_r, a.data/complex(b))

+(a::SparseOperator, b::SparseOperator) = (operators.check_samebases(a,b); SparseOperator(a.basis_l, a.basis_r, a.data+b.data))

-(a::SparseOperator) = SparseOperator(a.basis_l, a.basis_r, -a.data)
-(a::SparseOperator, b::SparseOperator) = (operators.check_samebases(a,b); SparseOperator(a.basis_l, a.basis_r, a.data-b.data))

# Fast in-place multiplication implementations
operators.gemm!{T<:Complex}(alpha::T, M::SparseOperator, b::DenseOperator, beta::T, result::DenseOperator) = sparsematrix.gemm!(alpha, M.data, b.data, beta, result.data)
operators.gemm!{T<:Complex}(alpha::T, a::DenseOperator, M::SparseOperator, beta::T, result::DenseOperator) = sparsematrix.gemm!(alpha, a.data, M.data, beta, result.data)
operators.gemv!{T<:Complex}(alpha::T, M::SparseOperator, b::Ket, beta::T, result::Ket) = sparsematrix.gemv!(alpha, M.data, b.data, beta, result.data)
operators.gemv!{T<:Complex}(alpha::T, b::Bra, M::SparseOperator, beta::T, result::Bra) = sparsematrix.gemv!(alpha, b.data, M.data, beta, result.data)

operators.tensor(a::SparseOperator, b::SparseOperator) = SparseOperator(tensor(a.basis_l, b.basis_l), tensor(a.basis_r, b.basis_r), kron(a.data, b.data))

operators.dagger(x::SparseOperator) = SparseOperator(x.basis_r, x.basis_l, x.data')

operators.norm(op::SparseOperator, p) = norm(op.data, p)
operators.trace(op::SparseOperator) = trace(op.data)

Base.identity(b::Basis) = SparseOperator(b, b, speye(Complex128, length(b)))
Base.identity(b1::Basis, b2::Basis) = SparseOperator(b1, b2, speye(Complex128, length(b1), length(b2)))
Base.identity(op::SparseOperator) = identity(op.basis_l, op.basis_r)

function operators.embed(basis::CompositeBasis, indices::Vector{Int}, operators::Vector{SparseOperator})
    @assert length(indices) == length(operators)
    @assert length(basis.bases) >= maximum(indices)
    if length(basis.bases) == 1
        return SparseOperator(basis, basis, deepcopy(operators[1].data))
    end
    op_total = (1 in indices ? operators[findfirst(indices, 1)] : identity(basis.bases[1]))
    for i=2:length(basis.bases)
        op = (i in indices ? operators[findfirst(indices, i)] : identity(basis.bases[i]))
        op_total = tensor(op_total, op)
    end
    return op_total
end

end # module


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

JLLWrappers.@generate_main_file_header("snappy")
JLLWrappers.@generate_main_file("snappy", UUID("fe1e1685-f7be-5f59-ac9f-4ca204017dfd"))
end  # module snappy_jll


text:
# SimpleProjections.jl
#
# A small module for solving a distributionally robust optimization on
# the chi-square divergence ball on the simplex. The bisection
# algorithm takes a desired (relative) solution tolerance as a
# parameter, and is extremely quick for reasonable values of n and
# tolerance.
#
# Uses careful book-keeping and a number of sorted vectors to actually
# perform the projections, which require time
#
# O(n * log(n) + (log(n) * log(1 / epsilon)))
#
# where the first term is a sort, the second is a careful binary search.

module SProj

function RobustExpectation(z::Vector{Float64}, rho::Float64;
                           sol_tolerance::Float64 = 1e-13)
  return dot(z, MinimizeLinearSquareDivergence(z, rho; sol_tolerance = sol_tolerance));
end

# x = MinimizeLinearSquareDivergence(z, rho)
#
# Sets x to be the minimizer of
#
# minimize  sum(x .* z)
#   s.t.    (1/2n) * sum((x * n - 1)^2) <= rho
#           x >= 0, sum(x) == 1.
#
# Does this by performing a binary search on the partially dualized form
#
# maximize  inf_{x >= 0, sum(x) == 1}
#             lambda * ((1/2) * norm(x - 1/n)^2 - rho / n) + x' * z,
#
# by performing a binary search over lambda. With appropriate
# book-keeping, each successive step of the binary search requires
# time O(log(n)).
function MinimizeLinearSquareDivergence(z::Vector{Float64}, rho::Float64;
                                        sol_tolerance::Float64 = 1e-13)
  z_sort = sort(z - mean(z)); # just normalization
  nn = length(z_sort);
  if all(z_sort .== 0)
    return (1/nn) * ones(nn);
  end
  lambda_min = 0;
  lambda_init_max = max(maximum(abs(z_sort)) * nn,
                        sqrt(nn / (2 * rho)) * norm(z_sort));
  lambda_max = lambda_init_max;
  sum_vec = cumsum(z_sort);
  sum_squared_vec = cumsum(z_sort.^2);
  z_norm_squared = sum(z_sort.^2);

  while (abs(lambda_max - lambda_min) > sol_tolerance * lambda_init_max)
    lambda = (lambda_max + lambda_min) / 2;
    (eta, ind) = FindShiftVectorOntoSimplex(z_sort, lambda, sum_vec);

    # Now compute the derivative w.r.t. lambda, which is
    #
    # norm(x - (1/n))^2 - rho / n.
    #
    # We can do this in closed form using our book-keeping
    square_sum = (sum_squared_vec[ind] / lambda^2 + ind * eta^2
                  + 2 * sum_vec[ind] * eta / lambda + (nn - ind) / nn^2);
    if (square_sum / 2 - rho / nn > 0)
      lambda_min = lambda;
    else
      lambda_max = lambda;
    end
  end
  lambda = (lambda_max + lambda_min) / 2;
  (eta, ind) = FindShiftVectorOntoSimplex(z_sort, lambda, sum_vec);
  x = max(ones(nn) / nn - (z - mean(z)) / lambda - eta, 0);
  return x;
end

# (eta, ind) = FindShiftVectorOntoSimplex(z::Vector{Float64}, lambda::Float64,
#                                         s::Vector{Float64})
#
# Let v = (1/n) - z / lambda, where z and s are n-vectors and
# lambda > 0.  Finds the eta such that if
#
#   x_i = (1/n - z / lambda - eta)_+
#
# then sum_i x_i = 1. The vector z must be in ascending sorted order,
# so that z[1] <= z[2] <= ... <= z[n]. The vector s is the cumulative
# sum of the z vector, so that s[i] = sum_{j = 1}^i z[j].
#
# This is equivalent to finding the projection of the vector
#
#   v = 1/n - z / lambda
#
# onto the probability simplex. To do so, we find the index i such that
#
# sum_{j = 1}^i (v[j] - v[i]) < 1
#
# and
#
# sum_{j = 1}^{i + 1} (v[j] - v[i + 1]) >= 1.
#
# This method is explained in Figure 1 of Duchi et al., Efficient
# Projections onto the L1-Ball for Learning in High Dimensions, ICML
# 2008.
#
# Method runs in O(log n) time, and returns the pair (eta, ind) where
# ind satisfies the two preceding inequalities when i = ind.
function FindShiftVectorOntoSimplex(z::Vector{Float64}, lambda::Float64,
                                    sumvec::Vector{Float64})
  nn = length(z);
  low_ind = 1;
  high_ind = nn;
  if (((1/nn) - z[nn] / lambda + (1/nn) * sumvec[end] / lambda) > 0)
    return (-(1/nn) * sumvec[end] / lambda, nn);
  end

  while (low_ind != high_ind)
    ii = floor(Int64, (high_ind + low_ind) / 2);
    left_sum = (ii * z[ii] - sumvec[ii]) / lambda;
    right_sum = ((ii + 1) * z[ii + 1] - sumvec[ii + 1]) / lambda;
    if (right_sum >= 1 && left_sum < 1)
      eta = (1 / nn - 1 / ii) - (1 / ii) * sumvec[ii] / lambda;
      return (eta, ii);
    elseif (left_sum >= 1)
      # Need to decrease index, left sum is too large
      high_ind = ii - 1;
    else
      # Increase index, because right sum is too low
      low_ind = ii + 1;
    end
  end
  ii = low_ind;
  eta = (1 / nn - 1 / ii) - (1 / ii) * sumvec[ii] / lambda;
  return (eta, ii);
end

# TestProjections()
#
# A unit test to verify that a number of simple special cases of
# projections work.
function TestProjections()
  test_tol = 1e-6;
  rho = 1.2;
  z = repmat(linspace(0, 1, 4), 4, 1)[:];
  x = MinimizeLinearSquareDivergence(z, rho);
  num_passed = 0;
  total_to_pass = 0;
  x_desired = [ 0.229583
               0.0204175
               0.0
               0.0
               0.229583
               0.0204175
               0.0
               0.0
               0.229583
               0.0204175
               0.0
               0.0
               0.229583
               0.0204175
               0.0
               0.0  ];
  total_to_pass += 1;
  if (norm(x - x_desired) > test_tol)
    warn(string("Projection of 0, .33, .66, 1 vector failed."));
  else
    num_passed += 1;
  end
  z =  [-0.705131
        -1.6036
        1.6838
        -0.167923
        -0.514595
        0.488477 ];
  x = MinimizeLinearSquareDivergence(z, rho);
  x_desired = [  0.194745
               0.722411
               0.0
               0.0
               0.0828439
               0.0
               ];
  total_to_pass += 1;
  if (norm(x - x_desired) > test_tol)
    warn(string("Random-ish string of z failed."));
  else
    num_passed += 1;
  end
  x = MinimizeLinearSquareDivergence(z, 200.0);
  total_to_pass += 1;
  if (norm(x - [0; 1; 0; 0; 0; 0]) > test_tol)
    warn(string("Large rho test failed."));
  else
    num_passed += 1;
  end
  x = MinimizeLinearSquareDivergence(z, 1e-14);
  total_to_pass += 1;
  if (norm(x - ones(6) / 6) > test_tol)
    warn(string("Very small rho test failed."));
  else
    num_passed += 1;
  end
  println("Passed ", num_passed, " of ", total_to_pass, " tests.");
end

end  # module SProj

max_stars_repo_name: syadlowsky/robustopt

max_issues_repo_name: syadlowsky/robustopt

max_forks_repo_name: syadlowsky/robustopt


text:
export splitting_field, issubfield, isdefining_polynomial_nice,
       quadratic_field, islinearly_disjoint, rationals_as_number_field

################################################################################
#
#  Base field
#
################################################################################

base_field(K::AnticNumberField) = FlintQQ

################################################################################
#
#  Order type
#
################################################################################

order_type(::AnticNumberField) = NfAbsOrd{AnticNumberField, nf_elem}

order_type(::Type{AnticNumberField}) = NfAbsOrd{AnticNumberField, nf_elem}

################################################################################
#
#  Predicates
#
################################################################################

issimple(::Type{AnticNumberField}) = true

issimple(::AnticNumberField) = true

################################################################################
#
#  Field constructions
#
################################################################################

@doc Markdown.doc"""
    NumberField(S::Generic.ResRing{fmpq_poly}; cached::Bool = true, check::Bool = true) -> AnticNumberField, Map

 The number field $K$ isomorphic to the ring $S$ and the map from $K\to S$.
"""
function NumberField(S::Generic.ResRing{fmpq_poly}; cached::Bool = true, check::Bool = true)
  Qx = parent(modulus(S))
  K, a = NumberField(modulus(S), "_a", cached = cached, check = check)
  mp = MapFromFunc(y -> S(Qx(y)), x -> K(lift(x)), K, S)
  return K, mp
end

function NumberField(f::fmpq_poly; cached::Bool = true, check::Bool = true)
  return NumberField(f, "_a", cached = cached, check = check)
end

function NumberField(f::fmpz_poly, s::Symbol; cached::Bool = true, check::Bool = true)
  Qx = Globals.Qx
  return NumberField(Qx(f), String(s), cached = cached, check = check)
end

function NumberField(f::fmpz_poly, s::AbstractString; cached::Bool = true, check::Bool = true)
  Qx = Globals.Qx
  return NumberField(Qx(f), s, cached = cached, check = check)
end

function NumberField(f::fmpz_poly; cached::Bool = true, check::Bool = true)
  Qx = Globals.Qx
  return NumberField(Qx(f), cached = cached, check = check)
end

function radical_extension(n::Int, gen::Integer; cached::Bool = true, check::Bool = true)
  return radical_extension(n, fmpz(gen), cached = cached, check = check)
end

function radical_extension(n::Int, gen::fmpz; cached::Bool = true, check::Bool = true)
  x = gen(Globals.Qx)
  return number_field(x^n - gen, cached = cached, check = check)
end

@doc doc"""
    cyclotomic_field(n::Int) -> AnticNumberField, nf_elem

The cyclotomic field defined by the $n$-th cyclotomic polynomial.

# Examples

```jldoctest
julia> cyclotomic_field(10)
(Cyclotomic field of order 10, z_10)
```
"""
function cyclotomic_field(n::Int; cached::Bool = true)
  return CyclotomicField(n, "z_$n", cached = cached)
end

# TODO: Some sort of reference?
@doc doc"""
    wildanger_field(n::Int, B::fmpz) -> AnticNumberField, nf_elem

Returns the field with defining polynomial $x^n + \sum_{i=0}^{n-1} (-1)^{n-i}Bx^i$.
These fields tend to have non-trivial class groups.

# Examples

```jldoctest
julia> wildanger_field(3, ZZ(10), "a")
(Number field over Rational Field with defining polynomial x^3 - 10*x^2 + 10*x - 10, a)
```
"""
function wildanger_field(n::Int, B::fmpz, s::String = "_\$"; check::Bool = true, cached::Bool = true)
  x = gen(Globals.Qx)
  f = x^n
  for i=0:n-1
    f += (-1)^(n-i)*B*x^i
  end
  return NumberField(f, s, cached = cached, check = check)
end

function wildanger_field(n::Int, B::Integer, s::String = "_\$"; cached::Bool = true, check::Bool = true)
  return wildanger_field(n, fmpz(B), s, cached = cached, check = check)
end

@doc Markdown.doc"""
    quadratic_field(d::IntegerUnion) -> AnticNumberField, nf_elem

Returns the field with defining polynomial $x^2 - d$.

# Examples

```jldoctest
julia> quadratic_field(5)
(Real quadratic field defined by x^2 - 5, sqrt(5))
```
"""
function quadratic_field(d::IntegerUnion; cached::Bool = true, check::Bool = true)
end

function quadratic_field(d::fmpz; cached::Bool = true, check::Bool = true)
  x = gen(Globals.Qx)
  if nbits(d) > 100
    a = div(d, fmpz(10)^(ndigits(d, 10) - 4))
    b = mod(abs(d), 10^4)
    s = "sqrt($a..($(nbits(d)) bits)..$b)"
  else
    s = "sqrt($d)"
  end
  q, a = number_field(x^2-d, s, cached = cached, check = check)
  set_special(q, :show => show_quad)
  return q, a
end

function show_quad(io::IO, q::AnticNumberField)
  d = trailing_coefficient(q.pol)
  if d < 0
    print(io, "Real quadratic field defined by ", q.pol)
  else
    print(io, "Imaginary quadratic field defined by ", q.pol)
  end
end

function quadratic_field(d::Integer; cached::Bool = true, check::Bool = true)
  return quadratic_field(fmpz(d), cached = cached, check = check)
end

@doc doc"""
    rationals_as_number_field() -> AnticNumberField, nf_elem

Returns the rational numbers as the number field defined by $x - 1$.

# Examples

```jldoctest
julia> rationals_as_number_field()
(Number field over Rational Field with defining polynomial x - 1, 1)
```
"""
function rationals_as_number_field()
  x = gen(Globals.Qx)
  return number_field(x-1)
end

################################################################################
#
#  Predicates
#
################################################################################

@doc Markdown.doc"""
    isdefining_polynomial_nice(K::AnticNumberField)

Tests if the defining polynomial of $K$ is integral and monic.
"""
function isdefining_polynomial_nice(K::AnticNumberField)
  return Bool(K.flag & UInt(1))
end

function isdefining_polynomial_nice(K::NfAbsNS)
  pols = K.pol
  for i = 1:length(pols)
    d = denominator(pols[i])
    if !isone(d)
      return false
    end
    if !isone(leading_coefficient(pols[i]))
      return false
    end
  end
  return true
end

################################################################################
#
#  Class group
#
################################################################################

@doc Markdown.doc"""
    class_group(K::AnticNumberField) -> GrpAbFinGen, Map

Shortcut for `class_group(maximal_order(K))`: returns the class
group as an abelian group and a map from this group to the set
of ideals of the maximal order.
"""
function class_group(K::AnticNumberField)
  return class_group(maximal_order(K))
end

################################################################################
#
#  Class number
#
################################################################################

@doc Markdown.doc"""
    class_number(K::AnticNumberField) -> fmpz

Returns the class number of $K$.
"""
function class_number(K::AnticNumberField)
  return order(class_group(maximal_order(K))[1])
end

################################################################################
#
#  Relative class number
#
################################################################################

@doc Markdown.doc"""
    relative_class_number(K::AnticNumberField) -> fmpz

Returns the relative class number of $K$. The field must be a CM-field.
"""
function relative_class_number(K::AnticNumberField)
  if degree(K) == 2
    @req istotally_complex(K) "Field must be a CM-field"
    return class_number(K)
  end

  fl, c = iscm_field(K)
  @req fl "Field must be a CM-field"
  h = class_number(K)
  L, _ = fixed_field(K, c)
  hp = class_number(L)
  @assert mod(h, hp) == 0
  return divexact(h, hp)
end

################################################################################
#
#  Basis
#
################################################################################

function basis(K::AnticNumberField)
  n = degree(K)
  g = gen(K);
  d = Array{typeof(g)}(undef, n)
  b = K(1)
  for i = 1:n-1
    d[i] = b
    b *= g
  end
  d[n] = b
  return d
end

################################################################################
#
#  Torsion units and related functions
#
################################################################################

@doc Markdown.doc"""
    istorsion_unit(x::nf_elem, checkisunit::Bool = false) -> Bool

Returns whether $x$ is a torsion unit, that is, whether there exists $n$ such
that $x^n = 1$.

If `checkisunit` is `true`, it is first checked whether $x$ is a unit of the
maximal order of the number field $x$ is lying in.
"""
function istorsion_unit(x::nf_elem, checkisunit::Bool = false)
  if checkisunit
    _isunit(x) ? nothing : return false
  end

  K = parent(x)
  d = degree(K)
  c = conjugate_data_arb(K)
  r, s = signature(K)

  while true
    @vprint :UnitGroup 2 "Precision is now $(c.prec) \n"
    l = 0
    @vprint :UnitGroup 2 "Computing conjugates ... \n"
    cx = conjugates_arb(x, c.prec)
    A = ArbField(c.prec, cached = false)
    for i in 1:r
      k = abs(cx[i])
      if k > A(1)
        return false
      elseif isnonnegative(A(1) + A(1)//A(6) * log(A(d))//A(d^2) - k)
        l = l + 1
      end
    end
    for i in 1:s
      k = abs(cx[r + i])
      if k > A(1)
        return false
      elseif isnonnegative(A(1) + A(1)//A(6) * log(A(d))//A(d^2) - k)
        l = l + 1
      end
    end

    if l == r + s
      return true
    end
    refine(c)
  end
end

@doc Markdown.doc"""
    torsion_unit_order(x::nf_elem, n::Int)

Given a torsion unit $x$ together with a multiple $n$ of its order, compute
the order of $x$, that is, the smallest $k \in \mathbb Z_{\geq 1}$ such
that $x^`k` = 1$.

It is not checked whether $x$ is a torsion unit.
"""
function torsion_unit_order(x::nf_elem, n::Int)
  ord = 1
  fac = factor(n)
  for (p, v) in fac
    p1 = Int(p)
    s = x^divexact(n, p1^v)
    if isone(s)
      continue
    end
    cnt = 0
    while !isone(s) && cnt < v+1
      s = s^p1
      ord *= p1
      cnt += 1
    end
    if cnt > v+1
      error("The element is not a torsion unit")
    end
  end
  return ord
end

#################################################################################################
#
#  Normal Basis
#
#################################################################################################

function normal_basis(K::AnticNumberField)
  # First try basis elements of LLL basis
  # or rather not
  # n = degree(K)
  # Aut = automorphisms(K)

  # length(Aut) != n && error("The field is not normal over the rationals!")

  # A = zero_matrix(FlintQQ, n, n)
  # _B = basis(lll(maximal_order(K)))
  # for i in 1:n
  #   r = elem_in_nf(_B[i])
  #   for i = 1:n
  #     y = Aut[i](r)
  #     for j = 1:n
  #       A[i,j] = coeff(y, j - 1)
  #     end
  #   end
  #   if rank(A) == n
  #     return r
  #   end
  # end

  O = EquationOrder(K)
  Qx = parent(K.pol)
  d = discriminant(O)
  p = 1
  for q in PrimesSet(degree(K), -1)
    if divisible(d, q)
      continue
    end
    #Now, I check if p is totally split
    R = GF(q, cached = false)
    Rt, t = PolynomialRing(R, "t", cached = false)
    ft = Rt(K.pol)
    pt = powermod(t, q, ft)
    if degree(gcd(ft, pt-t)) == degree(ft)
      p = q
      break
    end
  end

  return _normal_basis_generator(K, p)
end

function _normal_basis_generator(K, p)
  Qx = parent(K.pol)

  #Now, I only need to lift an idempotent of O/pO
  R = GF(p, cached = false)
  Rx, x = PolynomialRing(R, "x", cached = false)
  f = Rx(K.pol)
  fac = factor(f)
  g = divexact(f, first(keys(fac.fac)))
  Zy, y = PolynomialRing(FlintZZ, "y", cached = false)
  g1 = lift(Zy, g)
  return K(g1)
end

################################################################################
#
#  Subfield check
#
################################################################################

function _issubfield(K::AnticNumberField, L::AnticNumberField)
  f = K.pol
  R = roots(f, L, max_roots = 1)
  if isempty(R)
    return false, L()
  else
    h = parent(L.pol)(R[1])
    return true, h(gen(L))
  end
end

function _issubfield_first_checks(K::AnticNumberField, L::AnticNumberField)
  f = K.pol
  g = L.pol
  if mod(degree(g), degree(f)) != 0
    return false
  end
  t = divexact(degree(g), degree(f))
  try
    OK = _get_maximal_order_of_nf(K)
    OL = _get_maximal_order_of_nf(L)
    if mod(discriminant(OL), discriminant(OK)^t) != 0
      return false
    end
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
  end
  # We could factorize the discriminant of f, but we only test small primes.
  cnt_threshold = 10*degree(K)
  p = 3
  cnt = 0
  while cnt < cnt_threshold
    F = GF(p, cached = false)
    Fx = PolynomialRing(F, "x", cached = false)[1]
    fp = Fx(f)
    gp = Fx(g)
    if !issquarefree(fp) || !issquarefree(gp)
      p = next_prime(p)
	    continue
    end
    cnt += 1
    fs = factor_shape(fp)
    gs = factor_shape(gp)
    if !divisible(lcm(collect(keys(gs))), lcm(collect(keys(fs))))
      return false
    end
    p = next_prime(p)
  end
  return true
end

function issubfield(K::AnticNumberField, L::AnticNumberField)
  fl = _issubfield_first_checks(K, L)
  if !fl
    return false, hom(K, L, zero(L), check = false)
  end
  b, prim_img = _issubfield(K, L)
  return b, hom(K, L, prim_img, check = false)
end

function _issubfield_normal(K::AnticNumberField, L::AnticNumberField)
  f = K.pol
  f1 = change_base_ring(L, f)
  r = roots(f1, max_roots = 1, isnormal = true)
  if length(r) > 0
    h = parent(L.pol)(r[1])
    return true, h(gen(L))
  else
    return false, L()
  end
end

@doc Markdown.doc"""
      issubfield_normal(K::AnticNumberField, L::AnticNumberField) -> Bool, NfToNfMor

Returns `true` and an injection from $K$ to $L$ if $K$ is a subfield of $L$.
Otherwise the function returns "false" and a morphism mapping everything to 0.

This function assumes that $K$ is normal.
"""
function issubfield_normal(K::AnticNumberField, L::AnticNumberField)
  fl = _issubfield_first_checks(K, L)
  if !fl
    return false, hom(K, L, zero(L), check = false)
  end
  b, prim_img = _issubfield_normal(K, L)
  return b, hom(K, L, prim_img, check = false)

end

################################################################################
#
#  Isomorphism
#
################################################################################

@doc Markdown.doc"""
    isisomorphic(K::AnticNumberField, L::AnticNumberField) -> Bool, NfToNfMor

Returns "true" and an isomorphism from $K$ to $L$ if $K$ and $L$ are isomorphic.
Otherwise the function returns "false" and a morphism mapping everything to 0.
"""
function isisomorphic(K::AnticNumberField, L::AnticNumberField)
  f = K.pol
  g = L.pol
  if degree(f) != degree(g)
    return false, hom(K, L, zero(L), check = false)
  end
  if fmpq[coeff(f, i) for i = 0:degree(f)] == fmpq[coeff(g, i) for i = 0:degree(g)]
    return true, hom(K, L, gen(L))
  end
  if signature(K) != signature(L)
    return false, hom(K, L, zero(L), check = false)
  end
  try
    OK = _get_maximal_order_of_nf(K)
    OL = _get_maximal_order_of_nf(L)
    if discriminant(OK) != discriminant(OL)
      return false, hom(K, L, zero(L), check = false)
    end
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
    t = discriminant(f)//discriminant(g)
    if !issquare(numerator(t)) || !issquare(denominator(t))
      return false, hom(K, L, zero(L), check = false)
    end
  end
  p = 10^5
  cnt = 0
  df = denominator(f)
  dg = denominator(g)
  while cnt < max(20, 2*degree(K))
    p = next_prime(p)
    if divisible(df, p) || divisible(dg, p)
      continue
    end
    F = GF(p, cached = false)
    Fx = PolynomialRing(F, "x", cached = false)[1]
    fp = Fx(f)
    if degree(fp) != degree(f) || !issquarefree(fp)
      continue
    end
    gp = Fx(g)
    if degree(gp) != degree(g) || !issquarefree(gp)
      continue
    end
    cnt += 1
    lf = factor_shape(fp)
    lg = factor_shape(gp)
    if lf != lg
      return false, hom(K, L, zero(L), check = false)
    end
  end
  b, prim_img = _issubfield(K, L)
  if !b
    return b, hom(K, L, zero(L), check = false)
  else
    return b, hom(K, L, prim_img, check = false)
  end
end

################################################################################
#
#  Compositum
#
################################################################################

@doc Markdown.doc"""
    compositum(K::AnticNumberField, L::AnticNumberField) -> AnticNumberField, Map, Map

Assuming $L$ is normal (which is not checked), compute the compositum $C$ of the
2 fields together with the embedding of $K \to C$ and $L \to C$.
"""
function compositum(K::AnticNumberField, L::AnticNumberField)
  lf = factor(K.pol, L)
  d = degree(first(lf.fac)[1])
  if any(x->degree(x) != d, keys(lf.fac))
    error("2nd field cannot be normal")
  end
  KK = NumberField(first(lf.fac)[1])[1]
  Ka, mKa = absolute_simple_field(KK)
  mK = hom(K, Ka, mKa\gen(KK))
  mL = hom(L, Ka, mKa\(KK(gen(L))))
  embed(mK)
  embed(mL)
  return Ka, mK, mL
end

################################################################################
#
#  Serialization
#
################################################################################

# This function can be improved by directly accessing the numerator
# of the fmpq_poly representing the nf_elem
@doc Markdown.doc"""
    write(io::IO, A::Vector{nf_elem}) -> Nothing

Writes the elements of `A` to `io`. The first line are the coefficients of
the defining polynomial of the ambient number field. The following lines
contain the coefficients of the elements of `A` with respect to the power
basis of the ambient number field.
"""
function write(io::IO, A::Vector{nf_elem})
  if length(A) == 0
    return
  else
    # print some useful(?) information
    print(io, "# File created by Hecke $VERSION_NUMBER, $(Base.Dates.now()), by function 'write'\n")
    K = parent(A[1])
    polring = parent(K.pol)

    # print the defining polynomial
    g = K.pol
    d = denominator(g)

    for j in 0:degree(g)
      print(io, coeff(g, j)*d)
      print(io, " ")
    end
    print(io, d)
    print(io, "\n")

    # print the elements
    for i in 1:length(A)

      f = polring(A[i])
      d = denominator(f)

      for j in 0:degree(K)-1
        print(io, coeff(f, j)*d)
        print(io, " ")
      end

      print(io, d)

      print(io, "\n")
    end
  end
end

@doc Markdown.doc"""
    write(file::String, A::Vector{nf_elem}, flag::ASCIString = "w") -> Nothing

Writes the elements of `A` to the file `file`. The first line are the coefficients of
the defining polynomial of the ambient number field. The following lines
contain the coefficients of the elements of `A` with respect to the power
basis of the ambient number field.

Unless otherwise specified by the parameter `flag`, the content of `file` will be
overwritten.
"""
function write(file::String, A::Vector{nf_elem}, flag::String = "w")
  f = open(file, flag)
  write(f, A)
  close(f)
end

# This function has a bad memory footprint
@doc Markdown.doc"""
    read(io::IO, K::AnticNumberField, ::Type{nf_elem}) -> Vector{nf_elem}

Given a file with content adhering the format of the `write` procedure,
this function returns the corresponding object of type `Vector{nf_elem}` such that
all elements have parent $K$.

**Example**

    julia> Qx, x = FlintQQ["x"]
    julia> K, a = NumberField(x^3 + 2, "a")
    julia> write("interesting_elements", [1, a, a^2])
    julia> A = read("interesting_elements", K, Hecke.nf_elem)
"""
function read(io::IO, K::AnticNumberField, ::Type{Hecke.nf_elem})
  Qx = parent(K.pol)

  A = Vector{nf_elem}()

  i = 1

  for ln in eachline(io)
    if ln[1] == '#'
      continue
    elseif i == 1
      # the first line read should contain the number field and will be ignored
      i = i + 1
    else
      coe = map(Hecke.fmpz, split(ln, " "))
      t = fmpz_poly(Array(slice(coe, 1:(length(coe) - 1))))
      t = Qx(t)
      t = divexact(t, coe[end])
      push!(A, K(t))
      i = i + 1
    end
  end

  return A
end

@doc Markdown.doc"""
    read(file::String, K::AnticNumberField, ::Type{nf_elem}) -> Vector{nf_elem}

Given a file with content adhering the format of the `write` procedure,
this function returns the corresponding object of type `Vector{nf_elem}` such that
all elements have parent $K$.

**Example**

    julia> Qx, x = FlintQQ["x"]
    julia> K, a = NumberField(x^3 + 2, "a")
    julia> write("interesting_elements", [1, a, a^2])
    julia> A = read("interesting_elements", K, Hecke.nf_elem)
"""
function read(file::String, K::AnticNumberField, ::Type{Hecke.nf_elem})
  f = open(file, "r")
  A = read(f, K, Hecke.nf_elem)
  close(f)
  return A
end

#TODO: get a more intelligent implementation!!!
@doc Markdown.doc"""
    splitting_field(f::fmpz_poly) -> AnticNumberField
    splitting_field(f::fmpq_poly) -> AnticNumberField

Computes the splitting field of $f$ as an absolute field.
"""
function splitting_field(f::fmpz_poly; do_roots::Bool = false)
  Qx = PolynomialRing(FlintQQ, parent(f).S, cached = false)[1]
  return splitting_field(Qx(f), do_roots = do_roots)
end

function splitting_field(f::fmpq_poly; do_roots::Bool = false)
  return splitting_field([f], do_roots = do_roots)
end

function splitting_field(fl::Vector{fmpz_poly}; coprime::Bool = false, do_roots::Bool = false)
  Qx = PolynomialRing(FlintQQ, parent(fl[1]).S, cached = false)[1]
  return splitting_field([Qx(x) for x = fl], coprime = coprime, do_roots = do_roots)
end

function splitting_field(fl::Vector{fmpq_poly}; coprime::Bool = false, do_roots::Bool = false)
  if !coprime
    fl = coprime_base(fl)
  end
  ffl = fmpq_poly[]
  for x = fl
    append!(ffl, collect(keys(factor(x).fac)))
  end
  fl = ffl
  r = []
  if do_roots
    r = [roots(x)[1] for x = fl if degree(x) == 1]
  end
  fl = fl[findall(x->degree(x) > 1, fl)]
  if length(fl) == 0
    if do_roots
      return FlintQQ, r
    else
      return FlintQQ
    end
  end
  K, a = number_field(fl[1])#, check = false, cached = false)

  @assert fl[1](a) == 0
  gl = [change_base_ring(K, fl[1])]
  gl[1] = divexact(gl[1], gen(parent(gl[1])) - a)
  for i=2:length(fl)
    push!(gl, change_base_ring(K, fl[i]))
  end

  if do_roots
    K, R = _splitting_field(gl, coprime = true, do_roots = Val{true})
    return K, vcat(r, [a], R)
  else
    return _splitting_field(gl, coprime = true, do_roots = Val{false})
  end
end

copy(f::fmpq_poly) = parent(f)(f)
gcd_into!(a::fmpq_poly, b::fmpq_poly, c::fmpq_poly) = gcd(b, c)

@doc Markdown.doc"""
    splitting_field(f::PolyElem{nf_elem}) -> AnticNumberField

Computes the splitting field of $f$ as an absolute field.
"""
splitting_field(f::PolyElem{nf_elem}; do_roots::Bool = false) = splitting_field([f], do_roots = do_roots)

function splitting_field(fl::Vector{<:PolyElem{nf_elem}}; do_roots::Bool = false, coprime::Bool = false)
  if !coprime
    fl = coprime_base(fl)
  end
  ffl = eltype(fl)[]
  for x = fl
    append!(ffl, collect(keys(factor(x).fac)))
  end
  fl = ffl
  r = []
  if do_roots
    r = [roots(x)[1] for x = fl if degree(x) == 1]
  end
  lg = [k for k = fl if degree(k) > 1]
  if length(lg) == 0
    if do_roots
      return base_ring(fl[1]), r
    else
      return base_ring(fl[1])
    end
  end

  K, a = number_field(lg[1])#, check = false)
  ggl = [map_coefficients(K, lg[1])]
  ggl[1] = divexact(ggl[1], gen(parent(ggl[1])) - a)

  for i = 2:length(lg)
    push!(ggl, map_coefficients(K, lg[i]))
  end
  if do_roots == Val{true}
    R = [K(x) for x = r]
    push!(R, a)
    Kst, t = PolynomialRing(Ks, cached = false)
    return _splitting_field(vcat(ggl, [t-y for y in R]), coprime = true, do_roots = Val{true})
  else
    return _splitting_field(ggl, coprime = true, do_roots = Val{false})
  end
end

function _splitting_field(fl::Vector{<:PolyElem{<:NumFieldElem}}; do_roots::Type{Val{T}} = Val{false}, coprime::Bool = false) where T
  if !coprime
    fl = coprime_base(fl)
  end
  ffl = eltype(fl)[]
  for x = fl
    append!(ffl, collect(keys(factor(x).fac)))
  end
  fl = ffl
  K = base_ring(fl[1])
  r = elem_type(K)[]
  if do_roots == Val{true}
    r = elem_type(K)[roots(x)[1] for x = fl if degree(x) == 1]
  end
  lg = eltype(fl)[k for k = fl if degree(k) > 1]
  if iszero(length(lg))
    if do_roots == Val{true}
      return K, r
    else
      return K
    end
  end

  K, a = number_field(lg[1])#, check = false)
  Ks, nk, mk = collapse_top_layer(K)

  ggl = [map_coefficients(mk, lg[1])]
  ggl[1] = divexact(ggl[1], gen(parent(ggl[1])) - preimage(nk, a))

  for i = 2:length(lg)
    push!(ggl, map_coefficients(mk, lg[i]))
  end
  if do_roots == Val{true}
    R = [mk(x) for x = r]
    push!(R, preimage(nk, a))
    Kst, t = PolynomialRing(Ks, cached = false)
    return _splitting_field(vcat(ggl, [t-y for y in R]), coprime = true, do_roots = Val{true})
  else
    return _splitting_field(ggl, coprime = true, do_roots = Val{false})
  end
end

function Base.:(^)(a::nf_elem, e::UInt)
  b = parent(a)()
  ccall((:nf_elem_pow, libantic), Nothing,
        (Ref{nf_elem}, Ref{nf_elem}, UInt, Ref{AnticNumberField}),
        b, a, e, parent(a))
  return b
end

@doc Markdown.doc"""
    normal_closure(K::AnticNumberField) -> AnticNumberField, NfToNfMor

The normal closure of $K$ together with the embedding map.
"""
function normal_closure(K::AnticNumberField)
  s = splitting_field(K.pol)
  r = roots(K.pol, s)[1]
  return s, hom(K, s, r, check = false)
end

function set_name!(K::AnticNumberField, s::String)
  Nemo.set_special(K, :name => s)
end

function set_name!(K::AnticNumberField)
  s = find_name(K)
  s === nothing || set_name!(K, string(s))
end

################################################################################
#
#  Is linearly disjoint
#
################################################################################

function islinearly_disjoint(K1::AnticNumberField, K2::AnticNumberField)
  if gcd(degree(K1), degree(K2)) == 1
    return true
  end
  d1 = numerator(discriminant(K1.pol))
  d2 = numerator(discriminant(K2.pol))
  if gcd(d1, d2) == 1
    return true
  end
  try
    OK1 = _get_maximal_order(K1)
    OK2 = _get_maximal_order(K2)
    if iscoprime(discriminant(K1), discriminant(K2))
      return true
    end
  catch e
    if !isa(e, AccessorNotSetError)
      rethrow(e)
    end
  end
  f = change_base_ring(K2, K1.pol)
  return isirreducible(f)
end

################################################################################
#
#  more general coercion, field lattice
#
################################################################################

Nemo.iscyclo_type(::NumField) = false

function force_coerce(a::NumField{T}, b::NumFieldElem, throw_error::Type{Val{S}} = Val{true}) where {T, S}
  if Nemo.iscyclo_type(a) && Nemo.iscyclo_type(parent(b))
    return force_coerce_cyclo(a, b, throw_error)::elem_type(a)
  end
  if absolute_degree(parent(b)) <= absolute_degree(a)
    c = find_one_chain(parent(b), a)
    if c !== nothing
      x = b
      for f = c
        @assert parent(x) == domain(f)
        x = f(x)
      end
      return x::elem_type(a)
    end
  end
  if throw_error === Val{true}
    throw(error("no coercion possible"))
  else
    return false
  end
end

@noinline function force_coerce_throwing(a::NumField{T}, b::NumFieldElem) where {T}
  if absolute_degree(parent(b)) <= absolute_degree(a)
    c = find_one_chain(parent(b), a)
    if c !== nothing
      x = b
      for f = c
        @assert parent(x) == domain(f)
        x = f(x)
      end
      return x::elem_type(a)
    else
      throw(error("no coercion possible"))
    end
  else
    throw(error("no coercion possible"))
  end
end

#(large) fields have a list of embeddings from subfields stored (special -> subs)
#this traverses the lattice downwards collecting all chains of embeddings
function collect_all_chains(a::NumField, filter::Function = x->true)
  s = get_special(a, :subs)
  if s === nothing
    return s
  end
  all_chain = Dict{UInt, Array{Any}}(objectid(domain(f)) => [f] for f = s if filter(f))
  if isa(base_field(a), NumField)
    all_chain[objectid(base_field(a))] = [MapFromFunc(x->a(x), base_field(a), a)]
  end
  new_k = Any[domain(f) for f = s]
  while length(new_k) > 0
    k = pop!(new_k)
    s = get_special(k, :subs)
    s === nothing && continue
    for f in s
      if filter(domain(f))
        o = objectid(domain(f))
        if haskey(all_chain, o)
          continue
        end
        @assert !haskey(all_chain, o)
        all_chain[o] = vcat([f], all_chain[objectid(codomain(f))])
        @assert !(o in new_k)
        push!(new_k, domain(f))
        if isa(base_field(domain(f)), NumField)
          b = base_field(domain(f))
          ob = objectid(b)
          if !haskey(all_chain, ob)
            g = MapFromFunc(x->domain(f)(x), b, domain(f))
            all_chain[ob] = vcat([g], all_chain[objectid(domain(f))])
            push!(new_k, b)
          end
        end
      end
    end
  end
  return all_chain
end

#tries to find one chain (array of embeddings) from a -> .. -> t
function find_one_chain(t::NumField, a::NumField)
  s = get_special(a, :subs)
  if s === nothing
    return s
  end
  ot = objectid(t)
  all_chain = Dict{UInt, Array{Any}}(objectid(domain(f)) => [f] for f = s)
  if isa(base_field(a), NumField)
    all_chain[objectid(base_field(a))] = [MapFromFunc(x->a(x), base_field(a), a)]
  end
  new_k = Any[domain(f) for f = s]
  if haskey(all_chain, ot)
    return all_chain[ot]
  end
  new_k = Any[domain(f) for f = s]
  while length(new_k) > 0
    k = pop!(new_k)
    s = get_special(k, :subs)
    s === nothing && continue
    for f in s
      o = objectid(domain(f))
      if o == ot
        return vcat([f], all_chain[objectid(codomain(f))])
      end
      if o in keys(all_chain)
        continue
      end
      @assert !haskey(all_chain, o)
      all_chain[o] = vcat([f], all_chain[objectid(codomain(f))])
      @assert !(o in new_k)
      push!(new_k, domain(f))
      if isa(base_field(domain(f)), NumField)
        b = base_field(domain(f))
        ob = objectid(b)
        if !haskey(all_chain, ob)
          g = MapFromFunc(x->domain(f)(x), b, domain(f))
          all_chain[ob] = vcat([g], all_chain[objectid(domain(f))])
          push!(new_k, b)
        end
        if ob == ot
          return all_chain[ob]
        end
      end
    end
  end
  return nothing
end

@doc Markdown.doc"""
    embed(f::Map{<:NumField, <:NumField})

Registers `f` as a canonical embedding from the domain into the co-domain.
Once this embedding is registered, it cannot be changed.
"""
function embed(f::Map{<:NumField, <:NumField})
  d = domain(f)
  c = codomain(f)
  if c == d
    return
  end
  @assert absolute_degree(d) <= absolute_degree(c)
  cn = find_one_chain(d, c)
  if cn !== nothing
    if issimple(d)
      cgend = force_coerce(c, gen(d))
      if cgend != f(gen(d))
        error("different embedding already installed")
        return
      end
    else
      if any(x->c(x) != f(x), gens(d))
        error("different embedding already installed")
      end
    end
  end
  s = get_special(c, :subs)
  if s === nothing
    s = Any[f]
  else
    push!(s, f)
  end
  set_special(c, :subs => s)
  s = get_special(c, :sub_of)

  if s === nothing
    s = Any[WeakRef(c)]
  else
    push!(s, WeakRef(c))
  end

  set_special(d, :sub_of => s)
end

@doc Markdown.doc"""
    hasembedding(F::NumField, G::NumField) -> Bool

Checks if an embedding from $F$ into $G$ is already known.
"""
function hasembedding(F::NumField, G::NumField)
  if F == G
    return true
  end
  if absolute_degree(G) % absolute_degree(F) != 0
    return false
  end
  cn = find_one_chain(d, c)
  return cn !== nothing
end

#in (small) fields, super fields are stored via WeakRef's
# special -> :sub_of
#this find all superfields registered
function find_all_super(A::NumField, filter::Function = x->true)
  s = get_special(A, :sub_of)
  s === nothing && return Set([A])

  ls = length(s)
  filter!(x->x.value !== nothing, s)
  if length(s) < ls #pruning old superfields
    set_special(A, :sub_of)
  end

  #the gc could(?) run anytime, so even after the pruning above
  #things could get deleted

  all_s = Set([x.value for x = s if x.value !== nothing && filter(x.value)])
  new_s = copy(all_s)
  while length(new_s) > 0
    B = pop!(new_s)
    s = get_special(B, :sub_of)
    s === nothing && continue
    ls = length(s)
    filter!(x->x.value !== nothing, s)
    if length(s) < ls
      set_special(B, :sub_of)
    end
    for x = s
      v = x.value
      if v !== nothing && filter(v)
        push!(new_s, v)
        push!(all_s, v)
      end
    end
  end
  return all_s
end

#finds a common super field for A and B, using the weak-refs
# in special -> :sub_of
function common_super(A::NumField, B::NumField)
  A === B && return A
  if Nemo.iscyclo_type(A) && Nemo.iscyclo_type(B)
    return cyclotomic_field(lcm(get_special(A, :cyclotomic_field), get_special(B, :cyclotomic_field)))[1]
  end

  c = intersect(find_all_super(A), find_all_super(B))
  first = true
  m = nothing
  for C = c
    if first
      m = C
      first = false
    else
      if absolute_degree(C) < absolute_degree(m)
        m = C
      end
    end
  end
  return m
end

function common_super(a::NumFieldElem, b::NumFieldElem)
  C = common_super(parent(a), parent(b))
  if C === nothing
    return C, C
  end
  return C(a), C(b)
end

#tries to find a common parent for all "a" and then calls op on it.
function force_op(op::T, throw_error::Type{Val{S}}, a::NumFieldElem...) where {T <: Function, S}
  C = parent(a[1])
  for b = a
    C = common_super(parent(b), C)
    if C === nothing
      if throw_error === Val{true}
        throw(error("no common parent known"))
      else
        return nothing
      end
    end
  end
  return op(map(C, a)...)
end

@doc Markdown.doc"""
    embedding(k::NumField, K::NumField) -> Map

Assuming $k$ is known to be a subfield of $K$, return the embedding map.
"""
function embedding(k::NumField, K::NumField)
  if issimple(k)
    return hom(k, K, K(gen(k)))
  else
    return hom(k, K, map(K, gens(k)))
  end
end

function force_coerce_cyclo(a::AnticNumberField, b::nf_elem, throw_error::Type{Val{T}} = Val{true}) where {T}
#  Base.show_backtrace(stdout, backtrace())
  fa = get_special(a, :cyclo)
  fb = get_special(parent(b), :cyclo)

  if fa % fb == 0 #coerce up, includes fa == fb
    #so a = p(z) for p in Q(x) and z = gen(parent(b))
    q = divexact(fa, fb)
    c = parent(a.pol)()
    if fb == 2 # if the field is linear, elem_length is not well-defined
      return a(coeff(b, 0))
    end
    for i=0:b.elem_length
      setcoeff!(c, i*q, coeff(b, i))
    end
    return a(c)
  end
  if isodd(fa) && (2*fa) % fb == 0
    #so a = cyclo(odd) = cyclo(2*odd) (and thus fb = 2*odd as well)
    #write b  = f(zeta_(2l))
    # z_2l = -z_l^invmod(2, l)
    # z_l = z_a^div(fa, l) so
    # z_2l = -z_a^(div(fa, l)*invmod(2, l))
    s = divexact(2*fa, fb)*invmod(2, divexact(fb, 2)) % fa
    si = 1
    c = parent(a.pol)()
    for i=0:b.elem_length
      setcoeff!(c, s*i, si*coeff(b, i))
      si *= -1
    end
    return a(c)
  end

  if fb % fa == 0 ||
    (isodd(fb) && (2*fb) % fa == 0)

    zb = gen(parent(b))

    ff = parent(parent(b).pol)(b)
    if isodd(fb) && (2*fb) % fa == 0
      ff = inflate(ff, 2)
      #f(x^2)
      zb = -zb^invmod(2, fb)
      fb *= 2
    end
    za = zb^divexact(fb, fa)
    q = divexact(fb, fa)

    cb = [i for i=1:fb if gcd(i, fb) == 1] # the "conjugates" in the large field
    ca = [[i for i = cb if i % fa == j] for j=1:fa if gcd(j, fa) == 1] #broken into blocks

    #in general one could test first if the evaluation is constant on a block
    #equivalently, if the element is Galois invariant under the fix group of a.
    #the result of the interpolation is supposed to be over Q, so we could
    #do this modulo deg-1-primes as well
    #using a fast(er) interpolation would be nice as well
    #but, this is avoiding matrices, so teh complexity is better
    #
    #Idea
    # b is a poly in Qx, evaluating at gen(a)... will produce the conjugates
    # b in a is also a poly, of smaller degree producing the same conjugates
    # so I compute the conjugates from the large field and interpolate them to get
    # the small degree poly
    #actually, since we're using roots of one, we proably should use FFT techniques

    ex = [x[1] for x = ca]
    ky = PolynomialRing(parent(b), cached = false)[1]
    f = interpolate(ky, [(za)^(i) for i=ex],
                        [ff(zb^(i)) for i=ex])
    g = parent(ff)()
    for i=0:length(f)
      c = coeff(f, i)
      if !isrational(c)
        if throw_error === Val{true}
          throw(error("no coercion possible"))
        else
          return false
        end
      end
      setcoeff!(g, i, FlintQQ(c))
    end
    return ba = a(g)
  end #missing: (?) b could be in a subfield and thus still in a
  if throw_error === Val{true}
    throw(error("no coercion possible"))
  end
end

(::FlintRationalField)(a::nf_elem) = (isrational(a) && return coeff(a, 0)) || error("not a rational")
(::FlintIntegerRing)(a::nf_elem) = (isinteger(a) && return numerator(coeff(a, 0))) || error("not an integer")


text:
# Tests for economic optimization problems. Tet only basic functionality and errors as other test are in specific DEA models.
@testset "EconOptimProblems" begin

    # ------------------
    # Cost Minimization
    # ------------------
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6	8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]
    W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1]

    Xtarget, clambda = deamincost(X, Y, W)

    @test Xtarget == 2 * ones(8,2)
    @test all(clambda[:,1] .== 1)

    @test_throws DimensionMismatch deamincost([1; 2 ; 3], [4 ; 5], [1; 1; 1]) #  Different number of observations
    @test_throws DimensionMismatch deamincost([1; 2; 3], [4; 5; 6], [1; 2; 3; 4]) # Different number of observation in prices
    @test_throws DimensionMismatch deamincost([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3]) # Different number of input prices and inputs
    @test_throws ArgumentError deamincost([1; 2; 3], [4; 5; 6], [1; 2; 3], rts = :Error) # Invalid returns to scale
    @test_throws ArgumentError deamincost([1; 2; 3], [4; 5; 6], [1; 2; 3], dispos = :Error) # Invalid disposability

    # ------------------
    # Revenue Maximization
    # ------------------
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]
    P = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1]

    Ytarget, rlambda = deamaxrevenue(X, Y, P)

    @test Ytarget == 7 * ones(8,2)
    @test all(clambda[:,1] .== 1)

    @test_throws DimensionMismatch  deamaxrevenue([1; 2 ; 3], [4 ; 5], [1; 1; 1]) #  Different number of observations
    @test_throws DimensionMismatch  deamaxrevenue([1; 2; 3], [4; 5; 6], [1; 2; 3; 4]) # Different number of observation in prices
    @test_throws DimensionMismatch  deamaxrevenue([1; 2; 3], [4 4; 5 5; 6 6], [4 4 4; 5 5 5; 6 6 6]) # Different number of output prices and outputs
    @test_throws ArgumentError deamaxrevenue([1; 2; 3], [4; 5; 6], [1; 2; 3], rts = :Error) # Invalid returns to scale
    @test_throws ArgumentError deamaxrevenue([1; 2; 3], [4; 5; 6], [1; 2; 3], dispos = :Error) # Invalid disposability

    # ------------------
    # Profit Maximization
    # ------------------
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]
    W = [1; 1; 1; 1; 1; 1; 1; 1]
    P = [2; 2; 2; 2; 2; 2; 2; 2]

    Xtarget, Ytarget, plambda = deamaxprofit(X, Y, W, P)

    @test Xtarget == 8 * ones(8,1)
    @test Ytarget == 8 * ones(8,1)
    @test all(plambda[:,3] .== 1)

    @test_throws DimensionMismatch deamaxprofit([1; 2 ; 3], [4 ; 5], [1; 1; 1], [4; 5]) #  Different number of observations
    @test_throws DimensionMismatch deamaxprofit([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], [4; 5; 6]) # Different number of observation in input prices
    @test_throws DimensionMismatch deamaxprofit([1; 2; 3], [4; 5; 6], [1; 2; 3], [4; 5; 6; 7]) # Different number of observation in output prices
    @test_throws DimensionMismatch deamaxprofit([1 1; 2 2; 3 3], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], [4; 5; 6]) # Different number of input prices and inputs
    @test_throws DimensionMismatch deamaxprofit([1; 2; 3], [4 4; 5 5; 6 6], [1; 2; 3], [4 4 4; 5 5 5; 6 6 6]) # Different number of oputput prices and outputs

end


text:
# Note: All the estimates and inequalities below only apply
# for our iterated integrals because the Levy area and the
# corresponding error terms are skew-symmetric!

"""
    abstract type AbstractErrorNorm end

Abstract type for different kind of errors one might consider.
The most important are `MaxL2` and `FrobeniusL2`. These are actually
special cases of the `MaxLp{p}` and `lqLp{p,q}` norms:
```julia
const MaxL2 = MaxLp{2}
const FrobeniusL2 = lqLp{2,2}
```

All currently defined norms:
```jldoctest; setup=:(using InteractiveUtils; using LevyArea)
julia> subtypes(LevyArea.AbstractErrorNorm)
3-element Vector{Any}:
 LevyArea.LpMax
 LevyArea.MaxLp
 LevyArea.lqLp
```
"""
abstract type AbstractErrorNorm end

# maximum of entry-wise Lp(Ω) norms
# this requires new Lp convergence proofs
# and should be defined for each new algorithm
struct MaxLp{p} <: AbstractErrorNorm end

function errcoeff(m, q_12, stepsize, alg::AbstractIteratedIntegralAlgorithm, ::MaxLp{p}) where {p}
    maxqq = maximum(q_12[i]*q_12[j] for i=1:m for j=1:i-1)
    return maxqq * errcoeff(m, stepsize, alg, MaxLp{p}())
end

# lq norm of all entry-wise Lp(Ω) norms
# = (Σᵢⱼ||Mᵢⱼ||_Lp(Ω)^q)^(1/q)
# we can get this from the MaxLp norms
# p=q=2 is the FrobeniusL2 norm
struct lqLp{p,q} <: AbstractErrorNorm end

function errcoeff(m, stepsize, alg::AbstractIteratedIntegralAlgorithm, ::lqLp{p,q}) where {p,q}
    return (m^2-m)^(1/q) * errcoeff(m, stepsize, alg, MaxLp{p}())
end
function errcoeff(m, q_12, stepsize, alg::AbstractIteratedIntegralAlgorithm, ::lqLp{p,q}) where {p,q}
    trQq_sq = abs2(sum(x->x^q, q_12))
    tr_Qqsq = sum(x->x^2q, q_12)
    return (trQq_sq-tr_Qqsq)^(1/q) * errcoeff(m, stepsize, alg, MaxLp{p}())
end

# Lp norm of maximum of matrix
# can be bounded by the lqLp norm with q=p
# (this relies on the fact that lp-Lp == Lp-lp)
# the factor 1/2^1/p comes from the antisymmetry
struct LpMax{p} <: AbstractErrorNorm end

function errcoeff(m, stepsize, alg::AbstractIteratedIntegralAlgorithm, ::LpMax{p}) where {p}
    return 1/(2^(1/p)) * errcoeff(m, stepsize, alg, lqLp{p,p}())
end
function errcoeff(m, q_12, stepsize, alg::AbstractIteratedIntegralAlgorithm, ::LpMax{p}) where {p}
    return 1/(2^(1/p)) * errcoeff(m, q_12, stepsize, alg, lqLp{p,p}())
end

# most important/common norms
const MaxL2 = MaxLp{2}
const FrobeniusL2 = lqLp{2,2}

max_stars_repo_name: mrbauff/MasterThesis

max_issues_repo_name: mrbauff/MasterThesis