text:
##############optimsamin
mutable struct IterBBState{TΘ⁺<:VolumePartsIter, T} <: AbstractIterState{T}
  Θ⁺::TΘ⁺
  λ⁺::T

  starttime::Float64
  lasttime::Float64
  besttime::Float64
  maxelapsedsincebest::Float64
  maxelapsed::Float64

  steps::Int
  beststeps::Int
  maxstepssincebest::Int
  maxsteps::Int

  #lastreporttime::Float64
  #minreporttime::Float64
  stop::Union{Nothing,Symbol}
end

IterBBState(Θ::TΘ, ::Type{T};
  Θ⁺::TΘ = Θ |> deepcopy,
  λ⁺::T = T(Inf),

  starttime::Float64=time(),
  lasttime::Float64=time(),
  besttime::Float64=time(),
  maxelapsedsincebest::Float64=PARAM[:bbmaxelapsedsincebest],
  maxelapsed::Float64=PARAM[:bbmaxelapsed],

  #NOTE: using steps as a limit only makes sense within a particular method
  steps::Int=0,
  beststeps::Int=0,
  maxstepssincebest::Int=PARAM[:bbmaxstepsincebest],
  maxsteps::Int=PARAM[:bbmaxsteps],

  #lastreporttime::Float64
  #minreporttime::Float64
  stop::Union{Nothing,Symbol} = nothing,
  ) where {TΘ,T} = IterBBState{TΘ,T}(Θ⁺, λ⁺,
    starttime, lasttime, besttime, maxelapsedsincebest, maxelapsed,
    steps, beststeps, maxstepssincebest, maxsteps,
    #lastreporttime, minreporttime,
    stop)

#updates the problem state given the solution for an iteration
function callbackbb!(bb, Θ::AbstractVolumeParts{<:Any, <:Any, elTΘ},
  rΨ::Ref{IterBBState{TΘ,elTΨ}}) where {elTΘ, TΘ, elTΨ}

  #get the current time
  t::Float64 = time()

  #this next code is useful for comparing methods
  if (bb.start_time > rΨ[].lasttime) &&  (rΨ[].steps>0)#reset the state if the op solver itself reset
    rΨ[] = IterBBState(Θ, elTΨ, starttime=bb.start_time)
  end

  #now update the state
  Ψ = rΨ[]
  elapsed::Float64 = t - Ψ.starttime
  Ψ.lasttime = t

  Ψ.steps::Int = BlackBoxOptim.num_steps(bb)
  stepsperssecond::Float64 = round(Ψ.steps/elapsed, digits=1)
  stepssincelast::Int = bb.num_steps_since_last_report
  Δnumbetter::Int = bb.num_better_since_last_report
  improvementpercent::Int = Δnumbetter ÷ stepssincelast*100

  Ψ.steps += stepssincelast

  λ::elTΨ = bb |> BlackBoxOptim.best_fitness |> elTΨ
  if λ < Ψ.λ⁺ #keep track of the best solution
    copyto!(Ψ.Θ⁺.G, bb |> best_candidate |> Vector{elTΘ})
    Ψ.λ⁺ = λ
    Ψ.besttime=t
    Ψ.beststeps=Ψ.steps
  end

  #use the solver for stopping criteria
  #run through the stopping criteria
  elapsedsincebest::Int = round(t - Ψ.besttime)
  stepssincebest::Int = Ψ.steps - Ψ.beststeps

  println("bestloss=$(Ψ.λ⁺)," *
    "time=$(round(elapsed))/$(Ψ.maxelapsed), " *
    "timesincebest=$elapsedsincebest/$(Ψ.maxelapsedsincebest), " *
    "steps=$(Ψ.steps), " *
    "improvements=$(Δnumbetter), " *
    "stepssincebest≈$stepssincebest, " *
    "steps/s=$stepsperssecond")
  bb.num_better_since_last_report = 0
  bb.num_steps_since_last_report = 0

  if !(Ψ.stop===nothing)
    throw("Previous shutdown call failed to cut the optimizer. Manual abort!")
  elseif elapsed > Ψ.maxelapsed
    @info "STOP: max time $(round(elapsed))/$(Ψ.maxelapsed) reached"
    Ψ.stop=:maxtime
  elseif elapsedsincebest>Ψ.maxelapsedsincebest
    @info "STOP: max time with minimal improvement $elapsedsincebest/$(Ψ.maxelapsedsincebest) reached"
    Ψ.stop=:maxtimesincebest
  elseif Ψ.steps ≥ Ψ.maxsteps
    @info "STOP: max steps $(Ψ.steps)/$(Ψ.maxsteps) reached." *
      "(NOTE- the meaning of this depends on the model)"
    Ψ.stop=:maxiter
  elseif stepssincebest  ≥ Ψ.maxstepssincebest
    @info "STOP: max steps since best $stepssincebest/$(Ψ.maxstepssincebest) reached. " *
      "(NOTE- the meaning of this depends on the model)"
    Ψ.stop=:maxstepssincebest
  end

  Ψ.stop===nothing || BlackBoxOptim.shutdown_optimizer!(bb)

  return nothing
end

function BBOptions(callback, Θ::VolumePartsIter;
    bbadditionaloptions = [:PopulationSize, :TraceMode])

  dims = (K=size(Θ.A,1), T=size(Θ.A,2))

  #construct the required parameters
  method = PARAM[:bbmethod]
  searchrange = (PARAM[:bblowerbound], PARAM[:bbupperbound])
  numdimensions=(dims.T-1)*dims.K
  callbackinterval=PARAM[:iterminreporttime]

  #make sure we only exit via the callback, setting others to be ignored
  maxfuncevals = 0
  maxsteps = 0 #some large number

  #construct the
  bboptionsdict = Dict(
    :Method=>method,
    :SearchRange=>searchrange,
    :NumDimensions=>numdimensions,
    :MaxFuncEvals=>maxfuncevals,
    :MaxSteps=>maxsteps,
    :CallbackFunction=>callback,
    :CallbackInterval=>callbackinterval) # we willl run the callback function each loop

  #load additional optional parameters
  loadedoptionalparams::Vector{Symbol} = Vector{Symbol}()
  for k ∈ bbadditionaloptions
    #parse the option name into its PARAM name
    paramname::Symbol = k |> string |> lowercase |> (s)->"bb$s" |> Symbol
    if !(PARAM[paramname] === nothing)
      bboptionsdict[k] = PARAM[paramname]
      push!(loadedoptionalparams, k)
    end
  end
  @info "Loaded bb optional params $loadedoptionalparams"

  return (; bboptionsdict...)
end

function solveA!(loss, Θ::TΘ,
    Xv::AbstractXYIter{TM, TV, T}, _Xv::XYIterAlloc{TM, TV, T},
    SolveType::Val{:bb};
    RegressionType::Val = Val(PARAM[:iterregressiontype]),
    GrowthType::Val = Val(PARAM[:itergrowthtype])) where {TΘ, TM, TV, T}

    dims = (K=size(Θ.A,1), T=size(Θ.A,2))

    rΨ::Ref{IterBBState{TΘ,Float64}} = Ref(IterBBState(Θ,Float64))
    callback(opt) = callbackbb!(opt,Θ,rΨ)

    bboptions = BBOptions(callback, Θ)

    x::Vector{T} = initializeΘgvector(Θ.g, T, GrowthType) #provides an initial vector for the optimizer

    #create an objective stub

    function obj(x::Vector)
      λ=updateA₁!(x, loss, Θ, Xv, _Xv, RegressionType, GrowthType) |> Float64
      return isfinite(λ) ? λ : Inf
    end

    #test the objective
    obj(x)
    local res
    try
      res = bboptimize(obj; bboptions...)
    catch err
      @warn "Aborted optimization with error $err. "
      Ψ.stop = :aborted
    finally

      t=time()
      Ψ=rΨ[]
      elapsed::Int = round(time() - Ψ.starttime)
      elapsedsincebest::Int = round(t - Ψ.besttime)
      stepssincebest::Int = Ψ.steps - Ψ.beststeps
      Θ.G .= Ψ.Θ⁺.G
      updateA₁!(loss, Θ, Xv, _Xv, RegressionType)
      updateAfromGAₜ!(Θ,Xv, _Xv,1) #update to the final values of A

      b = IOBuffer()
      write(b, "***Results info (BlackBoxOptim):\n" *
        "stop_reason: $(Ψ.stop)\n" *
        "iterations: $(Ψ.steps)/$(Ψ.maxsteps)\n" *
        "elapsed_at_completion: $(elapsed)/$(Ψ.maxelapsed)\n" *
        "elapsed since best: ~$(elapsedsincebest)/$(Ψ.maxelapsedsincebest)\n" *
        "iter since best: ~$(stepssincebest)/$(Ψ.maxstepssincebest)\n" *
        "Best resid: $(Ψ.λ⁺)\n")

      write(b, "\n\n***Other Simulation parameters:\n" *
        "growth type: $GrowthType\n" *
        "regression method: $(RegressionType)\n" *
        "report interval: $(PARAM[:iterminreporttime])\n" *
        "Types (Matrix, Vector, T): $TM, $TV, $T\n" *
        "bboptions:\n")
      for n ∈ propertynames(bboptions)
        occursin("callback", n |> string |> lowercase) && continue #callback options covered elsewhere
        write(b, "..$n=>$(bboptions[n])\n")
      end

      iterinfo = String(take!(b))
      println(iterinfo)
      open("$(PARAM[:outputpath])\\$(iterresultsname(Ψ.λ⁺))_info.txt", "w+") do f
        write(f, iterinfo)
      end

      return Ψ.λ⁺
    end

    @assert false #should never get here

  end

  #######BBOptim

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# Collection of sources required to build GMPBuilder
name = "GMP"
version = v"6.1.2"
sources = [
    "https://gmplib.org/download/gmp/gmp-6.1.2.tar.bz2" =>
    "5275bb04f4863a13516b2f39392ac5e272f5e1bb8057b18aec1c9b79d73d8fb2",

]

