text:
#TODO: Extract tails (top 10-10, bottom 10-10)

#these functions build the symbols for each percentile ranking covering "period" years' return
function buildpercentilesym(period::Int, lbenchmark::NSymbol = nothing;
  lretsym::NSymbol = :lreturn)::Tuple{Symbol, Symbol}

  lretstring = string(something(lretsym, "zero"))
  lbenchstring::String = isnothing(lbenchmark) ? "" : string("net", lbenchmark)

  return (Symbol("p$lretstring", lbenchstring, period, "yr"),
    Symbol("L$(period)p$lretstring", lbenchstring, period, "yr"))
end

function buildpercentilesymtype(period::Int, lbenchmark::NSymbol = nothing;
  lretsym::NSymbol = :lreturn)::Tuple{Symbol, Symbol}

  lretstring = string(something(lretsym, "zero"))
  lbenchstring::String = isnothing(lbenchmark) ? "" : string("net", lbenchmark)

  return (Symbol("p$lretstring", lbenchstring, period, "yrtype"),
    Symbol("L$(period)p$lretstring", lbenchstring, period, "yrtype"))
end

#consider making a pooled version of this
function persistenceplot(data::NCCSData, period::Int, endyear::Int, lbenchmark::NSymbol;
  outputPath::String = GRAPH_PATH, pooled::Bool=true,
  graphsuffix::String = GRAPH_SUFFIX)::Nothing

  #containers for the grpahs
  local plotNames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()
  local subdf::SubDataFrame

  #declare and build the symbols used to get the fields
  local percentilesym::Symbol
  local percentilesymtype::Symbol
  local Lpercentilesym::Symbol
  local Lpercentilesymtype::Symbol

  local percentilebenchsym::Symbol
  local percentilebenchsymtype::Symbol
  local Lpercentilebenchsym::Symbol
  local Lpercentilebenchsymtype::Symbol

  #get the field symbols given the period and benchmark
  (percentilesym, Lpercentilesym) = buildpercentilesym(period, lbenchmark)
  (percentilesymtype, Lpercentilesymtype) = buildpercentilesymtype(period, lbenchmark)

  local dotsize::Vector = [0.4mm]

  if pooled
    subdf = view(data.df, data.df[:fisyr] .≤ endyear, :)
    subdf = view(subdf, ((endyear .- subdf[:fisyr]) .% period) .== 0, :)
  else
    subdf = view(data.df, data.df[:fisyr] .== endyear, :)
  end
  subdf = view(subdf, (!ismissing).(subdf[percentilesym]), :)
  subdf = view(subdf, (!ismissing).(subdf[Lpercentilesym]), :)

  push!(plotNames, "$(something(lbenchmark,""))$(period)yrPerf$(endyear)$(pooled ? "-pooled" : "")")
  #=push!(plots, plot(subdf, x=Lpercentilesym, y=percentilesym,
    #size=:lnetassets,
    Guide.xlabel("$(period)yr return in $(endyear-period)"),
    Guide.ylabel("$(period)yr return in $(endyear)"),
    Geom.point, style(highlight_width=0mm, point_size=0.3mm),
    Coord.cartesian(xmin = 0.0, xmax=1.0, ymin=0.0, ymax = 1.0)))=#
  push!(plots, plot(subdf, x=Lpercentilesym, y=percentilesym,
    #size=:lnetassets,
    Guide.xlabel("Lagged $period year cumulative return percentiles"),
    Guide.ylabel("$period year cumulative return percentiles"),
    Geom.hexbin(xbincount=100, ybincount=100),
    #style(highlight_width=0mm, point_size=0.3mm),
    Coord.cartesian(xmin = 0.0, xmax=1.0, ymin=0.0, ymax = 1.0)))

  #NOTE: Extract tails here!!!!!!!!!

  for i ∈ 1:length(plots) #write the graphs
    draw(PDF("$outputPath\\$(plotNames[i])_$(graphsuffix).pdf", 9inch, 7inch), plots[i])
    #savefig(plots[i], "$outputPath\\$(plotNames[i])$(graphsuffix).$ext")
    #println("Graph $(plotNames[i]) written.")
  end

  return nothing
end

#table is complex enough to require own type
struct PersistenceTable
  persistencecolumns::Vector{Vector{Union{MFloat64,Int}}}
  rowindex::Dict
  colindex::Dict

  typefields::Vector{Symbol}
  periods::Vector{Int}
  statfields::Vector{Symbol}
  pooled::Bool
end

function PersistenceTable(data::NCCSData, endyear::Int, lbenchmark::NSymbol;
  periods::Vector{Int}=PERSISTENCE_PERIODS,  pooled::Bool = true,
  cutoffs=PERSISTENCE_CUTOFFS)::PersistenceTable

  #variables for views on the data
  local subdf::SubDataFrame
  local ssubdf::SubDataFrame

  #variables to hold the field names
  local percentilesym::Symbol
  local Lpercentilesym::Symbol

  #it is helpful to structure the rows and columns via an index
  local typefields::Vector{Symbol} = [:pc, :pf, :co, :all]

  #form a dictionary with the column index by type then period
  local colindex::Dict =
    Dict((typefields[i],periods[j])=>j + (i-1)*length(periods)
      for i ∈ 1:length(typefields)
      for j ∈ 1:length(periods))
  local colindexkeys::Vector{Tuple{Symbol, Int}} = collect(keys(colindex))
  local ncols::Int = length(colindexkeys)

  #now build the row index
  local cutofffields::Vector{Symbol} = (f::Float64 -> Symbol("p",Int(round(f*100)))).(cutoffs)
  local statfields::Vector{Symbol} = [cutofffields; :ppersistent; :N]
  local rowindex::Dict = Dict(statfields[r]=>r for r ∈ 1:length(statfields))
  local nrows = length(statfields)

  #preallocate the columns
  persistencecolumns::Vector{Vector{Union{MFloat64,Int}}} =
    (i->Vector{Union{MFloat64,Int}}(undef, nrows)).(1:ncols)

  #holds the top and bottom cutoff for each interval
  cutoffintervals::Vector{Tuple{Float64,Float64}} =
    (i::Int->i==1 ? (0.0, cutoffs[i]) : (cutoffs[i-1], cutoffs[i])).(1:length(cutoffs))

  if !pooled #the not-pooled scenario is more precise but much less powerful
    subdf = view(data.df, data.df[:fisyr] .== endyear,:) #ensure same dataset for all periods
    periodyears = (y->[endyear]).(1:(length(periods))) # same year for all periods
    for p ∈ periods
      (percentilesym, Lpercentilesym) = buildpercentilesym(p, lbenchmark)

      subdf = view(subdf, (!ismissing).(subdf[percentilesym]), :)
      subdf = view(subdf, (!ismissing).(subdf[Lpercentilesym]), :)
    end
  end

  #filter down to a single type
  for p::Int ∈ periods
    if pooled
      subdf = view(data.df, data.df[:fisyr] .≤ endyear, :)
      subdf = view(subdf, ((endyear .- subdf[:fisyr]) .% p) .== 0, :)

      (percentilesym, Lpercentilesym) = buildpercentilesym(p, lbenchmark)
      subdf = view(subdf, (!ismissing).(subdf[percentilesym]), :)
      subdf = view(subdf, (!ismissing).(subdf[Lpercentilesym]), :)
    end

    for t::NSymbol ∈ typefields

      #build the symbols referencing the appropriate data column
      if t==:all
        ssubdf = subdf
        (percentilesym, Lpercentilesym) = buildpercentilesym(p, lbenchmark)
      else
        ssubdf = view(subdf, subdf[:charitytype] .== t, :)
        (percentilesym, Lpercentilesym) = buildpercentilesymtype(p, lbenchmark)
      end

      N::Int = size(ssubdf,1)

      #will hold the total number in each quantile
      NStarts::Vector{Int} = Vector{Int}(undef, length(cutoffs))

      #will hold the total number of perisstent in eachF quantile
      NPersistents::Vector{Int} = Vector{Int}(undef, length(cutoffs))

      for c::Int ∈ 1:length(cutoffs)

        #number of NPs in the first period
        NStarts[c] = sum((ssubdf[Lpercentilesym] .≥ cutoffintervals[c][1]) .&
          (ssubdf[Lpercentilesym] .< cutoffintervals[c][2]))

        #number of NPs in the first and seocnd period
        NPersistents[c] = sum((ssubdf[Lpercentilesym] .≥ cutoffintervals[c][1]) .&
          (ssubdf[Lpercentilesym] .< cutoffintervals[c][2]) .&
          (ssubdf[percentilesym] .≥ cutoffintervals[c][1]) .&
          (ssubdf[percentilesym] .< cutoffintervals[c][2])
          )

        #fill in the values
        persistencecolumns[colindex[t,p]][rowindex[cutofffields[c]]] = NPersistents[c]/NStarts[c]
      end

      #sum the persistence values
      persistencecolumns[colindex[t,p]][rowindex[:ppersistent]] = sum(NPersistents)/N
      persistencecolumns[colindex[t,p]][rowindex[:N]] = N

    end
  end

  return PersistenceTable(persistencecolumns,
    rowindex,
    colindex,
    typefields,
    periods,
    statfields,
    pooled)

end

function persistencetable(P::PersistenceTable, endyear::Int, lbenchmark::NSymbol;
  outputpath::String = TABLE_PATH,
  decimals::Int = DECIMALS_RETURN,
  tablesuffix::String = TABLE_SUFFIX
  )::Nothing

  local lbenchmarkstring::String = isnothing(lbenchmark) ? "" : "$lbenchmark"

  local ncols::Int = length(P.persistencecolumns) #number of columns excluding the year
  local nrows::Int = length(P.persistencecolumns[1])

  #collect the keys reverse the dicitonaries for notational convenience
  local colindexkeys::Vector{Tuple{Symbol,Int}} = collect(keys(P.colindex))
  local rowindexkeys::Vector{Symbol} = collect(keys(P.rowindex))

  local revcolindex::Dict = Dict(c=>P.colindex[colindexkeys[c]]  for c ∈ 1:ncols)
  local revrowindex::Dict = Dict(r=>P.rowindex[rowindexkeys[r]]  for r ∈ 1:nrows)

  #column headers
  local periodfields::Vector{Symbol} = (i::Int->Symbol("R", i,"y")).(P.periods)
  local typenames = (s->s==:all ? "all non-profits" : "$(charitytype[s])").(P.typefields)
  #typenames = (s->"$s")(typenames)

  #include make the column names
  colNames::Vector{Vector{String}} = Vector{Vector{String}}()
  push!(colNames, typenames)
  push!(colNames, [((s->string(s)).(periodfields) for t ∈ 1:length(P.typefields))...;])

  widthColNames::Vector{Vector{Int}} = Vector{Vector{Int}}()
  push!(widthColNames, [length(periodfields) for i ∈ 1:length(colNames[1])])
  push!(widthColNames, ones(Int, ncols))

  #make the row names
  local descRowNames::Vector{String} =
    (s->s == :ppersistent ? "\\% Persistant" : string(s)).(P.statfields)

  #need to print the descriptive rows
  descContent::Vector{Vector{String}} = (
    (i::Int)->Vector{String}(undef, ncols)).(1:nrows)

  #NOTE: P has the data in COLUMNS while descContet absorbs data in ROWS (flipped indices)
  for r::Int ∈ 1:nrows
    for c::Int ∈ 1:ncols

      #Ns are ints
      if P.statfields[r] == :N
        #lookup the collumn number from the column index for the pth column symbol
        descContent[r][c] = num2Str(P.persistencecolumns[c][r], Ints=true)
      else
        descContent[r][c] = num2Str(P.persistencecolumns[c][r], decimals)
      end
    end
  end

  outtabletitle::String = (isnothing(lbenchmark) ? "Persistence of Non-Profit Returns $(endyear)"
    : "Persistence of Non-Profit Returns Net of $(lbenchmark) $(endyear)")
  outtable::String = texTable(titleCaption="Persistence of Non-Profit Returns $(endyear)",
      caption="""See Tex File""", #caption
      colNames= colNames, #colNames
      descRowNames=descRowNames, #descRowNames
      descContent = descContent, #descContent
      nakedTable=true,
      widthColNames=widthColNames,
      alignmentstring=string("l ", join((w->join([ " | "; ("r" for i ∈ 1:w)...;])
        ).(widthColNames[1])))
    )

  #periodstring::String = string((s->"-$s").(columnsymbols)...)
  pooledstring::String = P.pooled ? :"(pooled)" : ""
  outName::String = (
    "persistence$(lbenchmarkstring)$(endyear)$pooledstring-$(length(P.periods))_$tablesuffix.tex")
  writeNakedTable(outtable, path=outputpath,
    outName = outName)
  return nothing
