text:
"""
    Code for convenience functions for calculating CCFs for each chunk of an AbstractChunkListTimeseries (see RvSpectMLBase)
Author: Eric Ford
Created: August 2020
"""

"""  `calc_order_ccf_chunklist_timeseries( chunklist_timeseries, ccf_plan )`
Convenience function to compute separate CCFs for each chunk (potentially an order or view around one or two lines) of each spectrum in a timeseries.
    CCF is evaluated using line list and mask_shape provided by the ccf plan for each chunk.
Uses multiple threads if avaliable.
# Inputs:
- chunklist_timeseries
# Optional Arguments:
- ccf_plan (BasicCCFPlan())
# Return:
A 3-d array containing the CCF at each (velocity, chunk, spectrum)
Note that the ccf_plan provided is used as a template for creating a custom ccf_plan for each chunk that
    only includes lines that reliably appear in that order for all spectra in the chunklist_timeseries.
"""
function calc_order_ccf_chunklist_timeseries(clt::AbstractChunkListTimeseries,
                plan::PlanT = BasicCCFPlan(); verbose::Bool = false ,
                Δfwhm::AbstractVector{T} = zeros(0)
                ) where { PlanT<:AbstractCCFPlan, T<:Real  }

  #@assert issorted( plan.line_list.λ )
  #nvs = length(calc_ccf_v_grid(plan))
  nvs = calc_length_ccf_v_grid(plan)
  norders = length(clt.chunk_list[1].data)
  nobs =  length(clt.chunk_list)
  order_ccfs = zeros(nvs, norders, nobs)
  num_lines = length(plan.line_list.λ)
  plan_for_chunk = Vector{BasicCCFPlan}(undef,num_chunks(clt))
  for chid in 1:num_chunks(clt)
      order = clt[1].order[chid]
      if typeof(plan.line_list) <: BasicLineList2D
          start_order_idx = searchsortedfirst(plan.line_list.order,order)
          stop_order_idx = searchsortedlast(view(plan.line_list.order,start_order_idx:num_lines),order) + start_order_idx-1
      else
          start_order_idx = 1
          stop_order_idx = num_lines
      end
      if (1<=start_order_idx<=stop_order_idx<=num_lines)
          # find the maximum lower wavelength for the chunk, and the minumum upper wavelength, over all observations
          λmin = maximum(map(obsid->first(clt.chunk_list[obsid].data[chid].λ), 1:length(clt) ))
          λmax = minimum(map(obsid->last( clt.chunk_list[obsid].data[chid].λ), 1:length(clt) ))
          # extend λmin/λmax by the velocity range over which we don't want the mask to change
          λmin  = λmin / (calc_doppler_factor(plan.v_center)*calc_doppler_factor(-plan.v_range_no_mask_change))
          λmax  = λmax / (calc_doppler_factor(plan.v_center)*calc_doppler_factor(plan.v_range_no_mask_change))
          # extend λmin/λmax by the mask width
          upper_edge_of_mask_for_line_at_λmin = λ_max(plan.mask_shape,λmin)
          lower_edge_of_mask_for_line_at_λmax = λ_min(plan.mask_shape,λmax)
          # find the first and last mask entries to use in each chunk
          start_line_idx = searchsortedfirst(view(plan.line_list.λ,start_order_idx:stop_order_idx),upper_edge_of_mask_for_line_at_λmin) + start_order_idx-1
          #stop_line_idx  = num_lines+1 - searchsortedfirst(view(plan.line_list.λ,num_lines:-1:1),lower_edge_of_mask_for_line_at_λmax,rev=true)
          stop_line_idx  = start_order_idx-1 + searchsortedlast(view(plan.line_list.λ,start_order_idx:stop_order_idx),lower_edge_of_mask_for_line_at_λmax)
          if verbose
              flush(stdout)
              println("extrema(plan.line_list.λ) = ",extrema(plan.line_list.λ) )
              if (1 <= start_line_idx <= length(plan.line_list.λ)) && (1 <= stop_line_idx <= length(plan.line_list.λ))
                  println("start_order_idx = ", start_order_idx, "  stop_order_idx = ", stop_order_idx)
                  println("start_line_idx = ", start_line_idx, " λ= ", plan.line_list.λ[start_line_idx]) #, " order= ", plan.line_list.order[start_line_idx])
                  println("stop_line_idx = ", stop_line_idx, " λ= ", plan.line_list.λ[stop_line_idx]) #, " order= ", plan.line_list.order[stop_line_idx])
              end
              println("upper_edge_of_mask_for_line_at_λmin = ", upper_edge_of_mask_for_line_at_λmin, " Lower_edge_of_mask_for_line_at_λmax = ", lower_edge_of_mask_for_line_at_λmax)
          end

          if (1 <= start_line_idx <= length(plan.line_list.λ)) && (1 <= stop_line_idx <= length(plan.line_list.λ))
              if verbose && typeof(plan.line_list) <: BasicLineList2D
                  println("start_line_idx = ", start_line_idx, " λ= ", plan.line_list.λ[start_line_idx], " order= ", plan.line_list.order[start_line_idx])
                  println("stop_line_idx = ", stop_line_idx, " λ= ", plan.line_list.λ[stop_line_idx], " order= ", plan.line_list.order[stop_line_idx])
                  println(" Using ", length(start_line_idx:stop_line_idx), " lines for chunk ", chid, " order ", order)
              end
              if ! issorted( view(plan.line_list.λ, start_line_idx:stop_line_idx ) )
                  println("# View is not sorted.  Why?")
                  for i in start_line_idx:stop_line_idx
                      println("i= ", i, " λ= ", plan.line_list.λ[i], " order= ", plan.line_list.order[i])
                  end
              end
              @assert issorted( view(plan.line_list.λ, start_line_idx:stop_line_idx ) )
              # line_list_for_chunk = BasicLineList2D(view(plan.line_list.λ,start_line_idx:stop_line_idx), view(plan.line_list.weight,start_line_idx:stop_line_idx), fill(order,length(start_line_idx:stop_line_idx)) )
              line_list_for_chunk = BasicLineList(view(plan.line_list.λ,start_line_idx:stop_line_idx), view(plan.line_list.weight,start_line_idx:stop_line_idx)  )
          else  # No lines in this chunk!
              line_list_for_chunk = EmptyBasicLineList()
          end
      else
        line_list_for_chunk = EmptyBasicLineList()
      end

      #create a plan for this chunk that only includes the mask entries we want for this chunk
      plan_for_chunk[chid] = BasicCCFPlan( line_list=line_list_for_chunk, midpoint=plan.v_center, step=plan.v_step, max=plan.v_max, mask_shape=plan.mask_shape )
  end
  #@threaded mapreduce(obsid->calc_ccf_chunklist(clt.chunk_list[obsid], plan_for_chunk),hcat, 1:length(clt) )
  #Threads.@threads for i in 1:nobs
    #  order_ccfs[:,:,i] .= calc_order_ccfs_chunklist(clt.chunk_list[i], plan)
  #end
  Threads.@threads for i in 1:nobs
      this_Δfwhm = length(Δfwhm) == nobs ? Δfwhm[i] : 0.0
      order_ccfs[:,:,i] .= calc_order_ccfs_chunklist(clt.chunk_list[i], plan_for_chunk, Δfwhm=this_Δfwhm, assume_sorted=true )
  end

  return order_ccfs
end

"""  `calc_order_ccf_and_var_chunklist_timeseries( chunklist_timeseries, ccf_plan )`
Convenience function to compute separate CCFs for each chunk (potentially an order or view around one or two lines) of each spectrum in a timeseries.
    CCF is evaluated using line list and mask_shape provided by the ccf plan for each chunk.
Uses multiple threads if avaliable.
# Inputs:
- chunklist_timeseries
# Optional Arguments:
- ccf_plan (BasicCCFPlan())
# Return:
A 3-d array containing the CCF at each (velocity, chunk, spectrum)
Note that the ccf_plan provided is used as a template for creating a custom ccf_plan for each chunk that
    only includes lines that reliably appear in that order for all spectra in the chunklist_timeseries.
"""
function calc_order_ccf_and_var_chunklist_timeseries(clt::AbstractChunkListTimeseries,
                plan::PlanT = BasicCCFPlan(); verbose::Bool = false,
                ccf_var_scale::Real = 1.0, Δfwhm::AbstractVector{T} = zeros(0)
                 ) where {
                    PlanT<:AbstractCCFPlan, T<:Real }
  #@assert issorted( plan.line_list.λ )
  #nvs = length(calc_ccf_v_grid(plan))
  nvs = calc_length_ccf_v_grid(plan)
  norders = length(clt.chunk_list[1].data)
  nobs =  length(clt.chunk_list)
  order_ccfs = zeros(nvs, norders, nobs)
  order_ccf_vars = zeros(nvs, norders, nobs)
  num_lines = length(plan.line_list.λ)
  plan_for_chunk = Vector{BasicCCFPlan}(undef,num_chunks(clt))
  for chid in 1:num_chunks(clt)
      order = clt[1].order[chid]
      if typeof(plan.line_list) <: BasicLineList2D
          start_order_idx = searchsortedfirst(plan.line_list.order,order)
          stop_order_idx = searchsortedlast(view(plan.line_list.order,start_order_idx:num_lines),order) + start_order_idx-1
      else
          start_order_idx = 1
          stop_order_idx = num_lines
      end
      if (1<=start_order_idx<=stop_order_idx<=num_lines)
          # find the maximum lower wavelength for the chunk, and the minumum upper wavelength, over all observations
          λmin = maximum(map(obsid->first(clt.chunk_list[obsid].data[chid].λ), 1:length(clt) ))
          λmax = minimum(map(obsid->last( clt.chunk_list[obsid].data[chid].λ), 1:length(clt) ))
          # extend λmin/λmax by the velocity range over which we don't want the mask to change
          λmin  = λmin / (calc_doppler_factor(plan.v_center)*calc_doppler_factor(-plan.v_range_no_mask_change))
          λmax  = λmax / (calc_doppler_factor(plan.v_center)*calc_doppler_factor(plan.v_range_no_mask_change))
          # extend λmin/λmax by the mask width
          upper_edge_of_mask_for_line_at_λmin = λ_max(plan.mask_shape,λmin)
          lower_edge_of_mask_for_line_at_λmax = λ_min(plan.mask_shape,λmax)
          # find the first and last mask entries to use in each chunk
          start_line_idx = searchsortedfirst(view(plan.line_list.λ,start_order_idx:stop_order_idx),upper_edge_of_mask_for_line_at_λmin) + start_order_idx-1
          stop_line_idx  = start_order_idx-1 + searchsortedlast(view(plan.line_list.λ,start_order_idx:stop_order_idx),lower_edge_of_mask_for_line_at_λmax)
          if verbose
              flush(stdout)
              println("extrema(plan.line_list.λ) = ",extrema(plan.line_list.λ) )
              if (1 <= start_line_idx <= length(plan.line_list.λ)) && (1 <= stop_line_idx <= length(plan.line_list.λ))
                  println("start_order_idx = ", start_order_idx, "  stop_order_idx = ", stop_order_idx)
                  println("start_line_idx = ", start_line_idx, " λ= ", plan.line_list.λ[start_line_idx]) #, " order= ", plan.line_list.order[start_line_idx])
                  println("stop_line_idx = ", stop_line_idx, " λ= ", plan.line_list.λ[stop_line_idx]) #, " order= ", plan.line_list.order[stop_line_idx])
              end
              println("upper_edge_of_mask_for_line_at_λmin = ", upper_edge_of_mask_for_line_at_λmin, " Lower_edge_of_mask_for_line_at_λmax = ", lower_edge_of_mask_for_line_at_λmax)
          end

          if (1 <= start_line_idx <= length(plan.line_list.λ)) && (1 <= stop_line_idx <= length(plan.line_list.λ))
              if verbose && typeof(plan.line_list) <: BasicLineList2D
                  println("start_line_idx = ", start_line_idx, " λ= ", plan.line_list.λ[start_line_idx], " order= ", plan.line_list.order[start_line_idx])
                  println("stop_line_idx = ", stop_line_idx, " λ= ", plan.line_list.λ[stop_line_idx], " order= ", plan.line_list.order[stop_line_idx])
                  println(" Using ", length(start_line_idx:stop_line_idx), " lines for chunk ", chid, " order ", order)
              end
              if ! issorted( view(plan.line_list.λ, start_line_idx:stop_line_idx ) )
                  println("# View is not sorted.  Why?")
                  for i in start_line_idx:stop_line_idx
                      println("i= ", i, " λ= ", plan.line_list.λ[i], " order= ", plan.line_list.order[i])
                  end
              end
              @assert issorted( view(plan.line_list.λ, start_line_idx:stop_line_idx ) )
              # line_list_for_chunk = BasicLineList(view(plan.line_list.λ,start_line_idx:stop_line_idx), view(plan.line_list.weight,start_line_idx:stop_line_idx) )
              line_list_for_chunk = BasicLineList(view(plan.line_list.λ,start_line_idx:stop_line_idx), view(plan.line_list.weight,start_line_idx:stop_line_idx)  )
          else  # No lines in this chunk!
              line_list_for_chunk = EmptyBasicLineList()
          end
      else
        line_list_for_chunk = EmptyBasicLineList()
      end
      #create a plan for this chunk that only includes the mask entries we want for this chunk
      plan_for_chunk[chid] = BasicCCFPlan( line_list=line_list_for_chunk, midpoint=plan.v_center, step=plan.v_step, max=plan.v_max, mask_shape=plan.mask_shape )
  end
  #@threaded mapreduce(obsid->calc_ccf_chunklist(clt.chunk_list[obsid], plan_for_chunk),hcat, 1:length(clt) )
  #Threads.@threads for i in 1:nobs
    #  order_ccfs[:,:,i] .= calc_order_ccfs_chunklist(clt.chunk_list[i], plan)
  #end
  Threads.@threads for i in 1:nobs
      this_Δfwhm = length(Δfwhm) == nobs ? Δfwhm[i] : 0.0
      for j in 1:num_chunks(clt)
          ( ccf_tmp, ccf_var_tmp ) = calc_ccf_and_var_chunk(clt.chunk_list[i][j], plan_for_chunk[j],
                                            ccf_var_scale=ccf_var_scale, Δfwhm=this_Δfwhm, assume_sorted=true )
          order_ccfs[:,j,i] .= ccf_tmp
          order_ccf_vars[:,j,i] .= ccf_var_tmp
      end
  end

  return (ccf=order_ccfs, ccf_var=order_ccf_vars)
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
using DifferentialEquations, LinearAlgebra, Plots; pyplot()