# Bash recipe for building across all platforms
script = raw"""
cd gmp-6.1.2
flags="--enable-cxx"
# On x86_64 architectures, build fat binary
if [[ ${proc_family} == intel ]]; then
    flags="${flags} --enable-fat"
fi
export CC="${CC} ${LDFLAGS}"
export CXX="${CXX} ${LDFLAGS}"
./configure --prefix=$prefix --host=$target --enable-shared --disable-static ${flags}
make -j
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libgmp", :libgmp)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl

max_stars_repo_name: isuruf/GMPBuilder

max_issues_repo_name: isuruf/GMPBuilder

max_forks_repo_name: isuruf/GMPBuilder


text:
login_payload= JSON.json(Dict("i"=>"", "t" => "login", "c"=> Dict("l"=>"admin", "p"=>"ep4ever", "t" => ""), "r" => "websocket"))

function _construct_request(i, command, item, c = "*", v = "", u = "")
    return JSON.json(
    Dict(
    "i"=>i,
    "t" => "request",
    "c" => [
    Dict("c"=> command,
        "p" => Dict( #parameters
            "p"=> Dict(
                "l"=> "0",  # line (0)
    			"a"=> "*",   #address (every device in line 0)
    			"c"=> c
                ),
            "i" => item,
            "v" => v,
            "u" => u
        )
    )
    ],
    "r" => "websocket")
    )
end

function _issue_command(ip, command::String, item::String, channel::String = "*", value::String = "", unit::String = "")
    results = []
    HTTP.WebSockets.open(ip) do ws
           write(ws, login_payload)
           d = JSON.parse(String(readavailable(ws)))
           session_id = d["i"]
           #@show session_id
           write(ws, _construct_request(session_id, command, item, channel, value, unit))
           d = JSON.parse(String(readavailable(ws)))
           push!(results, d)
           sleep(1)
    end
   results
end

export get_measured_HV
function get_measured_HV(ip = ip)
    results = []
    HTTP.WebSockets.open(ip) do ws
               write(ws, login_payload)
               d = JSON.parse(String(readavailable(ws)))
               session_id = d["i"]
               #@show session_id
               write(ws, _construct_request(session_id, "getItem", "Status.voltageMeasure"))
               d = JSON.parse(String(readavailable(ws)))
               push!(results, d)
           end;
    d = results[1]
    return parse(Float64,d[1]["c"][1]["d"]["v"]), parse(Float64,d[1]["c"][2]["d"]["v"])
end

export get_set_HV
function get_set_HV(ip = ip)
    results = []
    HTTP.WebSockets.open(ip) do ws
               write(ws, login_payload)
               d = JSON.parse(String(readavailable(ws)))
               session_id = d["i"]
               #@show session_id
               write(ws, _construct_request(session_id, "getItem", "Control.voltageSet"))
               d = JSON.parse(String(readavailable(ws)))
               push!(results, d)
           end;
    d = results[1]
    return parse(Float64,d[1]["c"][1]["d"]["v"]), parse(Float64,d[1]["c"][2]["d"]["v"])
end

export set_HV
function set_HV(ip::String, c::Int, v::Real)
    results = []
    HTTP.WebSockets.open(ip) do ws
               write(ws, login_payload)
               d = JSON.parse(String(readavailable(ws)))
               session_id = d["i"]
               #@show session_id
               write(ws, _construct_request(session_id, "setItem", "Control.voltageSet", "$c", "$v", "V"))
               d = JSON.parse(String(readavailable(ws)))
               push!(results, d)
           end;
    d = results[1][1]
    return d["trigger"]
end

export ramp_up
function ramp_up(ip::String, c::Int)
    results = []
    HTTP.WebSockets.open(ip) do ws
           write(ws, login_payload)
           d = JSON.parse(String(readavailable(ws)))
           session_id = d["i"]
           #@show session_id
           write(ws, _construct_request(session_id, "setItem", "Control.on", "$c", "1", ""))
           d = JSON.parse(String(readavailable(ws)))
           push!(results, d)
    end;
    d = results[1][1]
    return d["trigger"]
end

export ramp_down
function ramp_down(ip::String, c::Int)
    results = []
    HTTP.WebSockets.open(ip) do ws
           write(ws, login_payload)
           d = JSON.parse(String(readavailable(ws)))
           session_id = d["i"]
           #@show session_id
           write(ws, _construct_request(session_id, "setItem", "Control.on", "$c", "0", ""))
           d = JSON.parse(String(readavailable(ws)))
           push!(results, d)
    end;
    d = results[1][1]
    return d["trigger"]
end

export get_rampspeedUp
function get_rampspeedUp(ip = ip)
    results = _issue_command(ip, "getItem", "Control.voltageRampspeedUp")
    d = results[1]
    return parse(Float64,d[1]["c"][1]["d"]["v"]), parse(Float64,d[1]["c"][2]["d"]["v"])
end

export set_rampspeedUp
function set_rampspeedUp(ip, channel, value)
    results = _issue_command(ip, "setItem", "Control.voltageRampspeedUp", "$channel", "$value", "V/s")
    d = results[1][1]
    return d["trigger"]
end

export get_rampspeedDown
function get_rampspeedDown(ip = ip)
    results = _issue_command(ip, "getItem", "Control.voltageRampspeedDown")
    d = results[1]
    return parse(Float64,d[1]["c"][1]["d"]["v"]), parse(Float64,d[1]["c"][2]["d"]["v"])
end

export set_rampspeedDown
function set_rampspeedDown(ip, channel, value)
    results = _issue_command(ip, "setItem", "Control.voltageRampspeedDown", "$channel", "$value", "V/s")
    d = results[1][1]
    return d["trigger"]
end

export voltage_goto
function voltage_goto(ip::String, channel::Int, value::Real)
    results = []
    HTTP.WebSockets.open(ip) do ws
           write(ws, login_payload)
           d = JSON.parse(String(readavailable(ws)))
           session_id = d["i"]
           #@show session_id
           write(ws, _construct_request(session_id, "setItem", "Control.voltageSet", "$channel", "$value", "V"))
           d = JSON.parse(String(readavailable(ws)))
           push!(results, d)
           sleep(1)
           @info "Channel $channel: set Voltage: $value V\nstarting ramp ..."
           write(ws, _construct_request(session_id, "setItem", "Control.on", "$channel", "1", ""))
           d = JSON.parse(String(readavailable(ws)))
           push!(results, d)
    end;
    #d = results[2][1]
    #println(results)
    #return d["trigger"]
end


text:
using Optim
using Cubature
using Roots
using StatsFuns

# import Statistics: median, pdf, cdf
"pre-process pdf and cdf, given fixed pdf and cdf at x0, compute estimated support and check if pdf is proper"
function pre_process(x0::Array{T,2}, Fx0::Array{T,2}, pdf::Function, cdf::Function, dpdf::Function, quantbound::Function; yshift = 0.) where T<:Float64
    quantbound_fixed(p) = quantbound(x0, Fx0, p)
    pdf_fixed(y) = pdf(x0, Fx0, y)
    cdf_fixed(y) = cdf(x0, Fx0, y)
    dpdf_fixed(y) = dpdf(x0, Fx0, y)
    support = [1e-2, 5.]
    function support_comp!(pdf, support)
      current = pdf(support[1])
      for i in 1:5
        next = pdf(support[1]/10)
        if next < current
          support[1] /= 10
        else
          break
        end
        current = next
      end
      while pdf(support[2]) > 1e-6
        support[2] *= 1.2
      end
      support[2] += yshift
      # should make sure CDF(support[2]) - CDF(support[1]) > .96 to make 95% CI possible
      INT = cdf_fixed(support[2]) - cdf_fixed(support[1])
      if INT < .96
        @warn "pdf integral $INT"
      end
      return INT
    end
    INT = support_comp!(pdf_fixed, support)
    return (pdf_fixed, cdf_fixed, dpdf_fixed, quantbound_fixed, support, INT)
end

"wrap up all statistics computation"
function summary_comp(pdf_fixed::Function, cdf_fixed::Function, dpdf_fixed::Function, quantbound_fixed::Function, support::Array{T,1};
                       px = .5, confidence_level = .95) where T<:Float64
    quant_p, error_quant = quantile(cdf_fixed, quantbound, support; p=px)
    med, error_med = median(cdf_fixed, quantbound, support)
    mod = mode(pdf_fixed, support)
    CI_equal, error_CI_eq = credible_interval(cdf_fixed, quantbound, support;
                                                mode=:equal, wp = confidence_level)
    # CI_narrow, error_CI_nr = credible_interval(cdf_fixed, quantbound, support;
    #                                             mode=:narrow, wp = confidence_level)
    quantileInfo = (level = px, value = quant_p, error = error_quant)
    medianInfo = (value = med, error = error_med)
    CIequalInfo = (equal = CI_equal, error = error_CI_eq)
    # CInarrowInfo = (equal = CI_narrow, error = error_CI_nr)
    CInarrowInfo = nothing
    DistributionInfo = (quantile = quantileInfo, median = medianInfo, mode = mod, CIequal = CIequalInfo, CInarrow = CInarrowInfo)
    return DistributionInfo
end

"""
Given pdf, cdf and maybe pdf_deriv,
compute median, quantile, mode, symmetric/narrowest credible interval.
Warning: only for normalized values
"""
function median(cdf::Function, quantbound::Function, support::Array{T,1}; pdf = nothing, pdf_deriv=nothing) where T<:Float64
    med, err = quantile(cdf, quantbound, support)
    return med, err, bound
end

function quantile(cdf::Function, quantbound::Function, support::Array{T,1}; pdf = nothing, pdf_deriv=nothing, p::T=.5) where T<:Float64
    bound = support
    # try
    #   bound = quantbound(p)
    # catch err
    # end
    quant = fzero(y0 -> cdf(y0) - p, bound[1], bound[2])
    err = abs(p-cdf(quant))/p
    # status = err < 1e-5 ? 1 : 0
    return quant, err, bound
end

function mode(pdf::Function, support::Array{T,1}; cdf = nothing, pdf_deriv=nothing) where T<:Float64
    # maximize the pdf
    routine = optimize(x -> -pdf(x), support[1], support[2])
    mod = Optim.minimizer(routine)
    return mod
end

function credible_interval(cdf::Function, quantbound::Function, support::Array{T,1};
                            pdf=nothing, pdf_deriv=nothing, wp::T=.95, mode=:equal) where T<:Float64
    return credible_interval(cdf, quantbound, support, Val(mode); pdf_deriv=pdf_deriv, wp=wp)
end

function credible_interval(cdf::Function, quantbound::Function, support::Array{T,1}, ::Val{:equal};
                            pdf=nothing, pdf_deriv=nothing, wp::T=.95) where T<:Float64
    lower_qp = (1 - wp) / 2
    upper_qp = 1 - lower_qp
    lower_quant = lower_qp # random initialization
    try
        lower_quant = quantile(cdf, quantbound, support; p=lower_qp)[1]
    catch err
        wp = .9
        lower_qp = (1 - wp) / 2
        upper_qp = 1 - lower_qp
        lower_quant = quantile(cdf, quantbound, support; p=lower_qp)[1]
        @warn "Can't compute 95% credible intervel, do $(Int(wp*100))% instead."
    end
    upper_quant = quantile(cdf, quantbound, support; p=upper_qp)[1]
    err = abs(cdf(upper_quant) -  cdf(lower_quant) - wp)/wp
    return [lower_quant, upper_quant], wp, err
end

function credible_interval(cdf::Function, quantbound::Function, support::Array{T,1}, ::Val{:narrow};
                            pdf=nothing, pdf_deriv=nothing, wp::T=.95) where T<:Float64
  #=
  Brief idea: bisection
    Suppose the target interval is [alpha*, beta*], i.e. integral of pdf on [alpha*, beta*] = wp
    and the corresponding height is pdf(alpha*) = pdf(beta*) = l_height*
    Say l* is the horizontal line intersecting with pdf curve at alpha* and beta*,
    we first find two horizontal l lines sitting on each side of l*
    One is lower with height l_height_low, and one is higher with height l_height_high
    Then we use a bisection way to gradually find the target l*
  =#

  #= Notations:
  [alpha, beta]:
      interval for integration, target credible interval
      s.t. 1) pdf(alpha) = pdf(beta) and 2) integral of pdf from alpha to beta = wp (like 0.95)
  bound:
      the bound of support of the pdf, currently assume this is available
      i.e. pdf(x) = 0 if z > bound
  l_height_low/l_height_high:
      height of the l line that is lower/higher than the target one
  alpha_low/alpha_high, beta_low/beta_high:
      the corresponding interval
  int_low/int_high:
      the corresponding integral value, i.e. integral of pdf on [alpha_low/high, beta_low/high]
  hquadrature:
      function from the package Cubature for numerical integration
      e.g. integral_heightue, error = hquadrature(f, a, b) gives integral of f on [a,b] with some error
      we mostly use this to compute integral of pdf on [alpha, beta] for different alpha and beta
  =#

 # helper functions
  #=
  given a height l_height, find the two intersections s.t. pdf(α) = pdf(β) = l_height
  Input:
    l_height: given height
    α_intvl/β_intvl: interval for root finding for alpha/beta
  Output:
    [α, β]: the two intersections
    int: integral value of pdf on [α, β]
  =#
  function find_αβ(l_height, α_intvl, β_intvl)
    # find α and β within ginve intervals α_intvl and β_intvl respectively
    temp_fun(x) = pdf(x) - l_height
    α = fzero(temp_fun,  α_intvl[1], α_intvl[2])
    β = fzero(temp_fun,  β_intvl[1], β_intvl[2])
    int = hquadrature(x -> pdf(x), α, β)[1]
    return α, β, int
  end

  #=
  Adjust the height in case initial choice of low/high is not proper
    Since we want int_low > wp and int_high < wp to do bisection,
    we have to adjust height if initial int_low < wp or int_high > wp
  Input:
    l: current height
    [α, β]: current interval
    int: current integral value
    MODE: indicates we are adjusting the low line or the high line
  Output:
    According values after adjustment
    s.t. int > wp if MODE == "low"; int < wp else.
 =#
  function adjust(l, α, β, int, MODE)
    # adjust if int < wp
    if MODE == "low"
        while int < wp
            l /= 2
            α, β, int = find_αβ(l, [support[1], α], [β, support[2]])
        end
    else
        while int > wp
            l = (l + pdf(mode_d))/2
            α, β, int = find_αβ(l, [α, mode_d], [mode_d, β])
        end
    end
    return l, α, β, int
  end
  mode_d = mode(pdf, cdf, support)
  l_height_low = pdf(support[1])
  α_low = support[1]
  β_low = fzero(x -> pdf(x) - l_height_low,  mode_d, support[2])
  int_low = hquadrature(x -> pdf(x), α_low, β_low)[1]
  # α_low, β_low, int_low = find_αβ(l_height_low, [support[1], mode_d], [mode_d, support[2]])
  # adjust height if l low is not lower than l* (i.e. int < wp)
  # l_height_low, α_low, β_low, int_low = adjust(l_height_low, α_low, β_low, int_low, "low")
  l_height_high = pdf(mode_d)*0.9
  α_high, β_high, int_high = find_αβ(l_height_high, [α_low, mode_d], [mode_d, β_low])
  # adjust height if l_high is not higher than l* (i.e. int > wp)
  l_height_high, α_high, β_high, int_high = adjust(l_height_high, α_high, β_high, int_high, "high")

  α_mid = 0.
  β_mid = 0.
  l_mid = 0.
  N = 0
  int = int_low
  while !isapprox(int, wp) && N < 50
  l_mid = (l_height_high + l_height_low)/2
  α_mid, β_mid, int = find_αβ(l_mid, [α_low, α_high], [β_high, β_low])
  if int > wp
      l_height_low = l_mid
      α_low = α_mid
      β_low = β_mid
  else
      l_height_high = l_mid
      α_high = α_mid
      β_high = β_mid
  end
  N += 1
  end
  err1 = abs(cdf(β_mid) -  cdf(α_mid) - wp)/wp
  err = max(err1, abs(pdf(α_mid)- pdf(β_mid))/abs(pdf(β_mid)))
  return ([α_mid, β_mid], err)

end

# @doc raw"""
# """
# function map_estimate(btg::BTG)
#     # TODO
# end

# @doc raw"""
# """
# function cross_validate(btg::BTG)
#     # TODO
# end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using NumCME
using Catalyst
using StaticArrays: @MVector
using Sundials: CVODE_BDF

