text:
using Revise

#Question- do residuals of the within model equal residuals of a FI model?
#=using FixedEffectModels, DataFrames
function testwithinresid( N::Int=100, K=10; NG = N ÷ 2K)
  #form the test data
  X = rand(N,K)
  X .*= collect(1:K)'
  Gval = (i->i% NG).(1:N)

  for (i,r) ∈ enumerate(eachrow(X))
    r .+= Gval[i]
  end
  β = sqrt.(collect(1:K))
  y = X*β .+ K*rand(N)

  df = DataFrame(X)
  xnames = propertynames(df)
  df.y = y
  numericnames = propertynames(df)
  df.G = (g->Symbol(:G,g)).(Gval)
  df.i = (i->Symbol(:i,i)).(1:N)

  ffixed = @eval @formula y ~ $(Meta.parse(join(xnames, "+"))) + fe(G)

  dfwi = deepcopy(df)
  dfwis = groupby(df, :G)
  demean(v::AbstractVector{T}) where T<:Real = v .- mean(v)
  demean(v::Any) = v
  dfwi = combine(dfwis, [numericnames...; :i] .=> demean)

  #fix the names
  rename!(dfwi,(n->replace(n, "_demean"=>"")).(names(dfwi)))

  sort!(df, :i)
  sort!(dfwi, :i)

  display(dfwi)

  fwithin = @eval @formula y ~ $(Meta.parse(join(xnames, "+"))) + 0

  rfixed = reg(df, ffixed, save=:residuals)
  rwithin = reg(dfwi, fwithin, save=:residuals)

  df.resid = residuals(rfixed, df)
  dfwi.resid = residuals(rwithin, dfwi)
  println(rfixed)
  println(rwithin)

  @info "Residual difference: $(sum(abs.(df.resid .- dfwi.resid)))"
  return nothing
end

testwithinresid(100,2)=#

#=using CUDA, LinearAlgebra,SparseArrays
function hybridmatrixmultmwe()
  #generate test data
  dX = rand(Float32, 1000,10) |> cu
  S = (x->ifelse(x>0.99f0,1f0,0f0)).(rand(Float32, 1000,1000)) |> SparseMatrixCSC
  dS = S |> CUSPARSE.CuSparseMatrixCSR
  dShyb = dS |> CUSPARSE.switch2hyb

  println("sum of product (CSR):", sum(dS*dX))
  println("sum of product (hybrid):", sum(dShyb*dX))
end
hybridmatrixmultmwe()=#

#=using Revise, CUDA, LinearAlgebra, BenchmarkTools
simplefunc(Aᵢ,Bᵢ) = CUDA.abs(Aᵢ - Bᵢ)
function mweptr(resultsalloc, As, Bs)

  #use the pre-allocation
  results = resultsalloc.results
  Threads.@threads for i ∈ 1:length(As)
    results[i] .= simplefunc.(As[i], Bs[i])
  end

  return reduce(vcat,results)
end

function mweptr(N)
  As::Vector{CuMatrix} = [rand(rand(1:N), 3) |> cu for i ∈ 1:N]
  Bs::Vector{CuMatrix} = [rand(size(As[i],1), 3) |> cu for i ∈ 1:N]
  results::Vector{CuMatrix} = [CUDA.zeros(size(As[i],1), 3) |> cu for i ∈ 1:N]

  resultssum = 0f0
  resultsalloc = (results=results,)
  for i ∈ 1:10
    resultssum += mweptr(resultsalloc, As, Bs) |> sum
  end

  @info "Completed test with sum $(resultssum)"
end

@btime mweptr(500)=#

#=using Revise, LinearAlgebra, CUDA, BenchmarkTools

CUDA.allowscalar(false)
function mwesum(N)
  cpuv= rand(Float32, N)
  print("\nStandard cpu sum: ")
  @btime sum($cpuv)

  cuv=CUDA.cu(cpuv)
  print("Standard cuda sum: ")
  @btime CUDA.@sync sum($cuv)

  print("Summing using allocating a vector of 1s and dot: ")
  onesum(v) = dot(cuv, CUDA.ones(N))
  @btime CUDA.@sync $onesum($cuv)

  cuvones = CUDA.ones(N)
  print("Summing just using dot and a pre-allocated vector of 1s: ")
  onesum2(v) = dot(cuv, cuvones)
  @btime CUDA.@sync $onesum2($cuv)
end

sleep(0.5)
mwesum(256^2)=#

###MWE zygote
using Zygote, CUDA
CUDA.allowscalar(false)

function zygotedevice2host(M::TM, ::Type{T}=eltype(TM)
  ) where {T<:Real, TM<:CuArray}
  buf = Zygote.Buffer(Array{T}(undef, size(M)))
  copyto!(buf, M)
  return copy(buf)
end

function mwe(::Type{T}=Float32; N=5) where T
  x = CUDA.rand(N)

  #verify functions work
  M=CUDA.rand(T, N,N^2)
  v=CUDA.rand(T, N)

  function cuones(::Type{T}, N)
    #v1 = CUDA.rand(T, N) .* T(0.0) .+ T(1.0)
    v1 = ones(N) |> CuVector{T}
  end

  brokenfunctions = (
    shouldwork=()->CUDA.sum(M),
    ones=()->CUDA.ones(N),
    cuones=()->cuones(T, N),
    zeros=()->CUDA.zeros(N),
    fill=()->CUDA.fill(2f0, N),
    cumprod=()->CUDA.cumprod(M, dims=2),
    device2host=()-> rand(Float64, N,N) |> CuMatrix{T} |> zygotedevice2host |>sum,
    hcat=()->CUDA.hcat(M,M),
    hcat2=()->CUDA.hcat(M,cuones(T, N)),
    reduce=()->CUDA.reduce(+,M),
    reduce2=()->CUDA.reduce(hcat,[M,M])
    )

  testcuda(f,x) = CUDA.sum(f() .+ x)
  for n ∈ propertynames(brokenfunctions)
    print("testing $n: ")
    f = brokenfunctions[n]
    #check function works
    @assert testcuda(f,x) > 0
    ∇testcuda(x) = gradient((x)->testcuda(f,x), x)
    try
      ∇testcuda(x) = gradient((x)->testcuda(f,x), x)
      grad = ∇testcuda(x)
      @info "passed with val $grad"
    catch err
      @warn " function $n is broken ($err)"
    end
  end
end
mwe(Float64)

#below copy/paste from https://github.com/FluxML/Zygote.jl/issues/730
#=function mweoriginal(N=5)
  x = CUDA.rand(N)

  M=CUDA.rand(N,N^2)

  brokenfunctions = (
    shouldwork=()->CUDA.sum(M),
    ones=()->CUDA.ones(N),
    zeros=()->CUDA.zeros(N),
    fill=()->CUDA.fill(2f0, N),
    cumprod=()->CUDA.cumprod(M, dims=2),
    device2host=()-> CUDA.rand(N,N) |> Matrix |>sum,
    reduce=()->CUDA.reduce(+,M),
    reduce2=()->CUDA.reduce(hcat,[M,M])
    )

  testcuda(f,x) = CUDA.sum(f() .+ x)
  for n ∈ propertynames(brokenfunctions)
    print("testing $n: ")
    f = brokenfunctions[n]
    #check function works
    @assert testcuda(f,x) > 0
    try
      ∇testcuda(x) = gradient((x)->testcuda(f,x), x)
      grad = ∇testcuda(x)
      @info "passed with val $grad"
    catch err
      @warn " function $n is broken ($err)"
    end
  end
end=#

#mweoriginal()

#=using Revise, Optim
mutable struct AdditionalState
  fevalctr
end

function f(x, Ψ)
  Ψ.fevalctr += 1

  return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

