text:
using JSON, Printf
import Base: merge, length
import Mill: reflectinmodel

abstract type AbstractExtractor end
abstract type JSONEntry end
StringOrNumber = Union{String,Number}
max_keys = 10000

# so I can pretty print symbols
length(s::Symbol) = length(string(s))

function updatemaxkeys!(n::Int)
	global max_keys = n
end

"""
	mutable struct Entry <: JSONEntry
		counts::Dict{Any,Int}
		updated::Int
	end

	Keeps statistics about scalar values of a one key and also about items inside a key
	`count` counts how many times given value appeared (at most max_keys is held)
	`updated` counts how many times the entry was updated
"""
mutable struct Entry <: JSONEntry
	counts::Dict{Any,Int}
	updated::Int
end

Entry() = Entry(Dict{Any,Int}(),0);
types(e::Entry) = unique(typeof.(collect(keys(e.counts))))
Base.keys(e::Entry) = sort(collect(keys(e.counts)))
function Base.show(io::IO, e::Entry;pad =[], key = "")
	key *= isempty(key) ? ""  : ": "
	paddedprint(io, @sprintf("%s[Scalar - %s], %d unique values, updated = %d\n",key,join(types(e)),length(keys(e.counts)),e.updated))
end

function suggestextractor(e::Entry, settings = NamedTuple())
	t = promote_type(unique(typeof.(keys(e.counts)))...)
	t == Any  && @error "JSON does not have a fixed type scheme, quitting"

	for (c, ex) in get(settings, :scalar_extractors, default_scalar_extractor())
		c(e) && return ex(e)
	end
end

function default_scalar_extractor()
	[(e -> (length(keys(e.counts)) / e.updated < 0.1  && length(keys(e.counts)) <= 10000),
		e -> ExtractCategorical(collect(keys(e.counts)))),
	(e -> true,
		e -> extractscalar(promote_type(unique(typeof.(keys(e.counts)))...))),]
end

"""
		function update!(a::Entry, v)

		updates the entry when seeing value v
"""
function update!(a::Entry, v)
	if length(keys(a.counts)) < max_keys
		a.counts[v] = get(a.counts,v,0) + 1
		# it's there because otherwise, after filling the count keys not even the existing ones are updates
	elseif haskey(a.counts, v)
		a.counts[v] += 1
	end
	a.updated +=1
end

function merge(es::Entry...)
	updates_merged = sum(map(x->x.updated, es))
	counts_merged = merge(+, map(x->x.counts, es)...)
	if length(counts_merged) > max_keys
		counts_merged_list = sort(collect(counts_merged), by=x->x[2], rev=true)
		counts_merged = Dict(counts_merged_list[1:max_keys])
	end
	Entry(counts_merged, updates_merged)
end

"""
		mutable struct ArrayEntry <: JSONEntry
			items
			l::Dict{Int,Int}
			updated::Int
		end

		keeps statistics about an array entry in JSON.
		`items` is typeof `Entry` or nothing and keeps statistics about the elements of the array
		`l` keeps histogram of message length
		`updated` counts how many times the struct was updated.
"""
mutable struct ArrayEntry <: JSONEntry
	items	# Tried to make type stable optional type, didn't work
	l::Dict{Int,Int}
	updated::Int
end

ArrayEntry(items) = ArrayEntry(items,Dict{Int,Int}(),0)

function Base.show(io::IO, e::ArrayEntry; pad = [], key = "")
  c = COLORS[(length(pad)%length(COLORS))+1]
  # paddedprint(io,"Vector with $(length(e.items)) items(s). (updated = $(e.updated))\n", color=c)
  if isnothing(e.items)
	paddedprint(io, "$(key): [Empty List] (updated = $(e.updated))\n", color=c)
	return
  end
  paddedprint(io,"$(key): [List] (updated = $(e.updated))\n", color=c)
  paddedprint(io, "  └── ", color=c, pad=pad)
  show(io, e.items, pad = [pad; (c, "      ")])
end

function update!(a::ArrayEntry, b::Vector)
	n = length(b)
	a.updated +=1
	a.l[n] = get(a.l,n,0) + 1
	n == 0 && return
	if isnothing(a.items)
		 a.items = newentry(b).items
	end
	# foreach(v -> update!(a.items,v), b)
	for v in b
		update!(a.items,v)
	end
end

function suggestextractor(node::ArrayEntry, settings = NamedTuple())
	if isnothing(node.items)
		throw(ArgumentError("empty array, can not suggest extractor"))
	end
	e = suggestextractor(node.items, settings)
	isnothing(e) ? e : ExtractArray(e)
end

Base.isempty(e::ArrayEntry) = isnothing(e.items)
# todo: implement merding empty arrays
# todo: implement isempty for arrayentry so it's semantic
function merge(es::ArrayEntry...)
	updates_merged = sum(map(x->x.updated, es))
	l_merged = merge(+, map(x->x.l, es)...)
	nonempty_items = map(x->x.items, filter(!isempty, collect(es)))
	items_merged = isempty(nonempty_items) ? nothing : merge(merge, nonempty_items...)
	ArrayEntry(items_merged, l_merged, updates_merged)
end

"""
		mutable struct DictEntry <: JSONEntry
			childs::Dict{String,Any}
			updated::Int
		end

		keeps statistics about an object in json
		`childs` maintains key-value statistics of childrens. All values should be JSONEntries
		`updated` counts how many times the struct was updated.
"""

mutable struct DictEntry <: JSONEntry
	childs::Dict{Symbol, Any}
	updated::Int
end

DictEntry() = DictEntry(Dict{Symbol,Any}(),0)
Base.getindex(s::DictEntry, k::Symbol) = s.childs[k]
Base.getindex(s::DictEntry, k::String) = s.childs[Symbol(k)]
Base.setindex!(s::DictEntry, i, k::Symbol) = s.childs[k] = i
Base.setindex!(s::DictEntry, i, k::String) = s.childs[Symbol(k)] = i
Base.get(s::Dict{Symbol, <:Any}, key::String, default) = get(s, Symbol(key), default)

function json(io::IO, e::DictEntry)
	println(io,"{")
	println(io,"\"updated\": $(e.updated)")
	for k in keys(e.childs)
		print(io, "\"$(k)\": ")
		json(io, e.childs[k])
		println(io)
	end
	println(io,"}")
end

function Base.show(io::IO, e::DictEntry; pad=[], key = "")
    c = COLORS[(length(pad)%length(COLORS))+1]
    k = sort(collect(keys(e.childs)))
    if isempty(k)
    	paddedprint(io, "$(key)[Empty Dict] (updated = $(e.updated))\n", color=c)
    	return
    end
    ml = maximum(length.(k))
    key *= ": "
	  paddedprint(io, "$(key)[Dict] (updated = $(e.updated))\n", color=c)

    for i in 1:length(k)-1
    	s = "  ├──"*"─"^(ml-length(k[i]))*" "
			paddedprint(io, s, color=c, pad=pad)
			show(io, e.childs[k[i]], pad=[pad; (c, "  │"*" "^(ml-length(k[i])+2))], key = string(k[i]))
    end
    s = "  └──"*"─"^(ml-length(k[end]))*" "
    paddedprint(io, s, color=c, pad=pad)
    show(io, e.childs[k[end]], pad=[pad; (c, " "^(ml-length(k[end])+4))], key = string(k[end]))
end

function update!(s::DictEntry, d::Dict)
	s.updated +=1
	for (k,v) in d
		v == nothing && continue
		i = get(s.childs, k, newentry(v))
		update!(i,v)
		s[k] = i
	end
end

"""
		suggestextractor(e::DictEntry, settings = NamedTuple())

		create convertor of json to tree-structure of `DataNode`

		`e` top-level of json hierarchy, typically returned by invoking schema
		`settings.mincount` contains minimum repetition of the key to be included into
		the extractor (if missing it is equal to zero)
		`settings` can be any container supporting `get` function
"""
function suggestextractor(e::DictEntry, settings = NamedTuple())
	mincount = get(settings, :mincount, 0)
	ks = Iterators.filter(k -> updated(e.childs[k]) > mincount, keys(e.childs))
	isempty(ks) && return(nothing)
	c = [(string(k),suggestextractor(e.childs[k], settings)) for k in ks]
	c = filter(s -> s[2] != nothing, c)
	isempty(c) && return(nothing)
	mask = map(i -> extractsmatrix(i[2]), c)
	ExtractBranch(Dict(c[mask]),Dict(c[.! mask]))
end

function merge(es::DictEntry...)
	updates_merged = sum(map(x->x.updated, es))
	childs_merged = merge(merge, map(x->x.childs, es)...)
	DictEntry(childs_merged, updates_merged)
end

"""
		newentry(v)

		creates new entry describing json according to the type of v
"""
newentry(v::Dict) = DictEntry()
newentry(v::A) where {A<:StringOrNumber} = Entry()
newentry(v::Vector) = isempty(v) ? ArrayEntry(nothing) : ArrayEntry(newentry(v[1]))

"""
		function schema(a::Vector{T}) where {T<:Dict}
		function schema(a::Vector{T}) where {T<:AbstractString}

		create schema from an array of parsed or unparsed JSONs
"""
function schema(a::Vector{T}) where {T<:Dict}
	schema = DictEntry()
	foreach(f -> update!(schema,f),a)
	schema
end

function schema(a::Vector{T}) where {T<:AbstractString}
	schema = DictEntry()
	foreach(f -> update!(schema,JSON.parse(f)), a)
	schema
end

updated(s::T) where {T<:JSONEntry} = s.updated
merge(combine::typeof(merge), es::JSONEntry...) = merge(es...)

sample_synthetic(e::Entry) = first(keys(e.counts))
sample_synthetic(e::ArrayEntry) = repeat([sample_synthetic(e.items)], 2)
sample_synthetic(e::DictEntry) = Dict(k => sample_synthetic(v) for (k, v) in e.childs)

Base.hash(e::Entry, h::UInt) = hash((e.counts, e.updated), h)
Base.hash(e::ArrayEntry, h::UInt) = hash((e.items, e.l, e.updated), h)
Base.hash(e::DictEntry, h::UInt) = hash((e.childs, e.updated), h)

Base.:(==)(e1::Entry, e2::Entry) = e1.updated === e2.updated && e1.counts == e2.counts
Base.:(==)(e1::ArrayEntry, e2::ArrayEntry) = e1.updated === e2.updated && e1.l == e2.l && e1.items == e2.items
Base.:(==)(e1::DictEntry, e2::DictEntry) = e1.updated === e2.updated && e1.childs == e2.childs

reflectinmodel(sch::JSONEntry, ex::AbstractExtractor, db, da=d->SegmentedMean(d); b = Dict(), a = Dict()) =
	reflectinmodel(ex(sample_synthetic(sch)), db, da, b=b, a=a)

function make_selector(s::Symbol)
	if s == Symbol("[]")
		return (d) -> d.items
	else
		return (d) -> d.childs[s]
	end
end

"""
		Deletes `field` at the specified `path` from the schema `sch`.
		For instance, the following:
			`delete!(schema, ".field.subfield.[]", "x")``
		deletes the field `x` from `schema` at:
			`schema.childs[:field].childs[:subfield].items.childs`
"""
function Base.delete!(sch::JsonGrinder.JSONEntry, path::AbstractString, field::AbstractString)
	@assert field != "[]"

	selectors = map(Symbol, split(path, ".")[2:end])

	item = reduce((s, f) -> f(s), map(make_selector, selectors), init=sch)

	delete!(item.childs, Symbol(field))
end


text:
differ = -(@view(states[1:end-1]), @view(states[2:end]))
emp_var = empirical_var(differ,n, k)#sum(sum([d .^2 for d in differ]))/(2*k)
epsilon = 5*q/√(n*k)
@test abs(emp_var-q) < epsilon

ext: jl


text:
function perf_hdindexing{T}(u_i::AbstractArray{T, 6}, u_ip1::AbstractArray{T, 6})
  d1min = 3
  d1max = size(u_i, 5) - 2
  d2min = 3
  d2max = size(u_i, 4) - 2
  d3min = 3
  d3max = size(u_i, 3) - 2
  d4min = 3
  d4max = size(u_i, 2) - 2
  d5min = 3
  d5max = size(u_i, 1) - 2

  for d1 = d1min:d1max
    for d2 = d2min:d2max
      for d3 = d3min:d3max
        for d4 = d4min:d4max
          for d5 = d5min:d5max
            idx = (d5, d4, d3, d2, d1)
            kernel5(idx, u_i, u_ip1)
          end
        end
      end
    end
  end

  return nothing
end