end

function persistencetablepresentation(P::PersistenceTable, endyear::Int, lbenchmark::NSymbol;
  outputpath::String = TABLE_PATH,
  decimals::Int = DECIMALS_RETURN,
  tablesuffix::String = TABLE_SUFFIX
  )::Nothing

  local lbenchmarkstring::String = isnothing(lbenchmark) ? "" : "$lbenchmark"

  local ncols::Int = length(P.persistencecolumns) #number of columns excluding the year
  local nrows::Int = length(P.persistencecolumns[1])

  #collect the keys reverse the dicitonaries for notational convenience
  local colindexkeys::Vector{Tuple{Symbol,Int}} = collect(keys(P.colindex))
  local rowindexkeys::Vector{Symbol} = collect(keys(P.rowindex))

  local revcolindex::Dict = Dict(c=>P.colindex[colindexkeys[c]]  for c ∈ 1:ncols)
  local revrowindex::Dict = Dict(r=>P.rowindex[rowindexkeys[r]]  for r ∈ 1:nrows)

  #column headers
  local periodfields::Vector{Symbol} = (i::Int->Symbol("R", i,"y")).(P.periods)
  local typenames = (s->s==:all ? "all non-profits" : "$(charitytype[s])").(P.typefields)

  #include make the column names
  colNames::Vector{Vector{String}} = Vector{Vector{String}}()
  push!(colNames, (string).(periodfields))

  #build the row names
  local blocknames::Vector{String} =
    (s::Symbol->s == :ppersistent ? "\\% Persistant" : string(s)).(P.statfields)
  local descRowNames = [(n->blocknames).(1:length(typenames))...;]
  #println("$descRowNames")

  #push!(colNames, [((s->string(s)).(periodfields) for t ∈ 1:length(P.typefields))...;])

  #need to print the descriptive rows
  descContent::Vector{Vector{String}} = (
    (i::Int)->Vector{String}(undef, length(periodfields))).(1:length(descRowNames))

  #NOTE: Indices are hacked together due to the change in table shape
  for s::Int ∈ 1:length(P.statfields)
    for p::Int ∈ 1:length(periodfields)
      for t::Int ∈ 1:length(typenames)

        outrow::Int = (t-1)*(length(P.statfields)) + s
        incol::Int = (t-1)*(length(periodfields)) + p

        if s == 1 && p==1
          descRowNames[outrow] = "\\textit{$(typenames[t])}\\\\$(descRowNames[outrow])"
          if t ≠ 1
            descRowNames[outrow] = "\\midrule $(descRowNames[outrow])"
          end
        end

        #Ns are ints
        if P.statfields[s] == :N
          #lookup the collumn number from the column index for the pth column symbol
          descContent[outrow][p] = num2Str(P.persistencecolumns[incol][s], Ints=true)
        else
          descContent[outrow][p] = num2Str(P.persistencecolumns[incol][s], decimals)
        end
      end
    end
  end

  outtabletitle::String = (isnothing(lbenchmark) ? "Persistence of Non-Profit Returns $(endyear)"
    : "Persistence of Non-Profit Returns Net of $(lbenchmark) $(endyear)")
  outtable::String = texTable(titleCaption="Persistence of Non-Profit Returns $(endyear)",
      caption="""See Tex File""", #caption
      colNames= colNames, #colNames
      descRowNames=descRowNames, #descRowNames
      descContent = descContent, #descContent
      nakedTable=true#,
      #widthColNames=widthColNames,
      #alignmentstring=string("l ", join((w->join([ " | "; ("r" for i ∈ 1:w)...;])
      # ).(widthColNames[1])))
    )
  pooledstring::String = P.pooled ? :"(pooled)" : ""
  outName::String = ("prespersistence$(lbenchmarkstring)$(endyear)$pooledstring-" *
    "$(length(P.periods))_$tablesuffix.tex")

  writeNakedTable(outtable, path=outputpath,
    outName = outName)
end

function persistencetables(data::NCCSData; persistenceyears::Vector{Int} = PERSISTENCE_YEARS,
  lbenchmarks::Vector{NSymbol} = LBENCHMARKS)::Nothing

  for s ∈ lbenchmarks
    for y ∈ persistenceyears

      P::PersistenceTable = PersistenceTable(data, y, s, pooled=false)
      persistencetable(P, y, s)

      P = PersistenceTable(data, y, s, pooled=true)
      persistencetable(P, y, s)
      persistencetablepresentation(P, y, s)
    end
  end

  return nothing
end

#creates all of the persistance plots
function persistenceplots(data::NCCSData,
  persistenceperiods::Vector{Int} = PERSISTENCE_PERIODS,
  persistenceyears::Vector{Int} = PERSISTENCE_YEARS,
  lbenchmarks::Vector{NSymbol} = LBENCHMARKS)

  for p::Int ∈ persistenceperiods, y::Int ∈ persistenceyears, s ∈ lbenchmarks
    persistenceplot(data, p, y, s)
  end

  println("completed persistence plots")
end

#writes out the data as a csv if desired
function writepercentiles(data::NCCSData;
  lbenchmark::NSymbol = nothing #=:bmFF3=#, period::Int = 4, endyear::Int = 2015,
  datapath::String = DATA_PATH, pooled::Bool=true,
  datasuffix::String = datasuffix)::Nothing

  local percentilesym::Symbol
  local Lpercentilesym::Symbol

  #get the percentile symbols
  (percentilesym, Lpercentilesym) = buildpercentilesym(period, lbenchmark)
  subdf = view(data.df, data.df[:fisyr] .≤ endyear, :)

  if pooled
    subdf = view(subdf, (((endyear .- subdf[:fisyr]) .% period) .== 0), :)
  else
    subdf = view(subdf, endyear .== subdf[:fisyr], :)
  end

  subdf = view(subdf, ((s::MFloat64, Ls::MFloat64)->
      (!ismissing(s)) && (!ismissing(Ls))).(
      subdf[percentilesym], subdf[Lpercentilesym]),:)

  #build the path and file name
  path =
    "$datapath\\percentilesout_$(percentilesym)$(pooled ? "-pooled" : "")_$(datasuffix).csv"
  CSV.write(path, subdf)

  return nothing
end

function analyzepersistence(data::NCCSData; percentiles2files::Bool = false)
  #make the graphs
  persistenceplots(data)
  #make the table
  persistencetables(data)
  percentiles2files && writepercentiles(data::NCCSData)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
## template

## binary search
function search_left_border(left::Int, right::Int, key::Function)
    ## e.g. search_left_border(1, 10, >(3)) | returns 4
    while left <= right
        mid = left + (right - left) >> 1
        if key(mid) ## lies in the right part
            right = mid - 1
        else
            left = mid + 1
        end
    end
    left
end

## segment tree
mutable struct SegmentTree{T<:Real}
    n::Int
    tree::Vector{T}

    function SegmentTree(nums::Vector{T}) where T
        n = length(nums)
        tree = append!(zeros(Int, n - 1), nums)
        for i in (n - 1):-1:1
            tree[i] = tree[i << 1] + tree[i << 1 | 1]
        end
        new{T}(n, tree)
    end
end

function update!(ST::SegmentTree, ind::Int, val::Int)::Nothing
    tree = ST.tree
    ind += ST.n - 1
    delta = val - tree[ind]
    while ind > 0
        tree[ind] += delta
        ind >>= 1
    end
end

function sum_range(ST::SegmentTree, left::Int, right::Int)::Int
    left += ST.n - 1
    right += ST.n - 1
    res, tree = 0, ST.tree
    while left <= right
        if isodd(left) ## right child
            res += tree[left] ## record outside value
            left += 1
        end
        left >>= 1
        if iseven(right) ## left child
            res += tree[right] ## record outside value
            right -= 1
        end
        right >>= 1
    end
    res
end

## merge sort
function merge_sort(arr::AbstractArray)::AbstractArray
    (n = length(arr)) == 1 && return arr
    mid = n >> 1
    lpart = merge_sort(@view(arr[1:mid]))
    rpart = merge_sort(@view(arr[mid + 1:end]))
    return merge_sorted(lpart, rpart)
end

function merge_sorted(lpart::AbstractArray{T}, rpart::AbstractArray{T})::AbstractArray{T} where T
    l1, l2 = length(lpart), length(rpart)
    combined = Array{T}(undef, l1 + l2)
    p1 = p2 = 1
    for pos in eachindex(combined)
        if p2 > l2 || p1 <= l1 && lpart[p1] <= rpart[p2]
            combined[pos] = lpart[p1]
            p1 += 1
        else
            combined[pos] = rpart[p2]
            p2 += 1
        end
    end
    combined
end

## 二叉树遍历
## 先序遍历 + 中序遍历
traversal(::Nothing) = Int[]
function traversal(root::TreeNode)::Vector{Int}
    res, stack = Int[], []
    while !isempty(stack) || !isnothing(root)
        while !isnothing(root)
            push!(stack, root)
            ## preorder traversal
            root = root.left
        end
        root = pop!(stack)
        ## inorder traversal
        root = root.right
    end
    res
end

ext: jl


text:
##########
##########
# Get all initial populations before perturbation
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-08-09
##########
##########

# set up =======================================================================
import Pkg
Pkg.instantiate()

using DrWatson
@quickactivate "trait-based-rewiring"

using DifferentialEquations, DelimitedFiles

# logger set up
#using Logging: global_logger
#using TerminalLoggers: TerminalLogger
#global_logger(TerminalLogger())

include(srcdir("00_species_list_creation.jl"))
include(srcdir("01_eqn_component_functions.jl"))
include(srcdir("03_initial_sims_functions.jl"))
include(srcdir("04_perturbation_sims_functions.jl"))

# set all positional values
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12

# read in files
all_species_lists = readdlm(datadir("intermediate-files/species_lists_for_simulations.csv"), ',', Any)

# little line to set the network id (this happens from the shell line)
network_id = parse(Int64, ARGS[1])

# set the species list to be worked on
species_list = all_species_lists[all_species_lists[:,network_id_pos] .== network_id,:]
network = create_interaction_matrix(species_list)

# set up the simulation
#  allocate u0
u0 = Float64[]
for i in 1:size(species_list, 1)
    append!(u0, species_list[i, biomass_pos])
end

# set the extinction callback
# set vector callback
extinction_callback =
    VectorContinuousCallback(extinction_threshold_pre_perturb,
                            extinction_affect_pre_perturb!,
                            size(species_list, 1),
                            save_positions = ntuple(x->true,
                                                    size(species_list, 1)))

# put species list and network in same object
p = [species_list, network]

# set  tspan
tspan = (0.0, 1000.0)

# set problem
prob_setup = ODEProblem(biomass_basic,
                            u0,
                            tspan,
                            p)
@time sol = solve(prob_setup,
                    Tsit5(),
                    callback = extinction_callback,
                    save_times = [0.0, 1000.0])

# put the values back in the species list
species_list[:, biomass_pos] = sol.u[end]

# change extinction values
for i in 1:size(species_list, 1)
    if species_list[i, biomass_pos] == 0
        species_list[i, alive_pos] = 0
    end
end

name = lpad(network_id, 4, "0")

# bind back the network_id and write the species list
writedlm(datadir("./intermediate-files/species-lists-pre-perturbation/species_list_$network_id.csv"), species_list, ",")

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
abstract type Closure end
@inline closure_relation(γ, br, u, r) = exp.(-u .+ γ .+ br) .- γ .- 1.0

struct HypernettedChain <: Closure end
@inline closure_relation(γ, ::HypernettedChain, u, r) = closure_relation(γ, 0.0, u, r)

struct PercusYevick <: Closure end
@inline function closure_relation(γ, ::PercusYevick, u, r)
    return closure_relation(γ, -γ .+ log.(1.0 .+ γ), u, r)
end

struct ModifiedVerlet{R<:Real} <: Closure
    α::R
    β::R
    δ::R
end
ModifiedVerlet() = ModifiedVerlet{Float64}(-0.5, 0.0, 0.8)
ModifiedVerlet(a::T, b::T, c::T) where {T<:Real} = ModifiedVerlet{T}(a, b, c)

function closure_relation(γ, c::ModifiedVerlet, u, r)
    new_br = @. (c.α * γ^2) * (1.0 + c.β * γ)
    @. new_br /= (1.0 + c.δ * abs(γ))

    return closure_relation(γ, new_br, u, r)
end

struct MeanSpherical <: Closure end

function closure_relation(γ, ::MeanSpherical, u, r)
    cr = similar(r)
    less_one = r .< 1.0
    larger_one = r .> 1.0

    @. cr[less_one] = -1.0 - γ[less_one]
    @. cr[larger_one] = -u[larger_one]

    return cr