function g!(G, x)
  G[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
  G[2] = 200.0 * (x[2] - x[1]^2)
end

function optimcallback(opt,Ψ)

  @info "fevalctr: $(Ψ.fevalctr)"
  return Ψ.fevalctr > 3
end

function mwe()
  Ψ=AdditionalState(0)

  #prep for the optimization and set parameters
  obj(x)=f(x,Ψ)
  callback(opt) = optimcallback(opt,Ψ)
  optimoptions = Optim.Options(iterations = 10, callback=callback)
  x = zeros(2)
  lower = [-100.,-100.]
  upper = [Inf, Inf]

  res=optimize(obj, g!, lower, upper, x, Fminbox(LBFGS()), optimoptions)
  println(res)
end

mwe()=#

###CUDA tests of prodmat
#=using Revise, CUDA, Zygote, LinearAlgebra, BenchmarkTools, Flux, Transducers, BangBang

prodmatbaseline(M::CuMatrix) = CUDA.hcat((i->CUDA.prod(M[:,1:i], dims=2)).(1:size(M,2))...)
prodmatbaseline(M::Matrix) = hcat((i->prod(M[:,1:i], dims=2)).(1:size(M,2))...)

function prodmattrans(M::TM,
  ::Type{TV} = CuVector{eltype(TM)},
  ::Type{T}=eltype(TM)) where {TM<:CuMatrix,TV<:CuVector, T}

  Nrow,Ncol = size(M)
  cumprodcol(c1)::TV = c1
  cumprodcol(c1,c2)::TV = c1 .* c2

  return foldl(hcat, Scan(cumprodcol), eachcol(M))
  #return reshape(eachcol(M) |> Scan(cumprodcol) |> Cat() |> collect, Nrow, Ncol)
end

function prodmat(M::TM,
  ::Type{T}=eltype(TM)) where {TM<:CuMatrix,T}

  Nrow,Ncol = size(M)

  cumprodcol(c1) = c1
  cumprodcol(c1,c2) = c1 .* c2

  cols = eachcol(M) |> collect
  accumcols = Iterators.accumulate(cumprodcol,cols) |> TCat(1)
  return foldl(hcat, Iterators.accumulate(cumprodcol,cols))
end

function mwe(N=10)

  #verify functions work
  M=CUDA.rand(Float64, N,N^2)
  Mcpu = M |> Matrix

  brokenfunctions = (
    #cumprod=(M)->CUDA.cumprod(M, dims=2),
    #prodmattrans=(M)->prodmattrans(M),
    #prodmat=(M)->prodmat(M),
    prodmatbaseline=(M)->prodmatbaseline(M),
    )

  testcuda(f,M) = CUDA.sum(f(M))
  for n ∈ propertynames(brokenfunctions)
    print("testing $n performance: ")
    f = brokenfunctions[n]
    #check function works
    testval = testcuda(f,M)
    verval = sum(cumprod(M, dims=2))
    @assert (testval ≈ verval) || (testval - N ≈ verval)
    try
      ∇testcuda() = gradient(()->testcuda(f,M), Flux.params(M))
      grad = @btime $∇testcuda()
      @btime $testcuda($f,$M)
      @info "passed with val $(sum(grad[M]))"
    catch err
      @warn " function $n is broken ($err)- performance of call-"
      @btime $testcuda($f,$M)
    end
  end
end

mwe()=#

###regression tests
#=using Revise, CUDA, Zygote, LinearAlgebra, BenchmarkTools, Flux, Transducers, BangBang

reg(X::AbstractMatrix{T},y::Vector{T}, ::Val{:glm}) where T = coef(lm(X, y))
reg(X::AbstractMatrix{T},y::AbstractArray{T}, ::Val{:cholesky}) where T = cholesky!(X'*X)\ (X' * y)
reg(X::AbstractMatrix{T},y::Vector{T}, ::Val{:qr}) where T = qr(X)\y
reg(X::AbstractMatrix{T},y::Vector{T}, ::Val{:svd}) where T = svd(X)\y
reg(X::CuMatrix{T}, y::CuArray{T}, ::Val{:cholesky}) where T = hyperreg(X,y) |> CUDA.vec
reg(X::AbstractMatrix{T}, y::AbstractArray{T}, ::Val{:choleskyzygote}) where T = cholesky(X'*X)\ (X' * y)
reg(X::CuMatrix{T}, y::CuVector{T}, ::Val{:choleskyzygote}) where T = CUDA.cholesky(X'*X)\ (X' * y)
reg(X::AbstractMatrix{T}, y::AbstractArray{T}, ::Val{:pinv}) where T = pinv(X'*X)*(X'*y)
function reg(X::CuMatrix{T}, y::CuArray{T}, ::Val{:pinv}) where T
  #@warn "No reasonable pinv function for CuArrays- falling back to cholesky"
  return reg(X, y, Val{:choleskyzygote}())
end

function mwe(N=10)

  #verify functions work
  X=CUDA.rand(Float64, N*10^3, N)
  y=CUDA.rand(Float64, size(X,1))
  Xcpu = X |> Matrix
  ycpu = y |> Vector

  brokenfunctions = (
    #cumprod=(M)->CUDA.cumprod(M, dims=2),
    reg=(X)->reg(X,y, Val{:choleskyzygote}()),
    #prodmat=(M)->prodmat(M),
    )

  testcuda(f,X) = CUDA.sum(f(X))
  for n ∈ propertynames(brokenfunctions)
    print("testing $n performance: ")
    f = brokenfunctions[n]
    #check function works
    testval = testcuda(f,X)
    verval = (cholesky(Xcpu'*Xcpu)\(Xcpu'*ycpu)) |> sum
    @assert testval ≈ verval
    try
      ∇testcuda() = gradient(()->testcuda(f,X), Flux.params(X))
      grad = @btime $∇testcuda()
      @btime $testcuda($f,$X)
      @info "passed with val $(sum(grad[X]))"
    catch err
      @warn " function $n is broken ($err)- performance of call-"
      @btime $testcuda($f,$X)
    end
  end
end

mwe()=#

#=
###zygote cumprod cpu tests
using Revise, Zygote, LinearAlgebra, BenchmarkTools, Transducers, BangBang, LoopVectorization, CUDA

#prodmatbaseline(M::CuMatrix) = CUDA.hcat((i->CUDA.prod(M[:,1:i], dims=2)).(1:size(M,2))...)
prodmatbaseline(M::Matrix) = hcat((i->prod(M[:,1:i], dims=2)).(1:size(M,2))...)

function prodmattrans(M::TM,
  ::Type{TV} = Vector{eltype(TM)},
  ::Type{T}=eltype(TM)) where {TM<:Matrix,TV<:Vector, T}

  Nrow,Ncol = size(M)
  cumprodcol(c1)::TV = c1
  cumprodcol(c1,c2)::TV = c1 .* c2

  return foldl(hcat, Scan(cumprodcol), eachcol(M))
  #return reshape(eachcol(M) |> Scan(cumprodcol) |> Cat() |> collect, Nrow, Ncol)
end

function prodmat(M::TM,
  ::Type{T}=eltype(TM)) where {TM<:Matrix,T}

  Nrow,Ncol = size(M)
  cumprodcol(c1) = c1
  cumprodcol(c1,c2) = c1 .* c2

  cols = eachcol(M) |> collect
  return foldl(hcat, Iterators.accumulate(cumprodcol,cols))
end

function modprodmat(M::TM, threshold=0.1, ::Type{T}=eltype(TM)) where {TM <: Matrix, T}
  #Mt = M' |> Matrix

  modprod(x) = ifelse(abs(x) > threshold, x, threshold*sign(x))
  function modprod(x,y)
    xy = x*y
    ifelse(abs(xy) > threshold, xy, threshold*sign(xy))
  end

  modprodrow(r) = foldl(modprod, r, init=T(1.0))
  #prodpart(Mpart) = foldl(vcat, modprodrow.(eachrow(Mpart)))
  prodpart(Mpart) = map(r->modprodrow(Mpart[r,:]), 1:size(M,1))
  #cols = (i->prodpart(view(M, :, 1:i))).(1:size(M,2))
  cols = map(i->prodpart(M[:,1:i]), 1:size(M,2))#(i->prodpart(view(M, :, 1:i))).(1:size(M,2))

  return hcat(cols...)
end

#NOTE NOTE NOTE - integrate the below?
#conditiongrowth(g, ::Val{:nozero}, nozerotol=PARAM[:iternozerotol]) = (
#  g .+ nozerotol .* sign.(g) .* exp.(-g .* g))
function boundedcumprod(M::TM, scale::TV = rand(size(M,1)), ::Type{T}=eltype(TM);
   threshold::T=0.1, init=ones(T,size(M,1))) where {TM<:Matrix,TV<:Vector, T}

  #NOTE: more idiomatically, this is x*y > threshold ? x*y : threshold*sign(x*y)
  thresholds = threshold .* scale
  function checkthreshold(x, y, threshold)
    xy = x * y
    xy + threshold * sign(xy) * exp(-xy^2/(2*threshold))
    #ifelse(abs(xy) > threshold, xy, threshold*sign(xy))
  end

  cumprodcol(c1,c2) = vmap(checkthreshold, c1, c2, thresholds)

  cols = eachcol(M) |> collect
  return foldl(hcat, Iterators.accumulate(cumprodcol,cols, init=init))
end
boundedcumprod(M::TM, ::Type{T}=eltype(TM)) where {TM<:CuMatrix,T} = (
  boundedcumprod(M |> Matrix{T}) |> TM)

function mwe(N=25)

  #verify functions work
  M=rand(Float64, N,N^2)# |> CuMatrix{Float64}

  brokenfunctions = (
    #cumprod=(M)->cumprod(M, dims=2),
    #prodmattrans=(M)->prodmattrans(M),
    #modprodmat=(M)->modprodmat(M),
    boundedcumprod=(M)->boundedcumprod(M),
    #prodmatbaseline=(M)->prodmatbaseline(M),
    )

  #(boundedcumprod(M) ≈ modprodmat(M)) || error("boundedcumprod(M) ≈/ modprodmat(M)
  #  Δ: $(mean(abs.(boundedcumprod(M) .- modprodmat(M)), dims=1))
  #  boundedcumprod: $(boundedcumprod(M)[:,1])")
  test(f,M) = sum(f(M))
  for n ∈ propertynames(brokenfunctions)
    sleep(0.1)
    print("testing $n performance: ")
    f = brokenfunctions[n]
    #check function works
    testval = test(f,M)
    #verval = sum(modprodmat(M))
    #(testval ≈ verval) || N > 10 || error("testval ≈/ verval: testval=$testval, verval=$verval")
    try
      ∇test(M) = gradient((M)->test(f,M), M)
      grad = @btime $∇test($M)
      #@btime $test($f,$M)
      println("passed with val $(sum(grad[1]))")
    catch err
      @warn " function $n is broken ($err)- performance of call-"
      @btime $test($f,$M)
    end
  end
end

mwe()
=#

###mcmc tests
#=
using Revise, Distributions, Random, Plots, UnicodePlots, LinearAlgebra

function testacceptreject(N=10^7, ; μ1=-1.0, μ2=2.0, σ1 = 2.0, σ2=3.0)
  unicodeplots()

  σstar = 1/(1/σ1^2+1/σ2^2)
  μstar = σstar*(μ1/σ1^2+μ2/σ2^2)
  @info "μstar: $μstar, σstar^0.5: $(σstar^0.5)"

  #pdf distribution (unkown)
  p(x) = pdf(Normal(μstar, σstar^0.5),x)
  #propp - known to a constant of proportionality, unsamplable
  propp(x) = exp(-(μ1-x)^2/(2σ1^2)-(μ2-x)^2/(2σ2^2))

  q(x) = pdf(Normal(μ2, σ2),x)
  Mq(x) = exp(-(μ2-x)^2/(2σ2^2))

  horiz = rand(Normal(μ2, σ2), N)
  vertu = rand(Uniform(), N)

  propps = (propp).(horiz)
  Mqs = Mq.(horiz)
  @assert all(Mqs .≥ propps)

  ratios = propps ./ Mqs
  keeps = vertu .≤ ratios

  @info "We drew $(sum(keeps)) samples out of a maximum $N"

  ardraws = horiz[keeps]
  @info "Drawn mean: $(mean(ardraws)), stdev = $(std(ardraws))"
  #Plots.histogram(ardraws, normalize=true, bins=-10.:0.01:10.0) |> display
  #Plots.histogram(rand(Normal(μstar, σstar), sum(keeps)), normalize=true, bins=-10.:0.01:10.0) |> display
end

function testacceptreject2(N=10^7, K=7; μ1=-1.0, μ2=2.0, σ1 = 2.0, σ2=3.0)
  unicodeplots()

  X = rand(Normal(0,σ2),N,K)# .+ rand(N)
  Σ2 = cov(X)
  #@info ""
  Λ2 = inv(Σ2)
  Σstar = (diagm(fill(1/σ1^2,K)) .+ Λ2)\I
  μstar = Σstar*(μ1/σ1^2 .+ Λ2*fill(μ2,K))
  @info "μstar: $μstar
    diag(Σ2)^0.5: $(diag(Σ2).^0.5)
    diag(Σstar) ^0.5: $(diag(Σstar).^0.5)"

  #pdf distribution (unkown)
  p(x) = pdf(MvNormal(μstar, Σstar),x)
  #propp - known to a constant of proportionality, unsamplable
  propp(x) = exp(-(μ1.-x)'*(μ1.-x)/(2σ1^2)-(μ2.-x)'*Λ2*(μ2.-x)/2)

  q(x) = pdf(MvNormal(μ2, Σ2),x)
  Mq(x) = exp(-(μ2 .- x)'*Λ2*(μ2 .- x)/2)

  horiz = rand(MvNormal(fill(μ2,K), Σ2), N)' |> Matrix
  vertu = rand(Uniform(), N)

  t0 = time()
    propps = (propp).(eachrow(horiz))
    Mqs = Mq.(eachrow(horiz))
    @assert all(Mqs .≥ propps)

    ratios = propps ./ Mqs
    keeps = (vertu .≤ ratios) |> vec
  Δt = time()-t0

  Nkeeps = sum(keeps)

  @info "We drew $(Nkeeps) samples out of a maximum $N ($(N/Nkeeps) draws per sample,"*
    "$(Δt/Nkeeps) seconds)"

  ardraws = horiz[keeps, :]
  @info "Drawn mean: $(mean(ardraws, dims=1)), stdev = $(std(ardraws, dims=1))"
  #Plots.histogram(ardraws, normalize=true, bins=-10.:0.01:10.0) |> display
  #Plots.histogram(rand(Normal(μstar, σstar), sum(keeps)), normalize=true, bins=-10.:0.01:10.0) |> display
end

#TODO- finish the below. Use a counter to compute rejected draws
#fill a matrix with accepted draws- compare performance with previous
function testacceptreject3(N=10^6, K=7; μ1=1.0, μ2=-2.0, σ1 = 2.0, σ2=3.0)
  unicodeplots()

  X = rand(Normal(0,σ2),N,K)# .+ rand(N)
  Σ2 = cov(X)
  #@info ""
  Λ2 = inv(Σ2)
  Σstar = (diagm(fill(1/σ1^2,K)) .+ Λ2)\I
  μstar = Σstar*(μ1/σ1^2 .+ Λ2*fill(μ2,K))
  @info "μstar: $μstar
    diag(Σstar) ^0.5: $(diag(Σstar).^0.5)"
    #diag(Σ2)^0.5: $(diag(Σ2).^0.5)

  #pdf distribution (unkown)
  p(x) = pdf(MvNormal(μstar, Σstar),x)
  #propp - known to a constant of proportionality, unsamplable
  propp(x) = exp(-(μ1.-x)'*(μ1.-x)/(2σ1^2)-(μ2.-x)'*Λ2*(μ2.-x)/2)

  q(x) = pdf(MvNormal(μ2, Σ2),x)
  Mq(x) = exp(-(μ2 .- x)'*Λ2*(μ2 .- x)/2)

  σk=1 ./ diag(Λ2).^0.5 #Schur complements
  Σk=[Σ2[k:k, 1:K .≠ k] for k in 1:K] #given Σ, selects row vector Σ_k,-k
  function drawqcomponent(x,k)
    #@info "2nd part: $((Σk[k]./Σ2[k,k]) * (x[1:K .≠ k] .- μ2))"
    μ = μ2 - ((Σk[k]./Σ2[k,k]) * (x[1:K .≠ k] .- μ2))[]
    rand(Normal(μ,σk[k]))
  end

  horiz = rand(MvNormal(fill(μ2,K), Σ2), N)' |> Matrix
  vertu = rand(Uniform(), N)

  ctr::Int = 0
  draws::Matrix{Float64} = Matrix{Float64}(undef, N, K)
  seed::Vector{Float64} = ones(K)
  x::Vector{Float64} = seed |> deepcopy
  t0 = time()
  for n ∈ 1:N
    for k ∈ 1:K #for each x component
      while true
        ctr += 1
        x[k] = drawqcomponent(x,k)
        cut = propp(x)/Mq(x)
        (rand(Uniform())<cut) && break #accept if true
      end
    end
    draws[n,:] .= x
  end

  Δt = time()-t0

  @info "We drew $(N) samples after $ctr attempts ($(round(ctr/N, digits=2)) draws per sample).
    Time per sample = $(round(Δt/N, digits=6))"

  @info "Drawn mean: $(mean(draws, dims=1)), stdev = $(std(draws, dims=1))"
  drawupperhalf = draws[(end÷2):end, :]
  @info "Drawn 2nd half mean: $(mean(drawupperhalf, dims=1)), 2nd half stdev = $(std(drawupperhalf, dims=1))"
  #Plots.histogram(ardraws, normalize=true, bins=-10.:0.01:10.0) |> display
  #Plots.histogram(rand(Normal(μstar, σstar), sum(keeps)), normalize=true, bins=-10.:0.01:10.0) |> display
end

function testMHC(N=10^6, K=7; μ1=-1.0, μ2=2.0, σ1 = 2.0, σ2=3.0)
  unicodeplots()

  X = rand(Normal(0,σ2),N,K)# .+ rand(N)
  Σ2 = cov(X)
  #@info ""
  Λ2 = inv(Σ2)
  Σstar = (diagm(fill(1/σ1^2,K)) .+ Λ2)\I
  μstar = Σstar*(μ1/σ1^2 .+ Λ2*fill(μ2,K))
  @info "μstar: $μstar
    diag(Σstar) ^0.5: $(diag(Σstar).^0.5)"
    #diag(Σ2)^0.5: $(diag(Σ2).^0.5)

  #pdf distribution (unkown)
  p(x) = pdf(MvNormal(μstar, Σstar),x)
  #propp - known to a constant of proportionality, unsamplable
  propp(x) = exp(-(μ1.-x)'*(μ1.-x)/(2σ1^2)-(μ2.-x)'*Λ2*(μ2.-x)/2)

  q(x) = pdf(MvNormal(μ2, Σ2),x)
  Mq(x) = exp(-(μ2 .- x)'*Λ2*(μ2 .- x)/2)

  σk=1 ./ diag(Λ2).^0.5 #Schur complements
  Σk=[Σ2[k:k, 1:K .≠ k] for k in 1:K] #given Σ, selects row vector Σ_k,-k
  function drawqcomponent(x,k)
    #@info "2nd part: $((Σk[k]./Σ2[k,k]) * (x[1:K .≠ k] .- μ2))"
    μ = μ2 - ((Σk[k]./Σ2[k,k]) * (x[1:K .≠ k] .- μ2))[]
    d=Normal(μ,σk[k])
    (x)->pdf(d,x), rand(d)
  end

  horiz = rand(MvNormal(fill(μ2,K), Σ2), N)' |> Matrix
  vertu = rand(Uniform(), N)

  ctr::Int = 0
  draws::Matrix{Float64} = Matrix{Float64}(undef, N, K)
  #seed::Vector{Float64} = ones(K)
  x::Vector{Float64} = ones(K)
  xcandidate = x |> deepcopy
  acceptancedraws = rand(Uniform(),N,K)
  t0 = time()
  for n ∈ 1:N
    for k ∈ 1:K #for each x component
      Mf,xcandidate[k] = drawqcomponent(x,k)
      #cut = (propp(xcandidate)/Mq(xcandidate))/(propp(x)/Mq(x)) #NOTE- this can be used as a test
      cut = (propp(xcandidate)/Mf(xcandidate[k]))/(propp(x)/Mf(x[k]))
      x[k] = ifelse(acceptancedraws[n,k]<cut, xcandidate[k], x[k])
      xcandidate[k] = x[k]
    end
    draws[n,:] .= x
  end

  Δt = time()-t0

  @info "We drew $(N) samples after $ctr attempts ($(round(ctr/N, digits=2)) draws per sample).
    Time per sample = $(round(Δt/N, digits=6))"

  @info "Drawn mean: $(mean(draws, dims=1)), stdev = $(std(draws, dims=1))"
  drawupperhalf = draws[(end÷2):end, :]
  @info "Drawn 2nd half mean: $(mean(drawupperhalf, dims=1)), 2nd half stdev = $(std(drawupperhalf, dims=1))"
  @info "uniques: $((c->length(unique(c))).(eachcol(draws)))"
  #Plots.histogram(ardraws, normalize=true, bins=-10.:0.01:10.0) |> display
  #Plots.histogram(rand(Normal(μstar, σstar), sum(keeps)), normalize=true, bins=-10.:0.01:10.0) |> display
end

function testMH(N=10^6, K=7; μ1=-1.0, μ2=2.0, σ1 = 2.0, σ2=3.0)
  unicodeplots()

  X = rand(Normal(0,σ2),N,K)# .+ rand(N)
  Σ2 = cov(X)
  #@info ""
  Λ2 = inv(Σ2)
  Σstar = (diagm(fill(1/σ1^2,K)) .+ Λ2)\I
  μstar = Σstar*(μ1/σ1^2 .+ Λ2*fill(μ2,K))
  @info "μstar: $μstar
    diag(Σstar) ^0.5: $(diag(Σstar).^0.5)"
    #diag(Σ2)^0.5: $(diag(Σ2).^0.5)

  #pdf distribution (unkown)
  p(x) = pdf(MvNormal(μstar, Σstar),x)
  #propp - known to a constant of proportionality, unsamplable
  propp(x) = exp(-(μ1.-x)'*(μ1.-x)/(2σ1^2)-(μ2.-x)'*Λ2*(μ2.-x)/2)

  q(x) = pdf(MvNormal(μ2, Σ2),x)
  Mq(x) = exp(-(μ2 .- x)'*Λ2*(μ2 .- x)/2)

  ctr::Int = 0
  draws::Matrix{Float64} = Matrix{Float64}(undef, N, K)
  drawqs::Matrix{Float64} = rand(MvNormal(fill(μ2,K), Σ2), N)' |> Matrix
  #seed::Vector{Float64} = ones(K)
  x::Vector{Float64} = ones(K)
  xcandidate = x |> deepcopy
  acceptancedraws = rand(Uniform(),N)
  t0 = time()
  for n ∈ 1:N
    xcandidate .= drawqs[n,:]
    cut = (propp(xcandidate)/Mq(xcandidate))/(propp(x)/Mq(x))
    x .= ifelse(acceptancedraws[n]<cut, xcandidate, x)
    draws[n,:] .= x
  end

  Δt = time()-t0

  @info "We drew $(N) samples after $ctr attempts ($(round(ctr/N, digits=2)) draws per sample).
    Time per sample = $(round(Δt/N, digits=6))"

  @info "Drawn mean: $(mean(draws, dims=1)), stdev = $(std(draws, dims=1))"
  drawupperhalf = draws[(end÷2):end, :]
  @info "Drawn 2nd half mean: $(mean(drawupperhalf, dims=1)), 2nd half stdev = $(std(drawupperhalf, dims=1))"
  @info "uniques: $((c->length(unique(c))).(eachcol(draws)))"
  #Plots.histogram(ardraws, normalize=true, bins=-10.:0.01:10.0) |> display
  #Plots.histogram(rand(Normal(μstar, σstar), sum(keeps)), normalize=true, bins=-10.:0.01:10.0) |> display
end

function testacceptreject2abs(N=10^7, K=7; μ1=-1.0, μ2=2.0, σ1 = 2.0, σ2=3.0)
  unicodeplots()

  X = rand(Normal(0,σ2),N,K)# .+ rand(N)
  Σ2 = cov(X)
  #@info ""
  Λ2 = inv(Σ2)
  Σstar = (diagm(fill(1/σ1^2,K)) .+ Λ2)\I
  μstar = Σstar*(μ1/σ1^2 .+ Λ2*fill(μ2,K))
  @info "μstar: $μstar
    diag(Σ2)^0.5: $(diag(Σ2).^0.5)
    diag(Σstar) ^0.5: $(diag(Σstar).^0.5)"

  #pdf distribution (unkown)
  p(x) = pdf(MvNormal(μstar, Σstar),x)
  #propp - known to a constant of proportionality, unsamplable
  propp(x) = exp(-(μ1.-x)'*(μ1.-x)/(2σ1^2)-(μ2.-abs.(x))'*Λ2*(μ2.-abs.(x))/2)

  #q(x) = pdf(MvNormal(μ2, Σ2),x)
  #Mq(x) = exp(-(μ2 .- x)'*Λ2*(μ2 .- x)/2)
  q(x) = rand(Uniform(μ2-5*σ2, μ2+5*σ2), K)
  Mq(x)=1.0

  horiz = rand(Uniform(μ2-5*σ2, μ2+5*σ2), K, N)' |> Matrix
  vertu = rand(Uniform(), N)

  t0 = time()
    propps = (propp).(eachrow(horiz))
    Mqs = Mq.(eachrow(horiz))
    @assert all(Mqs .≥ propps)

    ratios = propps ./ Mqs
    keeps = (vertu .≤ ratios) |> vec
  Δt = time()-t0

  Nkeeps = sum(keeps)

  @info "We drew $(Nkeeps) samples out of a maximum $N ($(N/Nkeeps) draws per sample,"*
    "$(Δt/Nkeeps) seconds)"

  ardraws = horiz[keeps, :]
  @info "Drawn mean: $(mean(ardraws, dims=1)), stdev = $(std(ardraws, dims=1))"
  #Plots.histogram(ardraws, normalize=true, bins=-10.:0.01:10.0) |> display
  #Plots.histogram(rand(Normal(μstar, σstar), sum(keeps)), normalize=true, bins=-10.:0.01:10.0) |> display
end

function testMHCabs(N=10^6, K=7; μ1=-1.0, μ2=2.0, σ1 = 2.0, σ2=3.0)
  unicodeplots()

  X = rand(Normal(0,σ2),N,K)# .+ rand(N)
  Σ2 = cov(X)
  #@info ""
  Λ2 = inv(Σ2)
  Σstar = (diagm(fill(1/σ1^2,K)) .+ Λ2)\I
  μstar = Σstar*(μ1/σ1^2 .+ Λ2*fill(μ2,K))
  @info "μstar: $μstar
    diag(Σstar) ^0.5: $(diag(Σstar).^0.5)"
    #diag(Σ2)^0.5: $(diag(Σ2).^0.5)

  #pdf distribution (unkown)
  p(x) = pdf(MvNormal(μstar, Σstar),x)
  #propp - known to a constant of proportionality, unsamplable
  propp(x) = exp(-(μ1.-x)'*(μ1.-x)/(2σ1^2)-(μ2.-abs.(x))'*Λ2*(μ2.-abs.(x))/2)
  #propp(x) = exp(-(μ1.-x)'*(μ1.-x)/(2σ1^2)-(μ2.-x)'*Λ2*(μ2.-x)/2)

  q(x) = pdf(MvNormal(μ2, Σ2),x)
  Mq(x) = exp(-(μ2 .- x)'*Λ2*(μ2 .- x)/2)

  σk=1 ./ diag(Λ2).^0.5 #Schur complements
  Σk=[Σ2[k:k, 1:K .≠ k] for k in 1:K] #given Σ, selects row vector Σ_k,-k
  ΣMkinv = (k->Σ2[1:K .≠ k, 1:K .≠ k]\I).(1:K)
  function drawqcomponent(x,k)
    #@info "2nd part: $((Σk[k]./Σ2[k,k]) * (x[1:K .≠ k] .- μ2))"
    #μ = μ2 - ((Σk[k]./Σ2[k,k]) * (x[1:K .≠ k] .- μ2))[]
    μ = μ2 - ((Σk[k]*ΣMkinv[k]) * (x[1:K .≠ k] .- μ2))[]
    #d=Normal(μ,σk[k])
    d = Normal(-40.,3.)
    (x)->pdf(d,x), rand(d)
  end

  horiz = rand(MvNormal(fill(μ2,K), Σ2), N)' |> Matrix
  vertu = rand(Uniform(), N)

  ctr::Int = 0
  draws::Matrix{Float64} = Matrix{Float64}(undef, N, K)
  #seed::Vector{Float64} = ones(K)
  x::Vector{Float64} = ones(K) .* -39.0
  xcandidate = x |> deepcopy
  acceptancedraws = rand(Uniform(),N,K)
  t0 = time()
  for n ∈ 1:N
    for k ∈ 1:K #for each x component
      Mf,xcandidate[k] = drawqcomponent(x,k)
      #cut = (propp(xcandidate)/Mq(xcandidate))/(propp(x)/Mq(x)) #NOTE- this can be used as a test
      cut = (propp(xcandidate)/Mf(xcandidate[k]))/(propp(x)/Mf(x[k]))
      x[k] = ifelse(acceptancedraws[n,k]<cut, xcandidate[k], x[k])
      xcandidate[k] = x[k]
    end
    draws[n,:] .= x
  end

  Δt = time()-t0

  @info "We drew $(N) samples after $ctr attempts ($(round(ctr/N, digits=2)) draws per sample).
    Time per sample = $(round(Δt/N, digits=6))"

  @info "Drawn mean: $(mean(draws, dims=1)), stdev = $(std(draws, dims=1))"
  drawupperhalf = draws[(end÷2):end, :]
  @info "Drawn 2nd half mean: $(mean(drawupperhalf, dims=1)), 2nd half stdev = $(std(drawupperhalf, dims=1))"
  @info "uniques: $((c->length(unique(c))).(eachcol(draws)))"
  #Plots.histogram(ardraws, normalize=true, bins=-10.:0.01:10.0) |> display
  #Plots.histogram(rand(Normal(μstar, σstar), sum(keeps)), normalize=true, bins=-10.:0.01:10.0) |> display
end

#NOTE: OK
#testacceptreject(10^6) #1D only
#testacceptreject2(10^6,1) #slow but looks good- might be a nice check for few vars
#testMHC(10^5,3) #trust this a little more, looks great!

testacceptreject2abs(10^7,3)
testMHCabs(10^7,3)

#rejected
#testacceptreject3(10^6,10) #slow but looks good- might be a nice check for few vars
#testMH(10^4,10) #slow anbd not very good
=#

###regression zygote
#=using Zygote, CUDA, LinearAlgebra, Distributions, StatsBase, Statistics, GLM
CUDA.allowscalar(false)

zygotedevice2host(M::TM, ::Type{T}=eltype(TM)
  ) where {T<:Real, TM<:CuArray} = copyto!(Zygote.Buffer(Array{T}(undef, size(M))), M) |> copy

function regtest(::Type{TM}, ::Type{TV}, ::Type{T}; N=100, K=5) where {TM, TV, T}

  y = rand(N^2) |> TV
  A = rand(N,K)
  idx = rand([1:N...;],N*N)

  function runreg(A, f)
    X = A[idx, :] |> TM
    b=f(X,y)
    return sum((X*b .- y).^2)
  end

  regchol(X,y) = cholesky(X'*X)\(X'*y)
  regsvd(X,y) = svd(X)\y
  regqr(X,y) = qr(X'*X)\(X'*y)
  regpinv(X,y) = pinv(X'*X)*(X'*y)
  reglm(X,y) = coef(lm(X,y))

  regfunctions = (
    shouldwork=(A)->runreg(A, regchol),
    #svd=(A)->runreg(A, regsvd),
    #qr=(A)->runreg(A, regqr),
    pinv=(A)->runreg(A, regpinv),
    lm=(A)->runreg(A, reglm),
    )

  for n ∈ propertynames(regfunctions)
    print("\ntesting $n: ")
    f = regfunctions[n]
    #check function works

    try
      @assert f(A) > 0
      ∇test(A) = gradient((A)->f(A), A)
      grad = ∇test(A)[1]

      @info("passed with $(sum(grad .≠ 0.0)) non-zero gradient values")
    catch err
      errmessage = "$err"
      if length(errmessage) > 1000
        errmessage=errmessage[1:1000]
      end
      @warn " function $n is broken ($errmessage)"
    end
  end
end

#regtest(CuMatrix{Float64}, CuVector{Float64}, Float64, N=10, K=5)
regtest(Matrix{Float64}, Vector{Float64}, Float64, N=10, K=5)=#

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using DelimitedFiles: readdlm

read_data(fn::AbstractString) = vec(readdlm(fn, skipstart=1))
read_cpts(fn::AbstractString) = vec([0; readdlm(fn, Int, skipstart=1)])

dir = "Test_Files"

@testset "Test Normal mean change..." begin
    data = read_data("$(dir)/Normal_change_in_mean_data")

    r_cpts = read_cpts("$(dir)/Normal_change_in_mean_cpts")
    cpts, cost = @PELT data Normal(:?, 1.0) 2*log(length(data))

    @test r_cpts == cpts
end

# Normal var change
@testset "Test Normal variance change..." begin
    data = read_data("$(dir)/Normal_change_in_var_data")

    r_cpts = read_cpts("$(dir)/Normal_change_in_var_cpts")
    cpts, cost = @PELT data Normal(1.0, :?) 2*log(length(data))

    @test r_cpts == cpts
end

# Normal meanvar change
@testset "Test Normal mean and variance chaange..." begin
    data = read_data("$(dir)/Normal_change_in_meanvar_data")

    r_cpts = read_cpts("$(dir)/Normal_change_in_meanvar_cpts")
    cpts, cost = @PELT data Normal(:?, :?) 4*log(length(data))

    @test r_cpts == cpts
end

# Poisson change
@testset "Test Poisson change..." begin
    data = read_data("$(dir)/Poisson_data")

    r_cpts = read_cpts("$(dir)/Poisson_cpts")
    cpts, cost = @PELT data Poisson(:?) 2*log(length(data))

    @test r_cpts == cpts
end

# Exponential change
@testset "Test Exponential  change..." begin
    data = read_data("$(dir)/Exponential_data")

    r_cpts = read_cpts("$(dir)/Exponential_change_in_mean_cpts")
    cpts, cost = @PELT data Exponential(:?) 2*log(length(data))

    @test r_cpts == cpts
end

# Gamma rate changes
@testset "Test Gamma rate change..." begin
    data = read_data("$(dir)/Gamma_data")

    r_cpts = read_cpts("$(dir)/Gamma_cpts")
    cpts, cost = @PELT data Gamma(1.0,:?) 2*log(length(data))

    @test r_cpts == cpts
end

@testset "Test Nonparametric change..." begin
    data = read_data("$(dir)/Nonparametric_data")

    r_cpts = read_cpts("$(dir)/Nonparametric_cpts")
    cpts, cost = @PELT data Nonparametric(10) 4*log(length(data))

    @test r_cpts == cpts
end

ext: jl


text:
using Test
using LinearAlgebra
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary:
    prob_ode_linear, prob_ode_2Dlinear, prob_ode_lotkavoltera, prob_ode_fitzhughnagumo

prob = prob_ode_lotkavoltera
prob = remake(prob, tspan=(0.0, 10.0))
prob = ProbNumDiffEq.remake_prob_with_jac(prob)

@testset "Condition numbers of A,Q" begin
    h = 0.1 * rand()
    σ = rand()

    d, q = 2, 3

    A!, Q! = ProbNumDiffEq.vanilla_ibm(d, q)
    Ah = diagm(0 => ones(d * (q + 1)))
    Qh = zeros(d * (q + 1), d * (q + 1))
    A!(Ah, h)
    Q!(Qh, h, σ^2)

    A_p, Q_p = ProbNumDiffEq.ibm(d, q)
    Ah_p = A_p
    Qh_p = Q_p * σ^2

    # First test that they're both equivalent
    D = d * (q + 1)
    P, PI = ProbNumDiffEq.init_preconditioner(d, q)
    ProbNumDiffEq.make_preconditioner!(P, h, d, q)
    ProbNumDiffEq.make_preconditioner_inv!(PI, h, d, q)
    @test Qh_p ≈ P * Qh * P'
    @test Ah_p ≈ P * Ah * PI

    # Check that the preconditioning actually helps
    # @info "Condition numbers" cond(Qh) cond(Matrix(Qh_p)) cond(Ah) cond(Ah_p)
    @test cond(Qh) > cond(Matrix(Qh_p))
    @test cond(Qh) > cond(Matrix(Qh_p))^2
end


text:
using JuMP
using Hungarian
using LinearAlgebra
import MathOptInterface as MOI

"""
Approximately solve the QAP by linearizing the problem and solving
the resulting linear integer program.
"""
function qap_linearization(A, B, obj)
	P, _ = adams_johnson_linearization(A, B, obj, integer=true)
	matching = [findfirst(==(1), P[i,:]) for i in 1:size(A,1)]
	return P, matching
end

"""
Approximately solve the QAP by linearizing and relaxing the problem
and solving the resulting LIP with linear programming and rounding
via the Hungarian algorithm.
"""
function qap_lprounding(A, B, obj)
	x, _ = adams_johnson_linearization(A, B, obj, integer=false)

	W = maximum(x) .- x
	P = Hungarian.munkres(W)

	P = Matrix(P .== Hungarian.STAR)
	matching = [findfirst(P[i,:]) for i in 1:size(A,1)]

	return P, matching
end

"""
QAP linearization from "Improved linear programming-based
lower bounds for the quadratic assignment problem" by Adams and Johnson.
"""
function adams_johnson_linearization(A, B, obj; integer::Bool=true)
	N = size(A, 1)

	optimizer = get_optimizer()
	model = Model(optimizer)
	set_silent(model)

	@variable(model, x[1:N,1:N], binary=integer)
	@variable(model, y[1:N,1:N,1:N,1:N], binary=integer)

	@constraint(model, [u=1:N], sum(x[u,:]) == 1)
	@constraint(model, [p=1:N], sum(x[:,p]) == 1)
	@constraint(model, [v=1:N,p=1:N,q=1:N], sum(y[:,p,v,q]) == x[v,q])
	@constraint(model, [u=1:N,v=1:N,q=1:N], sum(y[u,:,v,q]) == x[v,q])
	@constraint(model, [u=1:N,p=1:N,v=1:N,q=1:N], y[u,p,v,q] == y[v,q,u,p])

	if !integer
		@constraint(model, [u=1:N,p=1:N], 0 ≤ x[u,p] ≤ 1)
		@constraint(model, [u=1:N,p=1:N,v=1:N,q=1:N], 0 ≤ y[u,p,v,q] ≤ 1)
	end

	model_obj = (obj == :min) ? MOI.MIN_SENSE : MOI.MAX_SENSE
	@objective(model, model_obj, sum(A[u,v] * B[p,q] * y[u,p,v,q] for u=1:N, p=1:N, v=1:N, q=1:N))

	optimize!(model)

	status = Int(termination_status(model))
	@assert status == 1

	x_sol = value.(x)
	y_sol = value.(y)

	clamp!(x_sol, 0, 1)
	clamp!(y_sol, 0, 1)

	if integer
		x_sol = x_sol .> 0.5
		y_sol = y_sol .> 0.5
	end

	return x_sol, y_sol
end


text:
include("../src/WCSLIB.jl")
using WCSLIB
using .Test

w = wcsprm(2; alt     = 'B',
              cdelt   = [-0.066667, 0.066667],
              colnum  = 2,
              ctype   = ["RA---AIR", "DEC--AIR"],
              crpix   = [-234.75, 8.3393],
              crval   = [0., -90],
              obsgeo  = [1., 2, 3],
              pc      = [1. 0; 0 1],
              pv      = [pvcard(2, 1, 45.0)],
              velangl = 3,
              wcsname = "WCSLIB.jl")
@test wcsset(w) == 0

pixcrd = [0. 0; 24 38; 45 98]'
world = wcsp2s(w, pixcrd)
expected_world = [ 267.96547027  -73.73660749;
 276.53931377  -71.97412809;
 287.77080792  -69.67813884]'
@test maximum(abs(world .- expected_world)) < 5e-9
pixcrd2 = wcss2p(w, world)
@test maximum(abs(pixcrd .- pixcrd2)) < 1e-9
wcsfree(w)

header = "SIMPLE  =                    T / file does conform to FITS standard             BITPIX  =                  -64 / number of bits per data pixel                  NAXIS   =                    2 / number of data axes                            NAXIS1  =                 3636 / length of data axis 1                          NAXIS2  =                 1939 / length of data axis 2                          EXTEND  =                    T / FITS dataset may contain extensions            COMMENT   FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H WCSAXES =                    2 / Number of coordinate axes                      CRPIX1  =               1818.0 / Pixel coordinate of reference point            CRPIX2  =                969.5 / Pixel coordinate of reference point            CDELT1  =             -0.00825 / [deg] Coordinate increment at reference point  CDELT2  =              0.00825 / [deg] Coordinate increment at reference point  CUNIT1  = 'deg'                / Units of coordinate increment and value        CUNIT2  = 'deg'                / Units of coordinate increment and value        CTYPE1  = 'RA---CEA'           / Right ascension, cylindrical equal area projectCTYPE2  = 'DEC--CEA'           / Declination, cylindrical equal area projection CRVAL1  =                357.5 / [deg] Coordinate value at reference point      CRVAL2  =                  0.0 / [deg] Coordinate value at reference point      PV2_1   =                  1.0 / CEA projection parameter                       LONPOLE =                  0.0 / [deg] Native longitude of celestial pole       LATPOLE =                 90.0 / [deg] Native latitude of celestial pole        RADESYS = 'ICRS'               / Equatorial coordinate system                   RESTFRQ =                    0 / [Hz] Line rest frequency                       RESTWAV =                    0 / [Hz] Line rest wavelength                      END                                                                             "

(ws,nreject) = wcspih(header)
@test length(ws) == 1
@test eltype(ws) <: wcsprm
@test nreject == 0
w = ws[1]
@test wcsset(w) == 0
@test w.naxis == 2

header2 = wcshdo(w)
@test header[641:end-240] == header2

wcsfree(w)


text:
# Calculate the integral of s=0 to t of exp(T(t-s)) beta alpha exp(T s) ds
# using the uniformized version of T (specified by P and r). See equation (44)
# of Okamura et al (2012), Improvement of expectation-maximization algorithm
# for phase-type distributions with grouped and truncated data.
function conv_int_unif(r, P, t, beta, alpha)
    p = size(P, 1)
    ϵ = 1e-3
    R = quantile(Poisson(r * t), 1-ϵ)

    betas = Array{Float64}(undef, p, R+1)
    betas[:,1] = beta

    for u = 1:R
        betas[:,u+1] = P * betas[:,u]
    end

    poissPDFs = pdf.(Poisson(r*t), 1:R+1)

    alphas = Array{Float64}(undef, p, R+1)
    alphas[:, R+1] = poissPDFs[R+1] .* alpha'

    for u = (R-1):-1:0
        alphas[:, u+1] = alphas[:, u+2]' * P + poissPDFs[u+1] .* alpha'
    end

    Υ = zeros(p, p)
    for u = 0:R
        Υ += betas[:, u+1] * alphas[:, u+1]' ./ r
    end

    Υ
end

function e_step_observed_uniform!(s::Sample, fit::PhaseType,
        Bs::AbstractArray{Float64}, Zs::AbstractArray{Float64},
        Ns::AbstractArray{Float64})
    p = fit.p

    K = length(s.obs)
    deltaTs = diff([0; s.obs])

    fs = zeros(p, K+1)
    bs = zeros(p, K+1)
    fs[:,1] = fit.π
    bs[:,1] = fit.t

    for k = 1:K
        expTDeltaT = exp(fit.T * deltaTs[k])
        fs[:,k+1] = fs[:,k]' * expTDeltaT
        bs[:,k+1] = expTDeltaT * bs[:,k]
    end

    cs = zeros(p, K)
    cs[:,K] = fit.π / dot(fit.π, bs[:,K+1])
    for k = (K-1):-1:1
        cs[:,k] = cs[:,k+1]' * exp(fit.T * deltaTs[k+1]) +
            fit.π' / dot(fit.π, bs[:,k+1])
    end

    H = zeros(p, p)
    r = 1.01 * maximum(abs.(diag(fit.T)))
    P = I + (fit.T ./ r)

    for k = 1:K
        H += s.obsweight[k] .* conv_int_unif(r, P, deltaTs[k], bs[:,k], cs[:,k])
    end

    for k = 1:K
        Bs[:] = Bs[:] + s.obsweight[k] .* (fit.π .* bs[:,k+1]) ./ dot(fit.π, bs[:,k+1])
        Ns[:,end] = Ns[:,end] + s.obsweight[k] .* (fs[:,k+1] .* fit.t) ./ dot(fit.π, bs[:,k+1])
    end

    Zs[:] = Zs[:] + diag(H)
    Ns[:,1:p] = Ns[:,1:p] + fit.T .* (H') .* (1 .- Matrix{Float64}(I, p, p))
end

function e_step_censored_uniform!(s::Sample, fit::PhaseType,
        Bs::AbstractArray{Float64}, Zs::AbstractArray{Float64},
        Ns::AbstractArray{Float64})
    p = fit.p

    int = s.int
    intweight = s.intweight

    # For some reason, we need an interval which starts at 0, so we just
    # add one in with zero weight attached to it.
    if length(s.int) > 0 && minimum(s.int) > 0
        int = vcat([0, minimum(int)]', int)
        intweight = vcat(0, intweight)
    end

    K = size(int, 1)
    deltaTs = int[:,2] - int[:,1]

    barfs = zeros(p, K+1)
    tildefs = zeros(p, K)
    barbs = zeros(p, K+1)
    tildebs = zeros(p, K)
    ws = zeros(K+1)
    N = sum(intweight)
    U = 0

    barfs[:,1] = fit.π' * inv(-fit.T)
    barbs[:,1] = ones(p)

    for k = 1:K
        expTDeltaT = exp(fit.T * deltaTs[k])

        barfs[:,k+1] = barfs[:,k]' * expTDeltaT
        tildefs[:,k] = barfs[:,k] - barfs[:,k+1]
        barbs[:,k+1] = expTDeltaT * barbs[:,k]
        tildebs[:,k] = barbs[:,k] - barbs[:,k+1]

        U += fit.π' * tildebs[:,k]

        ws[k] = intweight[k] / (fit.π' * tildebs[:,k])
    end
    U += fit.π' * barbs[:,K+1]
    ws[K+1] = 0

    cs = zeros(p, K)
    cs[:,K] = (ws[K+1] - ws[K]) .* fit.π'
    for k = (K-1):-1:1
        cs[:,k] = cs[:,k+1]' * exp(fit.T * deltaTs[k+1]) +
            (ws[k+1] - ws[k]) .* fit.π'
    end

    H = zeros(p, p)
    r = 1.01 * maximum(abs.(diag(fit.T)))
    P = I + (fit.T ./ r)

    for k = 1:K
        H += ws[k] .* ones(p) * tildefs[:,k]' +
            conv_int_unif(r, P, deltaTs[k], barbs[:,k], cs[:,k])
    end
    H += ws[K+1] .* ones(p) * barfs[:,K+1]'

    # Step 4
    for k = 1:K
        Bs[:] = Bs[:] + intweight[k] .* (fit.π .* tildebs[:,k]) ./
            (fit.π' * tildebs[:,k])
        Ns[:,end] = Ns[:,end] + intweight[k] .* (tildefs[:,k] .* fit.t) ./
            (tildefs[:,k]' * fit.t)
    end

    Zs[:] = Zs[:] + diag(H)
    Ns[:,1:p] = Ns[:,1:p] + fit.T .* (H') .* (1 .- Matrix{Float64}(I, p, p))
end


text:
x = range(0, stop=π/2, length=100)
d1 = plot(x, sin.(x), lw=0.05)
d2 = plot!(x, cos.(x), lw=0.05)
fill_between!(x, sin.(x), cos.(x), alpha=0.2)
x = x[1:5:end]
d3 = scatter!(x, sin.(x) .* cos.(x), msize=0.3, marker="fdiamond")
legend([d1, d2, d3], ["sine", "cosine", "scatter"])
xlim(0,π/2)


text:
using OpenQASM.Types
using OpenQASM.RBNF: Token

#todo support custom yao gates to qasm

"""
    convert_to_qasm(qc, ncreg)

Parses an `AbstractBlock` into code based on the `OpenQASM` spec

- `qc`: A `ChainBlock`(circuit that is to be run).
- `ncreg` (optional) : Number of classical registers
While performing operations like measuring, one can input desired number of classical regs(each size equal to number of qubits). Defaults to 1.
"""
function convert_to_qasm(qc::AbstractBlock{N}, ncreg::Int = 1) where {N}
    prog = generate_prog(qc, ncreg)
    MainProgram(v"2", prog)
end

"""
    generate_prog(qc)

Parses the YaoIR into a list of `QASM` instructions

- `qc`: A `ChainBlock`(circuit that is to be run).
"""
function generate_prog(qc::AbstractBlock{N}, ncreg::Int) where {N}
    prog = []
    generate_defs(prog, N, ncreg)
    generate_prog!(prog, basicstyle(qc), [0:N-1...], Int[])
    return prog
end

function generate_defs(prog, nq::Int, ncreg::Int)
    qregs = RegDecl(Token{:reserved}("qreg"), Token{:type}("q"), Token{:int}("$nq"))
    cregs = collect(
        RegDecl(Token{:reserved}("creg"), Token{:type}("c$i"), Token{:int}("$nq")) for
        i = 1:ncreg
    )
    push!(prog, Include(Bit("\"qelib1.inc\"")), qregs, cregs...)
end

function generate_prog!(prog, qc_simpl::ChainBlock, locs, controls)
    for block in subblocks(qc_simpl)
        generate_prog!(prog, block, locs, controls)
    end
end

function generate_prog!(prog, blk::PutBlock{N,M}, locs, controls) where {N,M}
    generate_prog!(prog, blk.content, sublocs(blk.locs, locs), controls)
end

function generate_prog!(prog, blk::ControlBlock{N,GT,C}, locs, controls) where {N,GT,C}
    any(==(0), blk.ctrl_config) && error("Inverse Control used in Control gate context")
    generate_prog!(
        prog,
        blk.content,
        sublocs(blk.locs, locs),
        [controls..., sublocs(blk.ctrl_locs, locs)...],
    )
end

function generate_prog!(prog, m::YaoBlocks.Measure{N}, locs, controls) where {N}
    mlocs = sublocs(m.locations isa AllLocs ? [1:N...] : [m.locations...], locs)
    (length(controls) == 0) || error("controlled measure is not supported")

    if !(m.operator isa ComputationalBasis) && m.operator isa AbstractBlock
        generate_prog!(prog, m.operator, mlocs, ())
    end

    # can be improved
    for i in mlocs
        cargs = Bit("c1", i)
        qargs = Bit("q", i)
        inst = Types.Measure(qargs, cargs)
        push!(prog, inst)
    end
end

# IBMQ Chip only supports ["id", "u1", "u2", "u3", "cx"]
# x, y, z and control x, y, z, id, t, swap and other primitive gates
for (GT, NAME, MAXC) in [
    (:XGate, "x", 2),
    (:YGate, "y", 2),
    (:ZGate, "z", 2),
    (:I2Gate, "id", 0),
    (:TGate, "t", 0),
    (:SWAPGate, "swap", 0),
]
    @eval function generate_prog!(prog, ::$GT, locs, controls)
        if length(controls) <= $MAXC
            # push!(prog, )
            name = "c"^(length(controls)) * $NAME
            if name != "cx" && name != "swap"
                cargs = [] #empty for now
                qargs =
                    isempty(controls) ? [Bit("q", locs...)] :
                    [Bit("q", controls...), Bit("q", locs...)]
                inst = Instruction(name, cargs, qargs)
            elseif name == "swap"
                cargs = [] #empty for now
                qargs = [Bit("q", i) for i in locs]
                inst = Instruction(name, cargs, qargs)
            else
                inst = CXGate(Bit("q", controls...), Bit("q", locs...))
            end
            push!(prog, inst)
        else
            error("too many control bits!")
        end
    end
end

for (GT, NAME, PARAMS, MAXC) in [
    (:(RotationGate{1,T,XGate} where {T}), "rx", :(b.theta), 0),
    (:(RotationGate{1,T,YGate} where {T}), "ry", :(b.theta), 0),
    (:(RotationGate{1,T,ZGate} where {T}), "rz", :(b.theta), 0),
    (:(ShiftGate), "p", :(b.theta), 1),
    (:(HGate), "h", :(nothing), 0),
]
    @eval function generate_prog!(prog, b::$GT, locs, controls)
        if length(controls) <= $MAXC
            name = "c"^(length(controls)) * $NAME
            if $PARAMS === nothing
                cargs = []
                qargs = [Bit("q", locs...)]
            else
                params = $PARAMS
                cargs = [Token{:float64}("$params")]
                qargs =
                    isempty(controls) ? [Bit("q", locs...)] :
                    [Bit("q", controls...), Bit("q", locs...)]
            end
            push!(prog, Instruction(name, cargs, qargs))
        else
            error("too many control bits! got $controls (length > $($(MAXC)))")
        end
    end
end

sublocs(subs, locs) = [locs[i] for i in subs]

function basicstyle(blk::AbstractBlock)
    YaoBlocks.Optimise.simplify(blk, rules = [YaoBlocks.Optimise.to_basictypes])
end

max_stars_repo_name: Sov-trotter/YaoBlocksQASM

max_issues_repo_name: Sov-trotter/YaoBlocksQASM

max_forks_repo_name: Sov-trotter/YaoBlocksQASM


text:
using HTTP
using HTTP.Test
using HTTP.Base64
using JSON
using MbedTLS: digest, MD_MD5, MD_SHA256

using HTTP.IOExtras
using HTTP.request

println("async tests")

stop_pool_dump = false

@async HTTP.listen() do http
    startwrite(http)
    write(http, """

            <title>HTTP.jl Connection Pool</title>
            <meta http-equiv="refresh" content="1">

    """)
    write(http, "
")
    buf = IOBuffer()
    HTTP.ConnectionPool.showpoolhtml(buf)
    write(http, take!(buf))
    write(http, "
")
end

@async begin
    sleep(5)
    try
        run(`open http://localhost:8081`)
    catch e
        while !stop_pool_dump
            HTTP.ConnectionPool.showpool(STDOUT)
            sleep(1)
        end
    end
end

# Tiny S3 interface...
s3region = "ap-southeast-2"
s3url = "https://s3.$s3region.amazonaws.com"
#s3(method, path, body=UInt8[]; kw...) =
#    request(method, "$s3url/$path", [], body; aws_authorization=true, kw...)
#s3get(path; kw...) = s3("GET", path; kw...)
#s3put(path, data; kw...) = s3("PUT", path, data; kw...)

#=
function create_bucket(bucket)
    s3put(bucket, """
        <CreateBucketConfiguration
                     xmlns="http://s3.amazonaws.com/doc/2006-03-01/">
            <LocationConstraint>$s3region</LocationConstraint>
        </CreateBucketConfiguration>""",
        statusexception=false)
end

create_bucket("http.jl.test")
=#

function dump_async_exception(e, st)
    buf = IOBuffer()
    write(buf, "==========\n@async exception:\n==========\n")
    show(buf, "text/plain", e)
    show(buf, "text/plain", st)
    write(buf, "==========\n\n")
    print(String(take!(buf)))
end

if haskey(ENV, "AWS_ACCESS_KEY_ID") ||
   (VERSION > v"0.7.0-DEV.2338" && haskey(ENV, "AWS_DEFAULT_PROFILE"))
@testset "async s3 dup$dup, count$count, sz$sz, pipw$pipe, $http, $mode" for
    count in [10, 100, 1000],
    dup in [0, 7],
    http in ["http", "https"],
    sz in [100, 10000],
    mode in [:request, :open],
    pipe in [0, 32]

if (dup == 0 || pipe == 0) && count > 100
    continue
end

global s3url
s3url = "$http://s3.$s3region.amazonaws.com"
println("running async s3 dup$dup, count$count, sz$sz, pipe$pipe, $http, $mode")

put_data_sums = Dict()
ch = 100
conf = [:reuse_limit => 90,
        :verbose => 0,
        :pipeline_limit => pipe,
        :connection_limit => dup + 1,
        :readtimeout => 120]

@sync for i = 1:count
    data = rand(UInt8(65):UInt8(75), sz)
    md5 = bytes2hex(digest(MD_MD5, data))
    put_data_sums[i] = md5
    @async try
        url = "$s3url/http.jl.test/file$i"
        r = nothing
        if mode == :open
            r = HTTP.open("PUT", url, ["Content-Length" => sz];
                    body_sha256=digest(MD_SHA256, data),
                    body_md5=digest(MD_MD5, data),
                    aws_authorization=true,
                    conf...) do http
                for n = 1:ch:sz
                    write(http, data[n:n+(ch-1)])
                    sleep(rand(1:10)/1000)
                end
            end
        end
        if mode == :request
            r = HTTP.request("PUT", url, [], data;
                    aws_authorization=true, conf...)
        end
        #println("S3 put file$i")
        @assert strip(HTTP.header(r, "ETag"), '"') == md5
    catch e
        dump_async_exception(e, catch_stacktrace())
        rethrow(e)
    end
end

get_data_sums = Dict()
@sync begin
    for i = 1:count
        @async try
            url = "$s3url/http.jl.test/file$i"
            buf = BufferStream()
            r = nothing
            if mode == :open
                r = HTTP.open("GET", url, ["Content-Length" => 0];
                              aws_authorization=true,
                              conf...) do http
                    buf = BufferStream() # in case of retry!
                    while !eof(http)
                        write(buf, readavailable(http))
                        sleep(rand(1:10)/1000)
                    end
                    close(buf)
                end
            end
            if mode == :request
                r = HTTP.request("GET", url; response_stream=buf,
                                             aws_authorization=true, conf...)
            end
            #println("S3 get file$i")
            bytes = read(buf)
            md5 = bytes2hex(digest(MD_MD5, bytes))
            get_data_sums[i] = (md5, strip(HTTP.header(r, "ETag"), '"'))
        catch e
            dump_async_exception(e, catch_stacktrace())
            rethrow(e)
        end
    end
end

for i = 1:count
    a, b = get_data_sums[i]
    @test a == b
    @test a == put_data_sums[i]
end

if !haskey(ENV, "HTTP_JL_TEST_LONG_ASYNC")
    break
end

end # testset
end # if haskey(ENV, "AWS_ACCESS_KEY_ID")

configs = [
    [:verbose => 0],
    [:verbose => 0, :reuse_limit => 200],
    [:verbose => 0, :reuse_limit => 50]
]

@testset "async $count, $num, $config, $http" for count in 1:1,
                                            num in [100, 1000, 2000],
                                            config in configs,
                                            http in ["http", "https"]

println("running async $count, 1:$num, $config, $http A")

    result = []
    @sync begin
        for i = 1:min(num,100)
            @async try
                r = HTTP.request("GET",
                 "$http://httpbin.org/headers", ["i" => i]; config...)
                r = JSON.parse(String(r.body))
                push!(result, r["headers"]["I"] => string(i))
            catch e
                dump_async_exception(e, catch_stacktrace())
                rethrow(e)
            end
        end
    end
    for (a,b) in result
        @test a == b
    end

    HTTP.ConnectionPool.showpool(STDOUT)
    HTTP.ConnectionPool.closeall()

    result = []

println("running async $count, 1:$num, $config, $http B")

    @sync begin
        for i = 1:min(num,100)
            @async try
                r = HTTP.request("GET",
                     "$http://httpbin.org/stream/$i"; config...)
                r = String(r.body)
                r = split(strip(r), "\n")
                push!(result, length(r) => i)
            catch e
                dump_async_exception(e, catch_stacktrace())
                rethrow(e)
            end
        end
    end

    for (a,b) in result
        @test a == b
    end

    HTTP.ConnectionPool.showpool(STDOUT)
    HTTP.ConnectionPool.closeall()

    result = []

#=
    asyncmap(i->begin
        n = i % 20 + 1
        str = ""
        r = HTTP.open("GET", "$http://httpbin.org/stream/$n";
                      retries=5, config...) do s
            str = String(read(s))
        end
        l = split(strip(str), "\n")
        #println("GOT $i $n")

        push!(result, length(l) => n)

    end, 1:num, ntasks=20)

    for (a,b) in result
        @test a == b
    end

    result = []
=#

println("running async $count, 1:$num, $config, $http C")

    @sync begin
        for i = 1:num
            n = i % 20 + 1
            @async try
                r = nothing
                str = nothing
                url = "$http://httpbin.org/stream/$n"
                if rand(Bool)
                    if rand(Bool)
                        for attempt in 1:4
                            try
                                #println("GET $i $n BufferStream $attempt")
                                s = BufferStream()
                                r = HTTP.request(
                                    "GET", url; response_stream=s, config...)
                                @assert r.status == 200
                                str = String(read(s))
                                break
                            catch e
                                if attempt == 10 ||
                                   !HTTP.RetryRequest.isrecoverable(e)
                                    rethrow(e)
                                end
                                buf = IOBuffer()
                                println(buf, "$i retry $e $attempt...")
                                write(STDOUT, take!(buf))
                                sleep(0.1)
                            end
                        end
                    else
                        #println("GET $i $n Plain")
                        r = HTTP.request("GET", url; config...)
                        @assert r.status == 200
                        str = String(r.body)
                    end
                else
                    #println("GET $i $n open()")
                    r = HTTP.open("GET", url; config...) do http
                        str = String(read(http))
                    end
                    @assert r.status == 200
                end

                l = split(strip(str), "\n")
                #println("GOT $i $n $(length(l))")
                if length(l) != n
                    @show r
                    @show str
                end
                push!(result, length(l) => n)
            catch e
                push!(result, e => n)
                dump_async_exception(e, catch_stacktrace())
                rethrow(e)
            end
        end
    end

    for (a,b) in result
        @test a == b
    end

    HTTP.ConnectionPool.showpool(STDOUT)
    HTTP.ConnectionPool.closeall()

    if !haskey(ENV, "HTTP_JL_TEST_LONG_ASYNC")
        break
    end

end # testset

stop_pool_dump=true

HTTP.ConnectionPool.showpool(STDOUT)

println("async tests done")


text:
function DormandPrince853IntegratorBuilder(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return DormandPrince853IntegratorBuilder((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function build_integrator(obj::DormandPrince853IntegratorBuilder, arg0::Orbit, arg1::OrbitType)
    return jcall(obj, "buildIntegrator", AbstractIntegrator, (Orbit, OrbitType), arg0, arg1)
end


text:
push!(LOAD_PATH, "$(homedir())/Development/")
using JSort

path = "$(homedir())/master/sortering/"*"28si.yaml"
parameters = Parameters(path)
#sortfile(parameters)
#savecoefficients_eΔe(parameters, coefficients...)
#@time events = loadlabr(parameters);
#@time makeeΔebin(events, parameters)
@time coeff_e, coeff_Δe = featurealign2d(parameters.readpath, plot=true)
@show coeff_e
p_e = ParticleCalibrator(coeff_e, var=:e)
p_Δe = ParticleCalibrator(coeff_e, var=:Δe)
write(joinpath(parameters.savepath, "coefficients_e.txt"), p_e)
write(joinpath(parameters.savepath, "coefficients_de.txt"), p_Δe)
#@time savecoefficients_eΔe(parameters, coefficients...)

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
using DataFrames

"""
Returns the version of Plink associated with the command `plink`
"""
plink_version() = parse(Float64,match(r"PLINK\D+v([\d\.]+)",readall(`plink --help`))[1])

abstract PlinkOps
abstract PlinkInput <: PlinkOps

type PlinkBinaryFile <: PlinkInput
	resolved::Bool
	prev::PlinkOps
	next::Array{PlinkOps,1}
	bed::ASCIIString
	bim::ASCIIString
	fam::ASCIIString
	PlinkBinaryFile(bed,bim,fam) = (
		x = new();
		x.resolved=false;
		x.prev=x;
		x.next=[];
		x.bed=bed;
		x.bim=bim;
		x.fam=fam;
		x)

PlinkBinaryFile(p,bed,bim,fam) = (
		x = new();
		x.resolved=false;
		x.prev=p;
		x.next=[];
		x.bed=bed;
		x.bim=bim;
		x.fam=fam;
		x)
end

type PlinkFile <: PlinkInput
	resolved::Bool
	prev::PlinkOps
	next::Array{PlinkOps,1}
	ped::ASCIIString
	map::ASCIIString
	PlinkFile(ped,map) = (
		x = new();
		x.resolved=false;
		x.prev=x;
		x.next=[];
		x.ped=ped;
		x.map=map;
		x)
end

type PlinkVCF <: PlinkInput
	resolved::Bool
	prev::PlinkOps
	next::Array{PlinkOps,1}
	vcf::ASCIIString
end

type PlinkOption <: PlinkOps
	resolved::Bool
	prev::PlinkOps
	next::Array{PlinkOps,1}
	op::Symbol
	args::Array{ASCIIString,1}
	files::Dict{Symbol,ASCIIString}
end

# Binary (or greater) PlinkOption with output files
function PlinkOption(prev::PlinkOps,op::Symbol,args::Array{ASCIIString,1},files::Dict{Symbol,ASCIIString})
	opt=PlinkOption(false,prev,Array{PlinkOps,1}(),op,args,files)
	push!(opt.prev.next,opt)
	opt
end

# Binary (or greater) PlinkOption without output files
PlinkOption(prev::PlinkOps,op::Symbol,args::Array{ASCIIString,1}) = PlinkOption(prev,op,args,Dict{Symbol,ASCIIString}())

# Unary PlinkOption with no extra tokens or output files
PlinkOption(prev::PlinkOps,op::Symbol) = PlinkOption(prev,op,Array{ASCIIString,1}(),Dict{Symbol,ASCIIString}())

# Unary PlinkOption with file output with output files
PlinkOption(prev::PlinkOps,op::Symbol,files::Dict{Symbol,ASCIIString}) = PlinkOption(prev,op,Array{ASCIIString,1}(),files)

# Checks first for a binary with given plinkstem
# If not found, checks for a non-binary plink file
function plink(plinkstem::ASCIIString)
	bfile = map(x -> string(plinkstem,".",x), ["bed","bim","fam"])
	file = map(x -> string(plinkstem,".",x), ["ped","map"])
	vcf = string(plinkstem, ".vcf")

	if all(isfile,bfile)
		PlinkBinaryFile(bfile...)
	elseif all(isfile,file)
		PlinkFile(file...)
	elseif isfile(vcf)
		PlinkVCF(vcf)
	else
		error(string("No suitable plink files found for ", plinkstem))
	end
end

#plink(plinkstem::ASCIIString) = plink(ASCIIString(plinkstem))

sym_to_arg(s::Symbol) = string("--", replace(string(s),"_","-"))

function plink_arguments(p::PlinkInput)
	fields = setdiff(fieldnames(p),[:prev,:next,:resolved])
	opt = map(sym_to_arg,fields)
	arg = map(i->getfield(p,i), fields)
	paired = map(x->[x[1],x[2]],zip(opt,arg))
	foldl(append!,[],paired)
end

plink_arguments(p::PlinkOption) = vcat(sym_to_arg(p.op),p.args)

collect_until_input(p::PlinkInput) = [ p ]
collect_until_input(p::PlinkOps) = [ p ; collect_until_input(p.prev) ]

function split_operations(p::PlinkOps)
	# A recode / make_bed operation always result
	# in a PlinkInput object, but it is not part of
	# the operations, but rather the result
	ptr = (typeof(p) <: PlinkInput)	? p.prev : p

	commands = Array{PlinkOps,1}()

	while ptr != ptr.prev
		opseq = collect_until_input(ptr)
		push!(commands, ptr)
		ptr = last(opseq).prev
	end

	return commands
end

build_plink_commands(p::PlinkInput) = `plink $(plink_arguments(p))`
build_plink_commands(p::PlinkOps) = `$(build_plink_commands(p.prev))  $(plink_arguments(p))`

plinkcmds(p::PlinkOps) = map(x->build_plink_commands(x),split_operations(p))

#build_command(pieces::Array{ASCIIString,1}) = foldl((x,y)->`$x $y`,``,pieces)

function get_file_name(p::PlinkOps, ext::ASCIIString)
	# Basecase: If there is no specification of out-name, so just assume "plink"
	if (typeof(p) <: PlinkInput)
		string("plink",ext)
	elseif (typeof(p)==PlinkOption && p.op == :out)
		string(p.args[1],ext)
	else
		get_file_name(p.prev,ext)
	end
end

function out(p::PlinkOps,dirpath::ASCIIString)
	# FIXME: Dirpath could be checked to be valid
	# FIXME: Check that a previous out was not allready specified
	PlinkOption(p,:out,[dirpath])
end

# Create "anonymous" plink option
function opt(p::PlinkOps,opt::Cmd)
	m=match(r"`--(\S+)\s?(.*)`",string(opt))
	println(m)
	if m[1]==nothing
		error(string("Invalid unchecked plink option: ", opt))
	elseif m[2]==""
		PlinkOption(p,Symbol(m[1]))
	else
		PlinkOption(p,Symbol(m[1]),map(ASCIIString,split(reverse(chomp(reverse(m[2])))," ")))
	end
end

extensions_dict(p,names::Array{Symbol,1}) = [ x => ASCIIString(get_file_name(p,string(".",x))) for x in names ]

##################################################################################
## Input filtering
##################################################################################

macro file_filter_function(name,namesym)
	eval(:(($name)(p::PlinkOps,f::ASCIIString ; range=false) =
		(isfile(f) ?
			PlinkOption(p,$namesym,range?["range",f]:[f])
			:
			error(string("No such file '", f, "' for plink option"))
		)
	))
end

## ID lists:
"""
`keep(p::PlinkOption,file)` accepts a space/tab-delimited text `file` with family
IDs in the first column and within-family IDs in the second column,
and removes all unlisted samples from the current analysis.
`remove` does the same for all listed samples.
Similarly `keep(p::PlinkOption,listed::Array{ASCIIString,2})` accepts a 2d array
and autogenerates `file` for plink.
"""
@file_filter_function(keep,:keep)

"""
`remove` accepts a space/tab-delimited text file with family IDs in the first column and within-family IDs in the second column,
and removes all listed samples from the current analysis. `keep` does the same for all unlisted samples.
"""
@file_filter_function(exclude,:exclude)

"""
`extract` normally accepts a text file with a list of variant IDs (usually one per line, but it's okay for them to just be separated by spaces), and removes all unlisted variants from the current analysis.
"""
@file_filter_function(extract,:extract)

"""`from` excludes all variants on different chromosomes than the named variant, as well as those with smaller base-pair position values."""
from(p::PlinkOps,variant) = PlinkOption(p,:from, [ variant ])

"""
--from excludes all variants on different chromosomes than the named variant, as well as those with smaller base-pair position values. --to is similar, excluding variants with larger position values instead. If they are used together but the --from variant is after the --to variant, they are automatically swapped.
"""
to(p::PlinkOps,variant::ASCIIString) = PlinkOption(p,:to, [ variant ])

snp(p::PlinkOps,variant::ASCIIString) = PlinkOption(p,:snp, [ variant ])

window(p::PlinkOps,windows_size_kb::Int32) = PlinkOption(p,:snp, [ string(windows_size_kb) ])

exclude_snp(p::PlinkOps,variant::ASCIIString) = PlinkOption(p,:exclude_snp, [ variant ])

"""
geno filters out all variants with missing call rates exceeding the provided value (default 0.1) to be removed
"""
geno(p::PlinkOps,threshold::Float64) = PlinkOption(p,:geno, [ string(threshold) ])

"""
mind filters out all individuals with missing call rates exceeding the provided value (default 0.1) to be removed
"""
mind(p::PlinkOps,threshold::Float64) = PlinkOption(p,:geno, [ string(threshold) ])

maf(p::PlinkOps,threshold::Float64) = PlinkOption(p,:maf, [ string(threshold) ])
max_maf(p::PlinkOps,threshold::Float64) = PlinkOption(p,:max_maf, [ string(threshold) ])

mac(p::PlinkOps,threshold::Int32) = PlinkOption(p,:mac, [ string(threshold) ])
max_mac(p::PlinkOps,threshold::Int32) = PlinkOption(p,:max_mac, [ string(threshold) ])

"""filters out all variants which have Hardy-Weinberg equilibrium exact test p-value below the provided threshold."""
hwe(p::PlinkOps,threshold::Float64) = PlinkOption(p,:hwe, [ string(threshold) ])

""" filters out all samples with missing phenotypes """
prune(p::PlinkOps) = PlinkOption(p,:prune)

filter_cases(p::PlinkOps) = PlinkOption(p,:filter_cases)
filter_controls(p::PlinkOps) = PlinkOption(p,:filter_controls)
filter_males(p::PlinkOps) = PlinkOption(p,:filter_males)
filter_females(p::PlinkOps) = PlinkOption(p,:filter_females)
filter_founders(p::PlinkOps) = PlinkOption(p,:filter_founders)
filter_nonfounders(p::PlinkOps) = PlinkOption(p,:filter_nonfounders)

##################################################################################
## Basic statistics
##################################################################################

freq(p::PlinkOps) = PlinkOption(p,:freq, Dict(:frq => get_file_name(p,".frq")))

hardy(p::PlinkOps) = PlinkOption(p,:hardy, Dict(:hwe => get_file_name(p,".hwe")))

ibc(p::PlinkOps) = PlinkOption(p,:ibc, Dict(:hwe => get_file_name(p,".ibc")))

check_sex(p::PlinkOps) = PlinkOption(p,:check_sex, Dict(:hwe => get_file_name(p,".sexcheck")))

#missing(p::PlinkOps,threshold::Float64) = PlinkOption(p,:missing, [ string(threshold) ])

##################################################################################
## Data management
##################################################################################

function make_bed(p::PlinkOps)
	po = PlinkOption(p,:make_bed, extensions_dict(p,[:bed,:bim,:fam]))
	PlinkBinaryFile(po,po.files[:bed],po.files[:bim],po.files[:fam])
end

##################################################################################
## Reading and writing of plink file formats
##################################################################################

# Basic types for representing markers
# Compared to the binary plinkformat, this is a quite wasteful representation.
# So it will take up more memory and it does when stored on disk
# We may consider different options, e.g., runlength coding

type BiallelicVariantInfo
	chr::ASCIIString
	name::ASCIIString
	genetic_distance::Int32
	pos::Int32
	a1::ASCIIString # major allele
	a2::ASCIIString # minor allele

	# These are necessary to determine minor allele ordering
	a1_count::Int32
	a2_count::Int32
end

type BiallelicVariant
	# This is going to require a full-size pointer, and will
	# cause significant overhead. I wonder if it can be done
	# better using parametric types
	info::BiallelicVariantInfo
	genotype::Int8
end

function Base.show(io::IO,x::BiallelicVariant)
	if (x.genotype==0)
		print(io, string(x.info.a1,"/",x.info.a1))
	elseif (x.genotype==1)
		print(io, string(x.info.a1,"/",x.info.a2))
	elseif (x.genotype==2)
		print(io, string(x.info.a2,"/",x.info.a2))
	else
		"??"
	end
end

type Sex
		code::Int8
end

unknown_sex() = Sex(0)
male() = Sex(1)
female() = Sex(2)

function Base.show(io::IO,x::SexCode)
		if (x.code == 1)
				print(io, "male")
		elseif(x.code == 2)
				print(io, "female")
		else(x.code == 0)
			print(io,"unknown")
		end
end

each_trimmed_nonempty_line(io::IOStream) = filter(l->l!="", map(l->chomp(lstrip(l)), eachline(io)))
getfields(line) = filter(x->x!="",split(chomp(lstrip(line)),[' ','\t']))

function map_readfile(mapfile)
	println(readall(mapfile))
	markers_info = Array{BiallelicVariantInfo,1}()
	open(mapfile) do mapf
		for line in each_trimmed_nonempty_line(mapf)
			chr,snpname,cm,pos = getfields(line)
			push!(markers_info,BiallelicVariantInfo(chr,snpname,parse(Int32,cm),parse(Int32,pos),"","",0,0))
		end
	end
	return markers_info
end

function bim_readfile(bimfile)
	markers_info = Array{BiallelicVariantInfo,1}()
	open(bimfile) do bimf
		for line in each_trimmed_nonempty_line(bimf)
			chr,snpname,cm,pos,a1,a2 = getfields(line)
			push!(markers_info,BiallelicVariantInfo(chr,snpname,parse(Int32,cm),parse(Int32,pos),a1,a2,0,0))
		end
	end
	return markers_info
end

# With a dict we temporary relax constraints of equal number of colums per row
fam_dict() = Dict(
	:FID =>  data(Array{ASCIIString,1}()),
	:IID =>  data(Array{ASCIIString,1}()),
	:PaternalID =>  data(Array{ASCIIString,1}()),
	:MaternalID =>  data(Array{ASCIIString,1}()),
	:Sex => data(Array{Sex,1}()),
	:Phenotype => data(Array{Int8,1}()))

function fam_push!(fam::Dict, fields)
	push!(fam[:FID], shift!(fields))
	push!(fam[:IID], shift!(fields))
	push!(fam[:PaternalID], shift!(fields))
	push!(fam[:MaternalID], shift!(fields))
	push!(fam[:Sex], Sex(parse(Int8,shift!(fields))))
	push!(fam[:Phenotype], parse(Int8,shift!(fields)))
end

function fam_dataframe(dict::Dict)
	df = DataFrame(dict)
	for id in keys(dict)
		df[id] = map(i -> (i==-9 ? NA : i), df[id])
	end
	return df
end

function fam_readfile(famfile)
	fam=fam_dict()
	open(famfile) do f
		for line in each_trimmed_nonempty_line(f)
			fields = getfields(line)
			@assert 6==length(fields) "$famfile malformed on line $line $famfile"
			fam_push!(fam,fields)
		end
	end
	return fam_dataframe(fam)
end

function readplink(p::PlinkFile)
	markers_info = map_readfile(p.map)

	markers = Array{DataArray{BiallelicVariant,1},1}()
	fam=fam_dict()

	open(p.ped) do ped
		for line in each_trimmed_nonempty_line(ped)
			fields = getfields(line)

			fam_push!(fam,fields)

			for i in 1:length(markers_info)
				a1=shift!(fields)
				a2=shift!(fields)

				if (a1!="0")
					if markers_info[i].a1 == ""
						markers_info[i].a1 = a1
					end

					if markers_info[i].a2 == ""  && (a1 != a2 || a1!=markers_info[i].a1)
						if markers_info[i].a1 == a1
							markers_info[i].a2 = a2
						else
							markers_info[i].a2 = a1
						end
					end
				end

				genoint = a1=="0" ? -9 : a1!=a2 ? 1 : a1==markers_info[i].a1 ? 0 : 2
				genotype = BiallelicVariant(markers_info[i],genoint)

				if length(markers) < i
					push!(markers,DataArray{BiallelicVariant,1}[])
				end

				push!(markers[i],genotype)
				if genotype.genotype == -9
					markers[i][end] = NA
				end
			end
		end
	end

	df = fam_dataframe(fam)

	for i in 1:length(markers)
		df[Symbol(markers_info[i].name)] = markers[i]
	end

	return df
end

function byte_to_geno(byte)
	map(0:3) do pos
		(byte >> (2*pos)) & 3
	end
end

function readplink(p::PlinkBinaryFile)
	markers_info = bim_readfile(p.bim)
	fam = fam_readfile(p.fam)
	markers = Array{DataArray{BiallelicVariant,1},1}()

	println(markers_info)

	# The first three bytes should be 0x6c, 0x1b, and 0x01 in that order.
	bed=open(p.bed,"r")
	magic = readbytes(bed,)

	@assert magic == [0x6c, 0x1b, 0x01] "magic header in .bed file does not look like a binary plink file"

	# The rest of the file is a sequence of V blocks of N/4 (rounded up) bytes each
	# where V is the number of variants and N is the number of samples.
	# The first block corresponds to the first marker in the .bim file, etc.
	for mkr in markers_info
		last_byte_samples = nrow(fam) % 4
		ablock = readbytes(bed,Integer(ceil(nrow(fam)/4)))
		four_sample_blocks = map(byte_to_geno,ablock)
		if last_byte_samples != 0
			for i in 1:(4-last_byte_samples)
				pop!(last(four_sample_blocks))
			end
		end

		genotypes=reduce(vcat,[],four_sample_blocks)

		push!(markers,DataArray{BiallelicVariant,1}[])
		for geno in genotypes
			if (geno == 0)
				push!(markers[end],BiallelicVariant(mkr,0))
			elseif (geno == 1)
				push!(markers[end],NA)
			elseif (geno == 2)
				push!(markers[end],BiallelicVariant(mkr,1))
			elseif (geno == 3)
				push!(markers[end],BiallelicVariant(mkr,2))
			end
		end
	end
	close(bed)

	for i in 1:length(markers)
		fam[Symbol(markers_info[i].name)] = markers[i]
	end

	return fam
end

"""
Write "plink dataframe" to a binary plink file
"""
#function writeplink(plink::DataFrame)
#	fam_only = plink[:,[:FID,:IID,:PaternalID,:MaternalID,:Sex,:Phenotype]]
#	println(fam_only)
#
#	nrow(plink)
#
#fam_only
#
#	fam_only
#end

#dataframe(p::PlinkBinaryFile) = error("Not implemented yet")

dataframe(p::PlinkVCF) = error("Not implemented yet")

function dataframe(p::PlinkOps)
	for cmd in reverse(plinkcmds(p))
		run(cmd)
	end

	println(p.files)

	if (length(p.files) == 1)
		readtable(p.files[first(keys(p.files))],header=true,separator=' ')
	else
		error(string("Ambigous set of files-type: ", keys(p.files)))
	end
end

#function dataframe(p::PlinkOps,fileid::Symbol)
#	run(cmd(p))
#	if(p.files
#end


text:
module AmplNLWriter

import MathOptInterface
const MOI = MathOptInterface

"""
    AbstractSolverCommand

An abstract type that allows over-riding the call behavior of the solver.

See also: [`call_solver`](@ref).
"""
abstract type AbstractSolverCommand end

"""
    call_solver(
        solver::AbstractSolverCommand,
        nl_filename::String,
        options::Vector{String},
        stdin::IO,
        stdout::IO,
    )::String

Execute the `solver` given the NL file at `nl_filename`, a vector of `options`,
and `stdin` and `stdout`. Return the filename of the resulting `.sol` file.

You can assume `nl_filename` ends in `model.nl`, and that you can write a `.sol`
file to `replace(nl_filename, "model.nl" => "model.sol")`.

If anything goes wrong, throw a descriptive error.
"""
function call_solver end

struct _DefaultSolverCommand{F} <: AbstractSolverCommand
    f::F
end

function call_solver(
    solver::_DefaultSolverCommand,
    nl_filename::String,
    options::Vector{String},
    stdin::IO,
    stdout::IO,
)
    solver.f() do solver_path
        ret = run(
            pipeline(
                `$(solver_path) $(nl_filename) -AMPL $(options)`,
                stdin = stdin,
                stdout = stdout,
            ),
        )
        if ret.exitcode != 0
            error("Nonzero exit code: $(ret.exitcode)")
        end
    end
    return replace(nl_filename, "model.nl" => "model.sol")
end

"""
    _solver_command(x::Union{Function,String})

Functionify the solver command so it can be called as follows:
```julia
foo = _solver_command(x)
foo() do path
    run(`\$(path) args...`)
end
```
"""
_solver_command(x::String) = _DefaultSolverCommand(f -> f(x))
_solver_command(x::Function) = _DefaultSolverCommand(x)
_solver_command(x::AbstractSolverCommand) = x

struct _NLResults
    raw_status_string::String
    termination_status::MOI.TerminationStatusCode
    primal_status::MOI.ResultStatusCode
    objective_value::Float64
    primal_solution::Dict{MOI.VariableIndex,Float64}
    dual_solution::Vector{Float64}
    zL_out::Dict{MOI.VariableIndex,Float64}
    zU_out::Dict{MOI.VariableIndex,Float64}
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::MOI.FileFormats.NL.Model
    solver_command::AbstractSolverCommand
    options::Dict{String,Any}
    stdin::Any
    stdout::Any
    results::_NLResults
    solve_time::Float64
end

"""
    Optimizer(
        solver_command::Union{String,Function},
        solver_args::Vector{String};
        stdin::Any = stdin,
        stdout:Any = stdout,
    )

Create a new Optimizer object.

## Arguments

 * `solver_command`: one of two things:
   * A `String` of the full path of an AMPL-compatible executable
   * A function that takes takes a function as input, initializes any
     environment as needed, calls the input function with a path to the
     initialized executable, and then destructs the environment.
 * `solver_args`: a vector of `String` arguments passed solver executable.
   However, prefer passing `key=value` options via `MOI.RawOptimizerAttribute`.
 * `stdin` and `stdio`: arguments passed to `Base.pipeline` to redirect IO. See
   the Julia documentation for more details by typing `? pipeline` at the Julia
   REPL.

## Examples

A string to an executable:
```julia
Optimizer("/path/to/ipopt.exe")
```

A function or string provided by a package:
```julia
Optimizer(Ipopt.amplexe)
# or
Optimizer(Ipopt_jll.amplexe)
```

A custom function
```julia
function solver_command(f::Function)
    # Create environment ...
    ret = f("/path/to/ipopt")
    # Destruct environment ...
    return ret
end
Optimizer(solver_command)
```

The following two calls are equivalent:
```julia
# Okay:
model = Optimizer(Ipopt_jll.amplexe, ["print_level=0"])
# Better:
model = Optimizer(Ipopt_jll.amplexe)
MOI.set(model, MOI.RawOptimizerAttribute("print_level"), 0
```
"""
function Optimizer(
    solver_command::Union{AbstractSolverCommand,String,Function} = "",
    solver_args::Vector{String} = String[];
    stdin::Any = stdin,
    stdout::Any = stdout,
)
    return Optimizer(
        MOI.FileFormats.NL.Model(),
        _solver_command(solver_command),
        Dict{String,String}(opt => "" for opt in solver_args),
        stdin,
        stdout,
        _NLResults(
            "Optimize not called.",
            MOI.OPTIMIZE_NOT_CALLED,
            MOI.NO_SOLUTION,
            NaN,
            Dict{MOI.VariableIndex,Float64}(),
            Float64[],
            Dict{MOI.VariableIndex,Float64}(),
            Dict{MOI.VariableIndex,Float64}(),
        ),
        NaN,
    )
end

Base.show(io::IO, ::Optimizer) = print(io, "An AMPL (.nl) model")

function MOI.empty!(model::Optimizer)
    MOI.empty!(model.inner)
    model.results = _NLResults(
        "Optimize not called.",
        MOI.OPTIMIZE_NOT_CALLED,
        MOI.NO_SOLUTION,
        NaN,
        Dict{MOI.VariableIndex,Float64}(),
        Float64[],
        Dict{MOI.VariableIndex,Float64}(),
        Dict{MOI.VariableIndex,Float64}(),
    )
    model.solve_time = NaN
    return
end

MOI.is_empty(model::Optimizer) = MOI.is_empty(model.inner)

MOI.get(model::Optimizer, ::MOI.SolverName) = "AmplNLWriter"

MOI.supports(::Optimizer, ::MOI.Name) = true

MOI.get(model::Optimizer, ::MOI.Name) = MOI.get(model.inner, MOI.Name())

function MOI.set(model::Optimizer, ::MOI.Name, name::String)
    MOI.set(model.inner, MOI.Name(), name)
    return
end

MOI.supports(::Optimizer, ::MOI.RawOptimizerAttribute) = true

function MOI.get(model::Optimizer, attr::MOI.RawOptimizerAttribute)
    return get(model.options, attr.name, nothing)
end

function MOI.set(model::Optimizer, attr::MOI.RawOptimizerAttribute, value)
    model.options[attr.name] = value
    return
end

const _SCALAR_FUNCTIONS = Union{
    MOI.VariableIndex,
    MOI.ScalarAffineFunction{Float64},
    MOI.ScalarQuadraticFunction{Float64},
}

const _SCALAR_SETS = Union{
    MOI.LessThan{Float64},
    MOI.GreaterThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64},
}

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{<:_SCALAR_FUNCTIONS},
    ::Type{<:_SCALAR_SETS},
)
    return true
end

function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VariableIndex},
    ::Type{<:Union{MOI.ZeroOne,MOI.Integer}},
)
    return true
end

MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{<:_SCALAR_FUNCTIONS}) = true

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

MOI.supports(::Optimizer, ::MOI.NLPBlock) = true

MOI.supports_incremental_interface(::Optimizer) = false

MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) = MOI.copy_to(dest.inner, src)

function MOI.supports(
    ::Optimizer,
    ::MOI.VariablePrimalStart,
    ::Type{MOI.VariableIndex},
)
    return true
end

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart, x, v::Real)
    model.inner.x[x].start = Float64(v)
    return
end

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart, x, ::Nothing)
    model.inner.x[x].start = nothing
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.VariablePrimalStart,
    x::MOI.VariableIndex,
)
    return model.inner.x[x].start
end

"""
    _interpret_status(solve_result_num::Int, raw_status_string::String)

Convert the `solve_result_num` and `raw_status_string` into MOI-type statuses.

For the primal status, assume a solution is present. Other code is responsible
for returning `MOI.NO_SOLUTION` if no primal solution is present.
"""
function _interpret_status(solve_result_num::Int, raw_status_string::String)
    if 0 <= solve_result_num < 100
        # Solved, and nothing went wrong. Even though we say `LOCALLY_SOLVED`,
        # some solvers like SHOT use this status to represent problems that are
        # provably globally optimal.
        return MOI.LOCALLY_SOLVED, MOI.FEASIBLE_POINT
    elseif 100 <= solve_result_num < 200
        # Solved, but the solver can't be sure for some reason. e.g., SHOT
        # uses this for non-convex problems it isn't sure is the global optima.
        return MOI.LOCALLY_SOLVED, MOI.FEASIBLE_POINT
    elseif 200 <= solve_result_num < 300
        return MOI.INFEASIBLE, MOI.UNKNOWN_RESULT_STATUS
    elseif 300 <= solve_result_num < 400
        return MOI.DUAL_INFEASIBLE, MOI.UNKNOWN_RESULT_STATUS
    elseif 400 <= solve_result_num < 500
        return MOI.OTHER_LIMIT, MOI.UNKNOWN_RESULT_STATUS
    elseif 500 <= solve_result_num < 600
        return MOI.OTHER_ERROR, MOI.UNKNOWN_RESULT_STATUS
    end
    # If we didn't get a valid solve_result_num, try to get the status from the
    # solve_message string. Some solvers (e.g. SCIP) don't ever print the
    # suffixes so we need this.
    message = lowercase(raw_status_string)
    if occursin("optimal", message)
        return MOI.LOCALLY_SOLVED, MOI.FEASIBLE_POINT
    elseif occursin("infeasible", message)
        return MOI.INFEASIBLE, MOI.UNKNOWN_RESULT_STATUS
    elseif occursin("unbounded", message)
        return MOI.DUAL_INFEASIBLE, MOI.UNKNOWN_RESULT_STATUS
    elseif occursin("limit", message)
        return MOI.OTHER_LIMIT, MOI.UNKNOWN_RESULT_STATUS
    elseif occursin("error", message)
        return MOI.OTHER_ERROR, MOI.UNKNOWN_RESULT_STATUS
    else
        return MOI.OTHER_ERROR, MOI.UNKNOWN_RESULT_STATUS
    end
end

function _readline(io::IO)
    if eof(io)
        error("Reached end of sol file unexpectedly.")
    end
    return strip(readline(io))
end
_readline(io::IO, T) = parse(T, _readline(io))

function _read_sol(filename::String, model::MOI.FileFormats.NL.Model)
    return open(io -> _read_sol(io, model), filename, "r")
end

"""
    _read_sol(io::IO, model::MOI.FileFormats.NL.Model)

This function is based on a Julia translation of readsol.c, available at
https://github.com/ampl/asl/blob/64919f75fa7a438f4b41bce892dcbe2ae38343ee/src/solvers/readsol.c
and under the following license:

Copyright (C) 2017 AMPL Optimization, Inc.; written by David M. Gay.
Permission to use, copy, modify, and distribute this software and its
documentation for any purpose and without fee is hereby granted,
provided that the above copyright notice appear in all copies and that
both that the copyright notice and this permission notice and warranty
disclaimer appear in supporting documentation.

The author and AMPL Optimization, Inc. disclaim all warranties with
regard to this software, including all implied warranties of
merchantability and fitness.  In no event shall the author be liable
for any special, indirect or consequential damages or any damages
whatsoever resulting from loss of use, data or profits, whether in an
action of contract, negligence or other tortious action, arising out
of or in connection with the use or performance of this software.
"""
function _read_sol(io::IO, model::MOI.FileFormats.NL.Model)
    raw_status_string = ""
    line = ""
    while !startswith(line, "Options")
        line = _readline(io)
        raw_status_string *= line
    end
    # Read through all the options. Direct copy of reference implementation.
    @assert startswith(line, "Options")
    num_options = _readline(io, Int)
    need_vbtol = false
    if num_options > 0
        if !(3 <= num_options <= 9)
            error("expected num_options between 3 and 9; " * "got $num_options")
        end
        _readline(io, Int)  # Skip this line
        if _readline(io, Int) == 3
            num_options -= 2
            need_vbtol = true
        end
        for _ in 3:num_options
            _readline(io, Int)  # Skip the rest of the option lines
        end
    end
    # Read number of constraints
    num_cons = _readline(io, Int)
    @assert(num_cons == length(model.g) + length(model.h))
    # Read number of dual solutions to read in
    num_duals_to_read = _readline(io, Int)
    @assert(num_duals_to_read == 0 || num_duals_to_read == num_cons)
    # Read number of variables
    num_vars = _readline(io, Int)
    @assert(num_vars == length(model.x))
    # Read number of primal solutions to read in
    num_vars_to_read = _readline(io, Int)
    @assert(num_vars_to_read == 0 || num_vars_to_read == num_vars)
    # Skip over vbtol line if present
    if need_vbtol
        _readline(io)
    end
    # Read dual solutions
    dual_solution = Float64[_readline(io, Float64) for i in 1:num_duals_to_read]
    # Read primal solutions
    primal_solution = Dict{MOI.VariableIndex,Float64}()
    if num_vars_to_read > 0
        for xi in model.order
            primal_solution[xi] = _readline(io, Float64)
        end
    end
    # Check for status code
    solve_result_num = -1
    while !eof(io)
        linevals = split(_readline(io), " ")
        if length(linevals) > 0 && linevals[1] == "objno"
            @assert parse(Int, linevals[2]) == 0
            solve_result_num = parse(Int, linevals[3])
            break
        end
    end
    zL_out = Dict{MOI.VariableIndex,Float64}()
    zU_out = Dict{MOI.VariableIndex,Float64}()
    while !eof(io)
        line = _readline(io)
        if startswith(line, "suffix")
            items = split(line, " ")
            n_suffix = parse(Int, items[3])
            suffix = _readline(io)
            if !(suffix == "ipopt_zU_out" || suffix == "ipopt_zL_out")
                continue
            end
            for i in 1:n_suffix
                items = split(_readline(io), " ")
                x = model.order[parse(Int, items[1])+1]
                dual = parse(Float64, items[2])
                if suffix == "ipopt_zU_out"
                    zU_out[x] = dual
                else
                    @assert suffix == "ipopt_zL_out"
                    zL_out[x] = dual
                end
            end
        end
    end
    termination_status, primal_status =
        _interpret_status(solve_result_num, raw_status_string)
    objective_value = NaN
    if length(primal_solution) > 0
        # .sol files don't seem to be able to return the objective
        # value. Evaluate it here instead.
        objective_value = MOI.FileFormats.NL._evaluate(model.f, primal_solution)
    end
    return _NLResults(
        raw_status_string,
        termination_status,
        length(primal_solution) > 0 ? primal_status : MOI.NO_SOLUTION,
        objective_value,
        primal_solution,
        dual_solution,
        zL_out,
        zU_out,
    )
end

function MOI.optimize!(model::Optimizer)
    start_time = time()
    temp_dir = mktempdir()
    nl_file = joinpath(temp_dir, "model.nl")
    open(io -> write(io, model.inner), nl_file, "w")
    options = String[isempty(v) ? k : "$(k)=$(v)" for (k, v) in model.options]
    try
        sol_file = call_solver(
            model.solver_command,
            nl_file,
            options,
            model.stdin,
            model.stdout,
        )
        model.results = _read_sol(sol_file, model.inner)
    catch err
        model.results = _NLResults(
            "Error calling the solver. Failed with: $(err)",
            MOI.OTHER_ERROR,
            MOI.NO_SOLUTION,
            NaN,
            Dict{MOI.VariableIndex,Float64}(),
            Float64[],
            Dict{MOI.VariableIndex,Float64}(),
            Dict{MOI.VariableIndex,Float64}(),
        )
    end
    model.solve_time = time() - start_time
    return
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    return model.results.objective_value
end

function MOI.get(
    model::Optimizer,
    attr::MOI.VariablePrimal,
    x::MOI.VariableIndex,
)
    MOI.check_result_index_bounds(model, attr)
    return model.results.primal_solution[x]
end

MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    return model.results.termination_status
end

function MOI.get(model::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index != 1
        return MOI.NO_SOLUTION
    end
    return model.results.primal_status
end

function MOI.get(model::Optimizer, attr::MOI.DualStatus)
    n_duals =
        length(model.results.dual_solution) +
        length(model.results.zL_out) +
        length(model.results.zU_out)
    if attr.result_index != 1 ||
       n_duals == 0 ||
       model.results.termination_status != MOI.LOCALLY_SOLVED
        return MOI.NO_SOLUTION
    end
    return MOI.FEASIBLE_POINT
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    return model.results.raw_status_string
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT ? 1 : 0
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{<:MOI.VariableIndex},
)
    MOI.check_result_index_bounds(model, attr)
    return model.results.primal_solution[MOI.VariableIndex(ci.value)]
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{<:MOI.ScalarAffineFunction},
)
    MOI.check_result_index_bounds(model, attr)
    return MOI.FileFormats.NL._evaluate(
        model.inner.h[ci.value].expr,
        model.results.primal_solution,
    )
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintPrimal,
    ci::MOI.ConstraintIndex{<:MOI.ScalarQuadraticFunction},
)
    MOI.check_result_index_bounds(model, attr)
    return MOI.FileFormats.NL._evaluate(
        model.inner.g[ci.value].expr,
        model.results.primal_solution,
    )
end

function MOI.get(model::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    # TODO(odow): replace this with the proper dual objective.
    return MOI.get(model, MOI.ObjectiveValue())
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.LessThan{Float64}},
)
    MOI.check_result_index_bounds(model, attr)
    dual = get(model.results.zU_out, MOI.VariableIndex(ci.value), 0.0)
    return model.inner.sense == MOI.MIN_SENSE ? dual : -dual
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}},
)
    MOI.check_result_index_bounds(model, attr)
    dual = get(model.results.zL_out, MOI.VariableIndex(ci.value), 0.0)
    return model.inner.sense == MOI.MIN_SENSE ? dual : -dual
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}},
)
    MOI.check_result_index_bounds(model, attr)
    x = MOI.VariableIndex(ci.value)
    dual = get(model.results.zL_out, x, 0.0) + get(model.results.zU_out, x, 0.0)
    return model.inner.sense == MOI.MIN_SENSE ? dual : -dual
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,MOI.Interval{Float64}},
)
    MOI.check_result_index_bounds(model, attr)
    x = MOI.VariableIndex(ci.value)
    dual = get(model.results.zL_out, x, 0.0) + get(model.results.zU_out, x, 0.0)
    return model.inner.sense == MOI.MIN_SENSE ? dual : -dual
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{<:MOI.ScalarAffineFunction},
)
    MOI.check_result_index_bounds(model, attr)
    dual = model.results.dual_solution[length(model.inner.g)+ci.value]
    return model.inner.sense == MOI.MIN_SENSE ? dual : -dual