function kernel5{T}(idx, u_i::AbstractArray{T,6}, u_ip1::AbstractArray{T,6})

  d1 = idx[1]
  d2 = idx[2]
  d3 = idx[3]
  d4 = idx[4]
  d5 = idx[5]

  delta_12 = 0.5
  delta_22 = 0.5
  delta_32 = 0.5
  delta_42 = 0.5
  delta_52 = 0.5

  u11_1 =  (1/12)*u_i[ d1 - 2, d2, d3, d4, d5, 1] +
           (-2/3)*u_i[ d1 - 1, d2, d3, d4, d5, 1] +
              (0)*u_i[ d1    , d2, d3, d4, d5, 1] +
            (2/3)*u_i[ d1 + 1, d2, d3, d4, d5, 1] +
          (-1/12)*u_i[ d1 + 2, d2, d3, d4, d5, 1]

  u22_1 =  (1/12)*u_i[ d1, d2 - 2, d3, d4, d5, 1] +
           (-2/3)*u_i[ d1, d2 - 1, d3, d4, d5, 1] +
              (0)*u_i[ d1, d2    , d3, d4, d5, 1] +
            (2/3)*u_i[ d1, d2 + 1, d3, d4, d5, 1] +
          (-1/12)*u_i[ d1, d2 + 2, d3, d4, d5, 1]

  u33_1 =  (1/12)*u_i[ d1, d2, d3 - 2, d4, d5, 1] +
           (-2/3)*u_i[ d1, d2, d3 - 1, d4, d5, 1] +
              (0)*u_i[ d1, d2, d3    , d4, d5, 1] +
            (2/3)*u_i[ d1, d2, d3 + 1, d4, d5, 1] +
          (-1/12)*u_i[ d1, d2, d3 + 2, d4, d5, 1]

  u44_1 =  (1/12)*u_i[ d1, d2, d3, d4 - 2, d5, 1] +
           (-2/3)*u_i[ d1, d2, d3, d4 - 1, d5, 1] +
              (0)*u_i[ d1, d2, d3, d4    , d5, 1] +
            (2/3)*u_i[ d1, d2, d3, d4 + 1, d5, 1] +
          (-1/12)*u_i[ d1, d2, d3, d4 + 2, d5, 1]

  u55_1 =  (1/12)*u_i[ d1, d2, d3, d4, d5 - 2, 1] +
           (-2/3)*u_i[ d1, d2, d3, d4, d5 - 1, 1] +
              (0)*u_i[ d1, d2, d3, d4, d5    , 1] +
            (2/3)*u_i[ d1, d2, d3, d4, d5 + 1, 1] +
          (-1/12)*u_i[ d1, d2, d3, d4, d5 + 2, 1]

  u_ip1[ d1, d2, d3, d4, d5, 2] = delta_12*u11_1 + delta_22*u22_1 + delta_32*u33_1 +
                                  delta_42*u44_1 + delta_52*u55_1

  u11_2 =  (1/12)*u_i[ d1 - 2, d2, d3, d4, d5, 2] +
           (-2/3)*u_i[ d1 - 1, d2, d3, d4, d5, 2] +
              (0)*u_i[ d1    , d2, d3, d4, d5, 2] +
            (2/3)*u_i[ d1 + 1, d2, d3, d4, d5, 2] +
          (-1/12)*u_i[ d1 + 2, d2, d3, d4, d5, 2]

  u22_2 =  (1/12)*u_i[ d1, d2 - 2, d3, d4, d5, 2] +
           (-2/3)*u_i[ d1, d2 - 1, d3, d4, d5, 2] +
              (0)*u_i[ d1, d2    , d3, d4, d5, 2] +
            (2/3)*u_i[ d1, d2 + 1, d3, d4, d5, 2] +
          (-1/12)*u_i[ d1, d2 + 2, d3, d4, d5, 2]

  u33_2 =  (1/12)*u_i[ d1, d2, d3 - 2, d4, d5, 2] +
           (-2/3)*u_i[ d1, d2, d3 - 1, d4, d5, 2] +
              (0)*u_i[ d1, d2, d3    , d4, d5, 2] +
            (2/3)*u_i[ d1, d2, d3 + 1, d4, d5, 2] +
          (-1/12)*u_i[ d1, d2, d3 + 2, d4, d5, 2]

  u44_2 =  (1/12)*u_i[ d1, d2, d3, d4 - 2, d5, 2] +
           (-2/3)*u_i[ d1, d2, d3, d4 - 1, d5, 2] +
              (0)*u_i[ d1, d2, d3, d4    , d5, 2] +
            (2/3)*u_i[ d1, d2, d3, d4 + 1, d5, 2] +
          (-1/12)*u_i[ d1, d2, d3, d4 + 2, d5, 2]

  u55_2 =  (1/12)*u_i[ d1, d2, d3, d4, d5 - 2, 2] +
           (-2/3)*u_i[ d1, d2, d3, d4, d5 - 1, 2] +
              (0)*u_i[ d1, d2, d3, d4, d5    , 2] +
            (2/3)*u_i[ d1, d2, d3, d4, d5 + 1, 2] +
          (-1/12)*u_i[ d1, d2, d3, d4, d5 + 2, 2]

  u_ip1[ d1, d2, d3, d4, d5, 1] = delta_12*u11_2 + delta_22*u22_2 + delta_32*u33_2 +
                                  delta_42*u44_2 + delta_52*u55_2

  return nothing
end


text:
# norm of off diagonal terms of a square matrix
function norm_off(M)
    if size(M[1],1)>1
        return sqrt(sum(abs2(M[i,j]) + abs2(M[j,i]) for i in 1:size(M,1) for j in i+1:size(M,1)))
    else
        return 0.0
    end
end

function diagonalization_iter(D)
    n = size(D[1],1)
    s = length(D)

    X = fill(zero(D[1][1,1]),n,n)
    Y = fill(zero(D[1][1,1]),n,n)

    A = fill(zero(D[1][1,1]),s,2)
    b = fill(zero(D[1][1,1]),s)
    for i in 1:n
        for j in 1:n
            if i != j
                for k in 1:s
                    A[k,1] = D[k][i,i]
                    A[k,2] = D[k][j,j]
                    b[k]   = -D[k][i,j]
                end
                v = A\b
                X[i,j] =  v[1]
                Y[i,j] =  v[2]
            end
        end
    end
    for i in 1:n
        X[i,i]=1
        Y[i,i]=1
    end
    return X, Y
end

function diagonalization(M::Vector{Matrix{C}},
                         Info = Dict{String,Any}(
                             "maxIter" => 10,
                             "epsIter" => 1.e-3)) where C
    n  = length(M)
    r  = size(M[1],1)

    N   = (haskey(Info,"maxIter") ? Info["maxIter"] : 10)
    eps = (haskey(Info,"epsIter") ? Info["epsIter"] : 1.e-3)

    M1 = sum(M[i]*randn(Float64) for i in 1:n)
    E  = eigvecs(M1)

    F  = inv(E)

    D  = vcat([Matrix{C}(I,r,r)],[F*M[i]*E for i in 1:length(M)])
    err = sum(norm_off.(D))
    delta = sum(norm.(D))
    #println("diag off: ", err)

    Info["d0"] = err
    nit = 0

    if err/delta > 5.e-2
        delta = err
        while nit < N && delta > eps
            err0 = err
            X,Y = diagonalization_iter(D)
            D = [Y*D[i]*X for i in 1:length(D)]
            E = E*X
            F = Y*F
            nit+=1
            err = sum(norm_off.(D))
            delta = err0-err
            #println("Off", nit,": ", err, "   delta: ", delta)
        end
        Info["d*"]= err
    end
    Info["nIter"] = nit

    Xi = fill(zero(E[1,1]),n,r)
    for i in 1:r
    	for j in 1:n
	    Xi[j,i] = D[j+1][i,i]/D[1][i,i]
            #Xi[j,i] =(E[:,i]\(M[j]*E[:,i]))[1]
	end
    end
    return Xi, E, Info
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

module UsageManagementClient

using Random
using Dates
using HTTP
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_RateCardApi.jl")
include("api_UsageAggregatesApi.jl")

# export models
export convert, ErrorResponse
export convert, InfoField
export convert, MeterInfo
export convert, OfferTermInfo
export convert, RateCardQueryParameters
export convert, ResourceRateCardInfo
export convert, UsageAggregation
export convert, UsageAggregationListResult
export convert, UsageSample
export convert, MonetaryCommitment
export convert, MonetaryCredit
export convert, RecurringCharge

# export operations
export convert, RateCardApi, UsageAggregatesApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
using LinearAlgebra: qr!
using Dates: @dateformat_str, now, format
using Printf: @sprintf

push!(LOAD_PATH, @__DIR__)
using MultiprecisionFloats
using RKTK2
using DZMisc: say, approx_norm, orthonormalize_columns!

const AccurateReal = Float64x2
setprecision(256)

################################################################################

const EVALUATOR = RKOCEvaluator{AccurateReal}(10, 16)
const ERROR_EVALUATOR = RKOCEvaluator{AccurateReal}(11, 16)
const NUM_VARS = 136
const NUM_CONSTRS = 1205

const CONSTR_IDXS = [
    1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 18, 19, 20, 21, 22, 23,
    24, 25, 26, 27, 28, 31, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 51, 52,
    58, 59, 62, 75, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99,
    100, 101, 106, 107, 110, 115, 116, 123, 134, 135, 201, 202, 203, 205, 206,
    207, 208, 210, 211, 214, 215, 221, 222, 225, 226, 227, 230, 231, 238, 249,
    250, 269, 270, 316, 404, 406, 487, 488, 489, 490, 491, 492, 493, 496, 497,
    498, 499, 500, 501, 502, 507, 508, 511, 516, 517, 524, 535, 536, 539, 544,
    555, 556, 602, 603, 650, 690, 984]

const CONSTR_DIM = length(CONSTR_IDXS)

const TEMP_RES = Vector{AccurateReal}(undef, NUM_CONSTRS)
const SHORT_RES = Vector{AccurateReal}(undef, CONSTR_DIM)
const TEMP_JAC = Matrix{AccurateReal}(undef, NUM_CONSTRS, NUM_VARS)
const SHORT_JAC = Matrix{AccurateReal}(undef, NUM_VARS, CONSTR_DIM)
const SHORT_JCT = Matrix{AccurateReal}(undef, CONSTR_DIM, NUM_VARS)
const APX_SHORT_JAC = Matrix{Float64}(undef, NUM_VARS, CONSTR_DIM)
const ERROR_COEFFS = Vector{AccurateReal}(undef, ERROR_EVALUATOR.num_constrs)
const ERROR_JAC = Matrix{AccurateReal}(undef, ERROR_EVALUATOR.num_constrs, NUM_VARS)

################################################################################

function approx_objective(x)
    evaluate_residual!(TEMP_RES, x, EVALUATOR)
    result = zero(Float64)
    for i = 1 : NUM_CONSTRS
        result += abs2(Float64(TEMP_RES[i]))
    end
    result
end

function compute_residual(x)
    evaluate_residual!(TEMP_RES, x, EVALUATOR)
    @simd ivdep for j = 1 : length(CONSTR_IDXS)
        @inbounds SHORT_RES[j] = TEMP_RES[CONSTR_IDXS[j]]
    end
end

function compute_jacobian(x)
    evaluate_jacobian!(TEMP_JAC, x, EVALUATOR)
    for j = 1 : length(CONSTR_IDXS)
        k = CONSTR_IDXS[j]
        @simd ivdep for i = 1 : NUM_VARS
            @inbounds SHORT_JCT[j, i] = TEMP_JAC[k, i]
        end
    end
end

function compute_orthonormalized_jacobian(x)
    evaluate_jacobian!(TEMP_JAC, x, EVALUATOR)
    for j = 1 : length(CONSTR_IDXS)
        k = CONSTR_IDXS[j]
        @simd ivdep for i = 1 : NUM_VARS
            @inbounds SHORT_JAC[i, j] = TEMP_JAC[k, i]
        end
    end
    orthonormalize_columns!(SHORT_JAC)
    for j = 1 : length(CONSTR_IDXS)
        @simd ivdep for i = 1 : NUM_VARS
            @inbounds APX_SHORT_JAC[i, j] = Float64(SHORT_JAC[i, j])
        end
    end
end

function compute_error_coefficients(point)
    evaluate_error_coefficients!(ERROR_COEFFS, point, ERROR_EVALUATOR)
end

function compute_error_jacobian(point)
    evaluate_error_jacobian!(ERROR_JAC, point, ERROR_EVALUATOR)
end

function approx_force(point::Vector{T}) where {T <: Real}
    compute_error_coefficients(point)
    compute_error_jacobian(point)
    force = transpose(ERROR_JAC[1206:end, :]) * ERROR_COEFFS[1206:end]
    compute_orthonormalized_jacobian(point)
    force - APX_SHORT_JAC * (transpose(APX_SHORT_JAC) * force)
end

################################################################################

function constrain(x)
    x_old, obj_old = x, approx_objective(x)
    while true
        compute_residual(x_old)
        compute_jacobian(x_old)
        direction = qr!(SHORT_JCT) \ SHORT_RES
        x_new = x_old - direction
        obj_new = approx_objective(x_new)
        if obj_new < obj_old
            x_old, obj_old = x_new, obj_new
        else
            break
        end
    end
    x_old, obj_old
end

const EPS_THREE_HALVES = Float64(BigFloat(2)^-130)

function perturb(x, direction, multiplier)
    x_new, obj = constrain(x + multiplier * direction)
    if obj < EPS_THREE_HALVES
        x_new, multiplier
    else
        multiplier /= 2
        perturb(x, direction, multiplier)
    end
end

################################################################################

const INPUT_FILENAME = maximum(filename
    for filename in readdir()
    if isfile(filename) && startswith(filename, "RKTK-ERROPT-")
                        && endswith(filename, ".txt"))

say("Reading initial point from data file: ", INPUT_FILENAME)
const INPUT_POINT = AccurateReal.(BigFloat.(split(read(INPUT_FILENAME, String))))