end

struct SoftMeanSpherical <: Closure end
function _softmeanspherical!(cr, u1, u2, gamma, idxs)
    @. cr[idxs] = exp(-u1[idxs]) * (1.0 + gamma[idxs] - u2[idxs])
    @. cr[idxs] += -gamma[idxs] - 1.0

    return nothing
end

function closure_relation(γ, ::SoftMeanSpherical, u, r)
    cr = similar(r)
    u₁ = similar(u)
    u₂ = similar(u)

    less_one = r .< 1.0
    larger_one = r .>= 1.0

    u₁[less_one] .= u[less_one]
    u₁[larger_one] .= 0.0
    u₂[less_one] .= 0.0
    u₂[larger_one] .= u[larger_one]

    _softmeanspherical!(cr, u₁, u₂, γ, less_one)
    _softmeanspherical!(cr, u₁, u₂, γ, larger_one)

    return cr
end

struct HMSA{T<:Real} <: Closure
    α::T
end

function _fmix!(fm, alpha, r, idxs)
    @. fm[idxs] = 1.0 - exp(-alpha * r[idxs])

    return nothing
end

function _hmsa(cr, u1, u2, gamma, fm, idxs)
    @. cr[idxs] = exp(-u1[idxs])
    @. cr[idxs] *= 1.0 + ((exp((gamma[idxs] - u2[idxs]) * fm[idxs]) - 1.0) / fm[idxs])
    @. cr[idxs] += -gamma[idxs] - 1.0

    return nothing
end

function closure_relation(γ, cls::HMSA, u, r)
    cr = similar(r)
    f_mix = similar(r)
    u₁ = similar(u)
    u₂ = similar(u)

    less_one = r .< 1.0
    larger_one = r .>= 1.0

    u₁[less_one] .= u[less_one]
    u₁[larger_one] .= 0.0
    u₂[less_one] .= 0.0
    u₂[larger_one] .= u[larger_one]

    _fmix!(f_mix, cls.α, r, less_one)
    _fmix!(f_mix, cls.α, r, larger_one)

    _hmsa(cr, u₁, u₂, γ, f_mix, less_one)
    _hmsa(cr, u₁, u₂, γ, f_mix, larger_one)

    return cr
end


text:
# make sure to include this file in src/GeophysOperator.jl
module FooPhysics
# completely independent of transD_GP
# you can use whatever packages you want here

# this is probably where you'll store all variables needed for your physics
# e.g., wavenumbers and so on
mutable struct Foo
    x         :: Int
    somefield :: Vector{Float64}
end

# this is what will return the computation from our physics engine
function returnphysics!(F::Foo, m::AbstractArray)
    n = length(F.somefield)
    F.somefield .= F.x*rand(n) .+ sum(m[:]) # pre-allocated assignment
end

end

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
include("jim.jl")

export jim


text:
macro git_cmd(ex)
    return _git_cmd(ex)
end

function git(f::Function;
             adjust_PATH::Bool = true,
             adjust_LIBPATH::Bool = true)
    git_path, env_mapping = _env_mapping(; adjust_PATH = adjust_PATH,
                                           adjust_LIBPATH = adjust_LIBPATH)
    return withenv(env_mapping...) do
        return f(git_path)
    end
end


text:
#####################################################
# Download data from UF Sparse Matrix Collection
#####################################################

using Base.Filesystem

# collect the keys from local database (MATRIXDICT or USERMATRIXDICT)
# provide a numerical id counting from 1 for either database.
function insertlocal(db::MatrixDatabase, T::Type{<:GeneratedMatrixData},
                     ldb::Dict{<:AbstractString,Function})

    cnt = 0
    ks = sort!(collect(keys(ldb)))
    for k in ks
        cnt += 1
        push!(db, T(k, cnt, ldb[k]))
    end
    cnt
end

dbpath(db::MatrixDatabase) = abspath(data_dir(), "db.data")
function readdb(db::MatrixDatabase)
    @info("reading database")
    cachedb = dbpath(db)
    open(cachedb, "r") do io
        dbx = deserialize(io)
        merge!(db.data, dbx.data)
        merge!(db.aliases, dbx.aliases)
    end
end

function writedb(db::MatrixDatabase)
    cachedb = dbpath(db)
    open(cachedb, "w") do io
        serialize(io, db)
    end
end

"""
    MatrixDepot.downloadindices(db)
download index files and store matrix data in `db`.
additionally generate aliases for local and loaded matrices.
Serialize db object in file `db.data`.
If a file `db.data` is present in the data directory, deserialize
the data instead of downloading and collecting data.
"""
function downloadindices(db::MatrixDatabase; ignoredb=false)
    added = 0
    if !ignoredb && isfile(dbpath(db))
        try
            readdb(db)
        catch ex
            println(ex)
            @warn("recreating database file")
            _downloadindices(db)
            added = 1
        end
    else
        @info("creating database file")
        _downloadindices(db)
        added = 1
    end
    @info("adding metadata...")
    added += addmetadata!(db)
    @info("adding svd data...")
    added += addsvd!(db)
    added += insertlocal(db, GeneratedMatrixData{:B}, MATRIXDICT)
    added += insertlocal(db, GeneratedMatrixData{:U}, USERMATRIXDICT)
    if added > 0
        @info("writing database")
        writedb(db)
    end
    nothing
end

function _downloadindices(db::MatrixDatabase)
    @info("reading index files")
    empty!(db)

    try
        readindex(preferred(SSRemoteType), db)
        readindex(preferred(MMRemoteType), db)
    finally
        for data in values(db.data)
            db.aliases[aliasname(data)] = data.name
        end
    end
    nothing
end

"""
    MatrixDepot.update()
update database index from the websites
"""
function update(db::MatrixDatabase=MATRIX_DB)
    cachedb = abspath(data_dir(), "db.data")
    isfile(cachedb) && rm(cachedb)
    uf_matrices = localindex(preferred(TURemoteType))
    isfile(uf_matrices) && rm(uf_matrices)
    mm_matrices = localindex(preferred(MMRemoteType))
    isfile(mm_matrices) && rm(mm_matrices)
    downloadindices(db, ignoredb=true)
end

function gunzip(fname)
    endswith(fname, ".gz") || error("gunzip: $fname: unknown suffix")

    destname = rsplit(fname, ".gz", limit=2)[1]
    BUFFSIZE = 1000000
    open(destname, "w") do f
        open(GzipDecompressorStream, fname) do g
            buffer = read(g, BUFFSIZE)
            while length(buffer) > 0
                write(f, buffer)
                buffer = read(g, BUFFSIZE)
            end
        end
    end
    destname
end

"""
    loadmatrix(data::RemoteMatrixData)

Download the files backing the data from a remote repository. That is currently
the TAMU site for the UFl collection and the NIST site for the MatrixMarket
collection. The files are uncompressed and un-tar-ed if necessary.
The data files containing the matrix data have to be in MatrixMarket format in
both cases. Note, that some of the files of the MM collection are not available
in MatrixMarket format. An error message results, if tried to load them.
"""
function loadmatrix(data::RemoteMatrixData)
    file = matrixfile(data)
    if isfile(file)
        addmetadata!(data)
        return 0
    end
    dirfn = localfile(data)
    dir = dirname(localdir(data))
    url = redirect(dataurl(data))

    isdir(dir) || mkpath(dir)
    wdir = pwd()
    try
        @info("downloading: $url")
        downloadfile(url, dirfn)
        tarfile = gunzip(dirfn)
        cd(dir)
        rfile = relpath(string(tarfile))
        if endswith(tarfile, ".tar")
            run(`tar -xf $rfile`)
            rm(tarfile; force=true)
        end
    finally
        cd(wdir)
        rm(dirfn, force=true)
    end
    addmetadata!(data)
    1
end
loadmatrix(data::GeneratedMatrixData) = 0

"""
    loadinfo(data::RemoteDate)
Download the first part of the data file. Stop reading, as soon as the initial
comment and the size values of the main matrix have been finished. Store this in
a file with extension `.info` in the same directory, where the `.mtx` file is.
If the complete file is already available, the download is not performed, because
the head of the `.mtx` file contains the same lines.
"""
function loadinfo(data::RemoteMatrixData)
    filemtx = matrixfile(data)
    file = matrixinfofile(data)
    if isfile(filemtx) || isfile(file)
        return 0
    end
    url = redirect(dataurl(data))
    pipe = downloadpipeline(url)
    out = IOBuffer()
    s = try
        @info("downloading head of $url")
        open(pipe, "r") do io
            skip = 0
            while ( s = readline(io) ) != ""
                skip = s[1] == '%'  || isempty(strip(s)) ? 0 : skip + 1
                skip <= 1 && println(out, s)
                if skip == 1 && length(split(s)) == 3
                    break
                end
            end;
            close(io)
        end
        String(take!(out))
    catch ex
        ex isa InterruptException && rethrow()
        String(take!(out))
    finally
        close(out)
    end
    if !isempty(s)
        mkpath(dirname(file))
        write(file, s)
        addmetadata!(data)
        1
    else
        0
    end
end
loadinfo(data::MatrixData) = 0

"""
    downloadpipeline(url)

Set up a command pipeline (external processes to download and expand data)
"""
function downloadpipeline(url::AbstractString)
    urls = rsplit(url, '.', limit=3)
    cmd = []
    push!(cmd, downloadcommand(url))
    if urls[end] == "gz"
        push!(cmd, `gzip -dc`)
        resize!(urls, length(urls)-1)
    end
    if urls[end] == "tar"
        push!(cmd, `tar -xOf -`)
    end
    pipeline(cmd...)
end

function downloadcommand(url::AbstractString, filename::AbstractString="-")
    `sh -c 'curl "'$url'" -Lso "'$filename'"'`
end

function data_warn(data::RemoteMatrixData, dn, i1, i2)
    @warn "$(data.name): header $dn = '$i1' file $dn = '$i2' $(data.properties[])"
    i1
end

function extremnnz(data::RemoteMatrixData, m, n, k)
    if issymmetric(data) || ishermitian(data) || isskew(data)
        if m != n
            @warn "$(data.name) needs to be quadratic but is ($mx$n)"
        end
        isskew(data) ? (2k, 2k) : (2k - max(m, n), 2k)
    else
        (k, k)
    end
end
function extremnnz(data::RemoteMatrixData)
    extremnnz(data, data.m, data.n, data.dnz)
end

function addmetadata!(db::MatrixDatabase)
    added = 0
    for data in values(db.data)
        if isunloaded(data)
            addmetadata!(data)
            added += isloaded(data)
        end
    end
    added
end

addmetadata!(::MatrixData) = 0
function addmetadata!(data::RemoteMatrixData)
    file = matrixfile(data)
    dir = dirname(file)
    empty!(data.metadata)
    isdir(dir) || return 0
    base = basename(file)
    name, ext = rsplit(base, '.', limit=2)
    filtop(x) = x == base || (startswith(x, string(name, '_')) && !endswith(x, "SVD.mat"))
    append!(data.metadata, filter(filtop, readdir(dir)))

    if !isfile(file)
        file = matrixinfofile(data)
    end
    if (finfo = mmreadheader(file)) !== nothing
        data.properties[] = MMProperties(MATRIX, finfo[:format], finfo[:field], finfo[:symmetry])
        m, n, k = finfo[:m], finfo[:n], get(finfo, :nz, 0)
        n1, n2 = extremnnz(data, m, n, k)
        hdr = data.header
        hdr.m = hdr.m in (0, m) ? m : data_warn(data, "m", hdr.m, m)
        hdr.n = hdr.n in (0, n) ? n : data_warn(data, "n", hdr.n, n)
        hdr.nnz = hdr.nnz == 0  ? n2 : hdr.nnz <= n2 ? hdr.nnz :
                                        data_warn(data, "nnz", hdr.nnz, k)
        hdr.dnz = k
        date = get(finfo, :date, 0)
        hdr.date = hdr.date in (0, date) ? date : data_warn(data, "date", hdr.date, date)
        kind = get(finfo, :kind, "")
        if kind != ""
            if hdr.kind == "" || hdr.kind == "other problem"
                hdr.kind = kind
            else
                if standardkind(kind) != standardkind(hdr.kind)
                    data_warn(data, "kind", hdr.kind, kind)
                end
            end
        end
        for s in (:title, :author, :ed, :fields, :notes)
            val = get(finfo, s, "")
            val != "" && setfield!(hdr, s, val)
        end
    end
    return 1
end

standardkind(s::AbstractString) = lowercase(replace(s, '-' => ' '))

function addsvd!(db::MatrixDatabase)
    added = 0
    for data in values(db.data)
        if data isa RemoteMatrixData && !issvdok(data)
            added += addsvd!(data)
        end
    end
    added