end

function MOI.get(
    model::Optimizer,
    attr::MOI.ConstraintDual,
    ci::MOI.ConstraintIndex{<:MOI.ScalarQuadraticFunction},
)
    MOI.check_result_index_bounds(model, attr)
    dual = model.results.dual_solution[ci.value]
    return model.inner.sense == MOI.MIN_SENSE ? dual : -dual
end

function MOI.get(model::Optimizer, attr::MOI.NLPBlockDual)
    MOI.check_result_index_bounds(model, attr)
    dual = model.results.dual_solution[1:model.inner.nlpblock_dim]
    return model.inner.sense == MOI.MIN_SENSE ? dual : -dual
end

end


text:
import Base.-, Base.+
# @inline -(a::WrappingInt32, b::WrappingInt32) = a.val - b.val
@inline +(a::WrappingInt32, b::UInt32) = WrappingInt32(a.val + b)
@inline -(a::WrappingInt32, b::UInt32) = a + -b
@inline +(a::WrappingInt32, b::Integer) = a + UInt32(b)
@inline -(a::WrappingInt32, b::Integer) = a + -UInt32(b)

function wrap(n::UInt, isn::WrappingInt32)
    tmp = UInt32((n << 32 >> 32))
    return isn + tmp
end

wrap(n::Integer, isn::WrappingInt32) = wrap(UInt(n), isn)

