text:
module ArtifactHelpers

using Pkg
using Pkg.Artifacts
using Pkg.PlatformEngines #Note: supplies unpack.
using Pkg.BinaryPlatforms
using SHA
using Downloads: download

using ZipFile

import Base: SHA1

export
    File, AutoDownloadable, Zip, Processed,
    build_artifact!, initialise_artifact, setup

const AUTODOWNLOADABLES = [".gz"]

TYPES = []

function __init__()
    global TYPES = [(autodownloadable, AutoDownloadableEntry), (iszip, Zip), (downloadable, File), (islocal, Processed) ] #Note: order is relevant.
    @static isdefined(PlatformEngines, :probe_platform_engines!) && probe_platform_engines!()
end

abstract type Entry end
abstract type DownloadableEntry <: Entry end
abstract type CustomEntry <: Entry end

function _url(url::String)
    return [url]
end

function _url(urls::Vector{String})
    return urls
end

mutable struct AutoDownloadableEntry <: DownloadableEntry
    artifact_name::String
    meta::Dict{String,Any}
    urls::Vector{String}
    platform::Union{Platform,Nothing}
    lazy::Bool
    new_download_info::Vector{Tuple}

    AutoDownloadableEntry(url; platform = nothing, lazy = true) = new(basename(first(_url(url))), Dict(), _url(url), platform, lazy, Vector{Tuple}())
    AutoDownloadableEntry(artifact_name, url; platform = nothing, lazy = true) = new(artifact_name, Dict(), _url(url), platform, lazy, Vector{Tuple}())
    AutoDownloadableEntry(artifact_name, dict::Dict; platform = nothing, lazy = true) = new(artifact_name, dict, Vector{String}(), platform, lazy, Vector{Tuple}())
end

const AutoDownloadable = AutoDownloadableEntry

mutable struct File <: DownloadableEntry
    artifact_name::String
    meta::Dict{String,Any}
    platform::Union{Platform,Nothing}

    File(url; platform = nothing) = new(basename(url), Dict("url"=>url), platform)
    File(artifact_name, url; platform = nothing) = new(artifact_name, Dict("url"=>url), platform)
    File(artifact_name, dict::Dict; platform = nothing) = new(artifact_name, dict, platform)
end

mutable struct Zip <: DownloadableEntry
    artifact_name::String
    meta::Dict{String,Any}
    platform::Union{Platform,Nothing}

    Zip(url; platform = nothing) = new(basename(url), Dict("url"=>url), platform)
    Zip(artifact_name, url; platform = nothing) = new(artifact_name, Dict("url"=>url), platform)
    Zip(artifact_name, dict::Dict; platform = nothing) = new(artifact_name, dict, platform)
end

mutable struct Processed <: CustomEntry
    artifact_name::String
    meta::Dict{String,Any}
    platform::Union{Platform,Nothing}

    Processed(artifact_name; platform = nothing) = new(artifact_name, Dict(), platform)
    Processed(artifact_name, dict::Dict; platform = nothing) = new(artifact_name, dict, platform)
end

function name(entry::Entry)
    return entry.artifact_name
end

function metadata(entry::Entry)
    return entry.meta
end

function Base.haskey(entry::E, key) where E <: Entry
    return haskey(metadata(entry), key)
end

function Base.get(entry::Entry, key)
    return metadata(entry)[key] #Note: default is not provided, we want an error if the key does not exist.
end

function Base.setindex!(entry::Entry, value, key...)
    return setindex!(metadata(entry), value, key...)
end

function hassha256(entry::E) where E <: Entry

    #Note: the sha256 will be in the meta.

    if haskey(entry, "sha256")
        return true
    end

    if haskey(entry, "download")
        for download_entry in get(entry, "download")

            if isa(download_entry, Dict) && haskey(download_entry, "sha256")
                return true
            end

            @error "Unaccounted download_entry" typeof(download_entry) download_entry #TODO: Handle Array of Dicts?

        end
    end

    return false

end

function record_sha256!(entry::DownloadableEntry, url, str_hash) #TODO: make consistent with possibility for vector of urls.

    setindex!(entry, str_hash, "sha256")

    return entry
end

function record_sha256!(entry::AutoDownloadableEntry, url, str_hash)
    push!(entry.new_download_info, (url, str_hash))

    return entry
end

function isurl(str::String)
    check = "http"
    len = length(check)
    return length(str) >= len && str[1:len] == check
end

function iszip(path::AbstractString)

    (_, ext) = splitext(path)

    return isurl(path) && ext == ".zip"
end

function iszip(meta::Dict)
    return haskey(meta, "url") && iszip(meta["url"])
end

function autodownloadable(entry::String)
    (_, ext) = splitext(entry)

    return isurl(entry) && in(ext, AUTODOWNLOADABLES)
end

function autodownloadable(meta::Dict)
    return haskey(meta, "download") #TODO: subject to change as custom types converge.
end

function downloadable(entry::String)
    return isurl(entry) #TODO: find and use official url checker.
end

function downloadable(meta::Dict)
    return haskey(meta, "url") #Note: needs to be called after autodownloadable and zip!
end

function islocal(entry::String)
    return !isurl(entry) #TODO: find and use official url checker.
end

function islocal(meta::Dict)
    return !haskey(meta, "url") #Note: needs to be called after autodownloadable and zip!
end

struct Verifiable{x} end
Verifiable(x::Bool) = Verifiable{x}()

function isverifiable(entry::Entry)
    return hassha256(entry)
end

function isverifiable(url::String)
    return false
end

function _download(url::AbstractString, dest::AbstractString; verbose::Bool = false)
    return download(url, dest; verbose = verbose) #Note: throws an error if unsuccessfull, otherwise returns path.
end

function _download(url::AbstractString, dest::AbstractString, tarball_hash; verbose::Bool = false)
    downloaded = _download(url, dest; verbose = verbose) #Note: throws an error if unsuccessfull, otherwise returns path.

    result = Pkg.PlatformEngines.verify(dest, tarball_hash, verbose = verbose)

    result || error("Unexpected download.")

    return downloaded
end

function _download(entry::Entry, dest::AbstractString; kwargs...)
    return _download(entry, dest, Verifiable(isverifiable(entry)); kwargs...)
end

function _download(entry::AutoDownloadableEntry, dest::AbstractString, ::Verifiable{false}; kwargs...) #

    meta = metadata(entry)

    # Attempt to download from all sources.
    for url in entry.urls
        try

            downloaded = _download(url, dest; kwargs...) #Note: throws erro if unsuccessfull.
            record_sha256!(entry, url, file_hash(dest))
            return downloaded

        catch e

            if isa(e, InterruptException) || isa(e, MethodError) || isa(e, UndefVarError)
                rethrow(e)
            end
            # If something went wrong during download, continue.
            continue
        end

    end

    error("Nothing downloaded.")

end

function _download(entry::AutoDownloadableEntry, dest::AbstractString, ::Verifiable{true}; kwargs...)

    meta = metadata(entry)

    # Attempt to download from all sources.
    for download_entry in meta["download"]

        url = download_entry["url"]
        tarball_hash = download_entry["sha256"]

        try
            downloaded = _download(url, dest, tarball_hash; kwargs...) #Note: throws erro if unsuccessfull.
            return downloaded
        catch e

            if isa(e, InterruptException) || isa(e, MethodError) || isa(e, UndefVarError)
                rethrow(e)
            end
            # If something went wrong during download, continue.
            continue
        end
    end

    error("Nothing downloaded.")

end

function _download(entry::DownloadableEntry, dest::AbstractString, ::Verifiable{false}; kwargs...)
    url = get(entry, "url") #TODO: vector of URLs.
    downloaded = _download(url, dest; kwargs...)
    record_sha256!(entry, url, file_hash(dest))
    return downloaded
end

function _download(entry::DownloadableEntry, dest::AbstractString, ::Verifiable{true}; kwargs...)
    return _download(get(entry, "url"), dest, get(entry, "sha256"); kwargs...)
end

function unpack(entry::AutoDownloadableEntry, src::AbstractString, dest::AbstractString; verbose::Bool = false)
    return Pkg.PlatformEngines.unpack(src, dest, verbose = verbose)
end

function unpack(entry::Zip, src::AbstractString, dest::AbstractString; verbose::Bool = false)
    return unzip(src, dest, verbose = verbose)
end

function unpack(entry::File, src::AbstractString, dest::AbstractString; verbose::Bool = false)
    verbose && @info "Copying." src dest
    return cp(src, joinpath(dest, name(entry))) #Note: src is expected for cleanup.
end

function unpack(entry::Entry, src::AbstractString, dest::AbstractString; verbose::Bool = false)
    verbose && @info "Nothing to unpack."
    return
end

function setup(str::String)

    for (test, type) in TYPES
        if test(str)
            return type(str)
        end
    end

    error("Unknown entry type.")

end

function setup(artifact_name, artifacts_toml)

    meta = artifact_meta(artifact_name, artifacts_toml) #TODO: platform.

    for (test, type) in TYPES
        if test(meta)
            return type(artifact_name, meta)
        end
    end

    error("Unknown entry type.")

end

function record!(artifacts_toml::AbstractString, entry::Entry) #TODO: type Entry actually captures custom toml entries.

    isfile(artifacts_toml) || error("Artifacts.toml does not exist at specified path:", artifacts_toml)

    # Load toml.
    dict_artifact = Pkg.Artifacts.parse_toml(artifacts_toml)

    # Obtain artifact's metadata.
    meta = get(dict_artifact, name(entry), Dict())::Dict{String,Any}

    # Merge with whatever is on file.
    meta = merge(meta, metadata(entry)) #TODO: what does this do with Dicts of Dicts?

    # Overwrite artifact's metadata.
    dict_artifact[name(entry)] = meta

    # Spit it out onto disk.
    open(artifacts_toml, "w") do io
        Pkg.TOML.print(io, dict_artifact, sorted = true)
    end

    return nothing
end

function record!(artifacts_toml::AbstractString, entry::AutoDownloadableEntry)
    return nothing
end

function process(entry::DownloadableEntry; force::Bool = false, verbose::Bool = false)

    tree_hash = create_artifact() do path_artifact #Note: this will create an artifact that is ready for use.

        mktemp() do file_tmp, io

            _download(entry, file_tmp, verbose = verbose)

            unpack(entry, file_tmp, path_artifact, verbose = verbose)

        end
    end

    return tree_hash
end

function build_artifact!(artifacts_toml::String, entry::Entry, process_func::Function = process; force::Bool = false, verbose::Bool = false)

    if !isfile(artifacts_toml)
        error("Artifacts.toml does not exist at specified path: ", artifacts_toml)
    end

    artifact_name = name(entry)

    # Obtain artifact's recorded hash.
    tree_hash = artifact_hash(artifact_name, artifacts_toml)

    if tree_hash === nothing
        @warn "Building and binding a new artifact." artifact_name

        tree_hash = process_func(entry, verbose = verbose)

        bind_artifact!(artifacts_toml, entry, tree_hash; verbose = verbose)

        @info "Built $artifact_name." tree_hash
        return tree_hash
    end

    if force === true
        @warn "Forcefully rebuilding artifact." artifact_name

        tree_hash = process_func(entry, force = force, verbose = verbose)

        bind_artifact!(artifacts_toml, entry, tree_hash; force = force, verbose = verbose)

        @info "Forcefully built $artifact_name." tree_hash
        return tree_hash
    end

    # Setup the artifact if it does not exist on disk.
    if !artifact_exists(tree_hash)
        @warn "Rebuilding artifact." artifact_name

        setup_hash = process_func(entry, verbose = verbose) #Note: not forcing during initialisation.
        setup_hash == tree_hash || error("Hash $setup_hash of setup artifact does not match the entry for \"$artifact_name\".")

        @info "Rebuilt $artifact_name." tree_hash
        return tree_hash
    end

    @info "Skipped build of $artifact_name." tree_hash

    return tree_hash

end

function build_artifact!(artifacts_toml::String, str::String, process_func::Function = process; kwargs...)
    return build_artifact!(artifacts_toml, setup(str), process_func; kwargs...)
end

function build_artifact!(kernel::Function, artifacts_toml::String, entry; kwargs...)

    function wrapped_kernel(entry::Entry; kwargs...) #Note: kwargs captures and diffuses.

        tree_hash = create_artifact() do path_artifact #Note: this will create an artifact that is ready for use.
            kernel(path_artifact)
        end

        return tree_hash
    end

    return build_artifact!(artifacts_toml, entry, wrapped_kernel; kwargs...)
end

function Pkg.Artifacts.bind_artifact!(artifacts_toml::AbstractString, entry::AutoDownloadableEntry, tree_hash::SHA1; force::Bool = false, verbose::Bool = false)

    # bind_artifact!(artifacts_toml::String, name::String, hash::SHA1; platform::Union{Platform,Nothing} = nothing, download_info::Union{Vector{<:Tuple},Nothing} = nothing, lazy::Bool = false, force::Bool = false)

    download_info = entry.new_download_info #TODO: merge with existing download_info?

    # Bind acquired artifact.
    bind_artifact!(artifacts_toml, name(entry), tree_hash, force = force, lazy = entry.lazy, download_info = download_info, platform = entry.platform)

    return nothing

end

