text:
export Result,
    ResultStatistics,
    statistics,
    seed,
    path_results,
    results,
    nresults,
    solutions,
    real_solutions,
    nonsingular,
    singular,
    at_infinity,
    failed,
    nfinite,
    nsolutions,
    nsingular,
    nat_infinity,
    nexcess_solutions,
    nfailed,
    nnonsingular,
    nreal,
    ntracked

######################
## MultiplicityInfo ##
######################

"""
    MultiplicityInfo

This contains information about the multiplicities of the solutions.
"""
struct MultiplicityInfo
    # indices of multiple solutions grouped by multiplicity
    multiplicities::Dict{Int,Vector{Vector{Int}}}
    # This set stores all path numbers which we want to ignore
    # if we don't show multiple results
    multiple_indicator::Set{Int}
end

function MultiplicityInfo(pathresults::Vector{PathResult})
    multiple_indicator = Set{Int32}()
    multiplicities = compute_multiplicities(pathresults)
    for clusters in values(multiplicities), cluster in clusters
        for i = 2:length(cluster)
            push!(multiple_indicator, path_number(pathresults[cluster[i]]))
        end
    end
    # the multiplicities are currently with respect to the indices of the pathresults
    # but we want to have them with respect to path numbers in case pathresults
    # is just a subset of all paths
    for k in keys(multiplicities)
        multiplicities[k] = map(multiplicities[k]) do cluster
            map(i -> path_number(pathresults[i]), cluster)
        end
    end
    MultiplicityInfo(multiplicities, multiple_indicator)
end

function compute_multiplicities(result::Vector{<:PathResult}; kwargs...)
    D = Dict{Int,Vector{Vector{Int}}}()
    for m in multiplicities(solution, result)
        if haskey(D, length(m))
            push!(D[length(m)], m)
        else
            D[length(m)] = [m]
        end
    end
    D
end

function assign_multiplicities!(path_results::Vector{<:PathResult}, I::MultiplicityInfo)
    #assign multiplicities
    for (k, clusters) in I.multiplicities, cluster in clusters, vᵢ in cluster
        path_results[vᵢ].multiplicity =
            max(k, something(path_results[vᵢ].winding_number, 1))
    end
    path_results
end

"""
    Result

The result of [`solve`](@ref). This is a wrapper around the results of each single path
([`PathResult`](@ref)) and it contains some additional information like a random seed to
replicate the result.
"""
struct Result
    path_results::Vector{PathResult}
    tracked_paths::Int
    seed::Union{Nothing,UInt32}
    start_system::Union{Nothing,Symbol}
    multiplicity_info::MultiplicityInfo
end
function Result(
    path_results::Vector{PathResult};
    seed::Union{Nothing,UInt32} = nothing,
    tracked_paths = length(path_results),
    start_system = nothing,
)
    multiplicity_info = MultiplicityInfo(filter(is_singular, path_results))
    assign_multiplicities!(path_results, multiplicity_info)
    Result(path_results, tracked_paths, seed, start_system, multiplicity_info)
end

Base.size(r::Result) = (length(r),)
Base.length(r::Result) = length(r.path_results)
Base.getindex(r::Result, I) = getindex(r.path_results, I)

Base.iterate(r::Result) = iterate(r.path_results)
Base.iterate(r::Result, state) = iterate(r.path_results, state)
Base.lastindex(r::Result) = lastindex(r.path_results)
Base.eltype(r::Type{Result}) = PathResult
Base.keys(r::Result) = 1:length(r.path_results)
Base.values(r::Result) = r.path_results

"""
    seed(::Result)

Returns the seed to replicate the result.
"""
seed(R::Result) = R.seed

"""
    path_results(::Result)

Returns the stored [`PathResult`](@ref)s.
"""
path_results(R::Result) = R.path_results

is_multiple_result(r::PathResult, R::Result) =
    path_number(r) ∈ R.multiplicity_info.multiple_indicator
is_multiple_result(r::PathResult, R::AbstractVector{PathResult}) = false

Base.@kwdef struct ResultStatistics
    total::Int
    nonsingular::Int
    singular::Int
    singular_with_multiplicity::Int = singular
    real::Int
    real_nonsingular::Int
    real_singular::Int
    real_singular_with_multiplicity::Int = real_singular
    at_infinity::Int
    excess_solution::Int = 0
    failed::Int = 0
end
Base.show(io::IO, stats::ResultStatistics) = print_fieldnames(io, stats)

function ResultStatistics(result::Result; real_tol::Float64 = 1e-6)
    failed = at_infinity = excess_solution = 0
    nonsingular = singular = real_nonsingular = real_singular = 0
    singular_with_multiplicity = real_singular_with_multiplicity = 0
    for r in result
        is_multiple_result(r, result) && continue
        if is_failed(r)
            failed += 1
        elseif is_at_infinity(r)
            at_infinity += 1
        elseif is_excess_solution(r)
            excess_solution += 1
        elseif is_singular(r)
            if is_real(r, real_tol)
                real_singular += 1
                real_singular_with_multiplicity += something(multiplicity(r), 1)
            end
            singular += 1
            singular_with_multiplicity += something(multiplicity(r), 1)
        else # finite, nonsingular
            if is_real(r, real_tol)
                real_nonsingular += 1
            end
            nonsingular += 1
        end
    end
    ResultStatistics(
        nonsingular = nonsingular,
        singular = singular,
        singular_with_multiplicity = singular_with_multiplicity,
        real_nonsingular = real_nonsingular,
        real_singular = real_singular,
        real_singular_with_multiplicity = real_singular_with_multiplicity,
        real = real_nonsingular + real_singular,
        at_infinity = at_infinity,
        excess_solution = excess_solution,
        failed = failed,
        total = result.tracked_paths,
    )
end

"""
    statistics(R::Result; real_tol = 1e-6)

Statistic about the number of (real) singular and non-singular solutions etc.
"""
statistics(r; kwargs...) = ResultStatistics(r; kwargs...)

const Results = Union{Result,AbstractVector{<:PathResult}}

"""
    results(
        result;
        only_real = false,
        real_tol = 1e-6,
        only_nonsingular = false,
        only_singular = false,
        only_finite = true,
        multiple_results = false,
    )
    results(f, result; options...)

Return all [`PathResult`](@ref)s for which satisfy the given conditions and apply,
if provided, the function `f`.
"""
results(R::Results; kwargs...) = results(identity, R; kwargs...)
function results(
    f::Function,
    R::Results;
    only_real::Bool = false,
    real_tol::Float64 = 1e-6,
    only_nonsingular::Bool = false,
    only_singular::Bool = false,
    only_finite::Bool = true,
    multiple_results::Bool = false,
)
    [
        f(r) for r in R if (!only_real || is_real(r, real_tol)) &&
            (!only_nonsingular || is_nonsingular(r)) &&
            (!only_singular || is_singular(r)) &&
            (!only_finite || is_finite(r)) &&
            (multiple_results || !is_multiple_result(r, R))
    ]
end

"""
    nresults(
        result;
        only_real = false,
        real_tol = 1e-6,
        only_nonsingular = false,
        only_singular = false,
        only_finite = true,
        multiple_results = false,
    )

Count the number of results which satisfy the corresponding conditions. See also
[`results`](@ref).
"""
function nresults(
    R::Results;
    only_real::Bool = false,
    real_tol::Float64 = 1e-6,
    only_nonsingular::Bool = false,
    only_singular::Bool = false,
    onlyfinite::Bool = true, # deprecated
    only_finite::Bool = onlyfinite,
    multiple_results::Bool = false,
)
    count(R) do r
        (!only_real || is_real(r, real_tol)) &&
            (!only_nonsingular || is_nonsingular(r)) &&
            (!only_singular || is_singular(r)) &&
            (!only_finite || isfinite(r)) &&
            (multiple_results || !is_multiple_result(r, R))
    end
end

"""
    solutions(result; only_nonsingular = true, conditions...)

Returns all solutions for which the given conditions apply, see [`results`](@ref) for the
possible conditions.

## Example
```julia-repl
julia> @var x y
julia> F = System([(x-2)y, y+x+3]);
julia> solutions(solve(F))
2-element Array{Array{Complex{Float64},1},1}:
 [2.0 + 0.0im, -5.0 + 0.0im]
 [-3.0 + 0.0im, 0.0 + 0.0im]
```
"""
solutions(result::Results; only_nonsingular = true, kwargs...) =
    results(solution, result; only_nonsingular = only_nonsingular, kwargs...)

"""
    real_solutions(result; tol=1e-6, conditions...)

Return all real solution for which the given conditions apply.
For the possible `conditions` see [`results`](@ref).
Note that `only_real` is always `true` and `real_tol` is now `tol`.

## Example
```julia-repl
julia> @var x y;
julia> F = System([(x-2)y, y+x+3]);
julia> real_solutions(solve(F))
2-element Array{Array{Float64,1},1}:
 [2.0, -5.0]
 [-3.0, 0.0]
```
"""
function real_solutions(result::Results; tol::Float64 = 1e-6, kwargs...)
    results(real ∘ solution, result; only_real = true, real_tol = tol, kwargs...)
end

"""
    nonsingular(result; conditions...)

Return all [`PathResult`](@ref)s for which the solution is non-singular.
This is just a shorthand for `results(R; only_nonsingular=true, conditions...)`.
For the possible `conditions` see [`results`](@ref).
"""
nonsingular(R::Results; kwargs...) = results(R; only_nonsingular = true, kwargs...)

"""
    singular(result; multiple_results=false, kwargs...)

Return all [`PathResult`]s for which the solution is singular.
If `multiple_results=false` only one point from each cluster of multiple solutions is
returned. If `multiple_results = true` all singular solutions in `R` are returned.
For the possible `kwargs` see [`results`](@ref).
"""
function singular(R::Results; kwargs...)
    results(R; only_singular = true, kwargs...)
end

"""
    real(result, tol=1e-6)

Get all results where the solutions are real with the given tolerance `tol`.
See [`is_real`](@ref) for details regarding the determination of 'realness'.
"""
Base.real(R::Results; tol::Float64 = 1e-6) = filter(r -> is_real(r, tol), path_results(R))

"""
    failed(result)
Get all results where the path tracking failed.
"""
failed(R::Results) = filter(is_failed, path_results(R))

"""
    at_infinity(result)

Get all results where the solutions is at infinity.
"""
at_infinity(R::Results) = filter(is_at_infinity, path_results(R))

"""
    nsolutions(result; only_nonsingular = true, options...)

The number of solutions. See [`results`](@ref) for the possible options.
"""
nsolutions(R::Results; only_nonsingular = true, options...) =
    nresults(R; only_nonsingular = only_nonsingular, options...)

"""
    nsingular(
        result;
        counting_multiplicities = false,
        kwargs...,
    )

The number of singular solutions. A solution is considered singular if its winding number is
larger than 1 or the condition number is larger than `tol`.
If `counting_multiplicities=true` the number of singular solutions times their
multiplicities is returned.
"""
function nsingular(R::Results; counting_multiplicities::Bool = false, kwargs...)
    S = results(R; only_singular = true, multiple_results = false, kwargs...)
    isempty(S) && return 0
    counting_multiplicities && return sum(multiplicity, S)
    length(S)
end

"""
    nat_infinity(result)

The number of solutions at infinity.
"""
nat_infinity(R::Results) = count(is_at_infinity, R)

"""
    nexcess_solutions(result)

The number of exess solutions. See also [`excess_solution_check`](@ref).
"""
nexcess_solutions(R::Results) = count(is_excess_solution, R)

"""
    nfailed(result)

The number of failed paths.
"""
nfailed(R::Results) = count(is_failed, R)

"""
    nnonsingular(result)

The number of non-singular solutions. See also [`is_singular`](@ref).
"""
nnonsingular(R::Result) = count(is_nonsingular, R)

"""
    nreal(result; tol=1e-6)

The number of real solutions. See also [`is_real`](@ref).
"""
nreal(R::Results; tol = 1e-6) = nresults(R, only_real = true, real_tol = tol)

"""
    ntracked(result)

Returns the total number of paths tracked.
"""
ntracked(R::Result) = R.tracked_paths

###
### Show
####
function Base.show(io::IO, x::Result)
    s = statistics(x)
    total = s.nonsingular + s.singular
    header = "Result with $total $(plural("solution", total))"
    println(io, header)
    println(io, "="^(length(header)))
    println(io, "• $(ntracked(x)) $(plural("path", ntracked(x))) tracked")
    println(
        io,
        "• $(s.nonsingular) non-singular $(plural("solution", s.nonsingular)) ",
        "($(s.real_nonsingular) real)",
    )
    if s.singular > 0
        println(
            io,
            "• $(s.singular) singular $(plural("solution", s.singular)) ",
            "($(s.real_singular) real)",
        )
    end
    # s.at_infinity > 0 &&
    #     println(io, "• $(s.at_infinity) $(plural("solution", s.at_infinity)) at infinity")
    s.excess_solution > 0 && println(
        io,
        "• $(s.excess_solution) excess $(plural("solution", s.excess_solution))",
    )
    println(io, "• random_seed: ", sprint(show, seed(x)))
    if !isnothing(x.start_system)
        println(io, "• start_system: :", x.start_system)
    end
    if s.singular > 0
        println(io, "• multiplicity table of singular solutions:")
        singular_multiplicities_table(io, x, s)
    end
end

function TreeViews.treelabel(io::IO, x::Result, ::MIME"application/prs.juno.inline")
    s = statistics(x)
    total = s.nonsingular + s.singular
    print(
        io,
        "Result" *
        " with $total $(plural("solution", total))",
    )
