text:
push!(LOAD_PATH, "/Users/sean/GitRepos/ParticleMapper/ParticleGridMapper.jl")
push!(LOAD_PATH, "/Users/sean/GitRepos/Chia-Yu_Hu_ParticleMapper/OctreeBH.jl")
using ParticleGridMapper
using Test

include("test_amr.jl")


text:
const CACHE_MODE = Symbol(uppercase(@load_preference("cache_mode", "DEFAULT")))
const CACHE_DIR = @load_preference("cache_dir", joinpath(DEPOT_PATH[1], "datadeps", "JuliaConSchedule"))
const TIMEOUT = parse(Float64, @load_preference("timeout", "5.0"))
const TERMINAL_LINKS = parse(Bool, @load_preference("terminal_links", "false"))

function set_cachemode(mode::Symbol)
    @assert mode in (:DEFAULT, :NEVER, :ALWAYS)
    @set_preferences!("cache_mode" => uppercase(string(mode)))
    @info("New cache mode set; restart your Julia session for this change to take effect!")
end
function set_cachedir(path::String)
    @set_preferences!("cache_dir" => string(path))
    @info("New cache dir set; restart your Julia session for this change to take effect!")
end
function set_timeout(val::Union{Float64, Int})
    @assert val > 0
    @set_preferences!("timeout" => string(val))
    @info("New timeout set; restart your Julia session for this change to take effect!")
end
function set_terminallinks(on::Bool)
    @set_preferences!("terminal_links" => string(on))
    @info("New terminal links preference set; restart your Julia session for this change to take effect!")
end

const PRETALX_JSON_URL = "https://pretalx.com/juliacon2020/schedule/export/schedule.json"
const DATA_ARCHIVE_JSON_URL = "https://raw.githubusercontent.com/JuliaCon/JuliaConDataArchive/master/juliacon2020_schedule/schedule.json"
const jcon = Ref{JuliaConSchedule}()

default_json_url() = DATA_ARCHIVE_JSON_URL
default_now() = Dates.now()

"""
    debugmode(on::Bool=true)

Simulates that we are live / in the middle of JuliaCon.
"""
function debugmode(on::Bool=true)
    if on
        @eval JuliaCon default_now() = Dates.DateTime("2020-07-29T16:30:00.000") # JuliaCon2020
    else
        @eval JuliaCon default_now() = Dates.now()
    end
    return nothing
end

"""
    set_json_source(src::Symbol)

Anticipated input: `:pretalx`, `:github` (JuliaConDataArchive)
"""
function set_json_source(src::Symbol)
    if src == :pretalx
        @eval JuliaCon default_json_url() = PRETALX_JSON_URL
    else
        @eval JuliaCon default_json_url() = DATA_ARCHIVE_JSON_URL
    end
    return nothing
end

ext: jl


text:
const γ_exact = 1.4
const μ_exact = 1 // 100
@noinline ρ_g(t, x, y, z, ::Val{2}) =
    sin(pi * x) * cos(pi * t) * cos(pi * y) + 3
@noinline U_g(t, x, y, z, ::Val{2}) =
    (sin(pi * x) * cos(pi * t) * cos(pi * y) + 3) *
    sin(pi * x) *
    cos(pi * t) *
    cos(pi * y)
@noinline V_g(t, x, y, z, ::Val{2}) =
    (sin(pi * x) * cos(pi * t) * cos(pi * y) + 3) *
    sin(pi * x) *
    cos(pi * t) *
    cos(pi * y)
@noinline W_g(t, x, y, z, ::Val{2}) = 0
@noinline E_g(t, x, y, z, ::Val{2}) =
    sin(pi * x) * cos(pi * t) * cos(pi * y) + 100
@noinline Sρ_g(t, x, y, z, ::Val{2}) =
    pi * (
        2 * sin(2 * pi * x) - 2 * sin(2 * pi * y) - sin(pi * (2 * t - 2 * x)) +
        sin(pi * (2 * t + 2 * x)) +
        sin(pi * (2 * t - 2 * y)) - sin(pi * (2 * t + 2 * y)) +
        2 * sin(pi * (2 * x + 2 * y)) +
        sin(pi * (-2 * t + 2 * x + 2 * y)) +
        sin(pi * (2 * t + 2 * x + 2 * y)) +
        10 * cos(pi * (-t + x + y)) - 2 * cos(pi * (t - x + y)) +
        2 * cos(pi * (t + x - y)) +
        14 * cos(pi * (t + x + y))
    ) / 8
@noinline SU_g(t, x, y, z, ::Val{2}) =
    pi * (
        -2.0 * sin(pi * t) * sin(pi * x)^2 * cos(pi * t) * cos(pi * y)^2 -
        3.0 * sin(pi * t) * sin(pi * x) * cos(pi * y) -
        3.0 * sin(pi * x)^3 * sin(pi * y) * cos(pi * t)^3 * cos(pi * y)^2 -
        6.0 * sin(pi * x)^2 * sin(pi * y) * cos(pi * t)^2 * cos(pi * y) +
        1.8 * sin(pi * x)^2 * cos(pi * t)^3 * cos(pi * x) * cos(pi * y)^3 +
        3.6 * sin(pi * x) * cos(pi * t)^2 * cos(pi * x) * cos(pi * y)^2 +
        0.0233333333333333 * pi * sin(pi * x) * cos(pi * t) * cos(pi * y) +
        0.00333333333333333 * pi * sin(pi * y) * cos(pi * t) * cos(pi * x) +
        0.4 * cos(pi * t) * cos(pi * x) * cos(pi * y)
    )
@noinline SV_g(t, x, y, z, ::Val{2}) =
    pi * (
        -2.0 * sin(pi * t) * sin(pi * x)^2 * cos(pi * t) * cos(pi * y)^2 -
        3.0 * sin(pi * t) * sin(pi * x) * cos(pi * y) -
        1.8 * sin(pi * x)^3 * sin(pi * y) * cos(pi * t)^3 * cos(pi * y)^2 -
        3.6 * sin(pi * x)^2 * sin(pi * y) * cos(pi * t)^2 * cos(pi * y) +
        3.0 * sin(pi * x)^2 * cos(pi * t)^3 * cos(pi * x) * cos(pi * y)^3 -
        0.4 * sin(pi * x) * sin(pi * y) * cos(pi * t) +
        6.0 * sin(pi * x) * cos(pi * t)^2 * cos(pi * x) * cos(pi * y)^2 +
        0.0233333333333333 * pi * sin(pi * x) * cos(pi * t) * cos(pi * y) +
        0.00333333333333333 * pi * sin(pi * y) * cos(pi * t) * cos(pi * x)
    )
@noinline SW_g(t, x, y, z, ::Val{2}) = 0
@noinline SE_g(t, x, y, z, ::Val{2}) =
    pi * (
        -1.0 * sin(pi * t) * sin(pi * x) * cos(pi * y) +
        1.6 * sin(pi * x)^4 * sin(pi * y) * cos(pi * t)^4 * cos(pi * y)^3 +
        3.6 * sin(pi * x)^3 * sin(pi * y) * cos(pi * t)^3 * cos(pi * y)^2 -
        1.6 * sin(pi * x)^3 * cos(pi * t)^4 * cos(pi * x) * cos(pi * y)^4 -
        0.0233333333333333 *
        pi *
        sin(pi * x)^2 *
        sin(pi * y)^2 *
        cos(pi * t)^2 -
        2.8 * sin(pi * x)^2 * sin(pi * y) * cos(pi * t)^2 * cos(pi * y) -
        3.6 * sin(pi * x)^2 * cos(pi * t)^3 * cos(pi * x) * cos(pi * y)^3 +
        0.0466666666666667 *
        pi *
        sin(pi * x)^2 *
        cos(pi * t)^2 *
        cos(pi * y)^2 +
        0.0133333333333333 *
        pi *
        sin(pi * x) *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y) - 140.0 * sin(pi * x) * sin(pi * y) * cos(pi * t) +
        2.8 * sin(pi * x) * cos(pi * t)^2 * cos(pi * x) * cos(pi * y)^2 -
        0.0233333333333333 *
        pi *
        cos(pi * t)^2 *
        cos(pi * x)^2 *
        cos(pi * y)^2 + 140.0 * cos(pi * t) * cos(pi * x) * cos(pi * y)
    )
@noinline ρ_g(t, x, y, z, ::Val{3}) =
    sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) + 3
@noinline U_g(t, x, y, z, ::Val{3}) =
    (sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) + 3) *
    sin(pi * x) *
    cos(pi * t) *
    cos(pi * y) *
    cos(pi * z)
@noinline V_g(t, x, y, z, ::Val{3}) =
    (sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) + 3) *
    sin(pi * x) *
    cos(pi * t) *
    cos(pi * y) *
    cos(pi * z)
@noinline W_g(t, x, y, z, ::Val{3}) =
    (sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) + 3) *
    sin(pi * x) *
    sin(pi * z) *
    cos(pi * t) *
    cos(pi * y)
@noinline E_g(t, x, y, z, ::Val{3}) =
    sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) + 100
@noinline Sρ_g(t, x, y, z, ::Val{3}) =
    pi * (
        -sin(pi * t) * sin(pi * x) * cos(pi * y) * cos(pi * z) -
        2 *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z)^2 -
        sin(pi * x)^2 * sin(pi * z)^2 * cos(pi * t)^2 * cos(pi * y)^2 +
        sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 * cos(pi * z)^2 -
        3 * sin(pi * x) * sin(pi * y) * cos(pi * t) * cos(pi * z) +
        2 *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        3 * sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) +
        3 * cos(pi * t) * cos(pi * x) * cos(pi * y) * cos(pi * z)
    )
@noinline SU_g(t, x, y, z, ::Val{3}) =
    pi * (
        -2.0 *
        sin(pi * t) *
        sin(pi * x)^2 *
        cos(pi * t) *
        cos(pi * y)^2 *
        cos(pi * z)^2 -
        3.0 * sin(pi * t) * sin(pi * x) * cos(pi * y) * cos(pi * z) -
        3.0 *
        sin(pi * x)^3 *
        sin(pi * y) *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z)^3 -
        2.0 *
        sin(pi * x)^3 *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * y)^3 *
        cos(pi * z) +
        1.0 * sin(pi * x)^3 * cos(pi * t)^3 * cos(pi * y)^3 * cos(pi * z)^3 -
        6.0 *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z)^2 -
        0.6 *
        sin(pi * x)^2 *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z) -
        3.0 * sin(pi * x)^2 * sin(pi * z)^2 * cos(pi * t)^2 * cos(pi * y)^2 +
        1.8 *
        sin(pi * x)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z)^3 +
        3.0 * sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 * cos(pi * z)^2 -
        1.2 *
        sin(pi * x) *
        sin(pi * z)^2 *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 +
        3.6 *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        0.0333333333333333 *
        pi *
        sin(pi * x) *
        cos(pi * t) *
        cos(pi * y) *
        cos(pi * z) +
        0.00333333333333333 *
        pi *
        sin(pi * y) *
        cos(pi * t) *
        cos(pi * x) *
        cos(pi * z) -
        0.00333333333333333 *
        pi *
        cos(pi * t) *
        cos(pi * x) *
        cos(pi * y) *
        cos(pi * z) +
        0.4 * cos(pi * t) * cos(pi * x) * cos(pi * y) * cos(pi * z)
    )
@noinline SV_g(t, x, y, z, ::Val{3}) =
    pi * (
        -2.0 *
        sin(pi * t) *
        sin(pi * x)^2 *
        cos(pi * t) *
        cos(pi * y)^2 *
        cos(pi * z)^2 -
        3.0 * sin(pi * t) * sin(pi * x) * cos(pi * y) * cos(pi * z) +
        0.6 *
        sin(pi * x)^3 *
        sin(pi * y) *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z) -
        1.8 *
        sin(pi * x)^3 *
        sin(pi * y) *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z)^3 -
        2.0 *
        sin(pi * x)^3 *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * y)^3 *
        cos(pi * z) +
        1.0 * sin(pi * x)^3 * cos(pi * t)^3 * cos(pi * y)^3 * cos(pi * z)^3 +
        1.2 *
        sin(pi * x)^2 *
        sin(pi * y) *
        sin(pi * z)^2 *
        cos(pi * t)^2 *
        cos(pi * y) -
        3.6 *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z)^2 -
        3.0 * sin(pi * x)^2 * sin(pi * z)^2 * cos(pi * t)^2 * cos(pi * y)^2 +
        3.0 *
        sin(pi * x)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z)^3 +
        3.0 * sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 * cos(pi * z)^2 -
        0.4 * sin(pi * x) * sin(pi * y) * cos(pi * t) * cos(pi * z) +
        0.00333333333333333 *
        pi *
        sin(pi * x) *
        sin(pi * y) *
        cos(pi * t) *
        cos(pi * z) +
        6.0 *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        0.0333333333333333 *
        pi *
        sin(pi * x) *
        cos(pi * t) *
        cos(pi * y) *
        cos(pi * z) +
        0.00333333333333333 *
        pi *
        sin(pi * y) *
        cos(pi * t) *
        cos(pi * x) *
        cos(pi * z)
    )