end

issvdok(data::RemoteMatrixData) = data.svdok
issvdok(::MatrixData) = false

"""
    downloadfile(url, out)

Copy file from remote or local url. Works around julia Downloads #69 and #36
"""
function downloadfile(url::AbstractString, out::AbstractString)
    run(downloadcommand(url, out))
    nothing
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# 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.
##########################################################################

# copy from src to dst
function Base.copy!(dst::Location, src::Location)
	dst.x = src.x
	dst.y = src.y
end

# return true if locations have same x and y values, false otherwise
function Base.:(==)(loc1::Location, loc2::Location)
	return (loc1.x == loc2.x && loc1.y == loc2.y)
end

# return squared distance between two locations
function squareDist(map::Map, loc1::Location, loc2::Location)
	return ((loc1.x - loc2.x)*map.xScale)^2 + ((loc1.y - loc2.y)*map.yScale)^2
end

# return distance between two locations
function normDist(map::Map, loc1::Location, loc2::Location)
	return sqrt(squareDist(map, loc1, loc2))
end

# calculate off road travel time between two locations for a given travel mode
function offRoadTravelTime(travelMode::TravelMode, map::Map, loc1::Location, loc2::Location)
	return normDist(map, loc1, loc2) / travelMode.offRoadSpeed
end

# calculate off road travel time for given distance
function offRoadTravelTime(travelMode::TravelMode, dist::Float)
	return dist / travelMode.offRoadSpeed
end

# # calculate off road travel time between two locations for a given speed
# function offRoadTravelTime(map::Map, speed::Float, loc1::Location, loc2::Location)
	# return normDist(map, loc1, loc2) / speed
# end

# # calculate off road travel time for a given distance and speed
# function offRoadTravelTime(speed::Float, dist::Float)
	# return dist / speed
# end

# function linearInterp(x1::Float, x::Float, x2::Float, y1::Float, y2::Float)
	# # given x1 <= x <= x2, return matching linear interpolation for y values
	# # @assert(x1 <= x && x <= x2 && x1 < x2) # have removed this line, for speed reasons
	# return (x - x1) / (x2 - x1) * (y2 - y1) + y1
# end

function linearInterpLocation(startLoc::Location, endLoc::Location, startTime::Float, endTime::Float, currentTime::Float)
	@assert(startTime <= currentTime && currentTime <= endTime && startTime < endTime)
	# will assume constant travel time
	p = (currentTime - startTime) / (endTime - startTime) # portion of travel complete
	currentLoc = Location()
	currentLoc.x = (1-p)*startLoc.x + p*endLoc.x
	currentLoc.y = (1-p)*startLoc.y + p*endLoc.y

	return currentLoc
end

# return random location from map
# use uniform distribution
# allow fractional trimming of border
function randLocation(map::Map; trim::Float = 0.0, rng::AbstractRNG = GLOBAL_RNG)
	@assert(trim >= 0 && trim <= 1)
	r = rand(rng, 2) .* (1-trim) .+ trim/2
	location = Location()
	location.x = map.xMin + map.xRange * r[1]
	location.y = map.yMin + map.yRange * r[2]
	return location
end

# return map with borders trimmed by fraction
function trimmedMap(map::Map, trim::Float = 0.0)
	@assert(trim >= 0 && trim <= 1)
	mapTrimmed = deepcopy(map)
	mapTrimmed.xMin = map.xMin + map.xRange * trim / 2
	mapTrimmed.xMax = map.xMax - map.xRange * trim / 2
	mapTrimmed.yMin = map.yMin + map.yRange * trim / 2
	mapTrimmed.yMax = map.yMax - map.yRange * trim / 2
	mapTrimmed.xRange = mapTrimmed.xMax - mapTrimmed.xMin
	mapTrimmed.yRange = mapTrimmed.yMax - mapTrimmed.yMin
	return mapTrimmed
end


text:
using DTMarkovChain

runs = [[(1,2),(1,1),(1,2),(2,1),(1,2),(2,2)],
  [(2,1),(1,1),(1,2),(2,2),(2,2),(1,1),(1,2),(2,1),(1,1)]]
states = [(1,1),(1,2),(2,1),(2,2)]
mc = fit(DTFiniteMarkovChain,runs,states)


text: Dict{String,Any}("valid key" => Dict{String,Any}())

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using Clang

#error("$(@__FILE__) is broken until Clang.jl supports wrapping C++")

# ROCBLAS_HEADERS are those headers to be wrapped.
const ROCBLAS_INCLUDE = joinpath(@__DIR__, "include") |> normpath

rocblas_headers_base = joinpath(@__DIR__, "include", "rocblas")
const ROCBLAS_HEADERS = String[]
push!(ROCBLAS_HEADERS, joinpath(rocblas_headers_base, "julia_rocblas.h"))

wc = init(; headers = ROCBLAS_HEADERS,
            output_file = joinpath(@__DIR__, "rocblas_interface.jl"),
            common_file = joinpath(@__DIR__, "rocblas_types.jl"),
            clang_includes = vcat(ROCBLAS_INCLUDE, CLANG_INCLUDE),
            #clang_args = ["-I", joinpath(ROCBLAS_INCLUDE, "..")],
            header_wrapped = (root, current)->root == current,
            header_library = x->"librocblas",
            clang_diagnostics = true,
            )

run(wc)

# Remove useless ROCBLAS_API macros
#run(Cmd(["sed", "-i", "/^.*ROCBLAS_API.*\$/d", "rocblas_types.jl"]))


text:
using .Networks

using Printf
using Dates
using Parameters
using DataStructures

@enum GeneratorState begin
    ON
    OFF
end

function Base.float(generator_state::GeneratorState)
    if generator_state == ON
        return 1.
    elseif generator_state == OFF
        return 0.
    else
        throw( error("Unknown GeneratorState value : ", generator_state) )
    end
end

function Base.parse(::Type{GeneratorState}, str::String)
    if lowercase(str) == "on"
        return ON
    elseif  lowercase(str) == "off"
        return OFF
    else
        throw( error("Unable to convert `", str, "` to a GeneratorState") )
    end
end

function Base.parse(::Type{GeneratorState}, val::Float64)
    if val > 1e-09
        return ON
    else
        return OFF
    end
end

##########################################
## Firmness
##########################################

@enum DecisionFirmness begin
    DECIDED # a firm decision was already taken
    TO_DECIDE # a firm decision needs to be taken
    FREE # by scenario decisions
end

@with_kw_noshow struct Firmness
    #gen,ts
    commitment::SortedDict{String, SortedDict{Dates.DateTime, DecisionFirmness} } =
        SortedDict{String, SortedDict{Dates.DateTime, DecisionFirmness} }()
    power_level::SortedDict{String, SortedDict{Dates.DateTime, DecisionFirmness} } =
        SortedDict{String, SortedDict{Dates.DateTime, DecisionFirmness} }()
end

Base.:(==)(a::Firmness, b::Firmness) = a.commitment==b.commitment && a.power_level==b.power_level

function get_commitment_firmness(firmness::Firmness)
    return firmness.commitment
end

function get_power_level_firmness(firmness::Firmness)
    return firmness.power_level
end

function set_commitment_firmness!(firmness::Firmness, gen_id::String, ts, decision_firmness::DecisionFirmness)
    get!(firmness.commitment, gen_id, SortedDict{Dates.DateTime, DecisionFirmness}())
    firmness.commitment[gen_id][ts] = decision_firmness
end

function set_power_level_firmness!(firmness::Firmness, gen_id::String, ts, decision_firmness::DecisionFirmness)
    get!(firmness.power_level, gen_id, SortedDict{Dates.DateTime, DecisionFirmness}())
    firmness.power_level[gen_id][ts] = decision_firmness
end

function get_commitment_firmness(firmness::Firmness, gen_id::String)
    if haskey(firmness.commitment,gen_id)
        return firmness.commitment[gen_id]
    else
        return missing
    end
end

function get_commitment_firmness(firmness::Firmness, gen_id::String, ts::Dates.DateTime)
    if haskey(firmness.commitment,gen_id)
        if haskey(firmness.commitment[gen_id], ts)
            return firmness.commitment[gen_id][ts]
        else
            return missing
        end
    else
        return missing
    end
end

function get_power_level_firmness(firmness::Firmness, gen_id::String)
    if haskey(firmness.power_level,gen_id)
        return firmness.power_level[gen_id]
    else
        return missing
    end
end

function get_power_level_firmness(firmness::Firmness, gen_id::String, ts::Dates.DateTime)
    if haskey(firmness.power_level,gen_id)
        if haskey(firmness.power_level[gen_id], ts)
            return firmness.power_level[gen_id][ts]
        else
            return missing
        end
    else
        return missing
    end
end

function Base.show(io::IO, firmness::Firmness)
    println(io, "commitment :")
    pretty_print(io, firmness.commitment)
    println(io, "power_level :")
    pretty_print(io, firmness.power_level)
end

##########################################
## Schedule
##########################################

struct GeneratorSchedule
    gen_id::String
    commitment::SortedDict{Dates.DateTime, UncertainValue{GeneratorState}}
    production::SortedDict{Dates.DateTime, UncertainValue{Float64}}
end
function GeneratorSchedule(gen_id::String)
    return  GeneratorSchedule(gen_id,
                            SortedDict{Dates.DateTime, UncertainValue{GeneratorState}}(),
                            SortedDict{Dates.DateTime, UncertainValue{Float64}}())
end

function get_commitment(generator_schedule::GeneratorSchedule)
    return generator_schedule.commitment
end

function get_production(generator_schedule::GeneratorSchedule)
    return generator_schedule.production
end

mutable struct Schedule <: AbstractSchedule
    decider_type::DeciderType
    decision_time::Dates.DateTime
    #gen_id => ts => uncertainValue(s)
    generator_schedules::SortedDict{String, GeneratorSchedule }

    #bus,ts,s
    loss_of_load_by_bus::SortedDict{Tuple{String,DateTime,String}, Float64 }
    #limitable_gen_id,ts,s
    capping::SortedDict{Tuple{String,DateTime,String}, Float64 }
end
function Schedule(decider_type, ech)
    return Schedule(decider_type, ech,
                    SortedDict{String, GeneratorSchedule}(),
                    SortedDict{Tuple{String,DateTime,String}, Float64 }(),
                    SortedDict{Tuple{String,DateTime,String}, Float64 }()
                    )
end
function Schedule(decider_type, ech, generator_schedules)
    return Schedule(decider_type, ech, generator_schedules,
                    SortedDict{Tuple{String,DateTime,String}, Float64 }(),
                    SortedDict{Tuple{String,DateTime,String}, Float64 }()
                    )
end

function init!(schedule::Schedule, network::Networks.Network,
            target_timepoints::Vector{Dates.DateTime}, scenarios::Vector{String})
    empty!(schedule.generator_schedules)
    for generator_l in Networks.get_generators(network)
        gen_id = Networks.get_id(generator_l)
        schedule.generator_schedules[gen_id] = GeneratorSchedule(gen_id)
        for ts in target_timepoints

            schedule.generator_schedules[gen_id].production[ts] = UncertainValue{Float64}(scenarios)

            #commitment values are only defined for generators that have a pmin > 0
            if (Networks.get_p_min(generator_l) > 1e-09)
                schedule.generator_schedules[gen_id].commitment[ts] = UncertainValue{GeneratorState}(scenarios)
            end
        end
    end
end

function get_sub_schedule(schedule::Schedule, gen_id::String)::Union{GeneratorSchedule,Missing}
    if haskey(schedule.generator_schedules, gen_id)
        return schedule.generator_schedules[gen_id]
    else
        throw( error("no generator schedule was initialized for generator ", gen_id))
        #return missing
    end
end

function get_prod_uncertain_value(sub_schedule::GeneratorSchedule, ts::Dates.DateTime)::Union{UncertainValue{Float64},Missing}
    if haskey(sub_schedule.production, ts)
        return sub_schedule.production[ts]
    else
        @warn("power level schedule is not defined for generator ", sub_schedule.gen_id)
        return missing
    end
end
function get_prod_uncertain_value(schedule::Schedule, gen_id::String, ts::Dates.DateTime)::Union{UncertainValue{Float64},Missing}
    return get_prod_uncertain_value(schedule.generator_schedules[gen_id], ts)
end

function get_commitment_uncertain_value(sub_schedule::GeneratorSchedule, ts::Dates.DateTime)::Union{UncertainValue{GeneratorState},Missing}
    if haskey(sub_schedule.commitment, ts)
        return sub_schedule.commitment[ts]
    else
        # @warn("commitment schedule is not defined for generator ", sub_schedule.gen_id)
        return missing
    end