end
TreeViews.hastreeview(::Result) = true
TreeViews.numberofnodes(::Result) = 10
function TreeViews.nodelabel(io::IO, x::Result, i::Int, ::MIME"application/prs.juno.inline")
    s = statistics(x)
    if i == 1
        print(io, "Paths tracked")
    elseif i == 2 && s.nonsingular > 0
        print(io, "$(s.nonsingular) non-singular ($(s.real_nonsingular) real)")
    elseif i == 3 && s.singular > 0
        print(io, "$(s.singular) singular ($(s.real_singular) real)")
    elseif i == 4 && (s.real_nonsingular + s.real_singular) > 0
        print(io, "$(s.real_nonsingular+s.real_singular) real")
    elseif i == 5 && s.excess_solution > 0
        print(io, "$(s.excess_solution) excess $(plural("solution", s.excess_solution))")
    elseif i == 6 && s.at_infinity > 0
        print(io, "$(s.at_infinity) at_infinity")
    elseif i == 7 && s.failed > 0
        print(io, "$(s.failed) failed")
    elseif i == 8
        print(io, "Random seed used")
    elseif i == 9 && !isnothing(x.start_system)
        print(io, "start_system")
    elseif i == 10 && s.singular > 0
        print(io, "  multiplicity table of singular solutions: \n")
        singular_multiplicities_table(io, x, s)
    end
end

function TreeViews.treenode(r::Result, i::Integer)
    s = statistics(r)
    if i == 1
        return ntracked(r)
    elseif i == 2 && s.nonsingular > 0
        return results(r, only_nonsingular = true)
    elseif i == 3 && s.singular > 0
        return results(r, only_singular = true)
    elseif i == 4 && (s.real_nonsingular + s.real_singular) > 0
        return results(r, only_real = true)
    elseif i == 5 && s.excess_solution > 0
        return filter(is_excess_solution, path_results(r))
    elseif i == 6 && s.at_infinity > 0
        return at_infinity(r)
    elseif i == 7 && s.failed > 0
        return failed(r)
    elseif i == 8
        return seed(r)
    elseif i == 9 && !isnothing(r.start_system)
        return r.start_system
    elseif i == 10 && s.singular > 0
        return missing
    end
    missing
end

function singular_multiplicities_table(io, result::Result, stats = statistics(result))
    M = result.multiplicity_info.multiplicities
    if isempty(M)
        n_higher_mults_total = 0
    else
        n_higher_mults_total = sum(length, values(M))
    end

    headers = ["mult.", "total", "# real", "# non-real"]
    mult_1_exists = n_higher_mults_total < stats.singular
    data = Matrix{Int}(undef, mult_1_exists + length(M), 4)

    n_higher_real = 0
    i = mult_1_exists + 1
    for k in sort(collect(keys(M)))
        data[i, 1] = k
        n_real_solsᵢ = count(Mᵢ -> is_real(result.path_results[Mᵢ[1]]), M[k])
        n_solsᵢ = length(M[k])
        data[i, 2] = n_solsᵢ
        data[i, 3] = n_real_solsᵢ
        data[i, 4] = n_solsᵢ - n_real_solsᵢ
        n_higher_real += n_real_solsᵢ
        i += 1
    end

    if mult_1_exists
        data[1, 1] = 1
        data[1, 2] = stats.singular - n_higher_mults_total
        data[1, 3] = stats.real_singular - n_higher_real
        data[1, 4] = data[1, 2] - data[1, 3]
    end

    PrettyTables.pretty_table(
        io,
        data,
        headers;
        tf = PrettyTables.tf_unicode_rounded,
        alignment = :c,
        header_crayon = PrettyTables.Crayon(bold = false),
        border_crayon = PrettyTables.Crayon(faint = true),
    )
end


text:
using Core: Argument
#
# Provides `create_array` methods for `StateVector` and `ParameterVector`
# instances!
#

@inline function SymbolicUtils.Code.create_array(A::Type{<:StateVector}, T, nd::Val, d::Val{dims}, elems...) where {dims}
    data = SymbolicUtils.Code.create_array(SArray, T, nd, d, elems...)
    return (A)(data.data)
end

ext: jl


text:
## --- 1D Histograms
    n = 10000
    a = rand(n)
    h = histcounts(a, 0:0.1:1)
    @test h == histcounts(a, 0, 1, 10)
    @test isa(h, Array{Int64,1})
    @test length(h) == 10
    @test sum(h) == n

    @test histcounts(1:100,0:10:100) == fill(10,10)
    @test histcounts(1:100,0,100,10) == fill(10,10)
    @test histcounts(1:100.,0.,100.,10) == fill(10,10)

    # Test results and warnings when N is too small to hold results
    w = "length(N) < nbins; any bins beyond length(N) will not be filled"
    @test (@test_logs (:warn, w) histcounts!(zeros(5), 1:100 ,0:10:100)) == fill(10,5)
    @test (@test_logs (:warn, w) histcounts!(zeros(5), 1:100. ,0:10:100)) == fill(10,5)

## --- 2D histograms

    x = y = 0.5:9.5
    xedges = yedges = 0:10
    N = histcounts(x,y,xedges,yedges)
    @test N == I(10)

    # Test results and warnings when N is too small to hold results
    w = "size(N,1) < nybins; any y bins beyond size(N,1) will not be filled"
    @test (@test_logs (:warn, w) histcounts!(zeros(5,10), x, y, xedges, yedges)) == I(10)[1:5,:]
    w = "size(N,2) < nxbins; any x bins beyond size(N,2) will not be filled"
    @test (@test_logs (:warn, w) histcounts!(zeros(10,5), x, y, xedges, yedges)) == I(10)[:,1:5]

## --- End of File


text:
miss_vec = [" ","0",""]
function quality_control(file,miss_vec; delim = ',', header = 0)
    """
    Zigui Section
    Remove duplicates
    Example: a8  a4  a6
             a8  a4  a6
    Return:  a8  a4  a6

    Find row that misses son
    Example:    a4  a6

    Find son that misses both sire and dam
    Example: a5
    """

    dat = CSV.read("pedigree_clean.txt", delim = ",", header = 1,missingstrings = miss_vec)
    #Potential risk, now read nothing will return "" instead of missing

    dat = unique!(dat) #remove duplicate
    miss_son_row=ismissing.(dat[:,1])  #find the row that miss son
    miss_par_row = (ismissing.(dat[:,2]).|ismissing.(dat[:,3])) #Find the row that miss both sire and dam
    miss_son_index = findall(miss_son_row)
    miss_par_index = findall(miss_par_row)
    miss_par_son = dat[miss_par,1]
    if sum(miss_son_row) >0 @warn "Find missing son record at row $miss_son_index" end
    if sum(miss_par_row) >0 @warn "Find son $miss_par_son that miss parents at row $miss_par_index" end

    println("Zigui section finish")

    dropmissing!(dat) #Remove all rows containing missing value (Do it for following code can be run)

    """
    Jiayi Section
    Remove individuals that are both parents and children
    Input:
    a5	a9	a7
    a6	a5	a4
    a7	a7	a8
    a9	a5	a4

    Output:
    a5	a9	a7
    a6	a5	a4
    a9	a5	a4

    """
    index = collect(1:1:nrow(dat))

    ConflictRole = index[dat[1] .== dat[2]]
    append!(ConflictRole, index[dat[1] .== dat[3]])
    Ind = String.(dat[ConflictRole, 1])
    deleterows!(dat, ConflictRole)

    if !isempty(ConflictRole)
        @warn "Conflict Role observed for individual $Ind"
    end

    println("Jiayi section finish")

    """
    Tianjing Section
    Goal:    warn duplicate son (this step is after removing all duplicate)
    Example: a4 a1 a2
             a4 a9 a10
             a4 a6 a2
    """
    root_signal = "0"
    dup_bool = nonunique(dat, 1)
    dup_son = unique(dat[1][dup_bool])
    #warning
    index_dup_son = findall(x -> x in dup_son, dat[1])
    dup_data = dat[index,:]
    @warn "Same son with different parents in row $index_dup_son
    $dup_data"

    """
    Goal:    warn element in both dam and sire
    Example: a8  a4 a6
             a9  a4 a6
             a13 a6 a2
    """
    dup = intersect(dat[2], dat[3])
    filter!(e -> e ≠ root_signal, dup)  #remove the dam and sire of root, eg. "0"
    dup_both = vcat(findall(x->x in dup, dat[2]),findall(x->x in dup, dat[3])) #find index of dup
    #warning
    unique!(dup_both)
    sort!(dup_both)
    dup_dat= dat[dup_both,:]
    @warn "Parents in row $dup_both appear in both dam and sire
    $dup_dat"

    println("Tianjing section finish")
    return dat

end

#Testing
quality_control("pedigree_clean.txt",miss_vec)


text:
# AddIDs

const SEL_H1234 = SelectOr(SelectTag.((:h1, :h2, :h3, :h4)))
const SEL_H234 = SelectOr(SelectTag.((:h2, :h3, :h4)))

"""
    AddID(selector = sel"h2, h3, h4"; idfn)

`Replacer` that adds an `id` attribute to every x-expression
selected by `selector`. It assumes that its first child is a `String`.
The id is created by applying `idfn` to that string. `idfn` defaults
to `CommonMark.slugify`.
"""
AddID(sel=SEL_H234; idfn=CommonMark.slugify) =
    Replacer(x -> addid(x; idfn=idfn), sel)

function addid(x; idfn=CommonMark.slugify)
    #(!isempty(children(x)) && children(x)[1] isa String) || error(
    #    "To add an ID, first child must be a `String`. Got xexpr\n$doc")
    text = gettext(x)
    #text = children(x)[1]
    id = idfn(text)
    return withattributes(x, Dict(:id => id))
end

# HTMLify

function HTMLify(sel=SelectNode(), htmltag=:div)
    return Replacer(sel) do x
        htmlify(x, htmltag)
    end
end

"""
    htmlify(doc, htmltag = :div)

If `doc.tag` is not a valid HTML tag, changes it into a :div and adds the attribute
`:class => doc.tag`.
"""
function htmlify(doc, htmltag=:div)
    if tag(doc) in HTMLTAGS
        return doc
    else
        return XNode(
            htmltag,
            # TODO: add to classes if exist
            merge(attributes(doc), Dict(:class => string(tag(doc)))),
            children(doc),
        )
    end
end

# ChangeLinkExtensions

function ChangeLinkExtension(ext, sel::Selector = SelectTag(:a); linkattr = :href)
    return Replacer(x -> changelinkextension(x, ext; attr = linkattr), sel)
end

function changelinkextension(doc::XNode, ext; attr = :href)
    if haskey(attributes(doc), attr)
        href = doc.attributes[attr]
        if startswith(href, "http") || startswith(href, "www")
            return doc
        else
            return withattributes(
                doc,
                merge(doc.attributes, Dict(attr => changehrefextension(href, ext))),
            )
        end
    else
        return doc
    end
end

# ChangeTag

function ChangeTag(t, sel)
    return Replacer(x -> withtag(x, t), sel)
end

# FormatCode

function FormatCode(codesel = SelectTag(:pre))
    return Replacer(x -> formatcodeblock(x), codesel)
end

function formatcodeblock(doc)
    if get(attributes(doc), :lang, "") == "julia"
        code = gettext(doc)
        try
            code = format_text(code)
        catch
            @info "Parsing error when formatting code snippet: \n\n$code"
        end
        return withchildren(doc, [XNode(:code, XLeaf(code))])
    else
        return doc
    end
end

#

const CSSLINKSELECTOR = SelectTag(:link) & SelectHasAttr(:href)

Base.@kwdef struct RelativeLinks <: Rewriter
    linktag::Symbol = :a
    linkattr::Symbol = :href
end

function rewritedoc(rewriter::RelativeLinks, p, doc)
    sel = SelectTag(rewriter.linktag) & SelectHasAttr(rewriter.linkattr)
    cata(doc, sel) do x
        href = attributes(x)[rewriter.linkattr]
        if startswith(href, '/')
            newhref = relpath(href, "/" * string(parent(p)))
            return withattributes(x, merge(attributes(x), Dict(rewriter.linkattr => newhref)))
        else
            return x
        end
    end
end

#

function createtitle(p, x)
    h1 = selectfirst(x, SelectTag(:h1))
    title = if !isnothing(h1)
        gettext(h1)
    else
        string(filename(p))
    end
    return XNode(:title, [XLeaf(title)])
end

struct HTMLRedirect <: Rewriter
    p::AbstractPath
end

function postbuild(redirect::HTMLRedirect, project, builder)
    builder isa FileBuilder || error("`HTMLRedirect` does not work with $builder")
    redirectpath = withext(redirect.p, formatextension(builder.format))
    content = """
    <!DOCTYPE html>
    <html>
    <head>
    <meta http-equiv = "refresh" content = "0; url = $redirectpath" />
    </head>
    </html>
    """
    open(joinpath(builder.dir, "index.html"), "w") do f
        write(f, content)
    end
end


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

mutable struct UpdateVirtualNetworkRuleProperties <: SwaggerModel
    subnetId::Any # spec type: Union{ Nothing, String } # spec name: subnetId

    function UpdateVirtualNetworkRuleProperties(;subnetId=nothing)
        o = new()
        validate_property(UpdateVirtualNetworkRuleProperties, Symbol("subnetId"), subnetId)
        setfield!(o, Symbol("subnetId"), subnetId)
        o
    end