k, b, M = 1.2, 0.3, 2.0
A = [0 1;
    -k/M -b/M]

initX = [8., 0.0]
tEnd = 50.0
tRange = 0:0.1:tEnd

manualSol = [exp(A*t)*initX for t in tRange]

linearRHS(x,Amat,t) = Amat*x
prob = ODEProblem(linearRHS, initX, (0,tEnd), A)
sol = solve(prob)

p1 = plot(first.(manualSol), last.(manualSol),
	c=:blue, label="Manual trajectory")
p1 = scatter!(first.(sol.u), last.(sol.u),
	c=:red, ms = 5, msw=0, label="DiffEq package")
p1 = scatter!([initX[1]], [initX[2]],
	c=:black, ms=10, label="Initial state",	xlims=(-7,9), ylims=(-9,7),
	ratio=:equal, xlabel="Displacement", ylabel="Velocity")
p2 = plot(tRange, first.(manualSol),
	c=:blue, label="Manual trajectory")
p2 = scatter!(sol.t, first.(sol.u),
	c=:red, ms = 5, msw=0, label="DiffEq package")
p2 = scatter!([0], [initX[1]],
	c=:black, ms=10, label="Initial state", xlabel="Time",
	ylabel="Displacement")
plot(p1, p2, size=(800,400), legend=:topright)

ext: jl

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
import ..Ipopt

DEFAULT_OPTIMIZER = Ipopt.Optimizer

function Ipopt.Optimizer(m::Model)
    _obj,_grad!,_con!,_jac!,jac_sparsity,_hess!,hess_sparsity = nlp_evaluator(m.obj,m.con)
    _jac_fill_sparsity! = @inline (I,J)->fill_sparsity!(I,J,jac_sparsity)
    _hess_fill_sparsity! = @inline (I,J)->fill_sparsity!(I,J,hess_sparsity)
    p = m.p

    prob =  Ipopt.createProblem(
        m.n,m.xl,m.xu,m.m,m.gl,m.gu,length(jac_sparsity),length(hess_sparsity),
        x->_obj(x,p), (x,g)->_con!(g,x,p), (x,f)->_grad!(f,x,p),
        (x, mode, I,J, jac)-> mode == :Structure ? _jac_fill_sparsity!(I,J) : _jac!(jac,x,p),
        (x, mode, I,J, sig, lag, hess)-> mode == :Structure ? _hess_fill_sparsity!(I,J) : _hess!(hess,x,p,lag,sig)
    )

    prob.x .= m.x

    for (name,option) in m.opt
        Ipopt.addOption(prob,string(name),option)
    end

    m.x = prob.x
    m.l = prob.mult_g
    m.zl= prob.mult_x_L
    m.zu= prob.mult_x_U

    return prob
end

optimize!(prob::Ipopt.IpoptProblem) = Ipopt.solveProblem(prob)


text:
## ---------------------------------------------------
function create_pkg_version(pkgdir::AbstractString=pwd();
        new_version::String="",
        up_major = false,
        up_minor = false,
        up_patch = false,
        registry::String=""
    )

    # find project
    !isdir(pkgdir) && error("Package directory '$pkgdir' not found")
    projfile = Base.current_project(pkgdir)
    (isnothing(projfile) || isempty(dirname(projfile)) || !isfile(projfile)) && error("Project file not found")
    projfile = abspath(projfile)
    pkgdir = dirname(projfile)

    # check Manifest
    unreg_pkgs = sort!(_find_unregistered_pkgs(pkgdir))
    !isempty(unreg_pkgs) && error("Some pkgs are in an 'unregistered' mode, pkgs: $(join(unreg_pkgs, ", "))")

    # up version
    projdict = TOML.parsefile(projfile)
    pkg_version = VersionNumber(projdict["version"])
    pkg_name = projdict["name"]
    pkg_uuid = projdict["uuid"]

    if !isempty(new_version)
        new_version = replace(new_version, "v"=>"")
        new_version = VersionNumber(new_version)
        (pkg_version > new_version) && error("New version $(new_version) is 'older' than current $(pkg_version)")
    elseif up_major
        new_version = VersionNumber(pkg_version.major + 1, 0, 0)
    elseif up_minor
        new_version = VersionNumber(pkg_version.major, pkg_version.minor + 1, 0)
    elseif up_patch
        new_version = VersionNumber(pkg_version.major, pkg_version.minor, pkg_version.patch + 1)
    else
        error("Not new version specify")
    end
    (new_version == pkg_version) && error("Equal new and old versions")

    _info("Package current status"; pkg_name, pkg_uuid, pkg_version)

    # write back Project
    projdict["version"] = string(new_version)
    _save_project(projfile, projdict)

    _info("Package new status"; pkg_name, pkg_uuid, new_version)

    # commit new project
    _info("Adding new Project.toml")
    run(Cmd(Cmd(["git", "add", projfile]); dir=pkgdir); wait=true)
    _info("Committing new Project.toml")
    run(_Cmd(["git", "commit", "-m", "up to $(new_version)"]; dir=pkgdir); wait=true)
    run(_Cmd(["git", "push"]; dir=pkgdir); wait=true)

    # tag and push
    _info("Tagging")
    tag = string("v", new_version)
    run(_Cmd(["git", "tag", tag]; dir=pkgdir); wait=true)
    run(_Cmd(["git", "push", "origin", tag]; dir=pkgdir); wait=true)

    # update registry
    if registry != "0"
        _info("Update registry"; registry, pkg_name, new_version)
        commit_to_registry(pkgdir;
            registry = registry == "1" ? "" : registry,
            verbose = false,
            push = true
        )
    else
        _warn("Ignoring registering"; pkg_name, new_version)
    end

    return pkgdir
end

## ---------------------------------------------------
function run_create_pkg_version(argv::Vector=ARGS)

    ## ---------------------------------------------------------
    argset = ArgParse.ArgParseSettings()
    ArgParse.@add_arg_table! argset begin
        "--pkgdir", "-d"
            help = "the package dir"
            arg_type = String
            default = pwd()
        "--version", "-v"
            help = "the new version to push"
            arg_type = String
            default = ""
        "--registry", "-r"
            help = "-r=Name specify the register to push. Use -r0 to ignore registering and r1 to use installed"
            arg_type = String
            default = "0"
        "--up-major", "-M"
            help = "Will push a new version with the major incremented by 1"
            action = :store_true
        "--up-minor", "-m"
            help = "Will push a new version with the minor incremented by 1"
            action = :store_true
        "--up-patch", "-p"
            help = "Will push a new version with the patch incremented by 1"
            action = :store_true
    end

    parsed_args = ArgParse.parse_args(argv, argset)
    new_version = parsed_args["version"]
    pkgdir = parsed_args["pkgdir"]
    registry = parsed_args["registry"]
    up_major = parsed_args["up-major"]
    up_minor = parsed_args["up-minor"]
    up_patch = parsed_args["up-patch"]

    ## ---------------------------------------------------------
    create_pkg_version(pkgdir;
        new_version, up_major, up_minor, up_patch, registry
    )

end


text:
@testset "model.jl" begin
    @testset "convenience functions" begin
        model = gdemo_default

        # sample from model and extract variables
        vi = VarInfo(model)
        s = vi[@varname(s)]
        m = vi[@varname(m)]

        # extract log pdf of variable object
        lp = getlogp(vi)

        # log prior probability
        lprior = logprior(model, vi)
        @test lprior ≈ logpdf(InverseGamma(2, 3), s) + logpdf(Normal(0, sqrt(s)), m)

        # log likelihood
        llikelihood = loglikelihood(model, vi)
        @test llikelihood ≈ loglikelihood(Normal(m, sqrt(s)), [1.5, 2.0])

        # log joint probability
        ljoint = logjoint(model, vi)
        @test ljoint ≈ lprior + llikelihood
        @test ljoint ≈ lp
    end

    @testset "rng" begin
        model = gdemo_default

        for sampler in (SampleFromPrior(), SampleFromUniform())
            for i in 1:10
                Random.seed!(100 + i)
                vi = VarInfo()
                model(Random.GLOBAL_RNG, vi, sampler)
                vals = DynamicPPL.getall(vi)

                Random.seed!(100 + i)
                vi = VarInfo()
                model(Random.GLOBAL_RNG, vi, sampler)
                @test DynamicPPL.getall(vi) == vals
            end
        end
    end

    @testset "defaults without VarInfo, Sampler, and Context" begin
        model = gdemo_default

        Random.seed!(100)
        s, m = model()

        Random.seed!(100)
        @test model(Random.GLOBAL_RNG) == (s, m)
    end

    @testset "nameof" begin
        @model function test1(x)
            m ~ Normal(0, 1)
            x ~ Normal(m, 1)
        end
        @model test2(x) = begin
            m ~ Normal(0, 1)
            x ~ Normal(m, 1)
        end

        @test nameof(test1(rand())) == :test1
        @test nameof(test2(rand())) == :test2
    end
end


text:
function prime_factors(n)
    found = false

    v = 2

    i = 2

    while i < n - 1 && !found
        if mod(n, i) == 0
            found = true

            v = i
        end

        i += 1
    end

    if found
        factors = vcat([v], prime_factors(div(n, v)))
    else
        factors = [n]
    end

    factors
end

println(maximum(prime_factors(600851475143)))

max_stars_repo_name: seaneshbaugh/rosetta-euler

max_issues_repo_name: seaneshbaugh/rosetta-euler

max_forks_repo_name: seaneshbaugh/rosetta-euler


text:
"""
    $(TYPEDEF)
GeneralBigDataEntry is a generic multipurpose data entry that creates a unique
reproducible key using userId_robotId_sessionId_variableId_key.
"""
mutable struct GeneralBigDataEntry <: AbstractBigDataEntry
    key::Symbol
    storeKey::Symbol # Could swap this to string, but using it as an index later, so better as a symbol I believe.
    createdTimestamp::DateTime
    lastUpdatedTimestamp::DateTime
    mimeType::String
end

"""
    $(SIGNATURES)
Internal function to generate a unique key for the entry - userId_robotId_sessionId_variable_key.
Simple symbol.
"""
function _uniqueKey(dfg::G, v::V, key::Symbol)::Symbol where {G <: AbstractDFG, V <: AbstractDFGVariable}
    key = join(String.([dfg.userId, dfg.robotId, dfg.sessionId, getLabel(v), String(key)]), "_")
    return Symbol(key)
end

GeneralBigDataEntry(key::Symbol, storeKey::Symbol;
                    mimeType::String="application/octet-stream") =
                    GeneralBigDataEntry(key, storeKey, now(), now(), mimeType)

function GeneralBigDataEntry(dfg::G, var::V, key::Symbol;
                             mimeType::String="application/octet-stream") where {G <: AbstractDFG, V <: AbstractDFGVariable}
    return GeneralBigDataEntry(key, _uniqueKey(dfg, var, key), mimeType=mimeType)
end

"""
    $(TYPEDEF)
BigDataEntry in MongoDB.
"""
struct MongodbBigDataEntry <: AbstractBigDataEntry
    key::Symbol
    oid::NTuple{12, UInt8} #mongodb object id
    #maybe other fields such as:
    #flags::Bool ready, valid, locked, permissions
    #MIMEType::Symbol
end

"""
    $(TYPEDEF)
BigDataEntry in a file.
"""
struct FileBigDataEntry <: AbstractBigDataEntry
    key::Symbol
    filename::String
end


text:
"""
A wrapper of `TransientFEOperator` that transforms it to `ODEOperator`, i.e.,
takes A(t,uh,∂tuh,∂t^2uh,...,∂t^Nuh,vh) and returns A(t,uF,∂tuF,...,∂t^NuF)
where uF,∂tuF,...,∂t^NuF represent the free values of the `EvaluationFunction`
uh,∂tuh,∂t^2uh,...,∂t^Nuh.
"""
struct ODEOpFromFEOp{C} <: ODEOperator{C}
  feop::TransientFEOperator{C}
end

get_order(op::ODEOpFromFEOp) = get_order(op.feop)