function unwrap(n::WrappingInt32, isn::WrappingInt32, checkpoint::Integer)
    ckpt = UInt(checkpoint)
    absolute_seqno_64 = UInt((n.val - isn.val))
    (ckpt <= absolute_seqno_64) && return absolute_seqno_64
    size_period = 0x100000000
    quotient = (ckpt - absolute_seqno_64) >> 32
    remainder = (ckpt - absolute_seqno_64) << 32 >> 32
    return absolute_seqno_64 + ((quotient + (remainder >= size_period ÷ 2)) << UInt32(32))
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
"""
    TaylorVector{N,T}

A data structure representing a vector of Taylor series with `N` terms.
Each element is an `NTuple{N,T}`.

    TaylorVector(T, N::Integer, n::Integer)

Create a vector of `n` `NTuple{N,T}`s.
"""
struct TaylorVector{N,T} <: AbstractVector{NTuple{N,T}}
    data::Matrix{T}
end
TaylorVector{N}(data::Matrix{T}) where {N,T} = TaylorVector{N,T}(data)
TaylorVector{N}(TV::TaylorVector{M,T}) where {N,M,T} = TaylorVector{N,T}(TV.data)
TaylorVector{N}(T, n::Integer) where {N} = TaylorVector{N}(zeros(T, N, n))