function Pkg.Artifacts.bind_artifact!(artifacts_toml::AbstractString, entry::Entry, tree_hash::SHA1; force::Bool = false, verbose::Bool = false)

    # bind_artifact!(artifacts_toml::String, name::String, hash::SHA1; platform::Union{Platform,Nothing} = nothing, download_info::Union{Vector{<:Tuple},Nothing} = nothing, lazy::Bool = false, force::Bool = false)

    # Bind acquired artifact.
    bind_artifact!(artifacts_toml, name(entry), tree_hash, force = force, lazy = false, platform = entry.platform)

    # Additionally record meta.
    record!(artifacts_toml, entry)

    return nothing

end

function initialise_artifact(artifacts_toml::String, artifact_name::String, process_func::Function = process; verbose::Bool = false)

    # Allow __init__ function to run when Artifacts.toml file does not exist.
    if !isfile(artifacts_toml)
        @warn "Artifacts.toml does not exist at specified path." artifacts_toml
        return nothing
    end

    # Obtain artifact's recorded hash.
    tree_hash = artifact_hash(artifact_name, artifacts_toml)

    # Allow __init__ function to run when the artifact entry does not exist in Artifacts.toml.
    if tree_hash == nothing
        @warn "An artifact entry for \"$artifact_name\" does not exist in Artifacts.toml."
        return nothing
    end

    # Setup the artifact if it does not exist on disk.
    if !artifact_exists(tree_hash)
        # setup_hash = setup_func(artifact_name, artifacts_toml, verbose = verbose)
        setup_hash = process_func(setup(artifact_name, artifacts_toml), force = false, verbose = verbose) #Note: not forcing during initialisation.
        setup_hash == tree_hash || error("Hash $setup_hash of setup artifact does not match the entry for \"$artifact_name\".")
    end

    verbose && @info "Initialised $artifact_name." tree_hash

    return tree_hash
end

function initialise_artifact(f::Function, args...; kwargs...)
    tree_hash = initialise_artifact(args...; kwargs...)

    if tree_hash === nothing
        return nothing
    end

    return f(tree_hash)
end

function file_hash(path)
    return open(path) do file
        bytes2hex(sha2_256(file))
    end
end

function unzip(src, dest; verbose::Bool = false)

    if !isfile(src)
        error("Source does not exist.")
    end

    if !isdir(dest)
        error("Destination does not exist.")
    end

    verbose && @info "Unzipping." src dest

    r = ZipFile.Reader(src)

    try
        for f in r.files

            verbose && (println("Filename: $(f.name)"))

            if f.method == ZipFile.Store
                mkpath(joinpath(dest, f.name))
            end

            if f.method == ZipFile.Deflate
                path = joinpath(dest, f.name)
                mkpath(dirname(path))
                !isdirpath(f.name) && write(path, read(f))
            end

        end
    finally
        close(r)
    end

    return dest
end

end # module


text:
@testset "grow" begin
    @testset "grow_last" begin
        for (m,f,s) in ((Axis(UnitMRange(1, 10)), Axis(1:10), Axis(UnitSRange(1, 10))),
                        (SimpleAxis(UnitMRange(1, 10)), SimpleAxis(1:10), SimpleAxis(UnitSRange(1, 10))))
            x = @inferred(grow_last(m, 2))
            @test m == 1:10
            @test x == 1:12

            x = @inferred(grow_last(f, 2))
            @test f == 1:10
            @test x == 1:12

            x = @inferred((s -> grow_last(s, 2))(s))

            @test s == 1:10
            @test x == 1:12
        end
    end

    @testset "grow_last!" begin
        for m in (Axis(UnitMRange(1, 10), UnitMRange(1, 10)), SimpleAxis(UnitMRange(1, 10)))
            x = @inferred(grow_last!(m, 2))
            @test m == 1:12
            @test x == 1:12
        end
    end

    # FIXME?
    #=
    @testset "grow_first" begin
        for (m,f,s) in ((Axis(UnitMRange(1, 10)), Axis(1:10), Axis(UnitSRange(1, 10))),
                        (SimpleAxis(UnitMRange(1, 10)), SimpleAxis(1:10), SimpleAxis(UnitSRange(1, 10))))
            m,f,s = UnitMRange(1, 10), 1:10, UnitSRange(1, 10)
            x = @inferred(grow_first(m, 2))
            @test m == 1:10
            @test x == -1:10

            x = @inferred(grow_first(f, 2))
            @test f == 1:10
            @test x == -1:10

            x = @inferred((s -> grow_first(s, 2))(s))
            @test s == 1:10
            @test x == -1:10
        end
    end

    @testset "grow_first!" begin
        for m in (Axis(UnitMRange(1, 10), UnitMRange(1, 10)), SimpleAxis(UnitMRange(1, 10)))
            x = @inferred(grow_first!(m, 2))
            @test m == -1:10
            @test x == -1:10
        end
    end
    =#
end

#=
@testset "shrink" begin
    @testset "shrink_last" begin
        for (m,f,s) in ((Axis(UnitMRange(1, 10)), Axis(1:10), Axis(UnitSRange(1, 10))),
                        (SimpleAxis(UnitMRange(1, 10)), SimpleAxis(1:10), SimpleAxis(UnitSRange(1, 10))))
            x = @inferred(shrink_last(m, 2))
            @test m == 1:10
            @test x == 1:8

            x = @inferred(shrink_last(f, 2))
            @test f == 1:10
            @test x == 1:8

            x = @inferred((s -> shrink_last(s, 2))(s))
            @test s == 1:10
            @test x == 1:8
        end
    end

    @testset "shrink_last!" begin
        for m in (Axis(UnitMRange(1, 10), UnitMRange(1, 10)), SimpleAxis(UnitMRange(1, 10)))
            x = @inferred(shrink_last!(m, 2))
            @test m == 1:8
            @test x == 1:8
        end
    end

    @testset "shrink_first" begin
        for (m,f,s) in (
            (Axis(UnitMRange(1, 10)), Axis(1:10), Axis(UnitSRange(1, 10))),
            (SimpleAxis(UnitMRange(1, 10)), SimpleAxis(1:10), SimpleAxis(UnitSRange(1, 10))))
            x = @inferred(shrink_first(m, 2))
            @test m == 1:10
            @test x == 3:10

            x = @inferred(shrink_first(f, 2))
            @test f == 1:10
            @test x == 3:10

            x = @inferred((s -> shrink_first(s, 2))(s))
            @test s == 1:10
            @test x == 3:10
        end
    end

    @testset "shrink_first!" begin
        for m in (Axis(UnitMRange(1, 10), UnitMRange(1, 10)), SimpleAxis(UnitMRange(1, 10)))
            x = @inferred(shrink_first!(m, 2))
            @test m == 3:10
            @test x == 3:10
        end
    end
end

=#

ext: jl


text:
function print_frame(H, W)
    println("#"^W)
    for i in 1:H-2
        println("#", "."^(W-2), "#")
    end
    println("#"^W)
end

while true
    H, W = parse.(Int, split(readline()))
    (H == 0 && W == 0) && break
    print_frame(H, W)
    println()
end


text:
# Checking gradients and Hessian vector products

function check_single_grad(loss::Function, grad!::Function, datapoint::Vector, theta::Vector, lambda::Float64, eps::Float64)
  p = length(theta)
  delta = randn(p)
  deriv_approx = (loss(datapoint, theta+eps*delta, lambda) - loss(datapoint, theta-eps*delta, lambda))/(2*eps)
  grad = zeros(p)
  grad!(datapoint, theta, lambda, grad)
  deriv_exact = dot(grad, delta)
  return norm(deriv_approx - deriv_exact)
end

function check_grad(loss::Function, grad!::Function, data::Matrix, theta::Vector, lambda::Float64, eps::Float64)
  n = size(data, 2)
  p = length(theta)
  delta = randn(p)
  deriv_approx = (loss(data, theta+eps*delta, lambda) - loss(data, theta-eps*delta, lambda))/(2*eps)
  grad = zeros(p)
  for i = 1:n
    grad!(vec(data[:,i]), theta, lambda, grad)
  end
  deriv_exact = dot(grad/n, delta)
  return norm(deriv_approx - deriv_exact)
end

# Compute the full gradient over the whole data matrix
function full_grad(grad!::Function, data::Matrix, theta::Vector, lambda::Float64)
  n = size(data, 2)
  result = zeros(length(theta))
  for i = 1:n
    grad!(vec(data[:,i]), theta, lambda, result)
  end
  return result
end

function single_grad(grad!::Function, datapoint::Vector, theta::Vector, lambda::Float64)
  result = zeros(length(theta))
  grad!(datapoint, theta, lambda, result)
  return result
end

# plase check with check_gradient if the gradient function is correct before you use this function
function check_hvp(grad!::Function, hvp!::Function, data::Matrix, theta::Vector, lambda::Float64, eps::Float64)
  n = size(data, 2)
  p = length(theta)
  delta = randn(p)
  hvp_exact = zeros(p)
  for i = 1:n
    hvp!(vec(data[:,i]), theta, delta, lambda, hvp_exact)
  end
  hvp_approx = (full_grad(grad!, data, theta+eps*delta, lambda) - full_grad(grad!, data, theta-eps*delta, lambda))/(2*eps)
  return norm(hvp_exact - hvp_approx)
end

function check_single_hvp(grad!::Function, hvp!::Function, datapoint::Vector, theta::Vector, lambda::Float64, eps::Float64)
  p = length(theta)
  delta = randn(p)
  hvp_exact = zeros(p)
  hvp!(datapoint, theta, delta, lambda, hvp_exact)
  hvp_approx = (single_grad(grad!, datapoint, theta+eps*delta, lambda) - single_grad(grad!, datapoint, theta-eps*delta, lambda))/(2*eps)
  return norm(hvp_exact - hvp_approx)
end

function ad_hoc_check(mode::String)
  println("Checking gradients and hessian vector products")
  data = randn(100, 100)
  if mode == "svm"
    data[100,:] = 2.0*randbool(100) - 1
  end
  if mode == "logistic regression"
    data[100,:] = 2.0*randbool(100) - 1
  end
  theta = randn(99)
  gnorm = norm(check_grad(loss, grad!, data, theta, 1e-3, 1e-4))
  hnorm = norm(check_hvp(grad!, hvp!, data, theta, 1e-3, 1e-4))
  println(gnorm)
  @assert gnorm <= 1e-4
  println(hnorm)
  @assert hnorm <= 1e-4
end

max_stars_repo_name: pcmoritz/slbfgs

max_issues_repo_name: pcmoritz/slbfgs

max_forks_repo_name: pcmoritz/slbfgs


text:
using RandomNumbers.MersenneTwisters
if !@isdefined RandomNumbers
    include("../common.jl")
end

@info "Testing MersenneTwisters"
stdout_ = stdout
pwd_ = pwd()
cd(dirname(@__FILE__))
rm("./actual"; force=true, recursive=true)
mkpath("./actual")

@test seed_type(MT19937) == NTuple{RandomNumbers.MersenneTwisters.N, UInt32}

for mt_name in (:MT19937, )

    outfile = open(string(
        "./actual/check-$(lowercase("$mt_name")).out"
    ), "w")
    redirect_stdout(outfile)

    @eval $mt_name()
    x = @eval $mt_name(123)
    @test copyto!(copy(x), x) == x

    for i in 1:100
        @printf "%.9f\n" rand(x)
    end

    close(outfile)
end
redirect_stdout(stdout_)

compare_dirs("expected", "actual")
cd(pwd_)


text:
@testset "Spaces" begin
  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x^2 - 2, "a1")
  OK = maximal_order(K)
  p = 3*OK
  Kt, t = K["t"]

  E, b = NumberField(t^2 + 3)
  s = involution(E)

  F = GF(3)

  Hecke.change_base_ring(::Hecke.NfRel, ::Hecke.gfp_mat) = error("asd")
  @test_throws ErrorException hermitian_space(E, F[1 2; 2 1])

  Hecke.change_base_ring(::Hecke.NfRel, x::Hecke.gfp_mat) = x
  @test_throws ErrorException hermitian_space(E, F[1 2; 2 1])

  V = @inferred hermitian_space(E, 3)
  @test fixed_field(V) == base_field(E)
  @test sprint(show, V) isa String
  @test !isquadratic(V)
  @test ispositive_definite(V)
  @test gram_matrix(V) == identity_matrix(E,3)
  @test V !== hermitian_space(E, 3, cached = false)

  @test_throws ArgumentError hermitian_space(E, E[1 b; b 1])
  @test_throws ArgumentError hermitian_space(E, FlintQQ[1 0; 1 1])

  V = @inferred hermitian_space(E, FlintQQ[1 2; 2 1])
  @test V === hermitian_space(E, FlintQQ[1 2; 2 1])
  @test V !== hermitian_space(E, FlintQQ[1 2; 2 1], cached = false)
  @test ishermitian(V)
  @test !isdefinite(V)
  @test involution(V) == s
  @test det(V) == -discriminant(V)

  U = hermitian_space(E, E[a 1; 1 1])
  V = hermitian_space(E, E[0 1; 1 0])
  W = hermitian_space(E, E[0 1 0; 1 0 0; 0 0 1])
  @test !Hecke.islocally_hyperbolic(U,p)
  @test !Hecke.isisotropic(U,p)
  @test Hecke.islocally_hyperbolic(V,p)
  @test Hecke.isisotropic(V,p)
  @test !Hecke.islocally_hyperbolic(W,p)
  @test Hecke.isisotropic(W,p)
  @test_throws AssertionError Hecke.islocally_hyperbolic(V, 2*OK)

  K, a = rationals_as_number_field()
  OK = maximal_order(K)
  Kt, t = K["t"]

  E, b = NumberField(t^2+17)

  p = 2*OK
  q = 17*OK

  V = hermitian_space(E, E[102 b; -b 0])
  H = hermitian_space(E, E[0 1; 1 0])
  W = hermitian_space(E, E[1 1 2; 1 2 3; 2 3 1])
  for r in [p,q]
    @test Hecke.islocally_represented_by(V,H,r)
  end
  @test Hecke.isrepresented_by(V,H)
  @test !Hecke.islocally_represented_by(W,V,p)
  @test !Hecke.isrepresented_by(W,V)