function allocate_cache(op::ODEOpFromFEOp)
  Ut = get_trial(op.feop)
  U = allocate_trial_space(Ut)
  Uts = (Ut,)
  Us = (U,)
  for i in 1:get_order(op)
    Uts = (Uts...,∂t(Uts[i]))
    Us = (Us...,allocate_trial_space(Uts[i+1]))
  end
  fecache = allocate_cache(op.feop)
  ode_cache = (Us,Uts,fecache)
  ode_cache
end

function allocate_cache(op::ODEOpFromFEOp,v::AbstractVector,a::AbstractVector)
  ode_cache = allocate_cache(op)
  (v,a, ode_cache)
end

function update_cache!(ode_cache,op::ODEOpFromFEOp,t::Real)
  _Us,Uts,fecache = ode_cache
  Us = ()
  for i in 1:get_order(op)+1
    Us = (Us...,evaluate!(_Us[i],Uts[i],t))
  end
  fecache = update_cache!(fecache,op.feop,t)
  (Us,Uts,fecache)
end

function allocate_residual(op::ODEOpFromFEOp,uhF::AbstractVector,ode_cache)
  Us,Uts,fecache = ode_cache
  uh = EvaluationFunction(Us[1],uhF)
  allocate_residual(op.feop,uh,fecache)
end

function allocate_jacobian(op::ODEOpFromFEOp,uhF::AbstractVector,ode_cache)
  Us,Uts,fecache = ode_cache
  uh = EvaluationFunction(Us[1],uhF)
  allocate_jacobian(op.feop,uh,fecache)
end

"""
It provides A(t,uh,∂tuh,...,∂t^Nuh) for a given (t,uh,∂tuh,...,∂t^Nuh)
"""
function residual!(
  b::AbstractVector,
  op::ODEOpFromFEOp,
  t::Real,
  xhF::Tuple{Vararg{AbstractVector}},
  ode_cache)
  Xh, = ode_cache
  dxh = ()
  for i in 2:get_order(op)+1
    dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i]))
  end
  xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh)
  residual!(b,op.feop,t,xh,ode_cache)
end

"""
It adds contribution to the Jacobian with respect to the i-th time derivative,
with i=0,...,N. That is, adding γ_i*[∂A/∂(∂t^iuh)](t,uh,∂tuh,...,∂t^Nuh) for a
given (t,uh,∂tuh,...,∂t^Nuh) to a given matrix J, where γ_i is a scaling coefficient
provided by the `ODESolver`, e.g., 1/Δt for Backward Euler; It represents
∂(δt^i(uh))/∂(uh), in which δt^i(⋅) is the approximation of ∂t^i(⋅) in the solver.
Note that for i=0, γ_i=1.0.
"""
function jacobian!(
  A::AbstractMatrix,
  op::ODEOpFromFEOp,
  t::Real,
  xhF::Tuple{Vararg{AbstractVector}},
  i::Integer,
  γᵢ::Real,
  ode_cache)
  Xh, = ode_cache
  dxh = ()
  for i in 2:get_order(op)+1
    dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i]))
  end
  xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh)
  jacobian!(A,op.feop,t,xh,i,γᵢ,ode_cache)
end

"""
Add the contribution of all jacobians ,i.e., ∑ᵢ γ_i*[∂A/∂(∂t^iuh)](t,uh,∂tuh,...,∂t^Nuh,vh)
"""
function jacobians!(
  J::AbstractMatrix,
  op::ODEOpFromFEOp,
  t::Real,
  xhF::Tuple{Vararg{AbstractVector}},
  γ::Tuple{Vararg{Real}},
  ode_cache)
  Xh, = ode_cache
  dxh = ()
  for i in 2:get_order(op)+1
    dxh = (dxh...,EvaluationFunction(Xh[i],xhF[i]))
  end
  xh=TransientCellField(EvaluationFunction(Xh[1],xhF[1]),dxh)
  jacobians!(J,op.feop,t,xh,γ,ode_cache)
end


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

"""
    ArrayShape{T,N} <: AbstractValueShape

Describes the shape of `N`-dimensional arrays of type `T` and a given size.

Constructor:

    ArrayShape{T}(dims::NTuple{N,Integer}) where {T,N}
    ArrayShape{T}(dims::Integer...) where {T}

e.g.

    shape = ArrayShape{Real}(2, 3)

See also the documentation of [`AbstractValueShape`](@ref).
"""
struct ArrayShape{T,N} <: AbstractValueShape
    dims::NTuple{N,Int}
end

export ArrayShape

ArrayShape{T}(dims::NTuple{N,Integer}) where {T,N} = ArrayShape{T,N}(map(Int, dims))
ArrayShape{T}(dims::Integer...) where {T} = ArrayShape{T}(dims)

@inline Base.size(shape::ArrayShape) = shape.dims
Base.length(shape::ArrayShape) = prod(size(shape))

@inline _valshapeoftype(T::Type{<:AbstractArray}) = throw(ArgumentError("Type $T does not have a fixed shape"))

@inline default_unshaped_eltype(shape::ArrayShape{T}) where {T} =
    default_unshaped_eltype(_valshapeoftype(T))

@inline shaped_type(shape::ArrayShape{T,N}, ::Type{U}) where {T,N,U<:Real} =
    Array{shaped_type(_valshapeoftype(T),U),N}

@inline function valshape(x::AbstractArray{T}) where T
    _valshapeoftype(T) # ensure T has a fixed shape
    ArrayShape{T}(size(x))
end

# Possible extension: valshape(x::AbstractArrayOfSimilarArrays{...})

totalndof(shape::ArrayShape{T}) where{T} =
    prod(size(shape)) * totalndof(_valshapeoftype(T))

(shape::ArrayShape{T,N})(::UndefInitializer) where {T,N} = Array{default_datatype(T),N}(undef, size(shape)...)

@static if VERSION < v"1.3"
    # Workaround for Julia issue #14919
    @inline (shape::ArrayShape)(data::AbstractVector{<:Real}) =
        _apply_shape_to_data(shape, data)
end

const ArrayAccessor{T,N} = ValueAccessor{ArrayShape{T,N}} where {T,N}

Base.@propagate_inbounds vs_getindex(data::AbstractVector{<:Real}, va::ArrayAccessor) = copy(view(data, va))

Base.@propagate_inbounds vs_unsafe_view(data::AbstractVector{<:Real}, va::ArrayAccessor{T,1}) where {T} =
    Base.unsafe_view(data, view_idxs(axes(data, 1), va))

Base.@propagate_inbounds vs_unsafe_view(data::AbstractVector{<:Real}, va::ArrayAccessor{T,N}) where {T,N} =
    reshape(Base.unsafe_view(data, view_idxs(axes(data, 1), va)), size(va.shape)...)

Base.@propagate_inbounds vs_setindex!(data::AbstractVector{<:Real}, v, va::ArrayAccessor) where {T} =
    setindex!(data, v, view_idxs(axes(data, 1), va))

Base.@propagate_inbounds function _bcasted_view(data::AbstractVectorOfSimilarVectors{<:Real}, va::ArrayAccessor{T,1}) where {T,N}
    flat_data = flatview(data)
    idxs = view_idxs(axes(flat_data, 1), va)
    fpview = view(flat_data, idxs, :)
    VectorOfSimilarVectors(fpview)
end

Base.@propagate_inbounds function _bcasted_view(data::AbstractVectorOfSimilarVectors{<:Real}, va::ArrayAccessor{T,N}) where {T,N}
    flat_data = flatview(data)
    idxs = view_idxs(axes(flat_data, 1), va)
    fpview = view(flat_data, idxs, :)
    VectorOfSimilarArrays(reshape(fpview, size(va.shape)..., :))
end

Base.copy(instance::VSBroadcasted2{typeof(getindex),AbstractVectorOfSimilarVectors{<:Real},Ref{<:ArrayAccessor}}) =
    copy(_bcasted_view(instance.args[1], instance.args[2][]))

Base.copy(instance::VSBroadcasted2{typeof(view),AbstractVectorOfSimilarVectors{<:Real},Ref{<:ArrayAccessor}}) =
    _bcasted_view(instance.args[1], instance.args[2][])

function _bcasted_view_unchanged(data::AbstractArray{<:AbstractVector{T}}, shape::ArrayShape{U,1}) where {T<:Real,U>:T}
    _checkcompat_inner(shape, data)
    data
end

# Specialize (::ArrayShape{T,1}).(::AbstractArray{<:AbstractVector{<:Real}}):
Base.copy(instance::VSBroadcasted1{ArrayShape{T,1},AbstractVector{<:AbstractVector{<:Real}}}) where T =
    _bcasted_view_unchanged(instance.args[1], instance.f)
Base.copy(instance::VSBroadcasted1{ArrayShape{T,1},AbstractArray{<:AbstractVector{<:Real}}}) where T =
    _bcasted_view_unchanged(instance.args[1], instance.f)

@inline _bcasted_unshaped(A::AbstractArrayOfSimilarArrays{<:Real,1}) = A

# Specialize unshaped.(::AbstractArray{<:AbstractVector{<:Real}}):
Base.copy(instance::VSBroadcasted1{typeof(unshaped),AbstractVector{<:AbstractVector{<:Real}}}) =
    _bcasted_unshaped(instance.args[1])
Base.copy(instance::VSBroadcasted1{typeof(unshaped),AbstractArray{<:AbstractVector{<:Real}}}) =
    _bcasted_unshaped(instance.args[1])

# TODO: Add support for StaticArray.

# Possible extension: variable/flexible array shapes?


text:
include("tunedpw.jl")
include("tunecbts.jl")


text:
module Databases
    import Base
    import ..Messages

    struct Database
        frame_id_index::Dict{UInt32, Ref{Messages.Message}}
        name_index::Dict{String, Ref{Messages.Message}}

        function Database(messages::Set{Messages.Message})
            v = [messages...]
            e = enumerate(v)
            l1 = [Messages.name(m1) == Messages.name(m2) for (i, m1) in e for (j, m2) in e if i < j]
            l2 = [Messages.frame_id(m1) == Messages.frame_id(m2) for (i, m1) in e for (j, m2) in e if i < j]
            a1 = any(l1)
            a2 = any(l2)

            if a1
                throw(DomainError(a1, "messages with the same name"))
            end

            if a2
                throw(DomainError(a2, "messages with the same frame_id"))
            end

            frame_id_index = Dict{UInt32, Ref{Messages.Message}}()
            name_index = Dict{String, Ref{Messages.Message}}()

            for message in messages
                m = Ref(message)
                frame_id_index[Messages.frame_id(message)] = m
                name_index[Messages.name(message)] = m
            end

            new(frame_id_index, name_index)
        end
    end

    function Database(messages::Messages.Message...)
        s = Set(messages)
        return Database(s)
    end

    function Base.getindex(db::Database, index::String)
        m_ref = db.name_index[index]
        return m_ref[]
    end

    function Base.getindex(db::Database, index::UInt32)
        m_ref = db.frame_id_index[index]
        return m_ref[]
    end

    function Base.getindex(db::Database, index::Integer)
        index = convert(UInt32, index)
        return db[index]
    end

    function Base.get(db::Database, key::String, default=nothing)
        try
            value = db[key]
            return value
        catch
            return default
        end
    end

    function Base.get(db::Database, key::UInt32, default=nothing)
        try
            value = db[key]
            return value
        catch
            return default
        end
    end

    function Base.get(db::Database, key::Integer, default=nothing)
        key = convert(UInt32, key)
        return get(db, key, default)
    end

    export Database
end


text:
# Michelle Medina
# July 7, 2020
# Functions for Hierarchical Bases

export basis, maxlevel, evaluate, H_2_Nodal, Nodal_2_H, calcNodal, Nodal_2_H_new

#************  INDEX  ********************************
"""
Function gives the highest level (l-1 because we are considering level 0)
"""
function maxlevel(H::HierarchicalBasis{T})::Int where {T}
	return length(H.levels) - 1
end

"""
Function gives the highest level (it's log2 because length must be a multiple of 2 + 1)
"""
function maxlevel(N::NodalBasis{T})::Int where {T}
	return log2(length(N.values) - 1)
end

"""
Function getindex now works with Hierarchical Bases
- converts the physical indices to sequential indices on the grid and returns the coefficients
"""
function Base.getindex(H::HierarchicalBasis{T}, l::Int, j::Int)::T where {T}
	@assert 0 <= l <= maxlevel(H)
	@assert 0 <= j <= (2^l)
	return H.levels[l+1].coefficients[j+1]
end

function Base.getindex(N::NodalBasis{T}, j::Int)::T where {T}
	@assert 0 <= j <= (2^maxlevel(N))
	return N.values[j+1]
end

function Base.setindex!(H::HierarchicalBasis{T}, X::T, l::Int, j::Int)::T where {T}
	@assert 0 <= l <= maxlevel(H)
	@assert 0 <= j <= (2^l)
	H.levels[l+1].coefficients[j+1] = X
end

function Base.setindex!(N::NodalBasis{T}, X::T, j::Int)::T where {T}
	@assert 0 <= j <= (2^maxlevel(N))
	N.values[j+1] = X
end

function Base.getindex(N::NodalBasis{T}, l::Int, j::Int)::T where {T}
	@assert 0 <= l <= maxlevel(N)
	@assert 0 <= j <= (2^maxlevel(N))
	J = 2^(maxlevel(N)-l)*j
	return N.values[J+1]