const r1 = 6.1e-3
const r2 = 6.9e-3
const η = 5.9
const Ahog = 9.3e9
const Mhog = 2.2e-2
function Hog1p(t)
    u = (1.0 - exp(-r1 * t)) * exp(-r2 * t)
    signal = Ahog * (u / (1.0 + u / Mhog))^η
    return signal
end
@parameters k01, k10, a, k12, k21, k23, k32, λ0, λ1, λ2, λ3, ktrans, γnuc, γcyt
rn = @reaction_network begin
    k01, G0 --> G1
    max(0, k10 - a*Hog1p(t)), G1 --> G0
    k12, G1 --> G2
    k21, G2 --> G1
    k23, G2 --> G3
    k32, G3 --> G2
    λ0, G0 --> G0 + RNAnuc
    λ1, G1 --> G1 + RNAnuc
    λ2, G2 --> G2 + RNAnuc
    λ3, G3 --> G3 + RNAnuc
    γnuc, RNAnuc --> ∅
    ktrans, RNAnuc --> RNAcyt
    γcyt, RNAcyt --> ∅
end k01 k10 a k12 k21 k23 k32 λ0 λ1 λ2 λ3 γnuc ktrans γcyt

param_values = Dict([
k01=> 2.6e-3,
k10=> 1.9e01,
a=> 0.0,
k12=> 7.63e-3,
k21=> 1.2e-2,
k23=> 4e-3,
k32=> 3.1e-3,
λ0=> 5.9e-4,
λ1=> 1.7e-1,
λ2=> 1.0,
λ3=> 3e-2,
ktrans=> 2.6e-1,
γnuc=> 2.2e-6,
γcyt=> 8.3e-3
])
# Simulate long-time behavior before MAPK signal
model = CmeModel(rn, collect(param_values))
p0 = FspVectorSparse([@MVector [1,0,0,0,0,0]], [1.0])
fsp = AdaptiveFspSparse(
    ode_method = CVODE_BDF(linear_solver=:GMRES),
    space_adapter = RStepAdapter(10,20,true)
)
sol = solve(
    model,
    p0,
    (0.0, 8*3600.0),
    fsp;
    saveat = [8*3600.0],
    fsptol = 1.0e-6,
    odeatol = 1.0e-14,
    odertol = 1.0e-6,
    verbose=false
);
pend = sol[end].p
# Simulate behavior after MAPK addition
param_values[a] = 3.2e04
model = CmeModel(rn, collect(param_values))
sol = solve(
    model,
    pend,
    (0.0, 3600.0),
    fsp;
    saveat = [t for t in 0.0:60.0:3600.0],
    fsptol = 1.0e-4,
    odeatol = 1.0e-14,
    odertol = 1.0e-6,
    verbose=true
);

#
using Plots

anim = @animate for i ∈ 1:length(sol)
    pfull = sol[i].p
    pnuc = sum(pfull, [1,2,3,4,6]) |> Array
    pcyt = sum(pfull, [1,2,3,4,5]) |> Array
    rnaplot = plot()
    plot!(rnaplot, 0:length(pnuc)-1, pnuc, color=:blue, label="Nuclear mRNA")
    plot!(rnaplot, 0:length(pcyt)-1, pcyt, color=:red, label="Cytoplasmic mRNA")
    ylims!(rnaplot, (0.0, 1.0))
    xlims!(rnaplot, (0, 20))
    xlabel!(rnaplot, "Molecule count")
    ylabel!(rnaplot,"Probability")
    title!(rnaplot, "t = $(sol.t[i]/60.0) minute")
end
gif(anim, "hog1p.gif", fps = 10)


text:
export ArcEagerSystem
export execute_transition, is_transition_valid, transitions_number, set_labels!

mutable struct ArcEagerSystem <: ParsingSystem
  labels::Vector{String}
  transitions::Vector{Transition}

  ArcEagerSystem() = new()
  function ArcEagerSystem(labels::Vector{String})
    system = new()
    set_labels!(system, labels)

    system
  end
end

function set_labels!(system::ArcEagerSystem, labels::Vector{String})
  transitions = Vector{Transition}()

  push!(transitions, Transition(Shift(), EMPTY_LABEL))
  push!(transitions, Transition(Reduce(), EMPTY_LABEL))

  foreach(labels) do label
    for move in [LeftArc(), RightArc()]
      push!(transitions, Transition(move, label))
    end
  end

  system.labels = labels
  system.transitions = transitions
end

function transitions_number(system::ArcEagerSystem)
  length(system.transitions)
end

function execute_transition(config::Configuration, trans::Transition, system::ArcEagerSystem)
  if !is_transition_valid(config, trans, system)
    return false
  end

  transition(config, trans.label, trans.move)

  true
end

function is_transition_valid(config::Configuration, transition::Transition, system::ArcEagerSystem)
  is_valid(config, transition.move)
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
struct SOneTo{n} <: AbstractUnitRange{Int} end

@pure SOneTo(n::Int) = SOneTo{n}()
function SOneTo{n}(r::AbstractUnitRange) where n
    ((first(r) == 1) & (last(r) == n)) && return SOneTo{n}()

    errmsg(r) = throw(DimensionMismatch("$r is inconsistent with SOneTo{$n}")) # avoid GC frame
    errmsg(r)
end
Base.Tuple(::SOneTo{N}) where N = ntuple(identity, Val(N))

@pure Base.axes(s::SOneTo) = (s,)
@pure Base.size(s::SOneTo{n}) where n = (n,)
@pure Base.length(s::SOneTo{n}) where n = n

# The axes of a Slice'd SOneTo use the SOneTo itself
Base.axes(S::Base.Slice{<:SOneTo}) = (S.indices,)
Base.unsafe_indices(S::Base.Slice{<:SOneTo}) = (S.indices,)
Base.axes1(S::Base.Slice{<:SOneTo}) = S.indices

@propagate_inbounds function Base.getindex(s::SOneTo, i::Int)
    @boundscheck checkbounds(s, i)
    return i
end
@propagate_inbounds function Base.getindex(s::SOneTo, s2::SOneTo)
    @boundscheck checkbounds(s, s2)
    return s2
end

@pure Base.first(::SOneTo) = 1
@pure Base.last(::SOneTo{n}) where n = n::Int
@pure Base.iterate(::SOneTo{n}) where n = n::Int < 1 ? nothing : (1, 1)
@pure function Base.iterate(::SOneTo{n}, s::Int) where {n}
    if s < n::Int
        s2 = s + 1
        return (s2, s2)
    else
        return nothing
    end
end

function Base.getproperty(::SOneTo{n}, s::Symbol) where {n}
    if s === :start
        return 1
    elseif s === :stop
        return n::Int
    else
        error("type SOneTo has no property $s")
    end
end

Base.show(io::IO, ::SOneTo{n}) where {n} = print(io, "SOneTo(", n::Int, ")")
Base.@pure function Base.checkindex(::Type{Bool}, ::SOneTo{n1}, ::SOneTo{n2}) where {n1, n2}
    return n1::Int >= n2::Int
end

Base.promote_rule(a::Type{Base.OneTo{T}}, ::Type{SOneTo{n}}) where {T,n} =
    Base.OneTo{promote_type(T, Int)}


text:
# canonical data version: 1.1.0

using Test

include("minesweeper.jl")

@testset "no rows" begin
    @test annotate([]) == []
end

@testset "no columns" begin
     @test annotate([""]) == [""]
end

@testset "no mines" begin
    minefield = ["   ",
                 "   ",
                 "   "]
    out = ["   ",
           "   ",
           "   "]
    @test annotate(minefield) == out
end

@testset "minefield only mines" begin
    minefield = ["***",
                 "***",
                 "***"]
    out = ["***",
           "***",
           "***"]
    @test annotate(minefield) == out
end

@testset "mine surrounded by spaces" begin
    minefield = ["   ",
                 " * ",
                 "   "]
    out = ["111",
           "1*1",
           "111"]
    @test annotate(minefield) == out
end

@testset "space surrounded by mines" begin
    minefield = ["***",
                 "* *",
                 "***"]
    out = ["***",
           "*8*",
           "***"]
    @test annotate(minefield) == out
end

@testset "horizontal line" begin
    minefield = [" * * "]
    out = ["1*2*1"]
    @test annotate(minefield) == out
end

@testset "horizontal line mines at edges" begin
    minefield = ["*   *"]
    out = ["*1 1*"]
    @test annotate(minefield) == out
end

@testset "vertical line" begin
    minefield = [" ",
                 "*",
                 " ",
                 "*",
                 " "]
    out = ["1",
           "*",
           "2",
           "*",
           "1"]
    @test annotate(minefield) == out
end

@testset "vertical line mines at edges" begin
    minefield = ["*",
                 " ",
                 " ",
                 " ",
                 "*"]
    out = ["*",
           "1",
           " ",
           "1",
           "*"]
    @test annotate(minefield) == out
end

@testset "cross" begin
    minefield = ["  *  ",
                 "  *  ",
                 "*****",
                 "  *  ",
                 "  *  "]
    out = [" 2*2 ",
           "25*52",
           "*****",
           "25*52",
           " 2*2 "]
    @test annotate(minefield) == out
end

@testset "large minefield" begin
    minefield = [" *  * ",
                 "  *   ",
                 "    * ",
                 "   * *",
                 " *  * ",
                 "      "]
    out = ["1*22*1",
           "12*322",
           " 123*2",
           "112*4*",
           "1*22*2",
           "111111"]
    @test annotate(minefield) == out
end


text:
export marginalrule

@marginalrule KernelGCV(:y_x) (m_y::MvNormalMeanCovariance, m_x::MvNormalMeanCovariance, q_z::MvNormalMeanCovariance, meta::KernelGCVMetadata) = begin
    kernelfunction = get_kernelfn(meta)
    Λ = approximate_kernel_expectation(get_approximation(meta), (z) -> cholinv(kernelfunction(z)), q_z)

    Λy = invcov(m_y)
    Λx = invcov(m_x)

    wy = Λy * mean(m_y)
    wx = Λx * mean(m_x)

    C = cholinv([ Λ + Λy -Λ; -Λ Λ + Λx ])
    m = C * [ wy ; wx ]

    return MvNormalMeanCovariance(m, C)
end

@marginalrule KernelGCV(:y_x) (m_y::MvNormalMeanPrecision, m_x::MvNormalMeanPrecision, q_z::MvNormalMeanPrecision, meta::KernelGCVMetadata) = begin
    kernelfunction = get_kernelfn(meta)
    C = approximate_kernel_expectation(get_approximation(meta), (z) -> cholinv(kernelfunction(z)), q_z)

    Cy = invcov(m_y)
    Cx = invcov(m_x)

    wy = Cy * mean(m_y)
    wx = Cx * mean(m_x)

    Λ = [ C + Cy -C; -C C + Cx ]
    μ = cholinv(Λ) * [ wy ; wx ]

    return MvNormalMeanPrecision(μ, Λ)
end


text:
function get_param_values(arg)
    key = ENV["USDA_QUICK_SURVEY_KEY"]

    request("GET", string(usda_url, "/api/get_param_values/?key=", key, "&param=", arg))
end


text:
"""
    eom_logistic4(u, p, t)

Equations of motion for a 4D transitive causal chain of unidirectionally
coupled logistic maps, where ``y_1 \\to y_2 \\to y_3 \\to y_4`` [1].

*Note: With the default parameters which are as in [1], for some initial conditions,
this system wanders off to ``\\pm \\infty`` for some of the variables. Make sure that
you have a good realisation before using the orbit for anything.**

## Equations of motion

```math
\\begin{aligned}
y_1(t+1) &= y_1(t)(r_1 - r_1 y_1) \\\\
y_2(t+1) &= y_2(t)(r_2 - c_{12} y_1 - r_2 y_2) \\\\
y_3(t+1) &= y_3(t)(r_3 - c_{23} y_2 - r_3 y_3) \\\\
y_4(t+1) &= y_4(t)(r_4 - c_{34} y_3 - r_4 y_4)
\\end{aligned}
```

## References

1. Ye, Hao, et al. "Distinguishing time-delayed causal interactions using
    convergent cross mapping." Scientific reports 5 (2015): 14750
"""
function eom_logistic4(u, p, t)
    r₁, r₂, r₃, r₄, c₁₂, c₂₃, c₃₄  = (p...,)
    y₁, y₂, y₃, y₄ = (u...,)

    dy₁ = y₁*(r₁ - r₁*y₁)
    dy₂ = y₂*(r₂ - c₁₂*y₁ - r₂*y₂)
    dy₃ = y₃*(r₃ - c₂₃*y₂ - r₃*y₃)
    dy₄ = y₄*(r₄ - c₃₄*y₃ - r₄*y₄)
    return SVector{4}(dy₁, dy₂, dy₃, dy₄)
end

function logistic4(u₀, r₁, r₂, r₃, r₄, c₁₂, c₂₃, c₃₄)
    p = [r₁, r₂, r₃, r₄, c₁₂, c₂₃, c₃₄]
    DiscreteDynamicalSystem(eom_logistic4, u₀, p)
end