end


text:
module problem104

description = """
Pandigital Fibonacci ends
Problem 104

The Fibonacci sequence is defined by the recurrence relation:

    Fn = Fn−1 + Fn−2, where F1 = 1 and F2 = 1.

It turns out that F541, which contains 113 digits, is the first Fibonacci number for which the last nine digits are 1-9 pandigital (contain all the digits 1 to 9, but not necessarily in order). And F2749, which contains 575 digits, is the first Fibonacci number for which the first nine digits are 1-9 pandigital.

Given that Fk is the first Fibonacci number for which the first nine digits AND the last nine digits are 1-9 pandigital, find k.
"""

# test if the string representation of a number starts with digits 1-9
is_start_pandigital(x :: String) = (length(x) >= 9) && all([in(c, x[1:9]) for c = '1':'9'])

# test if an integer ends with digits 1-9
# This has been optimized to avoid memory allocations as it is the inner loop
function is_end_pandigital(x :: Int, digits=zeros(Bool, 10))
  digits[:] = false
  for i = 1:9
    digits[x % 10  + 1] = true
    x = div(x, 10)
  end
  for i = 2:10
    if !digits[i]
      return false
    end
  end
  return true
end

# Find the first fibonacci that starts with and ends with digits 1-9
# Calculate fibonacci sequence in BigInts, and in Int64s (modulo 10^10)
# Test for end-pandigital on the native ints, and when found, test the BigInts for start-pandigital
function find_start_end_pandigital()
  a, b = 1, 1
  A, B = BigInt(a), BigInt(b)
  buf = zeros(Bool, 10)
  for k = 3:1000000
    a, b = b, (a+b) % 10000000000
    A, B = B, (A+B)
    if is_end_pandigital(b, buf)
      if is_start_pandigital(string(B))
        return k
      end
    end
  end
end

using Base.Test

@test is_start_pandigital("3219485766456")
@test !is_start_pandigital("10201049586723")
@test is_end_pandigital(958928534791826)

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
"""
like bloom_best2, but with continue in main loop instead of elseif
no significant change (as expected)

"""
function bloom_best end

function bloom_bestc(t::SearchSequence)
    n = sizeof(t)
    skip = n
    tlast = _nthbyte(t,n)
    bloom_mask = UInt64(_search_bloom_mask(tlast))
    bloom_bits = 1 # no. of bits set in bloom filter if <= 32
    bloom_skip = 1 # no. of bytes to skip if byte not in filter
    j = n
    while (j-=1)>=1  # j=n is already processed
        tj = _nthbyte(t,j)
        if tj == tlast && skip == n
            skip = n - j
        end
        if bloom_bits <= 32 # argument: is near max(p(bloom_skip)*bloom_skip)
            hash = _search_bloom_mask(tj)
            if hash&bloom_mask == 0
                if bloom_bits < 32
                    # put in bloom filter up to 32 bits
                    bloom_mask |= hash
                end
                bloom_bits += 1 # gets >32 to terminate bloom filter filling
            end
            if bloom_bits <= 32
                bloom_skip += 1
            end
        else
            if skip<n
                # we have finished adding hashes to bloom filter
                # and we have determined the skip distance if matching last byte
                # nothimg remains to be done in preprocessing - stop work.
                break
            end
        end
    end
    return t,bloom_mask,bloom_skip,bloom_bits,skip,tlast
end

function bloom_bestc(s::SearchSequence, t::SearchSequence, i::Integer)
    bloom_bestc(s,bloom_bestc(t),i)
end

function bloom_bestc(s::SearchSequence,p::Tuple,i::Integer,sv::MaybeVector=nothing)
    (t,bloom_mask,bloom_skip,bloom_bits,skip,tlast) = p
    DOSTATS = !(sv isa Nothing)
    n = sizeof(t)
    m = sizeof(s)
    if n == 0
        return 1 <= i <= m+1 ? max(1, i) : 0
    elseif m == 0
        return 0
    elseif n == 1
        return something(findnext(isequal(_nthbyte(t,1)), s, i), 0)
    end
    w = m - n
    if w < 0 || i - 1 > w
        return 0
    end

    tlast = _nthbyte(t,n)
    if DOSTATS loops = 0 end
    if DOSTATS bloomtests = 0 end
    if DOSTATS bloomskips = 0 end
    i +=n-1
    if bloom_bits <= 12 # best guess from benchmarks
        # do bloom test first in loop
        while i < m
            if DOSTATS loops += 1 end
            if DOSTATS bloomtests += 1 end
            if bloom_mask & _search_bloom_mask(_nthbyte(s,i)) == 0
                if DOSTATS bloomskips += 1 end
                i += bloom_skip
                continue
            end
            if _nthbyte(s,i) == tlast
                # check candidate
                j = 1
                while j < n
                    if _nthbyte(s,i-n+j) != _nthbyte(t,j)
                        break
                    end
                    j += 1
                    # match found?
                    if j == n
                        if DOSTATS sv[Int(SFloops)] = loops; sv[Int(SFtests)] = bloomtests; sv[Int(SFskips)] = bloomskips; sv[Int(SFbits)] = bitcount(bloom_mask) end
                        return i-n+1
                    end
                end
                # no match: skip and test bloom
                i += skip
                continue
            end
            i +=1
        end
    else
        # do byte test first
        while i < m
            if DOSTATS loops += 1 end
            if _nthbyte(s,i) == tlast
                # check candidate
                j = 1
                while j < n
                    if _nthbyte(s,i-n+j) != _nthbyte(t,j)
                        break
                    end
                    j += 1
                    # match found?
                    if j == n
                        if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
                        return i-n+1
                    end
                end
                # no match: skip and test bloom
                i += skip
                if DOSTATS bloomtests += 1 end
                if i>=w
                    break
                end
                continue
            end
            i += 1
            if DOSTATS bloomtests += 1 end
            if bloom_mask & _search_bloom_mask(_nthbyte(s,i)) == 0
                if DOSTATS bloomskips += 1 end
                i += bloom_skip
            end
        end
    end
    if i==m
        # test end match
        j = 1
        while j <= n
            if _nthbyte(s,i-n+j) != _nthbyte(t,j)
                break # not found
            end
            if j == n
                if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
                return i-n+1
            end # match at the very end
            j += 1
        end
    end
    if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
    0
end


text:
# Autogenerated wrapper script for libevdev_jll for powerpc64le-linux-gnu
export libevdev

JLLWrappers.@generate_wrapper_header("libevdev")
JLLWrappers.@declare_library_product(libevdev, "libevdev.so.2")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libevdev,
        "lib/libevdev.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Documenter, DiffEqParamEstim

makedocs(
    sitename="DiffEqParamEstim.jl",
    authors="Chris Rackauckas et al.",
    clean=true,
    doctest=false,
    modules=[DiffEqParamEstim],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://diffeqparamestim.sciml.ai/stable/"),

    pages=[
        "DiffEqParamEstim.jl: Parameter Estimation for Differential Equations" => "index.md",
        "Tutorials" => Any[
          "tutorials/ODE_inference.md",
          "tutorials/global_optimization.md",
          "tutorials/jump.md",
          "tutorials/generalized_likelihood.md",
          "tutorials/stochastic_evaluations.md"
        ],
        "Methods" => Any[
          "methods/recommended_methods.md",
          "methods/optimization_based_methods.md",
          "methods/collocation_loss.md",
          "methods/alternative_objectives.md"
        ],
    ]
)

deploydocs(
    repo="github.com/SciML/DiffEqParamEstim.jl";
    push_preview=true
)


text:
# -----------------------------------------------
# Generalized Linear Models routines
# -----------------------------------------------

# using Optim
# include("distribtools.jl")
# include("types.jl")

function glmregress{T<:Number}(X::Array{T}, Y::Array{T}, regtype; fast=false)
    if !(typeof(regtype)<:GLMRegressionType)
        error("glmregress() - Incorrect GLMRegressionType specified")
    else
        if isdefined(regtype, :intercept)
            X = regtype.intercept ? addintercept(X) : nointercept(X)
        end
        fitted_beta, residuals, converged, iter = regtype(X, Y)
        if regtype.robust == false || typeof(regtype.robust) == BasicVariance
            variance = BasicVariance(X, residuals, length(fitted_beta))
        elseif typeof(regtype.robust) == HCEVariance
            variance = HCEVariance(X, residuals)
        end
        if fast
            return (fitted_beta, variance)
        else
            return glmfitresult(regtype, X, Y, fitted_beta, residuals, variance, regtype.robust, converged, iter)
        end
    end
end

# #########   LOGIT   #########

function (logit_glmfit::Logit)(X::Array{T}, Y::Array{T}) where {T<:Number}
    theta0 = ones(size(X)[2])
    objective(theta) = -logLike(X, Y, theta, Logistic) + 0.05*norm(theta)
    optimresult = optimize(objective, theta0, NewtonTrustRegion())
    fitted_beta = Optim.minimizer(optimresult)
    residuals = Y - Alistair.logistic(aprod(X, fitted_beta))
    return fitted_beta, residuals, Optim.converged(optimresult), Optim.iterations(optimresult)
end


text:
# ------------------------------------------------------------------------------------------
# # 新手入门
#
# 主题：
# 1. 如何打印（print）
# 2. 如何给变量赋值
# 3. 如何注释（comment）
# 4. 基础数学语法
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## 如何打印（print）
#
# 在 Julia 中使用 `println()` 函数来打印输出（print）
# ------------------------------------------------------------------------------------------

println("I'm excited to learn Julia!")

# ------------------------------------------------------------------------------------------
# ## 如何给变量赋值
#
# 只需要变量名、变量值和一个等号！
#
# Julia 会自动判断类型。
# ------------------------------------------------------------------------------------------

my_answer = 42
typeof(my_answer)

my_pi = 3.14159
typeof(my_pi)

😺 = "开心的猫!"
typeof(😺)

# ------------------------------------------------------------------------------------------
# 要打出“开心的猫”，用 tab 键补全选择 emoji 的名字然后再按 tab 键即可
#
# > `\:smi` + `<tab>` ==> 上下键选定 + `<enter>` ==> `<tab>` + `<enter>` ==> 完成
# >
# > 即 `\:smile_cat:` + `<tab>`
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 在给一个变量赋值后，可以给它重新赋任意类型的值而不会报错。
# ------------------------------------------------------------------------------------------

😺 = 1

typeof(😺)

# ------------------------------------------------------------------------------------------
# 注意：Julia 允许用超广泛的字符，就比如😺
#
# 这就允许我们可以把代码写成这样：
# ------------------------------------------------------------------------------------------

😀 = 0
😞 = -1

😺 + 😞 == 😀

# ------------------------------------------------------------------------------------------
# ## 如何注释（comment）
# ------------------------------------------------------------------------------------------

# 单行注释用井号键（pound/hash key）

#=

多行注释，
用“#= =#”组合

=#

# ------------------------------------------------------------------------------------------
# ## 基础数学语法
# ------------------------------------------------------------------------------------------

sum = 3 + 7  # 加

difference = 10 - 3  # 减

product = 20 * 5  # 乘

quotient = 100 / 10  # 除

power = 10 ^ 2  # 乘方

modulus = 101 % 2  # 取余

# ------------------------------------------------------------------------------------------
# ## 练习
#
# ### 1.1
# 查看函数 `convert` 的帮助文档。
# ------------------------------------------------------------------------------------------

# 在此作答

# ------------------------------------------------------------------------------------------
# ### 1.2
# 给变量 `days` 赋值`365`，将 `days` 转为浮点数（float）并将其赋值给变量 `days_float`
# ------------------------------------------------------------------------------------------

days = 0
days_float = days

@assert days == 365
@assert days_float == 365.0

# ------------------------------------------------------------------------------------------
# ### 1.3
# 看看下面两种用法的区别
#
# ```julia
# convert(Int64, "1")
# ```
# 和
#
# ```julia
# parse(Int64, "1")
# ```
# ------------------------------------------------------------------------------------------

convert(Int64, "1")

parse(Int64, "1")