end

# function Base.abs(H::HierarchicalBasis{T})::T where {T}
# 	return abs.(H.levels.coefficients)
# end

#************ BASIS TRANSFORMATIONS ********************************

"""
Evaluating hat functions at each level l, at position j, at point x
- on the grid, j can't be larger than the total length
- j must be odd because those are the only places where there are peaks (for the hat functions)
"""
function basis(l::Int, j::Int, x::T)::T where {T}
	@assert 0 <= l
	@assert 0 <= j <= (2^l)
	h = 2.0^(-l)
	if (h*(j - 1) <= x <= h*(j + 1)) && (0 <= x <= 1)
		value = 1 - abs((x/h) - j)
	else
		value = 0
	end
	return value
end

function x(l::Int, j::Int)::Float64
	h = 1/(2^l)
	return h*j
end

"""
The coefficients of a hierarchical basis are added
"""
function evaluate(H::HierarchicalBasis{T}, x::T)::T where {T}
	value = T(0)
	for l in 0:maxlevel(H)
		for j in 0:2^l
			value = value + H[l, j]*basis(l, j, x)
		end
	end
	return value
end

function Base.zeros(H::HierarchicalBasis{T})::NodalBasis{T} where {T}
	return NodalBasis(zeros(T, 2^maxlevel(H) + 1))
end

function Base.zeros(N::NodalBasis{T})::HierarchicalBasis{T} where {T}
	H = []
	for l in 0:maxlevel(N)
		H = vcat(H, Level(zeros(T, 2^l + 1)))
	end
	return HierarchicalBasis{T}(H)
end

function H_2_Nodal(H::HierarchicalBasis{T})::NodalBasis{T} where {T}
  	N = zeros(H)
 	for j in 0:2^maxlevel(H)
 		N[j] = evaluate(H, x(maxlevel(H),j))
 	end
 	return N
 end

function Nodal_2_H(N::NodalBasis{T})::HierarchicalBasis{T} where {T}
 	H = zeros(N)
  	for l in 0:maxlevel(N), j in 0:2^l
 		H[l,j] = N[l,j] - evaluate(H, x(l,j))
 	end
 	return H
 end

 function Nodal_2_H_new(N::NodalBasis{T})::HierarchicalBasis{T} where {T}
 	 H = zeros(N)
 	 H[0,0] = N[0,0]
 	 H[0,1] = N[0,1]
 	 for l in 1:maxlevel(N), j in 1:2^l-1
 		 H[l,j] = N[l,j] - interpolate((x(l, j-1), N[l, j-1]), (x(l, j+1), N[l, j+1]), x(l,j))
 	 end
 	 return H
 end


text:
abstract type GradientAlgorithm <: MLAlgorithm end

# function gradient!(alg::GradientAlgorithm, model, data)
#     A_mul_B!(alg.nvec, data.x, model.β)          # nvec ← x * β
#     deriv!(alg.nvec, model.loss, data.y, alg.nvec) # nvec ← deriv(L, y, x * β)
#     multiply_by_weights!(alg.nvec, data.w)   # nvec .*= w ./ n
#     At_mul_B!(alg.pvec, o.x, alg.nvec)      # pvec ← x'nvec
# end

# multiply_by_weights!(nvec, w::Void) = scale!(nvec, 1 / length(nvec))
# function multiply_by_weights!(nvec, w)
#     wt = inv(length(nvec))
#     for i in eachindex(nvec)
#         @inbounds nvec[i] *= w[i] * wt
#     end
# end

mutable struct AdaptiveProxGrad <: GradientAlgorithm
    nvec::Vector{Float64}
    pvec::Vector{Float64}
    step::Vector{Float64}
    divisor::Float64
    function AdaptiveProxGrad(initstep = 1.0, divisor = 1.5)
        divisor > 1  || error("divisor must be > 1")
        initstep > 0 || error("step size must be > 0")
        new(zeros(0), zeros(0), fill(initstep, 1), divisor)
    end
end

function setup!(alg::AdaptiveProxGrad, mod::AbstractMLModel, data::MLData)
    _modelsetup!(mod, data)
    alg.nvec = zeros(nobs(data))
    alg.pvec = zeros(npredictors(data))
    alg.step = fill(alg.step[1], npredictors(data))
end


text:
using Statistics

include("common.jl")

example = "16,1,2,0,4,2,7,1,2,14"

get_locations(input::AbstractString) = map(x -> parse(Int64, x), split(input, ","))

fuel_cost(locations::AbstractVector{Int64}, x::Real) = sum(abs.(locations .- x))

function get_best_cost(locations::AbstractVector{Int64})
    optimal_loc = Int64(median(locations))
    return fuel_cost(locations, optimal_loc)
end

get_ans(input::AbstractString) = get_best_cost(get_locations(input))

@assert get_ans(example) == 37
get_ans(load_input(7))

# Part 2

function fuel_cost2(loc::Int64, x::Real)
    d = abs(loc - x)
    return Int(d * (d + 1) / 2)
end
fuel_cost2(locations::AbstractVector{Int64}, x::Real) = sum(fuel_cost2.(locations, x))

function get_best_cost2(locations::AbstractVector{Int64})
    # We're likely to be *close* to the mean, so start there and walk (since the objective
    # is convex).
    loc = Int64(round(mean(locations)))

    cost(x) = fuel_cost2(locations, x)

    while cost(loc + 1) < cost(loc)
        loc += 1
    end

    while cost(loc - 1) < cost(loc)
        loc -= 1
    end

    return cost(loc)
end

get_ans2(input::AbstractString) = get_best_cost2(get_locations(input))

@assert get_ans2(example) == 168
get_ans2(load_input(7))


text:
# For vectors
function permute_dims_expr(perm, D, C, B, R)
  s = Expr(:tuple) # size
  x = Expr(:tuple) # stride
  o = Expr(:tuple) # offsets
  Rnew = Expr(:tuple) # rank
  Dnew = Expr(:tuple) # dense
  Cnew = -1
  Bnew = -1
  N = length(perm)
  for n ∈ 1:N
    p = perm[n]
    push!(s.args, Expr(:ref, :s, p))
    push!(x.args, Expr(:ref, :x, p))
    push!(o.args, Expr(:ref, :o, p))
    push!(Rnew.args, R[p])
    push!(Dnew.args, D[p])
    if C == p
      Cnew = n
    end
    if B == p
      Bnew = n
    end
  end
  Dnew, Cnew, Bnew, Rnew, s, x, o
end
@generated function Base.permutedims(A::PtrArray{S,D,T,N,C,B,R}, ::Val{P}) where {S,D,T,N,C,B,R,P}
  Dnew, Cnew, Bnew, Rnew, s, x, o = permute_dims_expr(P, D, C, B, R)
  quote
    $(Expr(:meta,:inline))
    s = size(A)
    ptr = stridedpointer(A)
    x = strides(ptr)
    o = offsets(ptr)
    si = StrideIndex{$N,$Rnew,$Cnew}($x, $o)
    sp = stridedpointer(ptr.p, si, StaticInt{$B}())
    PtrArray(sp, $s, Val{$Dnew}())
  end
end

@inline Base.adjoint(A::AbstractStrideMatrix) = permutedims(A, Val{(2,1)}())
@inline Base.transpose(A::AbstractStrideMatrix) = permutedims(A, Val{(2,1)}())

@generated function Base.adjoint(a::PtrArray{S,D,T,1,C,B,R,X,O}) where {S,D,T,C,B,R,X,O}
  s = Expr(:tuple, :(One()), Expr(:ref, :s, 1))
  x₁ = Expr(:ref, :x, 1)
  x = Expr(:tuple, x₁, x₁)
  o = Expr(:tuple, :(One()), Expr(:ref, :o, 1))
  R1 = R[1]
  Rnew = Expr(:tuple, R1+1, R1)
  Dnew = Expr(:tuple, true, D[1])
  Cnew = C == 1 ? 2 : C
  quote
    $(Expr(:meta,:inline))
    s = size(a)
    ptr = stridedpointer(a)
    x = strides(ptr)
    o = offsets(ptr)
    si = StrideIndex{2,$Rnew,$Cnew}($x, $o)
    sp = stridedpointer(ptr.p, si, StaticInt{$B}())
    PtrArray(sp, $s, Val{$Dnew}())
  end
end

@inline Base.transpose(a::AbstractStrideVector) = adjoint(a)

@inline row_major(a::AbstractStrideVector) = a
@inline row_major(A::AbstractStrideMatrix) = A
@inline row_major(A::AbstractStrideArray{S,D,T,N}) where {S,D,T,N} = permutedims(A, Val(ntuple(Base.Fix1(-,N+1),Val(N))))


text:
#=
Copyright 2021 BlackRock, Inc.

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.
=#

module PiecewiseQuadratics

import Base: sum, length, iterate, getindex, setindex!, lastindex, push!, append!
import Base: copy, zero, convert, +, *, ==, ≈, -, ∈, ⊆, <, ≤, ≥, >, ∩, display, show,
             reverse
import Base: pop!, isempty, intersect, empty!
using Printf

# Basic types
export Interval, BoundedQuadratic, PiecewiseQuadratic, FixedMemoryPwq, PwqSumWorkspace

# Operations
export <, >, ≤, ≥, ≲, ≳, ==, ≈, ⊆, ∈
export isempty, is_point, is_almost_point, continuous_and_overlapping, is_convex
export +, -, *, ∩
export sum, intersect, minimize, derivative, _sum

export copy, copy!, display, show, domain
export zero, indicator, get_line, get_tangent, domain

export reverse, scale, perspective, shift, tilt, restrict_dom, extend_dom
export reverse!, scale!, perspective!, shift!, tilt!, restrict_dom!, extend_dom!

export envelope, append_envelope!

export getindex, setindex!, lastindex, iterate, length, pop!, push!, append!
export mult_scalar!, simplify, prox, get_plot

include("intervals.jl")
include("utils.jl")
include("bounded_quadratics.jl")
include("piecewise_quadratics.jl")
include("fixed_memory_pwq.jl")
include("envelope.jl")

end


text:
"""
rna.jl

Fit G state models (generalized telegraph models) to RNA abundance data
For single cell RNA (scRNA) technical noise is included as a yieldfactor
single molecule FISH (smFISH) is treated as loss less
"""

# Load data, model, and option structures
"""
transient_rna(nsets::Int,control::String,treatment::String,time::Float64,gene::String,r::Vector,decayprior::Float64,yieldprior::Float64,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp::Float64=10.,method::Int=1,warmupsteps=0,annealsteps=0)
Fit transient G model to time dependent mRNA data
"""

function transient_rna(path,name::String,time,gene::String,nsets::Int,r::Vector,decayprior::Float64,yieldprior::Float64,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp::Float64=10.,method::Int=1,warmupsteps=0,annealsteps=0)
    data = data_rna(path,name,time,gene,false)
    model = model_rna(r,G,nalleles,nsets,cv,fittedparam,decayprior,yieldprior,method)
    options = MHOptions(samplesteps,annealsteps,warmupsteps,maxtime,temp,temp)
    return data,model,options
end
function transient_fish(path,name::String,time,gene::String,r::Vector,decayprior::Float64,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp::Float64=10.,method::Int=1,warmupsteps=0,annealsteps=0)
    data = data_rna(path,name,time,gene,true)
    model,options = transient_fish(r,decayprior,G,nalleles,fittedparam,cv,maxtime,samplesteps,temp,method,warmupsteps,annealsteps)
    return data,model,options
end
function transient_fish(path,name::String,time,gene::String,r::Vector,decayprior::Float64,delayprior::Float64,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp::Float64=10.,method::Int=1,warmupsteps=0,annealsteps=0)
    data = data_rna(path,name,time,gene,true)
    model,options = transient_fish(r,decayprior,delayprior,G,nalleles,fittedparam,cv,maxtime,samplesteps,temp,method,warmupsteps,annealsteps)
    return data,model,options
end
function transient_fish(r::Vector,decayprior::Float64,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp::Float64=10.,method::Int=1,warmupsteps=0,annealsteps=0)
    model = model_rna(r,G,nalleles,2,cv,fittedparam,decayprior,method)
    options = MHOptions(samplesteps,annealsteps,warmupsteps,maxtime,temp,temp)
    return model,options
end
function transient_fish(r::Vector,decayprior::Float64,delayprior::Float64,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp::Float64=10.,method::Int=1,warmupsteps=0,annealsteps=0)
    model = model_delay_rna(r,G,nalleles,2,cv,fittedparam,decayprior,delayprior)
    options = MHOptions(samplesteps,annealsteps,warmupsteps,maxtime,temp,temp)
    return model,options
end
function transient_rnafish(path,name::String,time,gene::String,nsets::Int,r::Vector,decayprior::Float64,yieldprior::Float64,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp::Float64=10.,method::Int=1,warmupsteps=0,annealsteps=0)
    data = data_rna(path,name,time,gene)
    model,options = transient_rnafish(r,decayprior,yieldprior,G,nalleles,fittedparam,cv,maxtime,samplesteps,temp,method,warmupsteps,annealsteps)
    return data,model,options
end
function transient_rnafish(r::Vector,decayprior::Float64,yieldprior::Float64,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp::Float64=10.,method::Int=1,warmupsteps=0,annealsteps=0)
    model = model_rna(r,G,nalleles,nsets,cv,fittedparam,decayprior,yieldprior,method)
    options = MHOptions(samplesteps,annealsteps,warmupsteps,maxtime,temp,temp)
    return model,options
