text:
export PRIMITIVES,
DARTSModel,
Cell,
MixedOp,
DARTSEvalModel,
EvalCell,
Activations,
discretize,
DARTSEvalAuxModel,
DARTSModelBN,
parse_genotype

using Flux
using Base.Iterators
using StatsBase: mean
using Zygote
using LinearAlgebra
using CUDA
using Zygote: @adjoint, dropgrad
using Distributions

ReLUConvBN(channels_in, channels_out, kernel_size, stride, pad) =
    Chain(
        x -> relu.(x),
        Conv(kernel_size, channels_in => channels_out, bias = false),
        BatchNorm(channels_out),
    )

function FactorizedReduce(channels_in, channels_out, stride)
    odd =
        Conv(
            (1, 1),
            channels_in => channels_out ÷ stride,
            stride = (stride, stride),
            bias = false,
        ) |> gpu
    even =
        Conv(
            (1, 1),
            channels_in => channels_out ÷ stride,
            stride = (stride, stride),
            bias = false,
        ) |> gpu

    Chain(
        x -> relu.(x),
        x -> cat(odd(x), even(x[2:end, 2:end, :, :]), dims = 3), #TODO get rid of indexing
        BatchNorm(channels_out),
    )
end

SepConv(channels_in, channels_out, kernel_size, stride, pad) =
    Chain(
        x -> relu.(x),
        DepthwiseConv(
            kernel_size,
            channels_in => channels_in,
            stride = (stride, stride),
            pad = (pad, pad),
            bias = false,
        ),
        Conv(
            (1, 1),
            channels_in => channels_in,
            stride = (1, 1),
            pad = (0, 0),
            bias = false,
        ),
        BatchNorm(channels_in),
        x -> relu.(x),
        DepthwiseConv(
            kernel_size,
            channels_in => channels_in,
            pad = (pad, pad),
            stride = (1, 1),
            bias = false,
        ),
        Conv(
            (1, 1),
            channels_in => channels_out,
            stride = (1, 1),
            pad = (0, 0),
            bias = false,
        ),
        BatchNorm(channels_out),
    )

DilConv(channels_in, channels_out, kernel_size, stride, pad, dilation) =
    Chain(
        x -> relu.(x),
        DepthwiseConv(
            kernel_size,
            channels_in => channels_in,
            pad = (pad * dilation, pad * dilation),
            stride = (stride, stride),
            dilation = dilation,
            bias = false,
        ),
        Conv(
            (1,1),
            channels_in => channels_out,
            stride = (1, 1),
            pad = (0, 0),
            bias = false,
        ),
        BatchNorm(channels_out),
    )

SepConv_v(channels_in, channels_out, kernel_size, stride, pad) =
    Chain(
        x -> relu.(x),
        Conv(
            kernel_size,
            channels_in => channels_in,
            stride = (stride, stride),
            pad = (pad, pad),
            bias = false,
        ),
        Conv(
            (1, 1),
            channels_in => channels_in,
            stride = (1, 1),
            pad = (0, 0),
            bias = false,
        ),
        BatchNorm(channels_in),
        x -> relu.(x),
        Conv(
            kernel_size,
            channels_in => channels_in,
            pad = (pad, pad),
            stride = (1, 1),
            bias = false,
        ),
        Conv(
            (1, 1),
            channels_in => channels_out,
            stride = (1, 1),
            pad = (0, 0),
            bias = false,
        ),
        BatchNorm(channels_out),
    )

DilConv_v(channels_in, channels_out, kernel_size, stride, pad, dilation) =
    Chain(
        x -> relu.(x),
        Conv(
            kernel_size,
            channels_in => channels_in,
            pad = (pad * dilation, pad * dilation),
            stride = (stride, stride),
            dilation = dilation,
            bias = false,
        ),
        Conv(
            (1,1),
            channels_in => channels_out,
            stride = (1, 1),
            pad = (0, 0),
            bias = false,
        ),
        BatchNorm(channels_out),
    )

Identity() = x -> x
Zero(stride) =
    stride == 2 ? Chain(MeanPool((1, 1), stride = stride), x -> 0*x) :
    x -> 0*x

SkipConnect(channels_in, channels_out, stride, pad) =
    stride == 1 ? Identity() :
    FactorizedReduce(channels_in, channels_out, stride)

struct AdaptiveMeanPool{N}
    target_out::NTuple{N,Int}
end

AdaptiveMeanPool(target_out::NTuple{N,Integer}) where {N} =
    AdaptiveMeanPool(target_out)

function (m::AdaptiveMeanPool)(x)
    w = size(x, 1) - m.target_out[1] + 1
    h = size(x, 2) - m.target_out[2] + 1
    return meanpool(x, (w, h); pad = (0, 0), stride = (1, 1))
end

PRIMITIVES = [
    "none",
    "max_pool_3x3",
    "avg_pool_3x3",
    "skip_connect",
    "sep_conv_3x3",
    "sep_conv_5x5",
    #"sep_conv_7x7",
    "dil_conv_3x3",
    "dil_conv_5x5",
    #"conv_7x1_1x7"
]

OPS = Dict(
    "none" => (channels, stride, w) -> Chain(Zero(stride)),
    "avg_pool_3x3" =>
        (channels, stride, w) ->
            Chain(
                MeanPool((3, 3), stride = stride, pad = 1),
                BatchNorm(channels),
            ),
    "max_pool_3x3" =>
        (channels, stride, w) ->
            Chain(
                MaxPool((3, 3), stride = stride, pad = 1),
                BatchNorm(channels),
            ),
    "skip_connect" =>
        (channels, stride, w) ->
            Chain(SkipConnect(channels, channels, stride, 1)),
    "sep_conv_3x3" =>
        (channels, stride, w) ->
            SepConv(channels, channels, (3, 3), stride, 1),
    "sep_conv_5x5" =>
        (channels, stride, w) ->
            SepConv(channels, channels, (5, 5), stride, 2),
    "sep_conv_7x7" =>
        (channels, stride, w) ->
            SepConv(channels, channels, (7, 7), stride, 3),
    "dil_conv_3x3" =>
        (channels, stride, w) ->
            DilConv(channels, channels, (3, 3), stride, 1, 2),
    "dil_conv_5x5" =>
        (channels, stride, w) ->
            DilConv(channels, channels, (5, 5), stride, 2, 2),
    "conv_7x1_1x7" =>
        (channels, stride, w) ->
            Chain(
                x -> relu.(x),
                Conv(
                    (1, 7),
                    channels => channels,
                    pad = (0, 3),
                    stride = (1, stride),
                    bias = false,
                ),
                Conv(
                    (7, 1),
                    channels => channels,
                    pad = (3, 0),
                    stride = (stride, 1),
                    bias = false,
                ),
                BatchNorm(channels),
            ),
)

function my_softmax(xs::AbstractArray; dims = 1)
    softmax!(similar(xs), xs, dims=dims)
end

Zygote.@adjoint function my_softmax(xs::AbstractArray; dims = 1)
    softmax(xs, dims = dims), Δ -> begin
        (∇softmax(Δ, xs, dims = dims),)
    end
end

abstract type AbstractOp end

struct OpActs <: AbstractOp
    name::String
    op::Any
end

function showlayer(
    x::AbstractArray,
    layer,
    opname::String,
    outs::Array{Array{Float32,3}},
)
    out = layer(x)
    Zygote.ignore() do
        if !(typeof(layer) <: Flux.BatchNorm) && !occursin("none", opname)
            push!(outs, dropdims(mean(out, dims = (1, 2, 3)), dims = 3) |> cpu)
        end
    end
    out
end

function (opwrap::OpActs)(xin::AbstractArray, acts::Dict, cellid::Int64)
    outs = Array{Array{Float32,3}}(undef, 0)
    xout = foldl(
        (x, layer) -> showlayer(x, layer, opwrap.name, outs),
        opwrap.op,
        init = xin,
    )
    Zygote.ignore() do
        if length(outs) > 0
            if haskey(acts, opwrap.name)
                acts[opwrap.name] = vcat(acts[opwrap.name], hcat(outs...))
            else
                acts[opwrap.name] = hcat(outs...)
            end
        end #acts dims: cellid x oplayer x image
    end
    xout
end

Flux.@functor OpActs

struct Op <: AbstractOp
    name::String
    op::Chain
end

function (opwrap::Op)(
    xin::AbstractArray,
    acts::Union{Dict,Nothing},
    cellid::Int64 = 0,
)
    opwrap.op(xin)
end

Flux.@functor Op

struct MixedOp{T<:AbstractOp}
    cellid::Int64
    name::String
    ops::Vector{T}
end

function MixedOp(
    cellid::Int64,
    name::String,
    channels::Int64,
    stride::Int64,
    operations::Vector{String} = PRIMITIVES,
    track_acts::Bool = false,
)
    if track_acts
        op = OpActs
    else
        op = Op
    end
    MixedOp(
        cellid,
        name,
        [
            op(string(name, "-", operation), OPS[operation](channels, stride, 1) |> gpu)
            for operation in operations
        ],
    )
end

function (m::MixedOp)(
    x::CuArray{Float32, 4},
    #out::CuArray{Float32, 4},
    as::Vector{Float32} = [1f0],
    acts::Union{Nothing,Dict} = nothing,
)
    as = softmax(as)
    sum(as[i] * m.ops[i](x, acts, m.cellid) for i = 1:length(as))
end

"""
function (m::MixedOp)(
    x::CuArray{Float32, 4},
    precompute::AbstractArray,
    αs::Vector{Float32} = [1f0],
)
    as = my_softmax(copy(αs))
    if !precompute[1]
        precompute[2] = [m.ops[i](x, nothing, m.cellid) for i = 1:length(as)]
        precompute[1] = true
    end
    sum(as[i] * precompute[2][i] for i = 1:length(as))
end
"""

Flux.@functor MixedOp

struct Cell
    steps::Int64
    reduction::Bool
    multiplier::Int64
    prelayer1::Chain
    prelayer2::Chain
    mixedops::Vector{MixedOp}
end

function Cell(
    channels_before_last::Int64,
    channels_last::Int64,
    channels::Int64,
    reduce::Bool,
    reduce_previous::Bool,
    steps::Int64,
    multiplier::Int64,
    cellid::Int64,
    operations::Vector{String} = PRIMITIVES,
    track_acts::Bool = false,
)
    if reduce_previous
        prelayer1 = FactorizedReduce(channels_before_last, channels, 2) |> gpu
    else
        prelayer1 =
            ReLUConvBN(channels_before_last, channels, (1, 1), 1, 0) |> gpu
    end
    prelayer2 = ReLUConvBN(channels_last, channels, (1, 1), 1, 0) |> gpu
    mixedops = []
    for i = 3:steps+2 #op output
        for j = 1:i-1 #op input
            reduce && j < 3 ? stride = 2 : stride = 1
            mixedop =
                MixedOp(
                    cellid,
                    string(j, "-", i),
                    channels,
                    stride,
                    operations,
                    track_acts,
                )
            push!(mixedops, mixedop)
        end
    end
    Cell(steps, reduce, multiplier, prelayer1, prelayer2, mixedops)
end

function (m::Cell)(
    x1::CuArray{Float32, 4},
    x2::CuArray{Float32, 4},
    αs::Vector{Vector{Float32}},
    acts::Union{Nothing,Dict} = nothing,
)
    state1 = m.prelayer1(x1)

    states = Zygote.Buffer([state1], m.steps + 2)

    states[1] = state1
    states[2] = m.prelayer2(x2)
    offset = 0
    #precompute = Array{Array{Float32,4}}(undef, m.steps + 1)
    #precompute = Zygote.Buffer([], m.steps + 1)
    #for i = 1:m.steps+1
        #precompute[i] = [false, []]
    #end
    #@display(typeof(precompute))
    for step = 1:m.steps

        state = mapreduce(
            (mixedop, α, previous_state) ->
                mixedop(previous_state, α, acts),
            +,
            m.mixedops[offset+1:offset+step+1],
            αs[offset+1:offset+step+1],
            states,
        )

        offset += step + 1
        states[step+2] = state
    end
    reduce(ccat, copy(states[m.steps+2-m.multiplier+1:m.steps+2]))
end

Flux.@functor Cell

mutable struct Activations
    currentacts::Dict{String,Array{Float32,3}}
end

struct DARTSModel
    normal_αs::Vector{Vector{Float32}}
    reduce_αs::Vector{Vector{Float32}}
    stem::Chain
    cells::Array{Cell,1}
    global_pooling::AdaptiveMeanPool
    classifier::Dense
    operations::Vector{String}
    activations::Activations
end

function DARTSModel(;
    α_init = (num_ops -> 2f-3 * (rand(Float32, num_ops) .- 0.5f0)),
    num_classes::Int64 = 10,
    num_cells::Int64 = 8,
    channels::Int64 = 16,
    steps::Int64 = 4,
    mult::Int64 = 4,
    stem_mult::Int64 = 3,
    operations::Vector{String} = PRIMITIVES,
    track_acts::Bool = false,
)
    channels_current = channels * stem_mult
    stem =
        Chain(
            Conv((3, 3), 3 => channels_current, pad = (1, 1), bias = false),
            BatchNorm(channels_current),
        ) |> gpu

    channels_before_last = channels_current
    channels_last = channels_current
    channels_current = channels
    reduce_previous = false
    cells = []
    for i = 1:num_cells
        if i == num_cells ÷ 3 + 1 || i == 2 * num_cells ÷ 3 + 1
            channels_current = channels_current * 2
            reduce = true
        else
            reduce = false
        end
        cell =
            Cell(
                channels_before_last,
                channels_last,
                channels_current,
                reduce,
                reduce_previous,
                steps,
                mult,
                i,
                operations,
                track_acts
            )
        push!(cells, cell)

        reduce_previous = reduce
        channels_before_last = channels_last
        channels_last = mult * channels_current
    end
    global_pooling = AdaptiveMeanPool((1, 1)) |> gpu
    classifier = Dense(channels_last, num_classes) |> gpu
    k = floor(Int, steps^2 / 2 + 3 * steps / 2)
    num_ops = length(operations)
    α_normal = [α_init(num_ops) for _ = 1:k]
    α_reduce = [α_init(num_ops) for _ = 1:k]
    activations = Activations(Dict{String,Array{Float32,3}}())
    DARTSModel(
        α_normal,
        α_reduce,
        stem,
        cells,
        global_pooling,
        classifier,
        operations,
        activations
    )
end

function (m::DARTSModel)(x::CuArray{Float32, 4})
    normal_αs = m.normal_αs
    reduce_αs = m.reduce_αs
    s1 = m.stem(x)
    s2 = m.stem(x)
    m.activations.currentacts = Dict{String,Array{Float32,3}}()
    for (i, cell) in enumerate(m.cells)
        cell.reduction ? αs = reduce_αs : αs = normal_αs
        new_state = cell(s1, s2, αs, m.activations.currentacts)
        s1 = s2
        s2 = new_state
    end
    out = m.global_pooling(s2)
    m.classifier(squeeze(out))