# ------------------------------------------------------------------------------------------
# 请在完成练习运行 `@assert` 块以检测结果是否正确，若无报错即为正确。
# ------------------------------------------------------------------------------------------

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
@testset "Test 1D_B" begin
    # x = [0.0, 1.0, 2.0]                        # function nodes
    # x_b = [0.5, 1.5, 2.5]                      # function nodes
    # u = x^2
    # u_b = x_b^2                                    # function values in nodes
    # u_lb = u_b - 1.
    # u_ub = u_b + 1.
    # s = [2.0]                                  # function first derivative nodes
    # v = [4.0]                                  # function first derivative values

    @testset "Test 1D_B_RK_H0 kernel" begin
        x = [4.]
        x_b = [1., 2., 3.]
        u = [10.0]
        u_ub = [ 1., 3., 3.]
        u_lb = [ 0., 2., 0.]
        points = [0., 1., 2., 3., 4.]

        spl = approximate(x, u, x_b, u_lb, u_ub, RK_H0(0.1), 100)
        res = evaluate(spl, points)
        res_rounded = round.(res; digits=3)
        @test res_rounded ≈ [0.967, 1.0, 2.0, 3.0, 10.0]
        @test issetequal(spl._active, [1, 3, -2])
        cond = estimate_cond(spl)
        @test cond ≈ 100.0
#####
        x = [0., 5., 10.]
        x_b = [1., 2., 3., 4., 6., 7., 8., 9.]
        u = [0.0, 10.0, 0.]
        u_ub = [0.1, 1., 1., 10., 10., 1., 1., 0.1]
        u_lb = [0., 0., 0., 0., 0., 0., 0., 0.]
        points = [0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]

        spl = approximate(x, u, x_b, u_lb, u_ub, RK_H0(0.1), 100)
        res = evaluate(spl, points)
        res_rounded = round.(res; digits=3)
        @test res_rounded ≈ [0.0, 0.1, 0.55, 1.0, 5.5, 10.0, 5.5, 1.0, 0.55, 0.1, -0.0]
        @test issetequal(spl._active, [8, 1, 6, 3])
        cond = estimate_cond(spl)
        @test cond ≈ 10000.0

        spl_inter = interpolate(x, u, RK_H0(0.1))
        res_inter = evaluate(spl_inter, points)
        res_inter_rounded = round.(res_inter; digits=3)
        @test res_inter_rounded ≈ [0.0, 1.999, 3.999, 5.998, 7.999, 10.0, 7.999, 5.998, 3.999, 1.999, 0.0]
        cond_inter = estimate_cond(spl_inter)
        @test cond_inter ≈ 100.0
##
        x_b = [1., 2., 3., 4., 6., 7., 8., 9.]
        u_ub = [Inf, 1., Inf, 10., 10., 1., 1., 0.1]
        u_lb = [0., 0., 10., 0., 0., 0., 0., -Inf]
        points = [0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]

        spl = approximate(x_b, u_lb, u_ub, RK_H0(0.1), 100)
        res = evaluate(spl, points)
        res_rounded = round.(res; digits=2)
        @test res_rounded ≈ [0.98, 0.99, 1.0, 10.0, 7.75, 5.5, 3.25, 1.0, 0.55, 0.1, 0.1]
        @test issetequal(spl._active, [-3, 2, 6, 8])
        cond = estimate_cond(spl)
        @test cond ≈ 1000.0

####
        # σ = evaluate(spl, x)                # evaluate spline in nodes
        # @test σ ≈ u                         # compare with exact function values in nodes
        #
        # spl = prepare(x, RK_H0(0.1))
        # spl = construct(spl, u)
        # σ = evaluate(spl, x)
        # @test σ ≈ u
        # #
        # # Check that we get close when evaluating near the nodes
        # p = x .+ 1e-4*randn(size(x))   # evaluation points near the nodes
        # f = p.^2                       # exact function values in evaluation points
        # σ = evaluate(spl, p)           # evaluate spline in evaluation points
        # # compare spline values with exact function values in evaluation point
        # @test all(isapprox.(σ, f, atol = 0.05))
        #
        # σ = evaluate_one(spl, p[3])
        # @test σ ≈ f[3] atol = 0.05
    end

     @testset "Test 1D_B_H1 kernel" begin
         x = [0., 5., 10.]
         x_b = [1., 2., 3., 4., 6., 7., 8., 9.]
         u = [0.0, 10.0, 0.]
         u_ub = [0.1, 1., 1., 10., 10., 1., 1., 0.1]
         u_lb = [0., 0., 0., 0., 0., 0., 0., 0.]
         points = [0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]

         spl = approximate(x, u, x_b, u_lb, u_ub, RK_H1(0.1), 100)
         res = evaluate(spl, points)
         res_rounded = round.(res; digits=3)
         @test res_rounded ≈ [0.0, 0.039, 0.0, 1.0, 6.321, 10.0, 6.321, 1.0, 0.0, 0.039, 0.0]
         @test issetequal(spl._active, [6, 3, -2, -7])
         cond = estimate_cond(spl)
         @test cond ≈ 1.0e7

         spl_inter = interpolate(x, u, RK_H1(0.1))
         res_inter = evaluate(spl_inter, points)
         res_inter_rounded = round.(res_inter; digits=3)
         @test res_inter_rounded ≈ [0.0, 2.913, 5.627, 7.885, 9.428, 10.0, 9.428, 7.885, 5.627, 2.913, -0.0]
         cond_inter = estimate_cond(spl_inter)
         @test cond_inter ≈ 100000.0

##
         x = [0., 5., 10.]
         x_b = [1., 2., 3., 4., 6., 7., 8., 9.]
         u = [0.0, 10.0, 0.]
         u_ub = [Inf, 1., Inf, 10., 10., 1., 1., 0.1]
         u_lb = [0., 0., 0., 0., 0., 0., 0., -Inf]
         points = [0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]

         spl = approximate(x, u, x_b, u_lb, u_ub, RK_H1(0.1), 100)
         res = evaluate(spl, points)
         res_rounded = round.(res; digits=3)
         @test res_rounded ≈ [-0.0, 0.0, 1.0, 4.821, 9.028, 10.0, 5.77, 1.0, -0.0, -0.009, 0.0]
         @test issetequal(spl._active, [6, 2, -7, -1])
         cond = estimate_cond(spl)
         @test cond ≈ 1.0e7

###
        x_b = [1., 2., 3., 4., 6., 7., 8., 9.]
        u_ub = [Inf, 1., Inf, 10., 10., 1., 1., 0.1]
        u_lb = [0., 0., 10., 0., 0., 0., 0., -Inf]
        points = [0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]

        spl = approximate(x_b, u_lb, u_ub, RK_H1(0.1), 100)
        res = evaluate(spl, points)
        res_rounded = round.(res; digits=3)
        @test res_rounded ≈[1.691, -0.0, 1.0, 10.0, 10.0, 6.622, 3.343, 1.0, 0.0, -0.69, -1.362]
        @test issetequal(spl._active, [-3, 2, 6, -7, -1, 4])
        cond = estimate_cond(spl)
        @test cond ≈ 1.0e6

    #     spl = interpolate(x, u, RK_H1(0.1))  # create spline
    #     cond = estimate_cond(spl)              # get estimation of the gram matrix condition number
    #     @test cond ≈ 1.0e5
    #     σ = evaluate(spl, x)           # evaluate spline in nodes
    #     @test σ ≈ u                    # compare with exact function values in nodes
    #
    #     d1_σ = evaluate_derivative(spl, 1.0)     # evaluate spline first derivative in the node
    #     @test d1_σ ≈ 2.0 atol = 0.05              # compare with exact function first derivative value in the node
    #
    #     # Check that we get close when evaluating near the nodes
    #     p = x .+ 1e-4*randn(size(x))   # evaluation points near the nodes
    #     f = p.^2                       # exact function values in evaluation points
    #     σ = evaluate(spl, p)                # evaluate spline in evaluation points
    #     # compare spline values with exact function values in evaluation point
    #     @test all(isapprox.(σ, f, atol = 1e-2))
    #
    #     ###
    #     spl = interpolate(x, u, s, v, RK_H1(0.1)) # create spline by function and
    #                                         # first derivative values in nodes
    #     cond = estimate_cond(spl)              # get estimation of the gram matrix condition number
    #     @test cond ≈ 1.0e5
    #     σ = evaluate(spl, x)                # evaluate spline in nodes
    #     @test σ ≈ u                    # compare with exact function values in nodes
    #
    #     # Check that we get close when evaluating near the nodes
    #     p = x .+ 1e-3*randn(size(x))   # evaluation points near the nodes
    #     f = p.^2                       # exact function values in evaluation points
    #     σ = evaluate(spl, p)                # evaluate spline in evaluation points
    #     # compare spline values with exact function values in evaluation point
    #     @test all(isapprox.(σ, f, atol = 1e-2))
    #
    #     spl = prepare(x, s, RK_H1(0.1))
    #     spl = construct(spl, u, v)
    #     σ = evaluate_one(spl, p[3])
    #     @test σ ≈ f[3] atol = 0.05
     end
    # #
    @testset "Test 1D_B_H2 kernel" begin
        x = [4.]
        x_b = [1., 2., 3.]
        u = [10.0]
        u_ub = [ 1., 3., 3.]
        u_lb = [ 0., 2., 0.]
        points = [0., 1., 2., 3., 4.]

        spl = approximate(x, u, x_b, u_lb, u_ub, RK_H2(0.1), 100)
        res = evaluate(spl, points)
        res_rounded = round.(res; digits=3)
        @test res_rounded ≈ [-2.263, 1.0, 2.0, 3.0, 10.0]
        @test issetequal(spl._active, [1, 3, -2])
        cond = estimate_cond(spl)
        @test cond ≈ 1.0e7
##
        x = [0., 5., 10.]
        u = [0.0, 10.0, 0.]
        x_b = [1., 2., 3., 4., 6., 7., 8., 9.]
        u_ub = [0.1, 1., 1., 10., 10., 1., 1., 0.1]
        u_lb = [0., 0., 0., 0., 0., 0., 0., 0.]
        points = [0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]

        spl = approximate(x, u, x_b, u_lb, u_ub, RK_H2(0.1), 100)
        res = evaluate(spl, points)
        res_rounded = round.(res; digits=3)
        @test res_rounded ≈ [-0.0, 0.1, -0.0, 1.0, 6.298, 10.0, 6.298, 1.0, -0.0, 0.1, 0.0]
        @test issetequal(spl._active, [1, 8, 6, 3, -2, -7])
        cond = estimate_cond(spl)
        @test cond ≈ 1.0e12
##
        x_b = [1., 2., 3., 4., 6., 7., 8., 9.]
        u_ub = [Inf, 1., Inf, 10., 10., 1., 1., 0.1]
        u_lb = [0., 0., 10., 0., 0., 0., 0., -Inf]
        points = [0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]

        spl = approximate(x_b, u_lb, u_ub, RK_H2(0.1), 100)
        res = evaluate(spl, points)
        res_rounded = round.(res; digits=2)
        @test res_rounded ≈ [14.75, 0.0, 1.0, 10.0, 10.0, 4.45, 0.83, 0.0, 0.14, 0.1, -0.35]
        @test issetequal(spl._active, [-3, 2, -1, 4, 8, -6])
        cond = estimate_cond(spl)
        @test cond ≈ 1.0e10

    #     spl = interpolate(x, u, RK_H2(0.1))  # create spline
    #     cond = estimate_cond(spl)              # get estimation of the gram matrix condition number
    #     @test cond ≈ 1.0e7
    #     σ = evaluate(spl, x)           # evaluate spline in nodes
    #     @test σ ≈ u                    # compare with exact function values in nodes
    #
    #     d1_σ = evaluate_derivative(spl, 1.0)     # evaluate spline first derivative in the node
    #     @test d1_σ ≈ 2.0 atol = 0.005              # compare with exact function first derivative value in the node
    #
    #     # Check that we get close when evaluating near the nodes
    #     p = x .+ 1e-4*randn(size(x))   # evaluation points near the nodes
    #     f = p.^2                       # exact function values in evaluation points
    #     σ = evaluate(spl, p)                # evaluate spline in evaluation points
    #     # compare spline values with exact function values in evaluation point
    #     @test all(isapprox.(σ, f, atol = 1e-2))
    #
    #     ###
    #     spl = interpolate(x, u, s, v, RK_H2(0.1)) # create spline by function and
    #                                               # first derivative values in nodes
    #     σ = evaluate(spl, x)                # evaluate spline in nodes
    #     @test σ ≈ u                    # compare with exact function values in nodes
    #
    #     # Check that we get close when evaluating near the nodes
    #     p = x .+ 1e-3*randn(size(x))   # evaluation points near the nodes
    #     f = p.^2                       # exact function values in evaluation points
    #     σ = evaluate(spl, p)                # evaluate spline in evaluation points
    #     # compare spline values with exact function values in evaluation point
    #     @test all(isapprox.(σ, f, atol = 1e-2))
    end
end


text:
using BasicNT
using Test

primesieve(23)

@testset "BasicNT.jl" begin
    primesieve(23)
end


text:
#=

@authors : Spencer Lyon, John Stachurski

=#

using Optim