end

"""
steadystate_rna(nsets::Int,file::String,gene::String,r::Vector,decayprior::Float64,yieldprior::Float64,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp=10.,warmupsteps=0,annealsteps=0)
Fit G model to steady state data
"""
function steadystate_rna(path,name::String,gene::String,nsets,r::Vector,decayprior::Float64,yieldprior::Float64,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp=10.,warmupsteps=0,annealsteps=0)
    data = data_rna(path,name,gene,false)
    model = model_rna(r,G,nalleles,nsets,cv,fittedparam,decayprior,yieldprior,0)
    options = MHOptions(samplesteps,annealsteps,warmupsteps,maxtime,temp,temp)
    return data,model,options
end
function steadystate_rna(path,name::String,gene::String,nsets,r::Vector,decayprior::Float64,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp=10.,warmupsteps=0,annealsteps=0)
    data = data_rna(path,name,gene,false)
    model = model_rna(r,G,nalleles,nsets,cv,fittedparam,decayprior,0)
    options = MHOptions(samplesteps,annealsteps,warmupsteps,maxtime,temp,temp)
    return data,model,options
end
function steadystate_rnafish(path,name::String,gene::String,fish::Array,r::Vector,decayprior::Float64,noisepriors::Vector,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp=10.,method=1,warmupsteps=0,annealsteps=0)
    data = data_rna(path,name,gene,fish)
    model = model_rna(r,G,nalleles,cv,fittedparam,decayprior,noisepriors,method)
    options = MHOptions(samplesteps,annealsteps,warmupsteps,maxtime,temp,temp)
    return data,model,options
end
function thresholds_fish(path,name::String,gene::String,r::Vector,decayprior::Float64,noisepriors::Vector,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp=10.,warmupsteps=0,annealsteps=0)
    data = data_rna(path,name,gene,true)
    model,options = thresholds_fish(r,decayprior,noisepriors,G,nalleles,fittedparam,cv,maxtime,samplesteps,temp,warmupsteps,annealsteps)
    return data,model,options
end
function thresholds_fish(r::Vector,decayprior::Float64,noisepriors::Vector,G::Int,nalleles::Int,fittedparam::Vector,cv,maxtime::Float64,samplesteps::Int,temp=10.,warmupsteps=0,annealsteps=0)
    model = model_rna(r,G,nalleles,cv,fittedparam,decayprior,noisepriors,0)
    options = MHOptions(samplesteps,annealsteps,warmupsteps,maxtime,temp,temp)
    return model,options
end

#Prepare data structures
"""
data_rna(path,time,gene,time)
data_rna(path,time,gene)
Load data structure
"""
function data_rna(path,name,time,gene::String,fish::Bool)
    len,h = histograms_rna(path,gene,fish)
    TransientRNAData(name,gene,len,time,h)
end
function data_rna(path,name,gene::String,fish::Bool)
    len,h = histograms_rna(path,gene,fish)
    RNAData(name,gene,len,h)
end
function data_rna(path,name,gene::String,fish::Array{Bool,1})
    len,h = histograms_rna(path,gene,fish)
    RNAMixedData(name,gene,len,fish,h)
end
"""
histograms_rna(path,gene)
prepare mRNA histograms
"""
function histograms_rna(path::Array,gene::String,fish::Bool)
    n = length(path)
    h = Array{Array,1}(undef,n)
    lengths = Array{Int,1}(undef,n)
    for i in eachindex(path)
        lengths[i], h[i] = histograms_rna(path[i],gene,fish)
    end
    return lengths,h
end
function histograms_rna(path::String,gene::String,fish::Bool)
    if fish
        h = read_fish(path,gene,.98)
    else
        h = read_scrna(path,.99)
        if length(h) == 0
            throw("data not found")
        end
    end
    return length(h),h
end
function histograms_rna(path::Array,gene::String,fish::Array{Bool,1})
    n = length(path)
    h = Array{Array,1}(undef,n)
    lengths = Array{Int,1}(undef,n)
    for i in eachindex(path)
        lengths[i], h[i] = histograms_rna(path[i],gene,fish[i])
    end
    return lengths,h
end
# Prepare model structures
"""
model_rna(r,G,nalleles,nsets,propcv,fittedparam,decayprior,yieldprior,method)
model_rna(r,G,nalleles,nsets,propcv,fittedparam,decayprior,method)
model_rna(r,G,nalleles,nsets,propcv,fittedparam,decayprior,noisepriors,method)
Load model structure
"""
function model_rna(r::Vector,G::Int,nalleles::Int,nsets::Int,propcv,fittedparam::Array,decayprior::Float64,yieldprior::Float64,method::Int)
    # propcv = proposal_cv_rna(propcv,fittedparam)
    d = prior_rna(r,G,nsets,fittedparam,decayprior,yieldprior)
    GMlossmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,method)
end
function model_rna(r::Vector,G::Int,nalleles::Int,nsets::Int,propcv,fittedparam::Array,randomeffects::Tuple,decayprior::Float64,yieldprior::Float64,method::Int)
    # propcv = proposal_cv_rna(propcv,fittedparam)
    d = prior_rna(r,G,nsets,fittedparam,decayprior,yieldprior)
    GMfixedeffectslossmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,randomeffects,method)
end
function model_rna(r::Vector,G::Int,nalleles::Int,nsets::Int,propcv,fittedparam::Array,decayprior::Float64,method::Int)
    # propcv = proposal_cv_rna(propcv,fittedparam)
    d = prior_rna(r,G::Int,nsets,fittedparam,decayprior)
    GMmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,method)
end
function model_rna(r::Vector,G::Int,nalleles::Int,propcv,fittedparam::Array,decayprior::Float64,noisepriors::Array,method::Int)
    # propcv = proposal_cv_rna(propcv,fittedparam)
    d = prior_rna(r,G::Int,1,fittedparam,decayprior,noisepriors)
    if method == 1
        GMrescaledmodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,method)
    else
        GMmultimodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),typeof(method)}(G,nalleles,r,d,propcv,fittedparam,method)
    end
end

function model_delay_rna(r::Vector,G::Int,nalleles::Int,nsets::Int,propcv,fittedparam::Array,decayprior,delayprior)
    # propcv = proposal_cv_rna(propcv,fittedparam)
    d = prior_rna(r,G,nsets,fittedparam,decayprior,delayprior)
    GMdelaymodel{typeof(r),typeof(d),typeof(propcv),typeof(fittedparam),Int64}(G,nalleles,r,d,propcv,fittedparam,1)
end

"""
prior_rna(r::Vector,G::Int,nsets::Int,propcv,fittedparam::Array,decayprior,yieldprior)
prepare prior distribution
r[mod(1:2*G,nsets)] = rates for each model set  (i.e. rates for each set are stacked)
r[2*G*nsets + 1] == yield factor (i.e remaining after loss due to technical noise)
"""
function prior_rna(r::Vector,G::Int,nsets::Int,fittedparam::Array,decayprior::Float64,yieldprior::Float64,f=Gamma)
        ind = 2*G * nsets + 1
        if length(r) == ind
            rm,rcv = setpriorrate(G,nsets,decayprior,yieldprior)
            if in(ind,fittedparam)
                return distributionBeta_array(rm[fittedparam],rcv[fittedparam],findfirst(ind.==fittedparam),f)
            else
                return distribution_array(rm[fittedparam],rcv[fittedparam],f)
            end
        else
            throw("rates have wrong length")
        end
end
function prior_rna(r::Vector,G::Int,nsets::Int,fittedparam::Array,decayprior::Float64,f=Gamma)
        if length(r) == 2*G*nsets
            rm,rcv = setpriorrate(G,nsets,decayprior)
            return distribution_array(rm[fittedparam],rcv[fittedparam],f)
        else
            throw("rates have wrong length")
        end
end
"""
prior_rna(r::Vector,G::Int,nsets::Int,fittedparam::Array,decayprior::Float64,noisepriors::Array)

prior for multithresholded smFISH data
r[1:2G] = model rates
r[2G+1] = additive noise mean
r[2G + 1 + 1:length(noisepriors)] = remaining fraction after thresholding (i.e. yield)
"""
function prior_rna(r::Vector,G::Int,nsets::Int,fittedparam::Array,decayprior::Float64,noisepriors::Array)
        if length(r) == 2*G * nsets + length(noisepriors)
            rm,rcv = setpriorrate(G,nsets,decayprior,noisepriors)
            return distribution_array(rm[fittedparam],rcv[fittedparam])
        else
            throw("rates have wrong length")
        end
end
"""
function setpriorrate(G)
Set prior distribution for mean and cv of rates
"""
function setpriorrate(G::Int,nsets::Int,decayrate::Float64,yieldfactor::Float64)
    rm,rcv = setpriorrate(G,nsets,decayrate)
    return [rm;yieldfactor],[rcv;.25]
end
function setpriorrate(G::Int,nsets::Int,decayrate::Float64)
    r0 = [.01*ones(2*(G-1));.1;decayrate]
    rc = [.1*ones(2*(G-1));.01;0.01]
    rm = r0
    rcv = rc
    for i in 2:nsets
        rm = vcat(rm,r0)
        rcv = vcat(rcv,rc)
    end
    return rm,rcv
end

function setpriorrate(G::Int,nsets::Int,decayrate::Float64,noisepriors::Array)
    rm,rcv = setpriorrate(G,nsets,decayrate)
    for nm in noisepriors
        rm = vcat(rm,nm)
        rcv = vcat(rcv,.25)
    end
    return rm,rcv
end
"""
proposal_cv_rna(propcv)
set propcv to a vector or matrix
"""
proposal_cv_rna(cv) = typeof(cv) == Float64 ? propcv*ones(length(fittedparam)) : propcv

proposal_cv_rna(propcv,fittedparam) = typeof(propcv) == Float64 ? propcv*ones(length(fittedparam)) : propcv

"""
fittedparam_rna(G,nsets,loss)
select all parameters to be fitted except decay rates
yieldfactor is last parameter in loss models
"""
function fittedparam_rna(G,nsets,loss)
    fp = fittedparam_rna(G,nsets)
    if loss == 1
        return vcat(fp,nsets*2*G+1)
    else
        return fp
    end
end
function fittedparam_rna(G,nsets)
    nrates = 2*G  # total number of rates in GM model
    k = nrates - 1  # adjust all rate parameters except decay time
    fp = Array{Int,1}(undef,nsets*k)
    for j in 0:nsets-1
        for i in 1:k
            fp[k*j + i] = nrates*j + i
        end
    end
    return fp
end

"""
rescale_rate_rna(r,G,decayrate::Float64)

set new decay rate and rescale
transition rates such that steady state distribution is the same
"""
function rescale_rate_rna(r,G,decayrate::Float64)
    rnew = copy(r)
    if mod(length(r),2*G) == 0
        rnew *= decayrate/r[2*G]
    else
        stride = fld(length(r),fld(length(r),2*G))
        for i in 0:stride:length(r)-1
            rnew[i+1:i+2*G] *= decayrate/r[2*G]
        end
    end
    return rnew
end

# Read in data and construct histograms
"""
read_scrna(filename::String,yield::Float64=.99,nhistmax::Int=1000)
Construct mRNA count per cell histogram array of a gene
"""
function read_scrna(filename::String,threshold::Float64=.99,nhistmax::Int=500)
    if isfile(filename) && filesize(filename) > 0
        x = readdlm(filename)[:,1]
        x = truncate_histogram(x,threshold,nhistmax)
        if x == 0
            dataFISH = Array{Int,1}(undef,0)
        else
            dataFISH = x
        end
        return dataFISH
    else
        # println("data file not found")
        return Array{Int,1}(undef,0)
    end
end

"""
read_fish(path,gene,threshold)
Read in FISH data from 7timepoint type folders

"""
function read_fish(path::String,cond::String,threshold::Float64=.98)
    xr = zeros(1000)
    lx = 0
    for (root,dirs,files) in walkdir(path)
        for file in files
            target = joinpath(root, file)
            if occursin(cond,target) && occursin("cellular",target)
                println(target)
                x1 = readdlm(target)[:,1]
                x1 = truncate_histogram(x1,threshold,1000)
                lx = length(x1)
                println(lx)
                xr[1:min(lx,1000)] += x1[1:min(lx,1000)]
            end
        end
    end
    return truncate_histogram(xr,1.0,1000)
end

function read_fish(path1::String,cond1::String,path2::String,cond2::String,threshold::Float64=.98)
    x1 = read_fish(path1,cond1,threshold)
    x2 = read_fish(path2,cond2,threshold)
    combine_histogram(x1,x2)
end

"""
plot_histogram_rna()

functions to plot data and model predicted histograms

"""
function plot_histogram_rna(gene::String,cond::String,datapath::String)
    file = joinpath(datapath,gene * "_" * cond * ".txt")
    h = read_scrna(file)
    plot(normalize_histogram(h))
    return h
end

function plot_histogram_fish(gene::String,datapaths::Array,modelfile::String,time=[0.;30.;120.],fittedparam=[7;8;9;10;11])
    r = readrates(modelfile,1)
    data,model,_ = transient_fish(datapaths,"",time,gene,r,1.,3,2,fittedparam,1.,1.,10)
    h=likelihoodarray(r[fittedparam],data,model)
    figure(gene)
    for i in eachindex(h)
        plot(h[i])
        plot(normalize_histogram(data.histRNA[i]))
    end
    return h
