text:
@everywhere struct Constants
    C::Int
    P::Int
    J::Int
    Q::Int
    H::Int
    tt::Int
    delta::Int
end

@everywhere mutable struct Timer
    period::Int
    horizon_plan::Int
    available_cranes::Array{Int, 1}
end

@everywhere struct LTask
    p::Int
    b::Int
    c::Int
    t::Int
end

@everywhere mutable struct LoadingSequence
    order::Array{NamedTuple{(:task, :start_time, :qc),Tuple{LTask,Int,Int}}, 1}
    tasks_left::Int
    len::Int
    filled_pos::Array{Int, 1}
    loaded_cont::Array{Int, 1}
end

@everywhere mutable struct QuayCrane
    q::Int
    available_bays::Array{Int, 1}
    status::String
    time_left::Int
    current_bay::Number
    next_bay::Int
    task_buffer::Array{LTask, 1}
end

@everywhere function init_timer(CTS::Constants)
    TIME = Timer(0, CTS.H, collect(1:CTS.Q))
    return(TIME)
end

@everywhere function init_ls(CTS::Constants)
    LS = LoadingSequence([], CTS.P, 0, [], [])
    return(LS)
end

@everywhere function init_bay(q::Int, CTS::Constants)
    if q == 1
        return(1)
    elseif  q == CTS.Q
        return(CTS.J)
    else
        return(floor(CTS.J/(CTS.Q-1)*(q-1)))
    end
end

@everywhere function init_qc(CTS::Constants)
    QC=Array{QuayCrane, 1}()
    for q = 1:CTS.Q
        qc_pos = subset_bay(CTS, q)
        push!(QC, QuayCrane(q, qc_pos, "idle", 0, init_bay(q, CTS), init_bay(q, CTS), Array{LTask, 1}()))
    end
    return(QC)
end

@everywhere function Base.show(io::IO, LS::LoadingSequence)
    println(io, "LOADING SEQUENCE")
    println(io, "- tasks_left: ", LS.tasks_left)
    println(io, "- filled_pos: ", LS.filled_pos)
    println(io, "- loaded_cont: ", LS.loaded_cont)
    for i = 1:length(LS.order)
        print(io, "  "*string(i)*". ")
        println(io, LS.order[i])
    end
end

function Base.show(io::IO, ::MIME"text/plain", QC::QuayCrane)
    println(io, "QUAY CRANE - "*string(QC.q))
    println(io, " - status: ", QC.status)
    println(io, " - time_left: ", QC.time_left)
    println(io, " - current_bay: ", QC.current_bay)
    println(io, " - next_bay: ", QC.next_bay)
    println(io, " - task_buffer: ", QC.task_buffer)
    println(io, " - available_bays: ", QC.available_bays)
    println("\n")
end

function Base.show(io::IO, QC::QuayCrane)
    println(io, string(QC.q)*". QUAY CRANE")
    println(io, " - status: ", QC.status)
    println(io, " - time_left: ", QC.time_left)
    println(io, " - current_bay: ", QC.current_bay)
    println(io, " - next_bay: ", QC.next_bay)
    println(io, " - task_buffer: ", QC.task_buffer)
    println(io, " - available_bays: ", QC.available_bays)
    println("\n")
end

function Base.show(io::IO, TIME::Timer)
    println("--------------")
    println(io, "TIME PERIOD: ", TIME.period)
    println(io, " - available cranes: ", TIME.available_cranes)
    println("--------------")
end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
"""
    IS()

Importance sampling algorithm.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

Usage:

```julia
IS()
```

Example:

```julia
# Define a simple Normal model with unknown mean and variance.
@model gdemo(x) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt.(s))
    x[1] ~ Normal(m, sqrt.(s))
    x[2] ~ Normal(m, sqrt.(s))
    return s, m
end

sample(gdemo([1.5, 2]), IS(), 1000)
```
"""
struct IS{space} <: InferenceAlgorithm end

IS() = IS{()}()

mutable struct ISState{V<:VarInfo, F<:AbstractFloat} <: AbstractSamplerState
    vi                 ::  V
    final_logevidence  ::  F
end

ISState(model::Model) = ISState(VarInfo(model), 0.0)

function Sampler(alg::IS, model::Model, s::Selector)
    info = Dict{Symbol, Any}()
    state = ISState(model)
    return Sampler(alg, info, s, state)
end

function AbstractMCMC.step!(
    ::AbstractRNG,
    model::Model,
    spl::Sampler{<:IS},
    ::Integer,
    transition;
    kwargs...
)
    empty!(spl.state.vi)
    model(spl.state.vi, spl)

    return Transition(spl)
end

function AbstractMCMC.sample_end!(
    ::AbstractRNG,
    ::Model,
    spl::Sampler{<:IS},
    N::Integer,
    ts::Vector;
    kwargs...
)
    # Calculate evidence.
    spl.state.final_logevidence = logsumexp(map(x->x.lp, ts)) - log(N)
end

function DynamicPPL.assume(spl::Sampler{<:IS}, dist::Distribution, vn::VarName, vi)
    r = rand(dist)
    push!(vi, vn, r, dist, spl)
    return r, 0
end

function DynamicPPL.observe(spl::Sampler{<:IS}, dist::Distribution, value, vi)
    # acclogp!(vi, logpdf(dist, value))
    return logpdf(dist, value)
end

ext: jl


text:
using ForwardDiff

println("Testing Partials...")
t = @elapsed include("PartialsTest.jl")
println("done (took $t seconds).")

println("Testing Dual...")
t = @elapsed include("DualTest.jl")
println("done (took $t seconds).")

println("Testing derivative functionality...")
t = @elapsed include("DerivativeTest.jl")
println("done (took $t seconds).")

println("Testing gradient functionality...")
t = @elapsed include("GradientTest.jl")
println("done (took $t seconds).")

println("Testing jacobian functionality...")
t = @elapsed include("JacobianTest.jl")
println("done (took $t seconds).")

println("Testing hessian functionality...")
t = @elapsed include("HessianTest.jl")
println("done (took $t seconds).")

println("Testing perturbation confusion functionality...")
t = @elapsed include("ConfusionTest.jl")
println("done (took $t seconds).")

println("Testing miscellaneous functionality...")
t = @elapsed include("MiscTest.jl")
println("done (took $t seconds).")


text:
#=
    run_tests
    Copyright © 2019 Mark Wells <mwellsa@gmail.com>

    Distributed under terms of the MIT license.
=#

using InterpolatedRejectionSampling
using Test

import Random.seed!
seed!(1234)

@testset "irsample" begin

    X = range(0,π,length=5)
    Y = range(0,π/4,length=4)
    knots = (X,Y)
    prob = [sin(x)+tan(y) for x in X, y in Y]

    n = 100
    xy = irsample(knots, prob, n)
    @test isa(xy, Matrix{Float64})
    @test size(xy) == (2,n)

    xy = Matrix{Union{Float64,Missing}}(missing,2,n)
    irsample!(xy, knots, prob)
    @test isa(xy, Matrix{Union{Missing,Float64}})
    @test size(xy) == (2,n)
    @test iszero(count(ismissing.(xy)))

    x = irsample(X, sin.(X), n)
    @test isa(x, Vector{Float64})
    @test length(x) == n
end


text:
# using ComplexModeler, Modeler, Poly,
using LightXML
# using ComplexModeler
ROOT=pwd()
include(joinpath(ROOT,"..","src_PowSysMod", "PowSysMod_body.jl"))

filenames = sort([joinpath("instances", "IIDM", "etat_reseau", file) for file in readdir(joinpath("instances", "IIDM", "etat_reseau"))])

@time OPFpbs = build_OPFproblems(IIDMeodInput, filenames)

OPFpbs["BaseCase"].ds.bus["BUS_8"]

scenario = "BaseCase"
pb = build_Problem!(OPFpbs, scenario)
# print(pb)

voltage = Point()
for (busname, buselems) in OPFpbs["BaseCase"].ds.bus
  if !haskey(buselems, "VoltMeasure")
    warn("$busname has no volt measure")
  else
    voltage[OPFpbs["BaseCase"].mp.node_vars[busname]["Volt"]] = buselems["VoltMeasure"].V
  end
end

## Testing nodal power balance
println("$(length(pb.constraints)) constraints for initial OPF")
for cstrname in keys(pb.constraints)
  if ismatch(r"Cstr_", cstrname)
    delete!(pb.constraints, cstrname)
  end
end
println("$(length(pb.constraints)) constraints for power balance only OPF")

cstrs = get_constraints(pb, voltage)

slacks = get_slacks(pb, voltage)
println("minslack : ", get_minslack(pb, voltage))

ordered_slacks = sort([ [k abs(v)] for (k, v) in slacks], by=x->abs(x[2]))
# sort(collect(keys(slacks)), by=x->x[1], rev=true)
for (k, v) in ordered_slacks
  @printf("→ %18s  %.0e + %.0eim   %.0e + %.0e <  %.0e+%.0e im  < %.0e + %.0e\n", k, real(slacks[k]), imag(slacks[k]), real(pb.constraints[k[1]].lb), imag(pb.constraints[k[1]].lb), real(cstrs[k]), imag(cstrs[k]), real(pb.constraints[k[1]].ub), imag(pb.constraints[k[1]].ub))
  # println("→ $k   \t $(slacks[k]) \t$(pb.constraints[k[1]].lb) < \t $(cstrs[k]) < \t$(pb.constraints[k[1]].ub)")
end

## Power at orig and dest of lines.
# line_power = Dict{Link, Array{Complex128}}()
# for link in get_links(OPFpbs, "BaseCase")
#   link_elem = OPFpbs["BaseCase"].ds.link[link]
#   length(link_elem) == 1 || warn("link $link has several lines.")
#   for (line_name, line) in link_elem
#     line_power[link] = [line.S_orig line.S_dest]
#   end
# end

# ## Testing nodal power balance
# ds = OPFpbs["BaseCase"].ds
# gs = OPFpbs["BaseCase"].gs
# mc = OPFpbs["BaseCase"].mc
#
# pb_nodal_balance = Problem()
# for (busname, bus_elems) in ds.bus
#   bus_elems_formulations = mc.node_formulations[busname]
#   bus_elems_var = gs.node_vars[busname]
#
#   cstrname, Snode, lb, ub = get_Snodal(busname, bus_elems, bus_elems_formulations, bus_elems_var)
#   S_balance = Snode + get_Slinks_in(busname, ds, gs, mc) + get_Slinks_out(busname, ds, gs, mc)
#
#   if cstrname == "_"
#     cstrname = "LOAD_"
#   end
#   add_constraint!(pb_nodal_balance, cstrname*string(busname), lb << S_balance << ub)
# end
#
# get_constraints(pb_nodal_balance, voltage)
# slacks = get_slacks(pb_nodal_balance, voltage)
# minslack, _ = get_minslack(pb_nodal_balance, voltage)
# println("minslack on nodal powerbalance: $minslack")
# # for k in sort(collect(keys(slacks)), by=x->x[1], rev=true)
# #   println("→ $k   \t $(slacks[k]) \t$(pb_nodal_balance.constraints[k[1]].lb) < \t $(cstrs[k]) < \t$(pb_nodal_balance.constraints[k[1]].ub)")
# # end

# ## Testing line power balance
# line_power = Point()
# pb_line_balance = Problem()
# for link in get_links(OPFpbs, "BaseCase")
#   for (elementname, elem) in ds.link[link]
#     # elementname = collect(keys(ds.link[link]))[1]
#     # element = ds.link[link][elementname]
#     elem_formulation = OPFpbs["BaseCase"].mc.link_formulations[link][elementname]
#     link_vars = OPFpbs["BaseCase"].gs.link_vars[link]
#     Sor = Sorig(element, link, elementname, elem_formulation, link_vars)
#     Sde = Sdest(element, link, elementname, elem_formulation, link_vars)
#
#     add_constraint!("", element.S_orig << Sor << element.S_orig)
#     evaluate(Sor, voltage)
#     evaluate(Sde, voltage)
#
#     element.S_dest
#   end
# end

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
mutable struct Script
    instructions::Array{Union{UInt8, Array{UInt8, 1}}, 1}
    Script(instructions::Nothing) = new(Union{UInt8, Array{UInt8, 1}}[])
    Script(instructions) = new(instructions)
end

function show(io::IO, z::Script)
    for instruction in z.instructions
        if typeof(instruction) <: Integer
            if haskey(OP_CODE_NAMES, instruction)
                print(io, "\n", OP_CODE_NAMES[Int(instruction)])
            else
                print(io, "\n", string("OP_CODE_", Int(instruction)))
            end
        else
            print(io, "\n", bytes2hex(instruction))
        end
    end
end

"""
    scriptparse(::GenericIOBuffer) -> Script

Returns a Script object from an IOBuffer
"""
function scriptparse(s::Base.GenericIOBuffer{Array{UInt8,1}})
    length_ = read_varint(s)
    instructions = []
    count = 0
    while count < length_
        current = UInt8[]
        readbytes!(s, current, 1)
        count += 1
        current_byte = current[1]
        if current_byte >= 1 && current_byte <= 75
            n = current_byte
            instruction = UInt8[]
            readbytes!(s, instruction, n)
            push!(instructions, instruction)
            count += n
        elseif current_byte == 76
            # op_pushdata1
            n = UInt8[]
            readbytes!(s, n, 1)
            instruction = UInt8[]
            readbytes!(s, instruction, n[1])
            push!(instructions, instruction)
            count += n[1] + 1
        elseif current_byte == 77
            # op_pushdata2
            n = UInt8[]
            readbytes!(s, n, 2)
            n = reinterpret(Int16, n)[1]
            instruction = UInt8[]
            readbytes!(s, instruction, n)
            push!(instructions, instruction)
            count += n + 2
        else
            # op_code
            push!(instructions, current_byte)
        end
    end
    if count != length_
        error("Error: parsing Script failed")
    end
    return Script(instructions)