@assert length(INPUT_POINT) == NUM_VARS
say("Successfully read initial point.")

function main()
    x = INPUT_POINT[:]
    compute_error_coefficients(x)
    old_norm = approx_norm(ERROR_COEFFS[1206:end])
    println(old_norm)
    speed = 2500.0
    while true
        for _ = 1 : 1000
            force = approx_force(x)
            while true
                force *= inv(-speed * approx_norm(force))
                x_new, multiplier = perturb(x, force, 1.0)
                compute_error_coefficients(x_new)
                new_norm = approx_norm(ERROR_COEFFS[1206:end])
                if new_norm < old_norm
                    x = x_new
                    old_norm = new_norm
                    if multiplier == 1.0
                        speed *= 0.95
                    else
                        speed *= 2.0
                    end
                    println(new_norm, " ", speed, " ", multiplier)
                    break
                else
                    speed *= 2.0
                end
            end
        end
        filename = format(now(), dateformat"RKTK-\ERROPT-yyyymmdd-HHMMSS-sss.txt")
        say("Writing points to file: ", filename)
        write(filename, join(string.(BigFloat.(x)), "\n") * "\n")
    end
end

main()


text:
include(joinpath("src", "HNScraper.jl"))
using HTTP
using JSON

# see an example story
story_endpoint = "https://hacker-news.firebaseio.com/v0/item/25623858.json"
response = HTTP.get(story_endpoint)
response_body = String(response.body)
story = JSON.parse(response_body)
for key in keys(story)
    print(key, "\t", typeof(story[key]), "\n")
end

story["descendants"]

# H1: ranking has a direct influence on voting
# H2: influence of voting on ranking is not linear
# H3: distribution of votes over time is sinoidal

test = HNScraper.scrape_topstories(top_n = 10, n_timeslices = 2, interval = 2)

tss = unique(test["timestamp"])

Dates.unix2datetime(tss[1])

# scrape and write for preliminary analysis
prelim_data = HNScraper.scrape_topstories(top_n = 500, n_timeslices = 60, interval = 5)
if !("data" in readdir())
    mkdir("data")
end
Feather.write(joinpath("data", "five_hours.feather"), prelim_data)

begin
    i = 0
    cb(timer) = (global i += 1; println(i))
    t = Timer(cb, 4, interval=0.2)
    wait(t)
    sleep(0.5)
    close(t)
end

for i in 1:5
    t = Timer(5)
    @info "Now sleeping for 3 s."
    sleep(3)
    @info "Now waiting..."
    wait(t)
    print("Finally.\n")
end

HNScraper.scrape_topstories(n_timeslices = 2, interval_minutes = 2)

max_stars_repo_name: JohannesNakayama/hacker-news-analysis

max_issues_repo_name: JohannesNakayama/hacker-news-analysis

max_forks_repo_name: JohannesNakayama/hacker-news-analysis


text:
"""
    clone(env)

Create a clone of CommonEnv `env` at the current state.

Two clones are assumed to be completely independent of each other - no action applied to one will affect the other.
"""
function clone end

"""
    render(env)

Return a showable object that visualizes the environment at the current state.

Implementing this function will facilitate visualization through the [Julia Multimedia I/O system](https://docs.julialang.org/en/v1/base/io-network/#Multimedia-I/O-1), for example, calling `display(render(env))` will cause the visualization to pop up in a Jupyter notebook, the Juno plot pane, or an ElectronDisplay.jl window. `render` should return an object that has `Base.show` methods for many MIME types. Good examples include a Plots.jl plot, a Compose.jl `Context` graphic produced by Cairo, or a custom type with `show` methods.

# Example
```julia
using CommonRLInterface
using Plots

struct MyEnv <: AbstractEnv
    state::Tuple{Int, Int}
end

@provide CommonRLInterface.render(env::MyEnv) = scatter(env.state[1], env.state[2])
```
"""
function render end

"""
    state(env::AbstractEnv)

Return the state of the environment.

See setstate! for more information about the state.
"""
function state end

"""
    setstate!(env::AbstractEnv, s)

Set the state of the environment to s.

If two environments have the same state, future outputs (observations and rewards) should be statistically identical given the same sequence of actions.
"""
function setstate! end


text:
using HarwellRutherfordBoeing
using Krylov
using LinearOperators
# using ProfileView

# M = HarwellBoeingMatrix("data/illc1033.rra");
M = HarwellBoeingMatrix("data/illc1850.rra");
A = M.matrix;
(m, n) = size(A);
@printf("System size: %d rows and %d columns\n", n, m);  # We use A'.

# Define a linear operator with preallocation.
Ap = zeros(m);
Atq = zeros(n);
op = LinearOperator(m, n, false, false,
                    p -> A_mul_B!(1.0,  A, p, 0.0, Ap),
                    q -> At_mul_B!(1.0, A, q, 0.0, Atq),
                    q -> At_mul_B!(1.0, A, q, 0.0, Atq));