end

function plot_histogram_rna(gene,cond,G,nalleles,label,datafolder,folder,root)
    hn = get_histogram_rna(gene,cond,datafolder,root)
    ratepath = ratepath_Gmodel(gene,cond,G,nalleles,label,folder,root)
    println(ratepath)
    r = readrates(ratepath)
    m  = steady_state(r[1:2*G],r[end],G-1,length(hn),nalleles)
    plot(m)
    plot(hn)
    return r,hn,m,deviance(m,h),deviance(m,mediansmooth(h,3))
end

function get_histogram_rna(gene,cond,datafolder,root)
    datapath = scRNApath(gene,cond,datafolder,root)
    h = read_scrna(datapath,.99)
    normalize_histogram(h)
end

function plot_transient_rna(gene,cond,G,nalleles,label,datafolders::Vector,folder,root)
    ratepath = ratepath_Gmodel(gene,cond,G,nalleles,label,folder,root)
    println(ratepath)
    r = readrates(ratepath)
    maxdata = nhist_loss(10,r[end])
    h0 = steady_state_full(r[1:2*G],G-1,maxdata)
    h = transient([0.,30.,120.],r[2*G+1:4*G],r[end],G-1,nalleles,h0,1)
    i = 0
    for datafolder in datafolders
        i += 1
        figure()
        plot(h[i])
        plot(get_histogram_rna(gene,cond,datafolder,root))
    end
    return r,h
end

# functions to build paths to RNA/FISH data and results

function scRNApath(gene::String,cond::String,datapath::String,root::String)
    datapath = joinpath(root,datapath)
    if cond == ""
        joinpath(datapath,gene * ".txt")
    else
        joinpath(datapath,gene * "_" * cond * ".txt")
    end
end

scRNApath(gene,cond,datapath) = joinpath(datapath,gene * "_" * cond * ".txt")

FISHpath(gene,cond,datapath,root) = joinpath(joinpath(root,datapath),gene)

function ratepath_Gmodel(gene::String,cond::String,G::Int,nalleles::Int,label,folder,root)
    path_Gmodel("rates",gene,G,nalleles,label * "_" * cond,folder,root)
end

function path_Gmodel(type,gene::String,G::Int,nalleles::Int,label::String,folder,root)
    filelabel = label  * "_" * gene *  "_" * "$G" * "_" * "$nalleles" * ".txt"
    ratefile = type * "_" * filelabel
    joinpath(root, joinpath(folder,ratefile))
end

function stats_rna(genes::Vector,conds::Vector,datapath)
    z = Vector{Float64}(undef,0)
    h1 = Vector{Float64}(undef,0)
    h2 = Vector{Float64}(undef,0)
    for gene in genes
        t,m1,m2 = expression_rna(gene,conds,datapath)
        push!(z,t)
        push!(h1,m2)
        push!(h2,m1)
    end
    return z,h1,h2
end

function expression_rna(gene,conds::Vector,folder,threshold=.99)
    h = Vector{Vector}(undef,2)
    for i in eachindex(conds)
        datapath = scRNApath(gene,conds[i],folder)
        h[i] = read_scrna(datapath,threshold)
    end
    if length(h[1]) > 0 && length(h[1]) > 0
        return log_2sample(h[1],h[2]), mean_histogram(h[1]), mean_histogram(h[2])
    else
        return 0,0,0
    end
end

max_stars_repo_name: nih-niddk-mbs/StochasticGene

max_issues_repo_name: nih-niddk-mbs/StochasticGene

max_forks_repo_name: nih-niddk-mbs/StochasticGene


text:
# one parameter predicates

"""true iff midpoint(x) and radius(x) are zero"""
function iszero{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_zero), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

iszero{T<:Real}(x::T) = (x == zero(T))

"""true iff midpoint(x) or radius(x) are not zero"""
function notzero{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_zero), Int, (Ptr{ArbFloat},), &x)
    z == 0
end

notzero{T<:Real}(x::T) = (x != zero(T))

"""true iff zero is not within [upperbound(x), lowerbound(x)]"""
function nonzero{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_nonzero), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

nonzero{T<:Real}(x::T) = (x != zero(T))

"""true iff midpoint(x) is one and radius(x) is zero"""
function isone{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_one), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

isone{T<:Real}(x::T) = (x == one(T))

"""true iff midpoint(x) is not one or midpoint(x) is one and radius(x) is nonzero"""
function notone{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_one), Int, (Ptr{ArbFloat},), &x)
    z == 0
end

notone{T<:Real}(x::T) = (x != one(T))

"""true iff radius is zero"""
function isexact{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_exact), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

isexact{T<:Integer}(x::T) = true

"""true iff radius is nonzero"""
function notexact{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_exact), Int, (Ptr{ArbFloat},), &x)
    z == 0
end

notexact{T<:Integer}(x::T) = false

"""true iff midpoint(x) is an integer and radius(x) is zero"""
function isinteger{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_int), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

"""true iff midpoint(x) is not an integer or radius(x) is nonzero"""
function notinteger{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_int), Int, (Ptr{ArbFloat},), &x)
    z == 0
end

"""true iff lowerbound(x) is positive"""
function ispositive{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_positive), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

"""true iff upperbound(x) is negative"""
function isnegative{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_negative), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

"""true iff upperbound(x) is zero or negative"""
function notpositive{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_nonpositive), Int, (Ptr{ArbFloat},), &x)
    z != 0
end
excludesPositive{P}(x::ArbFloat{P}) = notpositive(x)

"""true iff lowerbound(x) is zero or positive"""
function notnegative{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_is_nonnegative), Int, (Ptr{ArbFloat},), &x)
    z != 0
end
excludesNegative{P}(x::ArbFloat{P}) = notnegative(x)

# two parameter predicates

"""true iff midpoint(x)==midpoint(y) and radius(x)==radius(y)"""
function isequal{P}(x::ArbFloat{P}, y::ArbFloat{P})
    z = ccall(@libarb(arb_equal), Int, (Ptr{ArbFloat}, Ptr{ArbFloat}), &x, &y)
    z != 0
end

"""true iff midpoint(x)!=midpoint(y) or radius(x)!=radius(y)"""
function notequal{P}(x::ArbFloat{P}, y::ArbFloat{P})
    z = ccall(@libarb(arb_equal), Int, (Ptr{ArbFloat}, Ptr{ArbFloat}), &x, &y)
    z == 0
end

"""true iff x and y have a common point"""
function overlap{P}(x::ArbFloat{P}, y::ArbFloat{P})
    z = ccall(@libarb(arb_overlaps), Int, (Ptr{ArbFloat}, Ptr{ArbFloat}), &x, &y)
    z != 0
end

"""true iff x and y have no common point"""
function donotoverlap{P}(x::ArbFloat{P}, y::ArbFloat{P})
    z = ccall(@libarb(arb_overlaps), Int, (Ptr{ArbFloat}, Ptr{ArbFloat}), &x, &y)
    z == 0
end

"""true iff x spans (covers) all of y"""
function contains{P}(x::ArbFloat{P}, y::ArbFloat{P})
    z = ccall(@libarb(arb_contains), Int, (Ptr{ArbFloat}, Ptr{ArbFloat}), &x, &y)
    z != 0
end

"""true iff y spans (covers) all of x"""
function iscontainedby{P}(x::ArbFloat{P}, y::ArbFloat{P})
    z = ccall(@libarb(arb_contains), Int, (Ptr{ArbFloat}, Ptr{ArbFloat}), &y, &x)
    z != 0
end

"""true iff x does not span (cover) all of y"""
function doesnotcontain{P}(x::ArbFloat{P}, y::ArbFloat{P})
    z = ccall(@libarb(arb_contains), Int, (Ptr{ArbFloat}, Ptr{ArbFloat}), &x, &y)
    z == 0
end

"""true iff y does not span (cover) all of x"""
function isnotcontainedby{P}(x::ArbFloat{P}, y::ArbFloat{P})
    z = ccall(@libarb(arb_contains), Int, (Ptr{ArbFloat}, Ptr{ArbFloat}), &y, &x)
    z == 0
end

"""true if it is quite likely that the arguments indicate the same value"""
function approxeq{P}(x::ArbFloat{P}, y::ArbFloat{P})
    if contains(x,y)
       true
    elseif contains(y,x)
       true
    else
       x = tidy(x)
       y = tidy(y)
       contains(x,y) || contains(y,x)
    end
end

(≊){P}(x::ArbFloat{P}, y::ArbFloat{P}) = approxeq(x,y)

"""true iff x contains an integer"""
function includesAnInteger{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_contains_int), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

"""true iff x does not contain an integer"""
function excludesIntegers{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_contains_int), Int, (Ptr{ArbFloat},), &x)
    z == 0
end

"""true iff x contains zero"""
function includesZero{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_contains_int), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

"""true iff x does not contain zero"""
function excludesZero{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_contains_int), Int, (Ptr{ArbFloat},), &x)
    z == 0
end

"""true iff x contains a positive value"""
function includesPositive{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_contains_positive), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

"""true iff x contains a negative value"""
function includesNegative{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_contains_negative), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

"""true iff x contains a nonpositive value"""
function includesNonpositive{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_contains_nonpositive), Int, (Ptr{ArbFloat},), &x)
    z != 0
end

"""true iff x contains a nonnegative value"""
function includesNonnegative{P}(x::ArbFloat{P})
    z = ccall(@libarb(arb_contains_nonnegative), Int, (Ptr{ArbFloat},), &x)
    z != 0
end


text:
# DATA_DIR = "simplecase"
DATA_DIR = "basecase"
# DATA_DIR = "basecase_lv" # same as basecase but with lower voltage limits
# DATA_DIR = "33buscase_pu"

include("data_handler.jl")

# Provides functions to read and manipulate case data
# Creates BUSES, LINES, GENERATORS as Dict
include("stochasticity.jl")
# Settings and inital calculations for load stochasticity
include("lindist.jl")
# Provides function Solve_LinDist(buses, lines, generators) for feasability testing
include("standard_cc.jl")
# Common CC OPF Formulation in JuMPChance
# provides function Solve_CCLinDist(buses, lines, generators, error_variances)
include("dd_dist_robust_cc.jl")
# provides function Solve_CCLinDist(buses, lines, generators, error_variances)
include("result_tester.jl")
# Provides function test_results(bus_results, line_results, samples)
include("soc_flow.jl")
# Provides function Solve_Dist_SOC(buses, lines, generators) with a SOC dist fow implementation
include("lin_soc_flow.jl")
# Provides the function function Solve_Linearized_Dist_SOC(buses, lines, generators, fp0, fq0, a0)
include("result_plotters.jl")
# Provedes some functiuons to vizualize results


text:
using Hwloc

mutable struct Topology
  topology::Union{Nothing,Hwloc.Object}
end
Topology() = Topology(nothing)
function safe_topology_load!()
  try
    TOPOLOGY.topology = Hwloc.gettopology()
  catch e
    @warn e
    @warn """
          Using Hwloc failed. Please file an issue with the above warning at: https://github.com/JuliaParallel/Hwloc.jl
          Proceeding with generic topology assumptions. This may result in reduced performance.
      """
  end
end

const TOPOLOGY = Topology()
safe_topology_load!()

function count_attr(topology::Hwloc.Object, attr::Symbol)
  count::Int = 0
  for t ∈ topology
    count += t.type_ === attr
  end
  if ((Sys.ARCH === :aarch64) && Sys.isapple()) && ((attr === :Core) | (attr === :PU))
    # count big cores only on Apple, as one process can only use big or small cores at a time.
    count >>>= one(count) # FIXME: how to actually identify big cores???
  end
  count
end
function count_attr(attr::Symbol)
  topology = TOPOLOGY.topology
  topology === nothing && return nothing
  count_attr(topology, attr)
end

@static if Sys.isapple() && Sys.ARCH === :aarch64 # detect M1
  redefine_attr_count() = nothing
  num_l1cache() = static(4)
  num_l2cache() = static(1)
  num_l3cache() = static(0)
  num_l4cache() = static(0)
  num_machines() = static(1)
  num_sockets() = static(1)
  num_cores() = static(4)
  sys_threads() = static(4)
  num_threads() = static(4)
else # not M1
  @noinline function define_attr_count(fname::Symbol, v)
    if v === nothing
      @eval $fname() = nothing
      return
    elseif v isa Integer
      @eval $fname() = StaticInt{$(convert(Int, v))}()
    elseif v isa Bool
      if v
        @eval $fname() = True()
      else
        @eval $fname() = False()
      end
    else
      @eval $fname() = $v
    end
    nothing
  end
  # define_attr(fname::Symbol, attr::Symbol) = define_attr_count(fname, )

  for (f, attr) ∈ [
    (:num_l1cache, :L1Cache),
    (:num_l2cache, :L2Cache),
    (:num_l3cache, :L3Cache),
    (:num_l4cache, :L4Cache),
    (:num_machines, :Machine),
    (:num_sockets, :Package),
    (:num_cores, :Core),
    (:sys_threads, :PU),
    (:num_threads, :PU),
  ]
    define_attr_count(f, count_attr(attr))
  end

  function redefine_attr_count()
    sys_thread = Int(sys_threads())::Int
    iter = [
      (num_l1cache(), :num_l1cache, :L1Cache),
      (num_l2cache(), :num_l2cache, :L2Cache),
      (num_l3cache(), :num_l3cache, :L3Cache),
      (num_l4cache(), :num_l4cache, :L4Cache),
      (num_machines(), :num_machines, :Machine),
      (num_sockets(), :num_sockets, :Package),
      (num_cores(), :num_cores, :Core),
      (sys_thread, :sys_threads, :PU),
    ]
    for (v, f, attr) ∈ iter
      ref = count_attr(attr)
      if ref ≠ v
        @debug "Redefining attr count $f = $ref."
        define_attr_count(f, ref)
      end
    end
    nothing
  end
  @noinline function _redefine_num_threads()
    @eval num_threads() = StaticInt{$(Threads.nthreads())}()
  end
  function redefine_num_threads()
    Int(num_threads()) > min(Threads.nthreads(), Int(sys_threads())) &&
      _redefine_num_threads()
  end