end # type UpdateVirtualNetworkRuleProperties

const _property_map_UpdateVirtualNetworkRuleProperties = Dict{Symbol,Symbol}(Symbol("subnetId")=>Symbol("subnetId"))
const _property_types_UpdateVirtualNetworkRuleProperties = Dict{Symbol,String}(Symbol("subnetId")=>"String")
Base.propertynames(::Type{ UpdateVirtualNetworkRuleProperties }) = collect(keys(_property_map_UpdateVirtualNetworkRuleProperties))
Swagger.property_type(::Type{ UpdateVirtualNetworkRuleProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_UpdateVirtualNetworkRuleProperties[name]))}
Swagger.field_name(::Type{ UpdateVirtualNetworkRuleProperties }, property_name::Symbol) =  _property_map_UpdateVirtualNetworkRuleProperties[property_name]

function check_required(o::UpdateVirtualNetworkRuleProperties)
    true
end

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


text:
module Absolute_Loss
    using LoopVectorization

    function func(output_matrix::Array{Float32}, sample_matrix::Array{Float32})
        loss_matrix = zeros(Float32, size(output_matrix))
        @avxt for i in axes(output_matrix, 1), j in axes(output_matrix, 2)
            loss_matrix[i,j] = abs(output_matrix[i,j]-sample_matrix[i,j])
        end
        return loss_matrix
    end

    function prop!(δ::Array{Float32}, output_matrix::Array{Float32}, sample_matrix::Array{Float32})
        @avxt for i in eachindex(δ)
            δ[i] = ifelse(output_matrix[i]-sample_matrix[i]>=0, 1, -1)
        end
    end
end


text:
#!/usr/local/bin/julia
# export JULIA_LOAD_PATH="/path/to/module/storage/folder"
#include("/Users/stewart/Projects/jlVCF/src/jlVCF.jl")
# ARGS=["/Users/stewart/Downloads/753_pair.MuTect2.call_stats.vcf","/Users/stewart/PycharmProjects/vcf2txt/out/test.tsv"]

include("jlVCF.jl")
include("VCFIterator.jl")
include("Variant.jl")

narg = 0
f=stdout
max_events=2^31
vcf_file1=ARGS[1]

narg=length(ARGS)
if (narg>1)
	a=ARGS[2]
	#if isnull(tryparse(Int32,a))
	if tryparse(Int32,a)==nothing
		println("vcf:\t",ARGS[1])
		println("tsv:\t",a)
		f = open(a,"w")
	else
	  max_events=parse(Int32,a)
	end
end

delim="\t"
c=fill("", 500)
num=Array{AbstractString,500}
if (max_events<1)
	delim="\n"
	c=Array{AbstractString,500}
	for i in 1:500
	  c[i]=string(i,".\t")
	end
end
#test loading
vc = VCFIterator(vcf_file1)

print(f,c[1],"CHRO",delim)
print(f,c[2],"POS",delim)
print(f,c[3],"ID",delim)
print(f,c[4],"REF",delim)
print(f,c[5],"ALT",delim)
print(f,c[6],"QUAL",delim)
print(f,c[8],"FILTER",delim)
global n1=8
for x in keys(vc.infoTypes)
    global n1+=1
	print(f,c[n1],x,delim)
end

NS=length(vc.samples)
for s in vc.samples
	for x in keys(vc.formatTypes)
		if vc.formatTypes[x].Number=="R"
		    global n1+=1
	  		print(f,c[n1],s,"_",x,"_REF",delim)
		    global n1+=1
	  		print(f,c[n1],s,"_",x,"_ALT",delim)
		else
		    global n1+=1
	  		print(f,c[n1],s,"_",x,delim)
		end
	end
end
print(f,"\n" )

count = 0
v = "dmy"
while !eof(vc)
    v = next1(vc)
    if isa(v,Bool)
    	continue
    end
    if sizeof(v)<2
    	continue
    end
    global count += 1
    if count>max_events
	  break
	end
	#print(f,count,"\t")
	print(f,v.CHROM,"\t")
	print(f,v.POS,"\t")
	print(f,join(v.ID,':'),"\t")
	print(f,v.REF,"\t")
	print(f,join(v.ALT,':'),"\t")
	print(f,v.QUAL,"\t")
	print(f,join(v.FILTER,':'),"\t")
	for x in keys(vc.infoTypes)
	  print(f,join((get(v.INFO,x,"")),":"),"\t")
	end
	for s in vc.samples
		for x in keys(vc.formatTypes)
		  if vc.formatTypes[x].Number=="R"
			print(f,v.FORMAT[s][x][1],"\t")
	  		print(f,v.FORMAT[s][x][2],"\t")
	 	  else
	  		q=get(v.FORMAT[s],x,"")
	  		print(f,join(q,":"),"\t")
	  	  end
	  	end
	end
	print(f,"\n" )
end

close(vc)


text:
# Evaluates the ITMIX runs by creating a `data` DataFrame containing summary info and, optionally, plots.
#
# NOTE: the run-plots show the mean thickness field not the one produced by using mtheta in the forward model.

error("remove Rhat")

# run this first to make sure all precompilation this through
include("itmix-setup.jl")
using FileIO, DataFrames, StatsBase
import Plots. PyPlot, Interpolations
global pl

const make_plots = false
const make_plots_thetas = make_plots
parallel = false

n_1d_btilde = 3 # number of spatial parameters
n_1d_fsl = 3 # number of spatial parameters
n_1d_temp = 1 # number of spatial parameters
subdir = ["",
          "res-7oct-2",
          "res-23-oct",
          "res-28-oct",
          "res-30-oct",
          "res-31-oct-rmiv",
          "res-03-nov-rmiv",
          "esa-14nov"
          ][end]

if !isdefined(:runit)
    runit = true
end

"""Construct an empty DataFrame"""
function make_data()
    error("remove Rhat")
    DataFrame(person = Symbol[],
              gl = Symbol[],
              prio = Int[],
              expnr = Int[],
              rm_iv = Bool[],
              nl = Int[],
              mean_h = Float64[],
              volume_area = Float64[],
              rmse = Float64[],
              rmse_fitting = Float64[],
              max = Float64[],
              max_fitting = Float64[],
              Rhat = Float64[],
              Rhat_max = Float64[],
              mtheta = Vector{Float64}[],
              err_h = VAWTools.Traj{Float64}[]
              )
end

function load_one(gid, expnr, rm_iv;
                  fit_target=BM.FitTarget.h_iv)
    gl,gb,pp,pm,pn,pl = gls[gid]
    priority = find([gid in pr for pr in I2.priorities])[1]
    fl, fls, fliv, flivs, flj = I2.filepaths(gl, pl, expnr, rm_iv, "results_Werder/"*subdir)
    dir_,fl_ = splitdir(fl)
    fl_ = splitext(fl_)[1]
    if !isfile(fl)
        println("File $fl not found")
        return nothing
    end

    jldloaded, jldd = try
        true, FileIO.load(flj)
    catch e
        println("File $flj not loading with $e")
        false, nothing
    end

    if jldloaded
        thetas, Rhat, sample_size = map(k->jldd[k], ("thetas","Rhat","sample_size"))
        mtheta = squeeze(mean(thetas,2),2)

        gln = BM.ITMIX2.make_ITMIX2_glacier(gl, expnr)
        gbn = BM.Bands(gb, gl=gln)
        th0d = BM.theta0_dict_defaults(pm, n_1d_btilde, n_1d_fsl, n_1d_temp)
        if length(mtheta)==8
            # make it compatible with older runs
            pop!(th0d, :iv_dist_exp2)
        end
        (theta0, logposterior, logposterior1d, logprior, logprior_debug,
         loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn,
         pmcmc_defaults, fit_target) =
             BM.init_inverse(gbn, pp, pm, pn,
                             theta0_dict=th0d,
                             fit_target=fit_target)
        fwdsol = fwdm_fn(theta0.th0)
    end

    r_lines_fitting = I2.get_rlines(gid, expnr)
    hs2d = VAWTools.read_agr(fl)
    hs2ds = VAWTools.read_agr(fls)
    iv2d = try VAWTools.read_agr(fliv) end
    iv2ds = try VAWTools.read_agr(flivs) end
    if gid==BITEModel.ITMIXGlacier(:Columbia, 2) || gid==BITEModel.ITMIXGlacier(:Elbrus, 2)
        # downsample agr
        hs2d = downsample(hs2d, 2)
        hs2ds = downsample(hs2ds, 2)
    end

    # tuples (control, fitting, all)
    RMSE, max_err, quantiles, smaller_20, smaller_std, smaller_2std, err_h = BM.compare_to_radar(gl, hs2d, hs2ds, r_lines_fitting)

    mean_ice_thickness, volume_area = VAWTools.meannan(hs2d.v), BM.volume_area_mean_h(gl)

    data_out = if jldloaded
        [:Werder, gid.id, priority, expnr, rm_iv, length(r_lines_fitting),
         map(x->signif(x,3), (mean_ice_thickness, volume_area, RMSE[1], RMSE[2], (max_err[1]), (max_err[2]), (mean(Rhat)), (maximum(Rhat))))...,
         mtheta, err_h]
    else
        [:Werder, gid.id, priority, expnr, rm_iv, length(r_lines_fitting),
         map(x->signif(x,3), (mean_ice_thickness, volume_area, RMSE[1], RMSE[2], (max_err[1]), (max_err[2]), 0.0, 0.0))...,
         Float64[], err_h]
    end

    fwdsol = fwdm_fn(mtheta);

    return data_out, hs2d, hs2ds, iv2d, iv2ds, r_lines_fitting, dir_, fl_, fwdsol, jldd, mtheta, thetas, gbn
end

# """
# Evaluate one ITMIX2 experiment and store/save relevant input.

# Returns the data to be appended to data-DataFrame
# """
@everywhere function evaluate_one(expnr, rm_iv, gb, pp, pm, pn, pl;
                                  n_1d_btilde=3, n_1d_fsl=3, n_1d_temp=1,
                                  make_plots = false,
                                  make_plots_thetas = make_plots,
                                  subdir="",
                                  fit_target=BM.FitTarget.h_iv )
    gl = gb.gl
    gid = gl.gid
    (data_out, hs2d, hs2ds, iv2d, iv2ds, r_lines_fitting, dir_,
     fl_, fwdsol, jldd, mtheta, thetas, gbn) = load_one(gid, expnr, rm_iv, fit_target=fit_target)

    # plot it
    if make_plots
        BM.plot2d_h(gl, hs2d.v; clims=:radar, plotline=I2.get_rlines(gid, expnr), reuse=true,
                    outlines=true)
        Plots.savefig(joinpath(dir_, "plots", "$(fl_)-hs2d.png"))
        BM.plot2d_h(gl, hs2ds.v; clims=:model, plotline=I2.get_rlines(gid, expnr), plottraj=false,
                    colorbar_title="h std (m)", reuse=true, outlines=true)
        Plots.savefig(joinpath(dir_, "plots", "$(fl_)-hs2d-std.png"))
        if iv2d!=nothing
            BM.plot2d_h(gl, iv2d.v; plotiv=true, reuse=true, clims=gl.iv.vals isa VAWTools.Traj ? :radar : :model,
                        outlines=true)
            Plots.savefig(joinpath(dir_, "plots", "$(fl_)-iv2d.png"))
            if !(gl.iv isa BM.SomeData{Void})
                if gl.iv.vals isa VAWTools.Gridded
                    err = iv2d.v-gl.iv.vals.v
                    err[.!gl.ivmask] = NaN
                    mi, ma = VAWTools.extremanan(err)
                    clim = max(abs(mi), abs(ma)) # set clims symmetric because of colorscale
                    BM.plot2d_h(gl, err; plotiv=true, colorbar_title="iv model-data (m)",
                                seriescolor=:pu_or, clims=(-clim,clim), reuse=true, outlines=true)
                    Plots.savefig(joinpath(dir_, "plots", "$(fl_)-iv2d-err.png"))
                elseif gl.iv.vals isa VAWTools.Traj
                    #BM.plot2d_h(gl, iv2d.v; clims=:radar, plotiv=true, colorbar_title="iv model-data (m)")
                else
                    error()
                end
            end
            BM.plot2d_h(gl, iv2ds.v; plotiv=true, colorbar_title="iv std (m)", reuse=true,
                        plottraj=false, clims=:model, outlines=true)
            Plots.savefig(joinpath(dir_, "plots", "$(fl_)-iv2d-std.png"))
        end
        # this isn't quite correct as, e.g., hs1d of mtheta is not equal to expected hs1d
        #
        # NOTE this plots the mean hs2d, but the others are from running the model with mean paras
        BM.plot_run(fwdsol,
                    reuse=true, plot2d=[:h,:iv,:ivlog][1], clims=[:model,:radar][2],
                    plotline=I2.get_rlines(gid, expnr),
                    radarline_gl=gl)
        Plots.savefig(joinpath(dir_, "plots", "$(fl_)-run.png"))
    end
    if make_plots_thetas
        names = jldd["theta0_names"]
        theta0 = BM.Theta0{Float64}(mtheta, names,
                                    OrderedDict{Symbol,LinSpace{Float64}}(), gbn, pp, pm, pn)
        BM.plottheta(thetas, theta0, fsize=(1000,1000), reuse=true)
        Plots.savefig(joinpath(dir_, "plots", "$(fl_)-thetas.png"))
    end
    # PyPlot.close("all") does not help
    return data_out
end