end

function rawserialize(s::Script)
    result = UInt8[]
    for instruction in s.instructions
        if typeof(instruction) == UInt8
            append!(result, instruction)
        else
            length_ = length(instruction)
            if length_ < 0x4b
                append!(result, UInt8(length_))
            elseif length_ > 0x4b && length_ < 0x100
                append!(result, 0x4c)
                append!(result, UInt8(length_))
            elseif length_ >= 0x100 && length_ <= 0x208
                append!(result, 0x4d)
                result += int2bytes(length_, 2)
            else
                error("too long an instruction")
            end
            append!(result, instruction)
        end
    end
    return result
end

function scriptserialize(s::Script)
    result = rawserialize(s)
    total = length(result)
    prepend!(result, encode_varint(total))
    return result
end

"""
    scriptevaluate(s::Script, z::Integer) -> Bool

Evaluate if Script is valid given the transaction signature hash
"""
function scriptevaluate(s::Script, z::Integer)
    instructions = copy(s.instructions)
    stack = Array{UInt8,1}[]
    altstack = Array{UInt8,1}[]
    while length(instructions) > 0
        instruction = popfirst!(instructions)
        if typeof(instruction) <: Integer
            operation = OP_CODE_FUNCTIONS[instruction]
            function badop(instruction::Integer)
                println("bad op: ", OP_CODE_NAMES[instruction])
            end
            if instruction in (99, 100)
                # op_if/op_notif require the  array
                if !operation(stack, instructions)
                    badop(instruction)
                    return false
                end
            elseif instruction in (107, 108)
                # op_toaltstack/op_fromaltstack require the altstack
                if !operation(stack, altstack)
                    badop(instruction)
                    return false
                end
            elseif instruction in (172, 173, 174, 175)
                if !operation(stack, z)
                    badop(instruction)
                    return false
                end
            elseif !operation(stack)
                badop(instruction)
                return false
            end
        else
            push!(stack, instruction)
        end
    end
    if length(stack) == 0
        return false
    end
    if pop!(stack) == Array{UInt8,1}[]
        return false
    end
    return true
end

"""
Takes a hash160 and returns the p2pkh scriptPubKey
"""
function p2pkh_script(h160::Array{UInt8,1})
    script = Union{UInt8, Array{UInt8,1}}[]
    pushfirst!(script, 0x76, 0xa9)
    push!(script, h160, 0x88, 0xac)
    return Script(script)
end

"""
Takes a hash160 and returns the p2sh scriptPubKey
"""
function p2sh_script(h160::Array{UInt8,1})
    script = Union{UInt8, Array{UInt8,1}}[]
    pushfirst!(script, 0xa9)
    push!(script, h160, 0x87)
    return Script(h160)
end


text:
"""
     add_singleton(x::AbstractArray, ::Val{dim}) where {dim}

Add an additional dimension `dim` of size 1 to array `x`.
"""
function add_singleton(x::AbstractArray, ::Val{dim}) where {dim}
    shape = ntuple(max(ndims(x) + 1, dim)) do i
        return i < dim ? size(x, i) : (i > dim ? size(x, i - 1) : 1)
    end
    return reshape(x, shape)
end

"""
    dot_matwise(x::AbstractArray, y::AbstractArray)

Compute the inner product of all matrices in `x` and `y`.

At least one of `x` and `y` has to be a matrix.
"""
dot_matwise(x::AbstractMatrix, y::AbstractMatrix) = dot(x, y)
function dot_matwise(x::AbstractArray, y::AbstractMatrix)
    xmat = reshape(x, size(x, 1) * size(x, 2), :)
    return reshape(reshape(y, 1, :) * xmat, size(x)[3:end])
end
dot_matwise(x::AbstractMatrix, y::AbstractArray) = dot_matwise(y, x)

"""
    checksize(μ::AbstractVecOrMat, ν::AbstractVecOrMat, C::AbstractMatrix)

Check that dimensions of source and target marginals `μ` and `ν` are consistent with cost
matrix `C`.
"""
function checksize(μ::AbstractVecOrMat, ν::AbstractVecOrMat, C::AbstractMatrix)
    size(C) == (size(μ, 1), size(ν, 1)) || throw(
        DimensionMismatch("cost matrix `C` must be of size `(size(μ, 1), size(ν, 1))`")
    )
    return nothing
end

"""
    checksize2(x::AbstractVecOrMat, y::AbstractVecOrMat)

Check if arrays `x` and `y` are compatible, then return a tuple of its broadcasted second
dimension.
"""
checksize2(::AbstractVector, ::AbstractVector) = ()
function checksize2(μ::AbstractVecOrMat, ν::AbstractVecOrMat)
    size_μ_2 = size(μ, 2)
    size_ν_2 = size(ν, 2)
    if size_μ_2 > 1 && size_ν_2 > 1 && size_μ_2 != size_ν_2
        throw(DimensionMismatch("size of source and target marginals is not compatible"))
    end
    return (max(size_μ_2, size_ν_2),)
end

"""
     isallapprox(x::AbstractVector)

Check that all entries of `x` are approximately equal
"""
function isallapprox(x::AbstractVecOrMat)
    return all(y -> isapprox(y, x[1]), x[2:end])
end

"""
     checkbalanced(μ::AbstractVecOrMat, ν::AbstractVecOrMat)

Check that source and target marginals `μ` and `ν` are balanced.
"""
function checkbalanced(μ::AbstractVector, ν::AbstractVector)
    sum(μ) ≈ sum(ν) || throw(ArgumentError("source and target marginals are not balanced"))
    return nothing
end
function checkbalanced(x::AbstractVecOrMat, y::AbstractVecOrMat)
    all(isapprox.(sum(x; dims=1), sum(y; dims=1))) ||
        throw(ArgumentError("source and target marginals are not balanced"))
    return nothing
end
function checkbalanced(x::AbstractMatrix)
    isallapprox(sum(x; dims=1)) ||
        throw(ArgumentError("source and target marginals are not balanced"))
    return nothing
end

"""
    A_batched_mul_B!(c::AbstractVector, A::AbstractMatrix, b::AbstractVector)

Compute the matrix-vector product `Ab` and write the result to `c`.
"""
A_batched_mul_B!(c::AbstractVector, A::AbstractMatrix, b::AbstractVector) = mul!(c, A, b)

"""
    A_batched_mul_B!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)

Compute the matrix-matrix product `AB` and write the result to `C`.
"""
A_batched_mul_B!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix) = mul!(C, A, B)

"""
    A_batched_mul_B!(C::AbstractMatrix, A::AbstractArray{<:Any,3}, B::AbstractMatrix)

Compute the matrix-vector products `A[:, :, i] B[:, i]` and write them to `C[:, i]`.
"""
function A_batched_mul_B!(C::AbstractMatrix, A::AbstractArray{<:Any,3}, B::AbstractMatrix)
    return NNlib.batched_mul!(add_singleton(C, Val(2)), A, add_singleton(B, Val(2)))
end

"""
    At_batched_mul_B!(c::AbstractVector, A::AbstractMatrix, b::AbstractVector)

Compute the matrix-vector product `transpose(A)b` and write the result to `c`.
"""
function At_batched_mul_B!(c::AbstractVector, A::AbstractMatrix, b::AbstractVector)
    return mul!(c, transpose(A), b)
end

"""
    At_batched_mul_B!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)

Compute the matrix-matrix product `transpose(A)B` and write the result to `C`.
"""
function At_batched_mul_B!(C::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix)
    return mul!(C, transpose(A), B)
end

"""
    At_batched_mul_B!(C::AbstractMatrix, A::AbstractArray{<:Any,3}, B::AbstractMatrix)

Compute the matrix-vector products `transpose(A[:, :, i]) B[:, i]` and write them to
`C[:, i]`.
"""
function At_batched_mul_B!(C::AbstractMatrix, A::AbstractArray{<:Any,3}, B::AbstractMatrix)
    return NNlib.batched_mul!(add_singleton(C, Val(1)), add_singleton(B, Val(1)), A)
end

struct FiniteDiscreteMeasure{X<:AbstractVector,P<:AbstractVector}
    support::X
    p::P

    function FiniteDiscreteMeasure{X,P}(support::X, p::P) where {X,P}
        length(support) == length(p) || error("length of `support` and `p` must be equal")
        isprobvec(p) || error("`p` must be a probability vector")
        return new{X,P}(support, p)
    end
end

"""
    discretemeasure(
        support::AbstractVector,
        probs::AbstractVector{<:Real}=fill(inv(length(support)), length(support))
    )

Construct a finite discrete probability measure with `support` and corresponding
`probabilities`. If the probability vector argument is not passed, then
equal probability is assigned to each entry in the support.

# Examples
```julia
using KernelFunctions
# rows correspond to samples
μ = discretemeasure(RowVecs(rand(7,3)), normalize!(rand(10),1))

# columns correspond to samples, each with equal probability
ν = discretemeasure(ColVecs(rand(3,12)))
```

!!! note
    If `support` is a 1D vector, the constructed measure will be sorted,
    e.g. for `mu = discretemeasure([3, 1, 2],[0.5, 0.2, 0.3])`, then
    `mu.support` will be `[1, 2, 3]` and `mu.p` will be `[0.2, 0.3, 0.5]`.
"""
function discretemeasure(
    support::AbstractVector{<:Real},
    probs::AbstractVector{<:Real}=fill(inv(length(support)), length(support)),
)
    return DiscreteNonParametric(support, probs)
end
function discretemeasure(
    support::AbstractVector,
    probs::AbstractVector{<:Real}=fill(inv(length(support)), length(support)),
)
    return FiniteDiscreteMeasure{typeof(support),typeof(probs)}(support, probs)
end

Distributions.support(d::FiniteDiscreteMeasure) = d.support
Distributions.probs(d::FiniteDiscreteMeasure) = d.p


text:
# C:\apps\Julia_1.5.0\bin\julia.exe

include("Ecl.jl")
#
#julia> @time a = EclFile("C:\\Test\\NORNE_ATW2013.SMSPEC");
#  0.037358 seconds (197.04 k allocations: 12.778 MiB)
#
#julia> @time a = EclFile("C:\\Test\\NORNE_ATW2013.UNSMRY");
#  1.472356 seconds (7.19 M allocations: 326.285 MiB, 25.26% gc time)
#

using Dates

function timedelta(delta, origin)
    # delta in days
    Dates.DateTime(origin) + Dates.Second(delta * 24 * 60 * 60)
end

function to_datetime(arr, origin)
    orig = Dates.DateTime(origin)
    #tmp = Array{Dates.DateTime}
    tmp = Dates.DateTime[]
    for delta in arr
        push!(tmp, orig + Dates.Second(delta * 24 * 60 * 60))
    end
    return tmp
end

#times = df["TIME::+:+:+:+"]
#
#[timedelta(i, "1997-11-06") for i in times]

insertcols!(df, 1, :DATETIME => to_datetime(df["TIME::+:+:+:+"], "1997-11-06") )
#to_datetime(df["TIME::+:+:+:+"], "1997-11-06")


text:
using NLsolve
using Distributions
using LinearAlgebra
using Printf
using BenchmarkTools

# Need to figure out how to a type system for solvers.
# E.g. Potential solvers inlcude forward solvers, splitting methods and ode solvers
abstract type Solver end
abstract type ForwardSolver <: Solver end
abstract type OperatorSplitting <: Solver end

# Forward Iteration
struct forward_iteration <: ForwardSolver
    tol
    maxIter
    verbose
    forward_iteration(;tol=1E-4, maxIter=1000, verbose=false) = new(tol, maxIter, verbose)
end
function (solver::forward_iteration)(f, z0)
    zp, z = z0, f(z0)
    for iter in range(1, stop=solver.maxIter)
        zp, z = z, f(z)
        err = norm(z - zp)
        if err < solver.tol
            if solver.verbose
                @printf("forward solve: \tIters = %d\terror=%1.2E\n", iter, err)
            end
            return z
        end
    end

    @printf("Solver did not converge\tError = %1.2E\n", norm(z - zp))
    return z
end

# Newton method
struct newton <: ForwardSolver
    tol
    maxIter
    verbose
    autodiff
    newton(;tol=1E-4, maxIter=1000, verbose=false, autodiff=:forward) = new(tol, maxIter, verbose, autodiff)
end
function (solver::newton)(f::Function, z0)
    function f!(F, z)
        F .= f(z) - z
    end
    res = nlsolve(f!, z0; ftol=solver.tol, xtol=solver.tol, method=:newton, autodiff=solver.autodiff)
    if solver.verbose
        println(res)
    end
    return res.zero
end

# Solver using Anderson Acceleration
struct anderson <: ForwardSolver
    tol
    maxIter
    verbose
    autodiff
    anderson(;tol=1E-4, maxIter=1000, verbose=false, autodiff=:forward) = new(tol, maxIter, verbose, autodiff)
end
function (solver::anderson)(f::Function, z0)
    function f!(F, z)
        F .= f(z) - z
    end
    res = nlsolve(f!, z0; ftol=solver.tol, xtol=solver.tol, method=:anderson, autodiff=solver.autodiff)
    if solver.verbose
        println(res)
    end
    if !res.x_converged && !res.f_converged
        println("Solver may not have converged.")
    end
    return res.zero