end # not M1

function dynamic_cache_inclusivity()::NTuple{4,Bool}
  @static if !((Sys.ARCH === :x86_64) || (Sys.ARCH === :i686))
    return (false, false, false, false)
  end
  function get_cache_edx(subleaf)
    # source: https://github.com/m-j-w/CpuId.jl/blob/401b638cb5a020557bce7daaf130963fb9c915f0/src/CpuInstructions.jl#L38
    # credit Markus J. Weber, copyright: https://github.com/m-j-w/CpuId.jl/blob/master/LICENSE.md
    Base.llvmcall(
      """
              ; leaf = %0, subleaf = %1, %2 is some label
              ; call 'cpuid' with arguments loaded into registers EAX = leaf, ECX = subleaf
              %2 = tail call { i32, i32, i32, i32 } asm sideeffect "cpuid",
                  "={ax},={bx},={cx},={dx},{ax},{cx},~{dirflag},~{fpsr},~{flags}"
                  (i32 4, i32 %0) #2
              ; retrieve the result values and return eax and edx contents
              %3 = extractvalue { i32, i32, i32, i32 } %2, 0
              %4 = extractvalue { i32, i32, i32, i32 } %2, 3
              %5  = insertvalue [2 x i32] undef, i32 %3, 0
              %6  = insertvalue [2 x i32]   %5 , i32 %4, 1
              ; return the value
              ret [2 x i32] %6
              """
      # llvmcall requires actual types, rather than the usual (...) tuple
      ,
      Tuple{UInt32,UInt32},
      Tuple{UInt32},
      subleaf % UInt32,
    )
  end
  # eax0, edx1 = get_cache_edx(0x00000000)
  t = (false, false, false, false)
  i = zero(UInt32)
  j = 0
  while (j < 4)
    eax, edx = get_cache_edx(i)
    i += one(UInt32)
    iszero(eax & 0x1f) && break
    iszero(eax & 0x01) && continue
    ci = ((edx & 0x00000002) != 0x00000000) & (eax & 0x1f != 0x00000000)
    t = Base.setindex(t, ci, (j += 1))
  end
  t
end

nothing_cache_summary() =
  (size = 0, linesize = 64, associativity = nothing, type = nothing, inclusive = nothing)
function dynamic_cache_summary(N)
  topology = TOPOLOGY.topology
  cache_name = (:L1Cache, :L2Cache, :L3Cache, :L4Cache)[N]
  if (topology === nothing) || (iszero(count_attr(cache_name)))
    return nothing_cache_summary()
  end
  c = first(t for t in topology if t.type_ == cache_name && t.attr.depth == N).attr
  (
    size = c.size,
    linesize = c.linesize,
    associativity = c.associativity,
    type = c.type_,
    inclusive = dynamic_cache_inclusivity()[N],
  )
end
cache_size(_) = StaticInt{0}()
cache_associativity(_) = nothing
cache_type(_) = nothing
cache_inclusive(_) = nothing

unwrap(::Val{S}) where {S} = S
@static if Sys.isapple() && Sys.ARCH === :aarch64
  cache_linesize(_) = StaticInt{128}() # assume...
  redefine_cache(_) = nothing
  cache_size(::Union{Val{1},StaticInt{1}}) = StaticInt{131072}()
  cache_linesize(::Union{Val{1},StaticInt{1}}) = StaticInt{128}()
  cache_associativity(::Union{Val{1},StaticInt{1}}) = StaticInt{0}()
  cache_type(::Union{Val{1},StaticInt{1}}) = Val{:Data}()
  cache_inclusive(::Union{Val{1},StaticInt{1}}) = False()

  cache_size(::Union{Val{2},StaticInt{2}}) = StaticInt{12582912}()
  cache_linesize(::Union{Val{2},StaticInt{2}}) = StaticInt{128}()
  cache_associativity(::Union{Val{2},StaticInt{2}}) = StaticInt{0}()
  cache_type(::Union{Val{2},StaticInt{2}}) = Val{:Unified}()
  cache_inclusive(::Union{Val{2},StaticInt{2}}) = False()

else # not M1
  cache_linesize(_) = StaticInt{64}() # assume...
  function define_cache(N, c = dynamic_cache_summary(N))
    c === nothing_cache_summary() || _define_cache(N, c)
  end
  @noinline function _define_cache(N, c)
    @eval begin
      cache_size(::Union{Val{$N},StaticInt{$N}}) = StaticInt{$(c.size)}()
      cache_linesize(::Union{Val{$N},StaticInt{$N}}) = StaticInt{$(c.linesize)}()
      cache_associativity(::Union{Val{$N},StaticInt{$N}}) = StaticInt{$(c.associativity)}()
      cache_type(::Union{Val{$N},StaticInt{$N}}) =
        Val{$(c.type === nothing ? nothing : QuoteNode(c.type))}()
      cache_inclusive(::Union{Val{$N},StaticInt{$N}}) =
        $(c.inclusive isa Bool && c.inclusive ? :True : :False)()
    end
    nothing
  end
  function redefine_cache(N)
    s = cache_size(StaticInt(N))
    l = cache_linesize(StaticInt(N))
    a = cache_associativity(StaticInt(N))
    t = cache_type(StaticInt(N))
    i = cache_inclusive(StaticInt(N))
    c = (
      size = Int(s)::Int,
      linesize = l === nothing ? nothing : Int(l)::Int,
      associativity = a === nothing ? nothing : Int(a)::Int,
      type = t === nothing ? nothing : unwrap(t)::Symbol,
      inclusive = i === nothing ? nothing : Bool(i)::Bool,
    )
    correct = dynamic_cache_summary(N)
    if c !== correct
      @debug "Redefining cache $N."
      _define_cache(N, correct)
    end
    nothing
  end
  foreach(define_cache, 1:4)
end

function __init__()
  Sys.isapple() && Sys.ARCH === :aarch64 && return # detect M1
  ccall(:jl_generating_output, Cint, ()) == 1 && return
  safe_topology_load!()
  if count_attr(:Core) ≢ BASELINE_CORES
    redefine_attr_count()
    foreach(redefine_cache, 1:4)
  end
  redefine_num_threads()
  return nothing
end


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/libim_convert_/include/im_convert.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static is_windows()? (const libim_convert_ = "im") : (const libim_convert_ = "libim")  # Name of IM shared lib.

function imImageCreate(width::Cint, height::Cint, color_space::Cint, data_type::Cint)
  ccall( (:imImageCreate, libim_convert_), Ptr{imImage}, (Cint, Cint, Cint, Cint), width, height, color_space, data_type)
end
function imImageInit(width::Cint, height::Cint, color_mode::Cint, data_type::Cint, data_buffer::Ptr{Void}, palette::Ptr{Clong}, palette_count::Cint)
  ccall( (:imImageInit, libim_convert_), Ptr{imImage}, (Cint, Cint, Cint, Cint, Ptr{Void}, Ptr{Clong}, Cint), width, height, color_mode, data_type, data_buffer, palette, palette_count)
end
function imImageCreateBased(image::Ptr{imImage}, width::Cint, height::Cint, color_space::Cint, data_type::Cint)
  ccall( (:imImageCreateBased, libim_convert_), Ptr{imImage}, (Ptr{imImage}, Cint, Cint, Cint, Cint), image, width, height, color_space, data_type)
end
function imImageDestroy(image::Ptr{imImage})
  ccall( (:imImageDestroy, libim_convert_), Void, (Ptr{imImage},), image)
end
function imImageAddAlpha(image::Ptr{imImage})
  ccall( (:imImageAddAlpha, libim_convert_), Void, (Ptr{imImage},), image)
end
function imImageSetAlpha(image::Ptr{imImage}, alpha::Cfloat)
  ccall( (:imImageSetAlpha, libim_convert_), Void, (Ptr{imImage}, Cfloat), image, alpha)
end
function imImageRemoveAlpha(image::Ptr{imImage})
  ccall( (:imImageRemoveAlpha, libim_convert_), Void, (Ptr{imImage},), image)
end
function imImageReshape(image::Ptr{imImage}, width::Cint, height::Cint)
  ccall( (:imImageReshape, libim_convert_), Void, (Ptr{imImage}, Cint, Cint), image, width, height)
end
function imImageCopy(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageCopy, libim_convert_), Void, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end
function imImageCopyData(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageCopyData, libim_convert_), Void, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end
function imImageCopyAttributes(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageCopyAttributes, libim_convert_), Void, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end
function imImageMergeAttributes(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imImageMergeAttributes, libim_convert_), Void, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end
function imImageCopyPlane(src_image::Ptr{imImage}, src_plane::Cint, dst_image::Ptr{imImage}, dst_plane::Cint)
  ccall( (:imImageCopyPlane, libim_convert_), Void, (Ptr{imImage}, Cint, Ptr{imImage}, Cint), src_image, src_plane, dst_image, dst_plane)
end
function imImageDuplicate(image::Ptr{imImage})
  ccall( (:imImageDuplicate, libim_convert_), Ptr{imImage}, (Ptr{imImage},), image)
end
function imImageClone(image::Ptr{imImage})
  ccall( (:imImageClone, libim_convert_), Ptr{imImage}, (Ptr{imImage},), image)
end
function imImageSetAttribute(image::Ptr{imImage}, attrib::Ptr{UInt8}, data_type::Cint, count::Cint, data::Ptr{Void})
  ccall( (:imImageSetAttribute, libim_convert_), Void, (Ptr{imImage}, Ptr{UInt8}, Cint, Cint, Ptr{Void}), image, attrib, data_type, count, data)
end
function imImageGetAttribute(image::Ptr{imImage}, attrib::Ptr{UInt8}, data_type::Ptr{Cint}, count::Ptr{Cint})
  ccall( (:imImageGetAttribute, libim_convert_), Ptr{Void}, (Ptr{imImage}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cint}), image, attrib, data_type, count)
end
function imImageGetAttributeList(image::Ptr{imImage}, attrib::Ptr{Ptr{UInt8}}, attrib_count::Ptr{Cint})
  ccall( (:imImageGetAttributeList, libim_convert_), Void, (Ptr{imImage}, Ptr{Ptr{UInt8}}, Ptr{Cint}), image, attrib, attrib_count)
end
function imImageClear(image::Ptr{imImage})
  ccall( (:imImageClear, libim_convert_), Void, (Ptr{imImage},), image)
end
function imImageIsBitmap(image::Ptr{imImage})
  ccall( (:imImageIsBitmap, libim_convert_), Cint, (Ptr{imImage},), image)
end
function imImageSetPalette(image::Ptr{imImage}, palette::Ptr{Clong}, palette_count::Cint)
  ccall( (:imImageSetPalette, libim_convert_), Void, (Ptr{imImage}, Ptr{Clong}, Cint), image, palette, palette_count)