# Good queries
# data[data[:gl].==:Unteraar,:]
# by(data, :gl, d->DataFrame( mean_rmse = mean(d[:rmse]), max_err = maximum(d[:max]), runs=length(d[:expnr])))

function summarize_data(data; expnrs=0:16, person=false)
    d = data[ in.(data[:expnr], (expnrs,)), :]
    bycols = person ? [:person, :gl] : [:gl]
    by(d, bycols, d->DataFrame(prio = round(Int, mean(d[:prio])),
                               median_mean_h =  median(d[:mean_h]),
                               quantiles_mean_h =  (summarystats(d[:mean_h]).q25, summarystats(d[:mean_h]).q75),
                               extrema_mean_h =  extrema(d[:mean_h]),
                               mean_rmse = mean(d[:rmse]),
                               rel_mean_rmse = mean(d[:rmse])/median(d[:mean_h]),
                               mean_max_err = mean(d[:max]),
                               maximum_max_err = maximum(d[:max]),
                               runs=length(d[:expnr]),
                               Rhat_median=median(d[:Rhat_max]),
                               Rhat_max=maximum(d[:Rhat_max])
                               ),
       sort = true)

    # Easier:
    # flds = [:gl, :rmse, :rmse_fitting, :max, :max_fitting, :Rhat, :Rhat_max]
    # sort(aggregate(d[flds], :gl, [mean, median, maximum]))
end

@everywhere begin
    include("itmix-setup.jl")
    using FileIO, DataFrames
    import Plots. PyPlot, Interpolations
    pl_kws = Dict{Symbol,Any}()
    PyPlot.matplotlib[:interactive](false) # memory leak https://github.com/JuliaPy/PyPlot.jl/issues/111

# global priority,gid,gl,gb,pp,pm,pn,pl,fl,fls,flj,hs2d,hs2ds,RMSE,max_err,smaller_20,thetas,theta0
# global fwdm_fn, jldd, gbn, mtheta, names

    gids = [BM.ITMIXv2_glaciers, BM.ITMIXGlacier.(BM.itmix_glaciers_iv,2)][2]
    if !isdefined(:gls) || length(gls)<length(gids)
        gls = try
            open("data/ITMIX2/results_Werder/$subdir/gls.jls", "r") do io
                deserialize(io)
            end
        catch
            gls = Dict()
            for gid in gids
                gls[gid] = BM.init_forward(gid; pl_kws...);
            end
            gls
        end
    end
end
# # load and save gls with
# open("data/ITMIX2/results_Werder/$subdir/gls.jls", "w") do io
#     serialize(io, gls)
# end
# gls = open("data/ITMIX2/results_Werder/$subdir/gls.jls", "r") do io
#     deserialize(io)
# end

if runit
    println("""
    ==================
    Running over directory $(joinpath(first(gls)[2][end].data_root,"results_Werder", subdir))
    ==================
    """)

    data = make_data()

    #gid = BITEModel.ITMIXGlacier(:Austfonna, 2)
    for priority=1:4, rm_iv in [false,true]
        println("\n\n++++++++++++++++ PRIORITY $priority ; rm_iv = $rm_iv ++++++++++++++++++++++++++\n\n")
        for gid in I2.priorities[priority] #[BM.ITMIXGlacier(:SouthGlacier,2)] #
            if !haskey(gls, gid)
                continue
            end
            gl,gb,pp,pm,pn,pl = gls[gid]

            dir_ = splitdir(I2.filepaths(gl, pl, 0, rm_iv, "results_Werder/"*subdir)[1])[1]
            !isdir(joinpath(dir_, "plots")) && mkdir(joinpath(dir_, "plots"))

            println("\n\n============================  Glacier $(BM.getname(gid))")
            #println("expnr, nl, rmse, rmse_fitting, max, Rhat, Rhat_max")
            if parallel
                out = pmap(expnr -> evaluate_one(expnr, rm_iv, gb, pp, pm, pn, pl;
                                                 make_plots=make_plots,
                                                 make_plots_thetas=make_plots,
                                                 subdir=subdir),
                           0:I2.nexps)
                [push!(data, o) for o in out if o!=nothing]
            else
                for expnr=0:I2.nexps

                println("\n\n                              Exp $expnr")
                    #try
                    data_out = evaluate_one(expnr, rm_iv, gb, pp, pm, pn, pl;
                                            make_plots=make_plots,
                                            make_plots_thetas=make_plots,
                                            subdir=subdir)
                    data_out!=nothing && push!(data, data_out)
                    # catch e
                    #     e isa InterruptException && rethrow(e)
                    #     print(e)
                # end
                end
            end
        end
    end
    # save data:
    FileIO.save(joinpath(pl.data_root,"results_Werder", subdir, "data.jld"), "data", data,
                "data_summarized", summarize_data(data)) #, "gls", gls) saving gls takes 3.5GiB!!!

    ## Load with
    # using DataFrames, FileIO, BITEModel
    # data, data_summarized = FileIO.load("data/ITMIX2/results_Werder/res-30-oct/data.jld", "data", "data_summarized")

end

# read others results

function read_others_results(gid, expnr, dir, person; prefix="thickness", thin=1, start=1)
    gl,gb,pp,pm,pn,pl = gls[gid]
    r_lines_fitting = I2.get_rlines(gid, expnr)

    expnr_ = expnr>9 ? Symbol("exp$expnr") : Symbol("exp0$expnr")
    fl = joinpath(dir, "$(prefix)_$(BM.getname(gid))_$(person)_$(expnr_).asc")
    hs2d = VAWTools.read_agr(fl)
    hs2d.v[hs2d.v.==0] = NaN
    if thin>1
        hs2d = VAWTools.downsample(hs2d, thin, start)
    end

    RMSE, max_err, quantiles, smaller_20, smaller_std, smaller_2std, err_h = BM.compare_to_radar(gl, hs2d, nothing, r_lines_fitting)
    mean_ice_thickness, volume_area = VAWTools.meannan(hs2d.v), BM.volume_area_mean_h(gl)

    data_out = [person, gid.id, -1, expnr, true, length(r_lines_fitting),
                map(x->signif(x,3), (mean_ice_thickness, volume_area, RMSE[1],
                                     RMSE[2], (max_err[1]), (max_err[2]), NaN, NaN))..., Float64[], err_h]
    return data_out
end


text:
using MPI
using PencilArrays
using Random
using Test

MPI.Initialized() || MPI.Init()

dims = (8, 12, 9)
comm = MPI.COMM_WORLD
pen = Pencil(dims, comm; permute = Permutation(2, 3, 1))

# This can be simplified to `redirect_stdout(devnull)` on Julia ≥ 1.6.
let dev_null = @static Sys.iswindows() ? "nul" : "/dev/null"
    MPI.Comm_rank(comm) == 0 || redirect_stdout(open(dev_null, "w"))
end

u = PencilArray{Float32}(undef, pen)
randn!(u)
ug = global_view(u)

@testset "Indices" begin
    for A in (u, ug)
        for (n, I) in zip(LinearIndices(A), CartesianIndices(A))
            @test A[n] === A[I]
        end

        for (n, v) in pairs(IndexLinear(), A)
            @test A[n] === v
        end

        for (I, v) in pairs(IndexCartesian(), A)
            @test A[I] === v
        end
    end
end


text:
module Helpers

export default_dev_dir, pkg_org, pkg_from_repo, pkg_url, pkg_giturl,
       pkg_names, pkg_parsebranch, fork_from_repo

const orgnames = Dict(
                     "AbstractAlgebra" => "Nemocas",
                     "Hecke" => "thofma",
                     "Nemo" => "Nemocas"
                )

const non_jl_repo = [ "libpolymake-julia", "libsingular-julia" ]

const default_dev_dir = "oscar-dev"

pkg_org(pkg::AbstractString) = get(orgnames, pkg, "oscar-system")

pkg_repo(pkg::AbstractString) = pkg in non_jl_repo ? pkg : "$pkg.jl"

pkg_from_repo(repo) = isnothing(repo) ? nothing : match(r"/([-_\w]+)(?:\.jl)?$",repo)[1]

fork_from_repo(repo) = isnothing(repo) ? nothing : match(r"^([-_\w]+)/",repo)[1]

function pkg_url(pkg::AbstractString; full=true, fork=nothing)
   org = isnothing(fork) ? pkg_org(pkg) : fork
   return (full ? "https://github.com/" : "") * org * "/" * pkg_repo(pkg)
end

function pkg_giturl(pkg::AbstractString; fork=nothing)
   org = isnothing(fork) ? pkg_org(pkg) : fork
   return "git@github.com:$org/" * pkg_repo(pkg)
end

function pkg_names(dir::AbstractString=default_dev_dir)
   return filter(x -> (isdir(joinpath(dir,x)) && x != "project"), readdir(dir))
end

function pkg_parsebranch(pkg::AbstractString, branch::AbstractString)
   fork = nothing
   if startswith(branch, "https://")
      urlmatch = match(r"https://github\.com/([-_\w]+)/[-_\w]+(?:\.jl)?#(.*)", branch)
      isnothing(urlmatch) &&
         @error "could not parse org and branch from $branch"
      if urlmatch[1] != pkg_org(pkg)
         fork = urlmatch[1]
      end
      branch = urlmatch[2]
   end
   if isnothing(fork)
      return (nothing, branch, "$pkg#$branch")
   else
      return (fork, branch, "$pkg@$fork#$branch")
   end
end

end


text:
"""
Measure Local sensitivity, Response coefficient and define fraction of activation
"""

module Sens

using Calculus
using PyCall
using SymPy
using Parameters

sp = pyimport("sympy")

include("Regression.jl")

export localSens, responseCoef, actFrac

@with_kw struct sensitivity_result
    f_x
    f_dx
    x
    f_response_coefficient
    activated_fraction
end

"""
Local Sensitivity Analysis

Compute the localsensitivity (df/ds) along the input s (number or vector).

# Arguments
- `f::function` : with single input f(s)
- `s`:: Single number (vector)
"""
function localSens(f, s)
    df = derivative(f)
    return df(s)
end

function localSens(f, s)
    df = derivative(f)
    return df.(s)
end

function localSens(f, s::Sym)
    df = sp.diff(f,s)
    return sp.simplify(df)
end

"""
Response Coefficient

Compute the response coefficient: \$R^{x}_{s} = \\frac{dx}{ds} \\frac{s}{x}\$.

where \$R^{x}_{s}\$ is the response coefficient of input \$s\$ and output \$x\$.

# Arguments
- `f::function` : with single input x = f(s)
- `s`:: Single number (can also be vector or sympy symbol).

# Warning

The result may be fluctuated due to the sinularity point.

# Reference:
1. Ferrell, J. E. & Ha, S. H. Ultrasensitivity part I: Michaelian responses and zero-order ultrasensitivity. Trends in Biochemical Sciences 39, 496–503 (2014).
2. Legewie, S., Blüthgen, N. & Herzel, H. Quantitative analysis of ultrasensitive responses. FEBS J. 272, 4071–4079 (2005).
"""
function responseCoef(f, s)
    df = derivative(f)
    return df(s) * s / f(s)
end

function responseCoef(f, s)
    df = Calculus.derivative(f)
    return df.(s) .* s ./ f.(s)
end

function responseCoef(f,s::Sym)
    df = sp.diff(f,s)
    return sp.simplify(df * s / f)
end

"""
Compute `activation fraction` of the stimulation process.

# Reference
1. Legewie, S., Blüthgen, N. & Herzel, H. Quantitative analysis of ultrasensitive responses. FEBS J. 272, 4071–4079 (2005).
"""
actFrac(x, x_basal, x_max) = (x-x_basal) / (x_max - x_basal)

function relative_amplification_approach(f,x)
    df = Calculus.derivative(f) # Differential function
    f_x = f.(x) # f(x)
    dfdx = df.(x) # df/dx
    f_resp_coef = dfdx .* x ./ f_x # response coefficient

    y_max = maximum(f_x)
    y_min = minimum(f_x)

    act_fraction = actFrac.(f_x, y_min, y_max)

    return sensitivity_result(
        f_x=f_x,
        f_dx= dfdx,
        x=x,
        f_response_coefficient= f_resp_coef,
        activated_fraction=act_fraction
    )

end

end


text:
using NCDatasets
using Test

new_cache_size = 32000000
new_cache_nelems = 2000
new_cache_preemption = .75

NCDatasets.nc_set_chunk_cache(new_cache_size,new_cache_nelems,new_cache_preemption)
new_cache_size2,new_cache_nelems2,new_cache_preemption2 = NCDatasets.get_chunk_cache()

@test new_cache_size2 == new_cache_size
@test new_cache_nelems2 == new_cache_nelems
@test new_cache_preemption2 ≈ new_cache_preemption

newest_cache_size = 32000001
NCDatasets.set_chunk_cache(size = newest_cache_size)
new_cache_size2,new_cache_nelems2,new_cache_preemption2 = NCDatasets.get_chunk_cache()

@test new_cache_size2 == newest_cache_size
# should not change
@test new_cache_nelems2 == new_cache_nelems
@test new_cache_preemption2 ≈ new_cache_preemption


text:
using HeartBeats, Test

ecg = example_ecg()
fs = 360

@test length(detect_heartbeats(ecg, fs)) == 478


text:
export ## Types
       CompositeQSystem,
       Control,
       Dissipation,
       Field,
       IndexSet,
       Interaction,
       ParametricInteraction,
       QSystem,
       ## Methods
       hamiltonian,
       dim,
       label,
       lowering,
       number,
       raising,
       strength,
       X,Y