function Base.show(io::IO, ::MIME"text/plain", X::TaylorVector)
    summary(io, X)
    isempty(X) && return
    println(io, ":")
    Base.print_array(io, map(i -> X[i], 1:length(X)))
end

Base.length(TV::TaylorVector) = size(TV.data, 2)
Base.size(TV::TaylorVector) = (length(TV),)
Base.eltype(::Type{TaylorVector{N,T}}) where {N,T} = NTuple{N,T}
Base.IndexStyle(::TaylorVector) = Base.IndexLinear()

Base.fill!(v::TaylorVector, x) = fill!(v.data, x)

"""
    vectors(TV::TaylorVec{N})

Return the Taylor series as `N` seperate vectors.
"""
@generated function vectors(TV::TaylorVector{N}) where {N}
    Expr(:tuple, (:(view(TV.data, $k, :)) for k = 1:N)...)
end

@generated function Base.getindex(TV::TaylorVector{N}, i::Integer) where {N}
    quote
        Base.@_propagate_inbounds_meta
        x = TV.data
        $(Expr(:tuple, (:(x[$k, i]) for k = 1:N)...))
    end
end
Base.getindex(TV::TaylorVector, i::Integer, j::Integer) = getindex(TV.data, j, i)

@generated function Base.setindex!(
    TV::TaylorVector{N,T},
    x::NTuple{M,S},
    i::Integer,
) where {N,T,M,S}
    tup = Expr(:tuple, (i ≤ M ? :(convert(T, x[$i])) : :(zero(T)) for i = 1:N)...)
    :(setindex!(TV, $tup, i))