end
function get_commitment_uncertain_value(schedule::Schedule, gen_id::String, ts::Dates.DateTime)::Union{UncertainValue{GeneratorState},Missing}
    return get_commitment_uncertain_value(schedule.generator_schedules[gen_id], ts)
end

function get_prod_value(sub_schedule::GeneratorSchedule, ts::Dates.DateTime)::Union{Float64, Missing}
    uncertain_value = get_prod_uncertain_value(sub_schedule, ts)
    if ismissing(uncertain_value)
        return missing
    end
    return get_value(uncertain_value)
end
function get_prod_value(schedule::Schedule, gen_id::String, ts::Dates.DateTime)::Union{Float64, Missing}
    return get_prod_value(schedule.generator_schedules[gen_id], ts)
end

function safeget_prod_value(schedule::Schedule, gen_id::String, ts::Dates.DateTime)::Float64
    return safeget_prod_value(schedule.generator_schedules[gen_id], ts)
end
function safeget_prod_value(sub_schedule::GeneratorSchedule, ts::Dates.DateTime)::Float64
    prod = get_prod_value(sub_schedule, ts)
    if ismissing(prod)
        msg = @sprintf("Missing production value for ts=%s in schedule %s",
                        ts, sub_schedule.gen_id)
        throw( error(msg) )
    else
        return prod
    end
end

function get_commitment_value(sub_schedule::GeneratorSchedule, ts::Dates.DateTime)::Union{GeneratorState, Missing}
    uncertain_value = get_commitment_uncertain_value(sub_schedule, ts)
    if ismissing(uncertain_value)
        return missing
    end
    return get_value(uncertain_value)
end
function get_commitment_value(schedule::Schedule, gen_id::String, ts::Dates.DateTime)::Union{GeneratorState, Missing}
    return get_commitment_value(schedule.generator_schedules[gen_id], ts)
end

function safeget_commitment_value(schedule::Schedule, gen_id::String, ts::Dates.DateTime)::GeneratorState
    return safeget_commitment_value(schedule.generator_schedules[gen_id], ts)
end
function safeget_commitment_value(sub_schedule::GeneratorSchedule, ts::Dates.DateTime)::GeneratorState
    commitment = get_commitment_value(sub_schedule, ts)
    if ismissing(commitment)
        msg = @sprintf("Missing commitment value for ts=%s in schedule %s",
                        ts, sub_schedule.gen_id)
        throw( error(msg) )
    else
        return commitment
    end
end

function get_prod_value(sub_schedule::GeneratorSchedule, ts::Dates.DateTime, scenario::String)::Union{Float64, Missing}
    uncertain_value = get_prod_uncertain_value(sub_schedule, ts)
    if ismissing(uncertain_value)
        return missing
    end
    return get_value(uncertain_value, scenario)
end
function get_prod_value(schedule::Schedule, gen_id::String, ts::Dates.DateTime, scenario::String)::Union{Float64, Missing}
    return get_prod_value(schedule.generator_schedules[gen_id], ts, scenario)
end
function safeget_prod_value(schedule::Schedule, gen_id::String, ts::Dates.DateTime, scenario::String)::Union{Float64, Missing}
    prod = get_prod_value(schedule.generator_schedules[gen_id], ts, scenario)
    if ismissing(prod)
        msg = @sprintf("Missing production value in schedule for (gen_id=%s,ts=%s,s=%s)",
                        gen_id,ts,scenario)
        throw( error(msg) )
    else
        return prod
    end
end
function safeget_prod_value(sub_schedule::GeneratorSchedule, ts::Dates.DateTime, scenario::String)::Float64
    prod = get_prod_value(sub_schedule, ts, scenario)
    if ismissing(prod)
        msg = @sprintf("Missing production value for (ts=%s,s=%s) in schedule %s",
                        ts, scenario, sub_schedule.gen_id)
        throw( error(msg) )
    else
        return prod
    end
end

function get_commitment_value(sub_schedule::GeneratorSchedule, ts::Dates.DateTime, scenario::String)::Union{GeneratorState, Missing}
    uncertain_value = get_commitment_uncertain_value(sub_schedule, ts)
    if ismissing(uncertain_value)
        return missing
    end
    return get_value(uncertain_value, scenario)
end
function get_commitment_value(schedule::Schedule, gen_id::String, ts::Dates.DateTime, scenario::String)::Union{GeneratorState, Missing}
    return get_commitment_value(schedule.generator_schedules[gen_id], ts, scenario)
end

function get_capping(schedule::Schedule, gen_id::String, ts::Dates.DateTime, scenario)::Union{Float64, Missing}
    key_l = (gen_id, ts, scenario)
    if !haskey(schedule.capping, key_l)
        return missing
    else
        return schedule.capping[key_l]
    end
end
function safeget_capping(schedule::Schedule, gen_id::String, ts::Dates.DateTime, scenario)::Float64
    capping = get_capping(schedule, gen_id, ts, scenario)
    if ismissing(capping)
        msg = @sprintf("Missing capping value for (gen_id=%s,ts=%s,s=%s) in schedule",
                        gen_id, ts, scenario)
        throw( error(msg) )
    else
        return capping
    end
end

function get_loss_of_load(schedule::Schedule, bus_id, ts::Dates.DateTime, scenario)::Union{Float64, Missing}
    key_l = (bus_id, ts, scenario)
    if !haskey(schedule.loss_of_load_by_bus, key_l)
        return missing
    else
        return schedule.loss_of_load_by_bus[key_l]
    end
end
function safeget_loss_of_load(schedule::Schedule, bus_id, ts::Dates.DateTime, scenario)::Float64
    loss_of_load = get_loss_of_load(schedule, bus_id, ts, scenario)
    if ismissing(loss_of_load)
        msg = @sprintf("Missing loss_of_load value for (gen_id=%s,ts=%s,s=%s) in schedule",
                        bus_id, ts, scenario)
        throw( error(msg) )
    else
        return loss_of_load
    end
end

function set_prod_value!(sub_schedule::GeneratorSchedule, ts::Dates.DateTime, scenario::String, value::Float64)
    uncertain_value = get_prod_uncertain_value(sub_schedule, ts)
    if ismissing(uncertain_value)
        msg = @sprintf("cannot set uninitialized schedule production value for generator `%s` at timestep %s and scenario %s",
                        sub_schedule.gen_id, ts, scenario)
        throw( error(msg) )
    end
    set_value!(uncertain_value, scenario, value)
end
function set_prod_value!(schedule, gen_id::String, ts::Dates.DateTime, scenario::String, value::Float64)
    return set_prod_value!(schedule.generator_schedules[gen_id], ts, scenario, value)
end

function set_commitment_value!(sub_schedule::GeneratorSchedule, ts::Dates.DateTime, scenario::String, value::GeneratorState)
    uncertain_value = get_commitment_uncertain_value(sub_schedule, ts)
    if ismissing(uncertain_value)
        msg = @sprintf("cannot set uninitialized schedule commitment value for generator `%s` at timestep %s and scenario %s",
                        sub_schedule.gen_id, ts, scenario)
        throw( error(msg) )
    end
    set_value!(uncertain_value, scenario, value)
end
function set_commitment_value!(schedule, gen_id::String, ts::Dates.DateTime, scenario::String, value::GeneratorState)
    return set_commitment_value!(schedule.generator_schedules[gen_id], ts, scenario, value)
end

function set_prod_definitive_value!(sub_schedule::GeneratorSchedule, ts::Dates.DateTime, value::Float64)
    uncertain_value = get_prod_uncertain_value(sub_schedule, ts)
    if ismissing(uncertain_value)
        msg = @sprintf("cannot set uninitialized schedule production value for generator `%s` at timestep %s",
                        sub_schedule.gen_id, ts)
        throw( error(msg) )
    end
    safeset_definitive_value!(uncertain_value, value)
end
function set_prod_definitive_value!(schedule, gen_id::String, ts::Dates.DateTime, value::Float64)
    return set_prod_definitive_value!(schedule.generator_schedules[gen_id], ts, value)
end

function set_commitment_definitive_value!(sub_schedule::GeneratorSchedule, ts::Dates.DateTime, value::GeneratorState)
    uncertain_value = get_commitment_uncertain_value(sub_schedule, ts)
    if ismissing(uncertain_value)
        msg = @sprintf("cannot set uninitialized schedule commitment value for generator `%s` at timestep %s",
                        sub_schedule.gen_id, ts)
        throw( error(msg) )
    end
    set_definitive_value!(uncertain_value, value)
end
function set_commitment_definitive_value!(schedule, gen_id::String, ts::Dates.DateTime, value::GeneratorState)
    return set_commitment_definitive_value!(schedule.generator_schedules[gen_id], ts, value)
end

function set_capping_value!(schedule::Schedule, gen_id::String, ts::Dates.DateTime, scenario::String, value::Float64)
    schedule.capping[gen_id, ts, scenario] = value
end
function set_loss_of_load_value!(schedule::Schedule, bus_id::String, ts::Dates.DateTime, scenario::String, value::Float64)
    schedule.loss_of_load_by_bus[bus_id, ts, scenario] = value
end

function reset_loss_of_load_by_bus!(schedule::Schedule)
    empty!(schedule.loss_of_load_by_bus)
end
function reset_capping!(schedule::Schedule)
    empty!(schedule.capping)
end

function get_commitment_sub_schedule(schedule::Schedule, gen_id::String)
    return get_commitment(get_sub_schedule(schedule, gen_id))
end

function get_production_sub_schedule(schedule::Schedule, gen_id::String)
    return get_production(get_sub_schedule(schedule, gen_id))
end

function Base.show(io::IO, gen_schedule::GeneratorSchedule)
    @printf("generator %s:\n", gen_schedule.gen_id)
    println("commitment:")
    pretty_print(io, gen_schedule.commitment)
    println("production:")
    pretty_print(io, gen_schedule.production)
end
function Base.show(io::IO, schedule::Schedule)
    @printf("schedule decided by %s at %s:\n", schedule.decider_type, schedule.decision_time)
    pretty_print(io, schedule.generator_schedules)
    println("capping:")
    pretty_print(io, schedule.capping)
    println("loss_of_load_by_bus:")
    pretty_print(io, schedule.loss_of_load_by_bus)
end

#######################
# Helpers
########################

function get_target_timepoints(gen_schedule::GeneratorSchedule)
    return union(keys(gen_schedule.commitment), keys(gen_schedule.production))
end

function get_start_value(preceding_state, current_state)
    if preceding_state==OFF && current_state==ON
        return 1
    else
        return 0
    end
end

function get_start_value(generator_reference_schedule, ts, preceding_ts, generator_initial_state)
    current_state = safeget_commitment_value(generator_reference_schedule, ts)
    preceding_state = ( isnothing(preceding_ts) ? generator_initial_state :
                            safeget_commitment_value(generator_reference_schedule, preceding_ts) )
    return get_start_value(preceding_state, current_state)
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using Fortran
using InteractiveUtils
using Test

@testset "Simple test" begin
    fun = FFunction(:add, FInteger, [Var(:x, FInteger), Var(:y, FInteger)], Var[],
                    Block([Assign(Var(:var"add.result", FInteger),
                                  Call(Fortran.add_integer, FExpr[Var(:x, FInteger), Var(:y, FInteger)]))]))

    cfun = clean_code(feval(fun))
    println("Compiled function:")
    println(cfun)
    eval(cfun)

    arg1 = FInteger(2)
    arg2 = FInteger(3)
    println("code_warntype:")
    @code_warntype add(arg1, arg2)
    println("code_llvm:")
    @code_llvm add(arg1, arg2)
    println("code_native:")
    @code_native add(arg1, arg2)
    println("result:")
    res = add(arg1, arg2)
    println(res)
    @test res ≡ FInteger(5)
end

@testset "Do loop" begin
    fun = FFunction(:loop, FInteger, [Var(:imin, FInteger), Var(:imax, FInteger)], [Var(:count, FInteger), Var(:i, FInteger)],
                    Block([Assign(Var(:count, FInteger), Const(FInteger(0))),
                           Do(Var(:i, FInteger), Var(:imin, FInteger), Var(:imax, FInteger), Const(FInteger(1)),
                              Block([Assign(Var(:count, FInteger),
                                            Call(Fortran.add_integer, FExpr[Var(:count, FInteger), Const(FInteger(1))]))])),
                           Assign(Var(:var"loop.result", FInteger), Var(:count, Integer))]))
    cfun = clean_code(feval(fun))
    println("Compiled function:")
    println(cfun)
    eval(cfun)
    println("code_warntype:")
    @code_warntype loop(FInteger(2), FInteger(5))
    println("code_llvm:")
    @code_llvm loop(FInteger(2), FInteger(5))
    println("code_native:")
    @code_native loop(FInteger(2), FInteger(5))
    println("result:")
    res = loop(FInteger(2), FInteger(5))
    println(res)
    @test res ≡ FInteger(5 - 2 + 1)