"""
    logistic4(;u₀ = rand(4), r₁ = 3.9, r₂ = 3.6, r₃ = 3.6, r₄ = 3.8,
        c₁₂ = 0.4, c₂₃ = 0.4, c₃₄ = 0.35) -> DiscreteDynamicalSystem

Initialise a system of a transitive chain of four unidirectionally coupled
logistic maps, where ``y_1 \\to y_2 \\to y_3 \\to y_4`` [1]. Default
parameters are as in [1].

*Note: With the default parameters which are as in [1], for some initial conditions,
this system wanders off to ``\\pm \\infty`` for some of the variables. Make sure that
you have a good realisation before using the orbit for anything.**

## Equations of motion

```math
\\begin{aligned}
y_1(t+1) &= y_1(t)(r_1 - r_1 y_1) \\\\
y_2(t+1) &= y_2(t)(r_2 - c_{12} y_1 - r_2 y_2) \\\\
y_3(t+1) &= y_3(t)(r_3 - c_{23} y_2 - r_3 y_3) \\\\
y_4(t+1) &= y_4(t)(r_4 - c_{34} y_3 - r_4 y_4)
\\end{aligned}
```

## References

1. Ye, Hao, et al. "Distinguishing time-delayed causal interactions using
    convergent cross mapping." Scientific reports 5 (2015): 14750
"""
logistic4(;u₀ = rand(4),
            r₁ = 3.9, r₂ = 3.6, r₃ = 3.6, r₄ = 3.8,
            c₁₂ = 0.4, c₂₃ = 0.4, c₃₄ = 0.35) =
    logistic4(u₀, r₁, r₂, r₃, r₄, c₁₂, c₂₃, c₃₄)


text:
function nlp_cvx_206_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - linear objective
    # - intersection convex quadratic constraints
    # - power cones
    # Variants
    #   010 - intersection set

    model = Model(optimizer)

    # NOTE, starting any of these at 0.0 will segfault libcoinmumps
    @variable(model, x >= 0, start = 0.1)
    @variable(model, y >= 0, start = 0.1)
    @variable(model, z >= 0, start = 0.1)

    @objective(model, Max, 2 * x + y + z)
    @NLconstraint(model, z <= x^(0.3) * y^(0.7))
    @NLconstraint(model, x <= z^(0.7) * y^(0.3))
    @NLconstraint(model, x^2 + y^2 <= z + 1)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, 4.0, tol = objective_tol)
    return check_solution([x, y, z], [1, 1, 1], tol = primal_tol)
end


text:
#=================================
Methods using Interval Arithmetic
=================================#

function _enclose(::NaturalEnclosure, f::Function, dom::Interval_or_IntervalBox; kwargs...)
    return f(dom...)
end


text:
"""
    PyIO(x; own=false, text=missing, buflen=4096)

Wrap the Python IO stream `x` as a Julia IO stream.

When this goes out of scope and is finalized, it is automatically flushed. If `own=true` then it is also closed.

If `text=false` then `x` must be a binary stream and arbitrary binary I/O is possible.
If `text=true` then `x` must be a text stream and only UTF-8 must be written (i.e. use `print` not `write`).
If `text` is not specified then it is chosen automatically.
If `x` is a text stream and you really need a binary stream, then often `PyIO(x.buffer)` will work.

For efficiency, reads and writes are buffered before being sent to `x`.
The size of the buffers is `buflen`.
The buffers are cleared using `flush`.
"""
mutable struct PyIO <: IO
    py::Py
    # true to close the file automatically
    own::Bool
    # true if `o` is text, false if binary
    text::Bool
    # true if we are definitely at the end of the file; false if we are not or don't know
    eof::Bool
    # input buffer
    ibuflen::Int
    ibuf::Vector{UInt8}
    # output buffer
    obuflen::Int
    obuf::Vector{UInt8}

    PyIO(::Val{:new}, py::Py, own::Bool, text::Bool, ibuflen::Int, obuflen::Int) = begin
        io = new(py, own, text, false, ibuflen, UInt8[], obuflen, UInt8[])
        finalizer(pyio_finalize!, io)
    end
end
export PyIO

function PyIO(o; own::Bool = false, text::Union{Missing,Bool} = missing, buflen::Integer = 4096, ibuflen::Integer = buflen, obuflen::Integer = buflen)
    if text === missing
        text = pyhasattr(o, "encoding")
    end
    buflen = convert(Int, buflen)
    buflen > 0 || error("buflen must be positive")
    ibuflen = convert(Int, ibuflen)
    ibuflen > 0 || error("ibuflen must be positive")
    obuflen = convert(Int, obuflen)
    obuflen > 0 || error("obuflen must be positive")
    PyIO(Val(:new), Py(o), own, text, ibuflen, obuflen)
end

pyio_finalize!(io::PyIO) = begin
    C.CTX.is_initialized || return
    io.own ? close(io) : flush(io)
    pydel!(io.py)
    return
end

ispy(io::PyIO) = true
getpy(io::PyIO) = io.py
pydel!(io::PyIO) = (finalize(io); nothing)

pyconvert_rule_io(::Type{PyIO}, x::Py) = pyconvert_return(PyIO(x))

"""
    PyIO(f, x; ...)

Equivalent to `f(PyIO(x; ...))` except the stream is automatically flushed or closed according to `own`.

For use in a `do` block, as in
```
PyIO(x, ...) do io
    ...
end
```
"""
function PyIO(f::Function, o; opts...)
    io = PyIO(o; opts...)
    try
        return f(io)
    finally
        pydel!(io)
    end
end

# If obuf is non-empty, write it to the underlying stream.
function putobuf(io::PyIO)
    if !isempty(io.obuf)
        data = io.text ? pystr_fromUTF8(io.obuf) : pybytes(io.obuf)
        pydel!(@py io.write(data))
        pydel!(data)
        empty!(io.obuf)
    end
    return
end

# If ibuf is empty, read some more from the underlying stream.
# After this call, if ibuf is empty then we are at EOF.
# TODO: in binary mode, `io.readinto()` to avoid copying data
function getibuf(io::PyIO)
    if isempty(io.ibuf)
        data = @py io.read(@jl io.ibuflen)
        if io.text
            append!(io.ibuf, pystr_asUTF8vector(data))
        else
            append!(io.obuf, pybytes_asvector(data))
        end
        pydel!(data)
    end
    return
end

function Base.flush(io::PyIO)
    putobuf(io)
    pydel!(@py io.flush())
    return
end

function Base.close(io::PyIO)
    flush(io)
    pydel!(@py io.close())
    return
end

function Base.eof(io::PyIO)
    if io.eof
        true
    elseif !isempty(io.ibuf)
        false
    else
        getibuf(io)
        io.eof = isempty(io.ibuf)
    end
end

Base.fd(io::PyIO) = pyconvert_and_del(Int, @py io.fileno())

Base.isreadable(io::PyIO) = pyconvert_and_del(Bool, @py io.readable())

Base.iswritable(io::PyIO) = pyconvert_and_del(Bool, @py io.writable())

Base.isopen(io::PyIO) = !pyconvert_and_del(Bool, @py io.closed)

function Base.unsafe_write(io::PyIO, ptr::Ptr{UInt8}, n::UInt)
    ntodo = n
    while true
        nroom = max(0, io.obuflen - length(io.obuf))
        if ntodo < nroom
            append!(io.obuf, unsafe_wrap(Array, ptr, ntodo))
            return n
        else
            append!(io.obuf, unsafe_wrap(Array, ptr, nroom))
            putobuf(io)
            ptr += nroom
            ntodo -= nroom
        end
    end
end

function Base.write(io::PyIO, c::UInt8)
    push!(io.obuf, c)
    if length(io.obuf) ≥ io.obuflen
        putobuf(io)
    end
    return
end

function Base.unsafe_read(io::PyIO, ptr::Ptr{UInt8}, n::UInt)
    ntodo = n
    while true
        navail = length(io.ibuf)
        if ntodo ≤ navail
            unsafe_wrap(Array, ptr, ntodo) .= splice!(io.ibuf, 1:ntodo)
            return
        else
            unsafe_wrap(Array, ptr, navail) .= io.ibuf
            ptr += navail
            ntodo -= navail
            empty!(io.ibuf)
            eof(io) && throw(EOFError())
        end
    end
    return
end

function Base.read(io::PyIO, ::Type{UInt8})
    eof(io) && throw(EOFError())
    popfirst!(io.ibuf)
end

function Base.seek(io::PyIO, pos::Integer)
    putobuf(io)
    empty!(io.ibuf)
    io.eof = false
    pydel!(@py io.seek(pos))
    return io
end

function Base.truncate(io::PyIO, pos::Integer)
    seek(io, position(io))
    pydel!(@py io.truncate(pos))
    return io
end

function Base.seekstart(io::PyIO)
    putobuf(io)
    empty!(io.ibuf)
    io.eof = false
    pydel!(@py io.seek(0))
    return io
end

function Base.seekend(io::PyIO)
    putobuf(io)
    empty!(io.ibuf)
    io.eof = false
    pydel!(@py io.seek(0, 2))
    return io
end

function Base.skip(io::PyIO, n::Integer)
    putobuf(io)
    if io.text
        if n ≥ 0
            read(io, n)
        else
            error("`skip(io, n)` for text PyIO streams only implemented for positive `n`")
        end
    else
        if 0 ≤ n ≤ io.ibuflen
            read(io, n)
        else
            pydel!(@py io.seek(@jl(n - length(io.ibuf)), 1))
            empty!(io.ibuf)
            io.eof = false
        end
    end
    return io
end

function Base.position(io::PyIO)
    putobuf(io)
    if io.text
        if isempty(io.ibuf)
            return pyconvert_and_del(Int, @py io.tell())
        else
            error("`position(io)` text PyIO streams only implemented for empty input buffer (e.g. do `read(io, length(io.ibuf))` first)")
        end
    else
        return pyconvert_and_del(Int, @py io.tell()) - length(io.ibuf)
    end
end


text:
"""
This struct contains the data for each scenario.
    name: Scenario name
    probability: Probability of the scenario
    parameter_multipliers: Parameter multipliers for each scenario - to be multiplied to updated parameters from kalman filters
    parameter_values: Dictionary containing parameter values for the scenario
"""
mutable struct Scenario
    name::String
    probability::Float64
    parameter_multipliers::Union{Nothing, Dict{String, Float64}}
    parameter_values::Vector{AxisArrays.AxisArray{Float64, 2}}
end

get_name(scenario::Scenario) = scenario.name
get_probability(scenario::Scenario) = scenario.probability
get_parameter_multipliers(scenario::Scenario) = scenario.parameter_multipliers
get_parameter_values(scenario::Scenario) = scenario.parameter_values

function set_parameter_values!(scenario::Scenario, iteration_year::Int64, parameter_values::AxisArrays.AxisArray{Float64, 2})
    scenario.parameter_values[iteration_year] = parameter_values
end


text:
# Data summary pass:
#
# Before we can construct some of the term types, we need to know additional
# invariants about the data, like the unique values of categorical data, or
# min/max values for splines.  This pass will consume the data once and compute
# summaries of variables where necessary.
#
# For now, this is just for categorical terms.

include("abstract.jl")
include("categorical.jl")
include("make.jl")
include("run.jl")

# TODO: use OnlineStatsBase.jl API: replace Summarizer with OnlineStat, which
# needs to provide a fit!, merge, and _value method.


text:
module BGZFStreams

export
    BGZFStream,
    VirtualOffset,
    BGZFDataError,
    virtualoffset

import CodecZlib

include("virtualoffset.jl")
include("bgzfstream.jl")

end # module


text:
## Occupation vector iteration and array properties.

function Base.iterate(basis::AbstractSzbasis, state=1)
    state > basis.D && return nothing
    @view(basis.vectors[:, state]), state+1
end

Base.eltype(::Type{AbstractSzbasis}) = Vector{Int}
Base.length(basis::AbstractSzbasis) = basis.D

function Base.in(v::AbstractVector{Int}, basis::Szbasis)
    length(v) == basis.K && sum(v) == basis.N
end

function Base.in(v::AbstractVector{Int}, basis::RestrictedSzbasis)
    length(v) == basis.K && sum(v) == basis.N && maximum(v) <= basis.M
end


text:
function add_celestial_body_loader(::Type{CelestialBodyFactory}, arg0::JString, arg1::CelestialBodyLoader)
    return jcall(CelestialBodyFactory, "addCelestialBodyLoader", void, (JString, CelestialBodyLoader), arg0, arg1)
end

function add_default_celestial_body_loader(::Type{CelestialBodyFactory}, arg0::JString)
    return jcall(CelestialBodyFactory, "addDefaultCelestialBodyLoader", void, (JString,), arg0)
end

function add_default_celestial_body_loader(::Type{CelestialBodyFactory}, arg0::JString, arg1::JString)
    return jcall(CelestialBodyFactory, "addDefaultCelestialBodyLoader", void, (JString, JString), arg0, arg1)
end

function clear_celestial_body_cache(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "clearCelestialBodyCache", void, ())
end

function clear_celestial_body_cache(::Type{CelestialBodyFactory}, arg0::JString)
    return jcall(CelestialBodyFactory, "clearCelestialBodyCache", void, (JString,), arg0)
end

function clear_celestial_body_loaders(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "clearCelestialBodyLoaders", void, ())
end