x_exact = A * ones(n);  # y = ones(n);
x_exact_norm = norm(x_exact);
x_exact /= x_exact_norm;
b = A' * x_exact;
(x, stats) = cgne(op', b);
# @profile (x, stats) = cgne(op', b);
@time (x, stats) = cgne(op', b);
show(stats);
resid = norm(A' * x - b) / norm(b);
@printf("CGNE: Relative residual: %7.1e\n", resid);
@printf("CGNE: ‖x - x*‖₂: %7.1e\n", norm(x - x_exact));

# ProfileView.view()


text:
using FactCheck

include("blapi.jl")

exitstatus()


text:
using Test
using NeXLCore
using NeXLSpectrum

@testset "Al2O3" begin
    al2o3 = loadspectrum(joinpath(@__DIR__, "K412 spectra", "Al2O3 std.msa"))
    @test isapprox(al2o3[:LiveTime], 1172.19288, atol = 0.00001)
    @test isapprox(al2o3[:BeamEnergy], 20.0e3, atol = 1.0)
    @test isapprox(al2o3[:ProbeCurrent], 1.10989, atol = 0.00001)
    @test isapprox(al2o3[:RealTime], 1491.4828, atol = 0.0001)
    @test isapprox(al2o3[:Elevation], deg2rad(35.0), atol = deg2rad(0.01))
    @test isapprox(material(al2o3[:Coating])[n"C"], 1.0, atol = 0.0001)
    @test isapprox(NeXLCore.thickness(al2o3[:Coating]), 10.0 * 1.e-7, atol = 1.0e-9)
    @test isapprox(al2o3[101], 22290.0, atol = 0.001)
    @test isapprox(energy(1, al2o3.energy), 1.63032, atol = 0.1)
    @test isapprox(energy(101, al2o3.energy), 100.0 * 9.99856 + 1.63032, atol = 0.1)
    @test isapprox(energy(1001, al2o3.energy), 1000.0 * 9.99856 + 1.63032, atol = 0.1)
    @test isapprox(energy(10001, al2o3.energy), 10000.0 * 9.99856 + 1.63032, atol = 0.1)
    @test channel(10000.0 * 9.99856 + 1.63032, al2o3.energy) == 10001
    @test channel(1.63032 - 1.0, al2o3.energy) == 0
    @test channel(100.0 * 9.99856 + 1.63032 + 1.0, al2o3.energy) == 101
    @test channel(100.0 * 9.99856 + 1.63032 - 1.0, al2o3.energy) == 100
    @test al2o3[:Name] == "Al2O3 std"
    @test isapprox(al2o3[:Composition][n"Al"], 0.529251, atol = 0.00001)
    @test isapprox(al2o3[:Composition][n"O"], 0.470749, atol = 0.00001)
    @test isapprox(al2o3[:WorkingDistance], 1.7, atol = 0.0001)
    @test al2o3[:Owner] == "Unknown"
end


text:
is_ch_type(::Val{:DateTime64}) = true

result_type(
    ::Val{:DateTime64},
    precision_str,
    timezone::Union{String, Nothing} = nothing,
) = Vector{DateTime}

function read_col_data(
    sock::ClickHouseSock,
    num_rows::VarUInt,
    ::Val{:DateTime64},
    precision_str::String,
    timezone::Union{String, Nothing} = nothing,
)
    precision = parse(Int, precision_str)
    data = chread(sock, Vector{Int64}, num_rows)

    # Maximum precision of DateTime is milliseconds
    datetime_muliplier =10^3 / 10^precision

    return DateTime(1970) .+ Millisecond.(round.(Int64, data .* datetime_muliplier))
end

function write_col_data(
    sock::ClickHouseSock,
    data::AbstractVector{DateTime},
    ::Val{:DateTime64},
    precision_str::String,
    timezone::Union{String, Nothing} = nothing,
)
    precision = parse(Int, precision_str)

    # Maximum precision of DateTime is milliseconds
    datetime_muliplier =10^3 / 10^precision

    d = round.(Int64,
        getproperty.(
            Millisecond.(data .- DateTime(1970)), :value
        ) ./ datetime_muliplier
    )
    chwrite(sock, d)
end


text:
using Distributions
using StatsPlots

plotpost(x; kwargs...) = _plotpost(x, false; kwargs...)
plotpost!(x; kwargs...) = _plotpost(x, true; kwargs...)

function _plotpost(x, add; cialpha=1, a=0.05, ms=8, truth=nothing, kwargs...)
  qs = [quantile(x, q) for q in [a/2, 1 - a/2]]
  if add
    histogram!(x, linealpha=0, grid=false, normalize=true, label=nothing; kwargs...)
  else
    histogram(x, linealpha=0, grid=false, normalize=true, label=nothing; kwargs...)
  end
  if cialpha > 0
    scatter!(qs, zero.(qs), ms=ms, alpha=cialpha, msc=kwargs[:color],
             color=kwargs[:color], label=nothing)
    # annotate!(qs[1], zero(qs[1]), text("<", 30, kwargs[:color]))
    # annotate!(qs[2], zero(qs[2]), text(">", 30, kwargs[:color]))
  end
  truth == nothing || vline!([truth], color=kwargs[:color], ls=:dot, lw=3)
  vline!([mean(x)], color=kwargs[:color])
end

# x = randn(10000)
# plotpost(x, color=:red, alpha=0.3, truth=0)
# plotpost!(x .+ 2, color=:blue, alpha=0.3, truth=2)

# scatter(randn(100), randn(100), msc=:blue, ms=10, color=:blue, alpha=0.5)

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
export maximum_radius_with_center, hchebyshevcenter, vchebyshevcenter, chebyshevcenter
using JuMP

"""
    maximum_radius_with_center(h::HRep, center)

Return the maximum radius `r` such that the Euclidean ball of center `center`
and radius `r` is included in the polyhedron `h`.
"""
function maximum_radius_with_center(h::HRep{T}, center) where T
    if any(hp -> !isapproxzero(coord(hp)), hyperplanes(h))
        return zero(T)
    else
        radius = nothing
        for hs in halfspaces(h)
            n = norm(hs.a, 2)
            if iszero(n)
                if hs.β < 0
                    # The polyhedron is empty
                    return zero(T)
                end
            else
                new_radius = (hs.β - center ⋅ hs.a) / n
                if radius === nothing
                    radius = new_radius
                else
                    radius = min(radius, new_radius)
                end
            end
        end
        if radius === nothing
            error("The polyhedron is the full space, its maximum radius is infinite.")
        elseif radius < 0
            error("The polyhedron does not contain the provided center $center.")
        end
        return radius
    end
end

_shrink(h::HyperPlane, radius, T::Type) = convert(similar_type(typeof(h), T), h)
_shrink(h::HalfSpace, radius, T::Type) = HalfSpace{T}(h.a, h.β - norm(h.a, 2) * radius)
function _shrink(p::HRep{Tin}, radius) where Tin
    T = MA.promote_operation(MA.add_mul, Tin, float(Tin), eltype(radius))
    f = (i, h) -> _shrink(h, radius, T)
    d = FullDim(p)
    return similar(p, d, T, hmap(f, d, T, p)...)
end

# TODO reference thesis where proper chebyshev center is defined
"""
    hchebyshevcenter(p::HRep[, solver]; linearity_detected=false, proper=true)

Return a tuple with the center and radius of the largest euclidean ball contained in the polyhedron `p`.
Throws an error if the polyhedron is empty or if the radius is infinite.
Linearity is detected first except if `linearity_detected`.

Note that a polytope may have several Chebyshev center.
In general, the set of Chebyshev center of a polytope `p` is a polytope which has a lower dimension than `p` if `p` has a positive dimension.
For instance, if `p` is the rectangle `[-2, 2] x [-1, 1]`, the Chebyshev radius of `p` is 1
and the set of Chebyshev centers is `[-1, 1] x {0}`.
The *proper* Chebyshev center is `(0, 0)`, the Chebyshev center of `[-1, 1] x {0}`.
If `!proper` then any Chebyshev center is returned (the one returned depends on the solver).
Otherwise the proper Chebyshev center is computed.
The proper Chebyshev center is defined by induction on the dimension of `p`.
If `p` has dimension 0 then it is a singleton and its proper Chebyshev center
  is the only element of `p`.
Otherwise, the dimension of the set `q` of Chebyshev centers of `p` is smaller than
the dimension of `p` and the proper Chebyshev center of `p` is the proper Chebyshev center of `q`.
"""
function hchebyshevcenter(p::HRepresentation, solver=default_solver(p; T=Float64); # Need Float64 for `norm(a, 2)`
                          linearity_detected=false, # workaround for https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/25
                          proper=true, verbose=1)
    if !linearity_detected
        p = detecthlinearity(p, solver)
    end
    model, T = layered_optimizer(solver)
    c = MOI.add_variables(model, fulldim(p))
    _constrain_in(model, hyperplanes(p), c, T)
    rs, _ = MOI.add_constrained_variables(model, MOI.Nonnegatives(1))
    r = rs[1]
    for (i, hs) in enumerate(halfspaces(p))
        func, set = _constrain_in_func_set(hs, c, T)
        push!(func.terms, MOI.ScalarAffineTerm{T}(norm(hs.a, 2), r))
        MOI.add_constraint(model, func, set)
    end
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(), MOI.SingleVariable(r))
    MOI.optimize!(model)
    term = MOI.get(model, MOI.TerminationStatus())
    if term ∉ [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
        if term == MOI.INFEASIBLE
            error("An empty polyhedron has no H-Chebyshev center.")
        elseif term == MOI.DUAL_INFEASIBLE
            error("The polyhedron contains euclidean ball of arbitrary large radius.")
        else
            _unknown_status(model, term, "computing the H-Chebyshev center.")
        end
    end
    radius = MOI.get(model, MOI.VariablePrimal(), r)
    if proper
        q = detecthlinearity(_shrink(p, radius), solver)
        p_dim = dim(p, true)
        q_dim = dim(q, true)
        if !iszero(q_dim)
            if q_dim >= p_dim
                error("The dimension of the set of Chebyshev centers `$q` is `$q_dim` while we expect it to have a smaller dimension than the original polyhedron which has dimension `$p_dim`.")
            end
            if verbose >= 1
                println("The set `q` of Chebyshev centers has dimension `$q_dim` which is nonzero but lower than the previous dimension `$p_dim`: we now compute the set of Chebyshev center of `q`.")
            end
            center, _ = hchebyshevcenter(q, solver, linearity_detected=true)
            return center, radius
        end
    end
    center = MOI.get(model, MOI.VariablePrimal(), c)
    return center, radius
end
function hchebyshevcenter(p::Polyhedron, solver=default_solver(p; T=Float64); linearity_detected=false, kws...)
    if !linearity_detected
        # We are going to detect it so we might as well do it at the level
        # of the `Polyhedron` so that it is saved and not recomputed if it is needed later.
        detecthlinearity!(p, solver) # FIXME `solver` was forced to be for `T=Float64` but it might not be the best choice here
    end
    return hchebyshevcenter(hrep(p), solver; linearity_detected=true, kws...)
end

# TODO solver here should not be VRepOptimizer
"""
    vchebyshevcenter(p::VRep[, solver])

Return a tuple with the center and radius of the smallest euclidean ball containing the polyhedron `p`.
Throws an error if the polyhedron is empty or if the radius is infinite (i.e. `p` is not a polytope, it contains rays).
"""
function vchebyshevcenter(p::VRep, solver=Polyhedra.default_solver(p))
    error("TODO")
end

"""
    chebyshevcenter(p::Rep[, solver])

If `p` is a H-representation or is a polyhedron for which the H-representation has already been computed, calls `hchebyshevcenter`, otherwise, call `vchebyshevcenter`.
"""
function chebyshevcenter(p::Polyhedron, solver=Polyhedra.default_solver(p; T=Float64); kws...)
    if hrepiscomputed(p)
        hchebyshevcenter(p, solver; kws...)
    else
        vchebyshevcenter(p, solver; kws...)
    end
end
chebyshevcenter(p::HRepresentation, solver=Polyhedra.default_solver(p; T=Float64)) = hchebyshevcenter(p, solver)
chebyshevcenter(p::VRepresentation, solver=Polyhedra.default_solver(p; T=Float64)) = vchebyshevcenter(p, solver)


text:
module TestWiringDiagramCore
using Test

using Catlab.WiringDiagrams.WiringDiagramCore
import Catlab.WiringDiagrams.WiringDiagramCore: validate_ports

# For testing purposes, check equality of port symbols.
function validate_ports(source_port::Symbol, target_port::Symbol)
  if source_port != target_port
    error("Source port $source_port does not equal target port $target_port")
  end
end

A, B, C, D = [:A], [:B], [:C], [:D]
f = Box(:f, A, B)
g = Box(:g, B, C)
h = Box(:h, C, D)

# Imperative interface
######################

# Operations on boxes
d = WiringDiagram(A, C)
@test nboxes(d) == 0
@test_throws KeyError box(d,input_id(d))
@test_throws KeyError box(d,output_id(d))

fv = add_box!(d, f)
@test nboxes(d) == 1
@test box(d, fv).value == :f
@test box(d, fv) == f
rem_box!(d, fv)
@test nboxes(d) == 0

fv = add_box!(d, f)
gv = add_box!(d, g)
hv = add_box!(d, h)
@test nboxes(d) == 3
@test [b.value for b in boxes(d)] == [:f,:g,:h]
@test boxes(d) == [f,g,h]
rem_boxes!(d, [fv,hv])
@test nboxes(d) == 1
@test boxes(d) == [g]

# Operations on ports
d = WiringDiagram(A, C)
fv = add_box!(d, f)
gv = add_box!(d, g)

@test input_ports(d, fv) == A
@test output_ports(d, fv) == B
@test input_ports(d, output_id(d)) == C
@test output_ports(d, input_id(d)) == A
@test_throws ErrorException input_ports(d, input_id(d))
@test_throws ErrorException output_ports(d, output_id(d))

@test port_value(d, Port(input_id(d),OutputPort,1)) == :A
@test port_value(d, Port(output_id(d),InputPort,1)) == :C
@test port_value(d, Port(fv,InputPort,1)) == :A
@test port_value(d, Port(fv,OutputPort,1)) == :B

# Operations on wires
@test nwires(d) == 0
@test !has_wire(d, fv, gv)
@test !has_wire(d, (fv,1) => (gv,1))
add_wire!(d, (input_id(d),1) => (fv,1))
add_wire!(d, (fv,1) => (gv,1))
add_wire!(d, (gv,1) => (output_id(d),1))
@test nwires(d) == 3
@test has_wire(d, fv, gv)
@test has_wire(d, (fv,1) => (gv,1))
@test_throws ErrorException add_wire!(d, (gv,1) => (fv,1))
@test wires(d) == map(Wire, [
  (input_id(d),1) => (fv,1),
  (fv,1) => (gv,1),
  (gv,1) => (output_id(d),1),
])

# Shallow copies.
d_copy = copy(d)
rem_boxes!(d_copy, [fv,gv])
@test nboxes(d) == 2
@test nwires(d) == 3

# Graph properties.
@test Set(all_neighbors(d, fv)) == Set([input_id(d),gv])
@test Set(all_neighbors(d, gv)) == Set([fv,output_id(d)])
@test neighbors(d, fv) == [gv]
@test outneighbors(d, fv) == [gv]
@test inneighbors(d, gv) == [fv]
@test wires(d, input_id(d)) == [ Wire((input_id(d),1) => (fv,1)) ]
@test wires(d, fv) == map(Wire, [
  ((input_id(d),1) => (fv,1)),
  ((fv,1) => (gv,1))
])
@test out_wires(d, fv) == [ Wire((fv,1) => (gv,1)) ]
@test out_wires(d, Port(fv,OutputPort,1)) == [ Wire((fv,1) => (gv,1)) ]
@test in_wires(d, gv) == [ Wire((fv,1) => (gv,1)) ]
@test in_wires(d, Port(gv,InputPort,1)) == [ Wire((fv,1) => (gv,1)) ]

rem_wires!(d, fv, gv)
@test nwires(d) == 2
@test !has_wire(d, fv, gv)
rem_wire!(d, (input_id(d),1) => (fv,1))
@test wires(d) == [ Wire((gv,1) => (output_id(d),1)) ]

# Induced subgraph.
d = WiringDiagram(A,D)
fv, gv, hv = add_box!(d, f), add_box!(d, g), add_box!(d, h)
add_wires!(d, Pair[
  (input_id(d),1) => (fv,1),
  (fv,1) => (gv,1),
  (gv,1) => (hv,1),
  (hv,1) => (output_id(d),1),
])
sub = WiringDiagram(A, D)
fv, gv = add_box!(sub, f), add_box!(sub, g)
add_wires!(sub, Pair[
  (input_id(sub),1) => (fv,1),
  (fv,1) => (gv,1),
])
@test induced_subdiagram(d, [fv, gv]) == sub

# Substitution
##############

sub = WiringDiagram(B,D)
gv = add_box!(sub, g)
hv = add_box!(sub, h)
add_wires!(sub, Pair[
  (input_id(sub),1) => (gv,1),
  (gv,1) => (hv,1),
  (hv,1) => (output_id(sub),1),
])
d0 = WiringDiagram(A,D)
fv = add_box!(d0, f)
subv = add_box!(d0, sub)
add_wires!(d0, Pair[
  (input_id(d0),1) => (fv,1),
  (fv,1) => (subv,1),
  (subv,1) => (output_id(d0),1),
])
@test boxes(d0) == [ f, sub ]
@test boxes(sub) == [ g, h ]
d = substitute(d0, subv)
@test nboxes(d) == 3
@test boxes(d) == [f, g, h]
box_map = Dict(box(d,v).value => v for v in box_ids(d))
@test nwires(d) == 4
@test Set(wires(d)) == Set(map(Wire, [
  (input_id(d),1) => (box_map[:f],1),
  (box_map[:f],1) => (box_map[:g],1),
  (box_map[:g],1) => (box_map[:h],1),
  (box_map[:h],1) => (output_id(d),1),
]))

# Encapsulation
###############

d0 = WiringDiagram(A,D)
fv = add_box!(d0, f)
gv = add_box!(d0, g)
hv = add_box!(d0, h)
add_wires!(d0, Pair[
  (input_id(d),1) => (fv,1),
  (fv,1) => (gv,1),
  (gv,1) => (hv,1),
  (hv,1) => (output_id(d),1)
])

d = encapsulate(d0, [fv,gv])
@test nboxes(d) == 2
@test sum(isa(b, WiringDiagram) for b in boxes(d)) == 1
@test nwires(d) == 3
sub = first(b for b in boxes(d) if isa(b, WiringDiagram))
@test nboxes(sub) == 2
@test boxes(sub) == [f, g]
box_map = Dict(box(sub,v).value => v for v in box_ids(sub))
@test Set(wires(sub)) == Set(map(Wire, [
  (input_id(sub),1) => (box_map[:f],1),
  (box_map[:f],1) => (box_map[:g],1),
  (box_map[:g],1) => (output_id(sub),1),
]))

d = encapsulate(d0, [fv,gv], discard_boxes=true, value=:e)
@test boxes(d) == [ Box(:e, A, C), h ]
box_map = Dict(box(d,v).value => v for v in box_ids(d))
@test Set(wires(d)) == Set(map(Wire, [
  (input_id(d),1) => (box_map[:e],1),
  (box_map[:e],1) => (box_map[:h],1),
  (box_map[:h],1) => (output_id(d),1),
]))

d0 = WiringDiagram(A,B)
v1 = add_box!(d0, f)
v2 = add_box!(d0, f)
add_wires!(d0, Pair[
  (input_id(d),1) => (v1,1),
  (input_id(d),1) => (v2,1),
  (v1,1) => (output_id(d),1),
  (v2,1) => (output_id(d),1),
])
d = encapsulate(d0, [v1,v2])
@test nboxes(d) == 1
@test nwires(d) == 2
sub = first(boxes(d))
@test sub == d0

end


text:
module QuadratureRules

    import FastGaussQuadrature
    import GenericLinearAlgebra
    import Polynomials
    import Polynomials: Polynomial

    include("utils.jl")

    include("quadrature_rule.jl")

    export QuadratureRule
    export nnodes, nodes, order, weights

    include("tabulated_quadratures.jl")

    export RiemannQuadratureLeft,
           RiemannQuadratureRight,
           MidpointQuadrature,
           TrapezoidalQuadrature

    include("chebyshev.jl")
    include("clenshaw_curtis.jl")
    include("gauss_legendre.jl")
    include("lobatto_legendre.jl")

    export ChebyshevQuadrature,
           GaussChebyshevQuadrature,
           LobattoChebyshevQuadrature,
           ClenshawCurtisQuadrature,
           GaussLegendreQuadrature,
           LobattoLegendreQuadrature

end


text:
module TestDiagrams
using Test

using Catlab.Theories, Catlab.Graphs, Catlab.CategoricalAlgebra
using Catlab.Graphs.BasicGraphs: TheoryGraph, TheorySymmetricGraph

const SchSGraph = TheorySymmetricGraph

# Diagrams
##########

# Diagram for paths of length 2.
C = FinCat(@acset Graph begin
  V = 3
  E = 2
  src = [1,2]
  tgt = [3,3]
end)
D = FinDomFunctor([:E,:E,:V], [:tgt,:src], C, FinCat(SchSGraph))
d = Diagram(D)
@test shape(d) == C
@test ob_map(d, 3) == SchSGraph[:V]
@test hom_map(d, 1) == SchSGraph[:tgt]
@test startswith(sprint(show, d), "Diagram{id}(")

# Diagram morphisms
###################

f = DiagramHom([(2,:inv), (1,:inv), 3], [2,1], d, d)
@test dom(f) == d
@test codom(f) == d
@test is_functorial(shape_map(f))
@test shape_map(f) == FinFunctor([2,1,3], [2,1], C, C)
ϕ = diagram_map(f)
@test is_natural(ϕ, check_equations=false)
@test ϕ[1] == SchSGraph[:inv]
@test ϕ[3] == id(SchSGraph[:V])
@test ob_map(f, 2) == (1, SchSGraph[:inv])
@test hom_map(f, 2) == Path(graph(C), 1)
@test collect_ob(f) == [(2, ϕ[1]), (1, ϕ[2]), (3, ϕ[3])]
@test collect_hom(f) == [Path(graph(C), 2), Path(graph(C), 1)]
f² = f⋅f
@test shape_map(f²) == FinFunctor(1:3, 1:2, C, C)

f = DiagramHom{op}([(2,:inv), (1,:inv), 3], [2,1], D, D)
ιV = FinDomFunctor([:V], FinCat(1), FinCat(SchSGraph))
g = DiagramHom{op}([(1,:src)], D, ιV)
@test dom(g) == Diagram{op}(D)
@test codom(g) == Diagram{op}(ιV)
@test compose(id(dom(g)), g) == g
@test compose(g, id(codom(g))) == g
@test is_natural(diagram_map(g))
@test startswith(sprint(show, f), "DiagramHom{op}(")

fg = f⋅g
@test ob_map(fg, 1) == (2, SchSGraph[:inv]⋅SchSGraph[:src])

d = dom(f)
@test op(op(d)) == d
@test op(op(f)) == f
@test dom(op(g)) == Diagram{co}(ιV)
@test codom(op(g)) == Diagram{co}(D)
@test op(g) == DiagramHom{co}([(1,:src)], ιV, D)
@test op(g)⋅op(f) == op(f⋅g)

# Monads of diagrams
####################

C = FinCat(TheoryGraph)
d = munit(Diagram{id}, C, :V)
@test is_discrete(shape(d))
@test only(collect_ob(diagram(d))) == TheoryGraph[:V]
f = munit(DiagramHom{id}, C, :src)
@test only(components(diagram_map(f))) == TheoryGraph[:src]

end


text:
function (tform::AffineMap)(a::SDFileAtom)
    return SDFileAtom(a.symbol, a.charge, a.multiplicity, a.mass, tform(a.coords), a.stereo)
end

function (tform::AffineMap)(mol::GraphMol)
    return GraphMol(mol.neighbormap, mol.edges, [tform(a) for a in nodeattrs(mol)], mol.edgeattrs, mol.cache, mol.attributes)
end

function affinetransform(x::MolGMM, tform::AffineMap)
    return MolGMM([tform(g) for g in x.gaussians], tform(x.graph), x.nodes, x.σfun, x.ϕfun)
end

function affinetransform(x::PharmacophoreGMM, tform::AffineMap)
    tformgmms = [tform(x.gmms[key]) for key in keys(x.gmms)]
    gmmdict = Dict{eltype(keys(x.gmms)),eltype(tformgmms)}()
    for (i,key) in enumerate(keys(x.gmms))
        push!(gmmdict, key=>tformgmms[i])
    end
    return PharmacophoreGMM(gmmdict, tform(x.graph), x.nodes, x.σfun, x.ϕfun, x.features, x.directional)
end


text:
mutable struct GLMCopulaVCObs{T <: BlasReal, D, Link}
    # data
    y::Vector{T}
    X::Matrix{T}
    V::Vector{Matrix{T}}
    n::Int
    p::Int          # number of mean parameters in linear regression
    m::Int          # number of variance components
    # working arrays
    ∇β::Vector{T}   # gradient wrt β
    ∇resβ::Matrix{T}# residual gradient matrix d/dβ_p res_ij (each observation has a gradient of residual is px1)
    ∇τ::Vector{T}   # gradient wrt τ
    ∇θ::Vector{T}   # gradient wrt θ
    Hβ::Matrix{T}   # Hessian wrt β
    Hθ::Matrix{T}   # Hessian wrt variance components θ
    Hτ::Matrix{T}   # Hessian wrt τ
    res::Vector{T}  # residual vector res_i
    t::Vector{T}    # t[k] = tr(V_i[k]) / 2
    q::Vector{T}    # q[k] = res_i' * V_i[k] * res_i / 2
    xtx::Matrix{T}  # Xi'Xi
    storage_n::Vector{T}
    m1::Vector{T}
    m2::Vector{T}
    storage_p1::Vector{T}
    storage_p2::Vector{T}
    storage_np::Matrix{T}
    storage_pp::Matrix{T}
    added_term_numerator::Matrix{T}
    added_term2::Matrix{T}
    η::Vector{T}    # η = Xβ systematic component
    μ::Vector{T}    # μ(β) = ginv(Xβ) # inverse link of the systematic component
    varμ::Vector{T} # v(μ_i) # variance as a function of the mean
    dμ::Vector{T}   # derivative of μ
    d::D            # distribution()
    link::Link      # link function ()
    wt::Vector{T}   # weights wt for GLM.jl
    w1::Vector{T}   # working weights in the gradient = dμ/v(μ)
    w2::Vector{T}   # working weights in the information matrix = dμ^2/v(μ)
end

function GLMCopulaVCObs(
    y::Vector{T},
    X::Matrix{T},
    V::Vector{Matrix{T}},
    d::D,
    link::Link) where {T <: BlasReal, D, Link}
    n, p, m = size(X, 1), size(X, 2), length(V)
    @assert length(y) == n "length(y) should be equal to size(X, 1)"
    # working arrays
    ∇β  = Vector{T}(undef, p)
    ∇resβ  = Matrix{T}(undef, n, p)
    ∇τ  = Vector{T}(undef, 1)
    ∇θ  = Vector{T}(undef, m)
    Hβ  = Matrix{T}(undef, p, p)
    Hθ  = Matrix{T}(undef, m, m)
    Hτ  = Matrix{T}(undef, 1, 1)
    res = Vector{T}(undef, n)
    t   = [tr(V[k])/2 for k in 1:m]
    q   = Vector{T}(undef, m)
    xtx = transpose(X) * X
    storage_n = Vector{T}(undef, n)
    m1        = Vector{T}(undef, m)
    m2        = Vector{T}(undef, m)
    storage_p1 = Vector{T}(undef, p)
    storage_p2 = Vector{T}(undef, p)
    storage_np = Matrix{T}(undef, n, p)
    storage_pp = Matrix{T}(undef, p, p)
    added_term_numerator = Matrix{T}(undef, n, p)
    added_term2 = Matrix{T}(undef, p, p)
    η = Vector{T}(undef, n)
    μ = Vector{T}(undef, n)
    varμ = Vector{T}(undef, n)
    dμ = Vector{T}(undef, n)
    wt = Vector{T}(undef, n)
    fill!(wt, one(T))
    w1 = Vector{T}(undef, n)
    w2 = Vector{T}(undef, n)
    # constructor
    GLMCopulaVCObs{T, D, Link}(y, X, V, n, p, m, ∇β, ∇resβ, ∇τ, ∇θ, Hβ, Hθ,
        Hτ, res, t, q, xtx, storage_n, m1, m2, storage_p1, storage_p2, storage_np,
        storage_pp, added_term_numerator, added_term2, η, μ, varμ, dμ, d, link, wt, w1, w2)
end

"""
GLMCopulaVCModel
GLMCopulaVCModel(gcs)
Gaussian copula variance component model, which contains a vector of
`GLMCopulaVCObs` as data, model parameters, and working arrays.
"""
struct GLMCopulaVCModel{T <: BlasReal, D, Link} <: MathProgBase.AbstractNLPEvaluator
    # data
    data::Vector{GLMCopulaVCObs{T, D, Link}}
    Ytotal::T
    ntotal::Int     # total number of singleton observations
    p::Int          # number of mean parameters in linear regression
    m::Int          # number of variance components
    # parameters
    β::Vector{T}    # p-vector of mean regression coefficients
    τ::Vector{T}    # inverse of linear regression variance parameter
    θ::Vector{T}    # m-vector: [θ12, ..., θm2]
    # working arrays
    ∇β::Vector{T}   # gradient from all observations
    ∇τ::Vector{T}
    ∇θ::Vector{T}
    XtX::Matrix{T}  # X'X = sum_i Xi'Xi
    Hβ::Matrix{T}    # Hessian from all observations
    Hθ::Matrix{T}
    Hτ::Matrix{T}
    Ainv::Matrix{T}
    Aevec::Matrix{T}
    M::Matrix{T}
    vcov::Matrix{T}
    ψ::Vector{T}
    TR::Matrix{T}         # n-by-m matrix with tik = tr(Vi[k]) / 2
    QF::Matrix{T}         # n-by-m matrix with qik = res_i' Vi[k] res_i
    storage_n::Vector{T}
    storage_m::Vector{T}
    storage_θ::Vector{T}
    d::Vector{D}
    link::Vector{Link}
    penalized::Bool
end

function GLMCopulaVCModel(gcs::Vector{GLMCopulaVCObs{T, D, Link}}; penalized::Bool = false) where {T <: BlasReal, D<:Union{Poisson, Bernoulli}, Link}
    n, p, m = length(gcs), size(gcs[1].X, 2), length(gcs[1].V)
    β       = Vector{T}(undef, p)
    τ       = [1.0]
    θ       = Vector{T}(undef, m)
    ∇β      = Vector{T}(undef, p)
    ∇τ      = Vector{T}(undef, 1)
    ∇θ      = Vector{T}(undef, m)
    XtX     = zeros(T, p, p) # sum_i xi'xi
    Hβ      = Matrix{T}(undef, p, p)
    Hθ      = Matrix{T}(undef, m, m)
    Hτ      = Matrix{T}(undef, 1, 1)
    Ainv    = zeros(T, p + m, p + m)
    Aevec   = zeros(T, p + m, p + m)
    M       = zeros(T, p + m, p + m)
    vcov    = zeros(T, p + m, p + m)
    ψ       = Vector{T}(undef, p + m)
    TR      = Matrix{T}(undef, n, m) # collect trace terms
    Ytotal  = 0.0
    ntotal  = 0.0
    d       = Vector{D}(undef, n)
    link    = Vector{Link}(undef, n)
    for i in eachindex(gcs)
        ntotal  += length(gcs[i].y)
        Ytotal  += sum(gcs[i].y)
        BLAS.axpy!(one(T), gcs[i].xtx, XtX)
        TR[i, :] = gcs[i].t
        d[i] = gcs[i].d
        link[i] = gcs[i].link
    end
    QF        = Matrix{T}(undef, n, m)
    storage_n = Vector{T}(undef, n)
    storage_m = Vector{T}(undef, m)
    storage_θ = Vector{T}(undef, m)
    GLMCopulaVCModel{T, D, Link}(gcs, Ytotal, ntotal, p, m, β, τ, θ,
        ∇β, ∇τ, ∇θ, XtX, Hβ, Hθ, Hτ, Ainv, Aevec, M, vcov, ψ, TR, QF,
        storage_n, storage_m, storage_θ, d, link, penalized)
end

"""
    loglikelihood!(gc, β, τ, θ)
Calculates the loglikelihood of observing `y` given mean `μ`, for the Poisson and Bernoulli base distribution using the GLM.jl package.
"""
function loglikelihood!(
    gc::Union{GLMCopulaVCObs, Poisson_Bernoulli_VCObs},
    β::Vector{T},
    θ::Vector{T},
    needgrad::Bool = false,
    needhess::Bool = false;
    penalized::Bool = false
    ) where {T <: BlasReal}
    # n, p, m = size(gc.X, 1), size(gc.X, 2), length(gc.V)
    needgrad = needgrad || needhess
    if needgrad
        fill!(gc.∇β, 0)
        fill!(gc.∇θ, 0)
    end
    needhess && fill!(gc.Hβ, 0)
    update_res!(gc, β)
    # @show gc.res
    standardize_res!(gc)
    # @show gc.res
    fill!(gc.∇resβ, 0.0) # fill gradient of residual vector with 0
    std_res_differential!(gc) # this will compute ∇resβ

    # evaluate copula loglikelihood
    @inbounds for k in 1:gc.m
        mul!(gc.storage_n, gc.V[k], gc.res) # storage_n = V[k] * res
        if needgrad
            BLAS.gemv!('T', θ[k], gc.∇resβ, gc.storage_n, 1.0, gc.∇β) # stores ∇resβ*Γ*res (standardized residual)
        end
        gc.q[k] = dot(gc.res, gc.storage_n) / 2
    end
    # loglikelihood
    logl = QuasiCopula.component_loglikelihood(gc)
    tsum = dot(θ, gc.t)
    logl += -log(1 + tsum)
    qsum  = dot(θ, gc.q)
    logl += log(1 + qsum)
    # add L2 ridge penalty
    if penalized
        logl -= 0.5 * dot(θ, θ)
    end
    if needgrad
        inv1pq = inv(1 + qsum)
        inv1pt = inv(1 + tsum) #
        # gc.∇θ .= inv1pq * gc.q .- inv1pt * gc.t
        gc.m1 .= gc.q
        gc.m1 .*= inv1pq
        gc.m2 .= gc.t
        gc.m2 .*= inv1pt
        gc.∇θ .= gc.m1 .- gc.m2
        if penalized
            gc.∇θ .-= θ
        end
        if needhess
            BLAS.syrk!('L', 'N', -abs2(inv1pq), gc.∇β, 1.0, gc.Hβ) # only lower triangular
            # does adding this term to the approximation of the hessian violate negative semidefinite properties?
            fill!(gc.added_term_numerator, 0.0) # fill gradient with 0
            fill!(gc.added_term2, 0.0) # fill hessian with 0
            @inbounds for k in 1:gc.m
                mul!(gc.added_term_numerator, gc.V[k], gc.∇resβ) # storage_n = V[k] * res
                BLAS.gemm!('T', 'N', θ[k], gc.∇resβ, gc.added_term_numerator, one(T), gc.added_term2)
            end
            gc.added_term2 .*= inv1pq
            gc.Hβ .+= gc.added_term2
            gc.Hβ .+= QuasiCopula.glm_hessian(gc)

            # hessian for vc
            fill!(gc.Hθ, 0.0)
            BLAS.syr!('U', one(T), gc.m2, gc.Hθ)
            BLAS.syr!('U', -one(T), gc.m1, gc.Hθ)
            copytri!(gc.Hθ, 'U')
            # gc.Hθ .= gc.m2 * transpose(gc.m2) - gc.m1 * transpose(gc.m1)
        end
        gc.storage_p2 .= gc.∇β .* inv1pq
        # @show gc.res
        gc.res .= gc.y .- gc.μ
        gc.∇β .= QuasiCopula.glm_gradient(gc)
        # @show gc.res
        gc.∇β .+= gc.storage_p2
    end
    logl
end

function loglikelihood!(
    gcm::Union{GLMCopulaVCModel, Poisson_Bernoulli_VCModel},
    needgrad::Bool = false,
    needhess::Bool = false
    )
    logl = 0.0
    if needgrad
        fill!(gcm.∇β, 0)
        fill!(gcm.∇θ, 0)
    end
    if needhess
        fill!(gcm.Hβ, 0)
        fill!(gcm.Hθ, 0)
    end
    logl = zeros(Threads.nthreads())
    Threads.@threads for i in eachindex(gcm.data)
        @inbounds logl[Threads.threadid()] += loglikelihood!(gcm.data[i], gcm.β,
         gcm.θ, needgrad, needhess; penalized = gcm.penalized)
     end
     @inbounds for i in eachindex(gcm.data)
         if needgrad
             gcm.∇β .+= gcm.data[i].∇β
             gcm.∇θ .+= gcm.data[i].∇θ
         end
         if needhess
             gcm.Hβ .+= gcm.data[i].Hβ
             gcm.Hθ .+= gcm.data[i].Hθ
         end
     end
    sum(logl)
end


text:
function updatevolumes!(system::CVSystem,n::Int64)
    system.arterialvolume = 0;
    system.peripheralvolume = 0;
    system.vcvolume = 0;
    system.heartvolume = 0;
    system.lungvolume = 0;
    system.livervolume = 0;

    # 1D arterial domain and peripheral compartments
    for i = 1:length(system.branches.ID)
        x = convert(Array{Float64,1},
            linspace(0,system.branches.lengthincm[i]*cmTom,
            system.solverparams.JL));
        system.arterialvolume += NumericalIntegration.integrate(x,system.branches.A[i][:,n+1]);
        if isempty(system.branches.children[i])
            system.arterialvolume += system.branches.term[i].V[n+1,1];
            system.peripheralvolume += sum(system.branches.term[i].V[n+1,2:5]);
        end
    end
     # liver, vena cava, heart, and lung compartments
    system.vcvolume += (system.svc.V[n+1] + system.ivc.V[n+1]);
    system.heartvolume += (system.heart.rv.V[n+1] + system.heart.ra.V[n+1] +
        system.heart.lv.V[n+1] + system.heart.la.V[n+1]);
    system.lungvolume += (sum(system.lungs.Va[n+1,:]) + sum(system.lungs.Vv[n+1,:]));
    system.livervolume += sum(system.liver.V[n+1,:]);

    if n == 0
        system.initialvolume = (system.arterialvolume +
            system.peripheralvolume + system.vcvolume + system.heartvolume +
            system.lungvolume + system.livervolume);
        println("Blood volume at t = $(system.t[n+1]) s: $(system.initialvolume/cm3Tom3) mL")
    else
        system.finalvolume = (system.arterialvolume +
            system.peripheralvolume + system.vcvolume + system.heartvolume +
            system.lungvolume + system.livervolume);
        println("Blood volume at t = $(system.t[n+1]) s: $(system.finalvolume/cm3Tom3) mL")
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
module VMLinux
    using CBinding

    vmlinux_iob = IOBuffer()
    run(pipeline(`bpftool btf dump file /sys/kernel/btf/vmlinux format c`; stdout=vmlinux_iob))
    vmlinux_str = String(take!(vmlinux_iob))

    vmlinux_path = joinpath(@__DIR__, "vmlinux.h")
    open(vmlinux_path, "w") do io
        write(io, vmlinux_str)
        flush(io)
    end

    c`-std=c99 -fparse-all-comments -I$(@__DIR__)`
    c"#include \"vmlinux.h\""
end


text:
## Primitives: Functions to be implemented for a newsvendor problems

"""
    underage_cost(anp::AbstractNewsvendorProblem)

Get the cost of underage of a newsvendor problem.
"""
underage_cost(anp::AbstractNewsvendorProblem) = error("Underage cost is not
 defined in the concrete type")

"""
    overage_cost(anp::AbstractNewsvendorProblem)

Get the cost of overage of a newsvendor problem.
"""
overage_cost(anp::AbstractNewsvendorProblem) = error("Overage cost is not defined
in the concrete type")

"""
	distr(anp::AbstractNewsvendorProblem)

Get the demand distribution of a newsvendor problem.
"""
distr(anp::AbstractNewsvendorProblem) = error("Demand distribution is not
defined in the concrete type")
##

## Optional primitive

"""
profit_shift(anp::AbstractNewsvendorProblem)

Define how profit is shifted because of fixed cost, penalty, etc.;
defautls to 0.0.
"""
profit_shift(anp::AbstractNewsvendorProblem) = zero(underage_cost(anp))

"""
	q_min(anp::AbstractNewsvendorProblem)

Define minimal feasible quanitity, e.g., due to contractual limitations;
defautls to 0.0.
"""
q_min(anp::AbstractNewsvendorProblem) = zero(underage_cost(anp))

"""
	q_max(anp::AbstractNewsvendorProblem)

Define maximal feasible quanitity, e.g., due to contractual limitations;
defautls to `Inf`.
"""
q_max(anp::AbstractNewsvendorProblem) = Inf
##

## Maximizing order quanitity

"""
	critical_fractile(anp::AbstractNewsvendorProblem)

Compute the critical fractile for a newsvendor problem:
```math
\\textrm{critical fractile} = \\frac{\\textrm{underage cost}}{\\textrm{underage cost} + \\textrm{overage cost}}
```
"""
function critical_fractile(anp::AbstractNewsvendorProblem)
    underage_cost(anp) / (underage_cost(anp) + overage_cost(anp))
end

"""
	q_opt(anp::AbstractNewsvendorProblem; rounded = true)

Compute the quanitity that maximizes expected profit for a newsvendor problem
(i.e., where critical fractile equals in-stock probability). Attempts to solve
```math
F(q_{opt}) = \\textrm{critical fractile}
```
where F is the c.d.f. of the demand distribution.
Returns closest next integer unless `rounded=false`; then, it returns exact real.
Clamps at q_min and q_max.
"""
function q_opt(anp::AbstractNewsvendorProblem; rounded = true)

    if critical_fractile(anp) <= 0
        return q_min(anp)
    end

    if critical_fractile(anp) >= 1
        return min(maximum(distr(anp)), q_max(anp))
    end

    q = clamp(quantile(distr(anp), critical_fractile(anp)), q_min(anp), q_max(anp))

    if rounded == false  return q end

    q = round(Int, q)

    if q_min(anp) < q < q_max(anp)
        return q
    else
        return clamp(q, q_min(anp), q_max(anp))
    end

end

@doc raw"""
	q_scarf(anp::AbstractNewsvendorProblem)

Compute the quanitity that maximizes the minimal expected profit among all
distributions with the same mean and variance. Worst-case solution (Scarf, 1958)
```math
q_{scarf} = \mu + \sigma/2 * (\sqrt{r} - \sqrt{1/r})
```
where
```math
r = \frac{\textrm{underage cost}}{ \textrm{overage cost}}
```
"""
function q_scarf(anp::AbstractNewsvendorProblem)
    μ = mean(distr(anp))
    σ = std(distr(anp))
    r = underage_cost(anp) / overage_cost(anp)
    return μ + σ / 2 * (√r - √(1 / r))
end

##

## Metrics: Expected values depending on order quanitity

@doc raw"""
	leftover(anp::AbstractNewsvendorProblem, q)

Compute expected leftover inventory when stocking quantity q.

```math
E[\textrm{leftover}] = \int_{-\infty}^q (q - x)f(x)dx
```
"""
function leftover(anp::AbstractNewsvendorProblem, q)
    if isinf(q) return Inf end
    return leftover(distr(anp), q)
end

"""
    sales(anp::AbstractNewsvendorProblem, q)

Compute expected sales when stocking quantity q.
"""
sales(anp::AbstractNewsvendorProblem, q) = q - leftover(anp, q)

"""
	lost_sales(anp::AbstractNewsvendorProblem, q)

Compute expected lost sales when stocking quantity q.
"""
lost_sales(anp::AbstractNewsvendorProblem, q) = mean(distr(anp)) - sales(anp, q)

"""
	mismatch_cost(anp::AbstractNewsvendorProblem, q)

Compute expected mismatch cost when stocking quantity q. It is given by
```math
E[\\textrm{mismatch cost}] = \\textrm{underage cost} \\times  E[\\textrm{lost sales}] + \\textrm{overage cost} \\times  E[\\textrm{leftover}]
```
"""
function mismatch_cost(anp::AbstractNewsvendorProblem, q)
    if isinf(q) return 0 end
    return (underage_cost(anp) * lost_sales(anp, q)
            +
            overage_cost(anp) * leftover(anp, q))
end

"""
    profit(anp::AbstractNewsvendorProblem)

Compute expected profit when stocking quantity q_opt.

	profit(anp::AbstractNewsvendorProblem, q)

Compute expected profit when stocking quantity q. It is given by
```math
E[\\textrm{profit}] = \\textrm{underage cost} \\times  \\mu - E[\\textrm{mismatch cost}] +  \\textrm{profit shift}
```
"""
profit(anp::AbstractNewsvendorProblem, q) = (underage_cost(anp) * mean(distr(anp))
                                             -
                                             mismatch_cost(anp, q) + profit_shift(anp))

profit(anp::AbstractNewsvendorProblem) = profit(anp, q_opt(anp))


text:
@testset "5.6.2 Inverse cosecant functions" begin
    (a, b, c, d, n, x, ) = @variables a b c d n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integration*Problems*Involving*Inverse*Cosecants=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*u*acsc(a*x^n)^p=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*acsc(a*x^n)=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [acsc(a*x^5)/x, x, 7, (1/10)*im*acsc(a*x^5)^2 - (1/5)*acsc(a*x^5)*log(1 - ℯ^(2*im*acsc(a*x^5))) + (1/10)*im*PolyLog(2, ℯ^(2*im*acsc(a*x^5)))]

    @test_int [x^3*acsc(sqrt(x)), x, 4, sqrt(-1 + x)/4 + (1/4)*(-1 + x)^(3/2) + (3/20)*(-1 + x)^(5/2) + (1/28)*(-1 + x)^(7/2) + (1/4)*x^4*acsc(sqrt(x))]
    @test_int [x^2*acsc(sqrt(x)), x, 4, sqrt(-1 + x)/3 + (2/9)*(-1 + x)^(3/2) + (1/15)*(-1 + x)^(5/2) + (1/3)*x^3*acsc(sqrt(x))]
    @test_int [x^1*acsc(sqrt(x)), x, 4, sqrt(-1 + x)/2 + (1/6)*(-1 + x)^(3/2) + (1/2)*x^2*acsc(sqrt(x))]
    @test_int [x^0*acsc(sqrt(x)), x, 3, sqrt(-1 + x) + x*acsc(sqrt(x))]
    @test_int [acsc(sqrt(x))/x^1, x, 7, im*acsc(sqrt(x))^2 - 2*acsc(sqrt(x))*log(1 - ℯ^(2*im*acsc(sqrt(x)))) + im*PolyLog(2, ℯ^(2*im*acsc(sqrt(x))))]
    @test_int [acsc(sqrt(x))/x^2, x, 5, -(sqrt(-1 + x)/(2*x)) - acsc(sqrt(x))/x - (1/2)*atan(sqrt(-1 + x))]
    @test_int [acsc(sqrt(x))/x^3, x, 6, -(sqrt(-1 + x)/(8*x^2)) - (3*sqrt(-1 + x))/(16*x) - acsc(sqrt(x))/(2*x^2) - (3/16)*atan(sqrt(-1 + x))]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [x^2*acsc(a/x), x, 5, (1/3)*a^3*sqrt(1 - x^2/a^2) - (1/9)*a^3*(1 - x^2/a^2)^(3/2) + (1/3)*x^3*asin(x/a)]
    @test_int [x^1*acsc(a/x), x, 4, (1/4)*a*x*sqrt(1 - x^2/a^2) - (1/4)*a^2*asin(x/a) + (1/2)*x^2*asin(x/a)]
    @test_int [x^0*acsc(a/x), x, 3, a*sqrt(1 - x^2/a^2) + x*asin(x/a)]
    @test_int [acsc(a/x)/x^1, x, 6, (-(1/2))*im*asin(x/a)^2 + asin(x/a)*log(1 - ℯ^(2*im*asin(x/a))) - (1/2)*im*PolyLog(2, ℯ^(2*im*asin(x/a)))]
    @test_int [acsc(a/x)/x^2, x, 5, -(asin(x/a)/x) - atanh(sqrt(1 - x^2/a^2))/a]
    @test_int [acsc(a/x)/x^3, x, 3, -(sqrt(1 - x^2/a^2)/(2*a*x)) - asin(x/a)/(2*x^2)]
    @test_int [acsc(a/x)/x^4, x, 6, -(sqrt(1 - x^2/a^2)/(6*a*x^2)) - asin(x/a)/(3*x^3) - atanh(sqrt(1 - x^2/a^2))/(6*a^3)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*acsc(a*x^n)*when*n*symbolic=#

    @test_int [acsc(a*x^n)/x, x, 7, (im*acsc(a*x^n)^2)/(2*n) - (acsc(a*x^n)*log(1 - ℯ^(2*im*acsc(a*x^n))))/n + (im*PolyLog(2, ℯ^(2*im*acsc(a*x^n))))/(2*n)]

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*x^m*acsc(a+b*x)^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*acsc(a+b*x)^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [x^4*acsc(a + b*x), x, 9, -((a*(20 + 53*a^2)*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(30*b^5)) - (11*a*x^2*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(60*b^3) + (x^3*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(20*b^2) + ((9 + 58*a^2)*(a + b*x)^2*sqrt(1 - 1/(a + b*x)^2))/(120*b^5) + (a^5*acsc(a + b*x))/(5*b^5) + (1/5)*x^5*acsc(a + b*x) + ((3 + 40*a^2 + 40*a^4)*atanh(sqrt(1 - 1/(a + b*x)^2)))/(40*b^5)]
    @test_int [x^3*acsc(a + b*x), x, 8, ((2 + 17*a^2)*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(12*b^4) + (x^2*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(12*b^2) - (a*(a + b*x)^2*sqrt(1 - 1/(a + b*x)^2))/(3*b^4) - (a^4*acsc(a + b*x))/(4*b^4) + (1/4)*x^4*acsc(a + b*x) - (a*(1 + 2*a^2)*atanh(sqrt(1 - 1/(a + b*x)^2)))/(2*b^4)]
    @test_int [x^2*acsc(a + b*x), x, 7, -((5*a*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(6*b^3)) + (x*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(6*b^2) + (a^3*acsc(a + b*x))/(3*b^3) + (1/3)*x^3*acsc(a + b*x) + ((1 + 6*a^2)*atanh(sqrt(1 - 1/(a + b*x)^2)))/(6*b^3)]
    @test_int [x^1*acsc(a + b*x), x, 6, ((a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(2*b^2) - (a^2*acsc(a + b*x))/(2*b^2) + (1/2)*x^2*acsc(a + b*x) - (a*atanh(sqrt(1 - 1/(a + b*x)^2)))/b^2]
    @test_int [x^0*acsc(a + b*x), x, 5, ((a + b*x)*acsc(a + b*x))/b + atanh(sqrt(1 - 1/(a + b*x)^2))/b]
    @test_int [acsc(a + b*x)/x^1, x, 14, acsc(a + b*x)*log(1 + (im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2))) + acsc(a + b*x)*log(1 + (im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2))) - acsc(a + b*x)*log(1 - ℯ^(2*im*acsc(a + b*x))) - im*PolyLog(2, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2)))) - im*PolyLog(2, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2)))) + (1/2)*im*PolyLog(2, ℯ^(2*im*acsc(a + b*x)))]
    @test_int [acsc(a + b*x)/x^2, x, 6, -((b*acsc(a + b*x))/a) - acsc(a + b*x)/x - (2*b*atan((a - tan((1/2)*acsc(a + b*x)))/sqrt(1 - a^2)))/(a*sqrt(1 - a^2))]
    @test_int [acsc(a + b*x)/x^3, x, 8, -((b*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(2*a*(1 - a^2)*x)) + (b^2*acsc(a + b*x))/(2*a^2) - acsc(a + b*x)/(2*x^2) + ((1 - 2*a^2)*b^2*atan((a - tan((1/2)*acsc(a + b*x)))/sqrt(1 - a^2)))/(a^2*(1 - a^2)^(3/2))]
    @test_int [acsc(a + b*x)/x^4, x, 9, -((b*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(6*a*(1 - a^2)*x^2)) + ((2 - 5*a^2)*b^2*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(6*a^2*(1 - a^2)^2*x) - (b^3*acsc(a + b*x))/(3*a^3) - acsc(a + b*x)/(3*x^3) - ((2 - 5*a^2 + 6*a^4)*b^3*atan((a - tan((1/2)*acsc(a + b*x)))/sqrt(1 - a^2)))/(3*a^3*(1 - a^2)^(5/2))]
    @test_int [acsc(a + b*x)/x^5, x, 10, -((b*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(12*a*(1 - a^2)*x^3)) + ((3 - 8*a^2)*b^2*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(24*a^2*(1 - a^2)^2*x^2) - ((6 - 17*a^2 + 26*a^4)*b^3*(a + b*x)*sqrt(1 - 1/(a + b*x)^2))/(24*a^3*(1 - a^2)^3*x) + (b^4*acsc(a + b*x))/(4*a^4) - acsc(a + b*x)/(4*x^4) + ((2 - 7*a^2 + 8*a^4 - 8*a^6)*b^4*atan((a - tan((1/2)*acsc(a + b*x)))/sqrt(1 - a^2)))/(4*a^4*(1 - a^2)^(7/2))]

    @test_int [x^3*acsc(a + b*x)^2, x, 20, -((a*x)/b^3) + (a + b*x)^2/(12*b^4) + ((a + b*x)*sqrt(1 - 1/(a + b*x)^2)*acsc(a + b*x))/(3*b^4) + (3*a^2*(a + b*x)*sqrt(1 - 1/(a + b*x)^2)*acsc(a + b*x))/b^4 - (a*(a + b*x)^2*sqrt(1 - 1/(a + b*x)^2)*acsc(a + b*x))/b^4 + ((a + b*x)^3*sqrt(1 - 1/(a + b*x)^2)*acsc(a + b*x))/(6*b^4) - (a^4*acsc(a + b*x)^2)/(4*b^4) + (1/4)*x^4*acsc(a + b*x)^2 - (2*a*acsc(a + b*x)*atanh(ℯ^(im*acsc(a + b*x))))/b^4 - (4*a^3*acsc(a + b*x)*atanh(ℯ^(im*acsc(a + b*x))))/b^4 + log(a + b*x)/(3*b^4) + (3*a^2*log(a + b*x))/b^4 + (im*a*PolyLog(2, -ℯ^(im*acsc(a + b*x))))/b^4 + (2*im*a^3*PolyLog(2, -ℯ^(im*acsc(a + b*x))))/b^4 - (im*a*PolyLog(2, ℯ^(im*acsc(a + b*x))))/b^4 - (2*im*a^3*PolyLog(2, ℯ^(im*acsc(a + b*x))))/b^4]
    @test_int [x^2*acsc(a + b*x)^2, x, 17, x/(3*b^2) - (2*a*(a + b*x)*sqrt(1 - 1/(a + b*x)^2)*acsc(a + b*x))/b^3 + ((a + b*x)^2*sqrt(1 - 1/(a + b*x)^2)*acsc(a + b*x))/(3*b^3) + (a^3*acsc(a + b*x)^2)/(3*b^3) + (1/3)*x^3*acsc(a + b*x)^2 + (2*acsc(a + b*x)*atanh(ℯ^(im*acsc(a + b*x))))/(3*b^3) + (4*a^2*acsc(a + b*x)*atanh(ℯ^(im*acsc(a + b*x))))/b^3 - (2*a*log(a + b*x))/b^3 - (im*PolyLog(2, -ℯ^(im*acsc(a + b*x))))/(3*b^3) - (2*im*a^2*PolyLog(2, -ℯ^(im*acsc(a + b*x))))/b^3 + (im*PolyLog(2, ℯ^(im*acsc(a + b*x))))/(3*b^3) + (2*im*a^2*PolyLog(2, ℯ^(im*acsc(a + b*x))))/b^3]
    @test_int [x^1*acsc(a + b*x)^2, x, 11, ((a + b*x)*sqrt(1 - 1/(a + b*x)^2)*acsc(a + b*x))/b^2 - (a^2*acsc(a + b*x)^2)/(2*b^2) + (1/2)*x^2*acsc(a + b*x)^2 - (4*a*acsc(a + b*x)*atanh(ℯ^(im*acsc(a + b*x))))/b^2 + log(a + b*x)/b^2 + (2*im*a*PolyLog(2, -ℯ^(im*acsc(a + b*x))))/b^2 - (2*im*a*PolyLog(2, ℯ^(im*acsc(a + b*x))))/b^2]
    @test_int [x^0*acsc(a + b*x)^2, x, 8, ((a + b*x)*acsc(a + b*x)^2)/b + (4*acsc(a + b*x)*atanh(ℯ^(im*acsc(a + b*x))))/b - (2*im*PolyLog(2, -ℯ^(im*acsc(a + b*x))))/b + (2*im*PolyLog(2, ℯ^(im*acsc(a + b*x))))/b]
    @test_int [acsc(a + b*x)^2/x^1, x, 17, acsc(a + b*x)^2*log(1 + (im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2))) + acsc(a + b*x)^2*log(1 + (im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2))) - acsc(a + b*x)^2*log(1 - ℯ^(2*im*acsc(a + b*x))) - 2*im*acsc(a + b*x)*PolyLog(2, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2)))) - 2*im*acsc(a + b*x)*PolyLog(2, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2)))) + im*acsc(a + b*x)*PolyLog(2, ℯ^(2*im*acsc(a + b*x))) + 2*PolyLog(3, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2)))) + 2*PolyLog(3, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2)))) - (1/2)*PolyLog(3, ℯ^(2*im*acsc(a + b*x)))]
    @test_int [acsc(a + b*x)^2/x^2, x, 12, -((b*acsc(a + b*x)^2)/a) - acsc(a + b*x)^2/x - (2*im*b*acsc(a + b*x)*log(1 + (im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) + (2*im*b*acsc(a + b*x)*log(1 + (im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) - (2*b*PolyLog(2, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2)))))/(a*sqrt(1 - a^2)) + (2*b*PolyLog(2, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2)))))/(a*sqrt(1 - a^2))]

    @test_int [x^2*acsc(a + b*x)^3, x, 25, ((a + b*x)*acsc(a + b*x))/b^3 - (3*im*a*acsc(a + b*x)^2)/b^3 - (3*a*(a + b*x)*sqrt(1 - 1/(a + b*x)^2)*acsc(a + b*x)^2)/b^3 + ((a + b*x)^2*sqrt(1 - 1/(a + b*x)^2)*acsc(a + b*x)^2)/(2*b^3) + (a^3*acsc(a + b*x)^3)/(3*b^3) + (1/3)*x^3*acsc(a + b*x)^3 + (acsc(a + b*x)^2*atanh(ℯ^(im*acsc(a + b*x))))/b^3 + (6*a^2*acsc(a + b*x)^2*atanh(ℯ^(im*acsc(a + b*x))))/b^3 + atanh(sqrt(1 - 1/(a + b*x)^2))/b^3 + (6*a*acsc(a + b*x)*log(1 - ℯ^(2*im*acsc(a + b*x))))/b^3 - (im*acsc(a + b*x)*PolyLog(2, -ℯ^(im*acsc(a + b*x))))/b^3 - (6*im*a^2*acsc(a + b*x)*PolyLog(2, -ℯ^(im*acsc(a + b*x))))/b^3 + (im*acsc(a + b*x)*PolyLog(2, ℯ^(im*acsc(a + b*x))))/b^3 + (6*im*a^2*acsc(a + b*x)*PolyLog(2, ℯ^(im*acsc(a + b*x))))/b^3 - (3*im*a*PolyLog(2, ℯ^(2*im*acsc(a + b*x))))/b^3 + PolyLog(3, -ℯ^(im*acsc(a + b*x)))/b^3 + (6*a^2*PolyLog(3, -ℯ^(im*acsc(a + b*x))))/b^3 - PolyLog(3, ℯ^(im*acsc(a + b*x)))/b^3 - (6*a^2*PolyLog(3, ℯ^(im*acsc(a + b*x))))/b^3]
    @test_int [x^1*acsc(a + b*x)^3, x, 16, (3*im*acsc(a + b*x)^2)/(2*b^2) + (3*(a + b*x)*sqrt(1 - 1/(a + b*x)^2)*acsc(a + b*x)^2)/(2*b^2) - (a^2*acsc(a + b*x)^3)/(2*b^2) + (1/2)*x^2*acsc(a + b*x)^3 - (6*a*acsc(a + b*x)^2*atanh(ℯ^(im*acsc(a + b*x))))/b^2 - (3*acsc(a + b*x)*log(1 - ℯ^(2*im*acsc(a + b*x))))/b^2 + (6*im*a*acsc(a + b*x)*PolyLog(2, -ℯ^(im*acsc(a + b*x))))/b^2 - (6*im*a*acsc(a + b*x)*PolyLog(2, ℯ^(im*acsc(a + b*x))))/b^2 + (3*im*PolyLog(2, ℯ^(2*im*acsc(a + b*x))))/(2*b^2) - (6*a*PolyLog(3, -ℯ^(im*acsc(a + b*x))))/b^2 + (6*a*PolyLog(3, ℯ^(im*acsc(a + b*x))))/b^2]
    @test_int [x^0*acsc(a + b*x)^3, x, 10, ((a + b*x)*acsc(a + b*x)^3)/b + (6*acsc(a + b*x)^2*atanh(ℯ^(im*acsc(a + b*x))))/b - (6*im*acsc(a + b*x)*PolyLog(2, -ℯ^(im*acsc(a + b*x))))/b + (6*im*acsc(a + b*x)*PolyLog(2, ℯ^(im*acsc(a + b*x))))/b + (6*PolyLog(3, -ℯ^(im*acsc(a + b*x))))/b - (6*PolyLog(3, ℯ^(im*acsc(a + b*x))))/b]
    @test_int [acsc(a + b*x)^3/x^1, x, 20, acsc(a + b*x)^3*log(1 + (im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2))) + acsc(a + b*x)^3*log(1 + (im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2))) - acsc(a + b*x)^3*log(1 - ℯ^(2*im*acsc(a + b*x))) - 3*im*acsc(a + b*x)^2*PolyLog(2, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2)))) - 3*im*acsc(a + b*x)^2*PolyLog(2, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2)))) + (3/2)*im*acsc(a + b*x)^2*PolyLog(2, ℯ^(2*im*acsc(a + b*x))) + 6*acsc(a + b*x)*PolyLog(3, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2)))) + 6*acsc(a + b*x)*PolyLog(3, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2)))) - (3/2)*acsc(a + b*x)*PolyLog(3, ℯ^(2*im*acsc(a + b*x))) + 6*im*PolyLog(4, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2)))) + 6*im*PolyLog(4, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2)))) - (3/4)*im*PolyLog(4, ℯ^(2*im*acsc(a + b*x)))]
    @test_int [acsc(a + b*x)^3/x^2, x, 14, -((b*acsc(a + b*x)^3)/a) - acsc(a + b*x)^3/x - (3*im*b*acsc(a + b*x)^2*log(1 + (im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) + (3*im*b*acsc(a + b*x)^2*log(1 + (im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2))))/(a*sqrt(1 - a^2)) - (6*b*acsc(a + b*x)*PolyLog(2, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2)))))/(a*sqrt(1 - a^2)) + (6*b*acsc(a + b*x)*PolyLog(2, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2)))))/(a*sqrt(1 - a^2)) - (6*im*b*PolyLog(3, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 - sqrt(1 - a^2)))))/(a*sqrt(1 - a^2)) + (6*im*b*PolyLog(3, -((im*a*ℯ^(im*acsc(a + b*x)))/(1 + sqrt(1 - a^2)))))/(a*sqrt(1 - a^2))]

    #= ::Subsubsection:: =#
    #=n<0=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*u*acsc(a +b*x^n)=#

    @test_int [x^3*acsc(a + b*x^4), x, 6, ((a + b*x^4)*acsc(a + b*x^4))/(4*b) + atanh(sqrt(1 - 1/(a + b*x^4)^2))/(4*b)]

    @test_int [x^(n-1)*acsc(a + b*x^n), x, 6, ((a + b*x^n)*acsc(a + b*x^n))/(b*n) + atanh(sqrt(1 - 1/(a + b*x^n)^2))/(b*n)]

    #= ::Section::Closed:: =#
    #=Integrands*involving*inverse*cosecant*functions*of*exponentials=#

    @test_int [acsc(c*ℯ^(a + b*x)), x, 7, (im*acsc(c*ℯ^(a + b*x))^2)/(2*b) - (acsc(c*ℯ^(a + b*x))*log(1 - ℯ^(2*im*acsc(c*ℯ^(a + b*x)))))/b + (im*PolyLog(2, ℯ^(2*im*acsc(c*ℯ^(a + b*x)))))/(2*b)]

    #= ::Section::Closed:: =#
    #=Integrands*involving*exponentials*of*inverse*cosecant*functions=#

    @test_int [ℯ^acsc(a*x)*x^2, x, 6, ((4/5 - (12*im)/5)*ℯ^((1 + 3*im)*acsc(a*x))*HypergeometricFunctions._₂F₁(3/2 - im/2, 3, 5/2 - im/2, ℯ^(2*im*acsc(a*x))))/a^3 - ((8/5 - (24*im)/5)*ℯ^((1 + 3*im)*acsc(a*x))*HypergeometricFunctions._₂F₁(3/2 - im/2, 4, 5/2 - im/2, ℯ^(2*im*acsc(a*x))))/a^3]
    @test_int [ℯ^acsc(a*x)*x^1, x, 6, ((8/5 + (4*im)/5)*ℯ^((1 + 2*im)*acsc(a*x))*HypergeometricFunctions._₂F₁(1 - im/2, 2, 2 - im/2, ℯ^(2*im*acsc(a*x))))/a^2 - ((16/5 + (8*im)/5)*ℯ^((1 + 2*im)*acsc(a*x))*HypergeometricFunctions._₂F₁(1 - im/2, 3, 2 - im/2, ℯ^(2*im*acsc(a*x))))/a^2]
    @test_int [ℯ^acsc(a*x)*x^0, x, 5, -(((1 - im)*ℯ^((1 + im)*acsc(a*x))*HypergeometricFunctions._₂F₁(1/2 - im/2, 1, 3/2 - im/2, ℯ^(2*im*acsc(a*x))))/a) + ((2 - 2*im)*ℯ^((1 + im)*acsc(a*x))*HypergeometricFunctions._₂F₁(1/2 - im/2, 2, 3/2 - im/2, ℯ^(2*im*acsc(a*x))))/a]
    @test_int [ℯ^acsc(a*x)/x^1, x, 6, (-im)*ℯ^acsc(a*x) + 2*im*ℯ^acsc(a*x)*HypergeometricFunctions._₂F₁(-(im/2), 1, 1 - im/2, ℯ^(2*im*acsc(a*x)))]
    @test_int [ℯ^acsc(a*x)/x^2, x, 3, (-(1/2))*a*ℯ^acsc(a*x)*sqrt(1 - 1/(a^2*x^2)) - ℯ^acsc(a*x)/(2*x)]
    @test_int [ℯ^acsc(a*x)/x^3, x, 5, (1/5)*a^2*ℯ^acsc(a*x)*cos(2*acsc(a*x)) - (1/10)*a^2*ℯ^acsc(a*x)*sin(2*acsc(a*x))]
    @test_int [ℯ^acsc(a*x)/x^4, x, 6, (-(1/8))*a^3*ℯ^acsc(a*x)*sqrt(1 - 1/(a^2*x^2)) - (a^2*ℯ^acsc(a*x))/(8*x) + (1/40)*a^3*ℯ^acsc(a*x)*cos(3*acsc(a*x)) + (3/40)*a^3*ℯ^acsc(a*x)*sin(3*acsc(a*x))]
    @test_int [ℯ^acsc(a*x)/x^5, x, 6, (1/10)*a^4*ℯ^acsc(a*x)*cos(2*acsc(a*x)) - (1/34)*a^4*ℯ^acsc(a*x)*cos(4*acsc(a*x)) - (1/20)*a^4*ℯ^acsc(a*x)*sin(2*acsc(a*x)) + (1/136)*a^4*ℯ^acsc(a*x)*sin(4*acsc(a*x))]

    #= ::Section::Closed:: =#
    #=Miscellaneous*integrands*involving*inverse*cosecants=#

    @test_int [acsc(a + b*x)/((a*d)/b + d*x), x, 8, (im*acsc(a + b*x)^2)/(2*d) - (acsc(a + b*x)*log(1 - ℯ^(2*im*acsc(a + b*x))))/d + (im*PolyLog(2, ℯ^(2*im*acsc(a + b*x))))/(2*d)]
    end