end

# Trust Region Solver
struct trust_region <: ForwardSolver
    tol
    maxIter
    verbose
    autodiff
    trust_region(;tol=1E-4, maxIter=1000, verbose=false, autodiff=:forward) = new(tol, maxIter, verbose, autodiff)
end
function (solver::trust_region)(f::Function, z0)
    function f!(F, z)
        F .= f(z) - z
    end
    res = nlsolve(f!, z0; ftol=solver.tol, xtol=solver.tol, method=:trust_region, autodiff=solver.autodiff)
    if solver.verbose
        println(res)
    end
    return res.zero
end

## Splitting methods
mutable struct PeacemanRachford <: OperatorSplitting
    tol
    maxIter
    α
    verbose
    cg
    PeacemanRachford(;tol=1E-4, maxIter=1000, α=1, verbose=false, cg=true) = new(tol, maxIter, α, verbose, cg)
end
function (solver::PeacemanRachford)(RA, RB, z0)
    # Takes in the resolvent operators of the operators A and B
    error = 1E10
    uk, zk, u_12, z_12, zn, un = zero(z0), zero(z0), zero(z0), zero(z0), zero(z0), zero(z0)

    error_log = []
    for k in (1:solver.maxIter)
        u_12 = 2zk - uk
        z_12 = RB(u_12)
        un = 2z_12 - u_12
        zn = RA(un)
        error = norm(zn - zk)
        zk = zn
        uk = un
        # println(norm(error))

        if isnan(error)
            break
        end
        append!(error_log, error)
        if error < solver.tol
            if solver.verbose
                println(k, ",\t", error)
            end
            return zk
        end
    end
    println("PeacemanRachford did not converge. Error: ", error)
    return zk
end

mutable struct DouglasRachford <: OperatorSplitting
    tol
    maxIter
    α
    verbose
    DouglasRachford(;tol=1E-4, maxIter=1000, α=1, verbose=false) = new(tol, maxIter, α, verbose)
end
function (solver::DouglasRachford)(RA, RB, z0)
    # Takes in the resolvent operators of the operators A and B
    error = 1E10
    uk, zk = z0, z0
    for k in (1:solver.maxIter)
        x_12 =  RB(zk)
        println("x_12 = ", norm(x_12))
        z_12 = 2x_12 - zk
        xn = RA(z_12)
        println("x_n = ", norm(x_n))
        zn = zk + xn - x_12
        error = norm(zn - zk) / (norm(zk) + 1E-7)
        if error < solver.tol
            println("\tError ", error, "\tIter ", k)
            return zk
        end
        zk, xk = zn, xn
    end
    println("DouglasRachford did not converge. Error: ", error)
    return zk
end

function test_solvers()

    f(x) = sin.(x)
    solver = anderson()

    z0 = [0.1]
    res = solver(f, z0)

    z = res.zero
end

max_stars_repo_name: imanchester/REN

max_issues_repo_name: imanchester/REN

max_forks_repo_name: imanchester/REN