end

@testset "Do-while loop" begin
    fun = FFunction(:loop2, FInteger, [Var(:imin, FInteger), Var(:imax, FInteger)], [Var(:count, FInteger), Var(:i, FInteger)],
                    Block([Assign(Var(:count, FInteger), Const(FInteger(0))), Assign(Var(:i, FInteger), Var(:imin, FInteger)),
                           DoWhile(Call(Fortran.le_integer, FExpr[Var(:i, FInteger), Var(:imax, FInteger)]),
                                   Block([Assign(Var(:count, FInteger),
                                                 Call(Fortran.add_integer, FExpr[Var(:count, FInteger), Const(FInteger(1))])),
                                          Assign(Var(:i, FInteger),
                                                 Call(Fortran.add_integer, FExpr[Var(:i, FInteger), Const(FInteger(1))]))])),
                           Assign(Var(:var"loop2.result", FInteger), Var(:count, Integer))]))
    cfun = clean_code(feval(fun))
    println("Compiled function:")
    println(cfun)
    eval(cfun)
    println("code_warntype:")
    @code_warntype loop2(FInteger(2), FInteger(5))
    println("code_llvm:")
    @code_llvm loop2(FInteger(2), FInteger(5))
    println("code_native:")
    @code_native loop2(FInteger(2), FInteger(5))
    println("result:")
    res = loop2(FInteger(2), FInteger(5))
    println(res)
    @test res ≡ FInteger(5 - 2 + 1)
end

@testset "If statment loop" begin
    fun = FFunction(:cond, FInteger, [Var(:c, FLogical), Var(:x, FInteger), Var(:y, FInteger)], Var[],
                    If([Var(:c, FLogical) => Assign(Var(:var"cond.result", FInteger), Var(:x, FInteger))],
                       Assign(Var(:var"cond.result", FInteger), Var(:y, FInteger))))
    cfun = clean_code(feval(fun))
    println("Compiled function:")
    println(cfun)
    eval(cfun)
    println("code_warntype:")
    @code_warntype cond(FLogical(true), FInteger(2), FInteger(5))
    println("code_llvm:")
    @code_llvm cond(FLogical(true), FInteger(2), FInteger(5))
    println("code_native:")
    @code_native cond(FLogical(true), FInteger(2), FInteger(5))
    println("result:")
    res = cond(FLogical(true), FInteger(2), FInteger(5))
    println(res)
    @test cond(FLogical(true), FInteger(2), FInteger(5)) ≡ FInteger(2)
    @test cond(FLogical(false), FInteger(2), FInteger(5)) ≡ FInteger(5)
end


text:
using POMDPs
using POMDPModelTools
using POMDPXFiles
using POMDPModels
using Test

# write your own tests here

file_name = "tiger_test.pomdpx"
pomdp = TigerPOMDP()
pomdpx = POMDPXFile(file_name)
write(pomdp, pomdpx)
av, aa = read_pomdp("mypolicy.policy")

@test av ≈ [-81.5975 3.01448 24.6954 28.4025 19.3711; 28.4025 24.6954 3.01452 -81.5975 19.3711]
@test aa == [1,0,0,2,0]

@testset begin "a, sp observation warning"
    struct BadObsPOMDP <: POMDP{Int,Int,Int} end
    POMDPs.states(m::BadObsPOMDP) = 1:2
    POMDPs.actions(m::BadObsPOMDP) = 1:2
    POMDPs.observations(m::BadObsPOMDP) = 1:2
    POMDPs.transition(m::BadObsPOMDP, s, a) = Deterministic(clamp(s+a, 1, 2))
    POMDPs.reward(m::BadObsPOMDP, s, a, o) = s
    POMDPs.observation(m::BadObsPOMDP, s, a, sp) = Deterministic(sp)
    POMDPs.discount(m::BadObsPOMDP) = 0.99
    POMDPs.initialstate(m::BadObsPOMDP) = Deterministic(1)
    POMDPs.stateindex(m::BadObsPOMDP, s) = s
    POMDPs.actionindex(m::BadObsPOMDP, s) = s
    POMDPs.obsindex(m::BadObsPOMDP, s) = s

    @test_throws MethodError cd(mktempdir()) do
        write(BadObsPOMDP(), POMDPXFile("bad_obs_test.pomdpx"))
    end

    POMDPs.observation(m::BadObsPOMDP, a, sp) = Deterministic(1)
    cd(mktempdir()) do
        write(BadObsPOMDP(), POMDPXFile("bad_obs_test.pomdpx"))
    end
end


text:
using Colors
using Gadfly
using Distributions

"""
Plot a Gamma Distribution
"""
function drawGammaDistribution()
    as = [1. 1.5 2.]
    bs = ones(1, length(as))
    xs = linspace(0.1, 7, 40)

    colors = ["orange", "red", "purple"]
    titles = [string("a=", as[i], ",b=", bs[i]) for i=1:length(as)]

    γs = [Gamma(as[i], bs[i]) for i=1:length(as)]
    ps = [pdf(γ, x) for x=xs, γ=γs]

    μ = as./bs

    ls = [layer(x=xs, y=ps[:,i], Geom.line,
                Theme(default_color=parse(Colorant, colors[i]), line_width=3px),
                xintercept=[μ[i]], Geom.vline(color=parse(Colorant, colors[i]), size=3px))
          for i=1:length(as)]

    pp = plot(Guide.manual_color_key("Gamma parameters", titles, colors), ls...)
end


text:
"""
    opt,sol,data = cs_tssos_first(supp::Vector{Vector{Vector{Vector{UInt16}}}}, coe::Vector{Vector{ComplexF64}},
    n, d; numeq=0, foc=100, CS="MF", minimize=false, assign="first", TS="block", merge=false, md=3, solver="Mosek",
    QUIET=false, solve=true, MomentOne=true)

Compute the first step of the CS-TSSOS hierarchy for constrained complex polynomial optimization with
relaxation order `d`. Here the complex polynomial optimization problem is defined by `supp` and `coe`,
corresponding to the supports and coeffients of `pop` respectively.

# Arguments
- `supp`: the supports of the complex polynomial optimization problem.
- `coe`: the coeffients of the complex polynomial optimization problem.
- `d`: the relaxation order of the moment-SOHS hierarchy.
- `numeq`: the number of equality constraints.
"""
function cs_tssos_first(supp::Vector{Vector{Vector{Vector{UInt16}}}}, coe, n, d; numeq=0, foc=100,
    CS="MF", minimize=false, assign="first", TS="block", merge=false, md=3, solver="Mosek",
    QUIET=false, solve=true, tune=false, solution=false, ipart=true, MomentOne=false, Mommat=false)
    println("***************************TSSOS***************************")
    println("TSSOS is launching...")
    m = length(supp)-1
    ind = [supp[1][i][1]<=supp[1][i][2] for i=1:length(supp[1])]
    supp[1] = supp[1][ind]
    coe[1] = coe[1][ind]
    supp[1],coe[1] = resort(supp[1], coe[1])
    dg = zeros(Int, m)
    for i = 1:m
        dg[i] = maximum([length(supp[i+1][j][1]) + length(supp[i+1][j][2]) for j=1:length(supp[i+1])])
    end
    time = @elapsed begin
    cliques,cql,cliquesize = clique_decomp(n, m, dg, supp, order=d, alg=CS, minimize=minimize)
    end
    if CS != false && QUIET == false
        mc = maximum(cliquesize)
        println("Obtained the variable cliques in $time seconds. The maximal size of cliques is $mc.")
    end
    J,ncc = assign_constraint(m, supp, cliques, cql, cliquesize, assign=assign)
    rlorder = init_order(dg, J, cliquesize, cql, foc=foc, order=d)
    if TS != false && QUIET == false
        println("Starting to compute the block structure...")
    end
    time = @elapsed begin
    blocks,cl,blocksize,sb,numb,basis,status = get_cblocks_mix(dg, J, rlorder, m, supp, cliques, cql, cliquesize, TS=TS, merge=merge, md=md)
    end
    if TS != false && QUIET == false
        mb = maximum(maximum.(sb))
        println("Obtained the block structure in $time seconds. The maximal size of blocks is $mb.")
    end
    opt,ksupp,Mmatrix = blockcpop_mix(n, m, supp, coe, basis, cliques, cql, cliquesize, J, ncc, blocks, cl, blocksize, numeq=numeq, QUIET=QUIET, TS=TS, solver=solver, solve=solve, tune=tune, solution=solution, ipart=ipart, MomentOne=MomentOne, Mommat=Mommat)
    data = mcpop_data(n, 0, m, numeq, supp, coe, basis, rlorder, ksupp, cql, cliques, cliquesize, J, ncc, sb, numb, blocks, cl, blocksize, Mmatrix, solver, 1e-4, 1)
    return opt,nothing,data
end