"""
The approximate Bellman operator, which computes and returns the
updated value function Tw on the grid points.  An array to store
the new set of values Tw is optionally supplied (to avoid having to
allocate new arrays at each iteration).  If supplied, any existing data in
Tw will be overwritten.

#### Arguments

`w` : Vector
      The value of the input function on different grid points
`grid` : Vector
         The set of grid points
`β` : AbstractFloat
         The discount factor
`u` : Function
      The utility function
`f` : Function
      The production function
`shocks` : Vector
           An array of draws from the shock, for Monte Carlo integration (to
           compute expectations).
`Tw` : Vector, optional (default=similar(w))
       Array to write output values to
`compute_policy` : Bool, optional (default=false)
                   Whether or not to compute policy function

"""
function bellman_operator(w::Vector,
                          grid::Vector,
                          β::AbstractFloat,
                          u::Function,
                          f::Function,
                          shocks::Vector,
                          Tw::Vector = similar(w);
                          compute_policy::Bool = false)

    # === Apply linear interpolation to w === #
    w_func = LinInterp(grid, w)

    if compute_policy
        σ = similar(w)
    end

    # == set Tw[i] = max_c { u(c) + β E w(f(y  - c) z)} == #
    for (i, y) in enumerate(grid)
        objective(c) = - u(c) - β * mean(w_func.(f(y - c) .* shocks))
        res = optimize(objective, 1e-10, y)

        if compute_policy
            σ[i] = res.minimizer
        end
        Tw[i] = - res.minimum
    end

    if compute_policy
        return Tw, σ
    else
        return Tw
    end
end


text:
"""
    introduce_Sgenvariables!(OPFpbs)

Modify `mp` field of all `Scenario` in `OPFpbs` in order to create Sgen variables when the problem will be constructed.
"""
function introduce_Sgenvariables!(OPFpbs)
    for (scenario, OPFpb) in OPFpbs
        node_formulations = OPFpbs[scenario].mp.node_formulations

        for (bus, elems_form) in node_formulations
            for elem in keys(elems_form)
                if ismatch(r"Gen", elem)
                    if scenario != basecase_scenario_name()
                        elems_form[elem] = :GOCcoupling
                    else
                        elems_form[elem] = :NewVar
                    end
                elseif ismatch(r"Load", elem) || ismatch(r"Shunt", elem)
                    elems_form[elem] = :NbMinVar
                else
                    elems_form[elem] = :NewVar
                end
            end
        end
    end
    return OPFpbs
end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
function add_depot(path::String; default::Bool = true)
    if !isdir(path) || !isdir(joinpath(path, "registries"))
        error("Depot not found at: $path")
    end

    if path in Pkg3.depots()
        @warn "Depot already being used"
    else
        push!(Pkg3.depots(), path)
    end

    if default
        if DEFAULT_DEPOT[] != path
            @info "Setting default depot to $path"
            DEFAULT_DEPOT[] = path
        end
    end
end

"""
    create_depot(path; default = true)

Make a new
"""
function create_depot(path::String; default::Bool = true)
    regpath = joinpath(path, "registries")
    if isdir(path)
        if isdir(regpath)
            if !isempty(readdir(regpath))
                @warn "Depot already exists at: $path"
            end
        else
            mkdir(regpath)
        end
    else
        mkdir(path)
        mkdir(joinpath(path, "registries"))
    end

    add_depot(path; default = default)
end

function default_depot(path::String)
    if !isdir(path) || !isdir(joinpath(path, "registries"))
        error("Depot not found at: $path")
    end
    if !(path in Pkg3.depots())
        if DEFAULT_DEPOT[] != path
            @info "Setting default depot to $path"
            DEFAULT_DEPOT[] = path
        end
    end
end


text:
using FastCap

m = FastCap.Model()

sur = FastCap.add_rectangle!(m,0,0,1,1)
FastCap.add_physicalsurf!(m,"conductor_bot",sur)
@time newsur = FastCap.extrude!(m,sur,(0,0,0.5))
FastCap.add_physicalsurf!(m,"new_cond",newsur)

open(f->println(f,FastCap.print_geo_file(m)),"simple.geo","w")
#spawn(`gmsh simple.geo`)

max_stars_repo_name: jebej/FastCap

max_issues_repo_name: jebej/FastCap

max_forks_repo_name: jebej/FastCap


text:
# ------------ Reading parameters and choices ------------ #

read_parameter(ctx::K, addr::Address) where K <: BackpropagationContext = read_parameter(ctx, ctx.params, addr)
read_parameter(ctx::K, params::Store, addr::Address) where K <: BackpropagationContext = getindex(ctx.initial_params, addr)

Zygote.@adjoint function read_parameter(ctx, params, addr)
    ret = read_parameter(ctx, params, addr)
    fn = param_grad -> begin
        state_grad = nothing
        params_grad = Store(Dict{Address, Any}(addr => param_grad))
        (state_grad, params_grad, nothing)
    end
    return ret, fn
end

read_choice(ctx::K, addr::Address) where K <: BackpropagationContext = read_choice(ctx, ctx.choices, addr)
read_choice(ctx::K, choices, addr::Address) where K <: BackpropagationContext = getindex(ctx.call, addr)

Zygote.@adjoint function read_choice(ctx, call, addr)
    ret = read_choice(ctx, call, addr)
    fn = choice_grad -> begin
        state_grad = nothing
        choices_grad = Store(Dict{Address, Any}(addr => choice_grad))
        (state_grad, choices_grad, nothing)
    end
    return ret, fn
end

# ------------ Simulate function calls ------------ #

# Grads for learnable parameters.
simulate_parameter_pullback(sel, params, param_grads, cl::T, args...) where T <: CallSite = get_ret(cl)

# Grads for choices with differentiable logpdfs.
simulate_choice_pullback(fillables, params, choice_grads, choice_target, cl::T, args...) where T <: CallSite = get_ret(cl)

# ------------ get_choice_gradients ------------ #

function get_choice_gradients(cl::T, ret_grad...) where T <: CallSite
    choice_grads = Gradients()
    vals, arg_grads = accumulate_choice_gradients!(Empty(), Empty(), choice_grads, SelectAll(), cl, ret_grad...)
    return vals, arg_grads, choice_grads
end

function get_choice_gradients(ps::P, cl::T, ret_grad...) where {P <: AddressMap, T <: CallSite}
    choice_grads = Gradients()
    vals, arg_grads = accumulate_choice_gradients!(Empty(), ps, choice_grads, SelectAll(), cl, ret_grad...)
    return vals, arg_grads, choice_grads
end

function get_choice_gradients(fillables::S, ps::P, cl::T, ret_grad...) where {S <: AddressMap, P <: AddressMap, T <: CallSite}
    choice_grads = Gradients()
    vals, arg_grads = accumulate_choice_gradients!(fillables, ps, choice_grads, SelectAll(), cl, ret_grad...)
    return vals, arg_grads, choice_grads
end

function get_choice_gradients(sel::K, cl::T, ret_grad...) where {T <: CallSite, K <: Target}
    choice_grads = Gradients()
    vals, arg_grads = accumulate_choice_gradients!(Empty(), Empty(), choice_grads, sel, cl, ret_grad...)
    return vals, arg_grads, choice_grads
end

function get_choice_gradients(sel::K, ps::P, cl::T, ret_grad...) where {T <: CallSite, K <: Target, P <: AddressMap}
    choice_grads = Gradients()
    vals, arg_grads = accumulate_choice_gradients!(Empty(), ps, choice_grads, sel, cl, ret_grad...)
    return vals, arg_grads, choice_grads
end

function get_choice_gradients(sel::K, fillables::S, ps::P, cl::T, ret_grad...) where {T <: CallSite, K <: Target, S <: AddressMap, P <: AddressMap}
    choice_grads = Gradients()
    vals, arg_grads = accumulate_choice_gradients!(fillables, ps, choice_grads, sel, cl, ret_grad...)
    return vals, arg_grads, choice_grads
end

# ------------ get_learnable_gradients ------------ #

function get_learnable_gradients(ps::P, cl::DynamicCallSite, ret_grad...; scaler::Float64 = 1.0) where P <: AddressMap
    param_grads = Gradients()
    arg_grads = accumulate_learnable_gradients!(target(), ps, param_grads, cl, ret_grad...; scaler = scaler)
    return arg_grads, param_grads
end

function get_learnable_gradients(sel::K, ps::P, cl::DynamicCallSite, ret_grad...; scaler::Float64 = 1.0) where {K <: AddressMap, P <: AddressMap}
    param_grads = Gradients()
    arg_grads = accumulate_learnable_gradients!(sel, ps, param_grads, cl, ret_grad...; scaler = scaler)
    return arg_grads, param_grads
end

# Convenience utilities (used in implementations of accumulate_learnable_gradients! and accumulate_choice_gradients! for each type of call site).
function acc!(param_grads, ::Nothing, scaler) end
function acc!(param_grads, ps_grad, scaler)
    for (addr, grad) in ps_grad.params
        accumulate!(param_grads, addr, scaler .* grad)
    end
end

function filter_acc!(choice_grads, cl, grad_tr::Store, sel::K) where K <: Target
    _, choices = projection(cl, sel)
    for (k, v) in shallow_iterator(cl)
        k in sel && begin
            haskey(grad_tr, k) && accumulate!(choice_grads, k, grad_tr[k])
        end
    end
    return choices
end

@inline filter_acc!(choice_grads, cl, grad_tr, sel::K) where K <: Target = projection(cl, sel)[2]

# ------------ Choice sites ------------ #

@inline function (ctx::ParameterBackpropagateContext)(call::typeof(trace),
                                                      addr::T,
                                                      d::Distribution{K}) where {T <: Address, K}
    s = getindex(ctx.call, addr)
    increment!(ctx, logpdf(d, s))
    return s
end

@inline function (ctx::ChoiceBackpropagateContext)(call::typeof(trace),
                                                   addr::T,
                                                   d::Distribution{K}) where {T <: Address, K}
    haskey(ctx.target, addr) || return get_value(get_sub(ctx.call, addr))
    s = read_choice(ctx, addr)
    increment!(ctx, logpdf(d, s))
    return s
end

# ------------ Learnable ------------ #

@inline (ctx::ParameterBackpropagateContext)(fn::typeof(learnable), addr::Address) = read_parameter(ctx, addr)

@inline function (ctx::ChoiceBackpropagateContext)(fn::typeof(learnable), addr::Address)
    haskey(ctx.initial_params, addr) && return getindex(ctx.initial_params, addr)
    error("(learnable): parameter not provided at address $addr.")
end

# ------------ Fillable ------------ #

@inline function (ctx::ParameterBackpropagateContext)(fn::typeof(fillable), addr::Address)
    haskey(ctx.fillables, addr) && return getindex(ctx.fillables, addr)
    error("(fillable): parameter not provided at address $addr.")
end

@inline function (ctx::ChoiceBackpropagateContext)(fn::typeof(fillable), addr::Address)
    haskey(ctx.target, addr) && return getindex(ctx.target, addr)
    error("(fillable): parameter not provided at address $addr.")
end

# ------------ Call sites ------------ #

@inline function (ctx::ParameterBackpropagateContext)(c::typeof(trace), addr::A, call::Function, args...) where A <: Address
    param_grads = Gradients()
    ret = simulate_parameter_pullback(get_sub(ctx.fillables, addr),
                                      get_sub(ctx.initial_params, addr),
                                      param_grads,
                                      get_sub(ctx.call, addr),
                                      args)
    ctx.param_grads.tree[addr] = param_grads
    return ret
end

@inline function (ctx::ChoiceBackpropagateContext)(c::typeof(trace), addr::A, call::Function, args...) where A <: Address
    choice_grads = Gradients()
    ret = simulate_choice_pullback(get_sub(ctx.fillables, addr),
                                   get_sub(ctx.initial_params, addr),
                                   choice_grads,
                                   get_sub(ctx.target, addr),
                                   get_sub(ctx.call, addr),
                                   args)
    ctx.choice_grads.tree[addr] = choice_grads
    return ret
end

# ------------ Parameter gradients ------------ #

Zygote.@adjoint function simulate_parameter_pullback(sel,
                                                     params,
                                                     param_grads,
                                                     cl::DynamicCallSite,
                                                     args...)
    ret = simulate_parameter_pullback(sel, params, param_grads, cl, args...)
    fn = ret_grad -> begin
        arg_grads = accumulate_learnable_gradients!(sel, params, param_grads, cl, ret_grad...)
        (nothing, nothing, nothing, nothing, arg_grads)
    end
    return ret, fn
end

function accumulate_learnable_gradients!(sel, initial_params, param_grads, cl::DynamicCallSite, ret_grad...; scaler::Float64 = 1.0)
    fn = (args, params) -> begin
        ctx = ParameterBackpropagate(cl, sel, initial_params, params, param_grads)
        ret = ctx(cl.fn, args...)
        (ctx.weight, ret)
    end
    blank = Store()
    _, back = Zygote.pullback(fn, cl.args, blank)
    arg_grads, ps_grad = back((1.0, ret_grad...))
    acc!(param_grads, ps_grad, scaler)
    return arg_grads
end

# ------------ Choice gradients ------------ #

Zygote.@adjoint function simulate_choice_pullback(fillables,
                                                  params,
                                                  choice_grads,
                                                  choice_target,
                                                  cl::DynamicCallSite,
                                                  args...)
    ret = simulate_choice_pullback(fillables, params, choice_grads, choice_target, cl, args...)
    fn = ret_grad -> begin
        choice_vals, arg_grads = accumulate_choice_gradients!(fillables, params, choice_grads, choice_target, cl, ret_grad...)
        (nothing, nothing, nothing, nothing, (choice_vals, choice_grads), arg_grads)
    end
    return ret, fn