function clear_celestial_body_loaders(::Type{CelestialBodyFactory}, arg0::JString)
    return jcall(CelestialBodyFactory, "clearCelestialBodyLoaders", void, (JString,), arg0)
end

function get_body(::Type{CelestialBodyFactory}, arg0::JString)
    return jcall(CelestialBodyFactory, "getBody", CelestialBody, (JString,), arg0)
end

function get_celestial_bodies(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getCelestialBodies", LazyLoadedCelestialBodies, ())
end

function get_earth(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getEarth", CelestialBody, ())
end

function get_earth_moon_barycenter(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getEarthMoonBarycenter", CelestialBody, ())
end

function get_jupiter(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getJupiter", CelestialBody, ())
end

function get_mars(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getMars", CelestialBody, ())
end

function get_mercury(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getMercury", CelestialBody, ())
end

function get_moon(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getMoon", CelestialBody, ())
end

function get_neptune(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getNeptune", CelestialBody, ())
end

function get_pluto(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getPluto", CelestialBody, ())
end

function get_saturn(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getSaturn", CelestialBody, ())
end

function get_solar_system_barycenter(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getSolarSystemBarycenter", CelestialBody, ())
end

function get_sun(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getSun", CelestialBody, ())
end

function get_uranus(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getUranus", CelestialBody, ())
end

function get_venus(::Type{CelestialBodyFactory})
    return jcall(CelestialBodyFactory, "getVenus", CelestialBody, ())
end


text:
@testset "Coordinate" begin
    for i in 1:10
        @test typeof(Coordinate(rand(-10.0:10.0), rand(-10.0:10.0))) == Coordinate

    end
end

@testset "Coordinate System" begin
    for i in 1:10
        top_left = Coordinate(rand(30.0:50.0), rand(0:10.0))
        bottom_right = Coordinate(rand(0.0:30.0), rand(10.0:20.0))
        height = rand(1:1000)
        width = rand(1:1000)
        system = CoordinateSystem(top_left, bottom_right, height, width)
        @test typeof(system) == CoordinateSystem
    end
end

@testset "Conversion functions" begin
    @test lat_to_row(INDIA_COORDINATE_SYSTEM, 19.6) == 4296
    @test long_to_column(INDIA_COORDINATE_SYSTEM, 73.33) == 1299
    @test isapprox(row_to_lat(INDIA_COORDINATE_SYSTEM, 100), 37.083325)
    @test isapprox(column_to_long(INDIA_COORDINATE_SYSTEM, 100), 68.333326760)
    @test coordinate_to_image(INDIA_COORDINATE_SYSTEM, Coordinate(19.49907, 72.721252)) == [4320, 1153]
    @test isapprox(image_to_coordinate(INDIA_COORDINATE_SYSTEM, [4320, 1153]).latitude, Coordinate(19.49964, 72.720827).latitude)
    @test isapprox(image_to_coordinate(INDIA_COORDINATE_SYSTEM, [4320, 1153]).longitude, Coordinate(19.49964, 72.720827).longitude)
    @test translate_geometry(INDIA_COORDINATE_SYSTEM, Coordinate(19.49907, 72.721252), Coordinate(18.849475, 73.074187)) == MUMBAI_COORDINATE_SYSTEM
end


text:
using Documenter, ResultTypes

makedocs(;
    modules=[ResultTypes],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "API" => "pages/api.md",
    ],
    repo="https://github.com/iamed2/ResultTypes.jl/blob/{commit}{path}#L{line}",
    sitename="ResultTypes.jl",
    authors="Eric Davies",
    assets=[],
)

deploydocs(;
    repo="github.com/iamed2/ResultTypes.jl",
)


text:
# Copyright 2017-2019 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using SimilaritySearch:
    Sequential, KnnResult
import SimilaritySearch: optimize!

export NearNeighborClassifier, optimize!, predict, predict_proba

mutable struct NearNeighborClassifier{IndexType,LabelType}
    dist::Function
    X::IndexType
    y::Vector{Int}
    k::Int
    le::LabelEncoder{LabelType}
    weight

    function NearNeighborClassifier(dist::Function, X::AbstractVector{ItemType}, y::AbstractVector{LabelType}, k::Int=1, weight=:uniform, create_index=nothing) where {ItemType, LabelType}
        le = fit(LabelEncoder, y)
        y_ = [transform(le, l) for l in y]

        if create_index == nothing
            index = fit(Sequential, X)
        else
            index = create_index(X)
        end

        new{typeof(index), LabelType}(dist, index, y_, k, le, weight)
    end
end

function predict(nnc::NearNeighborClassifier{IndexType,LabelType}, vector) where {IndexType,LabelType}
    [predict_one(nnc, item) for item in vector]
end

function predict_proba(nnc::NearNeighborClassifier{IndexType,LabelType}, vector::AbstractVector; smoothing=0.0) where {IndexType,LabelType}
    [predict_one_proba(nnc, item, smoothing=smoothing) for item in vector]
end

function _predict_one(nnc::NearNeighborClassifier{IndexType,LabelType}, item) where {IndexType,LabelType}
    res = KnnResult(nnc.k)
    search(nnc.X, nnc.dist, item, res)
    w = zeros(Float64, length(nnc.le.labels))
    if nnc.weight == :uniform
        for p in res
            l = nnc.y[p.objID]
            w[l] += 1.0
        end
    elseif nnc.weight == :distance
        for p in res
            l = nnc.y[p.objID]
            w[l] += 1.0 / (1.0 + p.dist)
        end
    else
        throw(ArgumentError("Unknown weighting scheme $(nnc.weight)"))
    end

    w
end

function predict_one(nnc::NearNeighborClassifier{IndexType,LabelType}, item) where {IndexType,LabelType}
    score, i = findmax(_predict_one(nnc, item))
    inverse_transform(nnc.le, i)
end

function predict_one_proba(nnc::NearNeighborClassifier{IndexType,LabelType}, item; smoothing=0.0) where {IndexType,LabelType}
    w = _predict_one(nnc, item)
    t = sum(w)

    s = t * smoothing
    ss = s * length(w)

    for i in 1:length(w)
        w[i] = (w[i] + s) / (t + ss)  # overriding previous w
    end

    w
end

function _train_create_table(dist::Function, train_X, train_y, test_X, k::Int)
    index = fit(Sequential, train_X)
    tab = Vector{Vector{Tuple{Int,Float64}}}(undef, length(test_X))

    for i in 1:length(test_X)
        res = search(index, dist, test_X[i], KnnResult(k))
        tab[i] = [(train_y[p.objID], p.dist) for p in res]
    end

    tab
end

function _train_predict(nnc::NearNeighborClassifier{IndexType,LabelType}, table, test_X, k) where {IndexType,LabelType}
    A = Vector{LabelType}(undef, length(test_X))
    w = Vector{Float64}(undef, length(nnc.le.labels))

    for i in 1:length(test_X)
        w .= 0.0
        row = table[i]
        for j in 1:k
            label, _d = row[j]
            if nnc.weight == :uniform
                w[label] += 1
            else
                w[label] += 1.0 / (1.0 + _d)
            end
        end

        score, label = findmax(w)
        A[i] = inverse_transform(nnc.le, label)
    end

    A
end

function optimize!(nnc::NearNeighborClassifier, scorefun::Function; runs=3, trainratio=0.5, testratio=0.5, folds=0, folds_shuffle=true)
    mem = Dict{Tuple,Float64}()

    function f(train_X, train_y, test_X, test_y)
        _nnc = NearNeighborClassifier(nnc.dist, train_X, train_y)
        kmax = sqrt(length(train_y)) |> round |> Int
        table = _train_create_table(_nnc.dist, train_X, train_y, test_X, kmax)
        k = 2
        while k <= kmax
            for weight in (:uniform, :distance)
                _nnc.weight = weight
                _nnc.k = k - 1
                pred_y = _train_predict(_nnc, table, test_X, _nnc.k)
                score = scorefun(test_y, pred_y)
                key = (k - 1, weight)
                mem[key] = get(mem, key, 0.0) + score
            end
            k += k
        end
        0
    end

    if folds > 1
        kfolds(f, nnc.X.db, nnc.y, folds=folds, shuffle=folds_shuffle)
        bestlist = [(score/folds, conf) for (conf, score) in mem]
    else
        montecarlo(f, nnc.X.db, nnc.y, runs=runs, trainratio=trainratio, testratio=testratio)
        bestlist = [(score/runs, conf) for (conf, score) in mem]
    end

    sort!(bestlist, by=x -> (-x[1], x[2][1]))
    best = bestlist[1]
    nnc.k = best[2][1]
    nnc.weight = best[2][2]

    bestlist
end


text:
include("src/impedance/impedance.jl")
#include("src/constants.jl")
#using .constants
#include("docs/make.jl")

using .impedance

#f = 10; σ = 200_000; α = 5; Ω = 0.5; q = 1.3; sf=0.85; σs=0.3; s̄=46*10^-7

f = 1000         # Hz, frequency
σ = 200_0     # Pa.s.m-2, flow resistivity
Ω = 0.5         # porosity

Zc1 = impedanceM1(f,σ)      # tested!

Zc2 = impedanceM2(f,σ)      # tested!

Zc3 = impedanceM3(f,σ,Ω)    # tested!

Zc4 = impedanceM4(f,σ,Ω)    # tested!

Zc5 = impedanceM5(f,σ,Ω)    # tested!

Zc6 = impedanceM6(f,σ,Ω)    # tested!

Zc7 = impedanceM7(f,σ,Ω)    # tested!

Zc8 = impedanceM8(f,σ,Ω)    # tested!

Zc9b = impedanceM9b(f,σ,Ω)


text:
# Prediction
# ==========
#
# Sequence type predictor.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# predict sequence type based on character frequencies in `seq[start:stop]`
function predict(seq::Vector{UInt8}, start, stop)
    # count characters
    a = c = g = t = u = n = alpha = 0
    for i in start:stop
        x = seq[i]
        if x == 0x41 || x == 0x61
            a += 1
        elseif x == 0x43 || x == 0x63
            c += 1
        elseif x == 0x47 || x == 0x67
            g += 1
        elseif x == 0x54 || x == 0x74
            t += 1
        elseif x == 0x55 || x == 0x75
            u += 1
        elseif x == 0x4e || x == 0x6e
            n += 1
        end
        if 0x41 ≤ x ≤ 0x5a || 0x61 ≤ x ≤ 0x7a
            alpha += 1
        end
    end

    # the threshold (= 0.95) is somewhat arbitrary
    if (a + c + g + t + u + n) / alpha > 0.95
        if t ≥ u
            return DNASequence
        else
            return RNASequence
        end
    else
        return AminoAcidSequence
    end
end


text:
include("../utils.jl")
include("group_utils.jl")

@testset "Special Orthogonal group" begin
    G = SpecialOrthogonal(3)
    @test repr(G) == "SpecialOrthogonal(3)"
    M = base_manifold(G)
    @test M === Rotations(3)
    x = Matrix(I, 3, 3)

    types = [Matrix{Float64}]
    ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]]
    pts = [exp(M, x, hat(M, x, ωi)) for ωi in ω]
    vpts = [hat(M, x, [-1.0, 2.0, 0.5]), hat(M, x, [1.0, 0.0, 0.5])]

    ge = similar(pts[1])
    copyto!(ge, Identity(G))
    @test isapprox(ge, I; atol=1e-10)

    for T in types
        gpts = convert.(T, pts)
        vgpts = convert.(T, vpts)
        @test compose(G, gpts[1], gpts[2]) ≈ gpts[1] * gpts[2]
        @test translate_diff(G, gpts[2], gpts[1], vgpts[1], LeftAction()) ≈ vgpts[1]
        @test translate_diff(G, gpts[2], gpts[1], vgpts[1], RightAction()) ≈ transpose(gpts[2]) * vgpts[1] * gpts[2]
        test_group(G, gpts, vgpts; test_diff = true)
    end

    @testset "Decorator forwards to group" begin
        DM = NotImplementedGroupDecorator(G)
        @test Manifolds.is_decorator_group(DM) === Val(true)
        @test base_group(DM) === G
        @test Identity(DM) === Identity(G)
        @test_throws DomainError is_manifold_point(DM, Identity(TranslationGroup(3)), true)
        test_group(DM, pts, vpts; test_diff = true)
    end

    @testset "Group forwards to decorated" begin
        retraction_methods = [Manifolds.PolarRetraction(),
                              Manifolds.QRRetraction()]

        inverse_retraction_methods = [Manifolds.PolarInverseRetraction(),
                                      Manifolds.QRInverseRetraction()]

        test_manifold(G, pts;
            test_reverse_diff = false,
            test_injectivity_radius = false,
            test_project_tangent = true,
            test_musical_isomorphisms = false,
            retraction_methods = retraction_methods,
            inverse_retraction_methods = inverse_retraction_methods,
            exp_log_atol_multiplier = 20,
            retraction_atol_multiplier = 12,
        )

        @test injectivity_radius(G) == injectivity_radius(M)
        @test injectivity_radius(G, pts[1]) == injectivity_radius(M, pts[1])
        @test injectivity_radius(G, pts[1], PolarRetraction()) == injectivity_radius(M, pts[1], PolarRetraction())

        y = similar(pts[1])
        exp!(G, y, pts[1], vpts[1])
        @test isapprox(M, y, exp(M, pts[1], vpts[1]))

        y = similar(pts[1])
        retract!(G, y, pts[1], vpts[1])
        @test isapprox(M, y, retract(M, pts[1], vpts[1]))

        w = similar(vpts[1])
        inverse_retract!(G, w, pts[1], pts[2])
        @test isapprox(M, pts[1], w, inverse_retract(M, pts[1], pts[2]))

        w = similar(vpts[1])
        inverse_retract!(G, w, pts[1], pts[2], QRInverseRetraction())
        @test isapprox(M, pts[1], w, inverse_retract(M, pts[1], pts[2], QRInverseRetraction()))

        z = collect(reshape(1.0:9.0, 3, 3))
        @test isapprox(M, project_point(G, z), project_point(M, z))
        z2 = similar(pts[1])
        project_point!(G, z2, z)
        @test isapprox(M, z2, project_point(M, z))
    end
end


text:
export import_from_dat

"""
    pb, point = import_from_dat(instancepath::String, precondcstrspath::String)

Build the polynomial optimization problem described by the `instancepath` file,
with the dispensory preconditionning descritpion `precondcstrspath` along with
the initial point possibly provided in the file (defaults value is null).
"""
function import_from_dat(instancepath::String, precondfilename::String="")
    point = Point()
    variables = SortedDict{String, Variable}()
    exponents = SortedDict{String, Exponent}()
    pb = Problem()

    instance_str = open(instancepath)
    l = jump_comments!(instance_str)

    ## Collect and define variables
    # line = matchall(r"\S+", l)
    line = collect(( m.match for m=eachmatch(r"\S+", l) ))
    while line[1] == "VAR_TYPE" && !eof(instance_str)
        if line[2] == "REAL"
            var = Variable(line[3], Real)
        elseif line[2] == "BOOL"
            var = Variable(line[3], Bool)
        elseif line[2] == "C"
            var = Variable(line[3], Complex)
        else
            error(LOGGER, "import_to_dat(): Unknown variable type $(line[2]) for variable $(line[3]).")
        end

        variables[line[3]] = var
        add_variable!(pb, var)

        setindex!(point, convert(var.kind, parse.(Float64, line[5]) + im*parse.(Float64,line[6])), var)

        ## Mark where objective definition begins
        mark(instance_str)

        l = readline(instance_str)
        line = collect(( m.match for m=eachmatch(r"\S+", l) ))
    end

    ## Move forward to the monomial definition section
    while line[1] != "MONO_DEF" && !eof(instance_str)
        l = readline(instance_str)
        line = collect(( m.match for m=eachmatch(r"\S+", l) ))
    end

    ## Build all monomials
    while line[1] == "MONO_DEF" && !eof(instance_str)
        exponame = line[2]
        var = variables[line[3]]
        if !haskey(exponents, exponame)
            exponents[exponame] = Exponent()
        end
        exponents[exponame] = product(exponents[exponame], Exponent(SortedDict(var=>Degree(parse.(Float64,line[5]), parse.(Float64,line[6])))))
        l = readline(instance_str)
        line = collect(( m.match for m=eachmatch(r"\S+", l) ))
    end

    ## Add final monomial
    if line[1] == "MONO_DEF"
        exponame = line[2]
        var = variables[line[3]]
        if !haskey(exponents, exponame)
            exponents[exponame] = Exponent()
        end
        exponents[exponame] = product(exponents[exponame], Exponent(SortedDict(var=>Degree(parse.(Float64,line[5]), parse.(Float64,line[6])))))
    end

    ## Reset stream to the objective defintion
    reset(instance_str)
    l = readline(instance_str)
    line = collect(( m.match for m=eachmatch(r"\S+", l) ))

    ## Build polynomial objective
    p = Polynomial()
    while line[2] == "OBJ"
        λ = parse_λ(line[5], line[6])
        var1, var2 = line[3:4]
        if line[1] == "MONO"
            add!(p, λ * exponents[var1])
        else
            quad_expo = Exponent()
            (var1!="NONE") && product!(quad_expo, conj(variables[var1]))
            (var2!="NONE") && product!(quad_expo, variables[var2])
            add!(p, Polynomial(SortedDict{Exponent, Number}(quad_expo=>λ)))
        end
        l = readline(instance_str)
        line = collect(( m.match for m=eachmatch(r"\S+", l) ))
    end
    set_objective!(pb, p)

    ## Build constraints
    next_state = :AssembleCtr
    cur_ctr = line[2]
    var1, var2 = line[3:4]

    lb = -Inf
    ub = +Inf
    p = Polynomial()
    λ = parse_λ(line[5], line[6])
    while next_state!=:ReadLine || (!eof(instance_str) && line[1] != "MONO_DEF")

        state = next_state
        if state == :ReadLine
            ## Readline
            l = readline(instance_str)
            line = collect(( m.match for m=eachmatch(r"\S+", l) ))

            var1, var2 = line[3:4]
            λ = parse_λ(line[5], line[6])

            ## Determine whether current ctr should be completed or new ctr should be set
            if line[2] == cur_ctr
                next_state = :AssembleCtr
            else
                next_state = :SaveCtr
            end

            ## If monomial difnition section reached, exit loop
            if line[1] == "MONO_DEF"
                next_state = :ReadLine
            end

        elseif state == :AssembleCtr
            if line[1] == "MONO"
                add!(p, λ * exponents[var1])
            elseif line[1] ∈ SortedSet(["CONST", "LIN", "QUAD"])
                quad_expo = Exponent()
                (var1!="NONE") && product!(quad_expo, conj(variables[var1]))
                (var2!="NONE") && product!(quad_expo, variables[var2])
                add!(p, Polynomial(SortedDict{Exponent, Number}(quad_expo=>λ)))
            elseif line[1] == "UB"
                ub = λ
            elseif line[1] == "LB"
                lb = λ
            else
                error(LOGGER, "import_from_dat(): Unknown variable type $(line[1]) for constraint $(line[2]).")
            end

            next_state = :ReadLine

        elseif state == :SaveCtr
            add_constraint!(pb, String(cur_ctr), lb << p << ub)

            lb = -Inf
            ub = +Inf
            p = Polynomial()
            cur_ctr = line[2]
            next_state = :AssembleCtr
        else
            error(LOGGER, "import_from_dat(): Unknown state $state")
        end
    end

    # Add final constraint
    add_constraint!(pb, String(cur_ctr), lb << p << ub)

    close(instance_str)

    ## Set preconditioning flag
    if precondfilename != ""
        precond_str = open(joinpath(splitdir(instancepath)[1], precondfilename))

        l = jump_comments!(precond_str)
        line = collect(( m.match for m=eachmatch(r"\S+", l) ))
        while !eof(precond_str)
            if line[2] == "SQRT"
                pb.constraints[line[1]].precond = :sqrt
            else
                error(LOGGER, "import_from_dat(): Unknown preconditioning $(line[2]) for constraint $(line[1]).")
            end
            l = readline(precond_str)
            line = collect(( m.match for m=eachmatch(r"\S+", l) ))
        end
        close(precond_str)
    end

    return pb, point
end

"""
    l = jump_comments!(io::IOStream)

Jump comments, i.e. lines starting by '#', in the `io` stream, and return the
first non commented line.
"""
function jump_comments!(io::IOStream)
    l = ""
    while !eof(io)
        l = readline(io)
        occursin(r"\s*#", l) || break
    end
    return l
end

function parse_λ(realpart::T, imagpart::T) where T <: Union{String, SubString}
    λ = 0
    if realpart == "Inf"
        λ += Inf
    elseif realpart == "-Inf"
        λ += -Inf
    else
        λ += parse.(Float64,realpart)
    end

    if imagpart == "Inf"
        λ += Inf*im
    elseif imagpart == "-Inf"
        λ += -Inf*im
    else
        λ += parse.(Float64,imagpart)*im
    end
    return λ
end


text:
struct Interpolator
    interpolation::Interpolations.BSplineInterpolation

    function Interpolator(field::Array{<:Any,N}) where {N,D}
        new(Interpolations.interpolate(field, BSpline(Linear())))
    end
end

function interpolate(interpolator::Interpolator, continuous_indices)
    interpolator.interpolation(continuous_indices...)
end

max_stars_repo_name: lars-frogner/PIC

max_issues_repo_name: lars-frogner/PIC

max_forks_repo_name: lars-frogner/PIC


text:
export export_to_dat
export export_matpower_to_dat

function get_varsconj(exp::Exponent)
  varconj = Vector{Variable}()
  for (var, expo) in exp
    for i=1:expo.conjvar
      push!(varconj, var)
    end
  end
  varconj
end
function get_varslin(exp::Exponent)
  varlin = Vector{Variable}()
  for (var, expo) in exp
    for i=1:expo.explvar
      push!(varlin, var)
    end
  end
  varlin
end

function collect_requiredmonomials!(p::Polynomial, expos::SortedDict{Exponent, String})
  for expo in keys(p)
    if !is_oneline_printable(expo) && !haskey(expos, expo)
      expos[expo] = "MONO_$(length(expos))"
    end
  end
end

function print_string(io, strng, len; alignright = true, indentedprint=true)
  indentedprint || (len = length(strng))
  if alignright
    print(io, " "^(len-length(strng)), strng, " ")
  else
    print(io, strng, " "^(len-length(strng)), " ")
  end
end

function print_dat_line(io, linetype, cstrname, var1, var2, val1, val2, maxvarlen, maxcstrlen)
  @printf(io, "%8s ", linetype)
  print_string(io, cstrname, maxcstrlen)
  print_string(io, var1, maxvarlen)
  print_string(io, var2, maxvarlen)
  @printf(io, "% 23.16e % 23.16e\n", val1, val2)
end

function print_quad_expo(io, expo::Exponent, cat::String, coeff, maxvarlen, maxcstrlen)
  vars_conj = get_varsconj(expo)
  vars_lin = get_varslin(expo)
  if length(vars_conj) == 1 && length(vars_lin) == 1
    print_dat_line(io, "QUAD", cat, vars_conj[1].name, vars_lin[1].name, real(coeff), imag(coeff), maxvarlen, maxcstrlen)
  elseif length(vars_conj) == 0 && length(vars_lin) == 2 && vars_lin[1].kind <: Real && vars_lin[2].kind <: Real
    print_dat_line(io, "QUAD", cat, vars_lin[1].name, vars_lin[2].name, real(coeff), imag(coeff), maxvarlen, maxcstrlen)
  elseif length(vars_conj) == 1 && length(vars_lin) == 0
    print_dat_line(io, "LIN", cat, vars_conj[1].name, "NONE", real(coeff), imag(coeff), maxvarlen, maxcstrlen)
  elseif length(vars_conj) == 0 && length(vars_lin) == 1
    print_dat_line(io, "LIN", cat, "NONE", vars_lin[1].name, real(coeff), imag(coeff), maxvarlen, maxcstrlen)
  elseif length(vars_conj) == 0 && length(vars_lin) == 0
    print_dat_line(io, "CONST", cat, "NONE", "NONE", real(coeff), imag(coeff), maxvarlen, maxcstrlen)
  else
    warn(LOGGER, "print_quad_expo(): Exponent $expo not supported.")
  end
end

function print_constraint(io::IO, cstrname::String, cstr::Constraint, maxvarlen::Int, maxcstrlen::Int, expos::SortedDict{Exponent, String})
  print_poly!(io, cstr.p, cstrname, maxvarlen, maxcstrlen, expos)

  if occursin(r"_Re", cstrname)
    if real(cstr.lb) != -Inf
      print_dat_line(io, "LB", cstrname, "NONE", "NONE", real(cstr.lb), imag(cstr.lb), maxvarlen, maxcstrlen)
    end
    if real(cstr.ub) != Inf
      print_dat_line(io, "UB", cstrname, "NONE", "NONE", real(cstr.ub), imag(cstr.ub), maxvarlen, maxcstrlen)
    end
  elseif occursin(r"_Im", cstrname)
    if imag(cstr.lb) != -Inf
      print_dat_line(io, "LB", cstrname, "NONE", "NONE", real(cstr.lb), imag(cstr.lb), maxvarlen, maxcstrlen)
    end
    if imag(cstr.ub) != Inf
      print_dat_line(io, "UB", cstrname, "NONE", "NONE", real(cstr.ub), imag(cstr.ub), maxvarlen, maxcstrlen)
    end
  else
    if imag(cstr.lb) != -Inf-im*Inf
      print_dat_line(io, "LB", cstrname, "NONE", "NONE", real(cstr.lb), imag(cstr.lb), maxvarlen, maxcstrlen)
    end
    if imag(cstr.ub) != +Inf+im*Inf
      print_dat_line(io, "UB", cstrname, "NONE", "NONE", real(cstr.ub), imag(cstr.ub), maxvarlen, maxcstrlen)
    end
  end
end
nb_from_str(string::String) = match(r"\d+", string).match
get_scenario(string::String) = String(split(string, "_")[1])

function print_doc(io, filename)
  println(io, "# $filename - exported $(Dates.now())")
  println(io, "# Description of a QCQP optimization problem. 3 sections:")
  println(io, "#  - List of all variables in the problem, along with type, name, and possibly complex value.")
  println(io, "#    Grouped by \"VAR_TYPE\" tag (col 1). Type is \"C\" Complex, \"REAL\" Real, \"BOOL\" Boolean (col 2). ")
  println(io, "#    Name (col 3). Real and imag part of value (col 5 and 6 resp.).")
  println(io, "#  - Description of the objective function (one quadratic polynomial).")
  println(io, "#    Sum of monomials, grouped by \"OBJ\" (col 2), either:")
  println(io, "#                order 0 (\"CONST\" col 1) : real and imag part of coef. resp. at col 5 and 6.")
  println(io, "#                order 1 (\"LIN\" col 1) : real and imag part of coef. resp. at col 5 and 6, either:")
  println(io, "#                         variable name in col 3 hence conjugate variable in polynomial,")
  println(io, "#                         variable name in col 4 hence variable in polynomial.")
  println(io, "#                order 2 (\"QUAD\" col 1) : real and imag part of coef. resp. at col 5 and 6, ")
  println(io, "#                         monomial is product of conjugate variable at col 3 and variable at col 4.")
  println(io, "#  - Description of the constraints : one polynomial and two complex numeric bounds.")
  println(io, "#    Constraints are grouped by name (col 2). Quadratic body is described as the objective.")
  println(io, "#    Lower bound (\"LB\" col 1) and upper bound (\"UB\" col 1) have their complex value in col 5 and 6.")
  println(io, "#")
end

"""
  print_variables(io::IO, variables, pt::Point)

  Write a .dat description of `variables` variables to `io`.
"""
function print_variables(io::IO, variables, pt::Point, maxvarlen, maxcstrlen)
  for (varname, varkind) in variables
    var = Variable(varname, varkind)
    if iscomplex(var)
      var_type = "C"
    elseif isbool(var)
      var_type = "BOOL"
    elseif isreal(var)
      var_type = "REAL"
    else
      error(LOGGER, "Export_to_dat(): unsuported variable type $(var.kind)")
    end
    val = 0

    if haskey(pt, var)
      val = pt[var]
    end

    print_dat_line(io, "VAR_TYPE", var_type, var.name, "NONE", real(val), imag(val), maxvarlen, maxcstrlen)
  end
end

"""
  print_poly!(io::IO, p::AbstractPolynomial, cat::String, maxvarlen, maxcstrlen, expos::SortedDict{Exponent, String})

Print the `p` polynomial corresponding to the constraint or objective
`cat` (category) to `io`. Each of the polynomial's exponent is printed in a line,
either explicitly if its degree is 2 or less, or implicitly by defining an
exponent name in the `expos` dict if it required, and print the exponent name
with the coefficient (which essentially is a linear term).
"""
function print_poly!(io::IO, p::AbstractPolynomial, cat::String, maxvarlen, maxcstrlen, expos::SortedDict{Exponent, String})
  constval = 0

  for (expo, coeff) in p
    oneline = is_oneline_printable(expo)

    if length(expo) == 0  # const value
      constval = coeff
    elseif oneline         # oneline printable monomial
      print_quad_expo(io, expo, cat, coeff, maxvarlen, maxcstrlen)
    else                  # general monomial case
      print_dat_line(io, "MONO", cat, expos[expo], "NONE", real(coeff), imag(coeff), maxvarlen, maxcstrlen)
    end
  end
    if constval != 0
      print_dat_line(io, "CONST", cat, "NONE", "NONE", real(constval), imag(constval), maxvarlen, maxcstrlen)
  end
end

function is_oneline_printable(expo::Exponent)
  explsum, conjsum = get_sumdegs(expo)
  vars_deg = collect(expo.expo)

  oneline = (length(expo) == 0)
  oneline = (oneline || ((length(expo) == 2) && ((explsum, conjsum) == (1,1)))) # Hermitian product of two complex variables
  oneline = (oneline || ((length(expo) == 2) && ((explsum, conjsum) == (2,0)) && isreal(first(vars_deg)[1]) && isreal(last(vars_deg)[1]))) # Product of two different real variables
  oneline = (oneline || ((length(expo) == 1) && ((explsum, conjsum) == (2,0)) && isreal(first(vars_deg)[1]))) # One squared real variable
  var = ((length(expo) == 1) && (((explsum, conjsum) == (0,1)) || ((explsum, conjsum) == (1,0)))) # One real or complex variable
  oneline = (var || oneline)
  return oneline
end

"""
  export_to_dat(pb_optim::Problem, outpath; filename="real_minlp_instance.dat", pt=Point(), exportprecond=false)

Write the `pb_optim` problem to the `outpath` folder, with the dispensory initial
point `pt`. Output files are `real_minlp_instance.dat` and
`real_minlp_precond_cstrs.dat`.
"""
function export_to_dat(pb_optim::Problem, outpath::String; filename::String="real_minlp_instance.dat", point::Point=Point(), exportprecond=false)
  expos = SortedDict{Exponent, String}()
  precond_cstrs = SortedSet{String}()

  ## First pass on problem for defining all monomials necessary for export
  collect_requiredmonomials!(pb_optim.objective, expos)
  for (ctrname, ctr) in pb_optim.constraints
    collect_requiredmonomials!(ctr.p, expos)
  end

  ## Get max length varname
  maxvarlen = 4
  for var in pb_optim.variables
    occursin(r" ", var[1]) && error(LOGGER, "export_to_dat(): Variable '$var' is not suitable for dat export - no spaces allowed.")
    (length(var[1]) > maxvarlen) && (maxvarlen = length(var[1]))
  end

  for exponame in values(expos)
    occursin(r" ", exponame) && error(LOGGER, "export_to_dat(): Variable '$var' is not suitable for dat export - no spaces allowed.")
    (length(exponame) > maxvarlen) && (maxvarlen = length(exponame))
  end

  ## Get max length dat constraint name
  maxcstrlen = 9
  for cstrname in keys(pb_optim.constraints)
    occursin(r" ", cstrname) && error(LOGGER, "export_to_dat(): Constraint '$cstrname' is not suitable for dat export - no spaces allowed.")
    ((length(cstrname) > maxcstrlen)) && (maxcstrlen = length(cstrname))
  end

  isdir(outpath) || mkpath(outpath)
  filename = joinpath(outpath, filename)
  touch(filename)
  outfile = open(filename, "w")

  ## Comment section
  print_doc(outfile, filename)

  ## Print variables
  variables = get_variables(pb_optim)
  print_variables(outfile, variables, point, maxvarlen, maxcstrlen)

  ## Print objective
  print_poly!(outfile, pb_optim.objective, "OBJ", maxvarlen, maxcstrlen, expos)

  ## Print constraints
  for (cstrname, cstr) in pb_optim.constraints
    cstr = pb_optim.constraints[cstrname]
    print_constraint(outfile, cstrname, cstr, maxvarlen, maxcstrlen, expos)

    if cstr.precond != :none
      push!(precond_cstrs, cstrname)
    end
  end

  ## Print collected monomials definition
  for (expo, exponame) in expos
    for (var, deg) in expo
      print_dat_line(outfile, "MONO_DEF", exponame, var.name, "NONE", deg.explvar, deg.conjvar, maxvarlen, maxcstrlen)
    end
  end
  close(outfile)

  ## Print constraints with preconditionning
  filename = joinpath(outpath, "real_minlp_precond_cstrs.dat")
  touch(filename)
  outfile = open(filename, "w")

  print_string(outfile, "#cstrname", maxcstrlen)
  @printf(outfile, "%10s\n", "Precondtype")
  if exportprecond
    for cstrname in precond_cstrs
      if get_constraint(pb_optim, cstrname).precond == :sqrt
        print_string(outfile, cstrname, maxcstrlen)
        @printf(outfile, "%10s\n", "SQRT")
      else
        warn(LOGGER, "export_dat(): Unknown preconditionning for cstr $cstrname.")
      end
    end
  end
  close(outfile)

  return
end

##
##
function export_matpower_to_dat(QCQP::Problem, filename::String, pt::Point = Point())
  ## Get max length varname
  maxvarlen = -1
  for var in QCQP.variables
    if length(var[1]) > maxvarlen
      maxvarlen = length(var[1])
    end
  end

  ## Sort constraints by type (voltm, unit and rest), and build dat constraint name
  cstr_keys = [key for key in SortedSet(keys(QCQP.constraints))]
  voltm_keys = filter(x->occursin(r"VOLTM", x), cstr_keys)
  unit_keys = filter(x->occursin(r"UNIT", x), cstr_keys)
  load_keys = setdiff(cstr_keys, union(unit_keys, voltm_keys))
  println(load_keys)
  id_to_loadkey = SortedDict(nb_from_str(str)=> (str, "LOAD_$(string(nb_from_str(str)))") for str in load_keys)
  id_to_voltmkey = SortedDict(nb_from_str(str)=> (str, String(split(str, "_")[2])) for str in voltm_keys)
  id_to_unitkey = SortedDict(nb_from_str(str)=> (str, "UNIT_$(string(nb_from_str(str)))") for str in unit_keys)

  ## Get max length dat constraint name
  maxcstrlen = -1
  for (_, val) in union(id_to_loadkey, id_to_voltmkey, id_to_unitkey)
    if length(val[2]) > maxcstrlen
      maxcstrlen = length(val[2])
    end
  end

  touch(filename)
  outfile = open(filename, "w")

  ## Comment section
  print_doc(outfile, filename)

  ## Print variables
  variables = get_variables(QCQP)
  print_variables(outfile, variables, pt, maxvarlen, maxcstrlen)

  ## Print objective
  const_printed = false
  const_val = 0
  for (expo, coeff) in QCQP.objective
      if length(expo) != 0
        print_quad_expo(outfile, expo, "OBJ", coeff, maxvarlen, maxcstrlen)
      else
        const_val = coeff
      end
  end
  print_dat_line(outfile, "CONST", "OBJ", "NONE", "NONE", real(const_val), imag(const_val), maxvarlen, maxcstrlen)

  for (i, loadkey) in id_to_loadkey
    cstr = QCQP.constraints[loadkey[1]]
    print_constraint(outfile, loadkey[2], cstr, maxvarlen, maxcstrlen, expos)
  end
  for (i, voltmkey) in id_to_voltmkey
    cstr = QCQP.constraints[voltmkey[1]]
    print_constraint(outfile, voltmkey[2], cstr, maxvarlen, maxcstrlen, expos)
  end
  for (i, unitkey) in id_to_unitkey
    cstr = QCQP.constraints[unitkey[1]]
    print_constraint(outfile, unitkey[2], cstr, maxvarlen, maxcstrlen, expos)
  end
  close(outfile)
end


text:
#===================================================================================================
  Generic pairwisematrix functions for kernels consuming two vectors
===================================================================================================#

# These two functions may be used to define any kernel
pairwise{T}(κ::PairwiseKernel{T}, x::T, y::T) = phi(κ, x, y)
function unsafe_pairwise{T}(κ::PairwiseKernel{T}, x::AbstractArray{T}, y::AbstractArray{T})
    s = zero(T)
    @inbounds for i in eachindex(x)
        s += phi(κ, x[i], y[i])
    end
    s
end

function pairwise{T}(κ::PairwiseKernel{T}, x::AbstractArray{T}, y::AbstractArray{T})
    length(x) == length(y) || throw(DimensionMismatch("Arrays x and y must have the same length."))
    unsafe_pairwise(κ, x, y)
end

for (order, dimension) in ((:(:row), 1), (:(:col), 2))
    isrowmajor = order == :(:row)
    @eval begin

        @inline function subvector(::Type{Val{$order}}, X::AbstractMatrix,  i::Integer)
            $(isrowmajor ? :(slice(X, i, :)) : :(slice(X, :, i)))
        end

        @inline function init_pairwisematrix{T}(
                 ::Type{Val{$order}},
                X::AbstractMatrix{T}
            )
            Array(T, size(X,$dimension), size(X,$dimension))
        end

        @inline function init_pairwisematrix{T}(
                 ::Type{Val{$order}},
                X::AbstractMatrix{T},
                Y::AbstractMatrix{T}
            )
            Array(T, size(X,$dimension), size(Y,$dimension))
        end

        function checkpairwisedimensions{T}(
                 ::Type{Val{$order}},
                K::Matrix{T},
                X::AbstractMatrix{T}
            )
            n = size(K,1)
            if size(K,2) != n
                throw(DimensionMismatch("Kernel matrix K must be square"))
            elseif size(X, $dimension) != n
                errorstring = string("Dimensions of K must match dimension ", $dimension, "of X")
                throw(DimensionMismatch(errorstring))
            end
            return n
        end

        function checkpairwisedimensions(
                 ::Type{Val{$order}},
                K::Matrix,
                X::AbstractMatrix,
                Y::AbstractMatrix
            )
            n = size(X, $dimension)
            m = size(Y, $dimension)
            if n != size(K,1)
                errorstring = string("Dimension 1 of K must match dimension ", $dimension, "of X")
                throw(DimensionMismatch(errorstring))
            elseif m != size(K,2)
                errorstring = string("Dimension 2 of K must match dimension ", $dimension, "of Y")
                throw(DimensionMismatch(errorstring))
            end
            return (n, m)
        end

        function pairwisematrix!{T}(
                σ::Type{Val{$order}},
                K::Matrix{T},
                κ::PairwiseKernel{T},
                X::AbstractMatrix{T},
                symmetrize::Bool
            )
            n = checkpairwisedimensions(σ, K, X)
            for j = 1:n
                xj = subvector(σ, X, j)
                for i = 1:j
                    xi = subvector(σ, X, i)
                    @inbounds K[i,j] = unsafe_pairwise(κ, xi, xj)
                end
            end
            symmetrize ? LinAlg.copytri!(K, 'U', false) : K
        end

        function pairwisematrix!{T}(
                σ::Type{Val{$order}},
                K::Matrix{T},
                κ::PairwiseKernel{T},
                X::AbstractMatrix{T},
                Y::AbstractMatrix{T},
            )
            n, m = checkpairwisedimensions(σ, K, X, Y)
            for j = 1:m
                yj = subvector(σ, Y, j)
                for i = 1:n
                    xi = subvector(σ, X, i)
                    @inbounds K[i,j] = unsafe_pairwise(κ, xi, yj)
                end
            end
            K
        end
    end
end

#=
@inline function pairwisematrix{T}(
        σ::Union{Type{Val{:row}},Type{Val{:col}}},
        κ::PairwiseKernel{T},
        X::AbstractMatrix{T},
        symmetrize::Bool = true
    )
    pairwisematrix!(σ, init_pairwisematrix(σ, X), κ, X, symmetrize)
end

@inline function pairwisematrix{T}(
        σ::Union{Type{Val{:row}},Type{Val{:col}}},
        κ::PairwiseKernel{T},
        X::AbstractMatrix{T},
        Y::AbstractMatrix{T}
    )
    pairwisematrix!(σ, init_pairwisematrix(σ, X, Y), κ, X, Y)
end
=#

#===================================================================================================
  ScalarProduct and SquaredDistance using BLAS/Built-In methods
===================================================================================================#

function squared_distance!{T<:AbstractFloat}(G::Matrix{T}, xᵀx::Vector{T}, symmetrize::Bool)
    if !((n = length(xᵀx)) == size(G,1) == size(G,2))
        throw(DimensionMismatch("Gramian matrix must be square."))
    end
    @inbounds for j = 1:n, i = (1:j)
        G[i,j] = xᵀx[i] - 2G[i,j] + xᵀx[j]
    end
    symmetrize ? LinAlg.copytri!(G, 'U') : G
end

function squared_distance!{T<:AbstractFloat}(G::Matrix{T}, xᵀx::Vector{T}, yᵀy::Vector{T})
    if size(G,1) != length(xᵀx)
        throw(DimensionMismatch("Length of xᵀx must match rows of G"))
    elseif size(G,2) != length(yᵀy)
        throw(DimensionMismatch("Length of yᵀy must match columns of G"))
    end
    @inbounds for I in CartesianRange(size(G))
        G[I] = xᵀx[I[1]] - 2G[I] + yᵀy[I[2]]
    end
    G
end

for (order, dimension) in ((:(:row), 1), (:(:col), 2))
    isrowmajor = order == :(:row)
    @eval begin

        function dotvectors!{T<:AbstractFloat}(
                 ::Type{Val{$order}},
                xᵀx::Vector{T},
                X::Matrix{T}
            )
            if !(size(X,$dimension) == length(xᵀx))
                errorstring = string("Dimension mismatch on dimension ", $dimension)
                throw(DimensionMismatch(errorstring))
            end
            fill!(xᵀx, zero(T))
            for I in CartesianRange(size(X))
                xᵀx[I.I[$dimension]] += X[I]^2
            end
            xᵀx
        end

        @inline function dotvectors{T<:AbstractFloat}(σ::Type{Val{$order}}, X::Matrix{T})
            dotvectors!(σ, Array(T, size(X,$dimension)), X)
        end

        function gramian!{T<:AbstractFloat}(
                 ::Type{Val{$order}},
                G::Matrix{T},
                X::Matrix{T},
                symmetrize::Bool
            )
            LinAlg.syrk_wrapper!(G, $(isrowmajor ? 'N' : 'T'), X)
            symmetrize ? LinAlg.copytri!(G, 'U') : G
        end

        @inline function gramian!{T<:AbstractFloat}(
                 ::Type{Val{$order}},
                G::Matrix{T},
                X::Matrix{T},
                Y::Matrix{T}
            )
            LinAlg.gemm_wrapper!(G, $(isrowmajor ? 'N' : 'T'), $(isrowmajor ? 'T' : 'N'), X, Y)
        end

        @inline function pairwisematrix!{T}(
                σ::Type{Val{$order}},
                K::Matrix{T},
                κ::ScalarProductKernel{T},
                X::Matrix{T},
                symmetrize::Bool
            )
            gramian!(σ, K, X, symmetrize)
        end

        @inline function pairwisematrix!{T}(
                σ::Type{Val{$order}},
                K::Matrix{T},
                κ::ScalarProductKernel{T},
                X::Matrix{T},
                Y::Matrix{T},
            )
            gramian!(σ, K, X, Y)
        end

        function pairwisematrix!{T}(
                σ::Type{Val{$order}},
                K::Matrix{T},
                κ::SquaredDistanceKernel{T},
                X::Matrix{T},
                symmetrize::Bool
            )
            gramian!(σ, K, X, false)
            xᵀx = dotvectors(σ, X)
            squared_distance!(K, xᵀx, symmetrize)
        end

        function pairwisematrix!{T}(
                σ::Type{Val{$order}},
                K::Matrix{T},
                κ::SquaredDistanceKernel{T},
                X::Matrix{T},
                Y::Matrix{T},
            )
            gramian!(σ, K, X, Y)
            xᵀx = dotvectors(σ, X)
            yᵀy = dotvectors(σ, Y)
            squared_distance!(K, xᵀx, yᵀy)
        end
    end
end


text:
using Photon.Layers
using Photon.Losses: MSELoss
using Photon

function test_performance()

  model = Sequential(
    Conv2D(16, (3,3), :relu),
    Conv2D(16, (3,3), :relu),
    BatchNorm(),
    MaxPool2D(4),
    Conv2D(64, (3,3), :relu),
    Conv2D(64, (3,3), :relu),
    BatchNorm(),
    MaxPool2D(4),
    Conv2D(256, (3,3), :relu),
    Conv2D(256, (3,3), :relu),
    BatchNorm(),
    MaxPool2D(4),
    Dense(64, :relu),
    Dense(10, :relu)
  )

  # Create dummy data
  X = (randn(224,224,3,16) for _=1:100)
  Y = (randn(10,16) for _=1:100)
  data = collect(zip(X,Y))

  # Create the workout
  workout = Workout(model, MSELoss())

  # One run ensure all is compiled
  train!(workout, data)

  # 10 runs we going to measure
  @time train!(workout, data, epochs=10)
end

test_performance()


text:
function wordcount(sentence)

    sentence = lowercase(sentence)
    sentence = replace(sentence, r"[^\w\s,']'?" => " ")
    sentence = replace(sentence, r"\s'|'\s" =>  " ")

    words = split(sentence, r"[\s,]+")
    words = filter(x -> !(x in ["", "'"]), words)

    counter = Dict{String, Int}()
    for word in words
        if word in keys(counter)
            counter[word] += 1
        else
            counter[word] = 1
        end
    end
    return counter
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
function parse_script(script::AbstractString)
    optlines = String[]
    codelines = String[]
    open(script) do file
        while true
            l = readline(file)
            l === "# ```julia" && @goto juliafound
            l === "# ```jlcli" && break
            eof(file) && return
        end
        while true
            l = readline(file)
            l === "# ```" && break
            push!(optlines, l)
            eof(file) && return
        end
        while true
            l = readline(file)
            l === "# ```julia" && break
            eof(file) && return
        end
        @label juliafound
        while true
            l = readline(file)
            l === "# ```" && break
            push!(codelines, l)
            eof(file) && return
        end
    end
    if isempty(codelines)
        return ErrorException("No ```julia ... ``` block in script: $script")
    end
    if !all(startswith(l, "# ") for l in optlines)
        return ErrorException(string(
            "Invalid ```jlcli ... ``` block in script: $script\n",
            "Some lines do not start with `# `."
        ))
    end
    if !all(startswith(l, "# ") for l in codelines)
        return ErrorException(string(
            "Invalid ```julia ... ``` block in script: $script\n",
            "Some lines do not start with `# `."
        ))
    end
    strip2(l) = chop(l, head=2, tail=0)
    map!(strip2, codelines, codelines)
    map!(strip2, optlines, optlines)

    # TODO: use parser
    optsetup = join(optlines, "\n")
    m = Module()
    try
        include_string(m, optsetup)
    catch err
        return err
    end
    if :julia ∈ names(m)
        julia = m.julia
        julia isa AbstractString ||
            return ArgumentError("`julia` option must be a `String`. Got:\n$julia")
    else
        julia = "julia"
    end

    project = realpath(script)
    code = join(codelines, "\n")
    return (julia = julia, project = project, code = code)
end

function compilescript(
    io::IO,
    code,
    manifest::AbstractString;
    julia::Union{AbstractString, Nothing} = nothing,
)
    print(
        io,
        """
        #!/usr/bin/env jlcli
        # This file is created by JuliaCLI.jl
        """,
    )
    if julia !== nothing
        print(
            io,
            """
            #
            # ```jlcli
            # julia = $(repr(julia))
            # ```
            """,
        )
    end
    print(
        io,
        """
        #
        # ```julia
        """
    )
    for line in split(code, "\n")
        println(io, "# ", line)
    end
    print(
        io,
        """
        # ```
        #
        """
    )
    write(io, read(manifest))
end

function compilescript(
    scriptpath::AbstractString,
    code,
    manifest::AbstractString;
    julia::Union{AbstractString, Nothing} = nothing,
    overwrite::Bool = false,
)
    if isfile(scriptpath) && !overwrite
        error("File `$scriptpath` exists.")
    end
    open(scriptpath, write=true) do io
        compilescript(io, code, manifest; julia=julia)
    end
    Base.Filesystem.chmod(scriptpath, 0o777)
end


text:
struct SimWfc3D{T<:AbstractFloat}
  character::Tuple{T,Int64,Symbol,Symbol}
  unit_cell::Tuple{Int64,Int64,Int64}
  center::Point3D{T}
end

function SimWfc3D(n::Int64,character::Symbol,sub_character::Symbol,unit_cell,center::Point3D{T},Z) where T
  if character==:s
    return SimWfc3D{T}((Z,n,character,sub_character),(unit_cell...),center)
  elseif character==:p
    return SimWfc3D{T}((Z,n,character,sub_character),(unit_cell...),center)
  end
end

mutable struct SimModel3D{T<:AbstractFloat}
  cell_origins::Array{Point3D{T},3}
  prim_vec::Matrix{T}
  wfcs::Array{SimWfc3D{T},1}
end

function SimModel3D(supercell,characters::Array{Tuple{Int64,Symbol,Symbol}},centers,Z_s,unit_cells,prim_vec::Matrix{T}) where T
  #we create the supercell + a bounding box so that we don't get boudary effects in the calculations
  tot_cells = supercell+2
  model_origin = -(div(tot_cells-1,2))*Point3D{T}(prim_vec[1,:]+prim_vec[2,:]+prim_vec[3,:])
  #we store the origin of each of the unit cells
  #right now indexing runs from 1->tot_cell, would it be nicer to have -div(tot_cells,2):div(tot_cells,2)?
  cell_origins = [model_origin+Point3D{T}(a*prim_vec[1,:]+b*prim_vec[2,:]+c*prim_vec[3,:]) for a=0:tot_cells-1,b=0:tot_cells-1,c=0:tot_cells-1]
  #convert relative unit_cells to absolute ones
  abs_unit_cells = [cell + div(tot_cells-1,2) for cell in unit_cells]
  abs_centers = [model_origin + center + Point3D{T}(unit_cell[1]*prim_vec[1,:]+unit_cell[2]*prim_vec[2,:]+unit_cell[3]*prim_vec[3,:]) for (center,unit_cell) in zip(centers,abs_unit_cells)]
  wfcs = [SimWfc3D(character[1],character[2],character[3],unit_cell,center,Z) for (character,unit_cell,center,Z) in zip(characters,unit_cells,abs_centers,Z_s)]
  return SimModel3D(cell_origins,prim_vec,wfcs)
end


text:
module MultiAssayExperiments
include("class.jl")
include("miscellaneous.jl")
include("exampleobject.jl")
include("expandsampledata.jl")
include("getters.jl")
include("setters.jl")
include("multifilter.jl")
include("filtersamplemap.jl")
include("dropunused.jl")
end # module


text:
"""
Graph Laplacian of a `m×n` lattice.
"""
function gridlaplacian(T, m, n)
    S = sparse(T(0.0)I, n*m, n*m)
    linear = LinearIndices((1:m, 1:n))
    for i in 1:m
        for j in 1:n
            for (i2, j2) in ((i + 1, j), (i, j + 1))
                if i2 <= m && j2 <= n
                    S[linear[i, j], linear[i2, j2]] -= 1
                    S[linear[i2, j2], linear[i, j]] -= 1

                    S[linear[i, j], linear[i, j]] += 1
                    S[linear[i2, j2], linear[i2, j2]] += 1
                end
            end
        end
    end
    S
end