max_stars_repo_path: test/5 Inverse trig functions/5.6 Inverse cosecant/5.6.2 Inverse cosecant functions.jl


text:
#    Problem : GROUP A
#    *********
#	 Shanno's TRIDIA quadratic tridiagonal problem
#
#    Origonal SIF Source:problem 8 in
#    Ph.L. Toint,
#    "Test problems for partially separable optimization and results
#    for the routine PSPMIN",
#    Report 83/4, Department of Mathematics, FUNDP (Namur, B), 1983.
#
#    See also Buckley#40 (p.96)
#
#    TRIDIA.SIF classification QUR2-AN-V-0
#
#	 This problem is decomposed in n linear groups, the last n-1 of which
#    are 2 x 2 and singular.
#
#    Number of variables is variable

f = (x) -> begin
	fx   = (x[1]-1)^2
	for i in firstindex(x)+1:lastindex(x)
		item    = 2x[i]-x[i-1]
    	fx     += item^2*i
	end
    return fx
end

g! = (g, x) -> begin
	g[1] = 2(x[1]-1)
	for i in firstindex(x)+1:lastindex(x)
		item    = 2x[i]-x[i-1]
        g[i]   += 4item*i
	    g[i-1] -= 2item*i
	end
    return g
end

fg! = (g, x) -> begin
	fx   = (x[1]-1)^2
	g[1] = 2(x[1]-1)
	for i in firstindex(x)+1:lastindex(x)
		item    = 2x[i]-x[i-1]
    	fx     += item^2*i
        g[i]   += 4item*i
	    g[i-1] -= 2item*i
	end
    return fx, g