function blockcpop_mix(n, m, supp::Vector{Vector{Vector{Vector{UInt16}}}}, coe, basis, cliques,
    cql, cliquesize, J, ncc, blocks, cl, blocksize; numeq=0, nb=0, QUIET=false, TS="block", solver="Mosek",
    tune=false, solve=true, solution=false, MomentOne=false, ipart=true, Mommat=false)
    tsupp = Vector{Vector{UInt16}}[]
    for i = 1:cql
        for j = 1:cl[i][1], k = 1:blocksize[i][1][j], r = k:blocksize[i][1][j]
            @inbounds bi = [basis[i][1][blocks[i][1][j][k]], basis[i][1][blocks[i][1][j][r]]]
            if bi[1] <= bi[2]
                push!(tsupp, bi)
            else
                push!(tsupp, bi[2:-1:1])
            end
        end
        for (j, w) in enumerate(J[i])
            for l = 1:cl[i][j+1], t = 1:blocksize[i][j+1][l], r = t:blocksize[i][j+1][l], s = 1:length(supp[w+1])
                ind1 = blocks[i][j+1][l][t]
                ind2 = blocks[i][j+1][l][r]
                @inbounds bi = [sadd(basis[i][j+1][ind1], supp[w+1][s][1]), sadd(basis[i][j+1][ind2], supp[w+1][s][2])]
                if bi[1] <= bi[2]
                    push!(tsupp, bi)
                else
                    push!(tsupp, bi[2:-1:1])
                end
            end
        end
    end
    for i ∈ ncc, j = 1:length(supp[i+1])
        if supp[i+1][j][1] <= supp[i+1][j][2]
            push!(tsupp, supp[i+1][j])
        end
    end
    if (MomentOne == true || solution == true) && TS != false
        ksupp = copy(tsupp)
        for i = 1:cql, j = 1:cliquesize[i]
            push!(tsupp, [UInt16[], UInt16[cliques[i][j]]])
            for k = j+1:cliquesize[i]
                push!(tsupp, [UInt16[cliques[i][j]], UInt16[cliques[i][k]]])
            end
        end
    end
    sort!(tsupp)
    unique!(tsupp)
    if (MomentOne == true || solution == true) && TS != false
        sort!(ksupp)
        unique!(ksupp)
    else
        ksupp = tsupp
    end
    objv = nothing
    Mmatrix = nothing
    if solve == true
        ltsupp = length(tsupp)
        if QUIET == false
            println("Assembling the SDP...")
            println("There are $ltsupp affine constraints.")
        end
        if solver == "Mosek"
            model = Model(optimizer_with_attributes(Mosek.Optimizer))
            if tune == true
                set_optimizer_attributes(model,
                "MSK_DPAR_INTPNT_CO_TOL_MU_RED" => 1e-7,
                "MSK_DPAR_INTPNT_CO_TOL_INFEAS" => 1e-7,
                "MSK_DPAR_INTPNT_CO_TOL_REL_GAP" => 1e-7,
                "MSK_DPAR_INTPNT_CO_TOL_DFEAS" => 1e-7,
                "MSK_DPAR_INTPNT_CO_TOL_PFEAS" => 1e-7,
                "MSK_DPAR_INTPNT_CO_TOL_NEAR_REL" => 1e6,
                "MSK_IPAR_BI_IGNORE_NUM_ERROR" => 1,
                "MSK_DPAR_BASIS_TOL_X" => 1e-3,
                "MSK_DPAR_BASIS_TOL_S" => 1e-3,
                "MSK_DPAR_BASIS_REL_TOL_S" => 1e-5)
            end
        elseif solver == "COSMO"
            model = Model(optimizer_with_attributes(COSMO.Optimizer, "max_iter" => 10000))
        elseif solver == "SDPT3"
            model = Model(optimizer_with_attributes(SDPT3.Optimizer))
        else
            @error "The solver is currently not supported!"
            return nothing,nothing,nothing
        end
        set_optimizer_attribute(model, MOI.Silent(), QUIET)
        time = @elapsed begin
        rcons = [AffExpr(0) for i=1:ltsupp]
        if ipart == true
            icons = [AffExpr(0) for i=1:ltsupp]
        end
        for i = 1:cql
            if (MomentOne == true || solution == true) && TS != false
                bs = cliquesize[i]+1
                pos = @variable(model, [1:2*bs, 1:2*bs], PSD)
                for t = 1:bs, r = t:bs
                    @constraint(model, pos[t,r]==pos[t+bs,r+bs])
                    @constraint(model, pos[r,t+bs]+pos[t,r+bs]==0)
                    if t == 1 && r == 1
                        bi = [UInt16[], UInt16[]]
                    elseif t == 1 && r > 1
                        bi = [UInt16[], UInt16[cliques[i][r-1]]]
                    else
                        bi = [UInt16[cliques[i][t-1]], UInt16[cliques[i][r-1]]]
                    end
                    Locb = bfind(tsupp, ltsupp, bi)
                    if ipart == true
                        @inbounds add_to_expression!(icons[Locb], pos[t+bs,r])
                    end
                    @inbounds add_to_expression!(rcons[Locb], pos[t,r])
                end
            end
            for l = 1:cl[i][1]
                @inbounds bs = blocksize[i][1][l]
                if bs == 1
                    pos = @variable(model, lower_bound=0)
                    @inbounds bi = [basis[i][1][blocks[i][1][l][1]], basis[i][1][blocks[i][1][l][1]]]
                    Locb = bfind(tsupp, ltsupp, bi)
                    @inbounds add_to_expression!(rcons[Locb], pos)
                else
                    if ipart == true
                        pos = @variable(model, [1:2bs, 1:2bs], PSD)
                    else
                        pos = @variable(model, [1:bs, 1:bs], PSD)
                    end
                    for t = 1:bs, r = t:bs
                        if ipart == true
                            @constraint(model, pos[t,r]==pos[t+bs,r+bs])
                            @constraint(model, pos[r,t+bs]+pos[t,r+bs]==0)
                        end
                        @inbounds ind1 = blocks[i][1][l][t]
                        @inbounds ind2 = blocks[i][1][l][r]
                        @inbounds bi = [basis[i][1][ind1], basis[i][1][ind2]]
                        if bi[1] <= bi[2]
                            Locb = bfind(tsupp, ltsupp, bi)
                            if ipart == true
                                @inbounds add_to_expression!(icons[Locb], pos[t+bs,r])
                            end
                        else
                            Locb = bfind(tsupp, ltsupp, bi[2:-1:1])
                            if ipart == true
                                @inbounds add_to_expression!(icons[Locb], -1, pos[t+bs,r])
                            end
                        end
                        @inbounds add_to_expression!(rcons[Locb], pos[t,r])
                    end
                end
            end
        end
        for i in ncc
            if i <= m-numeq
                pos = @variable(model, lower_bound=0)
            else
                pos = @variable(model)
            end
            for j = 1:length(supp[i+1])
                if supp[i+1][j][1] <= supp[i+1][j][2]
                    Locb = bfind(tsupp, ltsupp, supp[i+1][j])
                    if ipart == true
                        @inbounds add_to_expression!(icons[Locb], imag(coe[i+1][j]), pos)
                    end
                    @inbounds add_to_expression!(rcons[Locb], real(coe[i+1][j]), pos)
                end
            end
        end
        for i = 1:cql, (j, w) in enumerate(J[i]), l = 1:cl[i][j+1]
            bs = blocksize[i][j+1][l]
            if bs == 1
                if w <= m-numeq
                    pos = @variable(model, lower_bound=0)
                else
                    pos = @variable(model)
                end
                ind = blocks[i][j+1][l][1]
                for s = 1:length(supp[w+1])
                    @inbounds bi = [sadd(basis[i][j+1][ind], supp[w+1][s][1]), sadd(basis[i][j+1][ind], supp[w+1][s][2])]
                    if bi[1] <= bi[2]
                        Locb = bfind(tsupp,ltsupp,bi)
                        if ipart == true
                            @inbounds add_to_expression!(icons[Locb], imag(coe[w+1][s]), pos)
                        end
                        @inbounds add_to_expression!(rcons[Locb], real(coe[w+1][s]), pos)
                    end
                end
            else
                if w <= m-numeq
                    if ipart == true
                        pos = @variable(model, [1:2bs, 1:2bs], PSD)
                    else
                        pos = @variable(model, [1:bs, 1:bs], PSD)
                    end
                else
                    if ipart == true
                        pos = @variable(model, [1:2bs, 1:2bs], Symmetric)
                    else
                        pos = @variable(model, [1:bs, 1:bs], Symmetric)
                    end
                end
                if ipart == true
                    for t = 1:bs, r = t:bs
                        @constraint(model, pos[t,r]==pos[t+bs,r+bs])
                        @constraint(model, pos[r,t+bs]+pos[t,r+bs]==0)
                    end
                end
                for t = 1:bs, r = 1:bs
                    ind1 = blocks[i][j+1][l][t]
                    ind2 = blocks[i][j+1][l][r]
                    for s = 1:length(supp[w+1])
                        @inbounds bi = [sadd(basis[i][j+1][ind1], supp[w+1][s][1]), sadd(basis[i][j+1][ind2], supp[w+1][s][2])]
                        if bi[1] <= bi[2]
                            Locb = bfind(tsupp,ltsupp,bi)
                            @inbounds add_to_expression!(rcons[Locb], real(coe[w+1][s]), pos[t,r])
                            if ipart == true
                                @inbounds add_to_expression!(icons[Locb], imag(coe[w+1][s]), pos[t,r])
                                @inbounds add_to_expression!(icons[Locb], real(coe[w+1][s]), pos[t+bs,r])
                                @inbounds add_to_expression!(rcons[Locb], -imag(coe[w+1][s]), pos[t+bs,r])
                            end
                        end
                    end
                end
            end
        end
        rbc = zeros(ltsupp)
        if ipart == true
            ibc = zeros(ltsupp)
        end
        for i = 1:length(supp[1])
            Locb = bfind(tsupp, ltsupp, supp[1][i])
            if Locb == 0
               @error "The monomial basis is not enough!"
               return nothing,ksupp,nothing
            else
               rbc[Locb] = real(coe[1][i])
               if ipart == true
                   ibc[Locb] = imag(coe[1][i])
               end
            end
        end
        @variable(model, lower)
        if Mommat == true
            rcons[1] += lower
            @constraint(model, rcon[i=1:ltsupp], rcons[i]==rbc[i])
            if ipart == true
                @constraint(model, icon[i=1:ltsupp], icons[i]==ibc[i])
            end
        else
            @constraint(model, rcons[2:end].==rbc[2:end])
            @constraint(model, rcons[1]+lower==rbc[1])
            if ipart == true
                @constraint(model, icons.==ibc)
            end
        end
        @objective(model, Max, lower)
        end
        if QUIET == false
            println("SDP assembling time: $time seconds.")
            println("Solving the SDP...")
        end
        time = @elapsed begin
        optimize!(model)
        end
        if QUIET == false
            println("SDP solving time: $time seconds.")
        end
        status = termination_status(model)
        objv = objective_value(model)
        if status != MOI.OPTIMAL
           println("termination status: $status")
           status = primal_status(model)
           println("solution status: $status")
        end
        println("optimum = $objv")
        if Mommat == true
            rmeasure = -dual.(rcon)
            imeasure = nothing
            if ipart == true
                imeasure = -dual.(icon)
            end
            Mmatrix = get_cmoment(rmeasure, imeasure, tsupp, cliques, cql, cliquesize, blocks, cl, blocksize, basis, ipart=ipart)
        end
    end
    return objv,ksupp,Mmatrix
end

function get_cblocks_mix(dg, J, rlorder, m, supp::Vector{Vector{Vector{Vector{UInt16}}}}, cliques, cql,
    cliquesize; tsupp=[], basis=[], blocks=[], cl=[], blocksize=[], sb=[], numb=[], TS="block", nb=0, merge=false, md=3)
    if isempty(tsupp)
        blocks = Vector{Vector{Vector{Vector{UInt16}}}}(undef, cql)
        cl = Vector{Vector{UInt16}}(undef, cql)
        blocksize = Vector{Vector{Vector{UInt16}}}(undef, cql)
        sb = Vector{Vector{UInt16}}(undef, cql)
        numb = Vector{Vector{UInt16}}(undef, cql)
        basis = Vector{Vector{Vector{Vector{UInt16}}}}(undef, cql)
        tsupp = copy(supp[1])
        for i = 2:m+1, j = 1:length(supp[i])
            if supp[i][j][1] <= supp[i][j][2]
                push!(tsupp, supp[i][j])
            end
        end
        sort!(tsupp)
        unique!(tsupp)
        flag = 1
    else
        flag = 0
    end
    status = ones(Int, cql)
    for i = 1:cql
        lc = length(J[i])
        nvar = cliquesize[i]
        ind = [issubset(union(tsupp[j][1], tsupp[j][2]), cliques[i]) for j=1:length(tsupp)]
        fsupp = tsupp[ind]
        if flag == 1
            basis[i] = Vector{Vector{Vector{UInt16}}}(undef, lc+1)
            basis[i][1] = get_basis(cliques[i], rlorder[i])
            for s = 1:lc
                basis[i][s+1] = get_basis(cliques[i], rlorder[i]-ceil(Int, dg[J[i][s]]/2))
            end
            blocks[i] = Vector{Vector{Vector{UInt16}}}(undef, lc+1)
            cl[i] = Vector{UInt16}(undef, lc+1)
            blocksize[i] = Vector{Vector{UInt16}}(undef, lc+1)
            sb[i] = Vector{UInt16}(undef, lc+1)
            numb[i] = Vector{UInt16}(undef, lc+1)
            blocks[i],cl[i],blocksize[i],sb[i],numb[i],status[i] = get_cblocks(lc, fsupp, supp[J[i].+1], basis[i], TS=TS, QUIET=true, merge=merge, md=md)
        else
            blocks[i],cl[i],blocksize[i],sb[i],numb[i],status[i] = get_cblocks(lc, fsupp, supp[J[i].+1], basis[i], blocks=blocks[i], cl=cl[i], blocksize=blocksize[i], sb=sb[i], numb=numb[i], TS=TS, QUIET=true, merge=merge, md=md)
        end
    end
    return blocks,cl,blocksize,sb,numb,basis,maximum(status)
end

function assign_constraint(m, supp::Vector{Vector{Vector{Vector{UInt16}}}}, cliques, cql, cliquesize; assign="first")
    J = [UInt32[] for i=1:cql]
    ncc = UInt32[]
    for i = 2:m+1
        rind = copy(supp[i][1][1])
        for j = 2:length(supp[i])
            append!(rind, supp[i][j][1])
        end
        unique!(rind)
        if assign == "first"
            ind = findfirst(k->issubset(rind, cliques[k]), 1:cql)
            if ind != nothing
                push!(J[ind], i-1)
            else
                push!(ncc, i-1)
            end
        else
            temp = UInt32[]
            for j = 1:cql
                if issubset(rind, cliques[j])
                    push!(temp, j)
                end
            end
            if !isempty(temp)
                if assign == "min"
                    push!(J[temp[argmin(cliquesize[temp])]], i-1)
                else
                    push!(J[temp[argmax(cliquesize[temp])]], i-1)
                end
            else
                push!(ncc, i-1)
            end
        end
    end
    return J,ncc
end

function get_basis(var::Vector{UInt16}, d)
    n = length(var)
    lb = binomial(n+d,d)
    basis = Vector{Vector{UInt16}}(undef, lb)
    basis[1] = UInt16[]
    i = 0
    t = 1
    while i < d+1
        t += 1
        if length(basis[t-1])>=i && basis[t-1][end-i+1:end] == var[n]*ones(UInt16, i)
           if i < d
               basis[t] = var[1]*ones(UInt16, i+1)
           end
           i += 1
        else
            j = bfind(var, n, basis[t-1][1])
            basis[t] = copy(basis[t-1])
            ind = findfirst(x->basis[t][x]!=var[j], 1:length(basis[t]))
            if ind == nothing
                ind = length(basis[t])+1
            end
            if j != 1
                basis[t][1:ind-2] = var[1]*ones(UInt16, ind-2)
            end
            basis[t][ind-1] = var[j+1]
        end
    end
    return basis
end