end

function (m::DARTSModel)(x::CuArray{Float32, 4}, αs::Vector{Vector{Vector{Float32}}})
    normal_αs = αs[1]
    reduce_αs = αs[2]
    s1 = m.stem(x)
    s2 = m.stem(x)
    m.activations.currentacts = Dict{String,Array{Float32,3}}()
    for (i, cell) in enumerate(m.cells)
        cell.reduction ? αs = reduce_αs : αs = normal_αs
        new_state = cell(s1, s2, αs, m.activations.currentacts)
        s1 = s2
        s2 = new_state
    end
    out = m.global_pooling(s2)
    m.classifier(squeeze(out))
end

Flux.@functor DARTSModel

struct MixedOpBN{T<:AbstractOp}
    cellid::Int64
    name::String
    ops::Vector{T}
    batchnorm::BatchNorm
end

function MixedOpBN(
    cellid::Int64,
    name::String,
    channels::Int64,
    stride::Int64,
    operations::Vector{String} = PRIMITIVES[2:end],
    track_acts::Bool = false,
)
    if track_acts
        op = OpActs
    else
        op = Op
    end
    MixedOpBN(
        cellid,
        name,
        [
            op(string(name, "-", prim), OPS[prim](channels, stride, 1) |> gpu)
            for prim in operations
        ],
        BatchNorm(channels) |> gpu
    )
end

function (m::MixedOpBN)(
    x::CuArray{Float32, 4},
    αs::Vector{Float32} = [1f0],
    acts::Union{Nothing,Dict} = nothing,
)
    if any(αs .> 0)
        m.batchnorm(sum(tanh(αs[i]) * m.ops[i](x, acts, m.cellid) for i = 1:length(αs) if αs[i] > 0))
    else
        0.0 * m.ops[1](x, acts, m.cellid)
    end
end

Flux.@functor MixedOpBN

struct CellBN
    steps::Int64
    reduction::Bool
    multiplier::Int64
    prelayer1::Chain
    prelayer2::Chain
    mixedops::Vector{MixedOpBN}
end

function CellBN(
    channels_before_last::Int64,
    channels_last::Int64,
    channels::Int64,
    reduce::Bool,
    reduce_previous::Bool,
    steps::Int64,
    multiplier::Int64,
    cellid::Int64,
    operations::Vector{String} = PRIMITIVES[2:end],
    track_acts::Bool = false,
)
    if reduce_previous
        prelayer1 = FactorizedReduce(channels_before_last, channels, 2) |> gpu
    else
        prelayer1 =
            ReLUConvBN(channels_before_last, channels, (1, 1), 1, 0) |> gpu
    end
    prelayer2 = ReLUConvBN(channels_last, channels, (1, 1), 1, 0) |> gpu
    mixedops = Vector{MixedOpBN}()
    for i = 3:steps+2 #op output
        for j = 1:i-1 #op input
            reduce && j < 3 ? stride = 2 : stride = 1
            mixedop =
                MixedOpBN(
                    cellid,
                    string(j, "-", i),
                    channels,
                    stride,
                    operations,
                    track_acts,
                )
            push!(mixedops, mixedop)
        end
    end
    CellBN(steps, reduce, multiplier, prelayer1, prelayer2, mixedops)
end

function (m::CellBN)(
    x1::CuArray{Float32, 4},
    x2::CuArray{Float32, 4},
    αs::Vector{Vector{Float32}},
    acts::Union{Nothing,Dict} = nothing,
)
    state1 = m.prelayer1(x1)

    states = Zygote.Buffer([state1], m.steps + 2)

    states[1] = state1
    states[2] = m.prelayer2(x2)
    offset = 0
    for step = 1:m.steps
        state = mapreduce(
            (mixedop, α, previous_state) ->
                mixedop(previous_state, α, acts),
            +,
            m.mixedops[offset+1:offset+step+1],
            αs[offset+1:offset+step+1],
            states,
        )
        offset += step + 1
        states[step+2] = state
    end
    reduce(ccat, copy(states[m.steps+2-m.multiplier+1:m.steps+2]))
end

Flux.@functor CellBN

struct DARTSModelBN
    normal_αs::Vector{Vector{Float32}}
    reduce_αs::Vector{Vector{Float32}}
    stem::Chain
    cells::Vector{CellBN}
    global_pooling::AdaptiveMeanPool
    classifier::Dense
    operations::Vector{String}
    activations::Activations
end

function DARTSModelBN(;
    α_init = (num_ops -> atanh.(ones(Float32, num_ops) ./ num_ops + 2f-3 * (rand(Float32, num_ops) .- 0.5f0))),
    num_classes::Int64 = 10,
    num_cells::Int64 = 8,
    channels::Int64 = 16,
    steps::Int64 = 4,
    mult::Int64 = 4,
    stem_mult::Int64 = 3,
    operations::Vector{String} = PRIMITIVES[2:end],
    track_acts::Bool = false,
)
    channels_current = channels * stem_mult
    stem =
        Chain(
            Conv((3, 3), 3 => channels_current, pad = (1, 1), bias = false),
            BatchNorm(channels_current),
        ) |> gpu
    channels_before_last = channels_current
    channels_last = channels_current
    channels_current = channels
    reduce_previous = false
    cells = Vector{CellBN}(undef, num_cells)
    for i = 1:num_cells
        if i == num_cells ÷ 3 + 1 || i == 2 * num_cells ÷ 3 + 1
            channels_current = channels_current * 2
            reduce = true
        else
            reduce = false
        end
        cells[i] =
            CellBN(
                channels_before_last,
                channels_last,
                channels_current,
                reduce,
                reduce_previous,
                steps,
                mult,
                i,
                operations,
                track_acts,
            )

        reduce_previous = reduce
        channels_before_last = channels_last
        channels_last = mult * channels_current
    end
    global_pooling = AdaptiveMeanPool((1, 1)) |> gpu
    classifier = Dense(channels_last, num_classes) |> gpu
    k = floor(Int, steps^2 / 2 + 3 * steps / 2)
    num_ops = length(operations)
    α_normal = [α_init(num_ops) for _ = 1:k]
    α_reduce = [α_init(num_ops) for _ = 1:k]
    activations = Activations(Dict{String,Array{Float32,3}}())
    DARTSModelBN(
        α_normal,
        α_reduce,
        stem,
        cells,
        global_pooling,
        classifier,
        operations,
        activations,
    )
end

function (m::DARTSModelBN)(x)
    s1 = m.stem(x)
    s2 = m.stem(x)
    m.activations.currentacts = Dict{String,Array{Float32,3}}()
    for (i, cell) in enumerate(m.cells)
        cell.reduction ? αs = m.reduce_αs : αs = m.normal_αs
        new_state = cell(s1, s2, αs, m.activations.currentacts)
        s1 = s2
        s2 = new_state
    end
    out = m.global_pooling(s2)
    m.classifier(squeeze(out))
end

Flux.@functor DARTSModelBN

struct EvalCell
    steps::Int64
    reduction::Bool
    multiplier::Int64
    prelayer1::Chain
    prelayer2::Chain
    ops::Vector{Vector{Chain}}
    inputindices::Vector{Vector{Int}}
    names::Vector{Tuple{String,String}}
end

function maxk(a::Vector{Vector{Float32}}, k::Int64)
    b = partialsortperm(a, 1:k, rev = true)
    return collect(b)
end

function discretize(
    αs::Vector{Vector{Float32}},
    channels::Int64,
    reduce::Bool,
    steps::Int64;
    operations::Vector{String} = PRIMITIVES,
    plotting::Bool = false,
    disclude_1::Bool = true,
)
    if operations == PRIMITIVES && size(αs[1],1) == 7
        operations = operations[2:end]
    end
    if "none" in operations
        disclude_1 = true
    else
        disclude_1 = false
    end
    ops = Vector{Vector{Chain}}(undef,steps)
    opnames = Vector{Tuple{String,String}}(undef,steps)
    inputindices = Vector{Vector{Int}}(undef,steps)
    rows = 0
    for i = 3:steps+2
        if disclude_1
            for j = 1:i-1
                αs[rows+j][1] = -Inf32
            end
        end
        options = [findmax(αs[rows+j]) for j = 1:i-1]
        top2 = partialsortperm(options, 1:2, by = x -> x[1], rev = true)
        top2names = Tuple(operations[options[i][2]] for i in top2)
        if !plotting
            top2ops = [
                OPS[operations[options[i][2]]](channels, reduce && i < 3 ? 2 : 1, 1) |> gpu
                for i in top2
            ]
            ops[i-2] = top2ops
        end
        inputindices[i-2] = top2
        opnames[i-2] = top2names
        rows += i - 1
    end
    (inputindices, ops, opnames)
end

function EvalCell(
    channels_before_last::Int64,
    channels_last::Int64,
    channels::Int64,
    reduce::Bool,
    reduce_previous::Bool,
    steps::Int64,
    multiplier::Int64,
    αs::Vector{Vector{Float32}},
    operations::Vector{String} = PRIMITIVES,
)
    if reduce_previous
        prelayer1 = FactorizedReduce(channels_before_last, channels, 2) |> gpu
    else
        prelayer1 = ReLUConvBN(channels_before_last, channels, (1, 1), 1, 0) |> gpu
    end
    prelayer2 = ReLUConvBN(channels_last, channels, (1, 1), 1, 0) |> gpu
    inputindices, ops, names = discretize(αs, channels, reduce, steps, operations)
    EvalCell(steps, reduce, multiplier, prelayer1, prelayer2, ops, inputindices, names)
end

DARTS_V2 = (
    normal = [
        ("sep_conv_3x3", 0),
        ("sep_conv_3x3", 1),
        ("sep_conv_3x3", 0),
        ("sep_conv_3x3", 1),
        ("sep_conv_3x3", 1),
        ("skip_connect", 0),
        ("skip_connect", 0),
        ("dil_conv_3x3", 2),
    ],
    reduce = [
        ("max_pool_3x3", 0),
        ("max_pool_3x3", 1),
        ("skip_connect", 2),
        ("max_pool_3x3", 1),
        ("max_pool_3x3", 0),
        ("skip_connect", 2),
        ("skip_connect", 2),
        ("max_pool_3x3", 1),
    ],
)

function parse_genotype(;
    genotype::NamedTuple=DARTS_V2,
    channels::Int64=1,
    reduce::Bool=false,
    steps::Int64=4,
)
    prim = [
        "none",
        "max_pool_3x3",
        "avg_pool_3x3",
        "skip_connect",
        "sep_conv_3x3",
        "sep_conv_5x5",
        #"sep_conv_7x7",
        "dil_conv_3x3",
        "dil_conv_5x5",
        #"conv_7x1_1x7"
    ]

    ops = Vector{Vector{Chain}}(undef,steps)
    opnames = Vector{Tuple{String,String}}(undef,steps)
    inputindices = Vector{Vector{Int}}(undef,steps)
    if reduce
        genocell = genotype.reduce
    else
        genocell = genotype.normal
    end
    for i = 1:steps
        top2 = [genocell[2*i-1][2]+1, genocell[2*i][2]+1]
        top2names = (genocell[2*i-1][1], genocell[2*i][1])
        top2ops = [
            OPS[top2names[j]](channels, reduce && top2[j] < 3 ? 2 : 1, 1) |> gpu
            for j in 1:2
        ]
        inputindices[i] = top2
        opnames[i] = top2names
        ops[i] = top2ops
    end
    (inputindices, ops, opnames)
end

function EvalCell(
    channels_before_last::Int64,
    channels_last::Int64,
    channels::Int64,
    reduce::Bool,
    reduce_previous::Bool,
    steps::Int64,
    multiplier::Int64,
)
    if reduce_previous
        prelayer1 = FactorizedReduce(channels_before_last, channels, 2) |> gpu
    else
        prelayer1 = ReLUConvBN(channels_before_last, channels, (1, 1), 1, 0)  |> gpu
    end
    prelayer2 = ReLUConvBN(channels_last, channels, (1, 1), 1, 0) |> gpu
    inputindices, ops, names = parse_genotype(genotype=DARTS_V2, channels=channels, reduce=reduce, steps=steps)
    display(collect(zip(inputindices,names)))
    EvalCell(steps, reduce, multiplier, prelayer1, prelayer2, ops, inputindices, names)
end

function droppath(x::CuArray{Float32, 4}, drop_prob::Float32)
    if drop_prob > 0.0
        mask = rand(Bernoulli(1 - drop_prob), 1, 1, size(x, 3), 1) |> gpu
        x = x .* mask / (1 - drop_prob)
    end
    x
end

function (m::EvalCell)(x1::CuArray{Float32, 4}, x2::CuArray{Float32, 4}, drop_prob::Float32)
    state1 = m.prelayer1(x1)

    states = Zygote.Buffer([state1], m.steps + 2)

    states[1] = state1
    states[2] = m.prelayer2(x2)
    for step = 1:m.steps
        state = states[m.inputindices[step][1]]
        in1 = m.ops[step][1](states[m.inputindices[step][1]])
        if m.names[step][1] != "skip_connect"
            in1 = droppath(in1, drop_prob)
        end
        in2 = m.ops[step][2](states[m.inputindices[step][2]])
        if m.names[step][2] != "skip_connect"
            in2 = droppath(in2, drop_prob)
        end
        states[step+2] = in1 + in2
    end
    reduce(ccat, copy(states[m.steps+2-m.multiplier+1:m.steps+2]))
end

Flux.@functor EvalCell

struct DARTSEvalModel
    normal_αs::Vector{Vector{Float32}}
    reduce_αs::Vector{Vector{Float32}}
    stem::Chain
    cells::Vector{EvalCell}
    global_pooling::AdaptiveMeanPool
    classifier::Dense
    operations::Vector{String}
end

function DARTSEvalModel(
    α_normal::Vector{Vector{Float32}},
    α_reduce::Vector{Vector{Float32}};
    num_cells::Int64 = 8,
    channels::Int64 = 16,
    num_classes::Int64 = 10,
    steps::Int64 = 4,
    mult::Int64 = 4,
    stem_mult::Int64 = 3,
    operations::Vector{String} = PRIMITIVES,
)
    channels_current = channels * stem_mult
    stem = Chain(
        Conv((3, 3), 3 => channels_current, pad = (1, 1), bias = false),
        BatchNorm(channels_current),
    ) |> gpu
    channels_before_last = channels_current
    channels_last = channels_current
    channels_current = channels
    reduce_previous = false
    cells = Vector{EvalCell}(undef,num_cells)
    for i = 1:num_cells
        if i == num_cells ÷ 3 + 1 || i == 2 * num_cells ÷ 3 + 1
            channels_current = channels_current * 2
            reduce = true
            αs = α_reduce
        else
            reduce = false
            αs = α_normal
        end
        cells[i] = EvalCell(
            channels_before_last,
            channels_last,
            channels_current,
            reduce,
            reduce_previous,
            steps,
            mult,
            αs,
            operations,
        )

        reduce_previous = reduce
        channels_before_last = channels_last
        channels_last = mult * channels_current
    end

    global_pooling = AdaptiveMeanPool((1, 1)) |> gpu
    classifier = Dense(channels_last, num_classes) |> gpu
    DARTSEvalModel(α_normal, α_reduce, stem, cells, global_pooling, classifier, operations)