abstract type QSystem end

label(q::QSystem) = q.label
dim(q::QSystem) = q.dim
raising(q::QSystem) = diagm(-1 => sqrt.(1:(dim(q)-1)))
lowering(q::QSystem) = diagm(1 => sqrt.(1:(dim(q)-1)))
number(q::QSystem) = raising(q) * lowering(q)
X(q::QSystem) = raising(q) + lowering(q)
Y(q::QSystem) = 1im*(raising(q) - lowering(q))
hamiltonian(q::QSystem, t) = hamiltonian(q)

abstract type Control end
label(c::Control) = c.label

abstract type Interaction end
strength(i::Interaction) = i.strength
hamiltonian(i::Interaction, t) = hamiltonian(i)

abstract type ParametricInteraction <: Interaction end

abstract type Dissipation end

const IndexSet = Vector{Int}
const ExpansionIndices = Tuple{Vector{IndexSet},Vector{IndexSet},Vector{IndexSet}}

# [todo] - Should COmpositeQSystem <: QSystem ?
mutable struct CompositeQSystem
    # [feature] - Use something like OrderedDict for component enumeration
    subSystems::Vector{QSystem}
    interactions::Vector{Interaction}
    parametericInteractions::Vector{ParametricInteraction}
    subSystemExpansions::Vector{ExpansionIndices}
    interactionExpansions::Vector{ExpansionIndices}
    dissipatorExpansions::Vector{ExpansionIndices}
    dissipators::Vector{Dissipation}
end


text:
@registerN DiscreteChebyshev AbstractCDOP2 α β
export DiscreteChebyshev

"""
    DiscreteChebyshev

This uses p22 of [Koepf and Schmersau](https://arxiv.org/pdf/math/9703217.pdf)  to define a two-parameter  family of *non* orthogonal polynomials.
See  the example  in [`DiscreteWeightFunction`](@ref) for implementing  the [DiscreteChebyshev](https://en.wikipedia.org/wiki/Discrete_Chebyshev_polynomials)  polynomials  from Wikipedia.

```jldoctest
julia> using Polynomials, SpecialPolynomials

julia> import SpecialPolynomials: Δₓ, ∇ₓ

julia> α,β = 1/2, 1
(0.5, 1)

julia> P  = DiscreteChebyshev{α,β}
DiscreteChebyshev{0.5, 1}

julia> i = 5
5

julia> yᵢ = basis(P, i)
DiscreteChebyshev(1.0⋅K⁽ᵅᵝ⁾₅(x))

julia> x = variable(P)
DiscreteChebyshev(- 2.0⋅K⁽ᵅᵝ⁾₀(x) + 2.0⋅K⁽ᵅᵝ⁾₁(x))

julia> a,b,c,d,e = SpecialPolynomials.abcde(P)
(a = 0, b = 0, c = 1, d = 0.5, e = 1)

julia> λᵢ  = -(a*i*(i-1)  + d*i)
-2.5

julia> Δₓ(∇ₓ(yᵢ)) +  (α*x + β) * Δₓ(yᵢ) ≈ -λᵢ*yᵢ # p22: "are not orthogonal, but satisfy the difference equation..."
true
```

"""
DiscreteChebyshev

basis_symbol(::Type{<:DiscreteChebyshev{α,β}}) where {α,β} = "K⁽ᵅᵝ⁾"
Polynomials.domain(::Type{<:DiscreteChebyshev{α,β}}) where {α,β} =
    Polynomials.Interval(-Inf, Inf)
abcde(::Type{<:DiscreteChebyshev{α,β}}) where {α,β} =
    NamedTuple{(:a, :b, :c, :d, :e)}((0, 0, 1, α, β))

function k1k0(P::Type{<:DiscreteChebyshev{α,β}}, n::Int) where {α,β}
    α
end


text:
module BooSTjl

export BooST, BooSTMore, estimate_derivatives, predictBooST

using StatsBase
using Base.LinAlg.BLAS

include("auxiliary.jl")
include("export.jl")

end # module

max_stars_repo_name: gabrielrvsc/BooSTjl

max_issues_repo_name: gabrielrvsc/BooSTjl

max_forks_repo_name: gabrielrvsc/BooSTjl


text:
__precompile__()
module TerminalMenus

terminal = nothing  # The user terminal

function __init__()
    global terminal
    terminal = Base.Terminals.TTYTerminal(get(ENV, "TERM", is_windows() ? "" : "dumb"), STDIN, STDOUT, STDERR)
end

include("util.jl")
include("config.jl")

include("AbstractMenu.jl")
include("RadioMenu.jl")
include("MultiSelectMenu.jl")

export
    RadioMenu,
    MultiSelectMenu,
    request

end # module


text:
# Solver.jl by Alex Van Mechelen
# Creation Date: 23/04/2020
# Made for Julia 1.2.0

using Main

global firstTime = true
global ignoreEasySteps = false

function fillPMs!(game::Board,allOptions::Bool=true)
	for i in 1:length(game.grid)
		game.pencilMarks[i] = (allOptions && game.grid[i] == 0) ? Array{Int,1}(game.dim:-1:1) : zeros(Int,(game.dim))
	end
end

function remPM!(game::Board,xPos::Int,yPos::Int,height::Int)
	count(game.pencilMarks[xPos,yPos].==height) == 0 && return false
	pos = findall(game.pencilMarks[xPos,yPos].==height)[1]
	game.pencilMarks[xPos,yPos][pos] = 0
	sort!(game.pencilMarks[xPos,yPos],rev=true)
	return true
end

function remPM!(pm::Int,PMs::Array{Int,1})
	try
		pos = findall(PMs.==pm)[1]
		PMs[pos] = 0
		sort!(PMs,rev=true)
		return true
	catch # Value already removed
		return false
	end
end

function resetPM!(PMs::Array{Int,1})
	PMs = zeros(Int,length(PMs))
end

function row(game::Board,xPos::Int,orientTLBR::Bool=true)
	return orientTLBR ? (game.grid[xPos,:],game.pencilMarks[xPos,:]) : (reverse(game.grid[xPos,:]),reverse(game.pencilMarks[xPos,:]))
end

function column(game::Board,yPos::Int,orientTLBR::Bool=true)
	return orientTLBR ? (game.grid[:,yPos],game.pencilMarks[:,yPos]) : (reverse(game.grid[:,yPos]),reverse(game.pencilMarks[:,yPos]))
end

function setRow!(row::Array{Int,1},PMs::Array{Array{Int,1},1},game::Board,xPos::Int,orientTLBR::Bool=true)
	if !orientTLBR
		row = reverse(row)
		PMs = reverse(PMs)
	end
	game.grid[xPos,:] = row
	game.pencilMarks[xPos,:] = PMs
end

function setColumn!(column::Array{Int,1},PMs::Array{Array{Int,1},1},game::Board,yPos::Int,orientTLBR::Bool=true)
	if !orientTLBR
		column = reverse(column)
		PMs = reverse(PMs)
	end
	game.grid[:,yPos] = column
	game.pencilMarks[:,yPos] = PMs
end

function explainShowStep(showStep::Bool,explanation::String,params::Tuple{Board,Int,Bool},outData::Tuple{Array{Int,1},Array{Array{Int,1},1}},rowHl::Bool,colHl::Bool,switchHl::Bool,convHlL::Int,lineIndex::Int=1)
	game = params[1]

	if switchHl
		lineIndex = game.dim-lineIndex+1
	end

	if rowHl
		game.highLight = [convHlL,lineIndex]
	elseif colHl
		game.highLight = [lineIndex,convHlL]
	end

	if showStep
		print(game)
		printstyled(explanation*".\n",color = COL_EXPL)
		reprompt = true
		while reprompt
			input = lowercase(readline())
			if input == "ignore" || input == "show"
				global ignoreEasySteps = (input == "ignore")
				print(game)
				if input == "ignore"
					printstyled("Succesfully enabled ignoreEasySteps.\n",color = COL_QUESTION)
				else
					printstyled("Succesfully disabled ignoreEasySteps.\n",color = COL_QUESTION)
				end
			else
				reprompt = false
			end
		end
		Calibrate() # Screen size could have changed while prompting user
	end
end

# Main function that solves the board
function mainSolver(game::Board,showStep::Bool=false)
	# Initialize variables
	fillPMs!(game)
	succes = true
	memory = Array{Union{Array{Int,1},Board},1}()

	# Loop while there are still errors
	while !getErrors(game)[3]
		# Perform simple logical steps
		errorInLine,succes,game = loopAroundGrid!(showStep,game) # CheckLine

		# If simple step does not have succes, perform complexStep (before guessing)
		if !succes
			errorInLine,succes,game = loopAroundGrid!(showStep,game,true) #ComplexStep
		end

		if (!getErrors(game)[3] && (errorInLine || !succes || count(getErrors(game)[1].==true)>0)) # Game not solved && ErrorInLine, no progress has been made or encountered an error inside the grid

			# If errorInLine, the grid is fully filled in or there occurred errors inside the grid
			if errorInLine || count(game.grid.==0) == 0 || count(getErrors(game)[1].==true)>0
				if length(memory) < 1 # Game unsolvable
					global firstTime = true # Reset for next solve
					return false,game
				end

				# Get data of last picked PM
				lastPicked,numLeftPMs,xPosPicked,yPosPicked = memory[end]

				while numLeftPMs < 1 # One PM left, but appears to be incorrect since !succes
					if length(memory) < 1 # No memory options left
						global firstTime = true # Reset for next solve
						return false,game
					end

					# Remove last two values in memory
					memory = memory[1:end-2]

					# If memory smaller than one, game unsolvable
					if length(memory) < 1 # Game unsolvable
						return false,game
					end

					# Reset the game to last (assumed to be) correct game-state (as far as it knows now)
					game = deepcopy(memory[end-1])

					# Get data of last picked PM
					lastPicked,numLeftPMs,xPosPicked,yPosPicked = memory[end]
				end

				# Reset the game to last (assumed to be) correct game-state (as far as it knows now)
				game = deepcopy(memory[end-1])

				# Remove last picked PM from last game-state in memory
				remPM!(game,xPosPicked,yPosPicked,lastPicked)
				memory[end-1] = deepcopy(game)

				# Get options for new pick
				copyLastPMs = deepcopy(game.pencilMarks[xPosPicked,yPosPicked])
				copyNewPickOptions = deleteat!(copyLastPMs,findall(copyLastPMs.==0))

				# Pick a random new option
				newPick = rand(copyNewPickOptions)

				# Replace previous pick-data by new pick-data
				numPMsLeft = count(game.pencilMarks[xPosPicked,yPosPicked].!=0) - 1 # Number of PM's in the cell - 1
				memory[end] = [newPick,numPMsLeft,xPosPicked,yPosPicked]

			else # Pick a new PM somewhere and store to memory

				# Make a hypothesis
				found = false
				minimizeOptions = 2

				while !found
					# Loop over the whole game board
					for xPos in 1:game.dim
						for yPos in 1:game.dim
							game.grid[xPos,yPos] != 0 && continue
							numPMs = count(game.pencilMarks[xPos,yPos].!=0)
							if 0 < numPMs <= minimizeOptions
								# Choose a random PM in this cell
								copyPMs = deepcopy(game.pencilMarks[xPos,yPos])
								pick = rand(deleteat!(copyPMs,findall(copyPMs.==0)))

								# Set highLight to the cell you chose a PM value for
								game.highLight = [xPos,yPos]

								# Make this cell a guessed value
								game.guessVals[xPos,yPos] = true

								# Save the current game to memory
								push!(memory,deepcopy(game))

								# Save picked PM, number of PMs left and location to memory
								push!(memory,[pick,numPMs-1,xPos,yPos])

								# Give game cell the picked PM value
								game.grid[xPos,yPos] = pick

								found = true
								break
							end
						end
						found && break
					end

					# If not found, increase minimizeOptions by 1
					if !found
						if minimizeOptions > game.dim
							global firstTime = true # Reset for next solve
							return false,game
						else
							minimizeOptions += 1
						end
					end
				end
				# Picked a clue, new game board -> perform standard steps (loopAroundGrid())
			end
		end
	end
	# Game solved
	global firstTime = true # Reset for next solve
	game.highLightLine = 0 # Reset highLightLine to highLight no line
	return true,game
end

# Loops around grid and solves in rows/columns
function loopAroundGrid!(showStep::Bool,game::Board,performComplexStep::Bool=false)
	succes = false
	errorInLine = false
	for i in 1:game.dim
		for dir in 1:2
			params = (game,i,dir==1)

			game.highLightLine = i+(dir+1)*game.dim
			errorInLine1,succes1,outLine,outPMs = performComplexStep ? complexStep(showStep,game.clues[i,dir],game.clues[i,dir%2+1],column(params...)...,params) : checkLine(showStep,game.clues[i,dir],game.clues[i,dir%2+1],column(params...)...,params)
			succes1 && setColumn!(outLine,outPMs,params...)

			game.highLightLine = i+(dir-1)*game.dim
			errorInLine2,succes2,outLine,outPMs = performComplexStep ? complexStep(showStep,game.clues[i,dir+2],game.clues[i,dir%2+3],row(params...)...,params) : checkLine(showStep,game.clues[i,dir+2],game.clues[i,dir%2+3],row(params...)...,params)
			succes2 && setRow!(outLine,outPMs,params...)

			errorInLine |= errorInLine1 | errorInLine2
			succes |= succes1 || succes2
		end
	end
	global firstTime = false
	return errorInLine,succes,game
end