end

init = (n::Int=5000) -> begin
    return n, ones(n)
end

TestSet["TRIDIA"] = UncProgram("TRIDIA",f, g!, fg!, init)


text:
function _precompile_()
    stackframe(func, file, line; C=false) = StackFrame(Symbol(func), Symbol(file), line, nothing, C, false, 0)
    backtraces = UInt64[0, 4, 3, 2, 1,   # order: calles then caller
                        0, 6, 5, 1,
                        0, 8, 7,
                        0, 4, 3, 2, 1,
                        0]
    if isdefined(Profile, :add_fake_meta)
        backtraces = Profile.add_fake_meta(backtraces)
    end
    lidict = Dict{UInt64,StackFrame}(1=>stackframe(:f1, :file1, 1),
                                     2=>stackframe(:f2, :file1, 5),
                                     3=>stackframe(:f3, :file2, 1),
                                     4=>stackframe(:f2, :file1, 15),
                                     5=>stackframe(:f4, :file1, 20),
                                     6=>stackframe(:f5, :file3, 1),
                                     7=>stackframe(:f1, :file1, 2),
                                     8=>stackframe(:f6, :file3, 10))
    g = flamegraph(backtraces; lidict=lidict)
    win, c, fdraw = viewgui(FlameGraphs.default_colors, g)
    for obs in c.preserved
        if isa(obs, Observable) || isa(obs, Observables.ObserverFunction)
            precompile(obs)
        end
    end
    precompile(fdraw)
    closeall()   # necessary to prevent serialization of stale references (including the internal `empty!`)
    precompile(viewprof_func, (FlameColors, GtkObservables.Canvas{GtkObservables.UserUnit}, FlameGraphs.LeftChildRightSiblingTrees.Node{FlameGraphs.NodeData}, Int))
    precompile(Tuple{typeof(save_as_cb),Ptr{GObject},Tuple{Gtk.GtkCanvas, Vector{UInt64}, Dict{UInt64, Vector{StackFrame}}}})   # time: 0.008177923
    precompile(Tuple{typeof(save_as_cb),Ptr{GObject},Tuple{Gtk.GtkCanvas, Nothing, Nothing}})
    precompile(warntype_last, ())