end

function (m::DARTSEvalModel)(
    x::CuArray{Float32, 4},
    drop_prob::Float32 = 0f0,
)
    s1 = m.stem(x)
    s2 = m.stem(x)
    for (i, cell) in enumerate(m.cells)
        new_state = cell(s1, s2, drop_prob)
        s1 = s2
        s2 = new_state
    end
    out = m.global_pooling(s2)
    m.classifier(squeeze(out))
end

Flux.@functor DARTSEvalModel

struct DARTSEvalAuxModel
    normal_αs::Vector{Vector{Float32}}
    reduce_αs::Vector{Vector{Float32}}
    stem::Chain
    cells::Vector{EvalCell}
    auxiliary::Chain
    global_pooling::AdaptiveMeanPool
    classifier::Dense
    operations::Vector{String}
end

function DARTSEvalAuxModel(
    α_normal::Vector{Vector{Float32}},
    α_reduce::Vector{Vector{Float32}};
    num_cells::Int64 = 8,
    channels::Int64 = 16,
    num_classes::Int64 = 10,
    steps::Int64 = 4,
    mult::Int64 = 4,
    stem_mult::Int64 = 3,
    operations::Vector{String} = PRIMITIVES,
)
    channels_current = channels * stem_mult
    stem = Chain(
        Conv((3, 3), 3 => channels_current, pad = (1, 1), bias = false),
        BatchNorm(channels_current),
    ) |> gpu
    channels_before_last = channels_current
    channels_last = channels_current
    channels_current = channels
    channels_aux = channels
    reduce_previous = false
    cells = Vector{EvalCell}(undef,num_cells)
    for i = 1:num_cells
        if i == num_cells ÷ 3 + 1 || i == 2 * num_cells ÷ 3 + 1
            channels_current = channels_current * 2
            reduce = true
            αs = α_reduce
        else
            reduce = false
            αs = α_normal
        end
        cells[i] = EvalCell(
            channels_before_last,
            channels_last,
            channels_current,
            reduce,
            reduce_previous,
            steps,
            mult,
            αs,
            operations,
        )

        reduce_previous = reduce
        channels_before_last = channels_last
        channels_last = mult * channels_current
        if i == 2 * num_cells ÷ 3 + 1
            channels_aux = channels_last
        end
    end
    auxiliary = Chain(
        x -> relu.(x),
        MeanPool((5, 5), pad = 0, stride = 3),
        Conv((1, 1), channels_aux => 128, bias = false),
        BatchNorm(128),
        x -> relu.(x),
        Conv((2, 2), 128 => 768, bias = false),
        BatchNorm(768),
        x -> relu.(x),
        x -> dropdims(x, dims = (1,2)),
        Dense(768, num_classes),
    ) |> gpu
    global_pooling = AdaptiveMeanPool((1, 1)) |> gpu
    classifier = Dense(channels_last, num_classes) |> gpu
    DARTSEvalAuxModel(
        α_normal,
        α_reduce,
        stem,
        cells,
        auxiliary,
        global_pooling,
        classifier,
        operations,
    )
end

function DARTSEvalAuxModel(
    ;
    num_cells::Int64 = 8,
    channels::Int64 = 16,
    num_classes::Int64 = 10,
    steps::Int64 = 4,
    mult::Int64 = 4,
    stem_mult::Int64 = 3,
)
    channels_current = channels * stem_mult
    stem = Chain(
        Conv((3, 3), 3 => channels_current, pad = (1, 1), bias = false),
        BatchNorm(channels_current),
    ) |> gpu
    channels_before_last = channels_current
    channels_last = channels_current
    channels_current = channels
    channels_aux = channels
    reduce_previous = false
    cells = []
    for i = 1:num_cells
        if i == num_cells ÷ 3 + 1 || i == 2 * num_cells ÷ 3 + 1
            channels_current = channels_current * 2
            reduce = true
        else
            reduce = false
        end
        cell = EvalCell(
            channels_before_last,
            channels_last,
            channels_current,
            reduce,
            reduce_previous,
            steps,
            mult,
        )
        push!(cells, cell)

        reduce_previous = reduce
        channels_before_last = channels_last
        channels_last = mult * channels_current
        if i == 2 * num_cells ÷ 3 + 1
            channels_aux = channels_last
        end
    end
    auxiliary = Chain(
        x -> relu.(x), #inplace?
        MeanPool((5, 5), pad = 0, stride = 3),
        Conv((1, 1), channels_aux => 128, bias = false),
        BatchNorm(128),
        x -> relu.(x),
        Conv((2, 2), 128 => 768, bias = false),
        BatchNorm(768),
        x -> relu.(x),
        x -> dropdims(x, dims = (1,2)),
        Dense(768, num_classes),
    ) |> gpu
    global_pooling = AdaptiveMeanPool((1, 1)) |> gpu
    classifier = Dense(channels_last, num_classes) |> gpu
    DARTSEvalAuxModel(
        [],
        [],
        stem,
        cells,
        auxiliary,
        global_pooling,
        classifier,
        PRIMITIVES
    )
end

function (m::DARTSEvalAuxModel)(
    x::CuArray{Float32, 4},
    is_training::Bool = false,
    drop_prob::Float32 = 0f0,
)
    s1 = m.stem(x)
    s2 = m.stem(x)
    for (i, cell) in enumerate(m.cells)
        new_state = cell(s1, s2, Float32(drop_prob))
        s1 = s2
        s2 = new_state
        if i == 2 * length(m.cells) ÷ 3 + 1 && is_training
            out_aux = m.auxiliary(s2)
        end
    end
    out = m.global_pooling(s2)
    out = m.classifier(squeeze(out))
    if !is_training
        out_aux = out
    end
    out, out_aux
end

Flux.@functor DARTSEvalAuxModel

max_stars_repo_name: SuReLI/DifferentiableNAS

max_issues_repo_name: SuReLI/DifferentiableNAS

max_forks_repo_name: SuReLI/DifferentiableNAS


text:
# arbres binaires de recherche (immuables)

struct BST
    left
    value
    right
end

function add(x, t)
    if t == nothing
	BST(nothing, x, nothing)
    elseif x < t.value
	BST(add(x, t.left), t.value, t.right)
    elseif x > t.value
	BST(t.left, t.value, add(x, t.right))
    else
        t
    end
end

function mem(x, a)
    if x == a.value
        return true
    end
    if x < a.value && a.left != nothing
        return mem(x, a.left)
    end
    if a.right != nothing
        return mem(x, a.right)
    end
    false
end

function print_tree(a)
    if a == nothing
        return
    end
    print("(")
    if a.left != nothing
        print_tree(a.left)
    end
    print(a.value)
    if a.right != nothing
        print_tree(a.right)
    end
    print(")")
end

function main()
    tree = nothing
    for i = 1 : 9
        x = (55 * i) % 34
        tree = add(x, tree)
        print_tree(tree)
        print("\n")
    end
    if mem(8, tree) && !mem(0, tree) && mem(32, tree) && !mem(22, tree)
        println("ok")
    end
    tree = add(42, tree)
    tree = add(-1, tree)
    print_tree(tree)
    print("\n")
end

main()

ext: jl

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Geotherm.Geometry: Point2D, Rectangle

@testset "Test whether a point is in the rectangle" begin
    @test in(Point2D(6, 4), Rectangle(1, 10, 5, 20)) == false
    @test in(Point2D(6, 15), Rectangle(1, 10, 5, 20)) == true
    @test in(Point2D(5, 5), Rectangle(1, 10, 5, 20)) == true
    @test in(Point2D(1, 15), Rectangle(1, 10, 5, 20)) == true
    @test in(Point2D(1, 5), Rectangle(1, 10, 5, 20)) == true
end


text:
# Polynomial covariance function

"""
    Poly <: Kernel

Polynomial kernel (covariance)
```math
k(x,x') = σ²(xᵀx' + c)ᵈ
```
with signal standard deviation ``σ``, additive constant ``c``, and degree ``d``.
"""
mutable struct Poly <: Kernel
    "Constant"
    c::Float64
    "Signal variance"
    σ2::Float64
    "Degree of polynomial"
    deg::Int
    "Priors for kernel parameters"
    priors::Array

    """
        Poly(lc::Float64, lσ::Float64, deg::Int)

    Create `Poly` with signal standard deviation `exp(lσ)`, additive constant `exp(lc)`,
    and degree `deg`.
    """
    Poly(lc::Float64, lσ::Float64, deg::Int) = new(exp(lc), exp(2 * lσ), deg, [])
end

function KernelData(k::Poly, X::MatF64)
    XtX=X'*X
    LinearAlgebra.copytri!(XtX, 'U')
    LinIsoData(XtX)
end
kernel_data_key(k::Poly, X::MatF64) = "LinIsoData"

_cov(poly::Poly, xTy) = poly.σ2*(poly.c.+xTy).^poly.deg
function Statistics.cov(poly::Poly, x::VecF64, y::VecF64)
    K = _cov(poly, dot(x,y))
end
function cov!(cK::MatF64, poly::Poly, X::MatF64, data::LinIsoData)
    cK .= _cov(poly, data.XtX)
end
@inline @inbounds function cov_ij(poly::Poly, X::MatF64, data::LinIsoData, i::Int, j::Int, dim::Int)
    return _cov(poly, data.XtX[i, j])
end
Statistics.cov(poly::Poly, X::MatF64, data::LinIsoData) = _cov(poly, data.XtX)

get_params(poly::Poly) = Float64[log(poly.c), log(poly.σ2) / 2]
get_param_names(poly::Poly) = [:lc, :lσ]
num_params(poly::Poly) = 2

function set_params!(poly::Poly, hyp::VecF64)
    length(hyp) == 2 || throw(ArgumentError("Polynomial function has two parameters"))
    poly.c = exp(hyp[1])
    poly.σ2 = exp(2 * hyp[2])
end

@inline dk_dlc(poly::Poly, xTy::Float64) = poly.c*poly.deg*poly.σ2*(poly.c+xTy).^(poly.deg-1)
@inline dk_dlσ(poly::Poly, xTy::Float64) = 2 * _cov(poly,xTy)
@inline function dKij_dθp(poly::Poly, X::MatF64, i::Int, j::Int, p::Int, dim::Int)
    if p==1
        return dk_dlc(poly, dotij(X,i,j,dim))
    else
        return dk_dlσ(poly, dotij(X,i,j,dim))
    end
end
@inline function dKij_dθp(poly::Poly, X::MatF64, data::LinIsoData, i::Int, j::Int, p::Int, dim::Int)
    if p==1
        return dk_dlc(poly, data.XtX[i,j])
    else
        return dk_dlσ(poly, data.XtX[i,j])
    end
end


text:
using Base.Test
using PGM

struct NameV
    id::String
    v::Int
end

@testset "Lists" begin

@testset "GroupedList" begin

    gl0 = grlist_from_counts(String[], [0])
    @test eltype(gl0) == String
    @test length(gl0) == 0
    @test ngroups(gl0) == 1
    @test collect(gl0) == String[]
    @test group_length(gl0, 1) == 0
    @test collect(group(gl0, 1)) == String[]

    gl1 = grlist_from_counts(["a", "b", "c", "d", "e", "f"], [2, 1, 0, 3])
    @test eltype(gl1) == String
    @test length(gl1) == 6
    @test ngroups(gl1) == 4
    @test collect(gl1) == ["a", "b", "c", "d", "e", "f"]
    @test [group_length(gl1, k) for k = 1:4] == [2, 1, 0, 3]
    @test collect(group(gl1, 1)) == ["a", "b"]
    @test collect(group(gl1, 2)) == ["c"]
    @test collect(group(gl1, 3)) == String[]
    @test collect(group(gl1, 4)) == ["d", "e", "f"]

    gl2 = grlist_by([(10, 1), (20, 3), (30, 2), (40, 1), (50, 3), (60, 3)],
                    3, v->v[2])
    @test eltype(gl2) == Tuple{Int, Int}
    @test length(gl2) == 6
    @test ngroups(gl2) == 3
    @test collect(gl2) == [(10, 1), (40, 1), (30, 2), (20, 3), (50, 3), (60, 3)]
    @test [group_length(gl2, k) for k = 1:3] == [2, 1, 3]
    @test collect(group(gl2, 1)) == [(10, 1), (40, 1)]
    @test collect(group(gl2, 2)) == [(30, 2)]
    @test collect(group(gl2, 3)) == [(20, 3), (50, 3), (60, 3)]
end

@testset "NamedList" begin

    e1 = NameV("a", 12)
    e2 = NameV("x", 23)
    e3 = NameV("c", 36)

    nl1 = namedlist([e1, e2, e3], e -> e.id)
    @test eltype(nl1) == NameV
    @test ndims(nl1) == 1
    @test length(nl1) == 3
    @test size(nl1) == (3,)
    @test collect(nl1) == [e1, e2, e3]
    @test collect(eachindex(nl1)) == [1, 2, 3]

    @test nl1[1] == e1
    @test nl1[2] == e2
    @test nl1[3] == e3
    @test nl1["a"] == e1
    @test nl1["x"] == e2
    @test nl1["c"] == e3

    @test indexof("a", nl1) == 1
    @test indexof("x", nl1) == 2
    @test indexof("c", nl1) == 3
    @test indexof(("c", "a"), nl1) == [3, 1]
    @test_throws KeyError indexof("b", nl1)
    @test_throws KeyError namedlist([e1, e1], e -> e.id)

    nl2 = namedlist_from_pairs([("a", 12), ("x", 23), ("c", 36)])
    @test eltype(nl2) == Int
    @test ndims(nl2) == 1
    @test length(nl2) == 3
    @test size(nl2) == (3,)
    @test collect(nl2) == [12, 23, 36]
    @test collect(eachindex(nl2)) == [1, 2, 3]

    @test nl2[1] == 12
    @test nl2[2] == 23
    @test nl2[3] == 36
    @test nl2["a"] == 12
    @test nl2["x"] == 23
    @test nl2["c"] == 36
end

end  # Lists