function checkVisible(line::Array{Int,1})
	visible = 0
	if count(line.!=0)>0 # If there are cells filled in
		maxCell = 0
		for cell in line
			if cell > maxCell
				visible += 1
				maxCell = cell
			end
		end
	end
	return visible
end

function permutations(possibilityLine::Array{Array{Int,1},1},permList::Array{Array{Int,1},1}=Array{Array{Int,1},1}(),index::Int=1)
    if index == length(possibilityLine)
        return permList
    end

    if index == 1
        for i in possibilityLine[index]
            push!(permList,[i])
        end
    end

    for partialPossibleLineIndex in 1:length(permList)
        front = permList[partialPossibleLineIndex]
        for i in possibilityLine[index+1]
            if front == permList[partialPossibleLineIndex]
                permList[partialPossibleLineIndex] = [front...,i]
            else
                push!(permList,[front...,i])
            end
        end
    end
    return permutations(possibilityLine,permList,index+1)
end

function complexStep(showStep::Bool,clue::Int,oppClue::Int,line::Array{Int,1},PMs::Array{Array{Int,1},1},params::Tuple{Board,Int,Bool})

	succes = false

	# If line is solved, return same line
	count(line.==0)==0 && return false,false,line,PMs

	game = params[1]
	lastParams = params[2:end]

	convHlL,rowLHl,rowRHl,colTHl,colBHl = convHlLine(game,game.highLightLine)
	rowHl = rowLHl || rowRHl
	colHl = colTHl || colBHl
	switchHl = rowRHl || colBHl

	outLine = line
	outPMs = PMs

	dim = game.dim

	# Construct possibilityLine
	possibilityLine = Array{Array{Int,1},1}(undef,dim)
	PMpositionsInLine = Array{Bool,1}(undef,dim)
	for i in 1:dim
		# If cell-value exists, this is the only possibility for that cell, else outPMs
		if line[i]!=0
			possibilityLine[i] = [line[i]]
			PMpositionsInLine[i] = false
		else
			possibilityLine[i] = filter!(x -> x != 0, deepcopy(outPMs[i]))
			PMpositionsInLine[i] = true
		end
	end

	# Get all permutations of this possibilityLine
	permList = permutations(possibilityLine)
	outPermList = Array{Array{Int,1},1}()

	# Exclude impossible permutations
	for i in 1:length(permList)
		if length(unique(permList[i]))==dim && (clue==0 || towersVisible(permList[i])==clue) && (oppClue==0 || towersVisible(reverse(permList[i]))==oppClue)
			push!(outPermList,permList[i])
		end
	end

	# If there are no possible permutations (last guess was incorrect)
	if isempty(outPermList)
		# Return errorInLine,succes,line,PMs
		return true,false,line,PMs
	end

	# If a cell is empty (no value nor PMs), last guess was incorrect
	for cellIndex in 1:dim
		if outPMs[cellIndex] == zeros(Int,dim) && outLine[cellIndex] == 0
			# Return errorInLine,succes,line,PMs
			return true,false,line,PMs
		end
	end

	# For all PMpositionsInLine
	for i in findall(PMpositionsInLine)
		onlyOnePossible = true
		possibleVals = [outPermList[1][i]]
		for permutation in outPermList
			# If other value occurred, not onlyOnePossible
			if !(permutation[i] in possibleVals)
				onlyOnePossible && (onlyOnePossible = false)
				push!(possibleVals,permutation[i])
			end
		end

		# If only clue possible, this should be the cell's value
		if onlyOnePossible && possibleVals[1] != 0
			outLine[i] = possibleVals[1]
			explainShowStep(showStep,"Only possible cell value in all permutations of this line is $(possibleVals[1])",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL,i)
			succes |= true
		end

		# If less possibleVals than PMs in that cell, remove the ones that are not possible
		if length(possibleVals) < count(outPMs[i].!=0)
			# Reset outPMs[i] to all zeros
			outPMs[i] = zeros(Int,dim)

			# Set outPMs[i] to possibleVals (sorted) (and the remaining places zero)
			outPMs[i][1:length(possibleVals)] = sort(possibleVals,rev=true)

			explainShowStep(showStep,"Removing Pencil-Marks that never occur in any permutation of this line",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL,i)
			succes |= true
		end
	end

	# Check if error in line
	errorInLine = count(getDoubles(outLine).!=0)>0

	return errorInLine,succes,outLine,outPMs
end

function convHlLine(game::Board,hlL::Int)
	dim = game.dim

	# Convert hlL index to line/col index
	convHlL = (hlL-1)%dim+1

	# Check what line/col to highLight
	rowLHl = (hlL in 1:dim)
	rowRHl = (hlL in dim+1:2*dim)
	colTHl = (hlL in 2*dim+1:3*dim)
	colBHl = (hlL in 3*dim+1:4*dim)

	return convHlL,rowLHl,rowRHl,colTHl,colBHl
end

function checkLine(showStep::Bool,clue::Int,oppClue::Int,line::Array{Int,1},PMs::Array{Array{Int,1},1},params::Tuple{Board,Int,Bool})

	succes = false
	pSucces = false

	# If line is solved, return same line
	count(line.==0)==0 && return false,false,line,PMs

	game = params[1]
	lastParams = params[2:end]

	convHlL,rowLHl,rowRHl,colTHl,colBHl = convHlLine(game,game.highLightLine)
	rowHl = rowLHl || rowRHl
	colHl = colTHl || colBHl
	switchHl = rowRHl || colBHl

	outLine = line
	outPMs = PMs

	dim = game.dim

	if firstTime  # Only run first time
		# Tallest tower cannot occur in the first (clue-1) positions, the second tallest tower cannot occur in the first (clue-2) positions, ect.
		for height in 1:dim
			forbiddenPositions = clue+height-dim-1
			forbiddenPositions < 1 && continue
			for cellIndex in 1:forbiddenPositions
				pSucces |= remPM!(height,outPMs[cellIndex])
				!ignoreEasySteps && pSucces && explainShowStep(showStep,"Tallest tower cannot occur in the first (clue-1) positions, the second tallest tower cannot occur in the first (clue-2) positions, ect.",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL,cellIndex) # Shorten
				succes |= pSucces
			end
		end

		# For clue of 1, first tower must be of height dim
		if clue == 1 && outLine[1] != dim
			outLine[1] = dim
			resetPM!(outPMs[1])
			explainShowStep(showStep,"Only one tower can be visible, so the first tower must be the tallest tower",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL)
			succes |= true
		end

		# For clue of 2, second tower cannot be the second tallest one (this would result in either 1 or three towers visible)
		if clue == 2
			pSucces = remPM!(dim-1,outPMs[2])
			pSucces && explainShowStep(showStep,"For clues of 2, second tallest tower cannot be the second tallest one",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL,2)
			succes |= pSucces
		end

		# If clue = dim, all towers are visible in ascending order
		if clue == dim
			outLine = Array{Int,1}(1:dim)
			for PM in outPMs
				resetPM!(PM)
			end
			explainShowStep(showStep,"If clue = dim, all towers are visible in ascending order",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL)
			succes |= true
		end
	end

	if clue > 1
		# For $(clue) of 2, if first tower is smallest one, second tower must be tallest one
		if clue == 2 && outLine[1] == 1 && outLine[2] != dim
			outLine[2] = dim
			resetPM!(outPMs[2])
			explainShowStep(showStep,"For clue of 2, if first tower is smallest one, second tower must be tallest one",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL,2)
			succes |= true
		end

		# If tallest tower at $(clue)-th position and all positions behind that tower are occupied by the towers from the tallest to the clue-th tallest tower in any order, then the cells before the tallest tower contain the sequence of shortest to clue-1-th tower
		if outLine[clue] == dim && sort(outLine[clue+1:dim]) == Array{Int,1}(clue:dim-1) && outLine[1:clue-1] != Array{Int,1}(1:clue-1)
			outLine[1:clue-1] = Array{Int,1}(1:clue-1)
			for PM in outPMs[1:clue-1]
				resetPM!(PM)
			end
			explainShowStep(showStep,"Tallest tower at clue-th position and all clue highest towers (in any order) filled in after the tallest one, means the leftover ones should form an upcounting sequence.",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL)
			succes |= true
		end

		# If the first clue-1 towers have ascending values starting with 1 (step 1), then the next tower must be the tallest tower
		if outLine[1:clue-1] == Array{Int,1}(1:clue-1) && outLine[clue] != dim
			outLine[clue] = dim
			resetPM!(outPMs[clue])
			explainShowStep(showStep,"If the first clue-1 towers have ascending values starting with 1, then the next tower must be the tallest tower",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL,clue)
			succes |= true
		end

		# If clue is 1 smaller than dim and second cell is a one, then there is only one possible solution for that line
		if clue == dim-1 && outLine[2] == 1
			pSucces = (outLine[1] != 2 || outLine[3:dim] != 3:dim)
			outLine[1] = 2
			outLine[3:dim] = 3:dim
			for PMs in outPMs
				resetPM!(PMs)
			end
			pSucces && explainShowStep(showStep,"If clue is 1 smaller than dim and second cell is a one, then there is only one possible solution for that line (up-counting)",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL)
			succes |= pSucces
		end

		# If the first cell in this line is still empty
		if outLine[1] == 0

			filledPackedFromBack = 0
			filledFromBack = 0
			towerHeight = dim
			countingPacked = true
			for i in dim:-1:1
				outLine[i] == 0 && (countingPacked = false)
				if outLine[i] == towerHeight
					countingPacked && (filledPackedFromBack += 1)
					filledFromBack += 1
					towerHeight -= 1
				end
			end

			# If first cell empty and filledPackedFromBack = clue-1, first tower must be highest remaining tower
			if filledPackedFromBack == clue-1
				for height in dim:-1:1
					height in outLine && continue
					pSucces = (outLine[1] != height)
					outLine[1] = height
					resetPM!(outPMs[1])
					pSucces && explainShowStep(showStep,"If first cell empty and filledPackedFromBack = clue-1, first tower must be the highest remaining tower",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL)
					succes |= pSucces
					break
				end
			end

			# If filledFromBack = clue-1 and the towers before the clue-1-th highest tower are empty, then those towers (except the closest one to the clue) can't have height dim-clue+1
			if filledFromBack == clue-1
				countIndex = 1
				for i in 1:dim
					outLine[i] != 0 && break
					countIndex += 1
				end
				if outLine[countIndex] == dim-clue+2
					for i in 2:countIndex-1
						pSucces |= remPM!(dim-clue+1,outPMs[i])
					end
					pSucces && explainShowStep(showStep,"If filledFromBack = clue-1 and the towers before the clue-1-th highest tower are empty, then those towers (except the closest one to the clue) can't have height dim-clue+1",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL,countIndex-1)
					succes |= pSucces
				end
			end

			# If first cell is zero, second cell and tallest tower are filled in, and there is a visible sequence of clue number of towers, then the first tower must have a higher value than the second one
			if outLine[2] != 0 && dim in outLine && towersVisible(outLine) == clue
				for PM in 1:outLine[2]-1
					pSucces |= remPM!(PM,outPMs[1])
				end
				pSucces && explainShowStep(showStep,"If first cell is zero, second cell and tallest tower are filled in, and there is a visible sequence of clue number of towers, then the first tower must have a higher value than the second one",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL)
				succes |= pSucces
			end

			upperPresent = true
			for i in dim-clue+2:dim
				upperPresent &= i in outLine
			end

			if upperPresent
				counting = false
				currentMaxHeight = dim-clue+2
				inOrder = true
				for i in 2:dim
					currentMaxHeight == dim && break
					if counting && outLine[i] == currentMaxHeight + 1
						currentMaxHeight = outLine[i]
						continue
					end
					outLine[i] == dim-clue+2 && (counting = true)
				end

				if currentMaxHeight != dim
					succes |= remPM!(dim-clue+1,outPMs[1])
					succes && explainShowStep(showStep,"If the clue-1 highest towers are present, not in ascending order, then the first tower cannot be the clue-th highest tower",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL)
				end
			end
		end
	end

	# Only one PM on a cell means it should have that value
	for i in 1:dim # Loop over cells
		outLine[i] != 0 && continue
		if count(outPMs[i].!=0) == 1 && outLine[i] != outPMs[i][findall(outPMs[i].!=0)[1]]
			outLine[i] = outPMs[i][findall(outPMs[i].!=0)[1]]
			resetPM!(outPMs[i])
			explainShowStep(showStep,"Only one Pencil-Mark left, so this tower must be of height $(outPMs[i][findall(outPMs[i].!=0)[1]])",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL,i)
			succes |= true
		end
	end

	# Only cell in a line containing a certain PM gets that value
	for i in 1:dim # Loop over possible values
		i in outLine && continue
		index = 0 # Reset (zero -> no match found)
		for j in 1:dim # Loop over cells in line
			outLine[j] != 0 && continue
			if i in outPMs[j]
				if index == 0
					index = j
				else
					index = 0
					break
				end
			end
		end
		if index != 0 && outLine[index] != i
			outLine[index] = i
			resetPM!(outPMs[index])
			explainShowStep(showStep,"Only tower in this line with Pencil-Mark $(i)",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL,index)
			succes |= true
		end
	end

	# Remove PMs with the value of a nonzero value in that line
	for i in 1:dim # Loop over all the cells in this line
		outLine[i] == 0 && continue
		pSucces = false
		for cellPMs in outPMs
			pSucces |= remPM!(outLine[i],cellPMs)
			succes |= pSucces
		end
		!ignoreEasySteps && pSucces && explainShowStep(showStep,"Removing Pencil-Marks of each cell in this line with value $(outLine[i])",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL,i)
	end

	# When two cells in the line have matching Pencil-Marks (and two of them), those values can be excluded in Pencil-Marks of the other cells in that line
	for i in 1:dim
		outLine[i] != 0 && continue
		pSucces = false
		if count(outPMs.==outPMs[i])==2 && count(outPMs[i].!=0)==2
			for cellPMs in outPMs
				for PM in findall(outPMs[i].!=0)
					pSucces |= remPM!(PM,cellPMs)
				end
			end
		end
		pSucces && explainShowStep(showStep,"When two cells in the line have matching Pencil-Marks (and two of them), those values can be excluded in Pencil-Marks of the other cells in that line",(game,lastParams...),(outLine,outPMs),rowHl,colHl,switchHl,convHlL,i)
	end

	# Check if error in line
	errorInLine = count(getDoubles(outLine).!=0)>0

	return errorInLine,succes,outLine,outPMs