text:
fp = dirname(@__FILE__)
@testset "Creating lagged data" begin
    # lag
    matdata = load(joinpath(fp, "../../reference/test_lag.jld2"))
    nonzero = findfirst(vec(sum(matdata["lag_y"], dims = 2)) .!== 0.)
    lags = Int(matdata["lags"])
    @test DSGE.lag(matdata["y"], lags) ≈ matdata["lag_y"][nonzero:end, :]
    @test Matrix(DSGE.lag(Matrix(matdata["y"]'), lags;
              T_by_n = false)') ≈ matdata["lag_y"][nonzero:end, :]
    @test DSGE.lag(matdata["y"], lags;
              drop_obs = lags ÷ 2) ≈
                  matdata["lag_y"][nonzero + (lags ÷ 2):end, :]
    @test Matrix(DSGE.lag(Matrix(matdata["y"]'), lags;
                     drop_obs = lags ÷ 2,
                     T_by_n = false)') ≈ matdata["lag_y"][nonzero + (lags ÷ 2):end, :]

    # lag with padding
    padded = DSGE.lag(matdata["y"], lags, pad = true)
    padded_dropobs1 = DSGE.lag(matdata["y"], lags, pad = true,
                         drop_obs = lags ÷ 2)
    padded_dropobs2 = DSGE.lag(matdata["y"], lags, pad = true,
                          drop_obs = lags + 1)
    padded_dropobs3 = Matrix(DSGE.lag(Matrix(matdata["y"]'), lags, pad = true,
                         drop_obs = lags ÷ 2, T_by_n = false)')
    padded_dropobs4 = Matrix(DSGE.lag(Matrix(matdata["y"]'), lags, pad = true,
                          drop_obs = lags + 1, T_by_n = false)')
    padded[isnan.(padded)] .= 0.
    padded_dropobs1[isnan.(padded_dropobs1)] .= 0.
    padded_dropobs3[isnan.(padded_dropobs3)] .= 0.
    @test padded ≈ matdata["lag_y"]
    @test padded_dropobs1 ≈ matdata["lag_y"][(lags ÷ 2) + 1:end, :]
    @test padded_dropobs2 ≈ matdata["lag_y"][(lags + 2):end, :]
    @test padded_dropobs3 ≈ matdata["lag_y"][(lags ÷ 2) + 1:end, :]
    @test padded_dropobs4 ≈ matdata["lag_y"][(lags + 2):end, :]

    # lag_data
    @test DSGE.lag_data(Matrix(matdata["y"]'), lags) ≈ matdata["x"][lags + 1:end, :]
    @test DSGE.lag_data(Matrix(matdata["y"]'), lags,
                   pad = true)[lags + 1:end, :] ≈ matdata["x"][lags + 1:end, :]
    tmp = DSGE.lag_data(Matrix(matdata["y"]'), lags, pad = true)
    tmp[isnan.(tmp)] .= 0.
    @test tmp ≈ matdata["x"]
    @test DSGE.lag_data(Matrix(matdata["y"]'), lags,
                   use_intercept = false) ≈ matdata["x"][lags + 1:end, 2:end]
    @test DSGE.lag_data(Matrix(matdata["y"]'), lags, pad = true,
                   padding = zeros(lags, size(matdata["x"], 2)))[lags + 1:end, :] ≈
                       matdata["x"][lags + 1:end, :]
end

@testset "Check calculation of population moments" begin
    matdata = load(joinpath(fp, "../../reference/test_var_population_moments.jld2"))
    lags = Int(matdata["lags"])
    @test DSGE.lag_data(Matrix(matdata["YY"]'), Int(matdata["lags"])) ≈
        matdata["XX"][lags + 1:end, :]
    YYYY, XXYY, XXXX = DSGE.compute_var_population_moments(Matrix(matdata["YY"]'), lags)
    @test YYYY ≈ matdata["YY"][lags + 1:end, :]' * matdata["YY"][lags + 1:end, :]
    @test XXYY ≈ matdata["XX"][lags + 1:end, 2:end]' * matdata["YY"][lags + 1:end, :]
    @test XXXX ≈ matdata["XX"][lags + 1:end, 2:end]' * matdata["XX"][lags + 1:end, 2:end]

    # Use intercept now
    YYYY, XXYY, XXXX = DSGE.compute_var_population_moments(Matrix(matdata["YY"]'), lags;
                                                    use_intercept = true)
    @test YYYY ≈ matdata["YY"][lags + 1:end, :]' * matdata["YY"][lags + 1:end, :]
    @test XXYY ≈ matdata["XX"][lags + 1:end, 1:end]' * matdata["YY"][lags + 1:end, :]
    @test XXXX ≈ matdata["XX"][lags + 1:end, 1:end]' * matdata["XX"][lags + 1:end, 1:end]
end

fp = dirname(@__FILE__)
@testset "Draw stationary VAR coefficients and variance-covariance matrix" begin
    matdata = load(joinpath(fp, "../../reference/test_dsgevar_stationary_draws.jld2"))
    β_draw, Σ_draw = DSGE.draw_stationary_VAR(matdata["yyyyd"], matdata["xxyyd"],
                                         matdata["xxxxd"], Int(matdata["nobsc"]),
                                         Int(matdata["nvar"]), Int(matdata["nlags"]);
                                         testing = true, test_Σ_draw_shock = matdata["draw"],
                                         test_β_draw_shock = vec(matdata["beta_draw"]))
    @test matdata["cct_sim"] ≈ convert(Matrix{Float64}, β_draw')
    @test matdata["sig_sim"] ≈ Σ_draw
end

@testset "Draw VECM coefficients and variance-covariance matrix" begin
    matdata = load(joinpath(fp, "../../reference/draw_vecm.jld2"))
    β_draw, Σ_draw = DSGE.draw_VECM(matdata["yyyyc"], matdata["xxyyc"],
                                    matdata["xxxxc"], Int(matdata["Tbar"]),
                                    Int(matdata["n_obs"]), Int(matdata["nlags"]),
                                    Int(matdata["coint"]);
                                    testing = true, test_Σ_draw_shock = matdata["sigma_draw"],
                                    test_β_draw_shock = vec(matdata["beta_draw"]))
    @test maximum(abs.(matdata["beta_ans"] - β_draw)) < 1e-2
    @test @test_matrix_approx_eq matdata["sigma_ans"] Σ_draw
end


text:
using AnalogDiscovery

function analogOutGatedSineSettings(hdwf::Hdwf,ch::Int32)
  analogOutNodeEnableSet(hdwf,ch,ANALOGOUTNODECARRIER,true)
  analogOutNodeFunctionSet(hdwf,ch,ANALOGOUTNODECARRIER,FUNCSINE)
  analogOutNodeFrequencySet(hdwf,ch,ANALOGOUTNODECARRIER,10.0)
  analogOutNodeAmplitudeSet(hdwf,ch,ANALOGOUTNODECARRIER,0.05)
  analogOutNodeOffsetSet(hdwf,ch,ANALOGOUTNODECARRIER,0.0)

  analogOutNodeEnableSet(hdwf,ch,ANALOGOUTNODEAM,true)
  analogOutNodeFunctionSet(hdwf,ch,ANALOGOUTNODEAM,FUNCSQUARE)
  analogOutNodeFrequencySet(hdwf,ch,ANALOGOUTNODEAM,5.0)
  analogOutNodeAmplitudeSet(hdwf,ch,ANALOGOUTNODEAM,200.0)
  analogOutNodeOffsetSet(hdwf,ch,ANALOGOUTNODEAM,100.0)

  analogOutRunSet(hdwf,ch,0.2)
  analogOutRepeatSet(hdwf,ch,Int32(3))
end

function go()
  cDev = enumDevices()
  @show ("Available",cDev)
  myDev = Int32(0)
  if cDev > 0
    if enumDeviceIsOpened(myDev)
      error("Error opening device, already in use.")
    end
    hdwf = deviceOpen(myDev)
  else
    error("No devices available.")
  end

  ch = Int32(0)
  analogOutGatedSineSettings(hdwf,ch)

  for in in 1:10
    @time begin analogOutNodeFrequencySet(hdwf,ch,ANALOGOUTNODECARRIER,4000.0)
    analogOutConfigure(hdwf,ch,true) end
    sleep(3.0)
    @time begin analogOutNodeFrequencySet(hdwf,ch,ANALOGOUTNODECARRIER,7000.0)
    analogOutConfigure(hdwf,ch,true) end
    sleep(3.0)
  end

  deviceClose(hdwf)
end

go()


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module SwingFactorial

using Nemo

export ModuleSwingFactorial
export Sfactorial

"""

Basic implementation of the swing algorithm using no primes. Claims to be the most efficient simple algorithm to compute the factorial. An advanced version based on prime-factorization is available as the prime-swing factorial factorialPS.

* Sfactorial
"""
const ModuleSwingFactorial = ""

"""

Return the factorial of ``n``. Basic implementation of the swing algorithm using no primes. An advanced version based on prime-factorization is available as the prime-swing factorial factorialPS.
"""
function Sfactorial(n::Int)::fmpz

    smallOddFactorial =           fmpz[0x0000000000000000000000000000001,
    0x0000000000000000000000000000001, 0x0000000000000000000000000000001,
    0x0000000000000000000000000000003, 0x0000000000000000000000000000003,
    0x000000000000000000000000000000f, 0x000000000000000000000000000002d,
    0x000000000000000000000000000013b, 0x000000000000000000000000000013b,
    0x0000000000000000000000000000b13, 0x000000000000000000000000000375f,
    0x0000000000000000000000000026115, 0x000000000000000000000000007233f,
    0x00000000000000000000000005cca33, 0x0000000000000000000000002898765,
    0x00000000000000000000000260eeeeb, 0x00000000000000000000000260eeeeb,
    0x0000000000000000000000286fddd9b, 0x00000000000000000000016beecca73,
    0x000000000000000000001b02b930689, 0x00000000000000000000870d9df20ad,
    0x0000000000000000000b141df4dae31, 0x00000000000000000079dd498567c1b,
    0x00000000000000000af2e19afc5266d, 0x000000000000000020d8a4d0f4f7347,
    0x000000000000000335281867ec241ef, 0x0000000000000029b3093d46fdd5923,
    0x0000000000000465e1f9767cc5866b1, 0x0000000000001ec92dd23d6966aced7,
    0x0000000000037cca30d0f4f0a196e5b, 0x0000000000344fd8dc3e5a1977d7755,
    0x000000000655ab42ab8ce915831734b, 0x000000000655ab42ab8ce915831734b,
    0x00000000d10b13981d2a0bc5e5fdcab, 0x0000000de1bc4d19efcac82445da75b,
    0x000001e5dcbe8a8bc8b95cf58cde171, 0x00001114c2b2deea0e8444a1f3cecf9,
    0x0002780023da37d4191deb683ce3ffd, 0x002ee802a93224bddd3878bc84ebfc7,
    0x07255867c6a398ecb39a64b83ff3751, 0x23baba06e131fc9f8203f7993fc1495]

    function oddProduct(m::Int, len::Int)
        if len < 24
            p = fmpz(m)
            for k ∈ 2:2:2(len - 1)
                p *= (m - k)
            end
            return p
        end
        hlen = len >> 1
        oddProduct(m - 2 * hlen, len - hlen) * oddProduct(m, hlen)
    end

    function oddFactorial(n)
        if n < 41
            oddFact = smallOddFactorial[1+n]
            sqrOddFact = smallOddFactorial[1+div(n, 2)]
        else
            sqrOddFact, oldOddFact = oddFactorial(div(n, 2))
            len = div(n - 1, 4)
            (n % 4) != 2 && (len += 1)
            high = n - ((n + 1) & 1)
            oddSwing = div(oddProduct(high, len), oldOddFact)
            oddFact = sqrOddFact^2 * oddSwing
        end
        (oddFact, sqrOddFact)
    end

    n < 0 && ArgumentError("n must be ≥ 0")
    if n == 0
        return fmpz(1)
    end
    sh = n - count_ones(n)
    oddFactorial(n)[1] << sh
end

#START-TEST-########################################################

using Test

function test()
    @testset "SwingFactorial" begin
        for n ∈ 0:999
            S = Sfactorial(n)
            B = Base.factorial(BigInt(n))
            @test S == B
        end
    end
end

function demo() end

function perf()
    # n = 10^6
    # 0.372059 seconds     (   911 allocations:  48.628 MiB,  1.12% gc time)
    # 1.033867 seconds     (2.57 M allocations: 153.290 MiB, 13.31% gc time)

    # n = 7250000
    # 3.759344 seconds (8.86 k allocations: 483.946 MiB, 7.25% gc time)
    # 9.922223 seconds (22.30 M allocations: 1.316 GiB, 10.44% gc time)

    # n = 10^7
    #  6.256117 seconds   (11.34 k allocations: 682.953 MiB,  3.32% gc time)
    # 15.678733 seconds   (25.56 M allocations:   2.081 GiB, 10.72% gc time)

    # n = 10^8
    #  96.296497 seconds (381.72 k allocations:   9.686 GiB, 1.76% gc time)
    # 237.375884 seconds (258.40 M allocations:  28.133 GiB, 6.73% gc time)

    GC.gc()
    n = 1000000
    @time Base.factorial(BigInt(n))
    @time Sfactorial(n)
    println("n = $n")
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
abstract type BoundaryCondition end

"""
    FluidBC(momentum    = Impenetrable(NoSlip())
            temperature = Insulating())

The standard boundary condition for CNSEModel. The default options imply a "no flux" boundary condition.
"""
Base.@kwdef struct FluidBC{M, T} <: BoundaryCondition
    momentum::M = Impenetrable(NoSlip())
    temperature::T = Insulating()
end

abstract type StateBC end
abstract type MomentumBC <: StateBC end
abstract type MomentumDragBC <: StateBC end
abstract type TemperatureBC <: StateBC end

(bc::StateBC)(state, aux, t) = bc.flux(bc.params, state, aux, t)

"""
    Impenetrable(drag::MomentumDragBC) :: MomentumBC

Defines an impenetrable wall model for momentum. This implies:
  - no flow in the direction normal to the boundary, and
  - flow parallel to the boundary is subject to the `drag` condition.
"""
struct Impenetrable{D <: MomentumDragBC} <: MomentumBC
    drag::D
end

"""
    Penetrable(drag::MomentumDragBC) :: MomentumBC

Defines an penetrable wall model for momentum. This implies:
  - no constraint on flow in the direction normal to the boundary, and
  - flow parallel to the boundary is subject to the `drag` condition.
"""
struct Penetrable{D <: MomentumDragBC} <: MomentumBC
    drag::D
end

"""
    NoSlip() :: MomentumDragBC

Zero momentum at the boundary.
"""
struct NoSlip <: MomentumDragBC end

"""
    FreeSlip() :: MomentumDragBC

No surface drag on momentum parallel to the boundary.
"""
struct FreeSlip <: MomentumDragBC end

"""
    MomentumFlux(stress) :: MomentumDragBC

Applies the specified kinematic stress on momentum normal to the boundary.
Prescribe the net inward kinematic stress across the boundary by `stress`,
a function with signature `stress(problem, state, aux, t)`, returning the flux (in m²/s²).
"""
Base.@kwdef struct MomentumFlux{𝒯, 𝒫} <: MomentumDragBC
    flux::𝒯 = nothing
    params::𝒫 = nothing
end

"""
    Insulating() :: TemperatureBC

No temperature flux across the boundary
"""
struct Insulating <: TemperatureBC end

"""
    TemperatureFlux(flux) :: TemperatureBC

Prescribe the net inward temperature flux across the boundary by `flux`,
a function with signature `flux(problem, state, aux, t)`, returning the flux (in m⋅K/s).
"""
Base.@kwdef struct TemperatureFlux{𝒯, 𝒫} <: TemperatureBC
    flux::𝒯 = nothing
    params::𝒫 = nothing
end

function check_bc(bcs, label)
    bctype = FluidBC

    bc_ρu = check_bc(bcs, Val(:ρu), label)
    bc_ρθ = check_bc(bcs, Val(:ρθ), label)

    return bctype(bc_ρu, bc_ρθ)
end

function check_bc(bcs, ::Val{:ρθ}, label)
    if haskey(bcs, :ρθ)
        if haskey(bcs[:ρθ], label)
            return bcs[:ρθ][label]
        end
    end

    return Insulating()
end

function check_bc(bcs, ::Val{:ρu}, label)
    if haskey(bcs, :ρu)
        if haskey(bcs[:ρu], label)
            return bcs[:ρu][label]
        end
    end

    return Impenetrable(FreeSlip())
end

# these functions just trim off the extra arguments
function _cnse_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    bc,
    model,
    state⁺,
    aux⁺,
    n,
    state⁻,
    aux⁻,
    t,
    _...,
)
    return cnse_boundary_state!(nf, bc, model, state⁺, aux⁺, n, state⁻, aux⁻, t)
end

function _cnse_boundary_state!(
    nf::NumericalFluxSecondOrder,
    bc,
    model,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    return cnse_boundary_state!(
        nf,
        bc,
        model,
        state⁺,
        gradflux⁺,
        aux⁺,
        n,
        state⁻,
        gradflux⁻,
        aux⁻,
        t,
    )
end


text:
################################################################################
#
#  Integer functions
#
################################################################################

function rem(a::fmpz, b::UInt)
  return ccall((:fmpz_fdiv_ui, :libflint), UInt, (Ptr{fmpz}, UInt), &a, b)
end

function isless(a::BigFloat, b::Nemo.fmpz)
  if _fmpz_is_small(b)
    c = ccall((:mpfr_cmp_si, :libmpfr), Int32, (Ptr{BigFloat}, Int), &a, b.d)
  else
    c = ccall((:mpfr_cmp_z, :libmpfr), Int32, (Ptr{BigFloat}, UInt), &a, unsigned(b.d) << 2)
  end
  return c < 0
end

function mulmod(a::UInt, b::UInt, n::UInt, ni::UInt)
  ccall((:n_mulmod2_preinv, :libflint), UInt, (UInt, UInt, UInt, UInt), a, b, n, ni)
end

# TODO (CF):
# should be Bernstein'ed: this is slow for large valuations
# returns the maximal v s.th. z mod p^v == 0 and z div p^v
#   also useful if p is not prime....
#
# TODO: what happens to z = 0???

function remove(z::T, p::T) where T <: Integer
  z == 0 && return (0, z)
  v = 0
  while mod(z, p) == 0
    z = div(z, p)
    v += 1
  end
  return (v, z)
end

function remove(z::Rational{T}, p::T) where T <: Integer
  z == 0 && return(0, z)
  v, d = remove(denominator(z), p)
  w, n = remove(numerator(z), p)
  return w-v, n//d
end

function valuation(z::T, p::T) where T <: Integer
  z == 0 && return 0
  const v = 0
  while mod(z, p) == 0
    z = div(z, p)
    v += 1
  end
  return v
end

function valuation(z::Rational{T}, p::T) where T <: Integer
  z == 0 && error("Not yet implemented")
  v = valuation(denominator(z), p)
  w = valuation(numerator(z), p)
  return w-v
end

function remove!(a::fmpz, b::fmpz)
  v = ccall((:fmpz_remove, :libflint), Clong, (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), a, a, b)
  return v, a
end

function remove!(a::fmpq, b::fmpz)
  nr = ccall((:fmpq_numerator_ptr, :libflint), Ptr{fmpz}, (Ref{fmpq}, ), a)
  vn = ccall((:fmpz_remove, :libflint), Clong, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), nr, nr, b)
  #fmpq's are simplified: either num OR den will be non-trivial
  if vn != 0
    return vn, a
  end
  nr = ccall((:fmpq_denominator_ptr, :libflint), Ptr{fmpz}, (Ref{fmpq}, ), a)
  vn = ccall((:fmpz_remove, :libflint), Clong, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), nr, nr, b)
  return -vn, a
end

function valuation!(a::fmpq, b::fmpz)
  nr = ccall((:fmpq_numerator_ptr, :libflint), Ptr{fmpz}, (Ref{fmpq}, ), a)
  vn = ccall((:fmpz_remove, :libflint), Clong, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), nr, nr, b)
  #fmpq's are simplified: either num OR den will be non-trivial
  if vn != 0
    return vn
  end
  nr = ccall((:fmpq_denominator_ptr, :libflint), Ptr{fmpz}, (Ref{fmpq}, ), a)
  vn = ccall((:fmpz_remove, :libflint), Clong, (Ptr{fmpz}, Ptr{fmpz}, Ref{fmpz}), nr, nr, b)
  return -vn
end

function *(a::fmpz, b::BigFloat)
  return BigInt(a)*b
end

function BigFloat(a::fmpq)
  r = BigFloat(0)
  ccall((:fmpq_get_mpfr, :libflint), Void, (Ref{BigFloat}, Ref{fmpq}, Int32), r, a, __get_rounding_mode())
  return r
end

function isless(a::Float64, b::fmpq) return a<b*1.0; end
function isless(a::fmpq, b::Float64) return a*1.0<b; end

function isless(a::Float64, b::fmpz) return a<b*1.0; end
function isless(a::fmpz, b::Float64) return a*1.0<b; end

function (a::FlintIntegerRing)(b::fmpq)
  !isone(denominator(b)) && error("Denominator not 1")
  return deepcopy(numerator(b))
end

function //(a::fmpq, b::fmpz)
  return a//fmpq(b)
end

function //(a::fmpz, b::fmpq)
  return fmpq(a)//b
end

function *(a::fmpz, b::Float64)
  return BigInt(a)*b
end

function *(b::Float64, a::fmpz)
  return BigInt(a)*b
end

function *(a::fmpq, b::Float64)
  return Rational(a)*b
end

function *(b::Float64, a::fmpq)
  return Rational(a)*b
end

function Float64(a::fmpq)
  b = a*fmpz(2)^53
  Float64(div(numerator(b), denominator(b)))/(Float64(2)^53) #CF 2^53 is bad in 32bit
end

function convert(R::Type{Rational{Base.GMP.BigInt}}, a::Nemo.fmpz)
  return R(BigInt(a))
end

log(a::fmpz) = log(BigInt(a))
log(a::fmpq) = log(numerator(a)) - log(denominator(a))

function round(::Type{fmpz}, a::fmpq)
  s = sign(numerator(a))
  n = abs(numerator(a))
  d = denominator(a)
  if isone(d)
    return numerator(a)
  end
  q = div(n, d)
  r = mod(n, d)
  if r >= div(d, 2)
    return s*(q+1)
  else
    return s*q
  end
end

function round(a::fmpq)
  return round(fmpz, a)
end

function zero(::Type{Nemo.fmpz})
  return fmpz(0)
end

function one(::Type{Nemo.fmpz})
  return fmpz(1)
end

@doc """
  modord(a::fmpz, m::fmpz) -> Int
  modord(a::Integer, m::Integer)

  The multiplicative order of a modulo m (not a good algorithm).
""" ->
function modord(a::fmpz, m::fmpz)
  gcd(a,m)!=1 && throw("1st agrument not a unit")
  i = 1
  b = a % m
  while b != 1
    i += 1
    b = b*a % m
  end
  return i
end

function modord(a::Integer, m::Integer)
  gcd(a,m)!=1 && throw("1st agrument not a unit")
  i = 1
  b = a % m
  while b != 1
    i += 1
    b = b*a % m
  end
  return i
end

function isodd(a::fmpz)
  ccall((:fmpz_is_odd, :libflint), Int, (Ref{fmpz},), a) == 1
end

function iseven(a::fmpz)
  ccall((:fmpz_is_even, :libflint), Int, (Ref{fmpz},), a) == 1
end

function neg!(a::fmpz)
  ccall((:fmpz_neg, :libflint), Void, (Ref{fmpz}, Ref{fmpz}), a, a)
  return a
end

##
## to support rand(fmpz:fmpz)....
##
# note, we cannot get a UnitRange as this is only legal for subtypes of Real

function colon(a::Int, b::fmpz)
  return fmpz(a):b
end

function one(::fmpz)
  return fmpz(1)
end

function zero(::fmpz)
  return fmpz(0)
end

#the build-in show fails due to Integer(a::fmpz) not defined
# I don't think it would efficient to provide this for printing
#display() seems to NOT call my show function, why, I don't know
function Integer(a::fmpz)
  return BigInt(a)
end

function show(io::IO, a::StepRange{fmpz, fmpz})
  println(io, "2-element ", typeof(a), ":\n ", a.start, ",", a.stop)
end

#TODO
# need to be mapped onto proper Flint primitives
# flints needs a proper interface to randomness - I think
# currently one simply cannot use it at all
#
# should be tied(?) to the Julia rng stuff?
# similarly, all the derived rand functions should probably also do this
#
# inspired by/ copied from the Base/random.jl
#
function rand(rng::AbstractRNG, a::StepRange{fmpz, fmpz})
  m = Base.last(a) - Base.first(a)
  m < 0 && throw("range empty")
  nd = ndigits(m, 2)
  nl, high = divrem(nd, 8*sizeof(Base.GMP.Limb))
  if high>0
    mask = m>>(nl*8*sizeof(Base.GMP.Limb))
  end
  s = fmpz(0)
  while true
    s = fmpz(0)
    for i=1:nl
      s = s << (8*sizeof(Base.GMP.Limb))
      s += rand(Base.GMP.Limb)
    end
    if high >0
      s = s << high
      s += rand(0:Base.GMP.Limb(mask))
    end
    if s <= m break; end
  end
  return s + first(a)
end

function length(a::StepRange{fmpz, fmpz})
  return a.stop - a.start +1
end

struct RangeGeneratorfmpz# <: Base.Random.RangeGenerator
  a::StepRange{fmpz, fmpz}
end

function Base.Random.RangeGenerator(r::StepRange{fmpz,fmpz})
    m = last(r) - first(r)
    m < 0 && throw(ArgumentError("range must be non-empty"))
    return RangeGeneratorfmpz(r)
end

function rand(rng::AbstractRNG, g::RangeGeneratorfmpz)
  return rand(rng, g.a)
end

function Base.getindex(a::StepRange{fmpz,fmpz}, i::fmpz)
  a.start+(i-1)*Base.step(a)
end

################################################################################
#
#  Should go to Nemo?
#
################################################################################

one(::Type{fmpq}) = fmpq(1)

############################################################
# more unsafe function that Bill does not want to have....
############################################################

function divexact!(z::fmpz, x::fmpz, y::fmpz)
    y == 0 && throw(DivideError())
    ccall((:fmpz_divexact, :libflint), Void,
          (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
    return z
end

function lcm!(z::fmpz, x::fmpz, y::fmpz)
   ccall((:fmpz_lcm, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
   return z
end

function gcd!(z::fmpz, x::fmpz, y::fmpz)
   ccall((:fmpz_gcd, :libflint), Void,
         (Ref{fmpz}, Ref{fmpz}, Ref{fmpz}), z, x, y)
   return z
end

function mul!(z::fmpz, x::fmpz, y::Int)
  ccall((:fmpz_mul_si, :libflint), Void, (Ref{fmpz}, Ref{fmpz}, Int), z, x, y)
  return z
end

function mul!(z::fmpz, x::fmpz, y::UInt)
  ccall((:fmpz_mul_ui, :libflint), Void, (Ref{fmpz}, Ref{fmpz}, UInt), z, x, y)
  return z
end

function mul!(z::fmpz, x::fmpz, y::Integer)
  mul!(z, x, fmpz(y))
  return z
end

function one!(a::fmpz)
  ccall((:fmpz_one, :libflint), Void, (Ref{fmpz}, ), a)
  return a
end

################################################################################
#
#  power detection
#
################################################################################

doc"""
    ispower(a::fmpz) -> Int, fmpz
    ispower(a::Integer) -> Int, Integer
> Writes $a = r^e$ with $e$ maximal. Note: $1 = 1^0$.
"""
function ispower(a::fmpz)
  if iszero(a)
    error("must not be zero")
  end
  if isone(a)
    return 0, a
  end
  if a < 0
    e, r = ispower(-a)
    if isone(e)
      return 1, a
    end
    v, s = remove(e, 2)
    return s, -r^(2^v)
  end
  rt = fmpz()
  e = 1
  while true
    ex = ccall((:fmpz_is_perfect_power, :libflint), Int, (Ref{fmpz}, Ref{fmpz}), rt, a)
    if ex == 1 || ex == 0
      return e, a
    end
    e *= ex
    a = rt
  end
end

function ispower(a::Integer)
  e,r = ispower(fmpz(a))
  return e, typeof(a)(r)
end

doc"""
    ispower(a::fmpq) -> Int, fmpq
    ispower(a::Rational) -> Int, Rational
> Writes $a = r^e$ with $e$ maximal. Note: $1 = 1^0$.
"""
function ispower(a::fmpq)
  e, r = ispower(numerator(a))
  if e==1
    return e, a
  end
  f, s = ispower(denominator(a))
  g = gcd(e, f)
  return g, r^div(e, g)//s^div(f, g)
end

function ispower(a::Rational)
  T = typeof(denominator(a))
  e, r = ispower(fmpq(a))
  return e, T(numerator(r))//T(denominator(r))
end

doc"""
    ispower(a::fmpz, n::Int) -> Bool, fmpz
    ispower(a::fmpq, n::Int) -> Bool, fmpq
    ispower(a::Integer, n::Int) -> Bool, Integer
> Tests if $a$ is an $n$-th power. Return {{{true}}} and the root if successful.
"""
function ispower(a::fmpz, n::Int)
  b = root(a, n)
  return b^n==a, b
end

function ispower(a::Integer, n::Int)
  return ispower(fmpz(a), n)
end

function ispower(a::fmpq, n::Int)
  fl, nu = ispower(numerator(a), n)
  if !fl
    return fl, a
  end
  fl, de = ispower(denominator(a), n)
  return fl, fmpq(nu, de)
end

################################################################################
#
#  Chinese remaindering modulo UInts to fmpz
#
################################################################################

mutable struct fmpz_comb
  primes::Ptr{UInt}
  num_primes::Int
  n::Int
  comb::Ptr{Ptr{fmpz}}
  res::Ptr{Ptr{fmpz}}
  mod_n::UInt
  mod_ninv::UInt
  mod_norm::UInt

  function fmpz_comb(primes::Array{UInt, 1})
    z = new()
    ccall((:fmpz_comb_init, :libflint), Void, (Ptr{fmpz_comb}, Ptr{UInt}, Int),
            &z, primes, length(primes))
    finalizer(z, _fmpz_comb_clear_fn)
    return z
  end
end

function _fmpz_comb_clear_fn(z::fmpz_comb)
  ccall((:fmpz_comb_clear, :libflint), Void, (Ptr{fmpz_comb}, ), &z)
end

mutable struct fmpz_comb_temp
  n::Int
  comb_temp::Ptr{Ptr{fmpz}}
  temp::Ptr{fmpz}
  temp2::Ptr{fmpz}

  function fmpz_comb_temp(comb::fmpz_comb)
    z = new()
    ccall((:fmpz_comb_temp_init, :libflint), Void,
            (Ptr{fmpz_comb_temp}, Ptr{fmpz_comb}), &z, &comb)
    finalizer(z, _fmpz_comb_temp_clear_fn)
    return z
  end
end

function _fmpz_comb_temp_clear_fn(z::fmpz_comb_temp)
  ccall((:fmpz_comb_temp_clear, :libflint), Void, (Ptr{fmpz_comb_temp}, ), &z)
end

function fmpz_multi_crt_ui!(z::fmpz, a::Array{UInt, 1}, b::fmpz_comb, c::fmpz_comb_temp)
  ccall((:fmpz_multi_CRT_ui, :libflint), Void,
          (Ptr{fmpz}, Ptr{UInt}, Ptr{fmpz_comb}, Ptr{fmpz_comb_temp}, Cint),
          &z, a, &b, &c, 1)
  return z
end

function _fmpz_preinvn_struct_clear_fn(z::fmpz_preinvn_struct)
  ccall((:fmpz_preinvn_clear, :libflint), Void, (Ptr{fmpz_preinvn_struct}, ), &z)
end

function fdiv_qr_with_preinvn!(q::fmpz, r::fmpz, g::fmpz, h::fmpz, hinv::fmpz_preinvn_struct)
  ccall((:fmpz_fdiv_qr_preinvn, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}, Ptr{fmpz}, Ptr{fmpz}, Ptr{fmpz_preinvn_struct}), &q, &r, &g, &h, &hinv)
end

function submul!(z::fmpz, x::fmpz, y::fmpz)
  ccall((:fmpz_submul, :libflint), Void, (Ptr{fmpz}, Ptr{fmpz}, Ptr{fmpz}), &z, &x, &y)
end

################################################################################
#
#  Number of bits
#
################################################################################

function nbits(a::BigInt)
  return ndigits(a, 2)
end

doc"""
  nbits(a::Int) -> Int
  nbits(a::UInt) -> Int
  nbits(a::BigInt) -> Int

  Returns the number of bits necessary to represent a
"""
function nbits(a::Int)
  a==0 && return 0
  return Int(ceil(log(abs(a))/log(2)))
end
function nbits(a::UInt)
  a==0 && return 0
  return Int(ceil(log(a)/log(2)))
end

################################################################################
#
#  Modular reduction with symmetric residue system
#
################################################################################

function mod_sym(a::fmpz, b::fmpz)
  c = mod(a,b)
  @assert c>=0
  if b > 0 && 2*c > b
    return c-b
  elseif b < 0 && 2*c > -b
    return c+b
  else
    return c
  end
end

doc"""
    isinteger(a::fmpq) -> Bool

> Returns true iff the denominator of $a$ is one.
"""
function isinteger(a::fmpq)
  return isone(denominator(a))
end

################################################################################
#
#  sunit group
#
################################################################################

mutable struct MapSUnitGrpZFacElem{T} <: Map{T, FacElemMon{FlintRationalField}, HeckeMap, MapSUnitGrpZFacElem}
  header::MapHeader
  idl::Array{fmpz, 1}

  function MapSUnitGrpZFacElem{T}() where {T}
    return new{T}()
  end
end

function show(io::IO, mC::MapSUnitGrpZFacElem)
  println(io, "SUnits (in factored form) map of $(codomain(mC)) for $(mC.idl)")
end

mutable struct MapSUnitGrpZ{T} <: Map{T, FlintRationalField, HeckeMap, MapSUnitGrpZ}
  header::MapHeader
  idl::Array{fmpz, 1}

  function MapSUnitGrpZ{T}() where {T}
    return new{T}()
  end
end

function show(io::IO, mC::MapSUnitGrpZ)
  println(io, "SUnits map of $(codomain(mC)) for $(mC.idl)")
end

doc"""
***
    sunit_group_fac_elem(S::Array{fmpz, 1}) -> GrpAbFinGen, Map
    sunit_group_fac_elem(S::Array{Integer, 1}) -> GrpAbFinGen, Map
> The $S$-unit group of $Z$ supported at $S$: the group of
> rational numbers divisible only by primes in $S$.
> The second return value is the map mapping group elements to rationals
> in factored form or rationals back to group elements.
"""
function sunit_group_fac_elem(S::Array{T, 1}) where T <: Integer
  return sunit_group_fac_elem(fmpz[x for x=S])
end

function sunit_group_fac_elem(S::Array{fmpz, 1})
  S = coprime_base(S)  #TODO: for S-units use factor???
  G = DiagonalGroup(vcat([fmpz(2)], fmpz[0 for i=S]))
  S = vcat(fmpz[-1], S)

  mp = MapSUnitGrpZFacElem{typeof(G)}()
  mp.idl = S

  Sq = fmpq[x for x=S]

  function dexp(a::GrpAbFinGenElem)
    return FacElem(Sq, [a.coeff[1,i] for i=1:length(S)])
  end

  function dlog(a::fmpz)
    g = [a>=0 ? 0 : 1]
    g = vcat(g, [valuation(a, x) for x=S[2:end]])
    return G(g)
  end

  function dlog(a::Integer)
    return dlog(fmpz(a))
  end

  function dlog(a::Rational)
    return dlog(fmpq(a))
  end

  function dlog(a::fmpq)
    return dlog(numerator(a)) - dlog(denominator(a))
  end

  function dlog(a::FacElem)
    return sum([e*dlog(k) for (k,e) = a.fac])
  end

  mp.header = Hecke.MapHeader(G, FacElemMon(FlintQQ), dexp, dlog)

  return G, mp
end

doc"""
***
    sunit_group(S::Array{fmpz, 1}) -> GrpAbFinGen, Map
    sunit_group(S::Array{Integer, 1}) -> GrpAbFinGen, Map
> The $S$-unit group of $Z$ supported at $S$: the group of
> rational numbers divisible only by primes in $S$.
> The second return value is the map mapping group elements to rationals
> or rationals back to group elements.
"""
function sunit_group(S::Array{T, 1}) where T <: Integer
  return sunit_group(fmpz[x for x=S])
end

function sunit_group(S::Array{fmpz, 1})
  u, mu = sunit_group_fac_elem(S)

  mp = MapSUnitGrpZ{typeof(u)}()
  mp.idl = S

  function dexp(a::GrpAbFinGenElem)
    return evaluate(image(mu, a))
  end

  mp.header = Hecke.MapHeader(u, FlintQQ, dexp, mu.header.preimage)

  return u, mp
end

Hecke.gcd(a::fmpz, b::Integer) = Hecke.gcd(a, fmpz(b))
Hecke.gcd(a::Integer, b::fmpz) = Hecke.gcd(fmpz(a), b)
Hecke.lcm(a::fmpz, b::Integer) = Hecke.lcm(a, fmpz(b))
Hecke.lcm(a::Integer, b::fmpz) = Hecke.lcm(fmpz(a), b)

doc"""
    isprime_power(n::fmpz) -> Bool
    isprime_power(n::Integer) -> Bool
> Tests is $n$ is the exact power of a prime number.
"""
function isprime_power(n::fmpz)
  e, p = ispower(n)
  return isprime(p)
end

function isprime_power(n::Integer)
  return isprime_power(fmpz(n))
end

################################################################################
# random and factor
################################################################################

factor(a::RingElem) = Nemo.factor(a)
factor(a::Integer) = factor(fmpz(a))

mutable struct flint_rand_ctx_t
  a::Ptr{Void}
  function flint_rand_ctx_t()
    return new()
  end
end

function show(io::IO, A::flint_rand_ctx_t)
  println(io, "Flint random state\n")
end

function flint_rand_state()
  A = flint_rand_ctx_t()
  A.a = ccall((:flint_rand_alloc, :libflint), Ptr{Void}, (Int, ), 1)
  ccall((:flint_randinit, :libflint), Void, (Ptr{Void}, ), A.a)

  function clean_rand_state(A::flint_rand_ctx_t)
    ccall((:flint_randclear, :libflint), Void, (Ptr{Void}, ), A.a)
    ccall((:flint_rand_free, :libflint), Void, (Ptr{Void}, ), A.a)
    nothing
  end
  finalizer(A, clean_rand_state)
  return A
end

global flint_rand_ctx

function ecm(a::fmpz, B1::UInt, B2::UInt, ncrv::UInt, rnd = flint_rand_ctx)
  f = fmpz()
  r = ccall((:fmpz_factor_ecm, :libflint), Int32, (Ptr{fmpz}, UInt, UInt, UInt, Ptr{Void}, Ptr{fmpz}), &f, ncrv, B1, B2, rnd.a, &a)
  return r, f
end

function ecm(a::fmpz, B1::Int, B2::Int, ncrv::Int, rnd = flint_rand_ctx)
  return ecm(a, UInt(B1), UInt(B2), UInt(ncrv), rnd)
end

#data from http://www.mersennewiki.org/index.php/Elliptic_Curve_Method
B1 = [2, 11, 50, 250, 1000, 3000, 11000, 43000, 110000, 260000, 850000, 2900000];
nC = [25, 90, 300, 700, 1800, 5100, 10600, 19300, 49000, 124000, 210000, 340000];

function ecm(a::fmpz, max_digits::Int = div(ndigits(a), 2)+1, rnd = flint_rand_ctx)
  n = ndigits(a, 10)
  B1s = 15

  i = 1
  s = div(max_digits-15, 5)+2
  s = max(i, s)
  while i <= s
    e,f = ecm(a, B1[i]*1000, B1[i]*1000*100, nC[i], rnd)
    if e != 0
      return (e,f)
    end
    i += 1
    if i > length(B1)
      return (e, f)
    end
  end
  return (Int32(0), a)
end

function factor_trial_range(N::fmpz, start::Int=0, np::Int=10^5)
   F = Nemo.fmpz_factor()
   ccall((:fmpz_factor_trial_range, :libflint), Void, (Ptr{Nemo.fmpz_factor}, Ptr{fmpz}, UInt, UInt), &F, &N, start, np)
   res = Dict{fmpz, Int}()
   for i in 1:F.num
     z = fmpz()
     ccall((:fmpz_factor_get_fmpz, :libflint), Void,
           (Ptr{fmpz}, Ptr{Nemo.fmpz_factor}, Int), &z, &F, i - 1)
     res[z] = unsafe_load(F.exp, i)
   end
   return res, canonical_unit(N)
end

big_primes = fmpz[]

function factor(N::fmpz)
  N_in = N
  global big_primes
  r, c = factor_trial_range(N)
  for (p, v) = r
    N = divexact(N, p^v)
  end
#  if length(r)==0
#    r[1] = 1
#  end
#  @assert prod(p^v for (p, v) = r)*N == N_in
  if isunit(N)
#    @assert all(isprime, keys(r))
    @assert N == c
    return Nemo.Fac(c, r)
  end
  N *= c
  @assert N > 0

  for p = big_primes
    v, N = remove(N, p)
    if v > 0
      @assert !haskey(r, p)
      r[p] = v
    end
  end
#  @assert prod(p^v for (p, v) = r)*N == N_in
  fac, N = ispower(N)

 #TODO: problem(s)
 # Nemo.factor = mpqs is hopeless if > n digits, but asymptotically and practically
 # faster than ecm.
 # ecm is much better if there are "small" factors.
 # p-1 and p+1 methods are missing
 # so probably
 # if n is small enough -> Nemo
 # if n is too large: ecm
 # otherwise
 #  need ecm to find small factors
 # then recurse...

  if ndigits(N) > 60  # TODO: needs revision!
    e, f = ecm(N, div(ndigits(N), 3))
    #TODO: use coprime basis to refine stuff...
    while e != 0
      ee, f = ispower(f)
      ee = valuation(N, f) #careful, f does not need to be prime, so N/f^ee is not coprime to f
      if isprime(f)
        if haskey(r, f)
          r[f] += fac*ee
        else
          r[f] = fac*ee
        end
      else
        s = factor(f)
        for (p, ex) = s.fac
          if haskey(r, p)
            r[p] += fac*ex*ee
          else
            r[p] = fac*ex*ee
          end
        end
      end
  #    @assert N % f^ee == 0
  #    @assert N % f^(ee+1) != 0
      N = divexact(N, f^ee)
  #    @assert prod(p^v for (p, v) = r)*N == N_in
      if isone(N)
        break
      end
      e, f = ecm(N, div(ndigits(N), 3))
    end
  end
  s = Nemo.factor(N)
  for (p, ex) = s.fac
    if haskey(r, p)
      r[p] += fac*ex
    else
      r[p] = fac*ex
    end
  end
  for p = keys(r)
    if nbits(p) > 60 && !(p in big_primes)
      push!(big_primes, p)
    end
  end
  @assert all(isprime, keys(r))
  @assert prod(a^b for (a,b) = r) * c == N_in
  return Nemo.Fac(c, r)
end

function ceil(::Type{fmpz}, a::BigFloat)
  return fmpz(ceil(BigInt, a))
end

function floor(::Type{fmpz}, a::BigFloat)
  return fmpz(floor(BigInt, a))
end

function round(::Type{fmpz}, a::BigFloat)
  return fmpz(round(BigInt, a))
end

/(a::BigFloat, b::fmpz) = a/BigInt(b)


text:
export interpolate,solution

"""
    interpolate(mesh,vector;n=2^3)

Interpolate a vector `vector` defined on a uniform collocation grid defined by `mesh`.

Returns `(new_mesh,new_vector)` a new uniform mesh with `n` times as many values, and the vector of values on collocation points.

"""
function interpolate(mesh::Mesh,vector;n=2^3::Int)

    fourier=fft(vector)
    m=Int(mesh.N/2)
    new_fourier = [fourier[1:m] ;zeros((n-1)*2*m) ;fourier[m+1:end]]

    new_mesh = Mesh(mesh.xmin,mesh.xmax,n*mesh.N)
    new_vector=ifft(new_fourier)*n
    if all((vector[i] isa Real) for i in length(vector))
        new_vector=real.(new_vector)
    end
    return new_mesh,new_vector
end

"""
    interpolate(mesh,vector,x;fast)

Interpolate a vector `vector` defined on a uniform collocation grid defined by `mesh`, on collocation points given by `x`.

If the collocation points `x` are regularly spaced and the optional keyword argument `fast` is set to `true` (default is `false`),
then the algorithm is faster and uses less allocations, but is less precise.

Returns the vector of values on collocation points.

"""
function interpolate(mesh::Mesh,vector,x;fast=false)

    fourier=fft(vector)
    k = mesh.k
	x₀ = mesh.xmin
	if fast == false
		#same as new_vector=exp.(1im*(x.-x₀)*k')*fourier/length(k)
		new_vector = similar(x,Complex); z = complex.(zero(k));
		for i in 1:length(x)
			z .= exp.(-1im*(x[i].-x₀)*k)
			new_vector[i] = real(dot(z,fourier))
		end
		new_vector ./= length(k)
	else
		y=x[2:end]-x[1:end-1];y.-=y[1];
        if maximum(abs.(y))>8*eps(maximum(x))
            @error("Collocation points must be equally spaced.")
        end
		new_vector = similar(x,Complex); z = exp.(-1im*(x[1].-x₀)*k);eidk=exp.(-1im*(x[2]-x[1])*k);
		for i in 1:length(x)
			new_vector[i] = real(dot(z,fourier))
			z .*= eidk
		end
		new_vector ./= length(k)
	end

    if all((vector[i] isa Real) for i in length(vector))
        new_vector=real.(new_vector)
    end
    return new_vector
end

"""
    solution(p::Problem;t,x,interpolation)

Gives the solution of a solved initial-value at a given time `t`.

# Arguments
- Argument `p` is of type `Problem`.
- Keyword argument `t` is optional, the last computed time is returned by default.
- Keyword argument `x` is optional, if provided the solution is interpolated to the collocation vector `x`.
- Keyword argument `interpolation` is optional, if an integer is provided the solution is interpolated on as many collocation points (if `true`, then the default value `2^3` is chosen).

# Return values
Provides `(η,v,x,t)` where
- `η` is the surface deformation at collocation points;
- `v` is the tangential velocity (derivative of the trace of the velocity potential) at collocation points;
- `x` is the vector of collocation points;
- `t` the time (first computed time greater or equal to provided `t`).

"""
function solution(p::Problem; t=nothing, x=nothing, interpolation = false)
	if t == nothing t = p.times.tfin end
	t=min(max(t,0),p.times.tfin)
	index = indexin(false,p.times.ts.<t)[1]
	t=p.times.ts[index]
	if Symbol(typeof(p.model)) == :WaterWaves
		if x != nothing || interpolation != false
			@warn "cannot interpolate non-regularly spaced mesh."
		end
		(x,η,v) = (p.model.mapfro)(p.data.U[index])

	else
    	(η,v) = (p.model.mapfro)(p.data.U[index])
		mesh=p.mesh
		if x != nothing
			η = interpolate(mesh,η,x)
			v = interpolate(mesh,v,x)
		else
			x=mesh.x
		end
		if interpolation == true
			new_mesh,η = interpolate(mesh,η)
			new_mesh,v = interpolate(mesh,v)
			x = new_mesh.x
		elseif isa(interpolation,Int)
			new_mesh,η = interpolate(mesh,η;n=interpolation)
			new_mesh,v = interpolate(mesh,v;n=interpolation)
			x = new_mesh.x
		end
	end
	return η,v,x,t
end


text:
with(Paper("a4", "landscape")) do
    move(0.5, 0.5)
    draw(0.25, 0.25)
end

with(Paper("letter"), Pen("red", 0.05), Scale(100, 100)) do
    forwards(100)
    right(45)
    up()
    forwards(10)
    down()
end

with(File("drawing.ps"), Paper("a4", "landscape")) do
    move(0.5, 0.5)
    with(Pen("red")) do
        draw(0.1, 0.1)
    end
end

with(Paper("a4", orientation="landscape")) do
    move(0.5, 0.5)
    with(Pen("red")) do
        draw(0.1, 0.1)
    end
    write_to("drawing.ps")
end


text:
# # Chapter 5

WORD_SIZE

bitstring(3)

bitstring(-3)

isbitstype(Int64)

isbitstype(String)

myadd(x, y) = x + y
@code_native myadd(1,2)

typemax(Int64)

bitstring(typemax(Int32))

typemin(Int64)

bitstring(typemin(Int32))

9223372036854775806 + 1

9223372036854775806 + 1 + 1

2^64

2^65

# ##BigInt

big(9223372036854775806) + 1 + 1

big(2)^64

# ## Floating Point

bitstring(2.5)
bitstring(-2.5)

function floatbits(x::Float64)
   b = bitstring(x)
   b[1:1]*"|"*b[2:12]*"|"*b[13:end]
end

floatbits(2.5)

floatbits(-2.5)

0.1 > 1//10

Rational(0.1)

float(big(Rational(0.1)))

bitstring(0.10000000000000001) == bitstring(0.1)

eps(0.1)

nextfloat(0.1)

floatbits(0.1)

floatbits(nextfloat(0.1))

# ##Unchecked conversions

UInt64(1)
UInt32(4294967297)
4294967297 % UInt32
@btime UInt32(1)
@btime 1 % UInt32

# ## FastMath

function sum_diff(x)
    n = length(x); d = 1/(n-1)
    s = zero(eltype(x))
    s = s +  (x[2] - x[1]) / d
    for i = 2:length(x)-1
        s =  s + (x[i+1] - x[i+1]) / (2*d)
    end
    s = s + (x[n] - x[n-1])/d
end

function sum_diff_fast(x)
    n=length(x); d = 1/(n-1)
    s = zero(eltype(x))
    @fastmath s = s +  (x[2] - x[1]) / d
    @fastmath for i = 2:n-1
        s =  s + (x[i+1] - x[i+1]) / (2*d)
    end
    @fastmath s = s + (x[n] - x[n-1])/d
end

t=rand(2000);
@btime sum_diff($t)
@btime sum_diff_fast($t)

macroexpand(Main, :(@fastmath a + b / c))

# ## KBN Summation

sum([1 1e-100 -1])

using KahanSummation

sum_kbn([1 1e-100 -1])

t=[1, -1, 1e-100];

@btime sum($t)

@btime sum_kbn($t)

# ## Subnormal numbers

issubnormal(1.0)

issubnormal(1.0e-308)

3e-308 - 3.001e-308

issubnormal(3e-308 - 3.001e-308)

set_zero_subnormals(true)

3e-308 - 3.001e-308

3e-308 == 3.001e-308

get_zero_subnormals()

function timestep( b, a, dt )
    n = length(b)
    b[1] = 1
    two = eltype(b)(2)
    for i=2:n-1
        b[i] = a[i] + (a[i-1] - two*a[i] + a[i+1]) * dt
    end
    b[n] = 0
end

function heatflow( a, nstep )
    b = similar(a)
    o = eltype(a)(0.1)
    for t=1:div(nstep,2)
        timestep(b,a,o)
        timestep(a,b,o)
    end
end

set_zero_subnormals(false)

t=rand(1000);

@btime heatflow($t, 1000)

set_zero_subnormals(true)

@btime heatflow($t, 1000)


text:
begin

    N = 100

    A = rand(N)

    #stack up some functions
    @et res_et = exp(cos(sin(A)))
    res_native = exp(cos(sin(A)))

    @test res_et == res_native
end


text:
module DebugTest

using Test
using Parametron

import Parametron: mock_model, setdirty!, findallocs

@testset "findallocs" begin
    model = mock_model()
    x = Variable(model)
    p = Parameter{Int}(() -> 3, model)
    expr = @expression p * 4 + x
    setdirty!(model); findallocs(devnull, expr)
end

end


text:
using SparseArrays, DelimitedFiles, Random
using Test, LinearAlgebra
using IterativeSolvers, AlgebraicMultigrid
import AlgebraicMultigrid: Pinv, Classical
using JLD2
using FileIO

using Random: seed!

include("sa_tests.jl")

@testset "AlgebraicMultigrid Tests" begin

graph = include("test.jl")
ref_S = include("ref_S_test.jl")
ref_split = readdlm("ref_split_test.txt")

@testset "Strength of connection" begin

# classical strength of connection
A = poisson(5)
A = float.(A)
S, T = Classical(0.2)(A)
@test Matrix(S) == [ 1.0  0.5  0.0  0.0  0.0
                   0.5  1.0  0.5  0.0  0.0
                   0.0  0.5  1.0  0.5  0.0
                   0.0  0.0  0.5  1.0  0.5
                   0.0  0.0  0.0  0.5  1.0 ]
S, T = Classical(0.25)(graph)
diff = S - ref_S
@test maximum(diff) < 1e-10

end

@testset "Splitting" begin

# Ruge-Stuben splitting
S = poisson(7)
@test RS()(S) == [0, 1, 0, 1, 0, 1, 0]

a = include("thing.jl")
S, T = Classical(0.25)(a)
@test RS()(S) == [0, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0,
0, 0, 1, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0,
1, 0]

@test RS()(ref_S) == Int.(vec(ref_split))

end

@testset "Interpolation" begin

# Direct Interpolation
using AlgebraicMultigrid
A = poisson(5)
A = Float64.(A)
splitting = [1,0,1,0,1]
P, R = AlgebraicMultigrid.direct_interpolation(A, copy(A), splitting)
@test P ==  [ 1.0  0.0  0.0
              0.5  0.5  0.0
              0.0  1.0  0.0
              0.0  0.5  0.5
              0.0  0.0  1.0 ]
A = include("thing.jl")
ml = ruge_stuben(A)
@test size(ml.levels[2].A, 1) == 19
end

@testset "Coarse Solver" begin
A = float.(poisson(10))
b = A * ones(10)
@test sum(abs2, Pinv(A)(similar(b), b) - ones(10)) < 1e-6
end

@testset "Multilevel" begin
A = poisson(1000)
A = float.(A) #FIXME
ml = AlgebraicMultigrid.ruge_stuben(A)
@test length(ml) == 8
s = [1000, 500, 250, 125, 62, 31, 15]
n = [2998, 1498, 748, 373, 184, 91, 43]
for i = 1:7
    @test size(ml.levels[i].A, 1) == s[i]
    @test nnz(ml.levels[i].A) == n[i]
end
@test size(ml.final_A, 1) == 7
@test nnz(ml.final_A) == 19

A = include("randlap.jl")
ml = ruge_stuben(A)
@test length(ml) == 3
s = [100, 17]
n = [2066, 289]
for i = 1:2
    @test size(ml.levels[i].A, 1) == s[i]
    @test nnz(ml.levels[i].A) == n[i]
end
@test size(ml.final_A, 1) == 2
@test nnz(ml.final_A) == 4
@static if VERSION < v"0.7-"
    @test round(AlgebraicMultigrid.operator_complexity(ml), 3) ≈ 1.142
    @test round(AlgebraicMultigrid.grid_complexity(ml), 3) ≈ 1.190
else
    @test round(AlgebraicMultigrid.operator_complexity(ml), digits=3) ≈ 1.142
    @test round(AlgebraicMultigrid.grid_complexity(ml), digits=3) ≈ 1.190
end

include("gmg.jl")

A = float.(poisson(10^6))
ml = multigrid(A)
@test length(ml) == 10

end

@testset "Solver" begin
fsmoother = GaussSeidel(ForwardSweep())

A = poisson(1000)
A = float.(A)
ml = ruge_stuben(A)
x = solve(ml, A * ones(1000))
@test sum(abs2, x - ones(1000)) < 1e-8

ml = ruge_stuben(A, presmoother = fsmoother,
                     postsmoother = fsmoother)
x = solve(ml, A * ones(1000))
@test sum(abs2, x - ones(1000)) < 1e-8

A = include("randlap.jl")

ml = ruge_stuben(A, presmoother = fsmoother,
                    postsmoother = fsmoother)
x = solve(ml, A * ones(100))
@test sum(abs2, x - zeros(100)) < 1e-8

ml = ruge_stuben(A)
x = solve(ml, A * ones(100))
@test sum(abs2, x - zeros(100)) < 1e-6

end

@testset "Preconditioning" begin
A = include("thing.jl")
n = size(A, 1)
smoother = GaussSeidel(ForwardSweep())
ml = ruge_stuben(A, presmoother = smoother,
                    postsmoother = smoother)
p = aspreconditioner(ml)
b = zeros(n)
b[1] = 1
b[2] = -1
x = solve(p.ml, A * ones(n), maxiter = 1, tol = 1e-12)
diff = x - [  1.88664780e-16,   2.34982727e-16,   2.33917697e-16,
         8.77869044e-17,   7.16783490e-17,   1.43415460e-16,
         3.69199021e-17,   9.70950385e-17,   4.77034895e-17,
         3.77491328e-17,   5.07592420e-18,   6.32131628e-19,
        -1.60361276e-18,  -1.36749626e-16,  -4.16651794e-17,
        -3.48207590e-17,  -4.19334783e-17,  -4.60500098e-17,
        -6.91113945e-17,  -1.35997904e-16,  -9.85940056e-17,
         8.99433377e-17,   1.63924842e-16,  -2.43048120e-16,
        -3.21888624e-16,  -1.56389252e-16,  -3.66495834e-17,
         1.70606350e-16,   1.66788345e-16,  -3.26736922e-16,
        -3.39591125e-16,  -3.67075849e-16,  -3.89891523e-16,
        -4.34537627e-16,  -4.42863579e-16,  -6.62433333e-16,
        -5.56056397e-16,  -5.70242981e-16,  -6.48075679e-16,
        -6.58200572e-16,  -7.24886474e-16,  -7.55973538e-16,
        -6.76965535e-16,  -7.00643227e-16,  -6.23581397e-16,
        -7.03016682e-16]
@test sum(abs2, diff) < 1e-8
x = solve(p.ml, b, maxiter = 1, tol = 1e-12)
diff = x - [ 0.76347046, -0.5498286 , -0.2705487 , -0.15047352, -0.10248021,
        0.60292674, -0.11497073, -0.08460548, -0.06931461,  0.38230708,
       -0.055664  , -0.04854558, -0.04577031,  0.09964325,  0.01825624,
       -0.01990265, -0.02866185, -0.03049521,  0.03310897, -0.01709034,
       -0.02038031, -0.01325201, -0.01051535,  0.02992818,  0.01493605,
       -0.00633922, -0.01285614, -0.01155069, -0.01095907,  0.04415807,
        0.02213755,  0.018686  ,  0.02625713,  0.02007781,  0.01898018,
        0.02107552,  0.01909623,  0.01874986,  0.01852736,  0.01844719,
        0.01841821,  0.01841695,  0.01953195,  0.01885713,  0.01864432,
        0.0185079 ]
@test sum(abs2, diff) < 1e-8
x = cg(A, b, Pl = p)
diff = x - [ 0.82365077, -0.537589  , -0.30632349, -0.19370186, -0.14773294,
        0.68489145, -0.15550115, -0.1278148 , -0.11197922,  0.45362483,
       -0.08577219, -0.08598307, -0.08477946,  0.12985118,  0.02805496,
       -0.03907565, -0.05950957, -0.06544269,  0.05446686, -0.047537  ,
       -0.05203899, -0.04685981, -0.04491762,  0.05639249,  0.02792704,
       -0.02282528, -0.04062864, -0.04321821, -0.0441893 ,  0.07593055,
        0.05212038,  0.04464215,  0.05835841,  0.05079815,  0.04830733,
        0.05272397,  0.05028666,  0.0494817 ,  0.04960952,  0.0496615 ,
        0.04968258,  0.04968737,  0.05105749,  0.05009268,  0.04972329,
        0.04970173]
@test sum(abs2, diff) < 1e-8

# Symmetric GaussSeidel Smoothing

ml = ruge_stuben(A)
p = aspreconditioner(ml)

x = cg(A, b, Pl = p, maxiter = 100_000, tol = 1e-6)
diff = x - [0.823762, -0.537478, -0.306212, -0.19359, -0.147621, 0.685002,
            -0.155389, -0.127703, -0.111867, 0.453735, -0.0856607, -0.0858715,
            -0.0846678, 0.129962, 0.0281662, -0.0389642, -0.0593981, -0.0653311,
            0.0545782, -0.0474255, -0.0519275, -0.0467483, -0.0448061, 0.056504,
            0.0280386, -0.0227138, -0.0405172, -0.0431067, -0.0440778, 0.076042,
            0.052232, 0.0447537, 0.05847, 0.0509098, 0.0484189, 0.0528356,
            0.0503983, 0.0495933, 0.0497211, 0.0497731, 0.0497942, 0.049799,
            0.0511691, 0.0502043, 0.0498349, 0.0498134]
@test sum(abs2, diff) < 1e-8

x = solve(ml, b, maxiter = 1, tol = 1e-12)
diff = x - [0.775725, -0.571202, -0.290989, -0.157001, -0.106981, 0.622652,
            -0.122318, -0.0891874, -0.0709834, 0.392621, -0.055544, -0.0507485,
            -0.0466376, 0.107175, 0.0267468, -0.0200843, -0.0282827, -0.0299929,
            0.0420468, -0.0175585, -0.0181318, -0.0121591, -0.00902523, 0.0394795,
            0.019981, -0.00270916, -0.0106855, -0.0093661, -0.00837619, 0.052532,
            0.0301423, 0.0248904, 0.0333098, 0.0262179, 0.0246211, 0.026778,
            0.0245746, 0.0238448, 0.0233892, 0.0231593, 0.0230526, 0.0229771,
            0.0247913, 0.0238555, 0.0233681, 0.023096]
@test sum(abs2, diff) < 1e-8

end

@testset "Precision" begin

a = poisson(100)
b = rand(size(a,1))

# Iterate through all types
for (T,V) in ((Float64, Float64), (Float32,Float32),
        (Float64,Float32), (Float32,Float64))
        a = T.(a)
        ml = smoothed_aggregation(a)
        b = V.(b)
        c = cg(a, b, maxiter = 10)
        @test eltype(solve(ml, b)) == eltype(c)
end

end

# Smoothed Aggregation
@testset "Smoothed Aggregation" begin

@testset "Symmetric Strength of Connection" begin

    test_symmetric_soc()
end

@testset "Standard Aggregation" begin

    test_standard_aggregation()
end

@testset "Fit Candidates" begin
    test_fit_candidates()
end

@testset "Approximate Spectral Radius" begin
    test_approximate_spectral_radius()
end
end

@testset "Gauss Seidel" begin
    test_gauss_seidel()
end

@testset "Jacobi Prolongation" begin
    test_jacobi_prolongator()
end

@testset "Int32 support" begin
    a = sparse(Int32.(1:10), Int32.(1:10), rand(10))
    @inferred smoothed_aggregation(a)
end

# Issue #24
nodes_not_agg()

# Issue #26
test_symmetric_sweep()

# Issue #31
for sz in [10, 5, 2]
    a = poisson(sz)
    ml = ruge_stuben(a)
    @test isempty(ml.levels)
    @test size(ml.final_A) == (sz,sz)
    @test AlgebraicMultigrid.operator_complexity(ml) == 1
    @test AlgebraicMultigrid.grid_complexity(ml) == 1

    a = poisson(sz)
    ml = smoothed_aggregation(a)
    @test isempty(ml.levels)
    @test size(ml.final_A) == (sz,sz)
    @test AlgebraicMultigrid.operator_complexity(ml) == 1
    @test AlgebraicMultigrid.grid_complexity(ml) == 1
end

# Issue #46
for f in (smoothed_aggregation, ruge_stuben)

    a = load("bug.jld2")["G"]
    ml = f(a)
    p = aspreconditioner(ml)
    b = zeros(size(a,1))
    b[1] = 1
    b[2] = -1
    @test sum(abs2, a * solve(ml, b) - b) < 1e-10
    @test sum(abs2, a * cg(a, b, Pl = p, maxiter = 1000) - b) < 1e-10

end

end


text:
export CuArray, CuVector, CuMatrix, CuVecOrMat
export cupointer, curand, curandn, cuzeros, cuones

mutable struct CuArray{T,N} <: AbstractArray{T,N}
    ptr::Ptr{T}
    dims::Dims{N}
    dev::Int
end

const CuVector{T} = CuArray{T,1}
const CuMatrix{T} = CuArray{T,2}
const CuVecOrMat{T} = Union{CuVector{T},CuMatrix{T}}

function CuArray{T}(dims::Dims{N}) where {T,N}
    CUDAMalloc()(T, dims)
end
CuArray{T}(dims::Int...) where T = CuArray{T}(dims)
CuArray(x::Array{T,N}) where {T,N} = copyto!(CuArray{T}(size(x)), x)

### AbstractArray interface
getdevice(x::CuArray) = x.dev
Base.size(x::CuArray) = x.dims
Base.size(x::CuArray, i::Int) = i <= ndims(x) ? x.dims[i] : 1
Base.getindex(x::CuArray, i::Int) = throw("getindex $i")
Base.getindex(x::CuArray, I...) = CuArray(view(x,I...))
Base.setindex!(y::CuArray, v::Number, i::Int) = throw("setindex!")
Base.setindex!(y::CuArray{T}, x::CuArray{T}, I...) where T = copyto!(view(y,I), x)
Base.IndexStyle(::Type{<:CuArray}) = IndexLinear()
Base.length(x::CuArray) = prod(x.dims)
Base.similar(x::CuArray, ::Type{T}, dims::Dims{N}) where {T,N} = CuArray{T}(dims)
Base.similar(::Type{<:CuArray}, ::Type{T}, dims::Dims) where T = CuArray{T}(dims)
Base.ndims(x::CuArray{T,N}) where {T,N} = N
Base.eltype(x::CuArray{T,N}) where {T,N} = T
Base.strides(x::CuArray{T,1}) where T = (1,)
Base.strides(x::CuArray{T,2}) where T = (1,size(x,1))
Base.strides(x::CuArray{T,3}) where T = (1,size(x,1),size(x,1)*size(x,2))
Base.strides(x::CuArray{T,4}) where T = (1,size(x,1),size(x,1)*size(x,2),size(x,1)*size(x,2)*size(x,3))
function Base.stride(x::CuArray, i::Int)
    d = 1
    for i = 1:i-1
        d *= size(x, i)
    end
    d
end

Base.cconvert(::Type{Ptr{T}}, x::CuArray) where T = Base.cconvert(Ptr{T}, pointer(x))
function Base.pointer(x::CuArray{T}, index::Int=1) where T
    @assert index > 0
    x.ptr + sizeof(T)*(index-1)
end
cupointer(x::CuArray, index=1) = CuPtr(pointer(x,index))
Base.Array(x::CuArray{T}) where T = copyto!(Array{T}(undef,size(x)), x)
Base.unsafe_wrap(::Type{A}, ptr::Ptr{T}, dims) where {A<:CuArray,T} = CuArray(ptr, dims)

##### reshape #####
Base.reshape(x::CuArray{T}, dims::Dims{N}) where {T,N} = CuArray(x.ptr, dims, x.dev)
Base.reshape(x::CuArray{T}, dims::Int...) where T = reshape(x, dims)
Base.vec(x::CuArray{T}) where T = ndims(x) == 1 ? x : reshape(x,length(x))
function Base.dropdims(x::CuArray; dims)
    isa(dims,Int) && (dims = (dims,))
    for i in dims
        size(x,i) == 1 || throw(ArgumentError("dropdims must all be size 1"))
    end
    d = ()
    for i = 1:ndims(x)
        if !in(i, dims)
            d = tuple(d..., size(x,i))
        end
    end
    reshape(x, d)
end

##### copy #####
Base.copy(src::CuArray) = copyto!(similar(src), src)
function Base.copyto!(dest::Array{T}, src::CuArray{T}; stream=C_NULL) where T
    @assert length(dest) == length(src)
    copyto!(dest, 1, src, 1, length(dest), stream=stream)
end
function Base.copyto!(dest::CuArray{T}, src::Array{T}; stream=C_NULL) where T
    @assert length(dest) == length(src)
    copyto!(dest, 1, src, 1, length(dest), stream=stream)
end
function Base.copyto!(dest::CuArray{T}, src::CuArray{T}; stream=C_NULL) where T
    @assert length(dest) == length(src)
    copyto!(dest, 1, src, 1, length(dest), stream=stream)
end
function Base.copyto!(dest::Array{T}, doffs::Int, src::CuArray{T}, soffs::Int, n::Int; stream=C_NULL) where T
    p_dest = pointer(dest, doffs)
    p_src = pointer(src, soffs)
    @apicall :cuMemcpyDtoH (Ptr{Cvoid},Ptr{Cvoid},Csize_t) p_dest p_src n*sizeof(T)
    dest
end
function Base.copyto!(dest::CuArray{T}, doffs::Int, src::Array{T}, soffs::Int, n::Int; stream=C_NULL) where T
    p_dest = pointer(dest, doffs)
    p_src = pointer(src, soffs)
    @apicall :cuMemcpyHtoDAsync (Ptr{Cvoid},Ptr{Cvoid},Csize_t,Ptr{Cvoid}) p_dest p_src n*sizeof(T) stream
    dest
end
function Base.copyto!(dest::CuArray{T}, doffs::Int, src::CuArray{T}, soffs::Int, n::Int; stream=C_NULL) where T
    p_dest = pointer(dest, doffs)
    p_src = pointer(src, soffs)
    @apicall :cuMemcpyDtoDAsync (Ptr{Cvoid},Ptr{Cvoid},Csize_t,Ptr{Cvoid}) p_dest p_src n*sizeof(T) stream
    dest
end

function Base.fill!(x::CuArray{T}, value; stream=C_NULL) where T
    if sizeof(T) == 4
        ui = reinterpret(Cuint, T(value))
        @apicall :cuMemsetD32Async (Ptr{Cvoid},Cuint,Csize_t,Ptr{Cvoid}) x ui length(x) stream
    elseif sizeof(T) == 2
        us = reinterpret(Cushort, T(value))
        @apicall :cuMemsetD16Async (Ptr{Cvoid},Cushort,Csize_t,Ptr{Cvoid}) x us length(x) stream
    elseif sizeof(T) == 1
        uc = reinterpret(Cuchar, T(value))
        @apicall :cuMemsetD8Async (Ptr{Cvoid},Cuchar,Csize_t,Ptr{Cvoid}) x us length(x) stream
    else
        throw("Not supported.")
    end
    x
end

@generated function fill2!(x::CuArray{T}, value) where T
    Ct = cstring(T)
    k = Kernel("""
    __global__ void fill($Ct *x, int n, $Ct value) {
        int idx = blockIdx.x * blockDim.x + threadIdx.x;
        if (idx >= n) return;
        x[idx] = value;
    }""")
    quote
        gdims, bdims = cudims(length(x))
        $k(gdims, bdims, pointer(x), length(x), T(value))
        x
    end
end

##### IO #####
Base.show(io::IO, x::CuArray) = print(io, Array(x))
Base.display(x::CuArray) = display(Array(x))

#####
function curand(::Type{T}, dims::Dims{N}) where {T,N}
    # TODO: use curand library
    CuArray(rand(T,dims))
end
curand(::Type{T}, dims::Int...) where T = curand(T, dims)
curand(dims::Int...) = curand(Float64, dims)
function curandn(::Type{T}, dims::NTuple{N,Int}) where {T,N}
    # TODO: use curand library
    CuArray(randn(T,dims))
end
curandn(::Type{T}, dims::Int...) where T = curandn(T, dims)
curandn(dims::Int...) = curandn(Float64, dims)


text:
#
# Below:
#
# x = trg is the target variable,
# y = src is the source variable,
# E is the energy, and
# derivs allows for calculation of partial derivatives as well.
#

function lhelmfs(trg::Union{Float64,Complex{Float64}},energies::Float64;derivs::Bool=false)
    stdquad = 400
    h = 0.25
    meth = 1
    x1,x2 = reim(trg)
    lhfs(x1,x2,energies,derivs,stdquad,h,meth)
end

function lhelmfs(trg::Union{Vector{Float64},Vector{Complex{Float64}}},energies::Vector{Float64};derivs::Bool=false)
    n = length(trg)
    @assert n == length(energies)
    stdquad = 400
    h = 0.25
    meth = 1
    x1,x2 = reim(trg)
    u = Vector{Complex{Float64}}(undef,n)
    if derivs
        ux = Vector{Complex{Float64}}(undef,n)
        uy = Vector{Complex{Float64}}(undef,n)
        lhfs!(u,ux,uy,x1,x2,energies,derivs,stdquad,h,meth,n)
        return u,ux,uy
    else
        lhfs!(u,x1,x2,energies,derivs,stdquad,h,meth,n)
        return u
    end
end

function lhelmfs(trg::Union{Matrix{Float64},Matrix{Complex{Float64}}},E::Matrix{Float64};derivs::Bool=false)
    sizetrg,sizeE = size(trg),size(E)
    @assert sizetrg == sizeE

    if derivs
        u,ux,uy = lhelmfs(vec(trg),vec(E);derivs=derivs)
        return reshape(u,sizetrg),reshape(ux,sizetrg),reshape(uy,sizetrg)
    else
        u = lhelmfs(vec(trg),vec(E);derivs=derivs)
        return reshape(u,sizetrg)
    end
end

lhelmfs(trg::Union{VecOrMat{Float64},VecOrMat{Complex{Float64}}},E::Float64;derivs::Bool=false) =
    lhelmfs(trg,fill(E,size(trg));derivs=derivs)

lhelmfs(trg::Union{T1,VecOrMat{T1}},src::Union{T2,VecOrMat{T2}},E::Float64;derivs::Bool=false) where {T1<:Union{Float64,Complex{Float64}},T2<:Union{Float64,Complex{Float64}}} =
    lhelmfs(trg.-src,E.+imag.(src);derivs=derivs)


text:
module PhaseSpaceTools

using Reexport
@reexport using Distributions
@reexport using Statistics
@reexport using LinearAlgebra
@reexport using Parameters

import GSL:sf_laguerre_n

export State, Coherent, Squeezed, Fock, Thermal, Crescent, Bogoliubov
export glauberP, husimiQ, wigner, positiveP, positiveW
export randnc, realnoise, realbridge
export reject, plaguerre, plaguerre_asymptotic, laguerren

include("types.jl")
include("helpers.jl")
include("methods.jl")

end


text:
using Documenter, Mimi

makedocs(
	doctest = false,
    modules = [Mimi],
	sitename = "Mimi.jl",
	pages = [
		"Home" => "index.md",
		"Tutorials" => Any[
			"Tutorials Intro" => "tutorials/tutorial_main.md",
			"1 Install Mimi" => "tutorials/tutorial_1.md",
			"2 Run an Existing Model" => "tutorials/tutorial_2.md",
			"3 Modify an Existing Model" => "tutorials/tutorial_3.md",
			"4 Create a Model" => "tutorials/tutorial_4.md",
			"5 Monte Carlo + Sensitivity Analysis" => "tutorials/tutorial_5.md",
			"6 Create a Model with Composite Components" => "tutorials/tutorial_6.md"
		],
		"How-to Guides" => Any[
			"How-to Guides Intro" => "howto/howto_main.md",
			"1 Construct + Run a Model" => "howto/howto_1.md",
			"2 Explore Results" => "howto/howto_2.md",
			"3 Monte Carlo + SA" => "howto/howto_3.md",
			"4 Timesteps, Params, and Vars" => "howto/howto_4.md",
			"5 Update Time Dimension" => "howto/howto_5.md",
			"6 Port to v0.5.0" => "howto/howto_6.md",
			"7 Port to v1.0.0" => "howto/howto_7.md"
		],
		"Advanced How-to Guides" => Any[
			"Advanced How-to Guides Intro" => "howto_advanced/howto_adv_main.md",
			"Build and Init Functions" => "howto_advanced/howto_adv_buildinit.md",
			"Using Datum References" => "howto_advanced/howto_adv_datumrefs.md"
		],
		"Reference Guides" => Any[
			"Reference Guides Intro" => "ref/ref_main.md",
			"Mimi API" => "ref/ref_API.md",
			"Structures: Classes.jl and Types" => "ref/ref_structures_classes_types.md",
			"Structures: Definitions" => "ref/ref_structures_definitions.md",
			"Structures: Instances" => "ref/ref_structures_instances.md"
		],
		"Explanations" => Any[
			"Explanations Intro" => "explanations/exp_main.md",
			"Models as Packages" => "explanations/exp_pkgs.md"
		],
		"FAQ" => "faq.md",
	],
	format = Documenter.HTML(prettyurls = get(ENV, "JULIA_NO_LOCAL_PRETTY_URLS", nothing) === nothing)
)

deploydocs(
    repo = "github.com/mimiframework/Mimi.jl.git",
)