@noinline SW_g(t, x, y, z, ::Val{3}) =
    pi *
    (
        -2.0 *
        sin(pi * t) *
        sin(pi * x)^2 *
        cos(pi * t) *
        cos(pi * y)^2 *
        cos(pi * z) - 3.0 * sin(pi * t) * sin(pi * x) * cos(pi * y) -
        3.0 *
        sin(pi * x)^3 *
        sin(pi * y) *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z)^2 -
        0.8 * sin(pi * x)^3 * sin(pi * z)^2 * cos(pi * t)^3 * cos(pi * y)^3 +
        2.8 * sin(pi * x)^3 * cos(pi * t)^3 * cos(pi * y)^3 * cos(pi * z)^2 -
        6.0 *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z) +
        3.0 *
        sin(pi * x)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z)^2 +
        7.2 * sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 * cos(pi * z) -
        0.00333333333333333 * pi * sin(pi * x) * sin(pi * y) * cos(pi * t) +
        6.0 *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z) - 0.4 * sin(pi * x) * cos(pi * t) * cos(pi * y) +
        0.0333333333333333 * pi * sin(pi * x) * cos(pi * t) * cos(pi * y) +
        0.00333333333333333 * pi * cos(pi * t) * cos(pi * x) * cos(pi * y)
    ) *
    sin(pi * z)
@noinline SE_g(t, x, y, z, ::Val{3}) =
    pi * (
        0.001171875 *
        (1 - cos(2 * pi * x))^2 *
        (1 - cos(4 * pi * z)) *
        (cos(2 * pi * t) + 1)^2 *
        (cos(2 * pi * y) + 1)^2 -
        1.0 * sin(pi * t) * sin(pi * x) * cos(pi * y) * cos(pi * z) +
        0.8 *
        sin(pi * x)^4 *
        sin(pi * y) *
        sin(pi * z)^2 *
        cos(pi * t)^4 *
        cos(pi * y)^3 *
        cos(pi * z)^2 +
        1.6 *
        sin(pi * x)^4 *
        sin(pi * y) *
        cos(pi * t)^4 *
        cos(pi * y)^3 *
        cos(pi * z)^4 +
        0.2 * sin(pi * x)^4 * sin(pi * z)^4 * cos(pi * t)^4 * cos(pi * y)^4 -
        0.4 * sin(pi * x)^4 * cos(pi * t)^4 * cos(pi * y)^4 * cos(pi * z)^4 +
        1.8 *
        sin(pi * x)^3 *
        sin(pi * y) *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z) +
        3.6 *
        sin(pi * x)^3 *
        sin(pi * y) *
        cos(pi * t)^3 *
        cos(pi * y)^2 *
        cos(pi * z)^3 -
        0.8 *
        sin(pi * x)^3 *
        sin(pi * z)^2 *
        cos(pi * t)^4 *
        cos(pi * x) *
        cos(pi * y)^4 *
        cos(pi * z)^2 +
        0.6 *
        sin(pi * x)^3 *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * y)^3 *
        cos(pi * z) -
        1.6 *
        sin(pi * x)^3 *
        cos(pi * t)^4 *
        cos(pi * x) *
        cos(pi * y)^4 *
        cos(pi * z)^4 -
        1.2 * sin(pi * x)^3 * cos(pi * t)^3 * cos(pi * y)^3 * cos(pi * z)^3 -
        0.01 *
        pi *
        sin(pi * x)^2 *
        sin(pi * y)^2 *
        sin(pi * z)^2 *
        cos(pi * t)^2 -
        0.0233333333333333 *
        pi *
        sin(pi * x)^2 *
        sin(pi * y)^2 *
        cos(pi * t)^2 *
        cos(pi * z)^2 -
        0.0233333333333333 *
        pi *
        sin(pi * x)^2 *
        sin(pi * y) *
        sin(pi * z)^2 *
        cos(pi * t)^2 *
        cos(pi * y) -
        2.8 *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z)^2 -
        0.01 *
        pi *
        sin(pi * x)^2 *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * y) *
        cos(pi * z)^2 -
        1.8 *
        sin(pi * x)^2 *
        sin(pi * z)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z) -
        1.4 * sin(pi * x)^2 * sin(pi * z)^2 * cos(pi * t)^2 * cos(pi * y)^2 +
        0.0133333333333333 *
        pi *
        sin(pi * x)^2 *
        sin(pi * z)^2 *
        cos(pi * t)^2 *
        cos(pi * y)^2 -
        3.6 *
        sin(pi * x)^2 *
        cos(pi * t)^3 *
        cos(pi * x) *
        cos(pi * y)^3 *
        cos(pi * z)^3 +
        0.0533333333333333 *
        pi *
        sin(pi * x)^2 *
        cos(pi * t)^2 *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        1.4 * sin(pi * x)^2 * cos(pi * t)^2 * cos(pi * y)^2 * cos(pi * z)^2 +
        0.0133333333333333 *
        pi *
        sin(pi * x) *
        sin(pi * y) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y) *
        cos(pi * z)^2 -
        140.0 * sin(pi * x) * sin(pi * y) * cos(pi * t) * cos(pi * z) +
        0.0233333333333333 *
        pi *
        sin(pi * x) *
        sin(pi * z)^2 *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 +
        0.01 *
        pi *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        2.8 *
        sin(pi * x) *
        cos(pi * t)^2 *
        cos(pi * x) *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        140.0 * sin(pi * x) * cos(pi * t) * cos(pi * y) * cos(pi * z) -
        0.01 *
        pi *
        sin(pi * z)^2 *
        cos(pi * t)^2 *
        cos(pi * x)^2 *
        cos(pi * y)^2 -
        0.0233333333333333 *
        pi *
        cos(pi * t)^2 *
        cos(pi * x)^2 *
        cos(pi * y)^2 *
        cos(pi * z)^2 +
        140.0 * cos(pi * t) * cos(pi * x) * cos(pi * y) * cos(pi * z)
    )


text:
"""
	maximal_to_minimal_jacobian(mechanism, z)

	Jacobian of mapping from maximal to minimal representation

	mechanism: Mechanism
	z: maximal state
"""
function maximal_to_minimal_jacobian(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}, z::AbstractVector{Tz}) where {T,Nn,Ne,Nb,Ni,Tz}
	J = zeros(minimal_dimension(mechanism), maximal_dimension(mechanism) - Nb)
	timestep= mechanism.timestep
	row_shift = 0
	for id in mechanism.root_to_leaves
		(id > Ne) && continue # only treat joints
		joint = mechanism.joints[id]
		c_shift = 0
		v_shift = input_dimension(joint)
		ichild = joint.child_id - Ne
		for element in (joint.translational, joint.rotational)
			nu_element = input_dimension(element)

			c_idx = row_shift + c_shift .+ (1:nu_element)
			v_idx = row_shift + v_shift .+ (1:nu_element)

			xb, vb, qb, ωb = unpack_maximal_state(z, ichild)

			xb_idx = collect((ichild-1)*12 .+ (1:3))
			vb_idx = collect((ichild-1)*12 .+ (4:6))
			qb_idx = collect((ichild-1)*12 .+ (7:9))
			ωb_idx = collect((ichild-1)*12 .+ (10:12))

			if joint.parent_id != 0
				iparent = joint.parent_id - Ne
				xa, va, qa, ωa = unpack_maximal_state(z, iparent)

				xa_idx = collect((iparent-1)*12 .+ (1:3))
				va_idx = collect((iparent-1)*12 .+ (4:6))
				qa_idx = collect((iparent-1)*12 .+ (7:9))
				ωa_idx = collect((iparent-1)*12 .+ (10:12))

				J[c_idx, [xa_idx; qa_idx]] = minimal_coordinates_jacobian_configuration(:parent, element, xa, qa, xb, qb)
				J[v_idx, [xa_idx; qa_idx]] = minimal_velocities_jacobian_configuration(:parent, element, xa, va, qa, ωa, xb, vb, qb, ωb, timestep)
				J[v_idx, [va_idx; ωa_idx]] = minimal_velocities_jacobian_velocity(:parent, element, xa, va, qa, ωa, xb, vb, qb, ωb, timestep)
			else
				xa, va, qa, ωa = current_configuration_velocity(mechanism.origin.state)
			end

			J[c_idx, [xb_idx; qb_idx]] = minimal_coordinates_jacobian_configuration(:child, element, xa, qa, xb, qb)
			J[v_idx, [xb_idx; qb_idx]] = minimal_velocities_jacobian_configuration(:child, element, xa, va, qa, ωa, xb, vb, qb, ωb, timestep)
			J[v_idx, [vb_idx; ωb_idx]] = minimal_velocities_jacobian_velocity(:child, element, xa, va, qa, ωa, xb, vb, qb, ωb, timestep)

			c_shift += nu_element
			v_shift += nu_element
		end
		row_shift += 2 * input_dimension(joint)
	end
	return J
end

"""
    get_maximal_gradients!(mechanism, z, u; opts)

    return maximal gradients for mechanism
    note: this requires simulating the mechanism for one time step

    mechanism: Mechanism
    z: state
    u: input
    opts: SolverOptions
"""
function get_maximal_gradients!(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}, z::AbstractVector{T}, u::AbstractVector{T};
    opts=SolverOptions()) where {T,Nn,Ne,Nb,Ni}

    step!(mechanism, z, u, opts=opts)
    jacobian_state, jacobian_control = get_maximal_gradients(mechanism)

    return jacobian_state, jacobian_control
end

function get_maximal_gradients(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}) where {T,Nn,Ne,Nb,Ni}
	timestep= mechanism.timestep
	nu = input_dimension(mechanism)

	for entry in mechanism.data_matrix.nzval # reset matrix
		entry.value .= 0.0
	end
	jacobian_data!(mechanism.data_matrix, mechanism)
	nodes = [mechanism.joints; mechanism.bodies; mechanism.contacts]
	dimrow = length.(nodes)
	dimcol = data_dim.(nodes)
	index_row = [1+sum(dimrow[1:i-1]):sum(dimrow[1:i]) for i in 1:length(dimrow)]
	index_col = [1+sum(dimcol[1:i-1]):sum(dimcol[1:i]) for i in 1:length(dimcol)]

	index_state = [index_col[body.id][[14:16; 8:10; 17:19; 11:13]] for body in mechanism.bodies] # ∂ x2 v15 q2 ω15
	index_control = [index_col[joint.id][1:input_dimension(joint)] for joint in mechanism.joints] # ∂ u

	datamat = full_matrix(mechanism.data_matrix, dimrow, dimcol)
	solmat = full_matrix(mechanism.system)

	# data Jacobian
	data_jacobian = solmat \ datamat #TODO: use pre-factorization

	# Jacobian
	jacobian_state = zeros(12Nb,12Nb)
	jacobian_control = zeros(12Nb,nu)
	for (i, body) in enumerate(mechanism.bodies)
		id = body.id
		# Fill in gradients of v25, ω25
		jacobian_state[12*(i-1) .+ [4:6; 10:12],:] += data_jacobian[index_row[id], vcat(index_state...)]
		jacobian_control[12*(i-1) .+ [4:6; 10:12],:] += data_jacobian[index_row[id], vcat(index_control...)]

		# Fill in gradients of x3, q3
		x2 = body.state.x2
		q2 = body.state.q2
		v25 = body.state.vsol[2]
		ω25 = body.state.ωsol[2]
		q3 = next_orientation(q2, ω25, timestep)
		jacobian_state[12*(i-1) .+ (1:3), :] += linear_integrator_jacobian_velocity(x2, v25, timestep) * data_jacobian[index_row[id][1:3], vcat(index_state...)]
		jacobian_state[12*(i-1) .+ (1:3), 12*(i-1) .+ (1:3)] += linear_integrator_jacobian_position(x2, v25, timestep)
		jacobian_state[12*(i-1) .+ (7:9), :] += LVᵀmat(q3)' * rotational_integrator_jacobian_velocity(q2, ω25, timestep) * data_jacobian[index_row[id][4:6], vcat(index_state...)]
		jacobian_state[12*(i-1) .+ (7:9), 12*(i-1) .+ (7:9)] += LVᵀmat(q3)' * rotational_integrator_jacobian_orientation(q2, ω25, timestep, attjac=true)

		jacobian_control[12*(i-1) .+ (1:3),:] += linear_integrator_jacobian_velocity(x2, v25, timestep) * data_jacobian[index_row[id][1:3], vcat(index_control...)]
		jacobian_control[12*(i-1) .+ (7:9),:] += LVᵀmat(q3)' * rotational_integrator_jacobian_velocity(q2, ω25, timestep) * data_jacobian[index_row[id][4:6], vcat(index_control...)]
	end

	return jacobian_state, jacobian_control
end