end

max_stars_repo_path: Towers Game/Solver.jl

max_issues_repo_path: Towers Game/Solver.jl

max_forks_repo_path: Towers Game/Solver.jl


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

# TODO/missing:
#  assertion on the lifting
#  filtering of the roots (maybe we don't want all?)
#  find all/ one/ few roots
#  break down in various modules:
#  - find powers of ideal
#  - find LLL/HNF basis
#  - lifting?
#  use ResRing(Poly) instead of doing % pgg?
#  an open variant where k is increased until we have a root?
#
#missing: different (better) handling for x^n-a
# possibly based on https://scicomp.stackexchange.com/questions/18846/newton-iteration-for-cube-root-without-division
#
# g = x^-n -a^(n-1), and r = root(g), then r = a^(-(n-1)/n), hence
#                                         ar = a^(-(n-1)/n + 1) = a^(1/n)
# Here, the Newton-iteration is easy: (b = a^(n-1))
# g = x^-n - b, => g' = -nx^(-n-1)
#Newton: r -> r-g(r)/g'(r) = r - (r^(-n) -b)/(-n r^(-n-1))
#                          = r - r/(-n) - b/n r^(n+1)
#                          = r(1+1/n) - b/n r^(n+1)
# no division! and no double iteration
# however, for reconstruction: we don't want r (which is non-integral)
# but ar which is...

################################################################################
#
#  Helper functions
#
################################################################################

Base.isless(x::arb, y::arb) = x < y

function lift(R::NmodPolyRing, a::fq_nmod)
  f = R()
  for i=0:degree(parent(a))-1
    setcoeff!(f, i, _get_coeff_raw(a, i))
  end
  return f
end

function lift(R::GFPPolyRing, a::fq_nmod)
  f = R()
  for i=0:degree(parent(a))-1
    setcoeff!(f, i, _get_coeff_raw(a, i))
  end
  return f
end

function (Zx::FmpzPolyRing)(a::nf_elem)
  b = Zx()
  @assert denominator(a) == 1
  if degree(parent(a)) == 1
    # If the number field is linear, then a.elem_length is not properly
    # initialized, that is, it could be anything.
    setcoeff!(b, 0, numerator(coeff(a, 0)))
  else
    for i=0:a.elem_length
      setcoeff!(b, i, numerator(coeff(a, i)))
    end
  end
  return b
end

################################################################################
#
#  Root computation via Hensel lifting
#
################################################################################

# General comments
# ispure means that f = X^deg(f) + coeff(f, 0)
# isnormal means that if there are either no or all roots

function _roots_hensel(f::Generic.Poly{nf_elem};
                       max_roots::Int = degree(f),
                       ispure::Bool = false,
                       isnormal::Bool = false,
                       root_bound::Vector{fmpz} = fmpz[])
  # f must be squarefree
  # I should check that
  K::AnticNumberField = base_ring(f)

  if iszero(trailing_coefficient(f))
    rs = nf_elem[zero(K)]
    f = div(f, gen(parent(f)))
  else
    rs = nf_elem[]
  end

  if length(rs) >= max_roots
    return rs
  end

  n = degree(K)
  deg = degree(f)

  # First we find a prime ideal such that f is squarefree modulo P
  # (The discriminant of f has only finitely many divisors).

  p = degree(f)+1
  deg_p = 1

  good_deg_p = 0

  good_p = p

  found = false

  local factor_of_g::gfp_poly

  local good_fp::fq_nmod_poly

  current_num_rmodp = degree(f) + 1

  not_better = 1

  coeffs_f = Vector{elem_type(base_ring(f))}(undef, length(f))

  if n >= 10
    STABILIZED = 2
  else
    STABILIZED = 0
  end

  if !ispure
    for i in 0:degree(f)
      coeffs_f[i + 1] = coeff(f, i)
    end
  end

  while !found
    p = next_prime(p)

    Rp = Nemo.GF(p, cached=false)
    Rpt, t = PolynomialRing(Rp, "t", cached=false)
    gp = Rpt(K.pol)
    if iszero(discriminant(gp))
      continue
    end

    lp = factor(gp).fac

    #set up the mod p data:
    #need FiniteField as I need to factor (roots)
    # I want to find a residue field with less roots
    for gp_factor in keys(lp)
      deg_p = degree(gp_factor)

      S = FqNmodFiniteField(gp_factor, :z, false)
      ST, T = PolynomialRing(S, "T", cached=false)

      if ispure
        red_coeff1 = Vector{fq_nmod}(undef, length(f))
        red_coeff1[length(f)] = S(1)
        for i = 2:length(f)-1
          red_coeff1[i] = zero(S)
        end
        red_coeff1[1] = S(Rpt(coeff(f, 0)))
        fp = ST(red_coeff1)
      else
        red_coeff = Vector{fq_nmod}(undef, length(f))
        for i in 1:length(f)
          red_coeff[i] = S(Rpt(coeffs_f[i]))
        end
        fp = ST(red_coeff)
      end

      if !issquarefree(fp)
        continue
      end
      num_rmodp = degree(gcd(fp, powmod(T, fmpz(p)^deg_p, fp) - T))

      if num_rmodp == 0
        return nf_elem[]
      end

      if isnormal && num_rmodp < degree(f)
        return nf_elem[]
      end

      if num_rmodp < current_num_rmodp || (num_rmodp == current_num_rmodp && deg_p > good_deg_p)
        good_p = p
        factor_of_g = gp_factor
        current_num_rmodp = num_rmodp
        good_fp = fp
        good_deg_p = deg_p
      else
        not_better = not_better + 1
      end

      if not_better >= STABILIZED
        found = true
        break
      end
    end
  end

  # compute the lifting exponent a la Friedrich-Fieker

  E = EquationOrder(K)
  r1, r2 = signature(E)

  gsa = derivative(K.pol)(gen(K))
  gsa_con = conjugates_arb(gsa, 32)

  if length(root_bound) > 0
    @assert length(root_bound) == r1 + r2
    bound_root = Vector{arb}(undef, r1 + r2)
    for i in 1:(r1 + r2)
      bound_root[i] = root_bound[i] * abs(gsa_con[i])
    end
  elseif ispure
    conj = __conjugates_arb(coeff(f, 0), 32)
    R = parent(conj[1])
    a_con = [R(abs_upper_bound(abs(e), fmpz)) for e in conj]

    bound_root = Vector{arb}(undef, r1 + r2)

    for i in 1:r1+r2
      bound_root[i] = root(abs(a_con[i]), degree(f)) * abs(gsa_con[i])
    end
  else
    bound_root = Vector{arb}(undef, r1 + r2)

    CC = AcbField(64, cached = false)
    CCt, t = PolynomialRing(CC, "t", cached=false)
    conjugates_of_coeffs = Vector{acb}[ conjugates_arb(c, 32) for c in coeffs_f ]

    for i in 1:r1+r2
      g = CCt(acb[ conjugates_of_coeffs[j + 1][i] for j in 0:degree(f) ])
      bound_root[i] = roots_upper_bound(g) * abs(gsa_con[i])
    end
  end
  ss = _lifting_expo(good_p, good_deg_p, E, bound_root)
  rts = _hensel(f, factor_of_g, good_fp, Int(ss), max_roots = max_roots - length(rs), ispure = ispure)
  return vcat(rs, rts)
end

################################################################################
#
#  Main lifting functions
#
################################################################################

#this is (or should be) the HNF basis for P^??
function _get_basis(pp::fmpz, n::Int, pgg::fmpz_mod_poly, Qt::FmpzModPolyRing)
  M = zero_matrix(FlintZZ, n, n)
  Q = base_ring(Qt)
  #the lattice for reco:
  #zero!(M)
  for j = 1:degree(pgg)
    M[j, j] = pp
  end
  coeffarr = Vector{elem_type(Q)}(undef, degree(pgg))
  for j = 1:degree(pgg)-1
    coeffarr[j] = zero(Q)
  end
  coeffarr[degree(pgg)] = one(Q)
  pt = Qt(coeffarr)
  for j=degree(pgg)+1:n
    pt = shift_left(pt, 1)
    rem!(pt, pt, pgg)
    M[j,j] = 1
    for k=0:degree(pt)
      M[j, k+1] = -lift(coeff(pt, k))
    end
  end
  return M
end

function _get_LLL_basis(Mold, Miold, dold, p, pr, i, gg)
  n = nrows(Mold)
  ctx = lll_ctx(0.5, 0.51)
  @show pr[i-1], pr[i]
  modu = fmpz(p)^25
  for j = (pr[i-1]+25):25:pr[i]
    pp = fmpz(p)^j
    Q = ResidueRing(FlintZZ, pp, cached=false)
    Qt, t = PolynomialRing(Q, "t", cached=false)
    pgg = Qt(gg)
    M = _get_basis(pp, n, pgg, Qt)
    mul!(M, M, Miold)
    divexact!(M, M, dold)
    @vtime :Saturate 1 hnf_modular_eldiv!(M, modu)
    @vtime :Saturate 1 lll!(M, ctx)
    @vtime :Saturate 1 lll!(M)
    mul!(M, M, Mold)
    @vtime :Saturate 1 lll!(M, ctx)
    @vtime :Saturate 1 lll!(M)
    Mi, d = pseudo_inv(M)
    Mold = M
    Miold = Mi
    dold = d
  end
  if !iszero(mod(pr[i]-pr[i-1], 25))
    modu = fmpz(p)^mod(pr[i]-pr[i-1], 25)
    pp = fmpz(p)^pr[i]
    Q = ResidueRing(FlintZZ, pp, cached=false)
    Qt, t = PolynomialRing(Q, "t", cached=false)
    pgg = Qt(gg)
    M = _get_basis(pp, n, pgg, Qt)
    mul!(M, M, Miold)
    divexact!(M, M, dold)
    @vtime :Saturate 1 hnf_modular_eldiv!(M, modu)
    @vtime :Saturate 1 lll!(M, ctx)
    @vtime :Saturate 1 lll!(M)
    mul!(M, M, Mold)
    @vtime :Saturate 1 lll!(M, ctx)
    @vtime :Saturate 1 lll!(M)
    Mi, d = pseudo_inv(M)
    Mold = M
    Miold = Mi
    dold = d
  end
  return Mold, Miold, dold
end