end
Base.setindex!(TV::TaylorVector, x, i::Integer) = setindex!(TV, make_ntuple(x), i)
function Base.setindex!(TV::TaylorVector{1}, x::Number, i::Integer)
    TV.data[i] = x
    x
end
@generated function Base.setindex!(
    TV::TaylorVector{N,T},
    x::NTuple{N,T},
    i::Integer,
) where {N,T}
    quote
        Base.@_propagate_inbounds_meta
        d = TV.data
        $(Expr(:tuple, (:(d[$k, i]) for k = 1:N)...)) = x
        x
    end
end
Base.setindex!(TV::TaylorVector, x, i::Integer, j::Integer) = setindex!(TV.data, x, j, i)


text:
# Test Julia's Ref mechanism and compare it to Ptr.
#

if !isdefined(:libn)
    const refc = Pkg.dir("CSDP", "test", "ref.c")
    const libn = Pkg.dir("CSDP", "deps", "usr", "lib", "ref.$(Libdl.dlext)")
    isdir(dirname(libn)) || mkdir(dirname(libn))
end
run(`gcc -fPIC -shared -o $libn -std=c99 $refc`)

const OP1 = UInt32(0)
const OP2 = UInt32(1)

"""Simple type containing an enum and Ptr"""
type S
    c::UInt32
    n::Cint
    e::Ptr{Cdouble}
end

"""Use Ref instead of Ptr"""
type R
    c::UInt32
    n::Cint
    e::Ref{Cdouble}
end

# Say hello, print sizeof(int)
println(ccall((:hello,libn), Cdouble, ()))

vec = Cdouble[1.0, 2.0]
println("&vec = $(pointer(vec))")

# Call with S --> should work
s = S(OP1, length(vec), pointer(vec))
ccall((:sum,libn), Cdouble, (S,), s)

# Call with R --> does not work
r = R(OP2, length(vec), Ref(vec))
ccall((:sum,libn), Cdouble, (R,), r)

# Now call with S, instead of R
Base.cconvert(::Type{R}, r::R) =
    S(r.c, r.n, Base.unsafe_convert(Ptr{Cdouble}, r.e))
ccall((:sum,libn), Cdouble, (S,), r)