"""
	minimal_to_maximal_jacobian(mechanism, x)

	Jacobian of mapping from minimal to maximal representation

	mechanism: Mechanism
	y: minimal state
"""
function minimal_to_maximal_jacobian(mechanism::Mechanism{T,Nn,Ne,Nb,Ni}, x::AbstractVector{Tx}) where {T,Nn,Ne,Nb,Ni,Tx}
	timestep= mechanism.timestep
	J = zeros(maximal_dimension(mechanism, attjac=true), minimal_dimension(mechanism))

	# Compute partials
	partials = Dict{Vector{Int}, Matrix{T}}()
	for cnode in mechanism.bodies
		for joint in parent_joints(mechanism, cnode)
			pnode = get_node(mechanism, joint.parent_id, origin=true)
			partials[[cnode.id, joint.id]] = set_minimal_coordinates_velocities_jacobian_minimal(joint, pnode, cnode, timestep) # 12 x 2nu (xvqω x Δxθvω)
			partials[[cnode.id, pnode.id]] = set_minimal_coordinates_velocities_jacobian_parent(joint, pnode, cnode, timestep) # 12 x 12 (xvqω x xvqω)
		end
	end

	# Index
	row = [12(i-1)+1:12i for i = 1:Nb]
	col = [] # ordering joints from root to tree
	col_idx = zeros(Int,Ne)
	cnt = 0
	for id in mechanism.root_to_leaves
		(id > Ne) && continue # only keep joints
		cnt += 1
		nu = input_dimension(get_joint(mechanism, id))
		if length(col) > 0
			push!(col, col[end][end] .+ (1:2nu))
		else
			push!(col, 1:2nu)
		end
		col_idx[id] = cnt
	end

	 # chain partials together from root to leaves
	for id in mechanism.root_to_leaves
		!(Ne < id <= Ne+Nb) && continue # only treat bodies
		cnode = get_node(mechanism, id)
		for joint in parent_joints(mechanism, cnode)
			pnode = get_node(mechanism, joint.parent_id, origin=true)
			J[row[cnode.id-Ne], col[col_idx[joint.id]]] += partials[[cnode.id, joint.id]] # ∂zi∂θp(i)
			(pnode.id == 0) && continue # avoid origin
			J[row[cnode.id-Ne], :] += partials[[cnode.id, pnode.id]] * J[row[pnode.id-Ne], :] # ∂zi∂zp(p(i)) * ∂zp(p(i))/∂θ
		end
	end
	return J
end

"""
    get_minimal_gradients!(mechanism, y, u; opts)

    return minimal gradients for mechanism
    note: this requires simulating the mechanism for one time step

    mechanism: Mechanism
    y: state
    u: input
	opts: SolverOptions
"""
function get_minimal_gradients!(mechanism::Mechanism{T}, y::AbstractVector{T}, u::AbstractVector{T};
	opts=SolverOptions()) where T
	# simulate next state
	step!(mechanism, y, u, opts=opts)
	return get_minimal_gradients!(mechanism)
end

function get_minimal_gradients!(mechanism::Mechanism{T}) where T
	# current maximal state
	z = get_maximal_state(mechanism)
	# next maximal state
	z_next = get_next_state(mechanism)
	# current minimal state
	x = maximal_to_minimal(mechanism, z)
	# maximal dynamics Jacobians
	maximal_jacobian_state, minimal_jacobian_control = get_maximal_gradients(mechanism)
	# minimal to maximal Jacobian at current time step (rhs)
	min_to_max_jacobian_current = minimal_to_maximal_jacobian(mechanism, x)
	# maximal to minimal Jacobian at next time step (lhs)
	max_to_min_jacobian_next = maximal_to_minimal_jacobian(mechanism, z_next)
	# minimal state Jacobian
	minimal_jacobian_state = max_to_min_jacobian_next * maximal_jacobian_state * min_to_max_jacobian_current
	# minimal control Jacobian
	minimal_jacobian_control = max_to_min_jacobian_next * minimal_jacobian_control
	return minimal_jacobian_state, minimal_jacobian_control
end


text:
#= The Computer Language Benchmarks Game
   https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

   contributed by Ben Lauwens

   derived from the Chapel version by Tom Hildebrandt, Brad Chamberlain, and Lydia Duncan
   and inspired by the Julia version of Luiz M. Faria
=#

using ResumableFunctions

using Base.GMP.MPZ: add_ui!, mul_ui!, add!, tdiv_q!
const mpz_t = Ref{BigInt}
addmul_ui!(x::BigInt,a::BigInt,b::UInt64) = ccall((:__gmpz_addmul_ui,"libgmp"), Cvoid, (mpz_t,mpz_t,Culong), x, a, b)
submul_ui!(x::BigInt,a::BigInt,b::UInt64) = ccall((:__gmpz_submul_ui,"libgmp"), Cvoid, (mpz_t,mpz_t,Culong), x, a, b)
get_ui(x::BigInt) = ccall((:__gmpz_get_ui,"libgmp"), Culong, (mpz_t,), x)

@resumable function pidigits()
    k = one(UInt64)
    d = zero(UInt64)
    numer = BigInt(1)
    denom = BigInt(1)
    accum = BigInt(0)
    tmp1 = BigInt(0)
    tmp2 = BigInt(0)
    tmp3 = BigInt(0)
    while true
        k2 = 2k + one(UInt64)
        addmul_ui!(accum, numer, UInt64(2))
        mul_ui!(accum, k2)
        mul_ui!(denom, k2)
        mul_ui!(numer, k)
        k += one(UInt64)
        if numer <= accum
            mul_ui!(tmp1, numer, UInt64(3))
            add!(tmp1, accum)
            tdiv_q!(tmp2, tmp1, denom)
            mul_ui!(tmp1, numer, UInt64(4))
            add!(tmp1, accum)
            tdiv_q!(tmp3, tmp1, denom)
            if tmp2 == tmp3
                d = get_ui(tmp2)
                @yield d
                submul_ui!(accum, denom, d)
                mul_ui!(accum, UInt64(10))
                mul_ui!(numer, UInt64(10))
            end
        end
    end
end

function main(n = parse(Int64, ARGS[1]))
    i = 0
    pid = pidigits()
    for i in 1:n
        d = pid()
        print(d)
        if i % 10 === 0
            println("\t:", i)
        end
    end
    if i % 10 !== 0
        println(" "^(10 - (i % 10)), "\t:", n)
    end
end

main(100)


text:
module BoteSalvatICX

using Requires

include("xione.jl")

export ionizationcrosssection  # Computes the ionization crosssection for electron impact on free atoms
export hasedge # Is data available for the specified element and shell
export edgeenergy  # Bote-Salvat's assumed edge energy in eV

function __init__()
    @require Gadfly="c91e804a-d5a3-530f-b6f0-dfbca275c004" include("gadflysupport.jl")
end

end # module


text:
module Utils

using ThreadedScans

struct WithKwargs
    f::Any
    overrides::NamedTuple
end

(f::WithKwargs)(args...; kwargs...) = f.f(args...; kwargs..., f.overrides...)
with_kwargs(f; overrides...) = WithKwargs(f, (; overrides...))
with_kwargs(f::WithKwargs; overrides...) = WithKwargs(f.f, (; f.overrides..., overrides...))

unwrap_function(f) = f
unwrap_function(f::WithKwargs) = unwrap_function(f.f)

function Base.show(io::IO, f::WithKwargs)
    print(io, parentmodule(f.f), '.')
    print(io, nameof(f.f), "(...;")
    isfirst = true
    for (k, v) in pairs(f.overrides)
        if isfirst
            isfirst = false
        else
            print(io, ',')
        end
        print(io, ' ', k, " = ", v)
    end
    print(io, ')')
end

const SCAN_FUNCTIONS = Any[
    ThreadedScans.simple!,
    ThreadedScans.dac!,
    ThreadedScans.linear!,
    ThreadedScans.partitioned_hillis_steele!,
    # ...
]

function scan_functions_with_kwargs()
    fns = copy(SCAN_FUNCTIONS)
    push!(
        fns,
        with_kwargs(ThreadedScans.simple!; ntasks = 7),
        with_kwargs(ThreadedScans.dac!; ntasks = 7),
        with_kwargs(ThreadedScans.dac!; ntasks = 8),
        with_kwargs(ThreadedScans.linear!; ntasks = 7),
        with_kwargs(ThreadedScans.linear!; ntasks = 8),
        with_kwargs(ThreadedScans.partitioned_hillis_steele!; ntasks = 7),
        with_kwargs(ThreadedScans.partitioned_hillis_steele!; ntasks = 8),
    )
    for f in copy(fns)
        unwrap_function(f) === ThreadedScans.simple! && continue
        push!(fns, with_kwargs(f, spin = 100))
    end
    return fns
end

end  # module


text:
"""
Obtain coefficients (A, B) of the function Se = B(x - A)^2 for
Se(1.2) = B(1.2 - A)^2 and Se(1.0) = B(1.0 - A)^2 as:
Se(1.0) = Se(1.2)/(1.2 - A)^2 * (1.0 - A)^2 that yields
(1.2 - A)^2 Se(1.0) = Se(1.2) * (1.0 - A)^2 or expanding:
(Se(1.2) - Se(1.0)) A^2 + (2.4 Se(1.0) - 2 Se(1.2)) A + (Se(1.2) - 1.44 Se(1.0)) = 0
and uses the negative solution of the quadratic equation
"""
function get_quadratic_saturation(Se::Tuple{Float64, Float64})
    if Se[1] == 0.0 || Se[2] == 0.0
        @warn "$(Se[1]) or $(Se[2]) equals to zero. Ignoring saturation."
        return (0.0, 0.0)
    end
    if Se[2] <= Se[1]
        throw(IS.ConflictingInputsError("Se(1.2) <= Se(1.0). Saturation data is inconsistent."))
    end
    E1 = 1.0
    E12 = 1.2
    Sat_Se1 = Se[1]
    Sat_Se12 = Se[2]
    Sat_a = sqrt(E1 * Sat_Se1 / (E12 * Sat_Se12)) * ((Sat_Se12 > 0) + (Sat_Se12 < 0))
    Sat_A = E12 - (E1 - E12) / (Sat_a - 1)
    Sat_B = E12 * Sat_Se12 * (Sat_a - 1)^2 * ((Sat_a > 0) + (Sat_a < 0)) / (E1 - E12)^2

    return (Sat_A, Sat_B)
end

"""
Obtain coefficients (A, B) of the function Se = Bx^A for
Se(1.2) = B(1.2)^A and Se(1.0) = B(1.0)^A as:
B = Se(1.0) and hence
(1.2)^A = Se(1.2)/B -> A = log(Se(1.2)/B) / log(1.2)
"""
function get_exponential_saturation(Se::Tuple{Float64, Float64})
    if Se[1] == 0.0 || Se[2] == 0.0
        @warn "$(Se[1]) or $(Se[2]) equals to zero. Ignoring saturation."
        return (0.0, 0.0)
    end
    if Se[2] <= Se[1]
        throw(IS.ConflictingInputsError("Se(1.2) <= Se(1.0). Saturation data is inconsistent."))
    end
    B = Se[1]
    A = log(Se[2] / B) / log(1.2)

    return (A, B)
end


text:
using Test
using ThreeBodyTB
using Suppressor

#include("../includes_laguerre.jl")
#include("../Ewald.jl")

TESTDIR=ThreeBodyTB.TESTDIR

function test_basics()

    @testset "test basic scf Hk TB manipulation" begin

         @suppress begin
#        if true

            units_old = ThreeBodyTB.set_units()
            ThreeBodyTB.set_units(both="atomic")

            types=["Al"];

            #positions, crystal units
            pos=zeros((1,3));

            #lattice vectors, in Bohr units currently
            A=[ [3.7962751610 3.7962751610 0]; [3.7962751610 0 3.7962751610 ]; [ 0 3.7962751610 3.7962751610]];

            #makes the crystal
            c=makecrys(A, pos, types, units="Bohr")

            energy, tbc, flag = scf_energy(c)
            energyS, tbcS, flagS = scf_energy(c, mixing_mode=:simple)

#            @test abs(energy - -0.28986347423986025) < 1e-2
            @test abs(energy - -0.28707042) < 1e-2
            @test abs(energyS - energy) < 1e-5

            energy_fft = ThreeBodyTB.TB.calc_energy_fft(tbc)

            directgap, indirectgap, gaptype, bandwidth = ThreeBodyTB.BandStruct.band_summary(tbc)
            @test gaptype == :metal
            @test abs(indirectgap) < 1e-5
            @test abs(directgap) < 1e-5
            @test abs(bandwidth - 0.8334073432114558) < 1e-5

            @test abs(energy - energy_fft) < 1e-5

            energy_calc = ThreeBodyTB.TB.calc_energy(tbc)

            @test abs(energy - energy_calc) < 1e-5

            energy2, force, stress, tbc2 = scf_energy_force_stress(tbc)

            @test abs(energy - energy2) < 1e-4
            @test sum(abs.(force)) < 1e-4
            @test sum(abs.(stress))  < 1e-2

            vects, vals, ham, S,e =  Hk(tbc, [0 0 0 ])

#            @test abs( (vals[2] - vals[1]) - 1.9352469029181867) < 1e-2   #basic eigen value testing
            @test abs( (vals[2] - vals[1]) - 1.9585164167567397) < 1e-2   #basic eigen value testing
            @test abs( (vals[3] - vals[2]) ) < 1e-4
            @test abs( (vals[4] - vals[2]) ) < 1e-4

            vects2, vals2, ham2, S2,e2 =  Hk(tbc.tb, [0 0 0 ])

            @test sum(abs.(ham2 - ham)) < 1e-7

            tbck = ThreeBodyTB.TB.read_tb_crys_kspace("$TESTDIR/data_forces/al_fcc_projham_K.xml.gz")

            directgap, indirectgap, gaptype, bandwidth = ThreeBodyTB.BandStruct.band_summary(tbck)
            @test gaptype == :metal
            @test abs(indirectgap) < 1e-5
            @test abs(directgap) < 1e-5
            @test abs(bandwidth - 0.8418566664739338) < 1e-2

            vects3, vals3, ham3, S3,e3 =  Hk(tbck, [0 0 0 ])

            @test sum(abs.(vals3[1] - vals[1])) < 0.05

            tb = deepcopy(tbc.tb)
            ThreeBodyTB.TB.trim(tb)

            vectsT, valsT, hamT, ST,eT =  Hk(tb, [0 0 0 ])

            @test sum(abs.(valsT - vals2)) < 2e-2