end

function accumulate_choice_gradients!(fillables::S, initial_params::P, choice_grads, choice_target::K, cl::DynamicCallSite, ret_grad...) where {S <: AddressMap, P <: AddressMap, K <: Target}
    fn = (args, choices) -> begin
        ctx = ChoiceBackpropagate(cl, fillables, initial_params, choices, choice_grads, choice_target)
        ret = ctx(cl.fn, args...)
        (ctx.weight, ret)
    end
    blank = Store()
    _, back = Zygote.pullback(fn, cl.args, blank)
    arg_grads, grad_ref = back((1.0, ret_grad...))
    choice_vals = filter_acc!(choice_grads, cl, grad_ref, choice_target)
    return choice_vals, arg_grads
end

# ------------ Documentation ------------ #

@doc(
"""
```julia
mutable struct ParameterBackpropagateContext{T <: Trace} <: BackpropagationContext
    tr::T
    weight::Float64
    initial_params::AddressMap
    params::Store
    param_grads::Gradients
end
```
`ParameterBackpropagateContext` is used to compute the gradients of parameters with respect to following objective:

Outer constructors:
```julia
ParameterBackpropagate(tr::T, params) where T <: Trace = ParameterBackpropagateContext(tr, 0.0, params, Gradients())
ParameterBackpropagate(tr::T, params, param_grads::Gradients) where {T <: Trace, K <: Target} = ParameterBackpropagateContext(tr, 0.0, params, param_grads)
```
""", ParameterBackpropagateContext)

@doc(
"""
```julia
mutable struct ChoiceBackpropagateContext{T <: Trace} <: BackpropagationContext
    tr::T
    weight::Float64
    initial_params::AddressMap
    params::Store
    param_grads::Gradients
end
```
`ChoiceBackpropagateContext` is used to compute the gradients of choices with respect to following objective:

Outer constructors:
```julia
ChoiceBackpropagate(tr::T, init_params, params, choice_grads) where {T <: Trace, K <: Target} = ChoiceBackpropagateContext(tr, 0.0, params, choice_grads, SelectAll())
ChoiceBackpropagate(tr::T, init_params, params, choice_grads, sel::K) where {T <: Trace, K <: Target} = ChoiceBackpropagateContext(tr, 0.0, params, choice_grads, sel)
```
""", ChoiceBackpropagateContext)

@doc(
"""
```julia
gradients = get_choice_gradients(params, cl::T, ret_grad) where T <: CallSite
gradients = get_choice_gradients(cl::T, ret_grad) where T <: CallSite
```

Returns a `Gradients` object which tracks the gradients with respect to the objective of random choices with differentiable `logpdf` in the program.
""", get_choice_gradients)

@doc(
"""
```julia
gradients = get_learnable_gradients(params, cl::T, ret_grad, scaler::Float64 = 1.0) where T <: CallSite
```

Returns a `Gradients` object which tracks the gradients of the objective with respect to parameters in the program.
""", get_learnable_gradients)


text:
using Main
using Base.Test

@test fib(10) == 55

max_stars_repo_name: vmchale/project-init

max_issues_repo_name: vmchale/pi-templates

max_forks_repo_name: vmchale/project-init


text:
using TrajectoryOptimization
using Altro
using RobotDynamics
using Attitude
using SatelliteDynamics
using StaticArrays
using LinearAlgebra
using MATLAB
const TO = TrajectoryOptimization
const RD = RobotDynamics

include(joinpath(dirname(dirname(@__FILE__)),"ks_functions.jl"))

struct KSopt <: TO.AbstractModel
    μ::Float64
    tscale::Float64
    dscale::Float64
    uscale::Float64