end


text:
using ModelingToolkit
using Test

@parameters t
@variables a b(t) c d

b = ParentScope(b)
c = ParentScope(ParentScope(c))
d = GlobalScope(d)

renamed(nss, sym) = ModelingToolkit.getname(foldr(ModelingToolkit.renamespace, nss, init=sym))

@test renamed([:foo :bar :baz], a) == :foo₊bar₊baz₊a
@test renamed([:foo :bar :baz], b) == :foo₊bar₊b
@test renamed([:foo :bar :baz], c) == :foo₊c
@test renamed([:foo :bar :baz], d) == :d

eqs = [
    0 ~ a
    0 ~ b
    0 ~ c
    0 ~ d
]
@named sub4 = NonlinearSystem(eqs, [a, b, c, d], [])
@named sub3 = NonlinearSystem(eqs, [a, b, c, d], [])
@named sub2 = NonlinearSystem([], [], [], systems=[sub3, sub4])
@named sub1 = NonlinearSystem([], [], [], systems=[sub2])
@named sys = NonlinearSystem([], [], [], systems=[sub1])

names = ModelingToolkit.getname.(states(sys))
@test :d in names
@test :sub1₊c in names
@test :sub1₊sub2₊b in names
@test :sub1₊sub2₊sub3₊a in names
@test :sub1₊sub2₊sub4₊a in names