#            println("VALST $valsT")
#            println("VALS2 $vals2")

            data_onsite, data_arr = ThreeBodyTB.TB.organizedata(tbc.crys, tbc.tb)

#            @test abs(data_onsite[1,5] - -0.677103) < 1e-2  #onsite s orbital
#            @test abs(data_onsite[1,7] - 5.36874) < 1e-2 #n.n. distance

            @test abs(data_onsite[1,5] - -0.696210333954235) < 1e-2  #onsite s orbital
            @test abs(data_onsite[1,7] - 5.368743819186305) < 1e-2 #n.n. distance
            @test abs(data_onsite[1,11] - 1.0) < 1e-3 # overlap onsite

            h1, dq = ThreeBodyTB.TB.get_h1(tbc)

            @test dq[1] < 1e-5

            en, eden, VECTS, VALS, err = ThreeBodyTB.TB.get_energy_electron_density_kspace(tbck)

            @test abs(en - -0.28962628722369677) < 1e-3

            tbck2 = ThreeBodyTB.SCF.remove_scf_from_tbc(tbck)
            tbc2 = ThreeBodyTB.SCF.remove_scf_from_tbc(tbc)

            en2, eden2, VECTS2, VALS2, err2 = ThreeBodyTB.TB.get_energy_electron_density_kspace(tbck2)

            @test abs(en - en2) < 1e-3

            H = ThreeBodyTB.Atomdata.atoms["H"]

            @test H.Z == 1

            energies, dos, pdos, names  =  ThreeBodyTB.DOS.gaussian_dos(tbc, do_display=false)

            @test length(energies) == length(dos)

            energies, dos, pdos, names  =  ThreeBodyTB.DOS.gaussian_dos(tbck, do_display=false)

            @test length(energies) == length(dos)
            @test abs(sum(dos) * (energies[2] - energies[1]) - 4.0) < 1e-2
            ThreeBodyTB.set_units(energy = units_old[1], length=units_old[2])

        end
    end
end

test_basics()


text:
# Autogenerated wrapper script for Entwine_jll for x86_64-w64-mingw32-cxx11
export entwine, libentwine