text:
function model(du,u,p,t)
    ts = [50.0, 100.0, 200.0, 300.0]
    freq = 0.2
    amp = p["amp"]
    if t < ts[1]
        p["maxA"] = 0
    elseif t < ts[3]
        p["maxA"] = 2.3
    else
        p["maxA"] = 2.3*0.65
        p["gA"] = 0.005
    end
    if t < ts[2]
        p["maxF"] = 0
    elseif t < ts[3]
        p["maxF"] = 2.1
    else
        p["maxF"] = 2.1*0.65
        p["gF"] = 0.005
    end
    if t < ts[4]
        p["muA"] = 0
        p["muF"] = 0
    else
        Random.seed!(floor(Int, (t-ts[4])*freq)*p["y"]*floor(Int, p["x"]))
        p["muA"] =   sign(sin((t-ts[4])*pi*freq)) * amp * randn()*u[5]
        p["muF"] =  sign(sin((t-ts[4])*pi*freq)) * amp * randn()*u[6]
    end
    P5, P8, R5, R8, A, F, r2, c15, c25, c35, r9, c18, c28, c38 = u
    x, gA, gF, maxA, maxF, s05, s08, K5A, K8F, n5A, n8F, K58, n58, K85, n85  = p["x"], p["gA"], p["gF"], p["maxA"], p["maxF"], p["s05"], p["s08"], p["K5A"], p["K8F"], p["n5A"], p["n8F"], p["K58"], p["n58"], p["K85"], p["n85"]
    kR5, kR8, l05, l08, kP5, kP8 = p["kR5"], p["kR8"], p["l05"], p["l08"], p["kP5"], p["kP8"]
    s5, s8  = p["s5"], p["s8"]
    muA, muF = p["muA"], p["muF"]
    du[1] = dP5 = l05 * R5 - kP5 * P5
    du[2] = dP8 = l08 * R8 - kP8 * P8
    du[3] = dR5 = s5*((A/K5A)^n5A)/(1+(P8/K58)^n58+(A/K5A)^n5A) + s05 - kR5 * R5
    du[4] = dR8 = s8*((F/K8F)^n8F)/(1+(F/K8F)^n8F) + s08 - kR8 * R8
    du[5] = dA = gA*(maxA * exp(-0.01 * x) - A) + muA
    du[6] = dF = gF*(maxF * exp(-0.01 * (99 - x)) - F) + muF
    du[7] = dr2 = -r2
    du[8] = dc15 = -c15
    du[9] = dc25 = -c25
    du[10] = dc35 = -c35
    du[11] = dr9 = -r9
    du[12] = dc18 = -c18
    du[13] = dc28 = -c28
    du[14] = dc38 = -c38
end

d_ics1 = Dict("R8"=> 0.00900023576123889, "c15"=> 1.761337400332505e-06, "P8"=> 2.9560034347325604e-06, "P5"=> 3.4020485285993286e-07, "c25"=> 9.118941517979407e-05, "c35"=> 0.004721137474792935, "A"=> 1.5e-323, "c18"=> 8.844673354103184e-06, "R5"=> 0.000500023861019239, "c38"=> 0.001867111128611335, "F"=> 1.5e-323, "r2"=> 0.3524443000060298, "r9"=> 0.33001432105772416, "c28"=> 0.0002687664938900425)
svlist = ["P5", "P8", "R5", "R8", "A", "F", "r2", "c15", "c25", "c35", "r9", "c18","c28", "c38"]
u0 = [d_ics1[x] for x in svlist]
tspan = (0.0,1000.0)
p=Dict("x"=>50
    ,"amp"=>0.001
    ,"kR5"=>1.0
    ,"l05"=>1.0
    ,"kP5"=>1.0
    ,"maxA"=>2.3
    ,"maxF"=>2.1
    ,"kR8"=>1.0
    ,"l08"=>1.0
    ,"kP8"=>1.0
    ,"gA"=>1.0
    ,"gF"=>1.0
    ,"gc"=>1.0
    ,"s05"=>0.0002
    ,"s5"=>1.0
    ,"s08"=>0.0002
    ,"s8"=>1.0
    ,"K5A"=>0.7
    ,"n5A"=>60
    ,"K58"=>0.05
    ,"n58"=>20
    ,"K85"=>1000
    ,"n85"=>6.0
    ,"K8F"=>0.75
    ,"n8F"=>60
    ,"y"=>1)


text:
using NPZ
using Convex
using ECOS
using LinearAlgebra
using LineSearches
using ProgressBars
using WoodburyMatrices
using Printf: @sprintf
using .Threads

G1 = npzread("Kernels/G1.npy")
G2 = npzread("Kernels/G2.npy")
G3 = npzread("Kernels/G3.npy")
nfactor = npzread("Testing/nfactor.npy")