function _hensel(f::Generic.Poly{nf_elem},
                 fac_pol_mod_p::gfp_poly,
                 fp::fq_nmod_poly, k::Int;
                 max_roots::Int = degree(f),
                 ispure::Bool = false,
                 isnormal::Bool = false)
  # Let f be a polynomial over a number field K with defining polynomial g.
  # This function expects
  # fac_pol_mod_p = factor of g mod p corresponding to a prime ideal P of K
  # fp = f mod P, the reduction of f modulo P
  # k = lifting exponent
  # This function lifts the roots of f mod P to P^k and reconstructs them.
  # f is pure if and only if f = x^deg(f) + coeff(f, 0)
  @assert max_roots > 0

  caching = degree(base_ring(f)) > 20
  if caching
    # Setup the caching
    _cache = _get_prime_data_lifting(base_ring(f))
    _p = Int(characteristic(base_ring(fac_pol_mod_p)))
    if haskey(_cache, _p)
      _cache_p = _cache[_p]
      @vprint :Saturate 4 "Hitting the cache for the prime $(_p)\n"
      fac_pol_mod_p.parent = parent(first(keys(_cache_p)))
      if haskey(_cache_p, fac_pol_mod_p)
        @vprint :Saturate 1 "  Hitting the cache for the prime ideal\n"
        _cache_lll = _cache_p[fac_pol_mod_p]
      else
        @vprint :Saturate 1 "  Not hitting the cache for the prime ideal\n"
        _cache_lll = Dict()
        _cache_p[fac_pol_mod_p] = _cache_lll
      end
    else
      @vprint :Saturate 4 "No hitting the cache for the prime $(_p)\n"
      _cache_p = Dict()
      _cache_lll = Dict()
      _cache_p[fac_pol_mod_p] = _cache_lll
      _cache[_p] = _cache_p
    end
  end

  k = max(k, 2)

  # We want normalized PR's
  if caching
    new_k = Int(round(fmpz(k)//10^flog(fmpz(k), 10))) * 10^flog(fmpz(k), 10)
    new_k < k ? new_k = (1 + Int(round(fmpz(new_k)//10^flog(fmpz(new_k), 10)))) * 10^flog(fmpz(new_k), 10) : new_k = new_k
    k = new_k
  end

  Rp = base_ring(fac_pol_mod_p)
  Rpt = parent(fac_pol_mod_p)
  t = gen(Rpt)
  p = Int(characteristic(Rp))

  #assumes f squarefree
  #assumes constant_coefficient(f) != 0

  ZX, X = PolynomialRing(FlintZZ, "X", cached = false)
  K = base_ring(f)

  #to avoid embarrasment...

  #find the prime ideal - as I don't want to use orders, this is
  #fun (computing a max order just for this is wasteful)
  #fun fact: if g = prod g_i mod p^k, then P_i^k = <p^k, g_i>
  #so instead of powering, and simplify and such, lets write it down
  if degree(fac_pol_mod_p) != degree(K)
    g1 = lift(ZX, fac_pol_mod_p)
    gg = hensel_lift(ZX(K.pol), g1, fmpz(p), k)
  else
    gg = ZX(K.pol)
  end
  # now for all i<= k, <p^i, K(gg)+p^i> is a normal presentation for
  #                                     the prime ideal power.
  #later we'll get the HNF matrix for selected powers as well

  #set up the mod p data:
  #need FiniteField as I need to factor (roots)

  rt = roots(fp)

  # Do something else in the pure case
  if ispure
    rt = eltype(rt)[inv(x^(degree(f)-1)) for x = rt]
  end

  # Now if we are in the normal case and want max_roots, we only have
  # to lift max_roots

  if isnormal
    rt = eltype(rt)[1:max(max_roots, degree(f))]
  end

  #we're going to lift the roots - and for efficiency 1/f'(rt) as well:
  fps = derivative(fp)
  irt = eltype(rt)[inv(fps(x)) for x in rt]

  # this is in the finite field, but now I need it in the res. ring.
  # in ZX for ease of transport.
  RT = fmpz_poly[lift(ZX, lift(Rpt, x)) for x = rt]
  IRT = fmpz_poly[lift(ZX, lift(Rpt, x)) for x = irt]

  #the den: ala Kronnecker:
  den = ZX(derivative(K.pol)(gen(K)))
  iden = inv(derivative(K.pol)(gen(K)))

  ##for the result, to check for stabilising

  res = nf_elem[zero(K) for x = RT]
  rs = nf_elem[]

  #we'll be working in (Z/p^k)[t]/gg

  #an "optimal" lifting chain:
  pr = Int[k]
  while k>1
    k = div(k+1, 2)
    push!(pr, k)
  end
  pr = reverse(pr)

  ##lets start:
  f_coeff_ZX = Vector{fmpz_poly}(undef, length(f))
  if !ispure
    for j in 0:degree(f)
      f_coeff_ZX[j + 1] = ZX(coeff(f, j))
    end
  end

  # Keep a bit array to keep track of what we have to still lift
  roots_to_lift = trues(length(rt))

  n = degree(K)
  M = zero_matrix(FlintZZ, n, n)
  local Mi::fmpz_mat
  local d::fmpz
  @vprint :Saturate 1 "Maximum number of steps: $(length(pr))\n"
  for i=2:length(pr)
    @vprint :Saturate 1 "Step number $i\n"
    pp = fmpz(p)^pr[i]
    Q = ResidueRing(FlintZZ, pp, cached=false)
    Qt, t = PolynomialRing(Q, "t", cached=false)

    #possibly this should be done with max precision and then adjusted down
    #the poly mod P^??
    if !ispure
      fpp = fmpz_mod_poly[Qt(f_coeff_ZX[k + 1]) for k=0:degree(f)]
    end

    #we need to evaluate fp and fp' at the roots (later)
    #given that we evaluate "by hand" we don't need polynomials

    pgg = Qt(gg) #we'll do the reductions by hand - possibly not optimal

    ctx_lll = lll_ctx(0.3, 0.51)
    if caching && haskey(_cache_lll, pr[i])
      M, Mi, d = _cache_lll[pr[i]]::Tuple{fmpz_mat, fmpz_mat, fmpz}
    elseif i > 10
      #This is getting bad. We try to apply the trick twice.
      Mold = M
      Miold = Mi
      dold = d
      pr_intermediate = pr[i-1] + div(pr[i] - pr[i-1], 2)
      ppint = fmpz(p)^pr_intermediate
      Qint = ResidueRing(FlintZZ, ppint, cached = false)
      Qintt = PolynomialRing(Qint, "t", cached = false)[1]
      pggQint = Qintt(gg)
      Mint = _get_basis(ppint, n, pggQint, Qintt)
      mul!(Mint, Mint, Miold)
      divexact!(Mint, Mint, dold)
      exp_mod = fmpz(p)^(pr_intermediate - pr[i-1])
      @vtime :Saturate 1 hnf_modular_eldiv!(Mint, exp_mod)
      @vtime :Saturate 1 lll!(Mint, ctx_lll)
      @vtime :Saturate 1 lll!(Mint)
      mul!(Mint, Mint, Mold)
      @vtime :Saturate 1 lll!(Mint, lll_ctx(0.6, 0.51))
      @vtime :Saturate 1 lll!(Mint)
      Miold, dold = pseudo_inv(Mint)
      Mold = Mint
      M = _get_basis(pp, n, pgg, Qt)
      mul!(M, M, Miold)
      divexact!(M, M, dold)
      exp_mod = divexact(pp, ppint)
      @vtime :Saturate 1 hnf_modular_eldiv!(M, exp_mod)
      @vtime :Saturate 1 lll!(M, ctx_lll)
      @vtime :Saturate 1 lll!(M)
      mul!(M, M, Mold)
      @vtime :Saturate 1 lll!(M, lll_ctx(0.6, 0.51))
      @vtime :Saturate 1 lll!(M)
      Mi, d = pseudo_inv(M)
      if caching
        _cache_lll[pr[i]] = (M, Mi, d)
      end
    elseif i > 3
      Mold = M
      Miold = Mi
      dold = d
      M = _get_basis(pp, n, pgg, Qt)
      mul!(M, M, Miold)
      divexact!(M, M, dold)
      exp_mod = div(pr[i], 2)
      if isodd(pr[i])
        exp_mod += 1
      end
      @vtime :Saturate 1 hnf_modular_eldiv!(M, fmpz(p)^exp_mod)
      @vtime :Saturate 1 lll!(M, ctx_lll)
      @vtime :Saturate 1 lll!(M)
      mul!(M, M, Mold)
      @vtime :Saturate 1 lll!(M, lll_ctx(0.6, 0.51))
      @vtime :Saturate 1 lll!(M)
      Mi, d = pseudo_inv(M)
      if caching
        _cache_lll[pr[i]] = (M, Mi, d)
      end
    else
      M = _get_basis(pp, n, pgg, Qt)
      @vtime :Saturate 1  lll!(M)
      Mi, d = pseudo_inv(M)
      if caching
        _cache_lll[pr[i]] = (M, Mi, d)
      end
    end

    if ispure
      ap = Qt((-coeff(f, 0)))
      bp = powmod(ap, degree(f)-1, pgg)
      minv = invmod(fmpz(degree(f)), pp)
    end

    for j=1:length(RT)
      if !roots_to_lift[j]
        continue
      end
      #to eval fp and the derivative, we pre-compute the powers of the
      #evaluation point - to save on large products...

      if !ispure
        pow = fmpz_mod_poly[Qt(1), Qt(RT[j])]
        while length(pow) <= degree(f)+1
          push!(pow, pow[2]*pow[end] % pgg)
        end

        eval_f = pow[1] * fpp[1]
        for k in 2:length(fpp)
          eval_f = eval_f +  pow[k] * fpp[k]
        end
        eval_f = eval_f % pgg

        #eval_f = sum(pow[k] * fp[k] for k=1:length(fp)) % pgg
        eval_fs = pow[1] * fpp[2]
        for k in 3:length(fpp)
          eval_fs = eval_fs + pow[k-1] *(k-1)*fpp[k]
        end
        eval_fs = eval_fs % pgg
        #eval_fs = sum(pow[k-1] *(k-1)*fp[k] for k=2:length(fp)) % pgg

        #double lift:
        #IRT = invmod(fp'(rt), p^k)
        # using x -> x(2-xy) to compute the inverse of y
        IRT[j] = lift(ZX, Qt(IRT[j])*(Qt(2-IRT[j]*eval_fs) % pgg) %pgg)
        #RT = rt mod p^k normal Newton
        # using x -> x-fp(x)//fp'(x) = x-fp(x) * IRT
        RT[j] = lift(ZX, Qt(pow[2] - eval_f*Qt(IRT[j])) % pgg)

        #before the reconstruction, we need to scale by den
        cf = lift(ZX, Qt(RT[j]*den) % pgg)
      else
        RTjp = Qt(RT[j])
        RT[j] = lift(ZX, (RTjp*(1+minv) - bp*minv* powmod(RTjp, degree(f)+1, pgg)) % pgg)

        #before the reconstruction, we need to scale by den and by a
        cf = lift(ZX, (Qt(RT[j]*den) % pgg)*ap % pgg)
      end

      ve = matrix(FlintZZ, 1, n, [coeff(cf, k) for k=0:n-1])
      _ve = ve*Mi
      mu = matrix(FlintZZ, 1, n,  [ round(fmpz, _ve[1, k]//d) for k=1:n])
      ve = ve - mu*M
      z = ZX()
      for kk=1:n
        setcoeff!(z, kk-1, ve[1, kk])
      end
      zz = K(z)*iden
      if res[j] == zz || i == length(pr)
        if (!ispure && iszero(f(zz))) || (ispure && zz^degree(f) == -coeff(f, 0))
          push!(rs, zz)
          if length(rs) >= max_roots
            return rs
          end
          roots_to_lift[j] = false
        end
        # I lifted to the end and one root did not lift to a root of f
        # and f is normal. Then there are no roots.
        if i == length(pr) && isnormal
          return nf_elem[]
        end
      else
        res[j] = zz
      end
    end
  end

  return rs
end

function _hensel(f::Generic.Poly{nf_elem}, p::Int, k::Int; max_roots::Int = degree(f))
  # Let f be a polynomial over a number field K
  # This function expects
  # p = prime number
  # k = lifting exponent
  # This function lifts the roots of f mod P to P^k and reconstructs them,
  # where P is a prime ideal with smallest degree over p

  K = base_ring(f)
  k = max(k, 2)
  Rp = GF(p, cached=false)
  Rpt, t = PolynomialRing(Rp, "t", cached=false)
  gp = Rpt(K.pol)
  lp = factor(gp).fac
  lpfac = first(keys(lp))

  for lpfac in keys(lp)
    if issquarefree(lpfac)
      break
    end
  end

  S = FqNmodFiniteField(lpfac, :z, false)
  ST, T = PolynomialRing(S,"T", cached=false)
  fp = ST([S(Rpt(coeff(f, i))) for i=0:degree(f)])

  return _hensel(f, lpfac, fp, k, max_roots = max_roots)
end

################################################################################
#
#  Computation of the lifting exponent
#
################################################################################

function _lifting_expo(p::Int, deg_p::Int, O::NfOrd, bnd::Array{arb, 1})
  # return _lifting_expo_using_logbound(p, deg_p, O, arb[log(a) for a in bnd])
  # compute the lifting exponent a la Friedrich-Fieker
  # bnd has upper bounds on |x^{(i)}| 1<= i <= r1+r2 as arbs
  # we're using a prime ideal above p of intertia degree deg_p
  # O is the order where the result will be reconstructed in

  (c1, c2) = norm_change_const(O)
  r1, r2 = signature(O)
  R = parent(bnd[1])
  bd = zero(R)
  n = degree(O)
  #so   |x|_mink  <= c_1 |x|_coeff
  #and  |x|_coeff <= c_2 |x|_mink

  for i in 1:r1
    bd += bnd[i]^2
  end

  for i=1:r2
    bd += 2*bnd[i+r1]^2
  end

  boundt2 = max(bd, one(R))

  # Tommy: log(...) could contain a ball, which contains zero
  tmp = R(abs_upper_bound(R(c1)*R(c2)*boundt2*exp((R(n*(n-1))//2 + 2)*log(R(2)))//n, fmpz))

  # CF: there is a prob, in the paper wrt LLL bounds on |x| or |x|^2....
  boundk = R(n)*log(tmp)//(2*deg_p*log(R(p)))

  ss = abs_upper_bound(boundk, fmpz)
  return ss
end

#identical to hasroot - which one do we want?
function ispower(a::NfOrdElem, n::Int)
  fl, r = ispower(nf(parent(a))(a), n, isintegral = true)
  return fl, parent(a)(r)
end


text: include("direct.jl")


text:
module Sync
# synchronous task runner

import ..ArbiterBase: task_loop, TaskResult

export run_tasks

"""
Run an iterable of tasks.

tasks: The iterable of tasks
"""
function run_tasks(tasks)
    return task_loop(tasks, execute)
end

"""
Execute a task, returning a Nullable{TaskResult}
"""
function execute(task)
    return Nullable(TaskResult(task.name, task.func()))
end

end


text:
module TestIterators
using Test

using Phylo
using Diversity

@testset "Iterators" begin
    species = 10
    sc = 5
    abund = rand(species, sc)
    abund ./= sum(abund)
    ru = rand(Ultrametric(species))
    m = Metacommunity(abund, PhyloBranches(ru))
    @test length(TypeIterator(m)) == counttypes(PhyloBranches(ru))
    @test length(SubcommunityIterator(m)) == sc
    @test all(reduce(+, SubcommunityIterator(m)) .≈ getmetaabundance(m))
    @test all(reduce(+, TypeIterator(m)) .≈ getweight(m))
end

end