function get_graph(tsupp::Vector{Vector{Vector{UInt16}}}, basis; nb=0)
    lb = length(basis)
    G = SimpleGraph(lb)
    ltsupp = length(tsupp)
    for i = 1:lb, j = i+1:lb
        bi = [basis[i], basis[j]]
        if bi[1] > bi[2]
            bi = bi[2:-1:1]
        end
        if bfind(tsupp, ltsupp, bi)!=0
            add_edge!(G, i, j)
        end
    end
    return G
end

function get_cgraph(tsupp::Vector{Vector{Vector{UInt16}}}, supp, basis; nb=0)
    lb = length(basis)
    ltsupp = length(tsupp)
    G = SimpleGraph(lb)
    for i = 1:lb, j = i+1:lb
        r = 1
        while r <= length(supp)
            bi = [sadd(basis[i], supp[r][1]), sadd(basis[j], supp[r][2])]
            if bi[1] > bi[2]
                bi = bi[2:-1:1]
            end
            if bfind(tsupp, ltsupp, bi)!=0
               break
            else
                r += 1
            end
        end
        if r <= length(supp)
            add_edge!(G, i, j)
        end
    end
    return G
end

function clique_decomp(n, m, dg, supp::Vector{Vector{Vector{Vector{UInt16}}}}; order="min", alg="MF", minimize=false)
    if alg==false
        cliques = [UInt16[i for i=1:n]]
        cql = 1
        cliquesize = [n]
    else
        G = SimpleGraph(n)
        for i = 1:m+1
            if order == "min" || i == 1 || order == ceil(Int, dg[i-1]/2)
                for j = 1:length(supp[i])
                    add_clique!(G, unique([supp[i][j][1];supp[i][j][2]]))
                end
            else
                temp = copy(supp[i][1][1])
                for j = 2:length(supp[i])
                    append!(temp, supp[i][j][1])
                end
                add_clique!(G, unique(temp))
            end
        end
        if alg == "NC"
            cliques,cql,cliquesize = max_cliques(G)
        else
            cliques,cql,cliquesize = chordal_cliques!(G, method=alg, minimize=minimize)
        end
    end
    uc = unique(cliquesize)
    sizes = [sum(cliquesize.== i) for i in uc]
    println("------------------------------------------------------")
    println("The clique sizes of varibles:\n$uc\n$sizes")
    println("------------------------------------------------------")
    return cliques,cql,cliquesize
end

function get_cmoment(rmeasure, imeasure, tsupp, cliques, cql, cliquesize, blocks, cl, blocksize, basis; ipart=true)
    moment = Vector{Vector{Matrix{Union{Float64,ComplexF64}}}}(undef, cql)
    ltsupp = length(tsupp)
    for i = 1:cql
        moment[i] = Vector{Matrix{Union{Float64,ComplexF64}}}(undef, cl[i][1])
        for l = 1:cl[i][1]
            bs = blocksize[i][1][l]
            rtemp = zeros(Float64, bs, bs)
            if ipart == true
                itemp = zeros(Float64, bs, bs)
            end
            for t = 1:bs, r = t:bs
                ind1 = blocks[i][1][l][t]
                ind2 = blocks[i][1][l][r]
                bi = [basis[i][1][ind1], basis[i][1][ind2]]
                if bi[1] <= bi[2]
                    Locb = bfind(tsupp, ltsupp, bi)
                    if ipart == true
                        itemp[t,r] = imeasure[Locb]
                    end
                else
                    Locb = bfind(tsupp, ltsupp, bi[2:-1:1])
                    if ipart == true
                        itemp[t,r] = -imeasure[Locb]
                    end
                end
                rtemp[t,r] = rmeasure[Locb]
            end
            rtemp = (rtemp + rtemp')/2
            if ipart == true
                itemp = (itemp - itemp')/2
                moment[i][l] = rtemp + itemp*im
            else
                moment[i][l] = rtemp
            end
        end
    end
    return moment
end

function resort(supp, coe; field="real")
    nsupp = copy(supp)
    sort!(nsupp)
    unique!(nsupp)
    l = length(nsupp)
    if field == "real"
        ncoe = zeros(l)
    else
        ncoe = zeros(ComplexF64, l)
    end
    for i = 1:length(supp)
        locb = bfind(nsupp, l, supp[i])
        ncoe[locb] += coe[i]
    end
    return nsupp,ncoe
end

# function sign_type(a::Vector{UInt16})
#     st = UInt8[]
#     if length(a) == 1
#         push!(st, a[1])
#     elseif length(a) > 1
#         r = 1
#         for i = 2:length(a)
#             if a[i] == a[i-1]
#                 r += 1
#             else
#                 if isodd(r)
#                     push!(st, a[i-1])
#                 end
#                 r = 1
#             end
#         end
#         if isodd(r)
#             push!(st, a[end])
#         end
#     end
#     return st
# end


text:
######### StanSample example  ###########

using StanSample, Test

bernoulli_model = "
data {
  int<lower=1> N;
  array[N] int<lower=0,upper=1> y;
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = [
  Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
  Dict("N" => 10, "y" => [0, 1, 0, 0, 1, 0, 0, 0, 0, 1]),
  Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 1, 0]),
  Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 1, 0, 0, 1]),
]

sm = SampleModel("bernoulli", bernoulli_model)
rc = stan_sample(sm; data=bernoulli_data, delta=0.85, num_threads=1)
if success(rc)

  # Fetch the same output in the `sdf` ChainDataFrame
  sdf = read_summary(sm)

  @test sdf[sdf.parameters .== :theta, :mean][1] ≈ 0.33 rtol=0.05

end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
# Weave readme
using Pkg
cd("c:/git/DataConvenience/")
Pkg.activate("c:/git/DataConvenience/readme-env")
upcheck()
# Pkg.update()

using Weave

weave("README.jmd", out_path = :pwd, doctype = "github")

if false
    tangle("README.jmd")
end

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
"""
    FiniteTempBasisSet

Type for holding IR bases and sparse-sampling objects.

An object of this type holds IR bases for fermions and bosons
and associated sparse-sampling objects.

# Fields
- basis_f::FiniteTempBasis: Fermion basis
- basis_b::FiniteTempBasis: Boson basis
- beta::Float64: Inverse temperature
- wmax::Float64: Cut-off frequency
- tau::Vector{Float64}: Sampling points in the imaginary-time domain
- wn_f::Vector{Int}: Sampling fermionic frequencies
- wn_b::Vector{Int}: Sampling bosonic frequencies
- smpl_tau_f::TauSampling: Sparse sampling for tau & fermion
- smpl_tau_b::TauSampling: Sparse sampling for tau & boson
- smpl_wn_f::MatsubaraSampling: Sparse sampling for Matsubara frequency & fermion
- smpl_wn_b::MatsubaraSampling: Sparse sampling for Matsubara frequency & boson
- sve_result::Tuple{PiecewiseLegendrePoly,Vector{Float64},PiecewiseLegendrePoly}: Results of SVE
"""
struct FiniteTempBasisSet{
    TSf<:TauSampling,MSf<:MatsubaraSampling,TSb<:TauSampling,MSb<:MatsubaraSampling
}
    basis_f::_DEFAULT_FINITE_TEMP_BASIS
    basis_b::_DEFAULT_FINITE_TEMP_BASIS
    smpl_tau_f::TSf # TauSampling{Float64,_DEFAULT_FINITE_TEMP_BASIS,Float64,Float64}
    smpl_tau_b::TSb # TauSampling{Float64,_DEFAULT_FINITE_TEMP_BASIS,Float64,Float64}
    smpl_wn_f::MSf # MatsubaraSampling{Int,_DEFAULT_FINITE_TEMP_BASIS,ComplexF64,Float64}
    smpl_wn_b::MSb # MatsubaraSampling{Int,_DEFAULT_FINITE_TEMP_BASIS,ComplexF64,Float64}
end

"""
    FiniteTempBasisSet(β, wmax, ε; sve_result=compute_sve(LogisticKernel(β * wmax); ε))

Create basis sets for fermion and boson and
associated sampling objects.
Fermion and bosonic bases are constructed by SVE of the logistic kernel.
"""
function FiniteTempBasisSet(
    β::AbstractFloat, wmax::AbstractFloat, ε;
    sve_result=compute_sve(LogisticKernel(β * wmax); ε),
)
    # Create bases using the given sve results
    basis_f = FiniteTempBasis(fermion, β, wmax, ε; sve_result)
    basis_b = FiniteTempBasis(boson, β, wmax, ε; sve_result)

    return FiniteTempBasisSet(
        basis_f, basis_b,
        TauSampling(basis_f), TauSampling(basis_b),
        MatsubaraSampling(basis_f), MatsubaraSampling(basis_b),
    )
end

getbeta(bset::FiniteTempBasisSet) = getbeta(bset.basis_f)
getwmax(bset::FiniteTempBasisSet) = getwmax(bset.basis_f)

function Base.getproperty(bset::FiniteTempBasisSet, d::Symbol)
    if d === :tau
        # return getfield(bset, :smpl_tau_f).sampling_points
        return bset.smpl_tau_f.sampling_points
    elseif d === :wn_f
        # return getfield(bset, :smpl_wn_f).sampling_points
        return bset.smpl_wn_f.sampling_points
    elseif d === :wn_b
        # return getfield(bset, :smpl_wn_b).sampling_points
        return bset.smpl_wn_b.sampling_points
    elseif d === :sve_result
        # return getfield(bset, :basis_f).sve_result
        return bset.basis_f.sve_result
    else
        return getfield(bset, d)
    end
end

function Base.propertynames(::FiniteTempBasisSet, private::Bool=false)
    return (:tau, :wn_f, :wn_b, :sve_result, propertynames(FiniteTempBasisSet, private)...)
end

function Base.show(io::IO, b::FiniteTempBasisSet)
    return print(io, "FiniteTempBasisSet: beta=$(getbeta(b)), wmax=$(getwmax(b))")
end


text:
function glsphere(;pos=Point3f0(0f0), r=0.1f0, quality=24)
    GLNormalMesh(HyperSphere(pos, r), quality)
end

"""
    draw(cubic::CubicBasis)
    draw(cubics::Vector{CubicBasis})

Shows a 3D representation of the CubicBasis. If a Vector of CUbicBasis is given,
each Basis will be drawn with a different, random color.
"""
draw(cubic::CubicBasis; kwargs...) = draw([cubic], kwargs...)

function draw(
        cubics::Vector{CubicBasis};

        rot_v::Signal{Float32} = Signal(0.05f0),
        trans_v::Signal{Float32} = Signal(4f0),
        cam_position::Vec3f0 = Vec3f0(2),
        lookat::Vec3f0 = Vec3f0(0),
        up::Vec3f0=Vec3f0(0., 0., 1.),

        color::RGBA{Float32} = RGBA{Float32}(0.2, 0.4, 0.8, 1.0),
        primitive::AbstractMesh = glsphere()
    )

    window = init_window()
    cam = camera(
        window,
        rot_v = rot_v,
        trans_v = trans_v,
        position = cam_position,
        lookat = lookat,
        up = up
    )

    hotkeys = Hotkeys(window)
    init_screenshots!(hotkeys)
    angles = init_angles(hotkeys)
    init_QWEASD_controls!(hotkeys, cam, angles)

    for c in cubics
        draw(
            window,
            expand(c),
            cam = cam,
            color = RGBA{Float32}(rand(3)..., 1.0)
        )
    end

    window, cam, hotkeys
end

function draw(
        window::GLWindow.Screen,
        pos::Vector{Point3f0};
        cam = :perspective,
        color::RGBA{Float32} = RGBA{Float32}(0.2, 0.4, 0.8, 1.0),
        primitive::AbstractMesh = glsphere()
    )

    _view(visualize(
        (primitive, pos),
        color = color
    ), window, camera = cam)
end


text:
import Serialization

function stack(io::IO,msg::Vector{UInt8})
    frontbytes = reinterpret(UInt8,Int16[length(msg)])
    item = UInt8[frontbytes...,msg...]
    write(io,item)
end

function unstack(io::IO)
    sizebytes = [read(io,UInt8),read(io,UInt8)]
    size = reinterpret(Int16,sizebytes)[1]

    msg = UInt8[]
    for i in 1:size
        push!(msg,read(io,UInt8))
    end
    return msg
end

function unstack(io::IOBuffer)
    bytes = take!(io)
    size = reinterpret(Int16,bytes[1:2])[1]
    msg = bytes[3:size+2]
    if length(bytes)>size+2
        write(io,bytes[size+3:end])
    end
    return msg
end

function serialize(socket,msg)
    io = IOBuffer()
    Serialization.serialize(io,msg)
    bytes = take!(io)
    stack(socket,bytes)
end

function deserialize(socket)
    bytes = unstack(socket)
    io = IOBuffer(bytes)
    Serialization.deserialize(io)
end