end
function RobotDynamics.dynamics(model::KSopt, x, u_dot)

    # unpack state
    p = SVector(x[1],x[2],x[3],x[4])
    p_prime = SVector(x[5],x[6],x[7],x[8])
    h = x[9]
    #r = x[10]
    #z = x[11]
    u = SVector(x[12],x[13],x[14])

    P_j2 = SVector(u[1],u[2],u[3])/model.uscale
    # Levi-Civita transformation matrix
    L = L_fx(p)

    # append 0 to end of P_j2 acceleration vector
    P_j2 = SVector{4}([P_j2;0]*model.tscale^2/model.dscale)

    p_dp = -(h/2)*p + ((dot(p,p))/(2))*L'*P_j2

    return SVector(p_prime[1],p_prime[2],p_prime[3],p_prime[4],
                   p_dp[1],p_dp[2],p_dp[3],p_dp[4],
                   -2*dot(p_prime,L'*P_j2),
                   2*p'*p_prime,
                   2*(p_prime[1]*p[3] + p[1]*p_prime[3] + p_prime[2]*p[4] + p[2]*p_prime[4]),
                   u_dot[1],u_dot[2],u_dot[3])

end

# scaling
dscale = 1e7 # m
tscale = 20000 # s
uscale = 10000.0
μ = 3.986004418e14 *(tscale^2)/(dscale^3)
model = KSopt(μ,tscale, dscale, uscale)
Base.size(::KSopt) = 14,3
n,m = size(model)

# initial conditions
r_eci = @SVector [6.578137000000001e6,0,0]
v_eci = @SVector [0,9111.210131355681,4642.393437617197]
R_0 = r_eci/dscale
V_0 = v_eci*tscale/dscale
p_0 = u_from_x(R_0)                          # u
p_prime_0 = uprime_from_xdot(V_0,p_0)     # du/ds
h_0 = μ/norm(R_0) - .5*norm(V_0)^2         # m^2/s^2
t_0 = 0                                    # seconds
x0 = SVector{n}([p_0; p_prime_0; h_0; norm(R_0); R_0[3];zeros(3)])

# time stuff
N = 1501
dt = 3e-2
tf = dt*(N-1)

# cost function
scaling_num = 10
r_Desired = 42e6/dscale#4.2#6.77814
R_scale = .01
Q = Diagonal(SVector{n}([0*ones(9);scaling_num;100;R_scale*ones(3)]))
R = 0.001*Diagonal(SVector(R_scale,R_scale,R_scale))
q = SVector{n}([zeros(9);-r_Desired*scaling_num;0;zeros(3)])
costfun = QuadraticCost(Q,R,q=q)
costfun_term = QuadraticCost(Q,R,q=q)
obj = Objective(costfun, costfun_term, N)

# constraints
cons = ConstraintList(n,m,N)
u_bnd = 10.0
maxThrust = NormConstraint(n, m, u_bnd, TO.SecondOrderCone(), 12:14)
TO.add_constraint!(cons, maxThrust, 1:N)

xf = randn(n)
xf[11] = 4.2
xf = SVector{n}(xf)
prob = Problem(model, obj,xf, tf, x0=x0, constraints=cons,integration=RD.RK4)
opts = SolverOptions(verbose=2,
                     show_summary = true,
                     iterations = 20000,
                     projected_newton = false,
                     constraint_tolerance = 1e-4, penalty_scaling = 2.0, penalty_initial = 1e-2)
solver = ALTROSolver(prob,opts)

solve!(solver)

# get time stuff
function get_time_transfer(solver,dt,tscale)
    """get the time in days for a run. This allows us to remove t from state"""
    X = states(solver)
    t = 0.0
    t_hist = zeros(length(X))
    for i = 1:(length(X)-1)

        p0 = X[i][1:4]
        R0 = dot(p0,p0)

        p1 = X[i+1][1:4]
        R1 = dot(p1,p1)

        t += (R0 + R1)*dt/2
        t_hist[i+1] = copy(t)
    end
    t_days = t*tscale/(24*3600)
    t_hist = t_hist*tscale/(24*3600)
    return t_days, t_hist
end

t_days,t_hist = get_time_transfer(solver,dt,tscale)

X = states(solver)
U = controls(solver)

r_eci_hist = dscale*mat_from_vec([x_from_u(X[i][1:4]) for i = 1:length(X)])
v_eci_hist = (dscale/tscale)*mat_from_vec([xdot_from_u(X[i][1:4],X[i][5:8]) for i = 1:length(X)])

mat"
figure
hold on
plot3($r_eci_hist(1,:),$r_eci_hist(2,:),$r_eci_hist(3,:))
%axis equal
hold off
"

Xm = mat_from_vec(X)
U_real = Xm[12:14,:]

mat"
figure
hold on
title('Thrust Vector')
plot($t_hist,$U_real')
ylabel('Thrust Acceleration (0.0001m/s^2)')
xlabel('Time (days)')
hold off
xlim([0,$t_hist(end)])
legend('u_x','u_y','u_z')
%saveas(gcf,'long_traj.png')
"

U_real_vec = [X[i][12:14] for i = 1:length(X)]
Unorm = [norm(U_real_vec[i]) for i = 1:length(U)]

mat"
figure
hold on
plot($Unorm)
hold off
"

println("Time in days:  $t_days")# using TrajectoryOptimization

# now let's get orbital elements
const SD = SatelliteDynamics

oe_hist = [SD.sCARTtoOSC([r_eci_hist[:,i];v_eci_hist[:,i]]) for i = 1:size(r_eci_hist,2)]

a_hist = [oe_hist[i][1] for i = 1:length(oe_hist)]
e_hist = [oe_hist[i][2] for i = 1:length(oe_hist)]
i_hist = [rad2deg(oe_hist[i][3]) for i = 1:length(oe_hist)]

mat"
figure
hold on
subplot(3,1,1)
hold on
plot($t_hist,$a_hist)
subplot(3,1,2)
hold on
plot($t_hist,$e_hist)
subplot(3,1,3)
hold on
plot($t_hist,$i_hist)
hold off
"

# using JLD2
# X = states(solver)
# U = controls(solver)
# @save "rate_control_transfers/30_v3_day_transfer.jld2" X U t_hist

max_stars_repo_name: RoboticExplorationLab/KSLowThrust

max_issues_repo_name: RoboticExplorationLab/KSLowThrust

max_forks_repo_name: RoboticExplorationLab/KSLowThrust


text:
# OpenCL.CmdQueue

type CmdQueue <: CLObject
    id::CL_command_queue

    function CmdQueue(q_id::CL_command_queue; retain=false)
        if retain
            @check api.clRetainCommandQueue(q_id)
        end
        q = new(q_id)
        finalizer(q, x -> begin
            retain || _deletecached!(q)
            if x.id != C_NULL
                @check api.clReleaseCommandQueue(x.id)
                x.id = C_NULL
            end
        end )
        return q
    end
end

Base.pointer(q::CmdQueue) = q.id

function Base.show(io::IO, q::CmdQueue)
    ptr_val = convert(UInt, Base.pointer(q))
    ptr_address = "0x$(hex(ptr_val, Sys.WORD_SIZE>>2))"
    print(io, "OpenCL.CmdQueue(@$ptr_address)")
end

Base.getindex(q::CmdQueue, qinfo::Symbol) = info(q, qinfo)

function CmdQueue(ctx::Context)
    devs  = devices(ctx)
    flags = cl_command_queue_properties(0)
    if isempty(devs)
        throw(ArgumentError("OpenCL.CmdQueue Context argument does not have any devices"))
    end
    return CmdQueue(ctx, first(devs), flags)
end

function CmdQueue(ctx::Context, prop::Symbol)
    devs = devices(ctx)
    if isempty(devs)
        throw(ArgumentError("OpenCL.CmdQueue Context argument does not have any devices"))
    end
    return CmdQueue(ctx, first(devs), prop)
end

function CmdQueue(ctx::Context, props::NTuple{2, Symbol})
    devs = devices(ctx)
    if isempty(devs)
        throw(ArgumentError("OpenCL.CmdQueue Context argument does not have any devices"))
    end
    return CmdQueue(ctx, first(devs), props)
end

function CmdQueue(ctx::Context, dev::Device)
    flags = cl_command_queue_properties(0)
    return CmdQueue(ctx, dev, flags)
end

function CmdQueue(ctx::Context, dev::Device, prop::Symbol)
    flags = cl_command_queue_properties(0)
    if prop == :out_of_order
        flags |= CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE
    elseif prop == :profile
        flags |= CL_QUEUE_PROFILING_ENABLE
    else
        throw(ArgumentError("Only :out_of_order and :profile flags are valid, recognized flag $prop"))
    end
    return CmdQueue(ctx, dev, flags)
end

function CmdQueue(ctx::Context, dev::Device, props::NTuple{2,Symbol})
    if !(:out_of_order in props && :profile in props)
        throw(ArgumentError("Only :out_of_order and :profile flags are vaid, unrecognized flags $props"))
    end
    flags = CL_QUEUE_PROFILING_ENABLE | CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE
    return CmdQueue(ctx, dev, flags)
end

function CmdQueue(ctx::Context, dev::Device, props::CL_command_queue_properties)
    err_code = Array(CL_int, 1)
    queue_id = api.clCreateCommandQueue(ctx.id, dev.id, props, err_code)
    if err_code[1] != CL_SUCCESS
        if queue_id != C_NULL
            api.clReleaseCommandQueue(queue_id)
        end
        throw(CLError(err_code[1]))
    end
    return CmdQueue(queue_id)
end

function flush(q::CmdQueue)
    @check api.clFlush(q.id)
    return q
end

function finish(q::CmdQueue)
    @check api.clFinish(q.id)
    return q
end

let
    context(q::CmdQueue) = begin
        ctx_id = Array(CL_context, 1)
        @check api.clGetCommandQueueInfo(q.id, CL_QUEUE_CONTEXT,
                                         sizeof(CL_context), ctx_id, C_NULL)
        Context(ctx_id[1], retain=true)
    end

    device(q::CmdQueue) = begin
        dev_id = Array(CL_device_id, 1)
        @check api.clGetCommandQueueInfo(q.id, CL_QUEUE_DEVICE,
                                         sizeof(CL_device_id), dev_id, C_NULL)
        Device(dev_id[1])
    end

    reference_count(q::CmdQueue) = begin
        ref_count = Array(CL_uint, 1)
        @check api.clGetCommandQueueInfo(q.id, CL_QUEUE_REFERENCE_COUNT,
                                         sizeof(CL_uint), ref_count, C_NULL)
        ref_count[1]
    end

    properties(q::CmdQueue) = begin
        props = Array(CL_command_queue_properties, 1)
        @check api.clGetCommandQueueInfo(q.id, CL_QUEUE_PROPERTIES,
                                         sizeof(CL_command_queue_properties),
                                         props, C_NULL)
        props[1]
    end

    const info_map = Dict{Symbol, Function}(
        :context => context,
        :device => device,
        :reference_count => reference_count,
        :properties => properties
    )

    function info(q::CmdQueue, qinfo::Symbol)
        try
            func = info_map[qinfo]
            func(q)
        catch err
            if isa(err, KeyError)
                throw(ArgumentError("OpenCL.CmdQueue has no info for: $qinfo"))
            else
                throw(err)
            end
        end
    end
end

context(queue::CmdQueue) = info(queue, :context)


text:
struct OptResults
    LB::Float64                         # Lower Bound on the objective
    UB::Float64                         # Upper Bound on the objective
    status::Symbol                      # Status
    solvetime::Float64                  # Solve time in seconds
    iterations::Union{Nothing,Int}      # Number of iterations (for local search)
    nodecount::Union{Nothing,Int}       # Node count of the Branch-and-Bound tree
    lazycount::Union{Nothing,Int}       # Number of lazy constraints added during solving

    function OptResults(lb::Float64, ub::Float64, status::Symbol, solvetime::Float64,
                        iterations::Union{Nothing,Int}, nodecount::Union{Nothing,Int}, lazycount::Union{Nothing,Int})
        if ~(lb <= ub)
            throw(ArgumentError("Lower Bound value must be less than or equal to Upper Bound."))
        end
        if (status == :Optimal) && (lb != ub)
            throw(ArgumentError("An optimal solution implies that the Lower Bound is equal to the Upper Bound."))
        end
        if solvetime < 0
            throw(DomainError("Solve time must be a non-negative number (in seconds)."))
        end
        if (iterations != nothing) && (iterations < 0)
            throw(DomainError("Number of iterations must be a positive integer number."))
        end
        if (nodecount != nothing) && (nodecount < 0)
            throw(DomainError("Number of nodes in B&B must be a non-negative integer number."))
        end
        if (lazycount != nothing) && (lazycount < 0)
            throw(DomainError("Number of lazy constraints added must be a non-negative integer number."))
        end

        return new(lb, ub, status, solvetime, iterations, nodecount, lazycount)
    end
end

function displayResults(opt_results::OptResults)
    """
    Method to display optimization results.

    Parameters
    ----------
    opt_results : OptResults
        Optimization results object.
    """
    print_string = string("\n",
          "--------------- Opt Results ---------------\n",
          "Obj. LB = ",round(opt_results.LB, digits=5),"\n",
          "Obj. UB = ",round(opt_results.UB, digits=5),"\n",
          "Status: $(opt_results.status)\n",
          "Solve time: ",round(opt_results.solvetime, digits=5)," seconds\n")
    if opt_results.iterations != nothing
        print_string = string(print_string, "Iterations: $(opt_results.iterations)\n")
    end
    if opt_results.nodecount != nothing
        print_string = string(print_string, "Number of nodes: $(opt_results.nodecount)\n")
    end
    if opt_results.lazycount != nothing
        print_string = string(print_string, "Number of lazy constraints: $(opt_results.lazycount)\n")
    end
    print(print_string)
end

function saveResults(opt_results::OptResults, path::String; oneline::Bool=false, header::Bool=false)
    throw("not implemented")
end


text:
using ..Atmos
using ..Atmos: MoistureModel

# Helpers to gather the thermodynamic variables across the DG grid
#

function vars_thermo(m::AtmosModel, FT)
    @vars begin
        T::FT
        θ_dry::FT
        θ_vir::FT
        e_int::FT
        h_tot::FT
        h_moi::FT

        moisture::vars_thermo(m.moisture, FT)
    end
end
vars_thermo(::MoistureModel, FT) = @vars()
function vars_thermo(m::EquilMoist, FT)
    @vars begin
        q_liq::FT
        q_ice::FT
        q_vap::FT
        θ_liq_ice::FT
    end
end
num_thermo(bl, FT) = varsize(vars_thermo(bl, FT))
thermo_vars(bl, array) = Vars{vars_thermo(bl, eltype(array))}(array)

# visitor function, to use with `@visitQ`
function compute_thermo!(atmos::AtmosModel, state, aux, thermo)
    e_tot = state.ρe / state.ρ
    ts = thermo_state(atmos, state, aux)
    e_int = internal_energy(ts)

    thermo.T = air_temperature(ts)
    thermo.θ_dry = dry_pottemp(ts)
    thermo.θ_vir = virtual_pottemp(ts)
    thermo.e_int = e_int

    # Moist and total henthalpy
    R_m = gas_constant_air(ts)
    thermo.h_tot = e_tot + R_m * thermo.T
    thermo.h_moi = e_int + R_m * thermo.T

    compute_thermo!(atmos.moisture, state, aux, ts, thermo)

    return nothing
end
function compute_thermo!(::MoistureModel, state, aux, ts, thermo)
    return nothing
end
function compute_thermo!(moist::EquilMoist, state, aux, ts, thermo)
    Phpart = PhasePartition(ts)

    thermo.moisture.q_liq = Phpart.liq
    thermo.moisture.q_ice = Phpart.ice
    thermo.moisture.q_vap = vapor_specific_humidity(ts)
    thermo.moisture.θ_liq_ice = liquid_ice_pottemp(ts)

    return nothing
end


text:
function ethmult2(a::Integer, b::Integer)
    A = [a]
    B = [b]
    while A[end] > 1
        push!(A, halve(A[end]))
        push!(B, double(B[end]))
    end
    return sum(B[map(!even, A)])
end

@show ethmult2(17, 34)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module enoc

include("analog.jl")
using .Analog

using LinearAlgebra
using Statistics
using Random
using Serialization

using Distributions
using NearestNeighbors
using PyCall

include("ssa.jl")
include("ens_forecast.jl")
include("ssa_varimax.jl")

using .SSA
using .ens_forecast
using .ssa_varimax

nanmean(x) = mean(filter(!isnan,x))
nanmean(x,y) = mapslices(nanmean,x,dims=y)

tsa = pyimport("statsmodels.tsa.api")
xskillscore = pyimport("xskillscore")
xarray = pyimport("xarray")

function optimal_ens(info)
    D, m, N = size(info.ens)
    errs = [[sqrt(mean((nanmean(info.ens[:, (sortperm(info.r_errs[j, :]))[1:i], j], 2) - info.x_trues[:, j]).^2)) for j=1:N] for i=1:m]
    errs_rand = [[sqrt(mean((nanmean(info.ens[:, shuffle(sortperm(info.r_errs[j, :]))[1:i], j], 2) - info.x_trues[:, j]).^2)) for j=1:N] for i=1:m]
    return mean(hcat(errs...), dims=1)', mean(hcat(errs_rand...), dims=1)'
end

function optimal_ens_crps(info)
    D, m, N = size(info.ens)
    errs = zeros(m, N)
    for mp=1:m
        for j=1:N
            E_mp = info.ens[:, (sortperm(info.r_errs[j, :]))[1:mp], j]
            E_mp_array = xarray.DataArray(data=E_mp, dims=["dim", "member"])
            errs[mp, j] = xskillscore.crps_ensemble(info.x_trues[:, j], E_corr_array).values[1]
        end
    end
    return errs
end

function setup(; model, Δt, outfreq, obs_err_pct, M, record_length, transient,
               y0, D, osc_vars, modes, integrator, pcs=nothing, varimax, da,
               window, k, k_r)
    y = integrator(model, y0, 0., record_length*outfreq*Δt, Δt; inplace=false)[1:outfreq:end, :][(transient + 1):end, :]
    if (obs_err_pct > 0)
        R = Symmetric(diagm(0 => obs_err_pct*std(y, dims=1)[1, :]))
        obs_err = MvNormal(zeros(D), R)
        y = y + rand(obs_err, size(y)[1])'
    end

    ssa_info = ssa_decompose(y[:, osc_vars], M)

    if varimax
        ssa_info.eig_vals, ssa_info.eig_vecs = varimax_rotate!(ssa_info.eig_vals, ssa_info.eig_vecs, M, ssa_info.D, maximum(modes))
    end

    r = ssa_reconstruct(ssa_info, modes, sum_modes=true)

    if pcs === nothing
        tree = KDTree(copy(y'))
        tree_r = KDTree(copy(r'))
    else
        _, _, v = svd(r)
        v = v[:, 1:pcs]
        tree = KDTree(copy(y'))
    end

    return tree, tree_r, ssa_info, y, r
end

function run(; model, model_err, integrator, m, M, D, k, k_r, modes, osc_vars,
             outfreq, Δt, cycles, window, record_length, obs_err_pct,
             ens_err_pct, transient, y0, mp, varimax,
             check_bounds, test_time=false, da, inflation, y_fcst, α=nothing,
             preload=nothing)

    if (preload === nothing) | !isfile(preload)
        tree, tree_r, ssa_info, y, r = setup(model=model, record_length=record_length,
                                             integrator=integrator, Δt=Δt, y0=y0,
                                             transient=transient, outfreq=outfreq,
                                             obs_err_pct=obs_err_pct, osc_vars=osc_vars,
                                             D=D, M=M, modes=modes, varimax=varimax,
                                             da=da, window=window, k=k, k_r=k_r)

        if y_fcst
            var_model = tsa.VAR(y).fit(maxlags=15, ic="aic")
        else
            var_model = nothing
        end

        serialize(preload, (tree, tree_r, ssa_info, y, r, var_model))
    else
        tree, tree_r, ssa_info, y, r, var_model = deserialize(preload)
    end

    if da
        R = error_cov(y, r, M, window, k, k_r, osc_vars)
    else
        R = false
    end

    E = integrator(model_err, y[end, :], 0.0, m*Δt*outfreq, Δt, inplace=false)[1:outfreq:end, :]'

    stds = std(y, dims=1)[:]
    means = mean(y, dims=1)
    bounds = (minimum(y, dims=1)[:], maximum(y, dims=1)[:])

    ens_info = forecast(E=copy(E), model=model, model_err=model_err,
                        integrator=integrator, m=m, Δt=Δt, window=window,
                        cycles=cycles, outfreq=outfreq, D=D, k=k, k_r=k_r, r=r,
                        tree=tree, tree_r=tree_r, osc_vars=osc_vars, means=means,
                        stds=stds, err_pct=ens_err_pct, mp=mp,
                        check_bounds=check_bounds, test_time=test_time,
                        bounds=bounds, da=da, R=R, inflation=inflation,
                        var_model=var_model, α=α)

    return ens_info, ssa_info
end
end


text:
@doc """ copy build file from destine dir to source dir
""" ->
function makebuild(dst, src)
    build_fl = joinpath(src, "BUILD")
    @assert isfile(build_fl)

    data = readlines(build_fl)

    srcname = basename(src)
    dstname = basename(dst)

    data = convert(Array{ASCIIString,1}, data)

    #rstr = Regex("name = \"$(srcname)\"")

    data = map(line->replace(line, srcname, dstname), data)

    dstbuild = joinpath(dst, "BUILD")

    write(dstbuild, data)
end


text:
using Random
using Plots

function simulate_group(number_of_people::Int, group::Array, number_of_simulations::Int)
    local numerator = 0.0
    for i in 1:number_of_simulations::Int
        for k in 1:number_of_people
            group[k] = rand(1:365)
        end
        #local number_of_turns = rand(1:365, number_of_people)
        if length(unique(group)) < number_of_people
            numerator += 1.0
        end
    end
    return(numerator)
end

number_of_simulations = 4000

probability_shared_birthday_results = Array{Float32}(undef, 499)
time_per_simulation = Array{Float32}(undef, 499)

for num_people in 1:499
    local group = Array{Int16}(undef, num_people+1)
    stats = @timed numer = simulate_group(num_people+1, group, number_of_simulations)
    probability_shared_birthday_results[num_people] = numer/number_of_simulations
    time_per_simulation[num_people] = stats.time
end

plot([2:1:500], probability_shared_birthday_results, xlim=[2, 500], dpi=600, xlabel="Number of People", ylabel="Probability of Shared Birthdays", legend = false)
png("birthdayproblem.png")

plot([2:1:500], time_per_simulation, xlim=[2, 500], dpi=600, xlabel="Number of People", ylabel="Simulation Time (s)", legend = false)
png("birthdayproblemtime.png")


text:
using Robotlib
using Base.Test
import Robotlib: ad, adi

function simulateCalibration1(N)
    n       = 6
    q       = 2π*rand(N,n)

    wn      = zeros(3,n)
    pn      = zeros(3,n)
    xin     = zeros(6,n+1)
    wn[:,1] = [0,0,1]
    pn[:,1] = [0,0,0]
    wn[:,2] = [0,1,0]
    pn[:,2] = [0,0,0.3]
    wn[:,3] = [0,1,0]
    pn[:,3] = [0,0,1]
    wn[:,4] = [1,0,0]
    pn[:,4] = [1,0,1]
    wn[:,5] = [0,1,0]
    pn[:,5] = [1,0,1]
    wn[:,6] = [1,0,0]
    pn[:,6] = [1,0,1]
    T0 = [  0 -1 0  1.2;
    0 0 -1 0;
    1 0  0  1;
    0 0  0  1]
    for i = 1:n
        xin[:,i] = [-skew(pn[:,i])*wn[:,i]; wn[:,i]]
    end
    xin[:,n+1] = twistcoords(logT(T0))
    xinmod   = deepcopy(xin)
    for i = 1:n+1
        xinmod[:,i]  += [0.001randn(3); 0.01π/180*randn(3)]
        #conformize(xinmod[:,i])
    end
    Ta  = zeros(4,4,N)
    for i = 1:N
        Ta[:,:,i] = fkinePOE(xin,q[i,:])
    end
    return q, xin, T0, xinmod, Ta
end

function simulateCalibration2(N)
    n       = 6
    q       = 2π*rand(N,n)
    dh      = DH7600()
    AAA,BBB,T0,Ti0,Tn0 = jacobian(zeros(6),dh, eye(4));
    xin     = DH2twistsPOE(dh)
    xinmod   = deepcopy(xin)
    for i = 1:n+1
        xinmod[:,i]  += [0.001randn(3); 0.01π/180*randn(3)]
        #conformize(xinmod[:,i])
    end
    Ta  = zeros(4,4,N)
    for i = 1:N
        Ta[:,:,i] = fkinePOE(xin,q[i,:])
    end
    return q, xin, T0, xinmod, Ta
end

function simulateCalibrationLPOE(N)
    srand(1)
    n       = 6
    q       = 2π*rand(N,n)
    dh      = DH7600()
    AAA,BBB,T0,Ti0,Tn0 = jacobian(zeros(6),dh, eye(4));
    xin = DH2twistsLPOE(Tn0)
    Tn0mod = deepcopy(Tn0)
    for i = 1:n+1
        Tn0mod[1:3,4,i] += 0.1randn(3)
    end
    Ta  = zeros(4,4,N)
    for i = 1:N
        Ta[:,:,i] = fkineLPOE(Tn0,xin,q[i,:])
        AAA,BBB,T = jacobian(q[i,:],dh, eye(4));
        @assert T ≈ Ta[:,:,i]
    end
    return q, xin, Tn0, Tn0mod, Ta
end

function simulateCalibration4(N)
    #srand(1)
    n       = 6
    qt      = 2π*rand(N,n)
    q1      = qt + 0.01*π/180*randn(N,n)
    q2      = qt + 0.01*π/180*randn(N,n)
    q       = cat(3,q1,q2)
    dh      = DH7600()
    AAA,BBB,T0,Ti0,Tn0 = jacobian(zeros(6),dh, eye(4));
    xin = DH2twistsLPOE(Tn0)
    xin = cat(3,xin,xin)
    Tn0mod = cat(4,deepcopy(Tn0),deepcopy(Tn0))
    for i = 1:n+1
        Tn0mod[1:3,4,i,1] += 0.1randn(3)
        Tn0mod[1:3,4,i,2] += 0.1randn(3)
    end
    Ta  = zeros(4,4,N)
    for i = 1:N
        Ta[:,:,i] = fkineLPOE(Tn0,xin,qt[i,:]')
    end

    return q, xin, Tn0, Tn0mod, Ta
end

function run_tests()
    # q, xin,T0, xinmod, Ta= simulateCalibration2(100)
    # xic,et,er = Robotlib.Calibration.calibPOE(xinmod,Ta,q,maxiter=150, λ = 10000.0)
    # display(norm(xin[:,1:7]-xinmod[:,1:7]))
    # display(norm(xin[:,1:7]-xic[:,1:7]))
    N = 100
    println("===== Testing calibLPOE =====")
    q, xin,Tn0,Tn0mod, Ta = simulateCalibrationLPOE(N)
    @time Tn0c,xic,et,er = Robotlib.Calibration.calibLPOE(xin,Tn0mod,Ta,q,maxiter=10, λ = 100.0)
    @test et[et .!= 0][end] < 1e-12
    @test er[er .!= 0][end] < 1e-12
    println("Initial error: ", sqrt(sum((Tn0-Tn0mod).^2)))
    println("Final error: ", sqrt(sum((Tn0-Tn0c).^2)))
    @test sqrt(sum((Tn0-Tn0c).^2)) < 0.8sqrt(sum((Tn0-Tn0mod).^2))

    println("===== Testing calibLPOEdual =====")
    q, xin,Tn0,Tn0mod, Ta = simulateCalibration4(N)
    @time Tn0c,xic,et,er = Robotlib.Calibration.calibLPOEdual(xin,Tn0mod,q,maxiter=6, λ = 0.01)
    println("Error between Tn0 and Tn0mod: ",
        sqrt(sum((Tn0[1:3,4,:]-Tn0mod[1:3,4,:,1]).^2+(Tn0[1:3,4,:]-Tn0mod[1:3,4,:,2]).^2)/N))
    println("Error between Tn0 and Tn0c  : ",
        sqrt(sum((Tn0[1:3,4,:]-Tn0c[1:3,4,:,1]).^2+(Tn0[1:3,4,:]-Tn0c[1:3,4,:,2]).^2)/N))
    ei = 0.0
    ec = 0.0
    for i = 1:N
        T1 = fkineLPOE(Tn0mod,xin,q[i,:,1]')
        T2 = fkineLPOE(Tn0mod,xin,q[i,:,2]')
        ei += norm(twistcoords(logm(Ta[:,:,i]*trinv(T1))))
        ei += norm(twistcoords(logm(Ta[:,:,i]*trinv(T2))))
        T1 = fkineLPOE(Tn0c,xin,q[i,:,1]')
        T2 = fkineLPOE(Tn0c,xin,q[i,:,2]')
        ec += norm(twistcoords(logm(Ta[:,:,i]*trinv(T1))))
        ec += norm(twistcoords(logm(Ta[:,:,i]*trinv(T2))))
    end
    println("Initial error: ",round(ei/N,5), " Calibrated error: ", round(ec/N,5))
    @test ec < ei
end

run_tests()

println("===== Testing frames =====")
include("testFrames.jl")


text:
"""
    shiftedkernel = centered(kernel)

Shift the origin-of-coordinates to the center of `kernel`. The
center-element of `kernel` will be accessed by `shiftedkernel[0, 0,
...]`.

This function makes it easy to supply kernels using regular Arrays,
and provides compatibility with other languages that do not support
arbitrary axes.

See also: [`imfilter`](@ref).
"""
function centered(A::AbstractArray)
    offsetted = first.(axes(A)) .- 1
    total_offsets = .-((size(A) .+ 1 ) .÷ 2)
    OffsetArray(A, total_offsets .- offsetted)
end

"""
    kernfft = freqkernel([T::Type], kern, sz=size(kern); rfft=false)

Return a frequency-space representation of `kern`.
This embeds `kern` in an array of size `sz`, in a manner that implicitly imposes periodic boundary
conditions, and then returns the fourier transform. This is sometimes called the optical transfer
function, and known in some frameworks as `psf2otf`. If `rfft` is `true`, the fft for real-valued
arrays (`rfft`) is returned instead and the first dimension size will be approximately half of `sz[1]`.

`kern` should be zero-centered, i.e., `kern[0, 0]` should reference the center of your kernel.
See [`centered`](@ref). Optionally specify the numeric type `T` (which must be one of the
types supported by FFTW, either `Float32` or `Float64`).

The inverse of `freqkernel` is [`spacekernel`](@ref).
"""
function freqkernel(::Type{T}, kern::AbstractArray, sz::Dims=size(kern); rfft=false) where T<:Union{Float32,Float64}
    wrapindex(i, s) = i<1 ? i+s : i
    all(size(kern) .<= sz) || throw(DimensionMismatch("kernel size $(size(kern)) is bigger than supplied size $sz"))
    kernw = zeros(T, sz...)
    for I in CartesianIndices(kern)
        J = CartesianIndex(map(wrapindex, Tuple(I), sz))
        kernw[J] = kern[I]
    end
    return rfft ? FFTW.rfft(kernw) : fft(kernw)
end
freqkernel(kern::AbstractArray{T}, args...; rfft=false) where T =
    freqkernel(ffteltype(T), kern, args...; rfft=rfft)

"""
    kern = spacekernel(kernfft, axs; rfftsz=0)

Return a real-space representation of `kernfft`, a frequency-space representation of a kernel.
This performs an inverse fourier transform, implicitly imposes periodic boundary conditions,
and then trims & truncates axes of the output to `axs`. By default `kernfft` is assumed to have
been generated by `fft`; if it was instead generated by `rfft`, specify the original size
of the first dimension. (If `kernfft` was generated by [`freqkernel`](@ref), this is just `sz[1]`.)

The inverse of `spacekernel` is [`freqkernel`](@ref).
"""
function spacekernel(kernfft::AbstractArray, axs::Indices; rfftsz=0)
    wrapindex(i, s) = i<1 ? i+s : i
    kernw = rfftsz > 0 ? irfft(kernfft, rfftsz) : ifft(kernfft)
    kern = zeros(eltype(kernw), axs...)
    sz = size(kernw)
    for I in CartesianIndices(kern)
        J = CartesianIndex(map(wrapindex, Tuple(I), sz))
        kern[I] = kernw[J]
    end
    return kern
end
function spacekernel(kernfft::AbstractArray; rfftsz=0)
    sz = size(kernfft)
    if rfftsz > 0
        sz = (rfftsz, Base.tail(sz)...)
    end
    upper = map(s->s>>1, sz)
    axs = map((u, s)-> u-s+1:u, upper, sz)
    return spacekernel(kernfft, axs; rfftsz=rfftsz)
end

ffteltype(::Type{T}) where T<:Union{Float32,Float64} = T
ffteltype(::Type{Float16}) = Float32
ffteltype(::Type{T}) where T<:Normed = ffteltype(floattype(T))
ffteltype(::Type{T}) where T = Float64

dummyind(::Base.OneTo) = Base.OneTo(1)
dummyind(::AbstractUnitRange) = 0:0

dummykernel(inds::Indices{N}) where {N} = fill(1, map(dummyind, inds))

nextendeddims(inds::Indices) = sum(ind->length(ind)>1, inds)
nextendeddims(a::AbstractArray) = nextendeddims(axes(a))

function checkextended(inds::Indices, n)
    dimstr = n == 1 ? "dimension" : "dimensions"
    nextendeddims(inds) != n && throw(ArgumentError("need $n extended $dimstr, got axes $inds"))
    nothing
end
checkextended(a::AbstractArray, n) = checkextended(axes(a), n)

_reshape(A::OffsetArray{<:Any,N}, ::Val{N}) where N = A
_reshape(A::OffsetArray, ::Val{N}) where {N} = OffsetArray(reshape(parent(A), Val(N)), fill_to_length(A.offsets, -1, Val(N)))
_reshape(A::AbstractArray, ::Val{N}) where {N} = reshape(A, Val(N))

_vec(a::AbstractVector) = a
_vec(a::AbstractArray) = (checkextended(a, 1); a)
_vec(a::OffsetArray{<:Any,1}) = a
function _vec(a::OffsetArray)
    inds = axes(a)
    checkextended(inds, 1)
    i = findall(ind->length(ind)>1, inds)
    OffsetArray(vec(parent(a)), inds[i])
end

samedims(::Val{N}, kernel) where {N} = _reshape(kernel, Val(N))
samedims(::Val{N}, kernel::Tuple) where {N} = map(k->_reshape(k, Val(N)), kernel)
samedims(::AbstractArray{<:Any,N}, kernel) where {N} = samedims(Val(N), kernel)

_tail(R::CartesianIndices{0}) = R
_tail(R::CartesianIndices) = CartesianIndices(tail(axes(R)))

# ensure that overflow is detected, by ensuring that it doesn't happen
# at intermediate stages of the computation
accumfilter(pixelval, filterval) = pixelval * filterval
const SmallInts = Union{UInt8,Int8,UInt16,Int16}
accumfilter(pixelval::SmallInts, filterval::SmallInts) = Int(pixelval)*Int(filterval)
# advice: don't use FixedPoint for the kernel
accumfilter(pixelval::N0f8, filterval::N0f8) = Float32(pixelval)*Float32(filterval)
accumfilter(pixelval::Colorant{N0f8}, filterval::N0f8) = float32(c)*Float32(filterval)

# In theory, the following might need to be specialized. For safety, make it a
# standalone function call.
safe_for_prod(x, ref) = oftype(ref, x)