function d_oed(X1::Matrix{T}, X2::Matrix{T}, X3::Matrix{T}, nfactor, B, ϵ, λ; maxiter=2500, ls=MoreThuente()) where T <: Number
    #initial vector
    X = [X3 zeros(T, size(X3));
         zeros(T, size(X3)) X3;
         nfactor*X1 nfactor*X2]
    Xn = maximum(sqrt.(sum(X.^2, dims=2)))
    X ./= Xn

    n = size(X3, 1)
    na = size(X, 1)
    μr = ones(T, n)./n
    μr_tmp = zeros(T, n)
    μp = Variable(n)
    μ = ones(T, na)
    μ[1:n] .= μr
    μ[(n+1):(2*n)] .= μr
    gvec = zeros(T, size(μ))
    s = zeros(T, size(μ))

    function project!(μ, μr, μp, B)
            μr .= μ[1:n]
            projection = minimize(sumsquares(μp-μr), [0<=μp, μp<=1, sum(μp)<=B])
            solve!(projection, () -> ECOS.Optimizer(verbose=false))
            μr .= @views μp.value[:,1]
            μ[1:n] .= μr
            μ[(n+1):(2*n)] .= μr
            return μ
    end

    project!(μ) = project!(μ, μr_tmp, μp, B)

    #setup linesearch functions
    #linesearch - this needs to make sure that any interior functions actually have function barriers around global variables
    function Ml(μ, X, ϵ)
        M = ϵ*I + X'*Diagonal(μ)*X
        Mc = cholesky(Hermitian(M))
        l = μ.^2 .*diag(X*(Mc\X'))
        return (M, Mc, l)
    end

    Ml(μ) = Ml(μ, X, ϵ)

    function f(μ, λ)
        M, Mc, l = Ml(μ)
        return -logdet(M) + λ*maximum(l)
    end

    function ff(μ, X, ϵ, λ)
        M = ϵ*I + X'*Diagonal(μ)*X
        l = (μ.^2).*diag(X*(M\X'))
        return -logdet(M) + λ*maximum(l)
    end

    f(μ) = f(μ, λ)

    ff(μ) = ff(μ, X, ϵ, λ)

    function gi(i, Mc, μ, u, lmind, X, λ, T)
        ld_term = @views -X[i,:]'*(Mc\X[i,:])
        tmp_A = SymWoodbury(Mc, X[i,:], -μ[i])
        tmp_b = @views tmp_A \ X[i,:]
        tmp_1 = (u'*tmp_b)^2
        tmp_d = @views X[i,:]'*tmp_b
        tmp_den = (1+μ[i]*tmp_d)^2
        lev_term_1 =  -μ[lmind]*((1+μ[i]*tmp_d)*tmp_1 - μ[i]*tmp_1*tmp_d) / tmp_den
        lev_term_2 = (i==lmind) ? u'*(Mc\u) : zero(T)
        return ld_term + λ*(lev_term_1 + lev_term_2)
    end

    gi(i, Mc, μ, u, lmind) = gi(i, Mc, μ, u, lmind, X, λ, T)

    function g!(gvec, μ, X, n)
        M, Mc, l = Ml(μ)
        lm, lmind = findmax(l)
        u = @views X[lmind, :]
        @threads for i = 1:n
            g1 = gi(i, Mc, μ, u, lmind) # in respect to x component
            g2 = gi(i+n, Mc, μ, u, lmind) # in respect to y component
            gvec[i] = g1+g2
            gvec[i+n] = g1+g2
        end
    end

    g!(gvec, μ) = g!(gvec, μ, X, n)

    function fg!(gvec, μ, X, n)
        M, Mc, l = Ml(μ)
        lm, lmind = findmax(l)
        u = @views X[lmind, :]
        @threads for i = 1:n
            g1 = gi(i, Mc, μ, u, lmind) # in respect to x component
            g2 = gi(i+n, Mc, μ, u, lmind) # in respect to y component
            gvec[i] = g1+g2
            gvec[i+n] = g1+g2
        end
        return -logdet(M) + λ*maximum(l)
    end

    fg!(gvec, μ) = fg!(gvec, μ, X, n)

    ϕ(α) = f(project!(μ.+α.*s))

    function dϕ(α)
        g!(gvec, project!(μ.+α.*s))
        return dot(gvec, s)
    end

    function ϕdϕ(α)
        phi = fg!(gvec, project!(μ.+α.*s))
        dphi = dot(gvec, s)
        return (phi, dphi)
    end

    #get initial gradient and iterate
    g!(gvec, μ)

    fμ = f(μ)
    fμt = []
    μrarr = []
    push!(fμt, fμ)
    iter = ProgressBar(1:maxiter)
    set_description(iter, string(@sprintf("Starting Up")))

    for i in iter
        #set step direction
        s .= -gvec
        dϕ_0 = dot(s, gvec)
        #linesearch gradient
        α, fμ = ls(ϕ, dϕ, ϕdϕ, 1.0, fμ, dϕ_0)
        @. μ = μ + α*s
        project!(μ, μr, μp, B)
        fμ = f(μ)
        #update gradient
        g!(gvec, μ)
        set_description(iter, string(@sprintf("Objective: %.5f Update: %.2e", fμ, fμ-fμt[i])))
        push!(fμt, fμ)
        push!(μrarr, copy(μr))
    end

    function pround(μ, μr, n, B)
        idx = 1:length(μr)
        rμr = copy(μr)
        v1 = copy(μr)
        v2 = copy(μr)
        μ1 = copy(μ)
        μ2 = copy(μ)
        not_zero_or_one(x) = !(isapprox(x,0) || isapprox(x,1))
        feasible(v) = (all(v.>=0) && all(v.<=1) && (sum(v)<=B))
        while length(idx[not_zero_or_one.(rμr)]) > 1
            # j, k = rand(idx[not_zero_or_one.(rμr)], 2)
            idr = idx[not_zero_or_one.(rμr)]
            j = idr[argmax(rμr[idr])]
            k = idr[argmin(rμr[idr])]
            println(j, " ", k)
            δ = min(1-rμr[j], rμr[k])
            τ = min(1-rμr[k], rμr[j])
            v1[j] += δ
            v1[k] -= δ
            v2[j] -= τ
            v2[k] += τ
            sv1 = sum(v1)
            sv2 = sum(v2)
            if !feasible(v1) && !feasible(v2)
                error("No Feasible solution")
            elseif !feasible(v1)
                rμr .= v2
                v1 .= v2
            elseif !feasible(v2)
                rμr .= v1
                v2 .= v1
            else
                μ1[1:n] .= v1
                μ1[(n+1):(2*n)] .= v1
                f1 = ff(μ1)
                μ2[1:n] .= v2
                μ2[(n+1):(2*n)] .= v2
                f2 = ff(μ2)
                if f1 < f2
                    rμr .= v1
                else
                    rμr .= v2
                end
            end
        end
        μ1[1:n] .= rμr
        μ1[(n+1):(2*n)] .= rμr
        fr = ff(μ1)
        return (rμr, fr)
    end

    rμr, fr = pround(μ, μr, n, B)

    return μr, fμt, μrarr, Xn, rμr, fr
end

μr, fμt, μrarr, Xn, rμr,  fr = d_oed(G1, G2, G3, nfactor, 32, 1e-1, 10.0, maxiter=250)

npzwrite("resmin.npy", μr)
npzwrite("resround.npy", rμr)
;

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
using CubicSplineNR
using Base.Test

x = collect(linspace(0, 2pi, 30))
fx = sin.(x)
fx2 = cos.(x)

sin_itp = CubicSpline(x, fx)
sin_mut = CubicSpline(x, fx, mutable=true)

sin_itp(3.14)
sin_mut(3.14)

@time update!(sin_mut, fx2)

sin_mut(3.14)


text:
using fci
using JLD2

config = [1,2,3]
sign = 1

count = 2 #from bubble sort
arr = [2,3,4] #from bubble sort
positions = [1,4]

@testset "excit config" begin
    spot = first(findall(x->x==positions[1], config))
    config[spot] = positions[2]
    config_org = deepcopy(config)
    count_test, arr_test = fci.bubble_sort(config)
    @test count_test == count
    @test arr_test == arr
    if iseven(count_test)
        sign_test = 1
    else
        sign_test = -1
    end
    @test sign == sign_test
end

max_stars_repo_name: nbraunsc/fci

max_issues_repo_name: nbraunsc/fci

max_forks_repo_name: nbraunsc/fci


text:
export rational_reconstruction, farey_lift, div, valence, leading_coefficient,
       trailing_coefficient, constant_coefficient, factor_mod_pk,
       factor_mod_pk_init, hensel_lift, resultant_sircana, rres, rresx,
       coefficients

function PolynomialRing(R::Ring; cached::Bool = false)
  return PolynomialRing(R, "_x", cached = cached)
end

function PolynomialRing(R::FlintRationalField, a::Symbol; cached::Bool = true)
  Qx = FmpqPolyRing(R, a, cached)
  return Qx, gen(Qx)
end

function PolynomialRing(R::FlintIntegerRing, a::Symbol; cached::Bool = true)
  Zx = FmpzPolyRing(a, cached)
  return Zx, gen(Zx)
end

function FlintFiniteField(p::Integer)
  @assert isprime(p)
  return ResidueRing(FlintZZ, p, cached=false)
end

function FlintFiniteField(p::fmpz)
  @assert isprime(p)
  return ResidueRing(FlintZZ, p, cached=false)
end

function fmpz(a::Generic.Res{Nemo.fmpz})
  return a.data
end

function fmpz(a::Nemo.nmod)
  return fmpz(a.data)
end

function lift(::FlintIntegerRing, a::Generic.Res{Nemo.fmpz})
  return a.data
end

function (::FlintIntegerRing)(a::Generic.Res{Nemo.fmpz})
  return a.data
end

function lift(::FlintIntegerRing, a::Nemo.nmod)
  return fmpz(a.data)
end

function (::FlintIntegerRing)(a::Nemo.nmod)
  return fmpz(a.data)
end

function div(f::PolyElem, g::PolyElem)
  q,r = divrem(f,g)
  return q
end

function ismonic(a::PolyElem)
  return isone(lead(a))
end

doc"""
***
  valence(f::PolyElem) -> RingElem

>  The last non-zero coefficient of f
"""
function valence(f::PolyElem)
  for i=0:degree(f)
    c = coeff(f, i)
    if !iszero(c)
      return c
    end
  end
  return c
end

doc"""
***
  leading_coefficient(f::PolyElem) -> RingElem

>  The last leading coefficient of f.
"""
leading_coefficient(f::PolyElem) = lead(f)

doc"""
***
  trailing_coefficient(f::PolyElem) -> RingElem
  constant_coefficient(f::PolyElem) -> RingElem

>  The constant coefficient of f.
"""
function trailing_coefficient(f::PolyElem)
  if iszero(f)
    return base_ring(f)(0)
  end
  return coeff(f, 0)
end

doc"""
    induce_rational_reconstruction(a::fmpz_poly, M::fmpz) -> fmpq_poly
> Apply {{{rational_reconstruction}}} to each coefficient of $a$, resulting
> in either a fail (return (false, s.th.)) or (true, g) for some rational
> polynomial $g$ s.th. $g \equiv a \bmod M$.
"""
function induce_rational_reconstruction(a::fmpz_poly, M::fmpz)
  b = PolynomialRing(FlintQQ, parent(a).S)[1]()
  for i=0:degree(a)
    fl, x,y = rational_reconstruction(coeff(a, i), M)
    if fl
      setcoeff!(b, i, x//y)
    else
      return false, b
    end
  end
  return true, b
end

constant_coefficient = trailing_coefficient

function resultant(f::fmpz_poly, g::fmpz_poly, d::fmpz, nb::Int)
  z = fmpz()
  ccall((:fmpz_poly_resultant_modular_div, :libflint), Void,
     (Ptr{fmpz}, Ptr{fmpz_poly}, Ptr{fmpz_poly}, Ptr{fmpz}, Int),
     &z, &f, &g, &d, nb)
  return z
end

##############################################################
#
# Hensel
#
##############################################################

mutable struct fmpz_poly_raw  ## fmpz_poly without parent like in c
  coeffs::Ptr{Void}
  alloc::Int
  length::Int

  function fmpz_poly_raw()
    error("should not get called")
    z = new()
    ccall((:fmpz_poly_init, :libflint), Void, (Ptr{fmpz_poly},), &z)
    finalizer(z, _fmpz_poly_raw_clear_fn)
    return z
  end
end

function _fmpz_poly_raw_clear_fn(a::fmpz_poly)
  ccall((:fmpz_poly_clear, :libflint), Void, (Ptr{fmpz_poly},), &a)
end

mutable struct fmpz_poly_factor
  c::Int   # really an fmpz  - but there is no fmpz_raw to be flint compatible
  poly::Ptr{fmpz_poly_raw}
  exp::Ptr{Int}
  _num::Int
  _alloc::Int

  function fmpz_poly_factor()
    z = new()
    ccall((:fmpz_poly_factor_init, :libflint), Void,
            (Ptr{fmpz_poly_factor}, ), &z)
    finalizer(z, _fmpz_poly_factor_clear_fn)
    return z
  end
end

function _fmpz_poly_factor_clear_fn(a::fmpz_poly_factor)
  ccall((:fmpz_poly_factor_clear, :libflint), Void,
          (Ptr{fmpz_poly_factor}, ), &a)
end

function factor_to_dict(a::fmpz_poly_factor)
  res = Dict{fmpz_poly,Int}()
  Zx,x = PolynomialRing(FlintZZ, "x")
  for i in 1:a._num
    f = Zx()
    ccall((:fmpz_poly_set, :libflint), Void, (Ptr{fmpz_poly}, Ptr{fmpz_poly_raw}), &f, a.poly+(i-1)*sizeof(fmpz_poly_raw))
    res[f] = unsafe_load(a.exp, i)
  end
  return res
end

function show(io::IO, a::fmpz_poly_factor)
  ccall((:fmpz_poly_factor_print, :libflint), Void, (Ptr{fmpz_poly_factor}, ), &a)
end

mutable struct HenselCtx
  f::fmpz_poly
  p::UInt

  LF :: fmpz_poly_factor
  link::Ptr{Int}
  v::Ptr{fmpz_poly_raw}
  w::Ptr{fmpz_poly_raw}
  N::UInt
  prev::UInt
  r::Int  #for the cleanup
  lf:: Nemo.nmod_poly_factor

  function HenselCtx(f::fmpz_poly, p::fmpz)
    a = new()
    a.f = f
    a.p = UInt(p)
    Zx,x = PolynomialRing(FlintZZ, "x", cached=false)
    Rx,x = PolynomialRing(ResidueRing(FlintZZ, p, cached=false), "x", cached=false)
    a.lf = Nemo.nmod_poly_factor(UInt(p))
    ccall((:nmod_poly_factor, :libflint), UInt,
          (Ptr{Nemo.nmod_poly_factor}, Ptr{nmod_poly}), &(a.lf), &Rx(f))
    r = a.lf.num
    a.r = r
    a.LF = fmpz_poly_factor()
    @assert r > 1  #flint restriction
    a.v = ccall((:flint_malloc, :libflint), Ptr{fmpz_poly_raw}, (Int, ), (2*r-2)*sizeof(fmpz_poly_raw))
    a.w = ccall((:flint_malloc, :libflint), Ptr{fmpz_poly_raw}, (Int, ), (2*r-2)*sizeof(fmpz_poly_raw))
    for i=1:(2*r-2)
      ccall((:fmpz_poly_init, :libflint), Void, (Ptr{fmpz_poly_raw}, ), a.v+(i-1)*sizeof(fmpz_poly_raw))
      ccall((:fmpz_poly_init, :libflint), Void, (Ptr{fmpz_poly_raw}, ), a.w+(i-1)*sizeof(fmpz_poly_raw))
    end
    a.link = ccall((:flint_calloc, :libflint), Ptr{Int}, (Int, Int), 2*r-2, sizeof(Int))
    a.N = 0
    a.prev = 0
    finalizer(a, HenselCtx_free)
    return a
  end

  function free_fmpz_poly_array(p::Ptr{fmpz_poly_raw}, r::Int)
    for i=1:(2*r-2)
      ccall((:fmpz_poly_clear, :libflint), Void, (Ptr{fmpz_poly_raw}, ), p+(i-1)*sizeof(fmpz_poly_raw))
    end
    ccall((:flint_free, :libflint), Void, (Ptr{fmpz_poly_raw}, ), p)
  end
  function free_int_array(a::Ptr{Int})
    ccall((:flint_free, :libflint), Void, (Ptr{Int}, ), a)
  end
  function HenselCtx_free(a::HenselCtx)
    free_fmpz_poly_array(a.v, a.r)
    free_fmpz_poly_array(a.w, a.r)
    free_int_array(a.link)
  end
end

function show(io::IO, a::HenselCtx)
  println("factorisation of $(a.f) modulo $(a.p)^$(a.N)")
  if a.N > 0
    d = factor_to_dict(a.LF)
    println("currently: $d")
  end
end

function start_lift(a::HenselCtx, N::Int)
  a.prev = ccall((:_fmpz_poly_hensel_start_lift, :libflint), UInt,
       (Ptr{fmpz_poly_factor}, Ptr{Int}, Ptr{fmpz_poly_raw}, Ptr{fmpz_poly_raw}, Ptr{fmpz_poly}, Ptr{Nemo.nmod_poly_factor}, Int),
       &a.LF, a.link, a.v, a.w, &a.f, &a.lf, N)
  a.N = N
end

function continue_lift(a::HenselCtx, N::Int)
  a.prev = ccall((:_fmpz_poly_hensel_continue_lift, :libflint), Int,
       (Ptr{fmpz_poly_factor}, Ptr{Int}, Ptr{fmpz_poly_raw}, Ptr{fmpz_poly_raw}, Ptr{fmpz_poly}, UInt, UInt, Int, Ptr{fmpz}),
       &a.LF, a.link, a.v, a.w, &a.f, a.prev, a.N, N, &fmpz(a.p))
  a.N = N
end

doc"""
***
  factor_mod_pk(f::fmpz_poly, p::Int, k::Int) -> Dict{fmpz_poly, Int}

>  For f that is square-free modulo p, return the factorisation modulo p^k.
"""
function factor_mod_pk(f::fmpz_poly, p::Int, k::Int)
  H = HenselCtx(f, fmpz(p))
  start_lift(H, k)
  return factor_to_dict(H.LF)
end

doc"""
***
  factor_mod_pk_init(f::fmpz_poly, p::Int) -> HenselCtx

>  For f that is square-free modulo p, return a structure that allows to compute
>  the factorisaion modulo p^k for any k
"""
function factor_mod_pk_init(f::fmpz_poly, p::Int)
  H = HenselCtx(f, fmpz(p))
  return H
end

doc"""
***
  factor_mod_pk(H::HenselCtx, k::Int) -> RingElem

>  Using the result of factor_mod_pk_init, return a factorisation modulo p^k
"""
function factor_mod_pk(H::HenselCtx, k::Int)
  @assert k>= H.N
  if H.N == 0
    start_lift(H, k)
  else
    continue_lift(H, k)
  end
  return factor_to_dict(H.LF)
end

#I think, experimentally, that p = Q^i, p1 = Q^j and j<= i is the condition to make it tick.
function hensel_lift!(G::fmpz_poly, H::fmpz_poly, A::fmpz_poly, B::fmpz_poly, f::fmpz_poly, g::fmpz_poly, h::fmpz_poly, a::fmpz_poly, b::fmpz_poly, p::fmpz, p1::fmpz)
  ccall((:fmpz_poly_hensel_lift, :libflint), Void, (Ptr{fmpz_poly}, Ptr{fmpz_poly},  Ptr{fmpz_poly},  Ptr{fmpz_poly},  Ptr{fmpz_poly},  Ptr{fmpz_poly},  Ptr{fmpz_poly}, Ptr{fmpz_poly}, Ptr{fmpz_poly}, Ptr{fmpz}, Ptr{fmpz}), &G, &H, &A, &B, &f, &g, &h, &a, &b, &p, &p1)
end

doc"""
***
  hensel_lift(f::fmpz_poly, g::fmpz_poly, h::fmpz_poly, p::fmpz, k::Int) -> (fmpz_poly, fmpz_poly)

>  Given f = gh modulo p for g, h coprime modulo p, compute G, H s.th. f = GH mod p^k and
>  G = g mod p, H = h mod p.
"""
function hensel_lift(f::fmpz_poly, g::fmpz_poly, h::fmpz_poly, p::fmpz, k::Int)
  Rx, x = PolynomialRing(ResidueRing(FlintZZ, p, cached=false), cached=false)
  fl, a, b = gcdx(Rx(g), Rx(h))
  @assert isone(fl)
  @assert k>= 2
  ## if one of the cofactors is zero, this crashes.
  ## this can only happen if one of the factors is one. In this case, the other
  ## is essentially f and f would be a legal answer. Probably reduced mod p^k
  ## with all non-negative coefficients
  ## for now:
  @assert !iszero(a) && !iszero(b)
  a = lift(parent(g), a)
  b = lift(parent(g), b)
  G = parent(g)()
  H = parent(g)()
  A = parent(g)()
  B = parent(g)()
  g = deepcopy(g)
  h = deepcopy(h)

  # the idea is to have a good chain of approximations, ie.
  # to reach p^10, one should do p, p^2, p^3, p^5, p^10
  # rather than p, p^2, p^4, p^8, p^10
  # the chain has the same length, but smaller entries.
  l = [k]
  while k>1
    k = div(k+1, 2)
    push!(l, k)
  end
  ll = []
  for i=length(l)-1:-1:1
    push!(ll, l[i] - l[i+1])
  end
  P = p
  for i in ll
    p1 = p^i
    hensel_lift!(G, H, A, B, f, g, h, a, b, P, p1)
    G, g = g, G
    H, h = h, H
    A, a = a, A
    B, b = b, B
    P *= p1
  end
  return g, h
end

doc"""
***
  hensel_lift(f::fmpz_poly, g::fmpz_poly, p::fmpz, k::Int) -> (fmpz_poly, fmpz_poly)

>  Given f and g such that g is a divisor of f mod p and g and f/g are coprime, compute a hensel lift of g modulo p^k.
"""
function hensel_lift(f::fmpz_poly, g::fmpz_poly, p::fmpz, k::Int)
  Rx, x = PolynomialRing(ResidueRing(FlintZZ, p, cached=false), cached=false)
  h = lift(parent(f), div(Rx(f), Rx(g)))
  return hensel_lift(f, g, h, p, k)[1]
end

function fmpq_poly_to_nmod_poly_raw!(r::nmod_poly, a::fmpq_poly)
  ccall((:_fmpz_vec_get_nmod_poly, :libhecke), Void, (Ptr{nmod_poly}, Ptr{Int}, Int), &r, a.coeffs, a.length)
  p = r.mod_n
  den = ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ptr{Int}, UInt), &a.den, p)
  if den != UInt(1)
    den = ccall((:n_invmod, :libflint), UInt, (UInt, UInt), den, p)
    mul!(r, r, den)
  end
end

function fmpq_poly_to_fmpz_mod_poly_raw!(r::fmpz_mod_poly, a::fmpq_poly, t1::fmpz_poly = fmpz_poly(), t2::fmpz = fmpz())
  ccall((:fmpq_poly_get_numerator, :libflint), Void, (Ptr{fmpz_poly}, Ptr{fmpq_poly}), &t1, &a)
  ccall((:fmpz_mod_poly_set_fmpz_poly, :libflint), Void, (Ptr{fmpz_mod_poly}, Ptr{fmpz_poly}), &r, &t1)
  ccall((:fmpq_poly_get_denominator, :libflint), Void, (Ptr{fmpz}, Ptr{fmpq_poly}), &t2, &a)
  if !isone(t2)
    res = ccall((:fmpz_invmod, :libflint), Cint, (Ptr{fmpz}, Ptr{fmpz}, Ptr{fmpz}), &t2, &t2, &base_ring(r).modulus)
    @assert res != 0
    ccall((:fmpz_mod_poly_scalar_mul_fmpz, :libflint), Void, (Ptr{fmpz_mod_poly}, Ptr{fmpz_mod_poly}, Ptr{fmpz}), &r, &r, &t2)
  end
end

function fmpq_poly_to_nmod_poly(Rx::Nemo.NmodPolyRing, f::fmpq_poly)
  g = Rx()
  fmpq_poly_to_nmod_poly_raw!(g, f)
  return g
end

function fmpz_poly_to_nmod_poly_raw!(r::nmod_poly, a::fmpz_poly)
  ccall((:fmpz_poly_get_nmod_poly, :libflint), Void,
                  (Ptr{nmod_poly}, Ptr{fmpz_poly}), &r, &a)

end

function fmpz_poly_to_nmod_poly(Rx::Nemo.NmodPolyRing, f::fmpz_poly)
  g = Rx()
  fmpz_poly_to_nmod_poly_raw!(g, f)
  return g
end

function fmpq_poly_to_fmpz_mod_poly(Rx::Nemo.FmpzModPolyRing, f::fmpq_poly)
  g = Rx()
  fmpq_poly_to_fmpz_mod_poly_raw!(g, f)
  return g
end

function fmpz_poly_to_fmpz_mod_poly_raw!(r::fmpz_mod_poly, a::fmpz_poly)
  ccall((:fmpz_poly_get_fmpz_mod_poly, :libflint), Void,
                  (Ptr{fmpz_mod_poly}, Ptr{fmpz_poly}), &r, &a)

end

function fmpz_poly_to_fmpz_mod_poly(Rx::Nemo.FmpzModPolyRing, f::fmpz_poly)
  g = Rx()
  fmpz_poly_to_fmpz_mod_poly_raw!(g, f)
  return g
end

doc"""
    deflate(f::PolyElem, n::Int64) -> PolyElem
> Given a polynomial $f$ in $x^n$, write it as a polynomial in $x$, ie. divide
> all exponents by $n$.
"""
function deflate(x::PolyElem, n::Int64)
  y = parent(x)()
  for i=0:div(degree(x), n)
    setcoeff!(y, i, coeff(x, n*i))
  end
  return y
end

doc"""
    inflate(f::PolyElem, n::Int64) -> PolyElem
> Given a polynomial $f$ in $x$, return $f(x^n)$, ie. multiply
> all exponents by $n$.
"""
function inflate(x::PolyElem, n::Int64)
  y = parent(x)()
  for i=0:degree(x)
    setcoeff!(y, n*i, coeff(x, i))
  end
  return y
end

doc"""
    deflate(x::PolyElem) -> PolyElem
> Deflate the polynomial $f$ maximally, ie. find the largest $n$ s.th.
> $f$ can be deflated by $n$, ie. $f$ is actually a polynomial in $x^n$.
"""
function deflate(x::PolyElem)
  g = 0
  for i=0:degree(x)
    if coeff(x, i) != 0
      g = gcd(g, i)
      if g==1
        return x, 1
      end
    end
  end
  return deflate(x, g), g
end

function _hensel(f::PolyElem{T}, g::PolyElem{T}, h::PolyElem{T}, s::PolyElem{T}, t::PolyElem{T}) where T <: RingElem #from von zur Gathen: h needs to be monic
  @assert ismonic(h)
#  @assert isnilpotent(content(f-g*h))  #this guarantees useful lifting
#  @assert isnilpotent(content(g*s+h*t-1))
  e = f-g*h
  q, r = divrem(s*e, h)
#  @assert s*e == q*h+r
  g = g+t*e+q*g
  h = h+r
  @assert ismonic(h)
  b = s*g+t*h-1
  c, d = divrem(s*b, h)
#  @assert s*b == c*h+d
  s = s-d
  t = t-t*b-c*g

  return g, h, s, t
end

#factors f as unit * monic
#seems to not like 3*x^2+3*x+1 mod 108! mod 27 it is fine.
# requires some coefficient of f to be a unit
function fun_factor(f::PolyElem{<:RingElem})
  local g0
  local h0
  if isunit(lead(f))
    l= lead(f)
    return parent(f)(l), f*inv(l)
  end
  if isunit(f)
    return f, parent(f)(1)
  end
  t = gen(parent(f))
  g0 = parent(f)(0)
  for i=degree(f):-1:0
    if isunit(coeff(f, i))
      h0 = f - g0
      break
    else
      setcoeff!(g0, i, deepcopy(coeff(f, i)))
    end
  end

  #co-factors:
  g0 = parent(f)(0)
  setcoeff!(g0, degree(f) - degree(h0), lead(f))
  setcoeff!(g0, 0, lead(h0))
  s = parent(f)(inv(lead(h0)))
  h0 *= lead(s)
  t = parent(f)(0)

  g, h, s,t = _hensel(f, g0, h0, s, t)
  i = 1
  while g!= g0 || h != h0
#      @show content(g*s-t*h-1)
    i += 1
    g0 = g
    h0 = h
    g, h, s, t = _hensel(f, g0, h0, s, t)
    if i>20 #in general: loop forever... one could check that the
      # contents of f-gh and s*g+t*h - 1 is nilpotent.
      # this SHOULD ensure convergence
      @show f
      @show parent(f)
      @show g, h, s, t
      @show content(f-g*h)
      @show content(g*s-t*h-1)
      error("too long")
    end
  end

  return g0, h0
end

doc"""
    Base.rand(Rt::PolyRing{T}, n::Int) where T <: ResElem{S} where S <: Union{fmpz, Integer} -> PolElem{T}
> Find a random polynomial of degree(n)
"""
function Base.rand(Rt::PolyRing{T}, n::Int) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  f = Rt()
  R = base_ring(Rt)
  for i=0:n
    setcoeff!(f, i, rand(R))
  end
  return f
end

doc"""
    resultant_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer} -> T
> The resultant of $f$ anf $g$ using a quadratic-time algorithm.
"""
function resultant_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Rt = parent(f)
  R = base_ring(Rt)
  m = fmpz(modulus(R))
  e, p = ispower(m)
  easy = isprime(p)

  Zx = PolynomialRing(FlintZZ)[1]

  res = R(1)

  while true
    if degree(f) < 1 && degree(g) < 1
      if iszero(f) || iszero(g)
        res *= R(0)
      else
        res *= R(1)
      end
      return res
    end

    if degree(f) < 1
      res *= lead(f)^degree(g)
      return res
    end

    if degree(g) < 1
      res *= lead(g)^degree(f)
      return res
    end

    c, f = primsplit(f)
    if !isone(c)
      res *= c^degree(g)
    end

    c, g = primsplit(g)
    if !isone(c)
      res *= c^degree(f)
    end

    if degree(f) < degree(g)
      res *= (-1)^(degree(f)*degree(g))
      f, g = g, f
    end

    #want f % g which works iff lead(g) | lead(f)

    if isunit(lead(g)) #accelerate a bit...possibly.
      q, r = divrem(f, g)
      res *= lead(g)^(degree(f) - degree(r))
      res *= R(-1)^(degree(g)*(degree(f) - degree(r)))
      f = r
      continue
    end

    if gcd(lift(lead(f)), m)  % gcd(lift(lead(g)), m) == 0
      q = divexact(lead(f), lead(g))
      ff = f - q*g*gen(Rt)^(degree(f) - degree(g))
      @assert degree(f) > degree(ff)
      res *= lead(g)^(degree(f) - degree(ff))
      res *= R(-1)^(degree(g)*(degree(f) - degree(ff)))
      f = ff
      continue
    end
    break
  end

  #factoring case, need to split the ring as well.
  #merde : we need a coprime factorisation: take
  # 6t^2+2t+3 mod 36....
  if easy
    cp = [m]
  else
    cp = [gcd(lift(coeff(g, i)), m) for i=0:degree(g)]
    push!(cp, m)
    cp = [x for x = cp if !iszero(x)]
    cp = coprime_base(cp)
    cp = [x for x = cp if !isunit(x)] #error: [1, 1, 3, 27] -> [1,3]
  end
  resp = fmpz[]
  pg = fmpz[]
  for p = cp
    lg = p^valuation(m, p)
    push!(pg, lg)

    if lg != m
      R1 = ResidueRing(FlintZZ, S(lg), cached=false)
      R1t = PolynomialRing(R1, cached=false)[1]
      #g is bad in R1, so factor it
      gR1 = R1t(lift(Zx, g))
      fR1 = R1t(lift(Zx, f))
    else
      gR1 = g
      fR1 = f
      R1 = R
      R1t = Rt
    end

    if degree(fR1) < degree(f) &&
       degree(gR1) < degree(g)
       res1 = R1(0)
    elseif degree(fR1) < degree(f)
        res1 = R1(-1)^(degree(g) * (degree(f) - degree(fR1))) *
               lead(gR1)^(degree(f) - degree(fR1))
    elseif degree(gR1) < degree(g)
        res1 = lead(fR1)^(degree(g) - degree(gR1))
    else
        res1 = R1(1)
    end

    if !isunit(lead(gR1))
      g1, g2 = fun_factor(gR1)

      #careful: lead(g) = 0 mod lg is possible, so the degree drops!
      #from Wiki:
      # phi: R -> S
      # phi(res(f, g)) = res(phi(f), phi(g)) if the degrees are the same
      #                = 0                   if both degrees drop (zero column in
      #                                         sylvester matrix)
      #                = (-1)^(delta(f)*deg(g))*
      #                  phi(l(g))^delta(f)  if only degree(f) drops
      #                = phi(l(f))^delta(g)  if only degree(g) drops
      #
      # we use res(f, gh) = res(f, g)res(f, h) which seems to be true in general
      #next: res(fR1, gR1) = res(phi(f), g1)
      #g1 is = 1 mod p
      #however, reverse(fR1) can have a smaller degree then fR1 (eg. x^2 -> 1)
      # res(f, g) = t(g)^delta(f) * (-1)^(deg(g)*delta(f)) * res(rev(f), rev(g))
      # there is a (-1)^deg(f)*deg(g) from Wiki for the reverse on top.

      res1 *= resultant_sircana(reverse(fR1), reverse(g1))
      v = valuation(fR1, gen(parent(fR1)))
           #should be delta(f) = degree(rev(f)) - degree(f)
      res1 *= constant_coefficient(g1)^v*R1(-1)^(v*degree(g1))
      res1 *= R1(-1)^(degree(g1)*degree(fR1))

      res1 *= resultant_sircana(fR1, g2)
      push!(resp, lift(res1))
    else
      #gR1 has a invertible leading coeff
      res1 *= resultant_sircana(fR1, gR1)
      push!(resp, lift(res1))
    end
  end
  res *= length(cp)==1 ? R(resp[1]) : R(crt(resp, pg))
  return res
end

#key idea (Carlo): if g = ab and a is a unit mod p, then it is actually a unit
# in Z/p^kZ, hence the ideal (f, g) = (f, b) where b is now monic.
#Thus rres(f,g ) = rres(f, b).... and the division can continue
doc"""
    rres(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer} -> T
> The reduced resultant of $f$ and $g$ using a quadratic-time algorithm.
> That is a generator for the $(f, g) \cap Z$
"""
function rres(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  return rres_sircana(f, g)
end

doc"""
    rres(f::fmpz_poly, g::fmpz_poly) -> fmpz
> The reduced resultant of $f$ and $g$,
> that is a generator for the ideal $(f, g) \cap Z$
"""
function rres(f::fmpz_poly, g::fmpz_poly)
  return rres_bez(f,g)
end

function rres_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Rt = parent(f)
  R = base_ring(Rt)
  m = fmpz(modulus(R))
  e, p = ispower(m)
  easy = isprime(p)

  Zx = PolynomialRing(FlintZZ)[1]

  while true
    if degree(f) < 1 && degree(g) < 1
      if iszero(f) || iszero(g)
        res = R(0)
      else
        res = R(gcd(lift(lead(f)), lift(lead(g))))
      end
      return res
    end

    if degree(f) < degree(g)
      f, g = g, f
    end

    if degree(g) < 1
      if !isunit(lead(f))
        #need the constant coeff * the annihilator of the others...
        a = coeff(f, 1)
        for i = 2:degree(f)
          a = gcd(a, coeff(f, i))
          if isone(a)
            break
          end
        end
        a = annihilator(a)
        if iszero(a)
          return lead(g)
        else
          return gcd(lead(g), a*constant_coefficient(f))
        end
      else
        return gcd(R(0), constant_coefficient(g))
      end
    end

    if !isunit(lead(g))
      c, g = primsplit(g)
      if easy
        cp = [m]
      else
        cp = [gcd(lift(coeff(g, i)), m) for i=0:degree(g)]
        push!(cp, m)
        cp = [x for x = cp if !iszero(x)]
        cp = coprime_base(cp)
        cp = [x for x = cp if !isunit(x)] #error: [1, 1, 3, 27] -> [1,3]
      end
      resp = fmpz[]
      pg = fmpz[]
      for p = cp
        lg = p^valuation(m, p)
        push!(pg, lg)

        if lg != m
          R1 = ResidueRing(FlintZZ, S(lg), cached=false)
          R1t = PolynomialRing(R1, cached=false)[1]
          #g is bad in R1, so factor it
          gR1 = R1t(lift(Zx, g))
          fR1 = R1t(lift(Zx, f))
        else
          gR1 = g
          fR1 = f
          R1 = R
          R1t = Rt
        end
        if isunit(lead(gR1))
          g2 = gR1
        else
          g1, g2 = fun_factor(gR1)
        end
        push!(resp, lift(lift(c)*rres_sircana(fR1, g2)))
      end
      res = length(cp)==1 ? R(resp[1]) : R(crt(resp, pg))
      return gcd(R(0), res)
    end

    q, f = divrem(f, g)
  end
end

doc"""
    rresx(f::PolyElem{ResElem{fmpz}}, g::PolyElem{ResElem{fmpz}}) -> r, u, v
> The reduced resultant $r$ and polynomials $u$ and $v$ s.th.
> $r = uf+vg$ and $\langle r\rangle = \langle f, g\rangle\cap \mathbb Z$.
"""
function rresx(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  return rresx_sircana(f, g)
end

function rresx_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  @assert isunit(lead(f)) || isunit(lead(g)) #can be weakened to invertable lead
  res::T, u::PolyElem{T}, v::PolyElem{T} = _rresx_sircana(f, g)::Tuple{T, PolyElem{T}, PolyElem{T}}
  if !iszero(res)
    cu = canonical_unit(res)
    cu = inv(cu)
    res *= cu
    u *= cu
    v *= cu
  end
  if isunit(lead(g))
    q::PolyElem{T}, r::PolyElem{T} = divrem(u, g)
    return res, r, v+q*f
  else
    q, r = divrem(v, f)
    return res, u+q*g, r
  end
end

function _rresx_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)

  Rt = parent(f)
  R = base_ring(Rt)
  Zx = FlintZZ["x"][1]
  m::fmpz = fmpz(modulus(R))
  e, p::fmpz = ispower(m)
  easy = isprime(p)

  u::PolyElem{T}, v::PolyElem{T} = Rt(0), Rt(1)
  U::PolyElem{T}, V::PolyElem{T} = Rt(1), Rt(0)
#  res::T, uu::T, vv::T

#  u, v = Rt(0), Rt(1)  #g = u f_in + v g_in
#  U, V = Rt(1), Rt(0)  #f = U f_in + V g_in

  Zx = PolynomialRing(FlintZZ)[1]

  while true
    if degree(f) < 1 && degree(g) < 1
      if iszero(f) || iszero(g)
        res = R(0)
        u = Rt(0)
        v = Rt(0)
      else
        res, uu, vv = gcdx(lead(f), lead(g))
        #res = uu*f + vv*g = uu*(U f_in + V g_in) + vv*(u f_in + v g_in)
        #    = uu*U + vv*u  uu*V + vv*v
        u, v = (uu*U + vv*u), (uu*V + vv*v)
      end
      return res, u, v
    end

    if degree(f) < degree(g)
      f, g = g, f
      U, u = u, U
      V, v = v, V
    end

    if degree(g) < 1
      if !isunit(lead(f))
        #need the constant coeff * the annihilator of the others...
        a = coeff(f, 1)
        for i = 2:degree(f)
          a = gcd(a, coeff(f, i))
          if isone(a)
            break
          end
        end
        a = annihilator(a)
        if iszero(a)
          return lead(g), u, v
        else
          res, uu, vv = gcdx(a*constant_coefficient(f), lead(g))
          u, v = (uu*a*U + vv*u), (uu*a*V + vv*v)

          return res, u, v
        end
      else
        return constant_coefficient(g), u, v
      end
    end

    if !isunit(lead(g))
      c, g = primsplit(g)
      if easy
        cp = [m]
      else
        cp = [gcd(lift(coeff(g, i)), m) for i=0:degree(g)]
        push!(cp, m)
        cp = [x for x = cp if !iszero(x)]
        cp = coprime_base(cp)
        cp = [x for x = cp if !isunit(x)] #error: [1, 1, 3, 27] -> [1,3]
      end
      resp = fmpz[]
      resB = []
      pg = fmpz[]
      for p = cp
        lg = p^valuation(m, p)
        push!(pg, lg)

        #gR1::PolyElem{T}, fR1::PolyElem{T}
        #c::T
        #g1::PolyElem{T}, g2::PolyElem{T}
        if lg != m
          R1 = ResidueRing(FlintZZ, S(lg), cached=false)
          R1t = PolynomialRing(R1, cached=false)[1]
          #g is bad in R1, so factor it
          gR1 = R1t(lift(Zx, g))
          fR1 = R1t(lift(Zx, f))
        else
          gR1 = g
          fR1 = f
          R1 = R
          R1t = Rt
        end
        if iszero(R1(lift(c)))
          push!(resp, fmpz(0))
          push!(resB, (R1t(0), R1t(1))) #relation need to be primitive
        else
          if isunit(lead(gR1))
            g2 = gR1
            g1 = R1t(1)
          else
            g1, g2 = fun_factor(gR1)
          end
          @assert isunit(g1)
          rr, uu, vv = rresx_sircana(fR1, g2)
          push!(resp, lift(lift(c)*rr))
          push!(resB, (uu*lift(c), inv(g1)*vv))
        end
      end
      if length(cp) == 1
        res, u_, v_ = R(resp[1]), Rt(resB[1][1]), Rt(resB[1][2])
      else
        ce = crt_env(pg)
        res = R(crt(resp, ce))
        u_ = Rt(induce_crt([x[1] for x = resB], ce))
        v_ = Rt(induce_crt([x[2] for x = resB], ce))
      end
      # f = U*f_in + V*g_in
      # g = u*f_in + v*g_in
      # r = u_ * f + v_ * g

      (u_*U + v_*u), (u_*V + v_*v)
      return res, (u_*U + v_*u), (u_*V + v_*v)
    end

    q, f = divrem(f, g)
    #f -> f-qg, so U*f_in + V * g_in -> ...
    U = U - q*u
    V = V - q*v
  end
end

function my_divrem(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  @show f, g
  @show g1, g2 = fun_factor(g)
  if degree(g2) < 1 # g2 is monic, so it's 1
    return parent(f)(0), g2
  end
  u = invmod(g1, g2)
  q, r = divrem(f, g2)
  q, r = divrem(r*u, g2)
  return g1*r, g2
end

#polynomial remainder sequence - almost
function prs_sircana(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Rt = parent(f)
  R = base_ring(Rt)
  m = fmpz(modulus(R))
  e, p = ispower(m)
  easy = isprime(p)
  @assert easy

  Zx = PolynomialRing(FlintZZ)[1]

  rs = []

  while degree(g) >0
    g *= inv(canonical_unit(lead(g)))
    c, gg = primsplit(g)
    @show f, (g, mu) = gg, my_divrem(f, gg)
    push!(rs, (c, gg, mu))
  end
  return rs, g
end

function Nemo.gcd(a::ResElem{T}, b::ResElem{T}) where T <: Union{Integer, fmpz}
  m = modulus(a)
  return parent(a)(gcd(gcd(a.data, m), b.data))
end

function Nemo.gcdx(a::ResElem{T}, b::ResElem{T}) where T <: Union{Integer, fmpz}
  m = modulus(a)
  R = parent(a)
  g, u, v = gcdx(fmpz(a.data), fmpz(b.data))
  G, U, V = gcdx(g, fmpz(m))
  return R(G), R(U)*R(u), R(U)*R(v)
end

doc"""
    annihilator(a::ResElem{fmpz}) -> r
    annihilator(a::ResElem{Integer}) -> r
> The annihilator of $a$, ie. a generator for the annihilator ideal
> $\{x | xa = 0\} = \langle r\rangle$
"""
function annihilator(a::ResElem{T}) where T <: Union{Integer, fmpz}
  R = parent(a)
  m = modulus(R)
  return R(divexact(m, gcd(m, a.data)))
end

doc"""
    isunit(f::Union{fmpz_mod_poly,nmod_poly}) -> Bool
> Tests if $f$ is a unit in the polynomial ring, ie. if
> $f = u + n$ where $u$ is a unit in the coeff. ring
> and $n$ is nilpotent.
"""
function Nemo.isunit(f::Union{fmpz_mod_poly,nmod_poly})
  if !isunit(constant_coefficient(f))
    return false
  end
  for i=1:degree(f)
    if !isnilpotent(coeff(f, i))
      return false
    end
  end
  return true
end

doc"""
    isnilpotent(a::ResElem{fmpz}) -> Bool
    isnilpotent(a::ResElem{Integer}) -> Bool
> Tests iff $a$ is nilpotent.
"""
function isnilpotent(a::ResElem{T}) where T <: Union{Integer, fmpz}
  #a is nilpontent if it is divisible by all primes divising the modulus
  # the largest exponent a prime can divide is nbits(m)
  l = nbits(modulus(a))
  return iszero(a^l)
end

function iszerodivisor(f::Union{fmpz_mod_poly,nmod_poly})
  c = content(f)
  return isnilpotent(c)
end

function Nemo.inv(f::Union{fmpz_mod_poly,nmod_poly})
  if !isunit(f)
    error("impossible inverse")
  end
  g = parent(f)(inv(trailing_coefficient(f)))
  #lifting: to invert a, start with an inverse b mod m, then
  # then b -> b*(2-ab) is an inverse mod m^2
  # starting with this g, and using the fact that all coeffs are nilpotent
  # we have an invers modulo s.th. nilpotent. Hence it works
  c = f*g
  while !isone(c)
    g = g*(2-c)
    c = f*g
  end
  return g
end

function Nemo.invmod(f::T, M::T) where T <: Union{fmpz_mod_poly,nmod_poly}
  if !isunit(f)
    error("impossible inverse")
  end
  if !isunit(lead(M))
    error("modulus must be monic")
  end
  g = parent(f)(inv(trailing_coefficient(f)))
  #lifting: to invert a, start with an inverse b mod m, then
  # then b -> b*(2-ab) is an inverse mod m^2
  # starting with this g, and using the fact that all coeffs are nilpotent
  # we have an invers modulo s.th. nilpotent. Hence it works
  c = f*g
  q, c = divrem(c, M)
  while !isone(c)
    g = g*(2-c)
    q, g = divrem(g, M)
    c = f*g
    q, c = divrem(c, M)
  end
  return g
end

#for testing: the obvious(?) naive method(s)
function rres_hnf(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  R = base_ring(f)
  Zx = FlintZZ["x"][1]
  s = Nemo.Generic.sylvester_matrix(lift(Zx, f), lift(Zx, g))
  h = hnf(s)
  return gcd(R(0), R(h[rows(h), cols(h)]))
end

function rres_bez(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  R = base_ring(f)
  Zx = FlintZZ["x"][1]
  Qx = FlintQQ["x"][1]
  g, q, w = gcdx(Qx(lift(Zx, f)), Qx(lift(Zx, g)))
  return gcd(R(0), R(lcm(denominator(q), denominator(w))))
end

function rres_hnf(f::fmpz_poly, g::fmpz_poly)
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  s = Nemo.Generic.sylvester_matrix(f, g)
  h = hnf(s)
  return abs(h[rows(h), cols(h)])
end

function rres_bez(f::fmpz_poly, g::fmpz_poly)
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Qx = FlintQQ["x"][1]
  g, q, w = gcdx(Qx(f), Qx(g))
  return lcm(denominator(q), denominator(w))
end

doc"""
    rresx(f::fmpz_poly, g::fmpz_poly) -> r, u, v
> The reduced resultant, ie. a generator for the intersection
> of the ideal generated by $f$ and $g$ with the integers.
> As well as polynomials $u$ and $v$ s.th. $r = uf+vg$,
> $\deg u < \deg g$ and $\deg v < \deg f$.
"""
function rresx(f::fmpz_poly, g::fmpz_poly)
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Qx = FlintQQ["x"][1]
  g, q, w = gcdx(Qx(f), Qx(g))
  l = lcm(denominator(q), denominator(w))
  Zx = parent(f)
  return l, Zx(l*q), Zx(l*w)
end

doc"""
    xxgcd(a::ResElem{fmpz}, b::ResElem{fmpz}) -> g, e, f, u, v
    xxgcd(a::ResElem{Integer}, b::ResElem{Integer}) -> g, e, f, u, v
> Coputes $g = \gcd(a, b)$, the Bezout coefficients, $e$, $f$ s.th.
> $g = ea+fb$ and $u$, $v$ s.th. $ev-fu = 1$, $gu = a$ and $gv = b$
"""
function xxgcd(a::ResElem{S}, b::ResElem{S}) where S <: Union{fmpz, Integer}
  g, e, f = gcdx(a, b)
  if iszero(g)
    return g, e, f, f, e
  else
    return g, e, f, divexact(a, g), divexact(b, g)
  end
end

doc"""
    primsplit!(f::PolyElem{ResElem{fmpz}}) -> c, f
    primsplit!(f::PolyElem{ResElem{Integer}}) -> c, f
> Computes the contents $c$ and the primitive part of $f$ destructively.
"""
function primsplit!(f::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}

  d = degree(f)
  if d == 0
    if iszero(f)
      return base_ring(parent(f))(1), f
    end
    c = canonical_unit(coeff(f, 0))
    setcoeff!(f, 0, 1)
    return inv(c)*coeff(f, 0), f
  end

  g = coeff(f, 0)
  setcoeff!(f, 0, 1)
  for i=1:d
    h, _, _, u, v = xxgcd(g, coeff(f, i))
    setcoeff!(f, i, v)
    if  g != h
      for j=0:i-1
        setcoeff!(f, j, u*coeff(f, j))
      end
    end
    g = h
    if isone(g)
      return g, f
    end
  end
  return g, f
end

doc"""
    primsplit(f::PolyElem{ResElem{fmpz}}}) -> c, f
    primsplit(f::PolyElem{ResElem{Integer}}}) -> c, f
> Computes the contents $c$ and the primitive part of $f$.
"""
function primsplit(f::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  g = deepcopy(f)
  return primsplit!(g)
end

immutable PolyCoeffs{T <: PolyElem}
  f::T
end

function coefficients(f::PolyElem)
  return PolyCoeffs(f)
end

function Base.start(a::PolyCoeffs)
  return 0
end

function Base.next(a::PolyCoeffs, b::Int)
  return coeff(a.f, b), b+1
end

function Base.done(a::PolyCoeffs, b::Int)
  return b > degree(a.f)
end

function Base.length(a::PolyCoeffs)
  return length(a.f)
end

function Base.endof(a::PolyCoeffs)
  return degree(a.f)
end

function Base.getindex(a::PolyCoeffs, i::Int)
  return coeff(a.f, i)
end

doc"""
    resultant_valuation(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer} -> T
> A generator for the ideal of the resultant of $f$ anf $g$ using a quadratic-time algorithm.
> One of the two polynomials must be monic.
"""
function resultant_valuation(f::PolyElem{T}, g::PolyElem{T}) where T <: ResElem{S} where S <: Union{fmpz, Integer}
  #The algorithm is the same as the resultant. We assume that one fo the 2 polynomials is monic. Under this assumption, at every
  #step the same is true and we can discard the unti obtained from the fun_factor function
  Nemo.check_parent(f, g)
  @assert typeof(f) == typeof(g)
  Rt = parent(f)
  R = base_ring(Rt)
  m = fmpz(modulus(R))
  e, p = ispower(m)
  easy = isprime(p)

  Zx = PolynomialRing(FlintZZ)[1]

  res = R(1)

  while true
    if degree(f) < 1 && degree(g) < 1
      if iszero(f) || iszero(g)
        res *= R(0)
      end
      return res
    end

    if degree(f) < 1
      res *= lead(f)^degree(g)
      return res
    end

    if degree(g) < 1
      res *= lead(g)^degree(f)
      return res
    end

    c, f = primsplit(f)
    if !isone(c)
      res *= c^degree(g)
    end

    c, g = primsplit(g)
    if !isone(c)
      res *= c^degree(f)
    end

    if degree(f) < degree(g)
      f, g = g, f
    end

    #want f % g which works iff lead(g) | lead(f)

    if isunit(lead(g)) #accelerate a bit...possibly.
      f = rem(f, g)
      continue
    end

    break
  end

  #factoring case, need to split the ring as well.
  #merde : we need a coprime factorisation: take
  # 6t^2+2t+3 mod 36....
  if easy
    cp = [m]
  else
    cp = [gcd(lift(coeff(g, i)), m) for i=0:degree(g)]
    push!(cp, m)
    cp = [x for x = cp if !iszero(x)]
    cp = coprime_base(cp)
    cp = [x for x = cp if !isunit(x)] #error: [1, 1, 3, 27] -> [1,3]
  end
  resp = fmpz[]
  pg = fmpz[]
  for p = cp
    lg = p^valuation(m, p)
    push!(pg, lg)

    if lg != m
      R1 = ResidueRing(FlintZZ, S(lg), cached=false)
      R1t = PolynomialRing(R1, cached=false)[1]
      #g is bad in R1, so factor it
      gR1 = R1t(lift(Zx, g))
      fR1 = R1t(lift(Zx, f))
    else
      gR1 = g
      fR1 = f
      R1 = R
      R1t = Rt
    end

    if degree(fR1) < degree(f) && degree(gR1) < degree(g)
      res1 = R1(0)
    elseif degree(fR1) < degree(f)
      res1 = lead(gR1)^(degree(f) - degree(fR1))
    else
      res1 = R1(1)
    end

    if !isunit(lead(gR1))
      g1, g2 = fun_factor(gR1)
      res1 *= resultant_valuation(fR1, g2)
      push!(resp, lift(res1))
    else
      #gR1 has a invertible leading coeff
      res1 *= resultant_valuation(fR1, gR1)
      push!(resp, lift(res1))
    end
  end
  res *= length(cp)==1 ? R(resp[1]) : R(crt(resp, pg))
  return res
end

################################################################################
#
#  fmpq_poly with denominator 1 to fmpz_poly
#
################################################################################

function (a::FmpzPolyRing)(b::fmpq_poly)
  (!isone(denominator(b))) && error("Denominator has to be 1")
  z = a()
  ccall((:fmpq_poly_get_numerator, :libflint), Void,
              (Ref{fmpz_poly}, Ref{fmpq_poly}), z, b)
  return z
end

################################################################################
#
#  Power sums
#
################################################################################

doc"""
    polynomial_to_power_sums(f::PolyElem{T}, n::Int=degree(f)) -> Array{T, 1}
> Uses Newton (or Newton-Girard) formulas to compute the first $n$
> power sums from the coefficients of $f$.
"""
function polynomial_to_power_sums(f::PolyElem{T}, n::Int=degree(f)) where T <: FieldElem
  d = degree(f)
  R = base_ring(f)
  S = PowerSeriesRing(R, n+1, "gen(S)")[1]
  #careful: converting to power series and derivative do not commute
  #I also don't quite get this: I thought this was just the log,
  #but it isn't
  A = S([coeff(reverse(derivative(f)), i) for i=0:d-1], d, n+1, 0)
  B = S([coeff(reverse(f), i) for i=0:d], d+1, n+1, 0)
  L = A*inv(B)
  s = [coeff(L, i) for i=1:n]
  return s
end

#plain vanilla recursion
function polynomial_to_power_sums(f::PolyElem{T}, n::Int=degree(f)) where T
  d = degree(f)
  R = base_ring(f)
  E = T[(-1)^i*coeff(f, d-i) for i=0:min(d, n)] #should be the elementary symm.
  while length(E) <= n
    push!(E, R(0))
  end
  P = T[]

  push!(P, E[1+1])
  for k=2:n
    push!(P, (-1)^(k-1)*k*E[k+1] + sum((-1)^(k-1+i)*E[k-i+1]*P[i] for i=1:k-1))
  end
  return P
end

doc"""
    power_sums_to_polynomial(P::Array{T, 1}) -> PolyElem{T}
> Uses the Newton (or Newton-Girard) identities to obtain the polynomial
> coefficients (the elementary symmetric functions) from the power sums.
"""
function power_sums_to_polynomial(P::Array{T, 1}) where T <: FieldElem
  d = length(P)
  R = parent(P[1])
  S = PowerSeriesRing(R, d, "gen(S)")[1]
  s = S(P, length(P), d, 0)
  if !false
    r = - integral(s)
    r = exp(r)
  end

  if false
    r = S(T[R(1), -P[1]], 2, 2, 0)
    la = [d+1]
    while la[end]>1
      push!(la, div(la[end]+1, 2))
    end
    n = length(la)-1
    while n > 0
      set_prec!(r, la[n])
      rr = derivative(r)*inv(r)
      md = -(rr+s)
      m = S([R(1)], 1, la[n], 0)+integral(md)
      r *= m
      n -= 1
    end
    println("new exp $r")
  end
  d = Nemo.pol_length(r)
  v = valuation(r)
  @assert v==0
  while d>=0 && iszero(Nemo.polcoeff(r, d))
    d -= 1
  end
  return PolynomialRing(R, cached = false)[1]([Nemo.polcoeff(r, d-i) for i=0:d])
end

function power_sums_to_polynomial(P::Array{T, 1}) where T
  E = T[1]
  R = parent(P[1])
  last_non_zero = 0
  for k=1:length(P)
    push!(E, divexact(sum((-1)^(i-1)*E[k-i+1]*P[i] for i=1:k), R(k)))
    if E[end] != 0
      last_non_zero = k
    end
  end
  E = E[1:last_non_zero+1]
  d = length(E) #the length of the resulting poly...
  for i=1:div(d, 2)
    E[i], E[d-i+1] = (-1)^(d-i)*E[d-i+1], (-1)^(i-1)*E[i]
  end
  return PolynomialRing(R, cached = false)[1](E)
end

function factor(f::fmpq_poly, R::NmodRing)
  Rt, t = PolynomialRing(R, "t", cached=false)
  return factor(Rt(f))
end

function factor(f::fmpz_poly, R::NmodRing)
  Rt, t = PolynomialRing(R, "t", cached=false)
  return factor(Rt(f))
end

function roots(f::fmpq_poly, R::Nemo.FqNmodFiniteField)
  Rt, t = PolynomialRing(R, "t", cached=false)
  fp = FlintZZ["t"][1](f*denominator(f))
  fpp = Rt(fp)
  return roots(fpp)
end

function roots(f::fmpq_poly, R::Nemo.NmodRing)
  Rt, t = PolynomialRing(R, "t", cached=false)
  fp = FlintZZ["t"][1](f*denominator(f))
  fpp = Rt(fp)
  return roots(fpp)
end

function roots(f::T) where T <: Union{fq_nmod_poly, fq_poly} # should be in Nemo and
                                    # made available for all finite fields I guess.
  q = size(base_ring(f))
  x = gen(parent(f))
  if degree(f) < q
    x = powmod(x, q, f)-x
  else
    x = x^q-x
  end
  f = gcd(f, x)
  l = factor(f).fac
  return elem_type(base_ring(f))[-trailing_coefficient(x) for x = keys(l) if degree(x)==1]
end

function roots(f::PolyElem)
  lf = factor(f)
  return elem_type(base_ring(f))[-trailing_coefficient(x) for x= keys(lf.fac) if degree(x)==1]
end

function setcoeff!(z::fq_nmod_poly, n::Int, x::fmpz)
   ccall((:fq_nmod_poly_set_coeff_fmpz, :libflint), Void,
         (Ptr{fq_nmod_poly}, Int, Ptr{fmpz}, Ptr{FqNmodFiniteField}),
         &z, n, &x, &base_ring(parent(z)))
     return z
end


text:
using AlgebraPDF

@testset "Parameter utils" begin
    @test nt(:d) == (d=0.0, )
    @test nt(:d, 3) == (d=3, )
    @test nt(:d, (1.1,0.1)) == (d = (1.1,0.1), )
end

@testset "Operations on NamedTuple" begin
    @test (a=1,b=2) + (d=1,c=2) == (a=1,b=2,d=1,c=2)
    @test (a=1,b=2,d=1,c=2) - (d=1,) == (a=1,b=2,c=2)
    @test (a=1,b=2,d=1,c=2) - :d == (a=1,b=2,c=2)
    @test (a=1,b=2,d=1,c=2) - [:d, :a] == (b=2,c=2)
    @test (a=1,b=2,d=1,c=2) - (:d, :a) == (b=2,c=2)
end

@testset "Parameter structure: NamedTuple" begin
    p0 = (a=1,b=2)
    @test pars(p0) == p0
    @test freepars(p0) == p0
    @test fixedpars(p0) == ∅
    @test AlgebraPDF.updatevalueorflag(p0, :a, true, 2) == (a=2,b=2)
    @test_throws ArgumentError AlgebraPDF.updatevalueorflag(p0, :a, false)
end

@testset "Parameter structure: FlaggedNamedTuple" begin
    #
    p0 = FlaggedNamedTuple(a=1,b=2)
    @test p0.a == 1
    @test p0.b == 2
    #
    p1 = AlgebraPDF.updatevalueorflag(p0, :a, true, 2)
    @test pars(p1) == (a=2,b=2)
    #
    p2 = AlgebraPDF.updatevalueorflag(p0, :a, false)
    @test pars(p0) == pars(p2)
    @test freepars(p2) == (b=2,)
    @test fixedpars(p2) == (a=1,)
    #
    @test p0 == FlaggedNamedTuple((a=1,b=2))
    @test p2 == FlaggedNamedTuple(p2)
    #
    @test keys(p1) == keys(p2)
    @test keys(p2, true) == (:b,)
    @test keys(p2, false) == (:a,)
end

@testset "Addind parameter sets: FNT+FNT, FNT+NT, NT+FNT" begin
    p1 = FlaggedNamedTuple(a1=1,b1=2)
    p2 = FlaggedNamedTuple(a2=1,b2=2)
    p12 = p1 + p2
    @test pars(p12) == pars(p1) + pars(p2)

    p2′ = AlgebraPDF.updatevalueorflag(p2, :a2, false)
    p12′ = p1 + p2′
    @test pars(p12′, true) == pars(p1, true) + pars(p2′, true)
    @test pars(p12′, false) == pars(p1, false) + pars(p2′, false)

    pd1′′ = p12′ + (d4=2.2, )
    pd2′′ = (d4=2.2, ) + p12′
    @test :d4 ∈ keys(freepars(pd1′′))
    @test :d4 ∈ keys(freepars(pd2′′))
end


text:
# Test mapping arrays directly to aesthetics without a data fram.

using Gadfly
using DataFrames

plot(x=collect(1:100), y=sort(rand(100)))


text:
include("CoupledSimulations/clock.jl")
include("CouplerState/cplstate_interface.jl")
# include("CoupledSimulations/cplsolver.jl")

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
e = read_elp(joinpath(dirname(@__FILE__), "..", "data", "test.elp"))

@test length(e) == 2
@test e[1].label == "Fpz"
@test e[2].label == "Fp2"


text:
# Structured dataset coordinates can be specified using either a 4D array
# (3, Ni, Nj, Nk), or a tuple (x, y, z), or a 3D array of vectors (typically
# SVector{3}).
const Array4 = AbstractArray{T, 4} where T
const Array3Tuple3 = Tuple{Vararg{AbstractArray{T,3}, 3}} where T
const Array3ofVec3 = AbstractArray{T, 3} where T
const StructuredCoords = Union{Array4, Array3Tuple3, Array3ofVec3}

structured_dims(xyz::Array4) = ntuple(d -> size(xyz, d + 1), 3)
structured_dims(xyz::Array3Tuple3) = size(first(xyz))
structured_dims(xyz::Array3ofVec3) = size(xyz)

function vtk_grid(
        dtype::VTKStructuredGrid, filename::AbstractString,
        xyz::StructuredCoords;
        extent = structured_dims(xyz), whole_extent = extent, kwargs...,
    )
    Ns = structured_dims(xyz)
    Npts = prod(Ns)
    Ncomp = num_components(xyz, Npts)
    Ncls = num_cells_structured(Ns)

    if Ncomp != 3  # three components (x, y, z)
        # TODO the error message is only accurate when xyz isa Array4
        msg = "coordinate array `xyz` has incorrect dimensions.\n" *
              "Expected dimensions: $((3, Ns...)).\n" *
              "Actual dimensions: $((Ncomp, Ns...))"
        throw(DimensionMismatch(msg))
    end

    xvtk = XMLDocument()
    vtk = DatasetFile(dtype, xvtk, filename, Npts, Ncls; kwargs...)

    # VTKFile node
    xroot = vtk_xml_write_header(vtk)

    # StructuredGrid node
    xGrid = new_child(xroot, vtk.grid_type)
    let ext = extent_attribute(whole_extent)
        set_attribute(xGrid, "WholeExtent", ext)
    end

    # Piece node
    xPiece = new_child(xGrid, "Piece")
    let ext = extent_attribute(extent)
        set_attribute(xPiece, "Extent", ext)
    end

    # Points node
    xPoints = new_child(xPiece, "Points")

    # DataArray node
    data_to_xml(vtk, xPoints, xyz, "Points", 3)

    return vtk::DatasetFile
end

# 3D variant of vtk_grid with 4D array xyz.
vtk_grid(filename::AbstractString, xyz::AbstractArray{T,4};
        kwargs...) where {T <: Number} =
    vtk_grid(VTKStructuredGrid(), filename, xyz; kwargs...)

# 3D variant of vtk_grid with 3D arrays x, y, z.
function vtk_grid(filename::AbstractString, x::AbstractArray{T,3},
                  y::AbstractArray{T,3}, z::AbstractArray{T,3};
                  kwargs...) where {T <: Number}
    if !(size(x) == size(y) == size(z))
        throw(DimensionMismatch("size of x, y and z arrays must be the same."))
    end
    vtk_grid(VTKStructuredGrid(), filename, (x, y, z); kwargs...)
end

# 2D variant of vtk_grid with 3D array xy
function vtk_grid(filename::AbstractString, xy::AbstractArray{T,3};
            kwargs...) where {T <: Number}
    Ncomp, Ni, Nj = size(xy)
    if Ncomp != 2
        msg = "coordinate array `xy` has incorrect dimensions.\n" *
              "Expected dimensions: (2, Ni, Nj).\n" *
              "Actual dimensions: $(size(xy))"
        throw(DimensionMismatch(msg))
    end
    xyz = (
        reshape(view(xy, 1, :, :), Ni, Nj, 1),
        reshape(view(xy, 2, :, :), Ni, Nj, 1),
        Zeros{T}(Ni, Nj, 1),
    )
    vtk_grid(VTKStructuredGrid(), filename, xyz; kwargs...)
end

# 2D variant of vtk_grid with 2D arrays x, y.
function vtk_grid(filename::AbstractString, x::AbstractArray{T,2},
                  y::AbstractArray{T,2}; kwargs...) where {T <: Number}
    if size(x) != size(y)
        throw(DimensionMismatch("size of x and y arrays must be the same."))
    end
    Ni, Nj = size(x)
    xyz = (
        reshape(x, Ni, Nj, 1),
        reshape(y, Ni, Nj, 1),
        Zeros{T}(Ni, Nj, 1),
    )
    vtk_grid(VTKStructuredGrid(), filename, xyz; kwargs...)
end

# Case in which !(T <: Number). In this case, one typically has T <: SVector{3}.
function vtk_grid(filename::AbstractString, xyz::AbstractArray; kwargs...)
    eltype(xyz) <: Number && throw(ArgumentError("expected element type to be SVector or similar"))
    N = ndims(xyz)
    points = if N === 3
        xyz
    elseif N === 2
        reshape(xyz, size(xyz)..., 1)
    elseif N === 1
        reshape(xyz, size(xyz)..., 1, 1)
    end
    vtk_grid(VTKStructuredGrid(), filename, points; kwargs...)
end


text:
"""
Spikes

Includes types and methods for handling spike trains.
"""
module Spikes

using JLD
import Base: show, isequal, hash
include("metadata.jl")

################################################################################
#### Include appropriate files
################################################################################
include("SpikeTrains.jl")
include("poissonproc.jl")
include("fileio.jl")
# include("CRCNS/CRCNS.jl")

################################################################################
#### Miscellaneous functions and constants
################################################################################

export show, isequal, hash, savespikes, loadspikes,
       SpikeTrains, histogram, raster, n_cells,
       inhomogeneous_poisson_process

end #module


text:
module ParametricMachinesDemos

using Base: front, setindex
using NNlib: conv!, ∇conv_filter!, DenseConvDims
using LinearAlgebra: mul!
using Functors: @functor

using ChainRulesCore: NoTangent, unthunk, derivatives_given_output, ChainRulesCore

import ChainRules

export DenseMachine, ConvMachine, RecurMachine

include("slice.jl")
include("solve.jl")
include("architectures.jl")

end