end
function imImageMatchSize(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchSize, libim_convert_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end
function imImageMatchColor(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchColor, libim_convert_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end
function imImageMatchDataType(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchDataType, libim_convert_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end
function imImageMatchColorSpace(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatchColorSpace, libim_convert_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end
function imImageMatch(image1::Ptr{imImage}, image2::Ptr{imImage})
  ccall( (:imImageMatch, libim_convert_), Cint, (Ptr{imImage}, Ptr{imImage}), image1, image2)
end
function imImageSetMap(image::Ptr{imImage})
  ccall( (:imImageSetMap, libim_convert_), Void, (Ptr{imImage},), image)
end
function imImageSetBinary(image::Ptr{imImage})
  ccall( (:imImageSetBinary, libim_convert_), Void, (Ptr{imImage},), image)
end
function imImageSetGray(image::Ptr{imImage})
  ccall( (:imImageSetGray, libim_convert_), Void, (Ptr{imImage},), image)
end
function imImageMakeBinary(image::Ptr{imImage})
  ccall( (:imImageMakeBinary, libim_convert_), Void, (Ptr{imImage},), image)
end
function imImageMakeGray(image::Ptr{imImage})
  ccall( (:imImageMakeGray, libim_convert_), Void, (Ptr{imImage},), image)
end
function imFileLoadImage(ifile::Ptr{Cint}, index::Cint, error::Ptr{Cint})
  ccall( (:imFileLoadImage, libim_convert_), Ptr{imImage}, (Ptr{Cint}, Cint, Ptr{Cint}), ifile, index, error)
end
function imFileLoadImageFrame(ifile::Ptr{Cint}, index::Cint, image::Ptr{imImage}, error::Ptr{Cint})
  ccall( (:imFileLoadImageFrame, libim_convert_), Void, (Ptr{Cint}, Cint, Ptr{imImage}, Ptr{Cint}), ifile, index, image, error)
end
function imFileLoadBitmap(ifile::Ptr{Cint}, index::Cint, error::Ptr{Cint})
  ccall( (:imFileLoadBitmap, libim_convert_), Ptr{imImage}, (Ptr{Cint}, Cint, Ptr{Cint}), ifile, index, error)
end
function imFileLoadImageRegion(ifile::Ptr{Cint}, index::Cint, bitmap::Cint, error::Ptr{Cint}, xmin::Cint, xmax::Cint, ymin::Cint, ymax::Cint, width::Cint, height::Cint)
  ccall( (:imFileLoadImageRegion, libim_convert_), Ptr{imImage}, (Ptr{Cint}, Cint, Cint, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint), ifile, index, bitmap, error, xmin, xmax, ymin, ymax, width, height)
end
function imFileLoadBitmapFrame(ifile::Ptr{Cint}, index::Cint, image::Ptr{imImage}, error::Ptr{Cint})
  ccall( (:imFileLoadBitmapFrame, libim_convert_), Void, (Ptr{Cint}, Cint, Ptr{imImage}, Ptr{Cint}), ifile, index, image, error)
end
function imFileSaveImage(ifile::Ptr{Cint}, image::Ptr{imImage})
  ccall( (:imFileSaveImage, libim_convert_), Cint, (Ptr{Cint}, Ptr{imImage}), ifile, image)
end
function imFileImageLoad(file_name::Ptr{UInt8}, index::Cint, error::Ptr{Cint})
  ccall( (:imFileImageLoad, libim_convert_), Ptr{imImage}, (Ptr{UInt8}, Cint, Ptr{Cint}), file_name, index, error)
end
function imFileImageLoadBitmap(file_name::Ptr{UInt8}, index::Cint, error::Ptr{Cint})
  ccall( (:imFileImageLoadBitmap, libim_convert_), Ptr{imImage}, (Ptr{UInt8}, Cint, Ptr{Cint}), file_name, index, error)
end
function imFileImageLoadRegion(file_name::Ptr{UInt8}, index::Cint, bitmap::Cint, error::Ptr{Cint}, xmin::Cint, xmax::Cint, ymin::Cint, ymax::Cint, width::Cint, height::Cint)
  ccall( (:imFileImageLoadRegion, libim_convert_), Ptr{imImage}, (Ptr{UInt8}, Cint, Cint, Ptr{Cint}, Cint, Cint, Cint, Cint, Cint, Cint), file_name, index, bitmap, error, xmin, xmax, ymin, ymax, width, height)
end
function imFileImageSave(file_name::Ptr{UInt8}, format::Ptr{UInt8}, image::Ptr{imImage})
  ccall( (:imFileImageSave, libim_convert_), Cint, (Ptr{UInt8}, Ptr{UInt8}, Ptr{imImage}), file_name, format, image)
end
function imConvertDataType(src_image::Ptr{imImage}, dst_image::Ptr{imImage}, cpx2real::Cint, gamma::Cfloat, abssolute::Cint, cast_mode::Cint)
  ccall( (:imConvertDataType, libim_convert_), Cint, (Ptr{imImage}, Ptr{imImage}, Cint, Cfloat, Cint, Cint), src_image, dst_image, cpx2real, gamma, abssolute, cast_mode)
end
function imConvertColorSpace(src_image::Ptr{imImage}, dst_image::Ptr{imImage})
  ccall( (:imConvertColorSpace, libim_convert_), Cint, (Ptr{imImage}, Ptr{imImage}), src_image, dst_image)
end
function imConvertToBitmap(src_image::Ptr{imImage}, dst_image::Ptr{imImage}, cpx2real::Cint, gamma::Cfloat, abssolute::Cint, cast_mode::Cint)
  ccall( (:imConvertToBitmap, libim_convert_), Cint, (Ptr{imImage}, Ptr{imImage}, Cint, Cfloat, Cint, Cint), src_image, dst_image, cpx2real, gamma, abssolute, cast_mode)
end
function imImageGetOpenGLData(image::Ptr{imImage}, glformat::Ptr{Cint})
  ccall( (:imImageGetOpenGLData, libim_convert_), Ptr{Void}, (Ptr{imImage}, Ptr{Cint}), image, glformat)
end
function imImageCreateFromOpenGLData(width::Cint, height::Cint, glformat::Cint, gldata::Ptr{Void})
  ccall( (:imImageCreateFromOpenGLData, libim_convert_), Ptr{imImage}, (Cint, Cint, Cint, Ptr{Void}), width, height, glformat, gldata)
end
function imConvertPacking(src_data::Ptr{Void}, dst_data::Ptr{Void}, width::Cint, height::Cint, src_depth::Cint, dst_depth::Cint, data_type::Cint, src_is_packed::Cint)
  ccall( (:imConvertPacking, libim_convert_), Void, (Ptr{Void}, Ptr{Void}, Cint, Cint, Cint, Cint, Cint, Cint), src_data, dst_data, width, height, src_depth, dst_depth, data_type, src_is_packed)
end
function imConvertMapToRGB(data::Ptr{Cuchar}, count::Cint, depth::Cint, packed::Cint, palette::Ptr{Clong}, palette_count::Cint)
  ccall( (:imConvertMapToRGB, libim_convert_), Void, (Ptr{Cuchar}, Cint, Cint, Cint, Ptr{Clong}, Cint), data, count, depth, packed, palette, palette_count)
end
function imConvertRGB2Map(width::Cint, height::Cint, red::Ptr{Cuchar}, green::Ptr{Cuchar}, blue::Ptr{Cuchar}, map::Ptr{Cuchar}, palette::Ptr{Clong}, palette_count::Ptr{Cint})
  ccall( (:imConvertRGB2Map, libim_convert_), Cint, (Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Clong}, Ptr{Cint}), width, height, red, green, blue, map, palette, palette_count)
end


text:
# Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors    #src
# This Source Code Form is subject to the terms of the Mozilla Public License   #src
# v.2.0. If a copy of the MPL was not distributed with this file, You can       #src
# obtain one at https://mozilla.org/MPL/2.0/.                                   #src

# # The Rosenbrock function

# A nonlinear example of the classical Rosenbrock function.

using JuMP
import Ipopt
import Test

function example_rosenbrock()
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @variable(model, x)
    @variable(model, y)
    @NLobjective(model, Min, (1 - x)^2 + 100 * (y - x^2)^2)
    optimize!(model)

    Test.@test termination_status(model) == MOI.LOCALLY_SOLVED
    Test.@test primal_status(model) == MOI.FEASIBLE_POINT
    Test.@test objective_value(model) ≈ 0.0 atol = 1e-10
    Test.@test value(x) ≈ 1.0
    Test.@test value(y) ≈ 1.0
    return
end

example_rosenbrock()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#!/usr/bin/env julia

const DATA_DIR   = "../11_data"
const INPUT_FILE = "input.txt"

const SEAT_EMPTY    = 'L'
const SEAT_FLOOR    = '.'
const SEAT_OCCUPIED = '#'

# ----------------------------------------
macro assert(bool_expr)
    message = string("Assertion: ", bool_expr, " failed")
    :($(esc(bool_expr)) || error($message))
end
# ----------------------------------------

function nbhd_1(seats, row, col)
    row_max = length(seats)
    col_max = length(seats[1])

    nbrs = Tuple{Int, Int}[]
    for row_adj in (row - 1:row + 1)
        for col_adj in (col - 1:col + 1)
            if (row_adj == row && col_adj == col)
                continue
            end
            if (1 <= row_adj <= row_max && 1 <= col_adj <= col_max)
               push!(nbrs, (row_adj, col_adj))
            end
        end
    end
    return nbrs
end

function nbhd_2(seats, row, col)
    row_max = length(seats)
    col_max = length(seats[1])
    nbrs = []
    dirs = [[dr, dc] for dr in -1:1 for dc in -1:1 if (dr != 0 || dc != 0)]
    for dir in dirs
        current = [row, col]
        while true
            current += dir
            if !is_in_range(row_max, col_max, current[1], current[2])
                break
            end
            if ( is_empty(seats, current[1], current[2])
               || is_occ(seats, current[1], current[2])
               )
                push!(nbrs, current)
                break
            end
        end
    end
    return nbrs
end

function nbhd_occ_count(seats, row, col, nbhd_rule)
    nbhd_occ_count = 0
    nbhd = nbhd_rule(seats, row, col)
    for nbr in nbhd
        if is_occ(seats, nbr[1], nbr[2])
            nbhd_occ_count += 1
        end
    end
    return nbhd_occ_count
end

function get_next_seats(seats, nbhd_rule, nbhd_thresh, verbose=false)
    result = deepcopy(seats)
    for row in 1:length(seats)
        for col in 1:length(seats[1])
            if ( (is_empty(seats, row, col)
                 && nbhd_occ_count(seats, row, col, nbhd_rule) == 0)
               )
                result[row][col] = SEAT_OCCUPIED
            elseif ( is_occ(seats, row, col)
                       && nbhd_occ_count(seats, row, col, nbhd_rule) >= nbhd_thresh
                   )
                result[row][col] = SEAT_EMPTY
            end
        end
    end
    return result
end

function is_empty(seats, row, col)
    result = seats[row][col] == SEAT_EMPTY
    return result
end

function is_in_range(row_max, col_max, row, col)
    return 1 <= row <= row_max && 1 <= col <= col_max
end

function is_occ(seats, row, col)
    result = seats[row][col] == SEAT_OCCUPIED
    return result
end

function occs(seats)
    result = 0
    for row in 1:length(seats)
        for col in 1:length(seats[1])
            if is_occ(seats, row, col)
                result += 1
            end
        end
    end
    return result
end

function show(seats)
    for row in seats
        for seat in row
            print(seat)
        end
        println()
    end
end

function test_show_versions(seats)
    version = 1
    while version <= 6
        println("VERSION $version:")
        show(seats)

        seats = get_next_seats(seats)
        version += 1
        version < 6 && println("====================")
    end
end

# ----------------------------------------

function get_occupied_count(seats, nbhd_rule, nbhd_thresh)
    round_num = 0
    while true
        round_num += 1
        # println("Round # = $(round_num)")
        next_seats = get_next_seats(seats, nbhd_rule, nbhd_thresh)
        # println("INFO: next_seats={next_seats}, seats={seats}")
        if all(next_seats .== seats)
            println("Count of stably occupied seats=$(occs(seats))")
            break
        end
        seats = next_seats
    end
end

function main()
    filename = "$DATA_DIR/$INPUT_FILE"
    seats = map(s->[c for c in s], readlines(filename))
    println("========== PART 1 ==========")
    get_occupied_count(seats, nbhd_1, 4)
    println()
    println("========== PART 1 ==========")
    get_occupied_count(seats, nbhd_2, 5)
end

main()

max_stars_repo_name: jaycoskey/AdventOfCode

max_issues_repo_name: jaycoskey/AdventOfCode

max_forks_repo_name: jaycoskey/AdventOfCode


text:
# This file includes modified source code from https://github.com/odow/SDDP.jl
# as at 5b3ba3910f6347765708dd6e7058e2fcf7d13ae5

#  Copyright 2017, Oscar Dowson and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################

using SDDiP, JuMP, GLPKMathProgInterface, Base.Test, SDDP
# using Ipopt  # uncomment for the Level Method

srand(11111)

# Demand for newspapers
# There are two equally probable scenarios in each stage
#   Demand[stage, noise]
Demand = [
    10.0 15.0;
    12.0 20.0;
    8.0  20.0
]

# Markov state purchase prices
PurchasePrice = [5.0, 8.0]

RetailPrice = 7.0

# Transition matrix
Transition = Array{Float64, 2}[
    [1.0]',
    [0.6 0.4],
    [0.3 0.7; 0.3 0.7]
  ]

# Initialise SDDP Model
m = SDDPModel(
        sense             = :Max,
        stages            = 3,
        objective_bound   = 1000,
        markov_transition = Transition,
        solver            = GLPKSolverMIP()
                                                ) do sp, stage, markov_state

    # ====================
    #   State variable: binary decomposition of stock
    @binarystate(sp, stock <= 100, stock0 == 5, Int)

    # ====================
    #   Other variables
    @variables(sp, begin
        buy >= 0, Int  # Quantity to buy
        sell>= 0, Int  # Quantity to sell
    end)

    # ====================
    #   Scenarios
    @rhsnoises(sp, D=Demand[stage,:], begin
        sell <= D
        sell >= 0.5D
    end)

    # ====================
    #   Objective
    @stageobjective(sp, sell * RetailPrice - buy * PurchasePrice[markov_state])

    # ====================
    #   Dynamics constraint
    @constraint(sp, stock == stock0 + buy - sell)

    # Call to solve via Lagrangians
    # setSDDiPsolver!(sp, method=LevelMethod(quadsolver=IpoptSolver(print_level=0)))
    # setSDDiPsolver!(sp, method=SubgradientMethod(wait=10))
    setSDDiPsolver!(sp, method=KelleyMethod())

end

@time solvestatus = SDDP.solve(m,
    max_iterations = 12
)

@test isapprox(getbound(m), 97.9, atol=1e-3)