using PDAL_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("Entwine")
JLLWrappers.@declare_executable_product(entwine)
JLLWrappers.@declare_library_product(libentwine, "libentwine.dll")
function __init__()
    JLLWrappers.@generate_init_header(PDAL_jll, OpenSSL_jll)
    JLLWrappers.@init_executable_product(
        entwine,
        "bin\\entwine.exe",
    )

    JLLWrappers.@init_library_product(
        libentwine,
        "bin\\libentwine.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# FITS methods

const VERBOSE_MODE = @compat Dict("r"=>"read-only",
                                  "w"=>"read-write",
                                  "r+"=>"append")

# helper function for show()
function show_ascii_table(io, names, cols, spaces=2, indent=0)
    ncols = length(cols)
    ncols >= 1 || error("No columns")
    nrows = length(cols[1])
    length(names) == ncols || error("length of cols and names must match")
    for i=1:ncols
        length(cols[i]) == nrows || error("column length mismatch")
    end

    lengths = [max(maximum(length, cols[i]), length(names[i])) + spaces
               for i=1:ncols]
    for i = 1:ncols
        print(io, rpad(names[i], lengths[i]))
    end
    print(io, "\n")
    for j = 1:nrows
        print(io, " "^indent)
        for i=1:ncols
            print(io, rpad(cols[i][j], lengths[i]))
        end
        print(io, "\n")
    end
end

function length(f::FITS)
    fits_assert_open(f.fitsfile)
    @compat Int(fits_get_num_hdus(f.fitsfile))
end

endof(f::FITS) = length(f)

# Iteration
start(f::FITS) = 1
next(f::FITS, state) = (f[state], state + 1)
done(f::FITS, state) = state > length(f)

function show(io::IO, f::FITS)
    fits_assert_open(f.fitsfile)

    print(io, """File: $(f.filename)
    Mode: $(repr(f.mode)) ($(VERBOSE_MODE[f.mode]))
    """)

    nhdu = length(f)

    if nhdu == 0
        print(io, "No HDUs.")
    else
        print(io, "HDUs: ")

        names = Array(Compat.ASCIIString, nhdu)
        vers = Array(Compat.ASCIIString, nhdu)
        types = Array(Compat.ASCIIString, nhdu)
        for i = 1:nhdu
            t = fits_movabs_hdu(f.fitsfile, i)
            types[i] = (t == :image_hdu ? "Image" :
                        t == :binary_table ? "Table" :
                        t == :ascii_table ? "ASCIITable" :
                        error("unknown HDU type"))
            nname = fits_try_read_extname(f.fitsfile)
            names[i] = get(nname, "")
            nver = fits_try_read_extver(f.fitsfile)
            vers[i] = isnull(nver) ? "" : string(get(nver))
        end

        nums = [string(i) for i=1:nhdu]

        # only display version info if present
        if maximum(length, vers) > 0
            dispnames = ["Num", "Name", "Ver", "Type"]
            dispcols = Vector{Compat.ASCIIString}[nums, names, vers, types]
        else
            dispnames = ["Num", "Name", "Type"]
            dispcols = Vector{Compat.ASCIIString}[nums, names, types]
        end

        show_ascii_table(io, dispnames, dispcols, 2, 6)
    end
end

# Returns HDU object based on extension number
function getindex(f::FITS, i::Integer)
    fits_assert_open(f.fitsfile)

    if haskey(f.hdus, i)
        return f.hdus[i]
    end

    if i > length(f)
        error("index out of bounds")
    end
    hdutype = fits_movabs_hdu(f.fitsfile, i)
    f.hdus[i] = (hdutype == :image_hdu ? ImageHDU(f.fitsfile, i) :
                 hdutype == :binary_table ? TableHDU(f.fitsfile, i) :
                 hdutype == :ascii_table ? ASCIITableHDU(f.fitsfile, i) :
                 error("bad HDU type"))
    return f.hdus[i]
end

# Returns HDU based on hduname, version
function getindex(f::FITS, name::AbstractString, ver::Int=0)
    fits_assert_open(f.fitsfile)
    fits_movnam_hdu(f.fitsfile, name, ver)
    i = fits_get_hdu_num(f.fitsfile)

    if haskey(f.hdus, i)
        return f.hdus[i]
    end

    hdutype = fits_get_hdu_type(f.fitsfile)
    f.hdus[i] = (hdutype == :image_hdu ? ImageHDU(f.fitsfile, i) :
                 hdutype == :binary_table ? TableHDU(f.fitsfile, i) :
                 hdutype == :ascii_table ? ASCIITableHDU(f.fitsfile, i) :
                 error("bad HDU type"))
    return f.hdus[i]
end

function close(f::FITS)
    fits_assert_open(f.fitsfile)
    fits_close_file(f.fitsfile)
    f.filename = ""
    f.mode = ""
    empty!(f.hdus)
    nothing
end


text:
abstract type AbstractSolarData end
struct SolarData{T1<:AF} <: AbstractSolarData
    wav::Dict{Tuple{Symbol,Symbol}, AbstractArray{T1,2}}
    bis::Dict{Tuple{Symbol,Symbol}, AbstractArray{T1,2}}
    dep::Dict{Tuple{Symbol,Symbol}, AbstractArray{T1,2}}
    wid::Dict{Tuple{Symbol,Symbol}, AbstractArray{T1,2}}
    len::Dict{Tuple{Symbol,Symbol}, Int}
    line::AbstractLineProperties
end

function clean_input(wavall::AA{T,2}, bisall::AA{T,2}, depall::AA{T,2}, widall::AA{T,2}) where T<:AF
    @assert size(wavall) == size(bisall)
    @assert size(depall) == size(widall)

    # make boolean array (column will be stripped if badcol[i] == true)
    badcols = zeros(Bool, size(wavall,2))

    # find spread of data
    wav_std = std(wavall, dims=2)
    wid_std = std(wavall, dims=2)

    # loop through checking for bad columns
    for i in 1:size(wavall,2)
        # check for monotinicity
        if !ismonotonic(widall[:,i])
            badcols[i] = true
        end

        # check for bad width measurements
        if all(iszero(widall[:,i]))
            badcols[i] = true
        end

        # check for excessive NaNs
        idx = findfirst(x->isnan(x), wavall[:, i])
        if bisall[idx,i] < 0.85
            badcols[i] = true
        end

        # remove data that is significant outlier (bisector)
        if any(abs.(wavall[5:50,1] - wavall[5:50,i]) .> (4.0 .* wav_std[5:50]))
            badcols[i] = true
        end
    end

    # strip the bad columns and return new arrays
    return strip_columns(wavall, badcols), strip_columns(bisall, badcols),
           strip_columns(depall, badcols), strip_columns(widall, badcols)
end

function extrapolate_bisector(wavall::AA{T,2}, bisall::AA{T,2}; top::T=0.8) where T<:AF
    for i in 1:size(wavall,2)
        # take a slice for one time snapshot
        wavt = view(wavall, :, i)
        bist = view(bisall, :, i)

        # find the last index for good data
        ind1 = searchsortedfirst(bist, top)
        ind2 = searchsortedfirst(bist, top-0.1)

        # calculate an average slope and project it up to continuum
        dydx = mean((wavt[(ind2+1:ind1)] .- wavt[(ind2:ind1-1)]) ./ (bist[ind2+1:ind1] .- bist[(ind2:ind1-1)]))
        wavall[ind1:end, i] .= (dydx .* (bist[ind1:end] .- bist[ind1]) .+ wavt[ind1])
    end
    return nothing
end

function extrapolate_width(depall::AA{T,2}, widall::AA{T,2}) where T<:AF
    for i in 1:size(widall,2)
        # take slice for one time snapshot
        widt = view(widall, :, i)
        dept = view(depall, :, i)

        # fix the last value if necessary
        if (widt[end-1] == widt[end])
            dydx = (widt[end-1]/widt[end-2])/(dept[end-1]/dept[end-2])
            widall[end, i] = widt[end-2] / dydx + widt[end-2]
        end
    end
    return nothing
end

function relative_bisector_wavelengths(wav::AA{T,2}, lp::LineProperties) where T<:AF
    # allocate array of size(wav)
    relwav = similar(wav)
    for i in 1:size(wav,2)
        λgrav = (635.0/c_ms) * lp.λrest
        relwav[:,i] = wav[:,i] .- lp.λrest .- λgrav
    end
    return relwav
end

"""
    SolarData(; dir=soldir, relative=true, ...)

Construct a `SpecParams` composite type instance.

# Arguments
- `dir::String=soldir`: Directory containing the pre-processed input data. Default directory is set in src/config.jl
- `relative::Bool=true`: Set whether wavelengths are on absolute scale or expressed relative to rest wavelength.
"""
function SolarData(;dir::String=soldir, relative::Bool=true,
                   fixed_width::Bool=false, fixed_bisector::Bool=false,
                   extrapolate::Bool=true, contiguous_only::Bool=false,
                   adjust_mean::Bool=true)
    # sort the data directories and such
    bdf = sort_bisector_data(dir=dir)
    wdf = sort_width_data(dir=dir)

    # fill in line properties
    # TODO this is hard-coded for now
    iron_lp = LineProperties("FeI", 55.845, 0.86, 5434.5232, 0.0, 550.0)

    # allocate memory for data dictionaries
    wavdict = Dict{Tuple{Symbol,Symbol}, AA{Float64,2}}()
    bisdict = Dict{Tuple{Symbol,Symbol}, AA{Float64,2}}()
    depdict = Dict{Tuple{Symbol,Symbol}, AA{Float64,2}}()
    widdict = Dict{Tuple{Symbol,Symbol}, AA{Float64,2}}()
    lengths = Dict{Tuple{Symbol,Symbol}, Int}()

    # loop over unique mu + axis pairs
    for ax in unique(bdf.axis)
        for mu in unique(bdf.mu)
            # get data for mu + axis pair
            bdf_temp = bdf[((bdf.axis .== ax) .& (bdf.mu .== mu)), :]
            wdf_temp = wdf[((wdf.axis .== ax) .& (wdf.mu .== mu)), :]

            # move on if no data for mu + axis pair
            if isempty(bdf_temp)
                continue
            end

            # stitch the time series
            wavall, bisall = stitch_time_series(bdf_temp, adjust_mean=adjust_mean, contiguous_only=contiguous_only)
            depall, widall = stitch_time_series(wdf_temp, adjust_mean=adjust_mean, contiguous_only=contiguous_only)

            # clean the input
            wavall, bisall, depall, widall = clean_input(wavall, bisall, depall, widall)

            # deal with NaNs in wavall measurements
            if extrapolate
                extrapolate_bisector(wavall, bisall)
                # extrapolate_width(depall, widall)
            end

            # assign key-value pairs to dictionary
            if relative
                relwav = relative_bisector_wavelengths(wavall, iron_lp)
                wavdict[(Symbol(ax), Symbol(mu))] = relwav
            else
                wavdict[(Symbol(ax), Symbol(mu))] = wavall
            end
            bisdict[(Symbol(ax), Symbol(mu))] = bisall

            # assign width
            if fixed_width
                dep, wid = calc_fixed_width()
                depdict[(Symbol(ax), Symbol(mu))] = zeros(size(depall))
                widdict[(Symbol(ax), Symbol(mu))] = zeros(size(depall))
                depdict[(Symbol(ax), Symbol(mu))] .= dep
                widdict[(Symbol(ax), Symbol(mu))] .= wid
            else
                depdict[(Symbol(ax), Symbol(mu))] = depall
                widdict[(Symbol(ax), Symbol(mu))] = widall
            end

            # set bisector to zero if fixed
            if fixed_bisector
                wavdict[(Symbol(ax), Symbol(mu))] = zeros(size(wavall))
            end

            # get length
            lengths[(Symbol(ax), Symbol(mu))] = size(wavall,2)
        end
    end
    return SolarData(wavdict, bisdict, depdict, widdict, lengths, iron_lp)
end

function calc_fixed_width(;λ::T1=5434.5232, M::T1=26.0, T::T1=5778.0, v_turb::T1=3.5e5) where T1<:AF
    wid = width_thermal(λ=λ, M=M, T=T, v_turb=v_turb)
    λs = collect(range(λ - 1.0, λ + 1.0, step=λ/7e5))
    prof = gaussian_line.(λs, mid=λ, width=wid, depth=0.86)
    depall, widall = calc_width_at_depth(λs, prof, len=100)
    widall[1] = 0.01
    return depall, widall
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(13)
    n = 3000
    data = randn(n)
    fig = Figure(resolution = (1200, 800), font = "sans", fontsize = 20)
    ax1 = Axis(fig, xlabel = "value" )
    ax2 = Axis(fig, xlabel = "value" )
    ax3 = Axis(fig, xlabel = "value")
    ax4 = Axis(fig, xlabel = "value" )
    hist!(ax1, data,normalization = :none, color = :black, label = "none")
    hist!(ax2, data,normalization = :pdf, color = :orange, label = "pdf")
    hist!(ax3, data,normalization = :density,color = :dodgerblue,label = "density")
    hist!(ax4, data,normalization = :probability, label = "probability")
    axislegend(ax1, position = :rt)
    axislegend(ax2, position = :rt)
    axislegend(ax3, position = :rt)
    axislegend(ax4, position = :rt)
    fig[1,1] = ax1
    fig[1,2] = ax2
    fig[2,1] = ax3
    fig[2,2] = ax4
    fig
    save(joinpath(@__DIR__, "output", "histogramsNorms.png"), current_figure(), px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using IntervalLinearAlgebra, StaticArrays, LazySets
using Test

const IA = IntervalArithmetic

include("test_classify.jl")
include("test_multiplication.jl")
include("test_utils.jl")

include("test_eigenvalues/test_interval_eigenvalues.jl")
include("test_eigenvalues/test_verify_eigs.jl")
include("test_solvers/test_enclosures.jl")
include("test_solvers/test_epsilon_inflation.jl")
include("test_solvers/test_precondition.jl")
include("test_solvers/test_oettli_prager.jl")


text: @time reg(df, @model(y ~ x1 + x2, fe = id1))


text:
module NoiseCC

# mode-1: fft
include("FFT.jl")

# mode-1: cc
include("CC.jl")

# mode-2: fft and cc
include("FFT_CC.jl")

# stack
include("STACK.jl")

end

max_stars_repo_name: OUCyf/NoiseCC

max_issues_repo_name: OUCyf/NoiseCC

max_forks_repo_name: OUCyf/NoiseCC


text:
# Ported from: http://www.math.unl.edu/~tshores1/Public/Math496S06/MatlabTools/
#                                                Examples/chap7/examp1/bvls.m
# TODO: take a critical look at speed some time
#
# x=bvls(A,b,l,u)
#
# Solves a bounded variables least squares problem
#
#   min || Ax-b ||
#          l <= x <= u
#
# The algorithm is based on the bounded variables least squares algorithm
# of Stark and Parker.  See
#
#  P.B. Stark and R. L. Parker, "Bounded-Variable Least-Squares: An Algorithm
#  and Applications", Computational Statistics 10:129-141, 1995.
#

function bvls(A::AbstractMatrix, b::AbstractVector, l::AbstractVector, u::AbstractVector)
    #
    # Get the size of A for future reference.
    #
    m, n = size(A)
    #
    # Initialize a bunch of variables.  This speeds things up by avoiding
    # realloation of storage each time an entry is added to the end of a vector.
    #
    oopslist = falses(n)
    state = fill(0, n)
    x = zeros(n)
    atbound = falses(n)
    between = falses(n)
    criti = 0
    crits = 0
    myeps = 1.0e-10
    #
    # setup an initial solution with all vars locked at a lower or upper bound.
    # set any free vars to 0.
    #
    for i in 1:n
        if ((u[i] >= +Inf) && (l[i] <= -Inf))
            x[i] = 0
            state[i] = 0
            between[i] = true
        elseif (u[i] >= +Inf)
            x[i] = l[i]
            state[i] = 1
            atbound[i] = true
        elseif (l[i] <= -Inf)
            x[i] = u[i]
            state[i] = 2
            atbound[i] = true
        elseif (abs(l[i]) <= abs(u[i]))
            x[i] = l[i]
            state[i] = 1
            atbound[i] = true
        else
            x[i] = u[i]
            state[i] = 2
            atbound[i] = true
        end
    end
    #
    # The main loop.  Stop after 10*n iterations if nothing else.
    #
    iter = 0
    while (iter < 10*n)
        iter += 1
        #
        # optimality test.
        #
        grad = A'*(A*x-b)
        #
        # We ignore any variable that has already been tried and failed.
        #
        grad[oopslist] = 0.
        #
        # Check for optimality.
        #
        done = true
        for i in 1:n
            if ((abs(grad[i]) > (1+norm(b))*myeps) && (state[i] == 0)) ||
                   ((grad[i] < 0) && (state[i] == 1)) ||
                   ((grad[i] > 0) && (state[i] == 2))
               done = false
               break
           end
        end
        done && return x
        #
        # Not optimal, so we need to free up some vars at bounds.
        #
        newi = 0
        newg = 0.0
        for i in 1:length(atbound)
            if atbound[i]
                #
                # Don't free up a variable that was just locked at its bound.
                #
                i == criti && continue
                #
                # Look for the locked variable with the biggest gradient.
                #
                if (grad[i] > 0) && (state[i] == 2) && (abs(grad[i]) > newg)
                    newi = i
                    newg = abs(grad[i])
                end
                if (grad[i] < 0) && (state[i] == 1) && (abs(grad[i]) > newg)
                    newi = i
                    newg = abs(grad[i])
                end
            end
        end
        #
        # Free the locked variable with the biggest gradient if there is one.
        #
        if (newi != 0)
            atbound[newi] = false
            state[newi] = 0
            between[newi] = true
        end
        #
        # Make sure the projected problem is nontrivial.
        #
        if (length(between)==0)
          # println("Empty projected problem")
          continue
        end
        #
        # Construct the new projected problem.
        #
        Aproj = A[:, between]
        An = A[:, atbound]
        if any(atbound)
            bproj = b - An*x[atbound]
        else
            bproj = b
        end
        #
        # Solve the projected problem.
        #
        z = Aproj\bproj
        xnew = copy(x)
        xnew[between] = z

        if newi != 0 && ((xnew[newi] <= l[newi] && x[newi] == l[newi]) ||
                (xnew[newi] >= u[newi] && x[newi] == u[newi]))
            #
            # Ooops- the freed variable wants to go beyond its bound.  Lock it down
            # and look for some other variable.
            #
            oopslist[newi] = true
            if xnew[newi] <= l[newi] && state[newi] == 1
                state[newi] = 1
                x[newi] = l[newi]
            end
            if xnew[newi] >= u[newi] && state[newi] == 2
                state[newi] = 2
                x[newi] = u[newi]
            end
            atbound[newi] = true
            between[newi] = false
            continue
        end
        #
        # We've got a good variable freed up.  Reset the oopslist.
        #
        oopslist = falses(n)
        #
        # Move as far as possible towards the optimal solution to the projected
        # problem.
        #
        alpha = 1.
        for i in 1:length(atbound)
            if between[i]
                if (xnew[i] > u[i])
                    newalpha = min(alpha, (u[i]-x[i])/(xnew[i]-x[i]))
                    if (newalpha < alpha)
                        criti = i
                        crits = 2
                        alpha = newalpha
                    end
                end
                if (xnew[i] < l[i])
                    newalpha = min(alpha, (l[i]-x[i])/(xnew[i]-x[i]))
                    if (newalpha < alpha)
                        criti = i
                        crits = 1
                        alpha = newalpha
                    end
                end
            end
        end
        #
        # Take the step.
        #
        x = x + alpha*(xnew-x)
        #
        # Update the state of variables.
        #
        if (alpha < 1)
            between[criti] = false
            atbound[criti] = true
            state[criti] = crits
        end

        for i in 1:n
            if (x[i] >= u[i])
                x[i] = u[i]
                state[i] = 2
                between[i] = false
                atbound[i] = true
            end
            if (x[i] <= l[i])
                x[i] = l[i]
                state[i] = 1
                between[i] = false
                atbound[i] = true
            end
        end
    end
end

function bvls{m,n}(A::SMatrix{m,n}, b::SVector{m}, l::SVector{n}, u::SVector{n})
    SVector{n}(bvls(Matrix(A), Vector(b), Vector(l), Vector(u)))
end


text:
using PackageExample, Test

@testset "foo" begin
    x, y = 5, 7
    @test foo(x, y) == 7
    x = "blah"
    @test_throws MethodError foo(x, y)
end

@testset "bar" begin
    z = 4.
    @test bar(z) == 1.
end


text:
################################################################################
#
#  Ringmorphisms of orders onto finite fields
#
################################################################################

# NfOrd -> FqNmod
mutable struct NfOrdToFqNmodMor <: Map{NfOrd, FqNmodFiniteField, HeckeMap, NfOrdToFqNmodMor}
  header::MapHeader{NfOrd, FqNmodFiniteField}
  poly_of_the_field::nmod_poly
  P::NfOrdIdl

  function NfOrdToFqNmodMor()
    r = new()
    r.header = MapHeader{NfOrd, FqNmodFiniteField}()
    return r
  end

  function NfOrdToFqNmodMor(O::NfOrd, F::FqNmodFiniteField, g::nmod_poly)
    # assume that F = F_p[X]/(g) and g is a factor of f mod p

    z = new()
    d = degree(nf(O))
    p = characteristic(F)
    a = gen(nf(O))
    tmp_nmod_poly = parent(g)()
    z.poly_of_the_field = g

    powers = Vector{nf_elem}(d)

    powers[1] = a
    for i in 2:d
      powers[i] = powers[i - 1] * a
    end

    function _image(x::NfOrdElem)
      u = F()
      gg = parent(nf(O).pol)(elem_in_nf(x))::fmpq_poly
      fmpq_poly_to_nmod_poly_raw!(tmp_nmod_poly, gg)
      ccall((:nmod_poly_rem, :libflint), Void,
            (Ptr{nmod_poly}, Ptr{nmod_poly}, Ptr{nmod_poly}, Ptr{Void}),
            &tmp_nmod_poly, &tmp_nmod_poly, &g, pointer_from_objref(F)+sizeof(fmpz))
      ccall((:fq_nmod_set, :libflint), Void,
            (Ptr{fq_nmod}, Ptr{nmod_poly}, Ptr{FqNmodFiniteField}),
            &u, &tmp_nmod_poly, &F)
      ccall((:fq_nmod_reduce, :libflint), Void,
            (Ptr{fq_nmod}, Ptr{FqNmodFiniteField}), &u, &F)
      return u
    end

    # The lift is even simpler!
    function _preimage(y::fq_nmod)
      zz = O()
      zz.elem_in_nf = nf(O)(coeff(y, 0))
      for i in 2:d
        add!(zz.elem_in_nf, zz.elem_in_nf, powers[i - 1] * coeff(y, i - 1))
      end
      _mod!(zz.elem_in_nf, p)
      return zz
    end

    z.header = MapHeader{NfOrd, FqNmodFiniteField}(O, F, _image, _preimage)
    return z
  end

  function NfOrdToFqNmodMor(O::NfOrd, P::NfOrdIdl)
    z = NfOrdToFqNmodMor()
    z.P = P
    p = minimum(P)
    a, g, b = get_residue_field_data(P)
    psmall = Int(p)
    R = ResidueRing(FlintZZ, psmall, cached = false)
    Rx, x = PolynomialRing(R, "_\$", cached = false)
    F = FqNmodFiniteField(Rx(g), Symbol("_\$"), false)
    d = degree(g)
    n = degree(O)
    imageofbasis = Vector{fq_nmod}(n)
    powers = Vector{nf_elem}(d)
    c = Rx()
    for i in 1:n
      ib = F()
      @assert d == cols(b[i])
      for j in 1:d
        setcoeff!(c, j - 1, b[i][1, j])
      end
      ccall((:fq_nmod_set, :libflint), Void, (Ref{fq_nmod}, Ref{nmod_poly}, Ref{FqNmodFiniteField}), ib, c, F)
      imageofbasis[i] = ib
    end

    powers[1] = a.elem_in_nf
    for i in 2:d
      powers[i] = powers[i - 1] * a.elem_in_nf
    end

    tempF = F()

    function _image(x::NfOrdElem)
      v = elem_in_basis(x, Val{false})
      zz = zero(F)
      for i in 1:n
        ccall((:fq_nmod_mul_fmpz, :libflint), Void,
              (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fmpz}, Ref{FqNmodFiniteField}),
              tempF, imageofbasis[i], v[i], F)
        add!(zz, zz, tempF)
      end
      return zz
    end

    function _preimage(y::fq_nmod)
      zz = O()
      zz.elem_in_nf = nf(O)(coeff(y, 0))
      for i in 2:d
        add!(zz.elem_in_nf, zz.elem_in_nf, powers[i - 1] * coeff(y, i - 1))
      end
      _mod!(zz.elem_in_nf, p)
      return zz
    end

    z.header = MapHeader{NfOrd, FqNmodFiniteField}(O, F, _image, _preimage)
    return z
  end
end

_mod!(x::nf_elem, y::Integer) = _mod!(x, fmpz(y))

function _mod!(x::nf_elem, y::fmpz)
  zden = denominator(x)
  zden2 = denominator(x)
  mul!(x, x, zden)
  mul!(zden, zden, y)
  mod!(x, zden)
  divexact!(x, x, zden2)
  return x
end

mutable struct NfOrdQuoMap <: Map{NfOrd, NfOrdQuoRing, HeckeMap, NfOrdQuoMap}
  header::MapHeader{NfOrd, NfOrdQuoRing}

  function NfOrdQuoMap(O::NfOrd, Q::NfOrdQuoRing)
    z = new()

    _image = function (x::NfOrdElem)
      return Q(x)
    end

    _preimage = function (x::NfOrdQuoRingElem)
      return x.elem
    end

    z.header = MapHeader(O, Q, _image, _preimage)
    return z
  end
end

function Mor(O::NfOrd, F::FqNmodFiniteField, y::fq_nmod)
  return NfOrdToFqNmodMor(O, F, y)
end

function Mor(O::NfOrd, F::FqNmodFiniteField, h::nmod_poly)
  return NfOrdToFqNmodMor(O, F, h)
end

function evaluate(f::fmpz_poly, r::fq_nmod)
  #Horner - stolen from Claus

  if length(f) == 0
    return parent(r)()
  end

  l = f.length-1
  s = parent(r)(coeff(f, l))
  for i =l-1:-1:0
    s = s*r + parent(r)(coeff(f, i))
  end
  return s
end

function evaluate!(z::fq_nmod, f::fmpz_poly, r::fq_nmod)
  #Horner - stolen from Claus

  zero!(z)

  if length(f) == 0
    return z
  end

  l = f.length-1
  set!(z, parent(r)(coeff(f, l)))
  #s = parent(r)(coeff(f, l))
  for i =l-1:-1:0
    mul!(z, z, r)
    add!(z, z, parent(r)(coeff(f, i)))
    #s = s*r + parent(r)(coeff(f, i))
  end
  return z
end

function Mor(K::AnticNumberField, L::AnticNumberField, y::nf_elem)
  z = NfToNfMor(K, L, y)
  return z
end

function _get_coeff_raw(x::fq_nmod, i::Int)
  u = ccall((:nmod_poly_get_coeff_ui, :libflint), UInt, (Ptr{fq_nmod}, Int), &x, i)
  return u
end

function _get_coeff_raw(x::fq, i::Int)
  t = FlintIntegerRing()
  ccall((:fmpz_poly_get_coeff_fmpz, :libflint), Void, (Ptr{fmpz}, Ptr{fq}, Int), &t, &x, i)
  return t
end

# this is super slow
# improve
function (f::NfOrdQuoMap)(I::NfOrdIdl)
  O = domain(f)
  Q = codomain(f)
  B = Q.ideal + I
  nB = norm(B)
  b = basis(B, Val{false})

  z = O()

  nQ = norm(Q.ideal)
  OnQ = ideal(O, nQ)
  range1nQ2 = fmpz(1):nQ^2

  while true
    z = rand!(z, b, range1nQ2)
    #z = sum(rand(range1nQ2) * b[i] for i in 1:degree(O))
    if norm(ideal(O, z) + OnQ) == nB
      break
    end
  end

  return Q(z)
end

function (f::NfOrdQuoMap)(p::PolyElem{NfOrdElem})
  F = codomain(f)
  Fx,_ = PolynomialRing(F, "_\$", cached = false)

  ar = NfOrdElem[ coeff(p, i) for i in 0:degree(p) ]

  z = Fx(map(f, ar))

  return z
end

# Assume A is mxd, B is mxl and there is a unique X of size lxd
# with A = B * X
# this function will find it!
function _solve_unique(A::nmod_mat, B::nmod_mat)
  X = zero_matrix(base_ring(A), cols(B), rows(A))

  #println("solving\n $A \n = $B * X")
  r, per, L, U = lufact(B) # P*M1 = L*U

  inv!(per)

  @assert B == per*L*U
  Ap = inv(per)*A
  Y = parent(A)()

  #println("first solve\n $Ap = $L * Y")

  for i in 1:cols(Y)
    for j in 1:rows(Y)
      s = Ap[j, i]
      for k in 1:j-1
        s = s - Y[k, i]*L[j, k]
      end
      Y[j, i] = s
    end
  end

  @assert Ap == L*Y

  #println("solving \n $Y \n = $U * X")

  YY = sub(Y, 1:r, 1:cols(Y))
  UU = sub(U, 1:r, 1:r)
  X = _inv(UU)*YY

  @assert Y == U * X

  @assert B*X == A
  return X
end

function _solve_unique(A::Generic.Mat{Generic.Res{fmpz}}, B::Generic.Mat{Generic.Res{fmpz}})
  X = zero_matrix(base_ring(A), cols(B), rows(A))

  #println("solving\n $A \n = $B * X")
  r, per, L, U = _lufact(B) # P*M1 = L*U
  inv!(per)

  @assert B == per*L*U
  Ap = inv(per)*A
  Y = zero_matrix(base_ring(A), rows(A), cols(A))

  #println("first solve\n $Ap = $L * Y")

  for i in 1:cols(Y)
    for j in 1:rows(Y)
      s = Ap[j, i]
      for k in 1:j-1
        s = s - Y[k, i]*L[j, k]
      end
      Y[j, i] = s
    end
  end

  @assert Ap == L*Y

  #println("solving \n $Y \n = $U * X")

  YY = sub(Y, 1:r, 1:cols(Y))
  UU = sub(U, 1:r, 1:r)
  X = _inv(UU)*YY

  @assert Y == U * X

  @assert B*X == A
  return X
end

function _lufact!(P::Generic.perm, A::Nemo.MatElem{T}) where {T}
   m = rows(A)
   n = cols(A)
   rank = 0
   r = 1
   c = 1
   R = base_ring(A)
   t = R()
   while r <= m && c <= n
      if A[r, c] == 0
         i = r + 1
         while i <= m
            if !iszero(A[i, c])
               for j = 1:n
                  A[i, j], A[r, j] = A[r, j], A[i, j]
               end
               P[r], P[i] = P[i], P[r]
               break
            end
            i += 1
         end
         if i > m
            c += 1
            continue
         end
      end
      rank += 1
      d = -inv(A[r, c])
      for i = r + 1:m
         q = A[i, c]*d
         for j = c + 1:n
            t = mul!(t, A[r, j], q)
            A[i, j] = addeq!(A[i, j], t)
         end
         A[i, c] = R()
         A[i, rank] = -q
      end
      r += 1
      c += 1
   end
   inv!(P)
   return rank
end

function _lufact(A::Nemo.MatElem{T}, P = PermGroup(rows(A))) where {T}
   m = rows(A)
   n = cols(A)
   P.n != m && error("Permutation does not match matrix")
   p = P()
   R = base_ring(A)
   U = deepcopy(A)
   L = similar(A, m, m)
   rank = _lufact!(p, U)
   for i = 1:m
      for j = 1:n
         if i > j
            L[i, j] = U[i, j]
            U[i, j] = R()
         elseif i == j
            L[i, j] = R(1)
         elseif j <= m
            L[i, j] = R()
         end
      end
   end
   return rank, p, L, U
end

mutable struct NfOrdToFqMor <: Map{NfOrd, FqFiniteField, HeckeMap, NfOrdToFqMor}
  header::MapHeader{NfOrd, FqFiniteField}
  poly_of_the_field::fmpz_mod_poly
  P::NfOrdIdl
  fastpath::Bool
  # Some temporary variables
  tmp_fmpz_mod_poly::fmpz_mod_poly
  t_fmpz_poly::fmpz_poly
  t_fmpz::fmpz
  a::nf_elem

  function NfOrdToFqMor()
    z = new()
    return z
  end

  function NfOrdToFqMor(O::NfOrd, F::FqFiniteField, g::fmpz_mod_poly)
    # assume that F = F_p[X]/(g) and g is a factor of f mod p

    z = new()
    z.fastpath = true
    p = characteristic(F)
    z.tmp_fmpz_mod_poly = parent(g)()
    z.t_fmpz_poly = fmpz_poly()
    z.t_fmpz = fmpz()

    z.a = gen(nf(O))
    z.poly_of_the_field = g

    z.header = MapHeader{NfOrd, FqFiniteField}(O, F)# _image, _preimage)

    return z
  end

#  function NfOrdToFqMor(O::NfOrd, P::NfOrdIdl)
#    z = new()
#
#    z.fastpath = false
#
#    p = abs(minimum(P))
#
#    R = ResidueRing(FlintZZ, p, cached=false)
#
#    OP = quoringalg(O, P, p)
#    x = quoelem(OP, zero(O))
#    f = PolynomialRing(R, "x", cached = false)[1]()
#    d = length(OP.basis)
#
#    while true
#      r = rand(0:p-1, d)
#      x = quoelem(OP, dot([ x for x in OP.basis], r))
#      f = minpoly(x)
#      if degree(f) == d
#        break
#      end
#    end
#
#    F = FqFiniteField(f, Symbol("_\$"), false)
#
#    M2 = zero_matrix(R, degree(O), d)
#
#    for i in 1:d
#      coords = elem_in_basis((x^(i-1)).elem)
#      for j in 1:degree(O)
#        M2[j, i] = coords[j]
#      end
#    end
#
#    M3 = zero_matrix(R, degree(O), degree(O))
#
#    for i in 1:degree(O)
#      coords = elem_in_basis(mod(basis(O)[i], OP.ideal))
#      for j in 1:degree(O)
#        M3[j, i] = coords[j]
#      end
#    end
#
#    X = _solve_unique(M3, M2)
#
#    #for i in 1:degree(O)
#    #  @assert quoelem(OP, basis(O)[i]) == quoelem(OP, dot([(x^j).elem for j in 0:d-1], _lift([ X[j, i] for j in 1:d ])))
#    #end
#
#    function _image(y::NfOrdElem)
#      co = zero_matrix(R, degree(O), 1)
#      coeff = elem_in_basis(mod(y, P))
#
#      for i in 1:degree(O)
#        co[i, 1] = coeff[i]
#      end
#
#      co = X*co
#
#      z = F(lift(co[1, 1])) # totally inefficient
#
#      for i in 2:rows(co)
#        z = z  + F(lift(co[i, 1]))*gen(F)^(i-1)
#      end
#
#      return z
#    end
#
#    function _preimage(y::fq)
#      zz = nf(O)()
#
#      for i in 0:degree(F)-1
#        zz = zz + coeff(y, i)*(x.elem.elem_in_nf)^i
#      end
#
#      return mod(O(zz, false), P)
#    end
#
#    z.header = MapHeader{NfOrd, FqFiniteField}(O, F, _image, _preimage)
#    z.P = P
#
#    return z
#  end
#
end

function NfOrdToFqMor(O::NfOrd, P::NfOrdIdl)#, g::fmpz_poly, a::NfOrdElem, b::Vector{fmpz_mat})
  z = NfOrdToFqMor()
  z.fastpath = false
  z.P = P
  a, g, b = get_residue_field_data(P)
  p = minimum(P)
  R = ResidueRing(FlintZZ, p, cached = false)
  Rx, x = PolynomialRing(R, "_\$", cached = false)
  F = FqFiniteField(Rx(g), Symbol("_\$"), false)
  d = degree(g)
  n = degree(O)
  imageofbasis = Vector{fq}(n)
  powers = Vector{nf_elem}(d)
  c = Rx()

  for i in 1:n
    ib = F()
    @assert d == cols(b[i])
    for j in 1:d
      setcoeff!(c, j - 1, b[i][1, j])
    end
    ccall((:fq_set, :libflint), Void, (Ref{fq}, Ref{fmpz_mod_poly}, Ref{FqFiniteField}), ib, c, F)
    imageofbasis[i] = ib
  end

  powers[1] = a.elem_in_nf
  for i in 2:d
    powers[i] = powers[i - 1] * a.elem_in_nf
  end

  tempF = F()

  function _image(x::NfOrdElem)
    v = elem_in_basis(x, Val{false})
    zz = zero(F)
    for i in 1:n
      ccall((:fq_mul_fmpz, :libflint), Void,
            (Ref{fq}, Ref{fq}, Ref{fmpz}, Ref{FqFiniteField}),
            tempF, imageofbasis[i], v[i], F)
      add!(zz, zz, tempF)
    end
    return zz
  end

  function _preimage(y::fq)
    zz = O()
    zz.elem_in_nf = nf(O)(coeff(y, 0))
    for i in 2:d
      add!(zz.elem_in_nf, zz.elem_in_nf, powers[i - 1] * coeff(y, i - 1))
    end
    _mod!(zz.elem_in_nf, p)
    return zz
  end

  z.header = MapHeader{NfOrd, FqFiniteField}(O, F, _image, _preimage)
  return z
end

function image(f::NfOrdToFqMor, x::NfOrdElem)
  if f.fastpath
    F = codomain(f)
    O = domain(f)
    u = F()
    gg = parent(nf(O).pol)(elem_in_nf(x))::fmpq_poly
    fmpq_poly_to_fmpz_mod_poly_raw!(f.tmp_fmpz_mod_poly, gg, f.t_fmpz_poly, f.t_fmpz)
    ccall((:fmpz_mod_poly_rem, :libflint), Void, (Ptr{fmpz_mod_poly}, Ptr{fmpz_mod_poly}, Ptr{fmpz_mod_poly}), &f.tmp_fmpz_mod_poly, &f.tmp_fmpz_mod_poly, &f.poly_of_the_field)
    ccall((:fq_set, :libflint), Void, (Ptr{fq}, Ptr{fmpz_mod_poly}, Ptr{FqFiniteField}), &u, &f.tmp_fmpz_mod_poly, &F)
    ccall((:fq_reduce, :libflint), Void, (Ptr{fq}, Ptr{FqFiniteField}), &u, &F)
    return u
  else
    return f.header.image(x)
  end
end

function preimage(f::NfOrdToFqMor, x::fq)
  if f.fastpath
    O = domain(f)
    F = codomain(f)
    zz = nf(O)()

    a = f.a
    # TODO: Do something more clever here
    for i in 0:degree(F)-1
      zz = zz + coeff(x, i)*a^i
    end

    return O(zz, false)
  else
    @assert isdefined(f.header, :preimage)
    return f.header.preimage(x)
  end
end

function Mor(O::NfOrd, F::FqFiniteField, h::fmpz_mod_poly)
  return NfOrdToFqMor(O, F, h)
end

function sub(M::Nemo.MatElem{T}, r::UnitRange{<:Integer}, c::UnitRange{<:Integer}) where {T}
  z = similar(M, length(r), length(c))
  for i in 1:length(r)
    for j in 1:length(c)
      z[i, j] = M[r[i], c[j]]
    end
  end
  return z
end

_inv(a::nmod_mat) = inv(a)

function _inv(a::MatElem{Generic.Res{fmpz}})
  b, d = inv(a)
  return divexact(b, d)
end

mutable struct NfToFqNmodMor <: Map{AnticNumberField, FqNmodFiniteField, HeckeMap, NfToFqNmodMor}
  header::MapHeader{AnticNumberField, FqNmodFiniteField}

  function NfToFqNmodMor()
    r = new()
    r.header = MapHeader{AnticNumberField, FqNmodFiniteField}()
    return r
  end
end

mutable struct NfToFqMor <: Map{AnticNumberField, FqFiniteField, HeckeMap, NfToFqMor}
  header::MapHeader{AnticNumberField, FqFiniteField}

  function NfToFqMor()
    r = new()
    r.header = MapHeader{AnticNumberField, FqFiniteField}()
    return r
  end
end

function extend(f::Union{NfOrdToFqNmodMor, NfOrdToFqMor}, K::AnticNumberField)
  nf(domain(f)) != K && error("Number field is not the number field of the order")

  if f isa NfOrdToFqNmodMor
    z = NfToFqNmodMor()
  elseif f isa NfOrdToFqMor
    z = NfToFqMor()
  end

  z.header.domain = K
  z.header.codomain = f.header.codomain

  p = characteristic(z.header.codomain)
  Zx = PolynomialRing(FlintIntegerRing(), "x", cached = false)[1]
  y = f(NfOrdElem(domain(f), gen(K)))
  pia = anti_uniformizer(f.P)
  O = domain(f)
  P = f.P

  #function _image(x::nf_elem)
  #  g = parent(K.pol)(x)
  #  u = inv(z.header.codomain(denominator(g)))

  #  g = Zx(denominator(g)*g)
  #  return u*evaluate(g, y)
  #end
  function _image(x::nf_elem)
    m = denominator(x, domain(f))
    l = valuation(m, P)
    if l == 0
      return f(O(m*x))//f(O(m))
    else
      return f(O(pia^l * m * x))//f(O(pia^l * m))
    end
  end

  function _preimage(x::Union{fq, fq_nmod})
    return elem_in_nf(preimage(f, x))
  end

  z.header.image = _image
  z.header.preimage = _preimage

  return z
end

function (f::NfOrdToFqNmodMor)(p::PolyElem{NfOrdElem})
  F = codomain(f)
  Fx,_ = PolynomialRing(F, "_\$", cached = false)

  ar = NfOrdElem[ coeff(p, i) for i in 0:degree(p) ]

  z = Fx(map(f, ar))

  return z
end

function (f::NfOrdToFqMor)(p::PolyElem{NfOrdElem})
  F = codomain(f)
  Fx,_ = PolynomialRing(F, "_\$", cached = false)

  ar = NfOrdElem[ coeff(p, i) for i in 0:degree(p) ]

  z = Fx(map(f, ar))

  return z
end

mutable struct NfOrdToAlgAssMor{T} <: Map{NfOrd, AlgAss{T}, HeckeMap, NfOrdToAlgAssMor}
  header::MapHeader

  function NfOrdToAlgAssMor{T}(O::NfOrd, A::AlgAss{T}, _image::Function, _preimage::Function) where T
    z = new{T}()
    z.header = MapHeader(O, A, _image, _preimage)
    return z
  end
end


text:
# Newton-like methods for NEPs

using LinearAlgebra
using Printf
using Random

export newton
export resinv
export augnewton
export quasinewton
export newtonqr
export implicitdet

#############################################################################
"""
    λ,v = newton([eltype],nep::NEP;[errmeasure,][tol,][maxit,][λ,][v,][c,][logger,][armijo_factor=1,][armijo_max])

Applies Newton-Raphsons method on the system of
nonlinear equations with `n+1` unknowns:
```math
M(λ)v=0
```
```math
c^Hv-1=0
```
The vector `c` is the
orthogonalization vector.  If `c=0` the current approximation will be used for the orthogonalization. See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> using LinearAlgebra
julia> nep=nep_gallery("dep0");
julia> λ,v=newton(nep);
julia> minimum(svdvals(compute_Mder(nep,λ)))
1.9997125567227177e-16
```

# References
* Nichtlineare Behandlung von Eigenwertaufgaben, Z. Angew. Math. Mech. 30 (1950) 281-282.
* A. Ruhe, Algorithms for the nonlinear eigenvalue problem, SIAM J. Numer. Anal. 10 (1973) 674-689

"""
    newton(nep::NEP;params...)=newton(ComplexF64,nep;params...)
    function newton(::Type{T},
                    nep::NEP;
                    errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
                    tol::Real=eps(real(T))*100,
                    maxit::Int=10,
                    λ::Number=zero(T),
                    v::Vector=randn(size(nep,1)),
                    c::Vector=v,
                    logger=0,
                    armijo_factor::Real=1,
                    armijo_max::Int=5) where {T<:Number}

        @parse_logger_param!(logger)

        # Ensure types λ and v are of type T
        λ=T(λ)
        v=Vector{T}(v)
        c=Vector{T}(c)

        err=Inf;
        v[:] = v/dot(c,v);

        for k=1:maxit
            err=estimate_error(errmeasure,λ,v)

            push_iteration_info!(logger,k,err=err,λ=λ,v=v,continues=true);
            if (err< tol)
                push_info!(logger,"")
                return (λ,v)
            end

            # Compute NEP matrix and derivative
            M = compute_Mder(nep,λ)
            Md = compute_Mder(nep,λ,1)

            # Create jacobian
            J = [M Md*v; c' 0];
            F = [M*v; c'*v-1];

            # Compute update
            delta=-J\F;  # Hardcoded backslash

            Δv=Vector{T}(delta[1:size(nep,1)]);
            Δλ=T(delta[size(nep,1)+1]);

            (Δλ,Δv,j,scaling)=armijo_rule(nep,errmeasure,err,
                                          λ,v,Δλ,Δv,real(T(armijo_factor)),armijo_max)
            if (j>0)
                push_info!(logger," Armijo scaling=$scaling")
            else
                push_info!(logger,"")
            end

            # Update eigenvalue and eigvec
            v[:] += Δv
            λ = λ+Δλ
        end

        msg="Number of iterations exceeded. maxit=$(maxit)."
        throw(NoConvergenceException(λ,v,err,msg))
    end

    ############################################################################
"""
    λ,v = resinv([eltype],nep::NEP;[errmeasure,][tol,][maxit,][λ,][v,][c,][logger,][armijo_factor=1,][armijo_max,][linsolvecreator])

Applies residual inverse iteration method for nonlinear eigenvalue problems.
The kwarg `linsolvecreator`
is a function which specifies how the linear system is created.
The function calls `compute_rf` for the computation
of the Rayleigh functional.
See [`augnewton`](@ref) for other parameters.

# Example
The example shows how to specify if the method should run in real
or complex mode (or any other `Number` type).
```julia-repl
julia> nep=nep_gallery("qdep0");
julia> λ,v=resinv(nep,λ=-2,v=ones(size(nep,1)))
julia> typeof(λ)
Complex{Float64}
julia> norm(compute_Mlincomb(nep,λ,v))
6.688224435370382e-12
julia> λ,v=resinv(Float64,nep,λ=-2,v=ones(size(nep,1)))
julia> typeof(λ)
Float64
julia> norm(compute_Mlincomb(nep,λ,v))
5.939894690000396e-12
```

# References
*  A. Neumaier, Residual inverse iteration for the nonlinear eigenvalue problem, SIAM J. Numer. Anal. 22 (1985) 914-923

"""
    resinv(nep::NEP;params...)=resinv(ComplexF64,nep;params...)
    function resinv(::Type{T},
                    nep::NEP;
                    errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
                    tol::Real=eps(real(T))*100,
                    maxit::Int=100,
                    λ::Number=zero(T),
                    v::Vector=randn(real(T),size(nep,1)),
                    c::Vector=v,
                    logger=0,
                    inner_solver= @default_compute_rf_inner_solver(nep),
                    linsolvercreator=DefaultLinSolverCreator(),
                    armijo_factor::Real=1,
                    armijo_max::Int=5) where T

        @parse_logger_param!(logger)

        # Ensure types λ and v are of type T
        λ::T=T(λ)
        v=Vector{T}(v)
        c=Vector{T}(c)
        n=size(v,1);

        push_info!(logger,"Precomputing linsolver")
        local linsolver::LinSolver=create_linsolver(linsolvercreator,nep,λ)

        # If c is zero vector we take eigvec approx as left vector in
        # generalized Rayleigh functional
        use_v_as_rf_vector=false;
        if (norm(c)==0)
            use_v_as_rf_vector=true;
        end

        push_info!(logger,2,
                   "use_v_as_rf_vector=$use_v_as_rf_vector");

        σ::T=λ;
        err=Inf;

        for k=1:maxit
            # Normalize
            v[:] = v/norm(v);

            err=estimate_error(errmeasure,λ,v)

            if (use_v_as_rf_vector)
                c[:]=v;
            end

            push_iteration_info!(logger,k,err=err,λ=λ,v=v,continues=true);

            if (err< tol)
                push_info!(logger,"")
                return (λ,v)
            end

            # Compute eigenvalue update
            λ_vec = compute_rf(T, nep, v, inner_solver, y=c, λ=λ, target=σ)
            local λ1::T = closest_to(λ_vec,  λ)
            Δλ=λ1-λ

            # Compute eigenvector update
            Δv = -lin_solve(linsolver,compute_Mlincomb(nep,λ1,reshape(v,n,1))) #M*v);

            (Δλ,Δv,j,scaling)=armijo_rule(nep,errmeasure,err,
                                          λ,v,Δλ,Δv,real(T(armijo_factor)),armijo_max)
            if (j>0)
                push_info!(logger," Armijo scaling=$scaling")
            else
                push_info!(logger,"")
            end

            # Update the eigenpair
            λ+=Δλ
            v[:] += Δv;

        end

        msg="Number of iterations exceeded. maxit=$(maxit)."
        throw(NoConvergenceException(λ,v,err,msg))
    end

    # New augnewton
"""
    augnewton([eltype], nep::NEP; [errmeasure,][tol,][maxit,][λ,][v,][c,][logger,][linsolvercreator,][armijo_factor,][armijo_max])

Run the augmented Newton method. The method is equivalent to `newton()`
in exact arithmetic,  but works only with operations on vectors of
length `n`.

The following keyword arguments are in common for many NEP-solvers:

* `logger` is either a [`Logger`](@ref) object or an `Int`. If it is an `Int`, a `PrintLogger(logger)` will be instantiated. `logger=0` prints nothing, `logger=1` prints more, etc.

* `errmeasure` determines how error is measured. It is either a function handle or an object of the type `Errmeasure`.  If it is a function handle, it should take `(λ,v)` as input and return a real scalar (the error). See [`Errmeasure`](@ref) and [`ErrmeasureType`](@ref) for further description.

* `tol` is a scalar which determines termination. If `errmeasure` is less than `tol` the eigenpair is marked as converged.

* The scalar `λ` and the vector `v` are starting approximations.

* `maxit` determines the maximum number of iterations. The error `NoConvergenceException` is thrown if this is exceeded.

*  The `linsolvecreator` specifies how the linear system should be solved. See [`LinSolver`](@ref) for further information.

* `armijo_factor` specifies if an Armijo rule should be applied, and its value specifies the scaling factor of the step length (per reduction step). The variable `armijo_max` specifies the maximum number of step length reductions.

# Example
This illustrates the equivalence between `newton` and `augnewton`.
```julia-repl
julia> nep=nep_gallery("dep1")
julia> λ1,v1=newton(nep,maxit=20,v=ones(size(nep,1)),λ=0)
julia> λ2,v2=augnewton(nep,maxit=20,v=ones(size(nep,1)),λ=0)
julia> λ1-λ2
0.0 + 0.0im
```
# References
* Nichtlineare Behandlung von Eigenwertaufgaben, Z. Angew. Math. Mech. 30 (1950) 281-282.
* A. Ruhe, Algorithms for the nonlinear eigenvalue problem, SIAM J. Numer. Anal. 10 (1973) 674-689
"""
    augnewton(nep::NEP;kwargs...)=augnewton(ComplexF64,nep::NEP;kwargs...)
    function augnewton(::Type{T},
                       nep::NEP;
                       errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
                       tol::Real=eps(real(T))*100,
                       maxit::Int=30,
                       λ::Number=zero(T),
                       v::Vector=randn(real(T),size(nep,1)),
                       c::Vector=v,
                       logger=0,
                       linsolvercreator=DefaultLinSolverCreator(),
                       armijo_factor::Real=one(real(T)),
                       armijo_max::Int=5) where {T<:Number}

        @parse_logger_param!(logger)

        # Ensure types λ and v are of type T
        λ=T(λ)
        v=Vector{T}(v)
        c=Vector{T}(c)

        err=Inf;
        # If c is zero vector we take eigvec approx as normalization vector
        use_v_as_normalization_vector=false;
        if norm(c) == 0
            use_v_as_normalization_vector=true;
            c[:] = v / norm(v)^2
        end
        v[:] = v/dot(c,v);
        local linsolver::LinSolver
        local tempvec = Vector{T}(undef, size(nep,1))

        push_info!(logger,2,
                   "use_v_as_normalization_vector=$use_v_as_normalization_vector");

        for k=1:maxit
            err=estimate_error(errmeasure,λ,v)
            push_iteration_info!(logger,k,err=err,λ=λ,v=v,continues=true);
            if (err< tol)
                push_info!(logger,"")
                return (λ,v)
            end
            # tempvec =  (M(λ_k)^{-1})*M'(λ_k)*v_k
            # α = 1/(c'*(M(λ_k)^{-1})*M'(λ_k)*v_k);

            z::AbstractVector=compute_Mlincomb(nep,λ,v,[T(1.0)],1)

            linsolver = create_linsolver(linsolvercreator,nep,λ)
            tempvec[:] = Vector{T}(lin_solve(linsolver, z, tol=tol));

            if (use_v_as_normalization_vector)
                c[:] = v /norm(v)^2
            end
            α = T(1)/ dot(c,tempvec);

            Δλ=-α
            Δv=α*tempvec-v;

            (Δλ,Δv,j,scaling)=armijo_rule(nep,errmeasure,err,
                                          λ,v,Δλ,Δv,real(T(armijo_factor)),armijo_max)

            if (j>0)
                push_info!(logger," Armijo scaling=$scaling")
            else
                push_info!(logger,"")
            end

            λ+=Δλ
            v[:]+=Δv

        end

        msg="Number of iterations exceeded. maxit=$(maxit)."
        throw(NoConvergenceException(λ,v,err,msg))
    end

"""
    quasinewton([T=ComplexF64],nep,[errmeasure,][tol,][maxit,][λ,][v][ws][logger][linsolvercreator,][armijo_factor,][armijo_max])

An implementation of the quasi-Newton approach referred to as quasi-Newton 2 in the reference.
The method involves one linear system solve per iteration corresponding with the
matrix ``M(λ)``, where ``λ`` is constant.
The vector `ws` is a representation of the normalization, in the sense that ``c^T=w_s^TM(λ)``,
where all iterates satisfy ``c^Tx_i=1``.
See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> nep=nep_gallery("pep0")
julia> λ,v=quasinewton(nep,λ=1.0,v=ones(size(nep,1)));
julia> norm(compute_Mlincomb(nep,λ,v))/norm(v)
5.448264607410413e-12
```

# References
* Jarlebring, Koskela, Mele, Disguised and new Quasi-Newton methods for nonlinear eigenvalue problems, Numer. Algorithms, 79:311-335, 2018. [preprint](https://arxiv.org/abs/1702.08492)
"""
    quasinewton(nep::NEP;params...)=quasinewton(ComplexF64,nep;params...)
    function quasinewton(::Type{T},
                         nep::NEP;
                         errmeasure = DefaultErrmeasure(nep),
                         tol::Real=eps(real(T))*100,
                         maxit::Int=100,
                         λ::Number=zero(T),
                         v::Vector=randn(real(T),size(nep,1)),
                         ws::Vector=v,
                         logger=0,
                         linsolvercreator=DefaultLinSolverCreator(),
                         armijo_factor::Real=1,
                         armijo_max::Int=5) where T

        @parse_logger_param!(logger)

        # Ensure types λ and v are of type T
        λ=T(λ)
        v=Vector{T}(v)
        ws=Vector{T}(ws) # Left vector such that c'=w'M(λ) where c normalization

        n = size(nep,1)
        u = zeros(T,n)
        w = zeros(T,n)

        err=Inf;

        local linsolver::LinSolver;
        push_info!(logger,"Precomputing linsolver")
        linsolver = create_linsolver(linsolvercreator,nep,λ)

        for k=1:maxit
            err=estimate_error(errmeasure,λ,v)
            push_iteration_info!(logger,k,err=err,λ=λ,v=v,continues=true);
            if (err< tol)
                push_info!(logger,"")
                return (λ,v)
            end

            # Compute u=M(λ)v and w=M'(λ)v
            u[:] = compute_Mlincomb(nep,λ,v,[T(1)],0);
            w[:] = compute_Mlincomb(nep,λ,v,[T(1)],1);

            # Intermediate quantities
            Δλ=-dot(ws,u)/dot(ws,w);
            z=Δλ*w+u;
            # Throws an error if lin_solve returns incorrect type.
            local Δv::Vector{T}= -lin_solve(linsolver, z, tol=tol)

            normΔv=norm(Δv);
            push_info!(logger,2," norm(Δv)=$normΔv",continues=true)

            (Δλ,Δv,j,scaling)=armijo_rule(nep,errmeasure,err,
                                          λ,v,Δλ,Δv,real(T(armijo_factor)),armijo_max)

            if (j>0)
                push_info!(logger," Armijo scaling=$scaling")
            else
                push_info!(logger,"");
            end

            # Update eigenpair
            λ += Δλ
            v[:] += Δv; # eigvec update

        end

        msg="Number of iterations exceeded. maxit=$(maxit)."
        throw(NoConvergenceException(λ,v,err,msg))
    end

"""
    λ,v = newtonqr([eltype],nep::NEP;[errmeasure,][tol,][maxit,][λ,][v,][c,][logger])

This function implements the Newton-QR method as formulated in the reference. The method ivolves the computation of a rank-revealing QR factorization
of ``M(λ)``, with the idea that on convergence the the last diagonal element ``R[n,n]`` of the upper-triangular matrix ``R`` becomes zero as a result of ``M(λ)``
becoming singular. Since the computation of a QR factorization is expensive, it is advisable to use this method for problems of small size or problems with
a certain structure that makes the QR computation less expensive.
See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> nep=nep_gallery("pep0")
julia> λ,v=newtonqr(nep,v=ones(size(nep,1)));
julia> norm(compute_Mlincomb(nep,λ,v))/norm(v)
8.440206093655014e-15
```

# References
* Kublanovskaya, V. N., (1970).  On an approach to the solution of the generalized latent value problem for λ-matrices, SIAM J. Numer. Anal. 7, 532–537
* Güttel, S., & Tisseur, F. (2017). The nonlinear eigenvalue problem. Acta Numerica, 26, 1-94. doi:10.1017/S0962492917000034
"""
    newtonqr(nep::NEP;params...)=newtonqr(ComplexF64,nep;params...)
    function newtonqr(::Type{T},
                      nep::NEP;
                      errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
                      tol::Real=eps(real(T))*100,
                      maxit::Int=100,
                      λ::Number=zero(T),
                      v::Vector=randn(real(T),size(nep,1)),
                      c::Vector=v,
                      logger=0) where T

        @parse_logger_param!(logger)

        # Ensure types λ and v are of type T
        λ=T(λ)
        v=Vector{T}(v)
        c=Vector{T}(c)

        n = size(nep,1);
        local err
        local w

        en = zeros(n);
        en[n] = 1;

        for k=1:maxit
            A = compute_Mder(nep,λ);
            Q,R,PI = qr(A, Val(true)) #QR factorization with pivoting.
            Q = Matrix(Q)

            P = Matrix{T}(I, n, n)[:,PI] #The permutation matrix corresponding to the pivoted QR.

            p = R[1:n-1,1:n-1]\R[1:n-1,n];
            v = P*[-p;T(1)];#Right eigenvector
            w = Q*en;#Left eigenvector

            #err = abs(R[n,n])/norm(compute_Mder(nep,λ),2); # Frobenius norm
            err=estimate_error(errmeasure,λ,v);

            push_iteration_info!(logger,k,err=err,λ=λ,v=v);
            if(err < tol)
                return λ,v,w;
            end

            d = dot(Q[:,n],compute_Mlincomb(nep,λ,reshape(v,n,1),[T(1)],1));
            #d = dot(Q[:,n],compute_Mder(nep,λ,1)*P*[-p;T(1.0)]);
            λ = λ - R[n,n]/d;
        end

        msg="Number of iterations exceeded. maxit=$(maxit)."
        throw(NoConvergenceException(λ,v,err,msg))
    end

"""
    λ,v = implicitdet([eltype],nep::NEP;[errmeasure,][tol,][maxit,][λ,][v,][c,][logger])

This function implements the Implicit determinant method as formulated Algorithm 4.3 in the reference. The method applies Newton-Raphson to the equation
``det(M(λ))/det(G(λ)) = 0``, where ``G(λ)`` is a saddle point matrix with ``M(λ)``
in the (1,1) block. The (2,1) and (1,2) blocks of ``G(λ)`` are set to
``c^H`` and ``c`` respectively. Note that ``G(λ) `` can be non-singular even when ``M(λ) ``
is singular. See reference for more information.
See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> nep=nep_gallery("pep0")
julia> λ,v=implicitdet(nep,v=ones(size(nep,1)));
julia> norm(compute_Mlincomb(nep,λ,v))/norm(v)
2.566371972986362e-14
```

# References
* Spence, A., & Poulton, C. (2005). Photonic band structure calculations using nonlinear eigenvalue techniques, J. Comput. Phys., 204 (2005), pp. 65–8
* Güttel, S., & Tisseur, F. (2017). The nonlinear eigenvalue problem. Acta Numerica, 26, 1-94. doi:10.1017/S0962492917000034
"""
    implicitdet(nep::NEP;params...)=implicitdet(ComplexF64,nep;params...)
    function implicitdet(::Type{T},
                         nep::NEP;
                         errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
                         tol=eps(real(T))*100,
                         maxit=100,
                         λ=zero(T),
                         v=randn(real(T),size(nep,1)),
                         c=v,
                         logger=0) where T

        @parse_logger_param!(logger)

        n = size(nep,1);
        v = Vector{T}(vcat(v,one(T)))
        vp = zeros(T,n+1)
        c = Vector{T}(c);
        b = c;
        P = Matrix{T}(I, n+1, n+1)

        local err

        for k=1:maxit

            A = compute_Mder(nep,λ);
            AA = [A b;c' 0];#The matrix G(λ)

            F = lu(AA);

            v[:] = F\([zeros(T,n);T(1)]);
            vp[:] = F\([-1*compute_Mder(nep,λ,1)*v[1:n];0]);

            #err = estimate_error(errmeasure,λ,v[1:n]);
            err = abs(v[n+1])/norm(compute_Mder(nep,λ),2); # Frobenius norm based error
            push_iteration_info!(logger,k,err=err,λ=λ,v=v);
            if(err < tol)
                return λ,v[1:n];
            end

            λ = λ - v[n+1]/vp[n+1];#Newton update for the equation det(M(λ))/det(G(λ)) = 0
        end

        msg="Number of iterations exceeded. maxit=$(maxit)."
        throw(NoConvergenceException(λ,v,NaN,msg))
    end

    # Armijo rule implementation
    function armijo_rule(nep,errmeasure,err0,λ,v,Δλ,Δv,armijo_factor,armijo_max)
        j=0
        if (armijo_factor<1)
            # take smaller and smaller steps until errmeasure is decreasing
            while (estimate_error(errmeasure,λ+Δλ,v+Δv)>err0 && j<armijo_max)
                j=j+1;
                Δv=Δv*armijo_factor;
                Δλ=Δλ*armijo_factor;
            end
        end
        return  (Δλ,Δv,j,armijo_factor^j)
    end