text:
using DiffEqBase, DiffEqOperators, ForwardDiff, LinearAlgebra, Test
const A = rand(300,300)
f(du,u) = mul!(du,A,u)
f(u) = A*u
x = rand(300)
v = rand(300)
du = similar(x)

cache1 = ForwardDiff.Dual{DiffEqOperators.JacVecTag}.(x, v)
cache2 = ForwardDiff.Dual{DiffEqOperators.JacVecTag}.(x, v)
@test DiffEqOperators.num_jacvec!(du, f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6
@test DiffEqOperators.num_jacvec!(du, f, x, v, similar(v), similar(v)) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6
@test DiffEqOperators.num_jacvec(f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v rtol=1e-6

@test DiffEqOperators.auto_jacvec!(du, f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v
@test DiffEqOperators.auto_jacvec!(du, f, x, v, cache1, cache2) ≈ ForwardDiff.jacobian(f,similar(x),x)*v
@test DiffEqOperators.auto_jacvec(f, x, v) ≈ ForwardDiff.jacobian(f,similar(x),x)*v

f(du,u,p,t) = mul!(du,A,u)
f(u,p,t) = A*u
L = JacVecOperator(f,x)
@test L*x ≈ DiffEqOperators.auto_jacvec(f, x, x)
@test L*v ≈ DiffEqOperators.auto_jacvec(f, x, v)
@test mul!(du,L,v) ≈ DiffEqOperators.auto_jacvec(f, x, v)
DiffEqBase.update_coefficients!(L,v,nothing,nothing)
@test mul!(du,L,v) ≈ DiffEqOperators.auto_jacvec(f, v, v)

L = JacVecOperator(f,x,autodiff=false)
DiffEqBase.update_coefficients!(L,x,nothing,nothing)
@test L*x ≈ DiffEqOperators.num_jacvec(f, x, x)
@test L*v ≈ DiffEqOperators.num_jacvec(f, x, v)
@test mul!(du,L,v) ≈ DiffEqOperators.num_jacvec(f, x, v)
DiffEqBase.update_coefficients!(L,v,nothing,nothing)
@test mul!(du,L,v) ≈ DiffEqOperators.num_jacvec(f, v, v)

L2 = JacVecOperator{Float64}(f)
DiffEqBase.update_coefficients!(L2,x,nothing,nothing)
@test L2*x ≈ DiffEqOperators.auto_jacvec(f, x, x)
@test L2*v ≈ DiffEqOperators.auto_jacvec(f, v, v)
@test mul!(du,L2,x) ≈ DiffEqOperators.auto_jacvec(f, x, x)

L2 = JacVecOperator{Float64}(f,autodiff=false)
DiffEqBase.update_coefficients!(L2,x,nothing,nothing)
@test L2*x ≈ DiffEqOperators.num_jacvec(f, x, x)
@test L2*v ≈ DiffEqOperators.num_jacvec(f, v, v) rtol=1e-6

using OrdinaryDiffEq
function lorenz(du,u,p,t)
 du[1] = 10.0(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
ff = ODEFunction(lorenz,jac_prototype=JacVecOperator{Float64}(lorenz,u0))
prob = ODEProblem(ff,u0,tspan)
sol = solve(prob,Rosenbrock23())
sol = solve(prob,Rosenbrock23(linsolve=LinSolveGMRES(tol=1e-10)))

ff = ODEFunction(lorenz,jac_prototype=JacVecOperator{Float64}(lorenz,u0,autodiff=false))
prob = ODEProblem(ff,u0,tspan)
sol = solve(prob,Rosenbrock23())
sol = solve(prob,Rosenbrock23(linsolve=LinSolveGMRES(tol=1e-10)))


text:
function single_stance!(du, u, p, t)
    z = u
    walker = p
    θ₁ = z[1]
    θ₂ = z[3]
    ω₁ = z[2]
    ω₂ = z[4]
    M = walker.M; m = walker.m; I=walker.I; l = walker.l; c = walker.c;
    g = walker.g; γ = walker.γ

    A11 = 2*I + M*l^2 + 2*c^2*m + 2*l^2*m - 2*c*l*m - 2*c*l*m*cos(θ₂);
    A12 = I + c^2*m - c*l*m*cos(θ₂);
    A21 = I + c^2*m - c*l*m*cos(θ₂);
    A22 = I + c^2*m;
    A_ss = [A11 A12; A21 A22];

    b1 = c*g*m*sin(γ - θ₁) - M*g*l*sin(γ - θ₁) - c*g*m*sin(θ₁ - γ + θ₂) - 2*g*l*m*sin(γ - θ₁) - c*l*m*ω₂^2*sin(θ₂) - 2*c*l*m*ω₁*ω₂*sin(θ₂);
    b2 = -c*m*(g*sin(θ₁ - γ + θ₂) - l*ω₁^2*sin(θ₂));
    b_ss = [b1; b2];

    if walker.use_controller
        (U, θ₂ᵣ, θ₂̇ᵣ, θ₂̈ᵣ) = controller(t, z, walker)
    else
        U = 0.0
    end

    B = [0; 1]
    α = A_ss\(b_ss+B*U);
    du[1] = ω₁
    du[2] = α[1]
    du[3] = ω₂
    du[4] = α[2]

end

function footstrike(t, z, walker)
    z=z[1]
    θ₁_n = z[1]
    θ₂_n = z[3]
    ω₁_n = z[2]
    ω₂_n = z[4]

    θ₁ = θ₁_n + θ₂_n
    θ₂ = -θ₂_n

    M = walker.M; m = walker.m; I = walker.I
    l = walker.l; c=walker.c

    J11 = 1;
    J12 = 0;
    J13 = l*(cos(θ₁_n + θ₂_n) - cos(θ₁_n));
    J14 = l*cos(θ₁_n + θ₂_n);
    J21 = 0;
    J22 = 1;
    J23 = l*(sin(θ₁_n + θ₂_n) - sin(θ₁_n));
    J24 = l*sin(θ₁_n + θ₂_n);
    J = [J11 J12 J13 J14; J21 J22 J23 J24];

    A11 = M + 2*m;
    A12 = 0;
    A13 = (m*(2*c*cos(θ₁_n + θ₂_n) - 2*l*cos(θ₁_n)))/2 + m*cos(θ₁_n)*(c - l) - M*l*cos(θ₁_n);
    A14 = c*m*cos(θ₁_n + θ₂_n);
    A21 = 0;
    A22 = M + 2*m;
    A23 = (m*(2*c*sin(θ₁_n + θ₂_n) - 2*l*sin(θ₁_n)))/2 - M*l*sin(θ₁_n) + m*sin(θ₁_n)*(c - l);
    A24 = c*m*sin(θ₁_n + θ₂_n);
    A31 = (m*(2*c*cos(θ₁_n + θ₂_n) - 2*l*cos(θ₁_n)))/2 + m*cos(θ₁_n)*(c - l) - M*l*cos(θ₁_n);
    A32 = (m*(2*c*sin(θ₁_n + θ₂_n) - 2*l*sin(θ₁_n)))/2 - M*l*sin(θ₁_n) + m*sin(θ₁_n)*(c - l);
    A33 = 2*I + M*l^2 + 2*c^2*m + 2*l^2*m - 2*c*l*m - 2*c*l*m*cos(θ₂_n);
    A34 = I + c^2*m - c*l*m*cos(θ₂_n);
    A41 = c*m*cos(θ₁_n + θ₂_n);
    A42 = c*m*sin(θ₁_n + θ₂_n);
    A43 = I + c^2*m - c*l*m*cos(θ₂_n);
    A44 = I + c^2*m;
    A_n_hs = [A11 A12 A13 A14; A21 A22 A23 A24; A31 A32 A33 A34; A41 A42 A43 A44];

    X_n_hs = [0 0 ω₁_n ω₂_n]';
    b_hs = [A_n_hs*X_n_hs; 0; 0];
    A_hs = [A_n_hs -J' ; J zeros((2,2))];
    X_hs = A_hs\b_hs;
    ω = zeros(2)
    ω[1] = X_hs[3]+X_hs[4]; ω[2] = -X_hs[4];

    zplus = [θ₁ ω[1] θ₂ ω[2]]

    return zplus
end

ext: jl


text:
ismutable(el) = ismutable(typeof(el))
ismutable(::Type) = Val(false)
ismutable(::Type{<:Array}) = Val(true)

import Base: getindex, setindex!, length, copy, vcat, keys, values, iterate
import Base: zero

import Base: valtype

const IndexedTime = Tuple{Int64,Float64}

"""
    ContinuousTimeProcess{T}

Types inheriting from the abstract type `ContinuousTimeProcess{T}` characterize
the properties of a `T`-valued stochastic process, play
a similar role as distribution types like `Exponential` in the package
`Distributions`.
"""
abstract type ContinuousTimeProcess{T} end
const ProcessOrCoefficients = Union{ContinuousTimeProcess,Tuple{Function,Function}}

"""
    a(t, x, P::ProcessOrCoefficients)

Fallback for `a(t, x, P)` calling `σ(t, x, P)*σ(t, x, P)'`.
"""
a(t, x, P::ProcessOrCoefficients) = outer(σ(t, x, P))
Γ(t, x, P::ProcessOrCoefficients) = inv(a(t, x, P))

abstract type AbstractPath{T} end

"""
    valtype(::ContinuousTimeProcess) -> T

Returns statespace (type) of a `ContinuousTimeProcess{T]`.
"""
valtype(::ContinuousTimeProcess{T}) where {T} = T
valtype(::AbstractPath{T}) where {T} = T

"""
    outertype(P::ContinuousTimeProcess) -> T

Returns the type of `outer(x)`, where `x` is a state of `P`
"""
outertype(::ContinuousTimeProcess{Float64}) = Float64
outertype(P::ContinuousTimeProcess{<:StaticArray}) = typeof(outer(zero(valtype(P))))

"""
    SamplePath{T} <: AbstractPath{T}

The struct
```
struct SamplePath{T}
    tt::Vector{Float64}
    yy::Vector{T}
    SamplePath{T}(tt, yy) where {T} = new(tt, yy)
end
```
serves as container for discretely observed `ContinuousTimeProcess`es and for the sample path returned
by direct and approximate samplers. `tt` is the vector of the grid points of the observation/simulation
and `yy` is the corresponding vector of states.

It supports `getindex, setindex!, length, copy, vcat`.
"""
struct SamplePath{T} <: AbstractPath{T}
    tt::Vector{Float64}
    yy::Vector{T}
    SamplePath{T}(tt, yy) where {T} = new(tt, yy)
end
SamplePath(tt, yy::Vector{T}) where {T} = SamplePath{T}(tt, yy)

samplepath(tt, v) = samplepath(tt, v, ismutable(v))

samplepath(tt, v, ::Val{false}) = SamplePath(tt, fill(v, length(tt)))
samplepath(tt, v, ::Val{true}) = SamplePath(tt, [copy(v) for t in tt])

copy(X::SamplePath{T}) where {T} = SamplePath{T}(copy(X.tt), copy(X.yy))
length(X::SamplePath) = length(X.tt)

getindex(V::SamplePath, I::AbstractArray) = SamplePath(V.tt[I], V.yy[I])
getindex(V::SamplePath, i::Integer) = V.tt[i] => V.yy[i]
vcat(Ys::SamplePath{T}...) where {T} = SamplePath{T}(vcat(map(Y -> Y.tt, Ys)...), vcat(map(Y -> Y.yy, Ys)...))

# iterator
# start(Y::SamplePath) = 1
# next(Y::SamplePath, state) = Y[state], state + 1
# done(V::SamplePath, state) = state > endof(V)

# lastindex(V::SamplePath) = lastindex(V.tt)

function setindex!(V::SamplePath, y, I)
    V.tt[I], V.yy[I] = y
end

"""
    endpoint!(X::SamplePath, v)

Convenience functions setting the endpoint of `X to `v`.
"""
function endpoint!(X::SamplePath, v)
    X.yy[end] = v
    X
end

keys(X::SamplePath) = X.tt
values(X::SamplePath) = X.yy
SamplePath(X::Vector{Pair{Float64,T}}) where {T} = SamplePath{T}(map(first, X), map(last, X))

zero(X::SamplePath) = SamplePath(X.tt, X.yy)

import Base.broadcast
broadcast(f, X::SamplePath, Y::SamplePath) = (X.tt != Y.tt) ? throw(ArgumentError("differing sample times")) : SamplePath(X.tt, f.(X.yy, Y.yy))

struct VSamplePath{T} <: Bridge.AbstractPath{T}
    tt::Vector{Float64}
    yy::Matrix{T}
    function VSamplePath(tt, yy::Matrix{T}) where {T}
        length(tt) != size(yy, 2) && throw(DimensionMismatch("length(tt) != size(yy, 2)"))
        new{T}(tt, yy)
    end
end

length(X::VSamplePath) = length(X.tt)

allsametime(xx) = all(x -> x.tt == xx[1].tt, xx)
function stack(args::SamplePath...)
    assert(allsametime(args))
    VSamplePath(args[1].tt, vcat((X.yy' for X in args)...))
end

# separate a zip of two Vectors
if VERSION < v"1.1-"
sep(Z::Base.Iterators.Zip2{Vector{T1},Vector{T2}}) where {T1,T2} =
    T1[z[1] for z in Z], T2[z[2] for z in Z] # takes into account the minimum of length
else
sep(Z::Base.Iterators.Zip{Tuple{Vector{T1},Vector{T2}}}) where {T1,T2} =
    T1[z[1] for z in Z], T2[z[2] for z in Z] # takes into account the minimum of length
end

copy(X::VSamplePath) = VSamplePath(copy(X.tt), copy(X.yy))

"""
Like `VSamplePath`, but with assumptions on `tt` and dimensionality.
Planned replacement for `VSamplePath`
"""
struct GSamplePath{S,T,P} <: Bridge.AbstractPath{T}
    tt::S
    yy::P
    GSamplePath(tt::S, yy::P) where {S, P<:AbstractArray{T}} where {T} = new{S,T,P}(tt,yy)
end

import Base: length
length(X::GSamplePath) = length(X.tt)

"""
    Increments{S<:AbstractPath{T}}

Iterator over the increments of an AbstractPath.
Iterates over `(i, tt[i], tt[i+1]-tt[i], yy[i+1]-y[i])`.
"""
mutable struct Increments{S<:AbstractPath}
    X::S
end

iterate(dX::Increments, i = 1) = i + 1 > length(dX.X.tt) ? nothing : ((i, dX.X.tt[i], dX.X.tt[i+1]-dX.X.tt[i], dX.X.yy[.., i+1]-dX.X.yy[.., i]), i + 1)

increments(X::AbstractPath) = Increments(X)

# Interoperatibility SDEs

b!(t, u, du, fg::Tuple{Function,Function}) = fg[1](t, u, du)
σ!(t, u, dw, dm, fg::Tuple{Function,Function}) = fg[2](t, u, dw, dm)
b(t, x, fg::Tuple{Function,Function}) = fg[1](t, x)
σ(t, x, fg::Tuple{Function,Function}) = fg[2](t, x)

_b!((i,t)::IndexedTime, u, du, fg::Tuple{Function,Function}) = fg[1](t, u, du)
_b((i,t)::IndexedTime, x, fg::Tuple{Function,Function}) = fg[1](t, x)

# Interoperatibility ODEs

F(t, x, f::Tuple{Function}) = f[1](t, x)
#@inline _F(t, x, F::Function) = F(t, x)
#@inline _F(t, x, F::Tuple{Function,Function}) = F[1](t)*x + F[2](t)


text:
# ---
# title: 676. Implement Magic Dictionary
# id: problem676
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Trie
# link: <https://leetcode.com/problems/implement-magic-dictionary/description/>
# hidden: true
# ---
#
# Design a data structure that is initialized with a list of **different**
# words. Provided a string, you should determine if you can change exactly one
# character in this string to match any word in the data structure.
#
# Implement the `MagicDictionary` class:
#
#   * `MagicDictionary()` Initializes the object.
#   * `void buildDict(String[] dictionary)` Sets the data structure with an array of distinct strings `dictionary`.
#   * `bool search(String searchWord)` Returns `true` if you can change **exactly one character** in `searchWord` to match any string in the data structure, otherwise returns `false`.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["MagicDictionary", "buildDict", "search", "search", "search", "search"]
#     [[], [["hello", "leetcode"]], ["hello"], ["hhllo"], ["hell"], ["leetcoded"]]
#     **Output**
#     [null, null, false, true, false, false]
#
#     **Explanation**
#     MagicDictionary magicDictionary = new MagicDictionary();
#     magicDictionary.buildDict(["hello", "leetcode"]);
#     magicDictionary.search("hello"); // return False
#     magicDictionary.search("hhllo"); // We can change the second 'h' to 'e' to match "hello" so we return True
#     magicDictionary.search("hell"); // return False
#     magicDictionary.search("leetcoded"); // return False
#
#
#
#
# **Constraints:**
#
#   * `1 <= dictionary.length <= 100`
#   * `1 <= dictionary[i].length <= 100`
#   * `dictionary[i]` consists of only lower-case English letters.
#   * All the strings in `dictionary` are **distinct**.
#   * `1 <= searchWord.length <= 100`
#   * `searchWord` consists of only lower-case English letters.
#   * `buildDict` will be called only once before `search`.
#   * At most `100` calls will be made to `search`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module DeepViz

using PlotUtils
using Plots
using Zygote
using Zygote: hook, @adjoint
using Flux
using StatsBase
using Augmentor
using Images
using Metalhead

include("utils.jl")
include("vanilla_gradient.jl")
include("gradcam.jl")
include("guided_backprop.jl")
include("smooth_grad.jl")
include("image_generation.jl")

end


text:
module MonomialBasesBenchs

using InplaceArrays.TensorValues
using InplaceArrays.Fields
using InplaceArrays.Polynomials

@inline function repeat(n,f,args...)
  for i in 1:n
    f(args...)
  end
  nothing
end

function bench1(n)

  xi = Point(2,3)
  np = 5
  x = fill(xi,np)

  orders = (1,2)
  V = Float64
  b = MonomialBasis{2}(V,orders)

  cb = field_cache(b,x)
  @time repeat(n,evaluate_field!,cb,b,x)

  ∇b = ∇(b)
  c∇b = field_cache(∇b,x)
  @time repeat(n,evaluate_field!,c∇b,∇b,x)

end

function bench2(n)

  xi = Point(2,3)
  np = 5
  x = fill(xi,np)

  orders = (1,2)
  V = VectorValue{2,Float64}
  b = MonomialBasis{2}(V,orders)

  cb = field_cache(b,x)
  @time repeat(n,evaluate_field!,cb,b,x)

  ∇b = ∇(b)
  c∇b = field_cache(∇b,x)
  @time repeat(n,evaluate_field!,c∇b,∇b,x)

end

for n in (1,1,10,1000,100000)
  @eval begin
    println("+++ runing suite for n = $($n) +++")
    bench1($n)
    bench2($n)
  end
end

end # module


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

module RunCommandsClient

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

include("modelincludes.jl")

include("api_VirtualMachineRunCommandsApi.jl")
include("api_VirtualMachineScaleSetVMsApi.jl")
include("api_VirtualMachinesApi.jl")

# export models
export convert, InstanceViewStatus
export convert, RunCommandDocument
export convert, RunCommandDocumentBase
export convert, RunCommandInput
export convert, RunCommandInputParameter
export convert, RunCommandListResult
export convert, RunCommandParameterDefinition
export convert, RunCommandResult

# export operations
export convert, VirtualMachineRunCommandsApi, VirtualMachineScaleSetVMsApi, VirtualMachinesApi

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

end


text:
#!/usr/bin/env julia

# A Julia script that reproduces expected results of the test module
# 'test/rationalTest.cpp'.
#
# From a shell, run the script as:
#   julia /path/to/rational.jl
#
# From Julia, the script can be run as:
#   include("path/to/rational.jl")

function rationalTest()
    s = Rational(-12)
    println("Rational(-12): ", s)
    s = 34824//1000
    println("Rational('34.824'): ", s)
    s = (-4285714+4) // 999999
    println("-4.|285714 = ", s)
    println("abs(", s, ") = ", abs(s))
    s = (35167-35) // (10000-10)
    println("3.5|167 = ", s)
    println("abs(", s, ") = ", abs(s))

    println()
    a = 15//20
    println("a = ", a, " = ", num(a)*20, "//", den(a)*20)
    b = a
    c = Rational(0)
    print("b = ", b, "    c = ", c)
    c = a
    c1 = Rational(-7)
    println("    c = ", c, "   c = ", c1)

    println()
    b = inv(a)
    c = 2//5
    println("a = ", a, "    b = ", b, "    c = ", c, "    c^(-1) = ", inv(c))

    a = 1//3
    b = 1//2
    c = +a + b
    println(a, " + ", b, " = ", c, " = ", float(c))
    b += a
    a = 1//4
    c = b - a
    println(b, " - ", a, " = ", c, " = ", float(c))
    a -= c
    c = a * b
    println(a, " * ", b, " = ", c)
    c /= b
    a = b / c
    println(b, " / ", c, " = ", a)

    c = 0//1
    println("a=", a, "   isZero: ", a==0, "   isPositive: ", a>0, "    isNegative: ", a<0)
    println("b=", b, "   isZero: ", b==0, "   isPositive: ", b>0, "    isNegative: ", b<0)
    println("c=", c, "   isZero: ", c==0, "   isPositive: ", c>0, "    isNegative: ", c<0)

    println(a, " + 2 = ", a + 2)
    println("5 + ", a, " = ", 5 + a)
    println(a, " - 3 = ", a - 3)
    println("1 - ", a, " = ", 1 - a)
    println(a, " * (-4) = ", a * (-4))
    println("2 * ", a, " = ", 2 * a)
    println(a, " / 8 = ", a / 8)
    println("10 / ", a, " = ", 10 / a)

    c = a
    println("a = ", a, "    b = ", b, "    c = ", c)

    if a==c
        println("a == c")
      else
        println("a != c")
    end

    if a==b
        println("a == b")
      else
        println("a != b")
    end

    if a!=c
        println("a != c")
      else
        println("a == c")
    end

    if a!=b
        println("a != b")
      else
        println("a == b")
    end

    println(a, " == 3: ", a == 3)
    println("2 == ", b, ": ", 2 == b)
    println(a, " != 1: ", a != 1)
    println("-2 != ", b, ": ", -2 != b)
    println(a, " > 1: ", a > 1)
    println("3 > ", b, ": ", 3 > b)
    println(a, " >= 2: ", a >= 2)
    println("5 >= ", b, ": ", 5 >= b)
    println(a, " < 6: ", a < 6)
    println("-2 < ", b, ": ", -2 < b)
    println(a, " <= 4: ", a <= 4)
    println("3 <= ", b, ": ", 3 <= b)

    print(a, " + 3 = ")
    a += 3
    println(a)

    print(a, " - 1 = ")
    a -= 1
    println(a)

    print(a, " * 5 = ")
    a *= 5
    println(a)

    print(a, " / 2 = ")
    a /= 2
    println(a)
end

rationalTest()

max_stars_repo_name: jkovacic/math

max_issues_repo_name: jkovacic/math

max_forks_repo_name: jkovacic/math


text:
function get_sub_ind_ranges(a::Union{
        <:Tuple{AbstractVector{<:CartesianPoint{T}}, CartesianTicksTuple{T}},
        <:Tuple{AbstractVector{<:CylindricalPoint{T}}, CylindricalTicksTuple{T}}
        }) where {T}
    # For Cylindrical Coordinates this only works for z right now
    # There, we will only use t_idx_range_ax3 in the painting (for now...)
    pts, t = a[1], a[2]
    ax1, ax2, ax3 = getindex.(pts, 1), getindex.(pts, 2), getindex.(pts, 3)
    # ± 2 because due to merging close ticks of an axis of a grid,
    # the grid points which should "belong" to corners of a primitive
    # might be a bit outside of the intervals defined by the points
    # of the `extreme_points`-method.
    # ± 1 might also be enough -> testing...
    ax1_min_idx = searchsortedfirst(t[1], minimum(ax1)) - 2
    ax2_min_idx = searchsortedfirst(t[2], minimum(ax2)) - 2
    ax3_min_idx = searchsortedfirst(t[3], minimum(ax3)) - 2
    ax1_max_idx = searchsortedfirst(t[1], maximum(ax1)) + 2
    ax2_max_idx = searchsortedfirst(t[2], maximum(ax2)) + 2
    ax3_max_idx = searchsortedfirst(t[3], maximum(ax3)) + 2
    ls = (length(t[1]), length(t[2]), length(t[3]))
    if ax1_min_idx < 1 ax1_min_idx = 1 end
    if ax2_min_idx < 1 ax2_min_idx = 1 end
    if ax3_min_idx < 1 ax3_min_idx = 1 end
    if ax1_max_idx > ls[1] ax1_max_idx = ls[1] end
    if ax2_max_idx > ls[2] ax2_max_idx = ls[2] end
    if ax3_max_idx > ls[3] ax3_max_idx = ls[3] end
    t_idx_range_ax1 = ax1_min_idx:ax1_max_idx
    t_idx_range_ax2 = ax2_min_idx:ax2_max_idx
    t_idx_range_ax3 = ax3_min_idx:ax3_max_idx
    t_idx_range_ax1, t_idx_range_ax2, t_idx_range_ax3
end

get_sub_ind_ranges(p::ConstructiveSolidGeometry.AbstractSurfacePrimitive{T}, grid::CartesianGrid3D{T}) where {N,T} =
    get_sub_ind_ranges((ConstructiveSolidGeometry.extreme_points(p), TicksTuple(grid)))
get_sub_ind_ranges(p::ConstructiveSolidGeometry.AbstractSurfacePrimitive{T}, grid::CylindricalGrid{T}) where {N,T} =
    get_sub_ind_ranges((CylindricalPoint.(ConstructiveSolidGeometry.extreme_points(p)), TicksTuple(grid)))

function paint!(point_types, potential, face::AbstractSurfacePrimitive{T}, geometry, pot_value, grid::CartesianGrid3D{T}) where {T}
    t_idx_r1, t_idx_r2, t_idx_r3 = get_sub_ind_ranges(face, grid)
    ticks = (grid.axes[1].ticks, grid.axes[2].ticks, grid.axes[3].ticks)
    eX = CartesianVector{T}(1,0,0);
    eY = CartesianVector{T}(0,1,0);
    eZ = CartesianVector{T}(0,0,1);
    widths_ax1 = diff(get_extended_ticks(grid[1]))
    widths_ax2 = diff(get_extended_ticks(grid[2]))
    widths_ax3 = diff(get_extended_ticks(grid[3]))
    Δw_max_factor = T(1e-5)
    #=
        Δw_max_factor is chosen by trying out different values for it.
        This value seems to be okay if the grid is not to unevenly spaced.
        But this can be secured via the `max_ratio`- and the `the max_tick_distance`-keywords.
        The critical parameter is `csgtol` inside the `in`-method, which is currently defined through
        the widths of the voxel of the grid point next to the calculated intersection point.
        It would be probably better to turn this into an `NTuple{3,T}` and pass the widths of the voxel instead
        of one single value (`csgtol`). However, inside the `in` methods this would have to be transformed
        due to rotations and maybe the shape of the primitive. This might not be straight forward.
        But this can be investigated and certainly improved in the future.
    =#
    for i2 in t_idx_r2
        for i1 in t_idx_r1
            l = ConstructiveSolidGeometry.Line(CartesianPoint{T}(ticks[1][i1], ticks[2][i2], zero(T)), eZ)
            pts = ConstructiveSolidGeometry.intersection(face, l)
            for pt in pts
                i3 = searchsortednearest(ticks[3], pt[3])
                csgtol = Δw_max_factor * max(
                    widths_ax1[i1], widths_ax1[i1+1],
                    widths_ax2[i2], widths_ax2[i2+1],
                )
                if in(pt, grid) && in(pt, geometry, csgtol)
                    point_types[i1, i2, i3] = zero(PointType)
                    potential[i1, i2, i3] = pot_value
                end
            end
        end
    end
    for i3 in t_idx_r3
        for i1 in t_idx_r1
            l = ConstructiveSolidGeometry.Line(CartesianPoint{T}(ticks[1][i1], zero(T), ticks[3][i3]), eY)
            pts = ConstructiveSolidGeometry.intersection(face, l)
            for pt in pts
                i2 = searchsortednearest(ticks[2], pt[2])
                csgtol = Δw_max_factor * max(
                    widths_ax1[i1], widths_ax1[i1+1],
                    widths_ax3[i3], widths_ax3[i3+1]
                )
                if in(pt, grid) && in(pt, geometry, csgtol)
                    point_types[i1, i2, i3] = zero(PointType)
                    potential[i1, i2, i3] = pot_value
                end
            end
        end
    end
    for i3 in t_idx_r3
        for i2 in t_idx_r2
            l = ConstructiveSolidGeometry.Line(CartesianPoint{T}(zero(T), ticks[2][i2], ticks[3][i3]), eX)
            pts = ConstructiveSolidGeometry.intersection(face, l)
            for pt in pts
                i1 = searchsortednearest(ticks[1], pt[1])
                csgtol = Δw_max_factor * max(
                    widths_ax2[i2], widths_ax2[i2+1],
                    widths_ax3[i3], widths_ax3[i3+1]
                )
                if in(pt, grid) && in(pt, geometry, csgtol)
                    point_types[i1, i2, i3] = zero(PointType)
                    potential[i1, i2, i3] = pot_value
                end
            end
        end
    end
    nothing
end

function paint!(point_types, potential, face::AbstractSurfacePrimitive{T}, geometry, pot_value, grid::CylindricalGrid) where {T}
    t_idx_r1, t_idx_r2, t_idx_r3 = get_sub_ind_ranges(face, grid)
    ticks = (grid.axes[1].ticks, grid.axes[2].ticks, grid.axes[3].ticks)
    eZ = CartesianVector{T}(0,0,1);

    widths_ax1 = diff(get_extended_ticks(grid[1]))
    widths_ax2 = diff(get_extended_ticks(grid[2]))
    widths_ax3 = diff(get_extended_ticks(grid[3]))
    Δw_max_factor = T(1e-5)
    #=
        Δw_max_factor is chosen by trying out different values for it.
        This value seems to be okay if the grid is not to unevenly spaced.
        But this can be secured via the `max_ratio`- and the `the max_tick_distance`-keywords.
        The critical parameter is `csgtol` inside the `in`-method, which is currently defined through
        the widths of the voxel of the grid point next to the calculated intersection point.
        It would be probably better to turn this into an `NTuple{3,T}` and pass the widths of the voxel instead
        of one single value (`csgtol`). However, inside the `in` methods this would have to be transformed
        due to rotations and maybe the shape of the primitive. This might not be straight forward.
        But this can be investigated and certainly improved in the future.
    =#
    for i2 in eachindex(ticks[2])
        for i1 in eachindex(ticks[1])
            l = ConstructiveSolidGeometry.Line(CartesianPoint(CylindricalPoint{T}(ticks[1][i1], ticks[2][i2], zero(T))), eZ)
            pts = ConstructiveSolidGeometry.intersection(face, l)
            for pt in pts
                i3 = searchsortednearest(ticks[3], pt[3])
                csgtol = Δw_max_factor * max(
                    widths_ax3[i3], widths_ax3[i3+1],
                )
                if in(pt, grid) && in(pt, geometry, csgtol)
                    point_types[i1, i2, i3] = zero(PointType)
                    potential[i1, i2, i3] = pot_value
                end
            end
        end
    end
    # For this we need a the function `intersection_with_φ_axis`...
    # This will be fun... I skip this for now.

    # for i3 in t_idx_r3 # z
    #     for i1 in t_idx_r1 # r
    #         l = ConstructiveSolidGeometry.Line(CartesianPoint{T}(ticks[1][i1], zero(T), ticks[3][i3]), eY)
    #         pt = ConstructiveSolidGeometry.intersection(plane, l)
    #         i2 = searchsortednearest(ticks[2], pt[2])
    #         if in(pt, geometry, abs(ticks[2][i2] - pt[2]))
    #             point_types[i1, i2, i3] = zero(PointType)
    #             potential[i1, i2, i3] = pot_value
    #         end
    #     end
    # end

    for i3 in t_idx_r3 # z;   Maybe switch loops so that the direction of `l` has to be calculated less times..
        o = CartesianPoint{T}(zero(T), zero(T), ticks[3][i3])
        for i2 in eachindex(ticks[2]) # φ
            dir = CartesianVector(CartesianPoint(CylindricalPoint{T}(one(T), ticks[2][i2], zero(T))))
            l = ConstructiveSolidGeometry.Line(o, dir) # dir should be normalized
            pts_car = ConstructiveSolidGeometry.intersection(face, l)
            for pt_car in pts_car
                pt_cyl = CylindricalPoint(pt_car)
                # Intersection must be in positive r direction of dir: pt_cyl[2] == ticks[2][i2]
                # Use `abs(pt_cyl[2] - ticks[2][i2]) < 0.1` to avoid rounding issues
                # if it differs, it would always differ by π = 3.141... -> 0.1 is fine
                i1 = searchsortednearest(ticks[1], pt_cyl[1])
                csgtol = Δw_max_factor * max(
                    widths_ax1[i1], widths_ax1[i1+1],
                    widths_ax2[i2], widths_ax2[i2+1],
                )
                if in(pt_cyl, grid) && abs(pt_cyl[2] - ticks[2][i2]) < T(0.1) && in(pt_car, geometry, csgtol)
                    point_types[i1, i2, i3] = zero(PointType)
                    potential[i1, i2, i3] = pot_value
                end
            end
        end
    end
    nothing
end


text:
quantile(daten["Volume"].values, 0.9)
mean(daten["Volume"].values)

tpers = (daten["Volume"].values/60) / (quantile(daten["Volume"].values, 0.99)/60)
tvsmax = (daten["Volume"].values/60) / (maximum(daten["Volume"].values)/60)

plot(1:length(tpers), tpers)
plot(1:length(tpers), tvsmax)

histogram(tpers)
histogram(daten["Volume"].values)

# To try/ to do
# find a way to first sort the volume and then get the array index of the top 20/ 4 volumes
# try to rename the x axis with strings of the timestamps,
# try and redude the background grid so that you can navigate through the graph

# working like this perhaps avoids writing it to a CSV file.

symbol = "GBL"
exch = "DTB"
expiry = "201809"
currency = "EUR"
duration = "14 D"
barSize = "30 mins"

getConnectIB()

#can also do this, perhaps right a funciton with default end TimeArray, which then is the real getIB

endTime = "20180724 18:00:00"

datenSatz = getIB2(symbol, "20180725 18:00:00", exch, expiry, currency, barSize, duration)

using Plots
pyplot()

plot(datenSatz[:TimeStamp], datenSatz[:Close])
plot!(datenSatz[:TimeStamp][5:10], datenSatz[:Open][5:10])

names(datenSatz)

quantile(daten2["Volume"].values, 0.9)
mean(daten2["Volume"].values)

tpers = (daten2["Volume"].values[(end-299:end)]) / (quantile(daten2["Volume"].values, 0.95))
tvsmax = (daten2["Volume"].values) / (maximum(daten2["Volume"].values))

plot(1:length(tpers), tpers)
plot(1:length(tpers), tvsmax)

histogram(tpers)
histogram(daten2["Volume"].values)

daten2["Volume"].values[end]
daten2.timestamp[end]

#convert to number

trialTime = DateTime(2018, 6, 22, 15, 57, 20)

#length of regular time interval
#convert(Float64, (DateTime(now()) - daten2.timestamp[end]))

mscPerBar = 5 * 60 * 1000
timeToGoInBar = (5 * 60 * 1000) - convert(Float64, (daten2.timestamp[end] - trialTime)) # potentially either reduce "end" or different difference

#extrapolation factor
extrPolFactor = 1 - timeToGoInBar / mscPerBar

#now extrapolate current volume as if it would grow in the same rate,
#if this is then bigger then the quantile, it should raise a red flag

(daten2["Volume"].values[end] / extrPolFactor) / (quantile(daten2["Volume"].values, 0.95))

writetimearray(daten2, "bund.csv")

neu = readtimearray("bund.csv")

symbol = "GBL"
exch = "DTB"
expiry = "201809"
currency = "EUR"
duration = "1 M"
barSize = "1 day"

data2 = getDailyIB(symbol, exch, expiry, currency, barSize, duration)

Dates.firstdayofweek(Date(now()-Dates.Day(5))) + Dates.Day(1)
Dates.firstdayofweek(Date(now()-Dates.Day(5)))
Dates.firstdayofweek(Date(now()))


text:
# functions designed to be used outside the optimization

# get functions
"""
    get_x0(optparam::OptimizationParameter)

Returns initial value(s) for `optparam`
"""
get_x0(optparam::OptimizationParameter) = optparam.x0

"""
    get_lb(optparam::OptimizationParameter)

Returns lower bound(s) for `optparam`
"""
get_lb(optparam::OptimizationParameter) = optparam.lb

"""
    get_ub(optparam::OptimizationParameter)

Returns upper bound(s) for `optparam`
"""
get_ub(optparam::OptimizationParameter) = optparam.ub

"""
    get_scaling(optparam::OptimizationParameter)

Returns scaling parameter(s) for `optparam`
"""
get_scaling(optparam::OptimizationParameter) = optparam.scaling

"""
    get_dv(optparam::OptimizationParameter)

Returns design variable flag(s) for `optparam`
"""
get_dv(optparam::OptimizationParameter) = optparam.dv

"""
    get_description(optparam::OptimizationParameter)

Returns the description for `optparam`
"""
get_description(optparam::OptimizationParameter) = optparam.description

"""
    get_x0(optparams, field)

Returns initial value(s) for `field` in `optparams`
"""
get_x0(optparams, field) = optparams[field].x0

"""
    get_lb(optparams, field)

Returns lower bound(s) for `field` in `optparams`
"""
get_lb(optparams, field) = optparams[field].lb

"""
    get_ub(optparams, field)

Returns upper bound(s) for `field` in `optparams`
"""
get_ub(optparams, field) = optparams[field].ub

"""
    get_scaling(optparams, field)

Returns scaling parameter(s) for `field` in `optparams`
"""
get_scaling(optparams, field) = optparams[field].scaling

"""
    get_dv(optparams, field)

Returns design variable flag(s) for `field` in `optparams`
"""
get_dv(optparams, field) = optparams[field].dv

"""
    get_description(optparams, field)

Returns the description for `field` in `optparams`
"""
get_description(optparams, field) = optparams[field].description

# set functions
"""
    set_x0(optparam::OptimizationParameter, x0)

Returns `optparam` with the initial value(s) given by `x0`
"""
set_x0(optparam::OptimizationParameter, x0)= OptimizationParameter(
    x0, lb=optparam.lb, ub=optparam.ub, scaling=optparam.scaling, dv=optparam.dv, description=optparam.description)

"""
    set_lb(optparam::OptimizationParameter, lb)

Returns `optparam` with the lower bound(s) given by `lb`
"""
set_lb(optparam::OptimizationParameter, lb) = OptimizationParameter(
    optparam.x0, lb=lb, ub=optparam.ub, scaling=optparam.scaling, dv=optparam.dv, description=optparam.description)

"""
    set_ub(optparam::OptimizationParameter, ub)

Returns `optparam` with the upper bound(s) given by `ub`
"""
set_ub(optparam::OptimizationParameter, ub) = OptimizationParameter(
    optparam.x0, lb=optparam.lb, ub=ub, scaling=optparam.scaling, dv=optparam.dv, description=optparam.description)

"""
    set_scaling(optparam::OptimizationParameter, scaling)

Returns `optparam` with the scaling parameter given by `scaling`
"""
set_scaling(optparam::OptimizationParameter, scaling) = OptimizationParameter(
    optparam.x0, lb=optparam.lb, ub=optparam.ub, scaling=scaling, dv=optparam.dv, description=optparam.description)

"""
    set_dv(optparam::OptimizationParameter, dv)

Returns `optparam` with the design variable flag(s) given by `dv`
"""
set_dv(optparam::OptimizationParameter, dv) = OptimizationParameter(
    optparam.x0, lb=optparam.lb, ub=optparam.ub, scaling=optparam.scaling, dv=dv, description=optparam.description)

"""
    set_description(optparam::OptimizationParameter, description)

Returns `optparam` with the description given by `description`
"""
set_description(optparam::OptimizationParameter, description) = OptimizationParameter(
    optparam.x0, lb=optparam.lb, ub=optparam.ub, scaling=optparam.scaling, dv=optparam.dv, description=description)

"""
    set_x0!(optparams::AbstractDict, field, x0)

Modifies the `field` in `optparams` to have the initial value(s) given by `x0`
"""
function set_x0!(optparams::AbstractDict, field, x0)
    optparams[field] = set_x0(optparams[field], x0)
    return optparams
end

"""
    set_lb!(optparams::AbstractDict, field, lb)

Modifies the `field` in `optparams` to have the lower bound(s) given by `lb`
"""
function set_lb!(optparams::AbstractDict, field, lb)
    optparams[field] = set_lb(optparams[field], lb)
    return optparams
end

"""
    set_ub!(optparams::AbstractDict, field, ub)

Modifies the `field` in `optparams` to have the upper bound(s) given by `ub`
"""
function set_ub!(optparams::AbstractDict, field, ub)
    optparams[field] = set_ub(optparams[field], ub)
    return optparams
end

"""
    set_scaling!(optparams::AbstractDict, field, scaling)

Modifies the `field` in `optparams` to have the scaling parameter(s)
given by `scaling`
"""
function set_scaling!(optparams::AbstractDict, field, scaling)
    optparams[field] = set_scaling(optparams[field], scaling)
    return optparams
end

"""
    set_dv!(optparams::AbstractDict, field, dv)

Modifies the `field` in `optparams` to have the design variable flag(s)
given by `dv`
"""
function set_dv!(optparams::AbstractDict, field, dv)
    optparams[field] = set_dv(optparams[field], dv)
    return optparams
end

"""
    set_description!(optparams::AbstractDict, field, description)

Modifies the `field` in `optparams` to have the description given by `description`
"""
function set_description!(optparams::AbstractDict, field, description)
    optparams[field] = set_description(optparams[field], description)
    return optparams
end

"""
    set_x0(optparams::AbstractDict, field, x0)

Modifies the `field` in `optparams` to have the initial value(s) given by `x0`
"""
set_x0(optparams::AbstractDict, field, x0) = set_x0!(copy(optparams), field, x0)

"""
    set_lb(optparams::AbstractDict, field, lb)

Modifies the `field` in `optparams` to have the lower bound(s) given by `lb`
"""
set_lb(optparams::AbstractDict, field, lb) = set_lb!(copy(optparams), field, lb)

"""
    set_ub(optparams::AbstractDict, field, ub)

Modifies the `field` in `optparams` to have the upper bound(s) given by `ub`
"""
set_ub(optparams::AbstractDict, field, ub) = set_ub!(copy(optparams), field, ub)

"""
    set_scaling(optparams::AbstractDict, field, scaling)

Modifies the `field` in `optparams` to have the scaling parameter(s)
given by `scaling`
"""
set_scaling(optparams::AbstractDict, field, scaling) = set_scaling!(copy(optparams), field, scaling)

"""
    set_dv(optparams::AbstractDict, field, dv)

Modifies the `field` in `optparams` to have the design variable flag(s)
given by `dv`
"""
set_dv(optparams::AbstractDict, field, dv) = set_dv!(copy(optparams), field, dv)

"""
    set_description(optparams::AbstractDict, field, description)

Modifies the `field` in `optparams` to have the description given by `description`
"""
set_description(optparams::AbstractDict, field, description) = set_description!(copy(optparams), field, description)

"""
    set_x0(optparams::NamedTuple, field, x0)

Modifies the `field` in `optparams` to have the initial value(s) given by `x0`
"""
set_x0(optparams::NamedTuple{names}, field, x0) where names =
    NamedTuple{names}(Tuple(name == field ? set_x0(optparam, x0) : optparam
    for (name, optparam) in pairs(optparams)))

"""
    set_lb(optparams::NamedTuple, field, lb)

Modifies the `field` in `optparams` to have the lower bound(s) given by `lb`
"""
set_lb(optparams::NamedTuple{names}, field, lb) where names =
    NamedTuple{names}(Tuple(name == field ? set_lb(optparam, lb) :
    optparam for (name, optparam) in pairs(optparams)))

"""
    set_ub(optparams::NamedTuple, field, ub)

Modifies the `field` in `optparams` to have the upper bound(s) given by `ub`
"""
set_ub(optparams::NamedTuple{names}, field, ub) where names =
    NamedTuple{names}(Tuple(name == field ? set_ub(optparam, ub) :
    optparam for (name, optparam) in pairs(optparams)))

"""
    set_scaling(optparams::NamedTuple, field, scaling)

Modifies the `field` in `optparams` to have the scaling parameter(s)
given by `scaling`
"""
set_scaling(optparams::NamedTuple{names}, field, scaling) where names =
    NamedTuple{names}(Tuple(name == field ? set_scaling(optparam, scaling) :
    optparam for (name, optparam) in pairs(optparams)))

"""
    set_dv(optparams::NamedTuple, field, dv)

Modifies the `field` in `optparams` to have the design variable flag(s)
given by `dv`
"""
set_dv(optparams::NamedTuple{names}, field, dv) where names =
    NamedTuple{names}(Tuple(name == field ? set_dv(optparam, dv) :
    optparam for (name, optparam) in pairs(optparams)))

"""
    set_description(optparams::NamedTuple, field, description)

Modifies the `field` in `optparams` to have the description given by `description`
"""
set_description(optparams::NamedTuple{names}, field, description) where names =
    NamedTuple{names}(Tuple(name == field ? set_description(optparam, description) :
    optparam for (name, optparam) in pairs(optparams)))

"""
    update_parameters(optparams, params)

Inserts parameters from `params` into `optparams` and returns the result.
`optparams` is a named tuple or dictionary filled with objects of type OptimizationParameter.
`params` is a named tuple or dictionary with the parameter values to be replaced or a vector of scaled design variables.
"""
function update_parameters(optparams::NamedTuple, params)
    optparams = named_tuple_to_dict(optparams)
    params = named_tuple_to_dict(params)
    new = update_parameters!(optparams, params)
    return dict_to_named_tuple(new)
end

function update_parameters(optparams::NamedTuple, x::AbstractVector)
    optparams = named_tuple_to_dict(optparams)
    new = update_parameters!(optparams, x)
    return dict_to_named_tuple(new)
end

function update_parameters(optparams::AbstractDict, params)
    params = named_tuple_to_dict(params)
    new = update_parameters!(copy(optparams), params)
    return new
end

function update_parameters(optparams::AbstractDict, x::AbstractVector)
    new = update_parameters!(copy(optparams), x)
    return new
end

"""
    update_parameters!(optparams, params)

Inserts parameters from `params` into `optparams` and returns the result.
`optparams` is a dictionary filled with objects of type OptimizationParameter
`params` is a dictionary with the parameter values to be replaced or a vector of scaled design variables
"""
function update_parameters!(optparams, params)
    for name in keys(params)
        optparams[name] = set_x0(optparams[name], params[name])
    end
    return optparams
end

function update_parameters!(optparams, x::AbstractVector)
    params = get_values(optparams, x)
    for name in keys(params)
        optparams[name] = set_x0(optparams[name], params[name])
    end
    return optparams
end

"""
    update_design_variables(original, optimized)

Copies design variables that are active in both `original` and `optimized` from
`optimized` into `original` and returns the result.
`original` is a named tuple or dictionary filled with objects of type OptimizationParameter
`optimized` is a named tuple or dictionary filled with objects of type OptimizationParameter
"""
function update_design_variables(original::NamedTuple, optimized)
    original = named_tuple_to_dict(original)
    optimized = named_tuple_to_dict(optimized)
    new = update_design_variables!(original, optimized)
    return dict_to_named_tuple(new)

end

function update_design_variables(original::AbstractDict, optimized)
    optimized = named_tuple_to_dict(optimized)
    new = update_design_variables!(original, optimized)
    return new
end

"""
    update_design_variables!(original, optimized)

Copies design variables that are active in both `original` and `optimized` from
`optimized` into `original` and returns the result.
`original` is a named tuple or dictionary filled with objects of type OptimizationParameter
`optimized` is a named tuple or dictionary filled with objects of type OptimizationParameter
"""
function update_design_variables!(original, optimized)

    # loop through all fields in original optimization parameters
    for name in keys(original)

        # check if the optimized solution has this field
        if name in keys(optimized)

            # extract old and new parameters
            old = original[name]
            opt = optimized[name]

            # check if the old and optimized parameters have the same length
            if length(old.x0) == length(opt.x0)

                # update design variables that are active in both sets of parameters
                if typeof(old.x0) <: AbstractArray
                    T = promote_type(eltype(old.x0), eltype(opt.x0))
                    x0 = T.(old.x0)
                    idv = findall(old.dv .& opt.dv)
                    for i in idv
                        x0[i] = opt.x0[i]
                    end
                else
                    x0 = opt.x0
                end

                # overwrite old parameter
                original[name] = set_x0(old, x0)
            end
        end
    end

    return original
end

"""
    assemble_input(optparams)

Assembles the input for the optimization problem: x0, lb, ub
"""
function assemble_input(optparams)

    # get indices and number of design variables
    indices, ndv = design_variable_indices(optparams)

    # initialize output arrays
    x0 = Array{Float64,1}(undef, ndv)
    lb = Array{Float64,1}(undef, ndv)
    ub = Array{Float64,1}(undef, ndv)

    # populate output arrays
    for (name, optparam) in pairs(optparams)
        if any(optparam.dv)
            c_lb = optparam.lb
            c_ub = optparam.ub
            scaling = optparam.scaling
            idv = findall(optparam.dv)
            for j = 1:length(indices[name])
                initial = optparam.x0[idv[j]]
                lower = optparam.lb[idv[j]]
                upper = optparam.ub[idv[j]]
                scaling = optparam.scaling[idv[j]]
                offset = design_variable_offset(lower, upper)
                lb[indices[name][j]] = (lower + offset) * scaling
                x0[indices[name][j]] = (initial + offset) * scaling
                ub[indices[name][j]] = (upper + offset) * scaling
            end
        end
    end

    return x0, lb, ub
end

"""
    named_tuple_to_dict(nt)

Converts a named tuple to an ordered dictionary
"""
named_tuple_to_dict(nt) = OrderedDict(pairs(nt))

"""
    dict_to_named_tuple(dict)

Converts a dictionary to a named tuple
"""
dict_to_named_tuple(dict) = (; dict...)


text:
module GaussQuadrature

using SpecialFunctions: gamma, lgamma

# October 2013 by Bill McLean, School of Maths and Stats,
# The University of New South Wales.
#
# Based on earlier Fortran codes
#
# gaussq.f original version 20 Jan 1975 from Stanford
# gaussq.f modified 21 Dec by Eric Grosse
# gaussquad.f95 Nov 2005 by Bill Mclean
#
# This module provides functions to compute the abscissae x[j] and
# weights w[j] for the classical Gauss quadrature rules, including
# the Radau and Lobatto variants.  Thus, the sum
#
#           n
#           ∑  w[j] f(x[j])
#          j=1
#
# approximates
#
#           hi
#          ∫  f(x) w(x) dx
#          lo
#
# where the weight function w(x) and interval lo < x < hi are as shown
# in the table below.
#
# Name                      Interval     Weight Function
#
# Legendre                 -1 < x < 1          1
# Chebyshev (first kind)   -1 < x < 1     1 / sqrt(1-x²)
# Chebyshev (second kind)  -1 < x < 1       sqrt(1-x²)
# Jacobi                   -1 < x < 1     (1-x)ᵅ (1+x)ᵝ
# Laguerre                  0 < x < ∞     xᵅ exp(-x)
# Hermite                  -∞ < x < ∞      exp(-x²)
#
# In addition to these classical rules, the module generates Gauss rules
# for logarithmic weights of the form
#
#    w(x) = x^ρ log(1/x)   for 0 < x < 1.
#
# For the Jacobi and Laguerre rules we require α > -1 and
# β > -1, so that the weight function is integrable.  Likewise, for
# log weight we require ρ > -1.
#
# Use the endpt argument to include one or both of the end points
# of the interval of integration as an abscissa in the quadrature
# rule, as follows.
#
# endpt = neither   Default      lo < x[j] < hi, j = 1:n.
# endpt = left      Left Radau   lo = x[1] < x[j] < hi, j = 2:n.
# endpt = right     Right Radau  lo < x[j] < x[n] = hi, j = 1:n-1.
# endpt = both      Lobatto      lo = x[1] < x[j] < x[n] = hi, j = 2:n-1.
#
# These labels make up an enumeration of type EndPt.
#
# The code uses the Golub and Welsch algorithm, in which the abscissae
# x[j] are the eigenvalues of a symmetric tridiagonal matrix whose
# entries depend on the coefficients in the 3-term recurrence relation
# for the othonormal polynomials generated by the weighted inner product.
#
# References:
#
#   1.  Golub, G. H., and Welsch, J. H., Calculation of Gaussian
#       quadrature rules, Mathematics of Computation 23 (April,
#       1969), pp. 221-230.
#   2.  Golub, G. H., Some modified matrix eigenvalue problems,
#       Siam Review 15 (april, 1973), pp. 318-334 (section 7).
#   3.  Stroud and Secrest, Gaussian Quadrature Formulas, Prentice-
#       Hall, Englewood Cliffs, N.J., 1966.

export neither, left, right, both
export legendre, legendre_coefs
export chebyshev, chebyshev_coefs
export jacobi, jacobi_coefs
export laguerre, laguerre_coefs
export hermite, hermite_coefs
export shifted_legendre_coefs
export logweight, logweight_coefs
export modified_moments, modified_chebyshev
export custom_gauss_rule, orthonormal_poly
export special_eigenproblem!

# Enumeration type used to specify which endpoints of the integration
# interval should be included amongst the quadrature points: neither,
# left, right or both.
struct EndPt
    label :: Char
end

const neither = EndPt('N')
const left    = EndPt('L')
const right   = EndPt('R')
const both    = EndPt('B')

#include("to_be_removed.jl")

# Maximum number of QL iterations used by steig!.
# You might need to increase this.
maxiterations = Dict(Float32 => 30, Float64 => 30, BigFloat => 40)

"""
    x, w = legendre(T, n, endpt=neither)

Returns points `x` and weights `w` for the `n`-point Gauss-Legendre rule
for the interval `-1 < x < 1` with weight function `w(x) = 1`.

Use `endpt=left`, `right` or `both` for the left Radau, right Radau or
Lobatto rules, respectively.
"""
function legendre(::Type{T}, n::Integer,
                  endpt::EndPt=neither) where {T<:AbstractFloat}
    @assert n ≥ 1
    a, b = legendre_coefs(T, n)
    return custom_gauss_rule(-one(T), one(T), a, b, endpt)
end

"""
Convenience function with type `T = Float64`:

    x, w = legendre(n, endpt=neither)
"""
legendre(n, endpt=neither) = legendre(Float64, n, endpt)

function legendre_coefs(::Type{T}, n::Integer) where {T<:AbstractFloat}
    a = zeros(T, n)
    b = zeros(T, n+1)
    b[1] = sqrt(convert(T, 2))
    for k = 2:n+1
        b[k] = (k-1) / sqrt(convert(T, (2k-1)*(2k-3)))
    end
    return a, b
end

"""
    x, w = chebyshev(T, n, kind=1, endpt=neither)

Returns points `x` and weights `w` for the `n`-point Gauss-Chebyshev rule
for the interval `-1 < x < 1` with weight function

    w(x) = 1 / sqrt(1-x²)   if kind=1
    w(x) = sqrt(1-x²)       if kind=2.

Use `endpt=left`, `right` or `both` for the left Radau, right Radau or
Lobatto rules, respectively.
"""
function chebyshev(::Type{T}, n::Integer, kind::Integer=1,
                   endpt::EndPt=neither) where {T<:AbstractFloat}
    @assert n ≥ 1
    a, b = chebyshev_coefs(T, n, kind)
    return custom_gauss_rule(-one(T), one(T), a, b, endpt)
end

"""
Convenience function with type `T = Float64`:

    x, w = chebyshev(n, kind=1, endpt=neither)

"""
chebyshev(n, kind=1, endpt=neither) = chebyshev(Float64, n, kind, endpt)

function chebyshev_coefs(::Type{T}, n::Integer,
                         kind::Integer) where {T<:AbstractFloat}
    half = convert(T, 1//2)
    a = zeros(T, n)
    b = fill(half, n+1)
    if kind == 1
        b[1] = sqrt(convert(T, pi))
	if n >= 2
            b[2] = sqrt(half)
        end
    elseif kind == 2
        b[1] = sqrt(half * convert(T, pi))
    else
        error("Unsupported value for kind")
    end
    return a, b
end

"""
    x, w = jacobi(n, α, β, endpt=neither)

Returns points `x` and weights `w` for the `n`-point Gauss-Jacobi rule
for the interval `-1 < x < 1` with weight function

    w(x) = (1-x)ᵅ (1+x)ᵝ.

Use `endpt=left`, `right` or `both` for the left Radau, right Radau or
Lobatto rules, respectively.
"""
function jacobi(n::Integer, α::T, β::T,
                endpt::EndPt=neither) where {T<:AbstractFloat}
    @assert n ≥ 1
    @assert α > -1.0 && β > -1.0
    a, b = jacobi_coefs(n, α, β)
    custom_gauss_rule(-one(T), one(T), a, b, endpt)
end

function jacobi_coefs(n::Integer, α::T, β::T) where {T<:AbstractFloat}
    a = zeros(T, n)
    b = zeros(T, n+1)
    ab = α + β
    abi = ab + 2
    b[1] = 2^((ab+1)/2) * exp(
             ( lgamma(α+1) + lgamma(β+1) - lgamma(abi) )/2 )
    a[1] = ( β - α ) / abi
    b[2] = sqrt( 4*(α+1)*(β+1) / ( (ab+3)*abi*abi ) )
    a2b2 = β*β - α*α
    for i = 2:n
        abi = ab + 2i
        a[i]   = a2b2 / ( (abi-2)*abi )
        b[i+1] = sqrt( 4i*(α+i)*(β+i)*(ab+i) /
                     ( (abi*abi-1)*abi*abi ) )
    end
    return a, b
end

"""
    x, w = laguerre(n, α, endpt=neither)

Returns points `x` and weights `w` for the `n`-point Gauss-Laguerre rule
for the interval `0 < x < ∞` with weight function

    w(x) = xᵅ exp(-x),   α > -1.

Use `endpt=left` for the left Radau rule.
"""
function laguerre(n::Integer, α::T,
                  endpt::EndPt=neither) where {T<:AbstractFloat}
    @assert n ≥ 1
    @assert α > -1.0
    @assert endpt in [neither, left]
    a, b = laguerre_coefs(n, α)
    custom_gauss_rule(zero(T), convert(T, Inf), a, b, endpt)
end

function laguerre_coefs(n::Integer, α::T) where {T<:AbstractFloat}
    a = zeros(T, n)
    b = zeros(T, n+1)
    b[1] = sqrt(gamma(α+1))
    for i = 1:n
        a[i] = 2i - 1 + α
        b[i+1] = sqrt( i*(α+i) )
    end
    return a, b
end

"""
    x, w = hermite(T, n)

Returns points `x` and weights `w` for the `n`-point Gauss-Laguerre rule
for the interval `-∞ < x < ∞` with weight function

    w(x) = exp(-x²).
"""
function hermite(::Type{T}, n::Integer) where {T<:AbstractFloat}
    @assert n ≥ 1
    a, b = hermite_coefs(T, n)
    custom_gauss_rule(convert(T, -Inf), convert(T, Inf), a, b, neither)
end

"""
Convenience function with type `T = Float64`:

    x, w = hermite(n)
"""
hermite(n) = hermite(Float64, n)

function hermite_coefs(::Type{T}, n::Integer) where {T<:AbstractFloat}
    a = zeros(T, n)
    b = zeros(T, n+1)
    b[1] = sqrt(sqrt(convert(T, pi)))
    for i = 1:n
        iT = convert(T, i)
        b[i+1] = sqrt(iT/2)
    end
    return a, b
end

"""
    x, w = logweight(T, n, r=0, endpt=neither)

Returns points `x` and weights `w` for the n-point Gauss rule on
the interval `0 < x < 1` with weight function

    w(x) = xʳ log(1/x),    r ≥ 0.
"""
function logweight(::Type{T}, n::Integer, r::Integer,
                   endpt::EndPt=neither) where {T<:AbstractFloat}
    @assert n ≥ 1
    @assert r ≥ 0
    α, β = logweight_coefs(T, n, r)
    custom_gauss_rule(zero(T), one(T), α, β, endpt)
end

"""
Convenience function with type `T = Float64`:

    x, w = logweight(n, r=0, endpt=neither)
"""
logweight(n, r=0, endpt=neither) = logweight(Float64, n, r, endpt)

"""
More general method works when `w(x) = x^ρ log(1/x)` for real `ρ > -1`:

    x, w = logweight(n, ρ, endpt=neither)
"""
function logweight(n::Integer, ρ::T,
                   endpt::EndPt=neither) where {T<:AbstractFloat}
    @assert n ≥ 1
    @assert ρ > -1
    α, β = logweight_coefs(n, ρ)
    custom_gauss_rule(zero(T), one(T), α, β, endpt)
end

function logweight_coefs(::Type{T}, n::Integer,
                         r::Integer) where {T<:AbstractFloat}
    a, b = shifted_legendre_coefs(T, 2n)
    ν = modified_moments(T, n, r)
    α, β, σ = modified_chebyshev(a, b, ν)
    return α, β
end

function logweight_coefs(n::Integer, ρ::T) where {T<:AbstractFloat}
    a, b = shifted_legendre_coefs(T, 2n)
    ν = modified_moments(n, ρ)
    α, β, σ = modified_chebyshev(a, b, ν)
    return α, β
end

function modified_moments(::Type{T}, n::Integer,
                          r::Integer) where {T<:AbstractFloat}
    @assert r ≥ 0
    @assert n ≥ 0
    m = min(2n, r+1)
    rrp1 = 1 / convert(T, r+1) # reciprocal r+1
    B = one(T)
    S = rrp1
    ν = zeros(T, 2n)
    ν[1] = rrp1 * S
    for l = 2:m
        j = l - 1
	rrpj = 1 / convert(T, r + 1 + j)
        rmj = convert(T, r + 1 - j)
	rrmj = 1 / rmj
        B *= rmj * rrpj
        S +=  rrpj - rrmj
        tlm1 = convert(T, 2l-1)
        ν[l] = rrp1 * B * S * sqrt(tlm1)
    end
    if 2n > r+1
        p = one(T)
        for j = 1:r
            p *= j / convert(T, 2(2j+1))
        end
        ν[r+2] = - sqrt(convert(T, 2r+3)) *p / ( 2(r+1) )
        for l = r+2:2n-1
            Tl = convert(T, l)
            ν[l+1] = - ν[l] * ( (Tl-r-1)/(Tl+r+1) ) * sqrt((2Tl+1)/(2Tl-1) )
        end
    end
    return ν
end

function modified_moments(n::Integer, ρ::T) where {T<:AbstractFloat}
    @assert ρ > -1
    @assert n ≥ 0
    r = (ρ<0) ? 0 : round(Integer, ρ)
    m = min(2n, r+1)
    S = 1 / ( 1 + ρ )
    B = one(T)
    ν = zeros(T, 2n)
    ν[1] = S / (1+ρ)
    for l = 2:m
        j = l - 1
        B *= (ρ+1-j) / (ρ+1+j)
	S += 1 / (ρ+1+j) - 1 / (ρ+1-j)
        ν[l] = B * S * sqrt(convert(T,2l-1))
    end
    if 2n > r+1
        l = r + 2
	j = l - 1
	X = ( (ρ-r) / (ρ+r+2) - 1 ) / (ρ+r+2)
        tlm1 = convert(T, 2l-1)
	ν[l] = ( B / (1+ρ) ) * ( X + ( (ρ-r)/(ρ+r+2) ) * S ) * sqrt(tlm1)
        for l = r+3:2n
	    j = l - 1
            B *= (ρ+1-j) / (ρ+1+j)
            S += 1 / (ρ+1+j) - 1 / (ρ+1-j)
        tlm1 = convert(T, 2l-1)
	    ν[l] = ( B / (1+ρ) ) * ( X + ( (ρ-r)/(ρ+r+2) ) * S ) * sqrt(tlm1)
        end
    end
    return ν
end

function shifted_legendre_coefs(T, n)
    a = zeros(T, n)
    b = zeros(T, n+1)
    fill!(a, 1/convert(T,2))
    b[1] = one(T)
    for k = 2:n+1
        b[k] = (k-1) / ( 2 * sqrt( convert(T, (2k-1)*(2k-3)) ) )
    end
    return a, b
end

"""
    x, w = custom_gauss_rule(lo, hi, a, b, endpt, maxits=maxiterations[T])

Generates the points `x` and weights `w` for a Gauss rule with weight
function `w(x)` on the interval `lo < x < hi`.

The arrays `a` and `b` hold the coefficients (as given, for instance, by
`legendre_coeff`) in the three-term recurrence relation for the monic
orthogonal polynomials `p(0,x)`, `p(1,x)`, `p(2,x)`, ... , that is,

    p(k, x) = (x-a[k]) p(k-1, x) - b[k]² p(k-2, x),    k ≥ 1,

where `p(0, x) = 1` and, by convention, `p(-1, x) = 0` with

              hi
    b[1]^2 = ∫  w(x) dx.
             lo

Thus, `p(k, x) = xᵏ + lower degree terms` and

     hi
    ∫  p(j, x) p(k, x) w(x) dx = 0 if j ≠ k.
    lo
"""
function custom_gauss_rule(lo::T, hi::T,
             a::Array{T,1}, b::Array{T,1}, endpt::EndPt,
             maxits::Integer=maxiterations[T]) where {T<:AbstractFloat}
    n = length(a)
    @assert length(b) == n+1
    if endpt == left
        if n == 1
            a[1] = lo
        else
            a[n] = solve(n, lo, a, b) * b[n]^2 + lo
        end
    elseif endpt == right
        if n == 1
            a[1] = hi
        else
            a[n] = solve(n, hi, a, b) * b[n]^2 + hi
        end
    elseif endpt == both
        if n == 1
            error("Must have at least two points for both ends.")
        end
        g = solve(n, lo, a, b)
        t1 = ( hi - lo ) / ( g - solve(n, hi, a, b) )
        b[n] = sqrt(t1)
        a[n] = lo + g * t1
    end
    w = zero(a)
    special_eigenproblem!(a, b, w, maxits)
    for i = 1:n
        w[i] = (b[1] * w[i])^2
    end
    idx = sortperm(a)
    # Ensure end point values are exact.
    if endpt in (left, both)
        a[idx[1]] = lo
    elseif endpt in (right, both)
        a[idx[n]] = hi
    end
    return a[idx], w[idx]
end

"""
    α, β, σ = modified_chebyshev(a, b, ν)

Implements the modified Chebyshev algorithm described in `doc/notes.tex`
and used in `logweight_coefs`.
"""
function modified_chebyshev(a::Vector{T}, b::Vector{T},
                            ν::Vector{T}) where {T<:AbstractFloat}
    m = length(ν)
    @assert m % 2 == 0 && m >= 2
    n = div(m, 2)
    @assert length(a) >= max(1, 2n-1)
    @assert length(b) >= max(1, 2n-1)
    α = zeros(T, n)
    β = zeros(T, n+1)
    β[1] = sqrt(ν[1])
    σ = zeros(T, 2n, n)
    for l = 1:2n
        σ[l,1] = ν[l] / β[1]
    end
    α[1] = a[1] + b[2] * ( σ[2,1] / σ[1,1] )
    if n >= 2
        # handle k=1 separately since we have omitted σ[l,0]=0.
        t = ( 1 + (b[3]/b[2]) * (σ[3,1]/σ[1,1])
                + ( (a[2]-α[1])/b[2] ) * (σ[2,1]/σ[1,1]) )
        if t < 0
            error("modified Chebyshev algorithm failed at k = 1")
        end
        β[2] = b[2] * sqrt(t)
        for l = 1:2n-2
            σ[l+1,2] = ( (b[l+2]/β[2])*σ[l+2,1]
                       + ((a[l+1]-α[1])/β[2]) * σ[l+1,1]
  	               + (b[l+1]/β[2]) * σ[l,1] )
        end
        α[2] = a[2] + b[3] * (σ[3,2]/σ[2,2]) - β[2] * (σ[2,1]/σ[2,2])
        # general k
        for k = 2:n-1
            t = ( 1 + (b[k+2]/b[k+1]) * (σ[k+2,k]/σ[k,k])
                    + ( (a[k+1]-α[k])/b[k+1] ) * (σ[k+1,k]/σ[k,k])
                    - (β[k]/b[k+1]) * (σ[k+1,k-1]/σ[k,k]) )
            β[k+1] = b[k+1] * sqrt(t)
            if t < 0
                error("modified Chebyshev algorithm failed at k = $k")
            end
            for l = k:2n-k-1
                σ[l+1,k+1] = ( (b[l+2]/β[k+1]) * σ[l+2,k]
                             + ( (a[l+1]-α[k])/β[k+1] ) * σ[l+1,k]
                             + (b[l+1]/β[k+1]) * σ[l,k]
                             - (β[k]/β[k+1]) * σ[l+1,k-1] )
            end
	    α[k+1] = ( a[k+1] + b[k+2] * (σ[k+2,k+1]/σ[k+1,k+1])
                              - β[k+1] * (σ[k+1,k]/σ[k+1,k+1]) )
        end
    end
    return α, β, σ
end

function solve(n::Integer, shift::T, a::Array{T,1},
               b::Array{T,1}) where {T<:AbstractFloat}
    #
    # Perform elimination to find the nth component s = delta[n]
    # of the solution to the nxn linear system
    #
    #     ( J_n - shift I_n ) delta = e_n,
    #
    # where J_n is the symmetric tridiagonal matrix with diagonal
    # entries a[i] and off-diagonal entries b[i], and e_n is the nth
    # standard basis vector.
    #
    t = a[1] - shift
    for i = 2:n-1
        t = a[i] - shift - b[i]^2 / t
    end
    return one(t) / t
end

function special_eigenproblem!(d::Array{T,1}, e::Array{T,1}, z::Array{T,1},
                               maxits::Integer) where {T<:AbstractFloat}
    #
    # Finds the eigenvalues and first components of the normalised
    # eigenvectors of a symmetric tridiagonal matrix by the implicit
    # QL method.
    #
    # d[i]   On entry, holds the ith diagonal entry of the matrix.
    #        On exit, holds the ith eigenvalue.
    #
    # e[i]   On entry, holds the [i,i-1] entry of the matrix for
    #        i = 2, 3, ..., n.  (The value of e[1] is not used.)
    #        On exit, e is overwritten.
    #
    # z[i]   On exit, holds the first component of the ith normalised
    #        eigenvector associated with d[i].
    #
    # maxits The maximum number of QL iterations.
    #
    # Martin and Wilkinson, Numer. Math. 12: 377-383 (1968).
    # Dubrulle, Numer. Math. 15: 450 (1970).
    # Handbook for Automatic Computation, Vol ii, Linear Algebra,
    #        pp. 241-248, 1971.
    #
    # This is a modified version of the Eispack routine imtql2.
    #
    n = length(z)
    z[1] = one(T)
    z[2:n] .= zero(T)
    e[n+1] = zero(T)

    if n == 1 # Nothing to do for a 1x1 matrix.
        return
    end
    for l = 1:n
        for j = 1:maxits
            # Look for small off-diagonal elements.
            m = n
            for i = l:n-1
                if abs(e[i+1]) <= eps(T) * ( abs(d[i]) + abs(d[i+1]) )
                    m = i
                    break
                end
            end
            p = d[l]
            if m == l
                continue
            end
            if j == maxits
                msg = string("No convergence after ", j, " iterations",
                             " (try increasing maxits)")
                error(msg)
            end
            # Form shift
            g = ( d[l+1] - p ) / ( 2 * e[l+1] )
            r = hypot(g, one(T))
            g = d[m] - p + e[l+1] / ( g + copysign(r, g) )
            s = one(T)
            c = one(T)
            p = zero(T)
            for i = m-1:-1:l
                f = s * e[i+1]
                b = c * e[i+1]
                if abs(f) <  abs(g)
                    s = f / g
                    r = hypot(s, one(T))
                    e[i+2] = g * r
                    c = one(T) / r
                    s *= c
                else
                    c = g / f
                    r = hypot(c, one(T))
                    e[i+2] = f * r
                    s = one(T) / r
                    c *= s
                end
                g = d[i+1] - p
                r = ( d[i] - g ) * s + 2 * c * b
                p = s * r
                d[i+1] = g + p
                g = c * r - b
                # Form first component of vector.
                f = z[i+1]
                z[i+1] = s * z[i] + c * f
                z[i]   = c * z[i] - s * f
            end # loop over i
            d[l] -= p
            e[l+1] = g
            e[m+1] = zero(T)
        end # loop over j
    end # loop over l
end

function orthonormal_poly(x::AbstractVector{T}, a::AbstractVector{T},
                          b::AbstractVector{T}) where {T<:AbstractFloat}
    # p[i,j] = value at x[i] of orthonormal polynomial of degree j-1.
    m = length(x)
    n = length(a)
    p = zeros(T, m, n+1)
    rb1 = one(T) / b[1]
    rb2 = one(T) / b[2]
    for i = 1:m
        p[i,1] = rb1
        p[i,2] = rb2 * ( x[i] - a[1] ) * p[i,1]
    end
    for j = 2:n
       rb = one(T) / b[j+1]
       for i = 1:m
           p[i,j+1] = rb * ( (x[i]-a[j]) * p[i,j]
                                - b[j] * p[i,j-1] )
       end
    end
    return p
end

end


text:
module IntegratedHU

include("./mass_hu.jl")

export
    # Export mass_hu functions
    mass_hu

end


text:
using Test

@info "Notice: These tests check paths some of which are supposed to throw warnings. Others of which are not. Some attention should be paid when reading the output."

tests = [
    "util",
    "locations",
    "main",
    "preupload"
]

for filename in tests
    @testset "$filename" begin
        include(filename * ".jl")
    end
end

examples =  [
    "examples.jl",
    "examples_manual.jl",
]

@testset "examples" for fn in examples
    @testset "$fn" begin
        tempdir = mktempdir()
        try
            @info("sending all datadeps to $tempdir")
            withenv("DATADEPS_LOAD_PATH"=>tempdir,
                    "DATADEPS_NO_STANDARD_LOADPATH"=>true) do
                @testset "download and use" begin
                    include(fn)
                end
                withenv("DATADEPS_DISABLE_DOWNLOAD"=>"true") do
                    @testset "use already downloaded" begin
                        include(fn)
                    end
                end
            end
        finally
    		try
    			@info("removing $tempdir")
                cd(@__DIR__)  # Ensure not currently in directory being deleted
                rm(tempdir, recursive=true, force=true)
    		catch err
    			@warn("Something went wrong with removing $tempdir")
    			@warn(err)
    		end
        end
    end
end


text:
"""
    splitobs(n::Int; at) -> Tuple

Pre-compute the indices for two disjoint subsets and return them
as a tuple of two ranges. The first range will span the first
`at` fraction of possible indices, while the second range will
cover the rest. These indices are applicable to any data
container of size `n`.

```julia
julia> splitobs(100, at=0.7)
(1:70, 71:100)
```

A tuple `at` can be passed for multiple splits:

```julia
julia> splitobs(100, at=(0.1, 0.4))
(1:10, 11:50, 51:100)
```
"""
splitobs(n::Int; at) = _splitobs(n, at)

_splitobs(n::Int, at::Integer) = _splitobs(n::Int, at / n)
_splitobs(n::Int, at::NTuple{N, <:Integer}) where {N} = _splitobs(n::Int, at ./ n)

_splitobs(n::Int, at::Tuple{}) = (1:n,)

function _splitobs(n::Int, at::AbstractFloat)
    0 <= at <= 1 || throw(ArgumentError("the parameter \"at\" must be in interval (0, 1)"))
    n1 = clamp(round(Int, at*n), 0, n)
    (1:n1, n1+1:n)
end

function _splitobs(n::Int, at::NTuple{N,<:AbstractFloat}) where N
    at1 = first(at)
    a, b = _splitobs(n::Int, at1)
    n1 = a.stop
    n2 = b.stop
    at2 = Base.tail(at) .* n ./ (n2 - n1)
    rest = map(x -> n1 .+ x, _splitobs(n2-n1, at2))
    return (a, rest...)
end

"""
    splitobs(data; at) -> Tuple

Split the `data` into multiple subsets proportional to the
value(s) of `at`.

Note that this function will perform the splits statically and
thus not perform any randomization. The function creates a
`NTuple` of data subsets in which the first N-1 elements/subsets
contain the fraction of observations of `data` that is specified
by `at`.

For example, if `at` is a `Float64` then the return-value will be
a tuple with two elements (i.e. subsets), in which the first
element contains the fracion of observations specified by `at`
and the second element contains the rest. In the following code
the first subset `train` will contain the first 70% of the
observations and the second subset `test` the rest.

```julia
train, test = splitobs(X, at=0.7)
```

If `at` is a tuple of `Float64` then additional subsets will be
created. In this example `train` will have the first 50% of the
observations, `val` will have next 30%, and `test` the last 20%

```julia
train, val, test = splitobs(X, at=(0.5, 0.3))
```

It is also possible to call `splitobs` with multiple data
arguments as tuple, which all must have the same number of total
observations. This is useful for labeled data.

```julia
train, test = splitobs((X, y), at=0.7)
(x_train,y_train), (x_test,y_test) = splitobs((X, y), at=0.7)
```

If the observations should be randomly assigned to a subset,
then you can combine the function with `shuffleobs`

```julia
# This time observations are randomly assigned.
train, test = splitobs(shuffleobs((X, y)), at=0.7)
```
"""
function splitobs(data; at)
    n = numobs(data)
    map(idx -> obsview(data, idx), splitobs(n; at))
end


text:
function drgenes(d::AbstractVector{T}) where T<:Real
    d²=d.^2
    fc=d²./mean(d²)
    χ² = Chisq(1)
    pVals = ccdf.(χ², fc)
    pAdjusted = MultipleTesting.adjust(pVals, BenjaminiHochberg())
    return fc,pVals,pAdjusted
end

max_stars_repo_name: jamesjcai/ScTenifold

max_issues_repo_name: jamesjcai/ScTenifold

max_forks_repo_name: jamesjcai/ScTenifold


text:
using LinearAlgebra

abstract type AbstractTrueOnlineTraceOptimizer <: AbstractTraceOptimizer end

mutable struct TrueOnlineTraceOptimizer{T, TF} <: AbstractTrueOnlineTraceOptimizer where {T<:Real,TF<:AbstractFuncApprox}
    f::TF
    α::T
    γ::T
    λ::T
    e::Array{T}
    g::Array{T}
    Δ::Array{T}
    function TrueOnlineTraceOptimizer(fun::AbstractFuncApprox, α::T, γ::T, λ::T) where {T<:Real}
        tp = eltype(get_params(fun))
        s = size(get_params(fun))
        new{tp, typeof(fun)}(fun, tp(α), tp(γ), tp(λ), zeros(tp, s), zeros(tp, s), zeros(tp, s))
    end
end

function update!(opt::TrueOnlineTraceOptimizer{T}, delta::T, ϕ::Array{T}) where {T<:Real}
    opt.e .= (opt.γ * opt.λ) .* opt.e .+ opt.α .* ϕ .- (opt.α * opt.γ * opt.λ * dot(vec(opt.e), vec(ϕ))) .* ϕ
    if !isnan(delta)
        opt.g .= opt.e .* delta .- (opt.α * dot(vec(opt.Δ), ϕ)) .* ϕ
        opt.Δ .= opt.g
        add_to_params!(opt.f, opt.g)
    end
end

function new_episode!(opt::AbstractTrueOnlineTraceOptimizer)
    fill!(opt.e, 0.)
    fill!(opt.g, 0.)
    fill!(opt.Δ, 0.)
end

function clone(opt::TrueOnlineTraceOptimizer)
    opt2 = TrueOnlineTraceOptimizer(opt.f, opt.α, opt.γ, opt.λ)
    opt2.e .= deepcopy(opt.e)
    opt2.g .= deepcopy(opt.g)
    return opt2
end

function clone(opt::TrueOnlineTraceOptimizer, f::AbstractFuncApprox)
    opt2 = TrueOnlineTraceOptimizer(f, opt.α, opt.γ, opt.λ)
    opt2.e .= deepcopy(opt.e)
    opt2.g .= deepcopy(opt.g)
    return opt2
end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
## FastForward (FF) delete-relaxation heuristic ##
export FFHeuristic

"Precomputed domain information for FF heuristic."
struct FFCache
    domain::Domain # Preprocessed domain
    axioms::Vector{Clause} # Preprocessed axioms
    preconds::Dict{Symbol,Vector{Vector{Term}}} # Preconditions in DNF
    additions::Dict{Symbol,Vector{Term}} # Action add lists
end

"FastForward (FF) delete-relaxation heuristic."
mutable struct FFHeuristic <: Heuristic
    cache::FFCache # Precomputed domain information
    pre_key::UInt64 # Key to check if information needs to be precomputed again
    FFHeuristic() = new()
end

Base.hash(heuristic::FFHeuristic, h::UInt) = hash(FFHeuristic, h)

function precompute!(h::FFHeuristic,
                     domain::Domain, state::State, spec::Specification)
    # Check if cache has already been computed
    if is_precomputed(h, domain, state, spec) return h end
    h.pre_key = objectid(domain) # Precomputed data is unique to each domain
    domain = copy(domain) # Make a local copy of the domain
    # Preprocess axioms
    axioms = regularize_clauses(domain.axioms) # Regularize domain axioms
    axioms = [Clause(ax.head, [t for t in ax.body if t.name != :not])
              for ax in axioms] # Remove negative literals
    domain.axioms = Clause[] # Remove axioms so they do not affect execution
    # Preprocess actions
    preconds = Dict{Symbol,Vector{Vector{Term}}}()
    additions = Dict{Symbol,Vector{Term}}()
    for (act_name, act_def) in domain.actions
        # Convert preconditions to DNF without negated literals
        conds = get_preconditions(act_def; converter=to_dnf)
        conds = [c.args for c in conds.args]
        for c in conds filter!(t -> t.name != :not, c) end
        preconds[act_name] = conds
        # Extract additions from each effect
        diff = effect_diff(act_def.effect)
        additions[act_name] = diff.add
    end
    h.cache = FFCache(domain, axioms, preconds, additions)
    return h
end

function is_precomputed(h::FFHeuristic,
                        domain::Domain, state::State, spec::Specification)
    return (isdefined(h, :cache) && objectid(domain) == h.pre_key)
end

function compute(h::FFHeuristic,
                 domain::Domain, state::State, spec::Specification)
    # Precompute if necessary
    if !is_precomputed(h, domain, state, spec)
        precompute!(h, domain, state, spec) end
    @unpack cache = h
    @unpack domain = cache
    @unpack types, facts = state
    goals = get_goal_terms(spec)
    # Initialize fact levels and achievers in a GraphPlan-style planning graph
    levels = Dict{Term,Int}(f => 1 for f in facts)
    achievers = Dict{Term,Vector{Term}}(f => Term[] for f in facts)
    cur_level = 1
    while true
        facts = Set(keys(levels))
        state = State(types, facts, Dict{Symbol,Any}())
        # Break out of loop once all goals are achieved
        if is_goal(spec, domain, state) break end
        cur_level += 1
        # Add all one-step derivations of domain axioms
        for ax in cache.axioms
            _, subst = resolve(ax.body, [Clause(f, []) for f in facts])
            for s in subst
                derived = substitute(ax.head, s)
                if haskey(levels, derived) continue end
                levels[derived] = cur_level
                achievers[derived] = Term[]
            end
        end
        # Add effects of available actions
        actions = available(state, domain)
        for act in actions
            act_vars = domain.actions[act.name].args
            subst = Subst(var => val for (var, val) in zip(act_vars, act.args))
            additions = [substitute(a, subst) for a in cache.additions[act.name]]
            for fact in additions
                if get(levels, fact, Inf) < cur_level continue end
                levels[fact] = cur_level # Update fact level
                push!(get!(achievers, fact, Term[]), act) # Add achieving actions
            end
        end
        # Terminate if we hit a fixed point before satisfying all goals
        if length(levels) == length(facts) && keys(levels) == facts
            return Inf end
    end
    # Backout relaxed partial-order plan from the planning graph
    plan = Set{Term}[]
    while cur_level > 1
        pushfirst!(plan, Set{Term}())
        next_goals = Set{Term}()
        for goal in goals
            if levels[goal] < cur_level
                push!(next_goals, goal) # Handle at next level
                continue
            end
            actions = achievers[goal] # Get all achieving actions
            if length(actions) == 0 continue end # No achievers, continue
            # TODO: match actions with their achieving preconditions
            act = actions[1] # Select the first achiever
            push!(plan[1], act) # Add to plan
            preconds = cache.preconds[act.name][1] # Select the first DNF clause
            act_vars = domain.actions[act.name].args
            s = Subst(var => val for (var, val) in zip(act_vars, act.args))
            preconds = [substitute(p, s) for p in preconds]
            union!(next_goals, preconds) # Handle preconds at next level
        end
        goals = next_goals
        cur_level -= 1
    end
    # Return length of relaxed plan
    # TODO : Return helpful actions
    return sum(length.(plan))
end


text:
struct IntervalCensoredSparseModel{T <: BlasReal}
    basemodel :: IntervalCensoredModel{T} # low-dimensional base model
    X         :: AbstractMatrix{T}        # predictors for variable selection
    βx        :: Vector{T}                # regression coefficients for X
    k         :: Int                      # ||βₓ||₀ = number of nonzero coefficients
    # working arrays
    ∇βx       :: Vector{T}
    ηx        :: Vector{T}
    storage_m :: Vector{T}
    storage_px :: Vector{T}
    storage_pz :: Vector{T}
end

function IntervalCensoredSparseModel(
    basemodel :: IntervalCensoredModel{T},
    X         :: AbstractMatrix{T},
    k         :: Integer
    ) where {T <: BlasReal}
    n, m = length(basemodel.C), length(basemodel.ts)
    px, pz = size(X, 2), length(basemodel.β)
    @assert size(X, 1) == n "rows of X does not match that of Z"
    βx         = Vector{T}(undef, px)
    ∇βx        = Vector{T}(undef, px)
    ηx         = Vector{T}(undef, n)
    storage_m  = Vector{T}(undef, m)
    storage_px = Vector{T}(undef, px)
    storage_pz = Vector{T}(undef, pz)
    # constructor
    IntervalCensoredSparseModel{T}(
        basemodel, X, βx, k,
        ∇βx, ηx, storage_m, storage_px, storage_pz)
end

function iht!(
    icsm    :: IntervalCensoredSparseModel{T};
    maxiter :: Integer = 100_000,
    reltol  :: Real = 1e-4,
    verbose :: Bool = false,
    printfreq :: Integer = 10,
    debiasing :: Bool = false) where T <: Real
    # fit the base model
    @info "Fit the base model:"
    icsm.basemodel.isfitted[1] || fit!(icsm.basemodel)
    display(icsm.basemodel)
    # ICM for debiasing model
    n, pz, px = length(icsm.basemodel.C), length(icsm.basemodel.β), length(icsm.βx)
    if debiasing
        icm_debias = IntervalCensoredModel(
            [icsm.basemodel.Z zeros(T, n, icsm.k)],
            icsm.basemodel.L,
            icsm.basemodel.R)
    end
    # start βₓ from all 0s
    fill!(icsm.βx, 0)
    logl = loglikelihood!(icsm.basemodel, true, true)
    verbose && @info("IHT iterations:")
    verbose && println("iter=0, logl=$logl")
    for iter in 1:maxiter
        # update λ₀
        stepsize = abs2(norm(icsm.basemodel.∇λ₀)) / dot(icsm.basemodel.∇λ₀,
        mul!(icsm.storage_m, icsm.basemodel.Hλ₀λ₀, icsm.basemodel.∇λ₀))
        for (j, λ₀j) in enumerate(icsm.basemodel.λ₀)
            icsm.basemodel.λ₀[j] =
            max(icsm.basemodel.λ₀[j] + stepsize * icsm.basemodel.∇λ₀[j], 0)
        end
        cumsum!(icsm.basemodel.Λ₀, icsm.basemodel.λ₀)
        # update βz
        stepsize = abs2(norm(icsm.basemodel.∇β)) / dot(icsm.basemodel.∇β,
        mul!(icsm.storage_pz, icsm.basemodel.Hββ, icsm.basemodel.∇β))
        icsm.basemodel.β .+= stepsize .* icsm.basemodel.∇β
        mul!(icsm.basemodel.η, icsm.basemodel.Z, icsm.basemodel.β)
        # update βx
        mul!(icsm.∇βx, transpose(icsm.X), icsm.basemodel.βres)
        mul!(icsm.basemodel.storage_n, icsm.X, icsm.∇βx)
        icsm.basemodel.storage_n .*= icsm.basemodel.glmwt
        stepsize = abs2(norm(icsm.∇βx) / norm(icsm.basemodel.storage_n))
        icsm.βx .+= stepsize .* icsm.∇βx
        # project βx to set sparse set {||βx||₀ ≤ k}
        icsm.storage_px .= abs.(icsm.βx)
        pivot = partialsort!(icsm.storage_px, icsm.k, rev=true)
        for (j, βxj) in enumerate(icsm.βx)
            if abs(βxj) < pivot; icsm.βx[j] = 0; end
        end
        # debiasing step
        βx_sp = sparsevec(icsm.βx)
        if debiasing
            # warm start debiased model
            copyto!(icm_debias.β, icsm.basemodel.β)
            copyto!(icm_debias.β, pz + 1, βx_sp.nzval, 1, icsm.k)
            @views icm_debias.Z[:, pz+1:end] .= icsm.X[:, βx_sp.nzind]
            mul!(icm_debias.η, icm_debias.Z, icm_debias.β)
            icm_debias.expη .= exp.(icm_debias.η)
            copyto!(icm_debias.λ₀, icsm.basemodel.λ₀)
            copyto!(icm_debias.Λ₀, icsm.basemodel.Λ₀)
            # fit debiased model
            fit!(icm_debias)
            # copy debiased estimates back
            copyto!(icsm.basemodel.β, 1, icm_debias.β, 1, pz)
            copyto!(βx_sp.nzval, 1, icm_debias.β, pz + 1, icsm.k)
            icsm.βx .= Vector{T}(βx_sp)
            copyto!(icsm.basemodel.η, icm_debias.η)
            copyto!(icsm.basemodel.expη, icm_debias.expη)
            copyto!(icsm.basemodel.λ₀, icm_debias.λ₀)
            copyto!(icsm.basemodel.Λ₀, icm_debias.Λ₀)
        else
            # update η = ηx + ηz
            icsm.basemodel.η   .+= mul!(icsm.ηx, icsm.X, βx_sp)
            icsm.basemodel.expη .= exp.(icsm.basemodel.η)
        end
        # update gradient and Hessian wrt λ₀, βz
        logl_old = logl
        logl = loglikelihood!(icsm.basemodel, true, true)
        verbose && iter%printfreq == 0 && println("iter=$iter, logl=$logl")
        # check convergence
        if abs(logl - logl_old) < reltol * (abs(logl_old) + 1)
            break
        elseif iter == maxiter
            @warn "maximum iterations reached!"
        end
    end
    icsm
end


text:
"""
Strings
"""
const rchar = "rml"
const rmchar = "rmean"
const rqchar = "rquant"
const wchar = "waic"
const routchar = "rout"
const runchar = "run"
const chichar = "chi"
const chioutchar = "chiout"
const rlastchar = "rlast"
const onchar = "LCON"
const offchar = "LCOFF"
const fishchar = "FISH"
const txtstr = ".txt"

"""
eig_decompose(M)
"""
function eig_decompose(M)
    Meiv = eigen(M)
    return Meiv.values, Meiv.vectors
end

"""
nonzero_rows(T)
Locate all rows that have at least one nonzero element
"""
function nonzero_rows(T)
    n = size(T)[2]
    nonzeros = Array{Int}(undef,0)
    for  a in 1:size(T)[1]
        if T[a,:] != zeros(n)
            push!(nonzeros,a)
        end
    end
    return nonzeros
end

normalize_histogram(hist) = hist/sum(hist)

"""
var_update(vartuple, newValue)
iteration step for online algorithm to compute
mean and variance
"""
function var_update(vartuple, newValue)
    count, mean, M2 = vartuple
    count += 1
    delta = newValue - mean
    mean .+= delta / count
    delta2 = newValue - mean
    M2 .+= delta .* delta2
    return count, mean, M2
end
"""
cov_update(covtuple,newValue1,newValue2)
"""
function cov_update(covtuple,newValue1,newValue2)
    count, meanx, meany, C = covtuple
    count += 1
    dx = newx - meanx
    meanx .+= dx / count
    meany .+= (newy - meany) / count
    C .+= dx * (newy - meany)
    return count, meanx, meany, C
end
"""
finalize(vartuple)
Retrieve the mean, variance and sample variance from an aggregate
"""
function finalize(vartuple)
    count, mean, M2 = vartuple
    if count < 2
        return mean, M2
    else
        return mean, M2 ./ (count - 1)
    end
end

"""
truncatehistogram(x::Array{Int,1},yield::Float64,nhistmax::Int)
"""
function truncate_histogram(x::Array,yield::Float64,nhistmax::Int)
    counts = sum(x)
    if counts > 0
        nhist = 1
        ratio = 0.
        # Only keep yield of all data to limit support of histogram
        while ratio < yield && nhist <= nhistmax
            ratio = sum(x[1:nhist])/counts
            nhist += 1
        end
        return x[1:min(nhist,length(x))]
    else
        return 0
    end
end

function combine_histogram(x1::Array,x2::Array)
    n = min(length(x1),length(x2))
    x1[1:n] + x2[1:n]
end

pooled_mean(means,counts) = counts'*means/sum(counts)
pooled_variance(vars,counts) = (counts .- 1)'*vars/sum(counts .- 1)
pooled_std(std,counts) = sqrt(pooled_variance(std.^2,counts))

var_ratio(mua,mub,vara,varb,cov) = mua^2/mub^2*(vara/mua^2 - 2*cov/(mua*mub) + varb/mub^2)

"""
decimal(x::Vector)
convert digits of base 3 number to base 10
"""
function decimal(x::Vector,base::Int=3)
    nr = length(x)
    x' * base .^ collect(0:nr-1)
end
"""
findbase(l,n,nr)
Find the number of G states for transition rate matrix of size l
"""
function findbase(l,n,nr)
    if l == (n+1)*3^nr
        return 3
    elseif l == (n+1)*2^nr
        return 2
    else
        throw("wrong length")
    end
end

"""
KL(x,y)

Kullback-Leibler distance
"""
KL(x,y) = sum(x .* (log.(max.(x,eps(Float64))) - log.(max.(y,eps(Float64)))))

"""
trim(h::Array,nh::Array)

"""
function trim(h::Array,nh::Array)
    for i in eachindex(h)
        h[i] = h[i][1:nh[i]]
    end
    return h
end

"""
LogNormal_array(param,cv)
Prior distribution arrays
"""
LogNormal_array(param,cv) = distributionarray(param,cv,LogNormal)
Gamma_array(param,cv) = distributionarray(param,cv,Gamma)

function LogNormalBeta_array(param,cv,ind)
    if ind == 1
        d = [setBeta(param[ind],cv[ind])]
    else
        barind = 1:ind-1
        d = LogNormalarray(param[barind],cv[barind])
        push!(d,setBeta(param[ind],cv[ind]))
    end
    return d
end
"""
distribution_array(param,cv,dist)

fill an array with dist(param,cv)
"""

function distributionBeta_array(param::Vector,cv::Vector,ind::Int,dist=LogNormal)
    if ind == 1
        d = [Beta_meancv(param[ind],cv[ind])]
    else
        barind = 1:ind-1
        d = distribution_array(param[barind],cv[barind],dist)
        push!(d,Beta_meancv(param[ind],cv[ind]))
    end
    return d
end
function distribution_array(param::Vector,cv,dist=LogNormal)
    d = []
    for i in eachindex(param)
        push!(d,dist(param[i],cv[i]))
    end
    return d
end

"""
Gamma_meancv(mean,cv)
LogNormal_meancv(mean,cv)
Beta_meancv(mean,cv)

Reparameterized distributions to mean and coefficient of variation arguments
"""
Gamma_meancv(mean,cv) = Gamma(1/cv^2,mean * cv^2)
LogNormal_meancv(mean,cv) = LogNormal(mulognormal(mean,cv), sigmalognormal(cv))
function Beta_meancv(m,cv)
    cv2 = cv^2
    fac = (1-m)/cv2/m - 1
    if fac <= 0.
        fac = 2/m
    end
    alpha = m*fac
    beta = (1-m)*fac
    Beta(alpha,beta)
end
"""
sigmalognormal(cv)
mulognormal(mean,cv)

compute LogNormal mu and sigma parameters
given desired mean and coefficient of variation
"""
sigmalognormal(cv) = sqrt.(log.(1 .+ cv .^ 2))

mulognormal(mean,cv) = log.(mean) - .5*log.(1 .+ cv .^ 2)

"""
Parameter information for GR models
"""
function num_rates(G,R)
    n = G - 1
    2*n + 2*R + 2
end
function Grange(G)
    n = G - 1
    1 : 2*n
end
function initiation(G)
    n = G - 1
    2*n + 1
end
function Rrange(G,R)
    n = G - 1
    2*n + 2 : 2*n + R + 1
end
function Srange(G,R)
    n = G - 1
    2*n + R + 2 : 2*n + 2*R + 1
end

mean_histogram(x) = (collect(1:length(x)) .- 1)' * x/sum(x)
m2_histogram(x) = ((collect(1:length(x)) .- 1).^2)' * x/sum(x)
function var_histogram(x)
    m2_histogram(x) - mean_histogram(x)^2
end
# mean_histogram(x) = ((collect(length(x)) .- 1)' * x)/sum(x)
tstat_2sample(x1,x2) = (mean_histogram(x1) - mean_histogram(x2))/(sqrt(var_histogram(x1)/length(x1) + var_histogram(x2)/length(x2)))

log_2sample(x1,x2) = log(mean_histogram(x1)) - log(mean_histogram(x2))
delta_2sample(x1,x2) = mean_histogram(x1)-mean_histogram(x2)
delta_2frac(x1,x2) = delta_2sample(x1,x2)/mean_histogram(x1)

function mediansmooth(xin,window)
    x = copy(xin)
    halfwin = div(window-1,2)
    for i in 1:length(x)-window+1
        x[i+halfwin] = median(x[i:i+window-1])
    end
    x
end

"""
plot_histogram()

functions to plot data and model predicted histograms

"""
function plot_histogram(data::RNAData{Vector{Int64}, Vector{Array}},model::GMlossmodel)
    h=likelihoodarray(model.rates,data,model)
    for i in eachindex(h)
        figure()
        plot(h[i])
        plot(normalize_histogram(data.histRNA[i]))
        savefig(string(i))
    end
    return h
end
function plot_histogram(data::AbstractRNAData{Array{Array,1}},model)
    h=likelihoodarray(model.rates,data,model)
    figure(data.gene)
    for i in eachindex(h)
        plot(h[i])
        plot(normalize_histogram(data.histRNA[i]))
        savefig(string(i))
    end
    return h
end
function plot_histogram(data::AbstractRNAData{Array{Float64,1}},model)
    h=likelihoodfn(get_param(model),data,model)
    figure(data.gene)
    plot(h)
    plot(normalize_histogram(data.histRNA))
    return h
end

function plot_histogram(data::RNALiveCellData,model)
    h=likelihoodtuple(model.rates,data,model)
    figure(data.gene)
    plot(h[1])
    plot(normalize_histogram(data.OFF))
    plot(h[2])
    plot(normalize_histogram(data.ON))
    figure("FISH")
    plot(h[3])
    plot(normalize_histogram(data.histRNA))
    return h
end

function plot_histogram(data::TransientRNAData,model::AbstractGMmodel)
    h=StochasticGene.likelihoodarray(model.rates,data,model)
    for i in eachindex(h)
        figure(data.gene *":T" * "$(data.time[i])")
        plot(h[i])
        plot(normalize_histogram(data.histRNA[i]))
    end
    return h
end

function plot_histogram(data::RNAData,model::AbstractGMmodel)
    h=StochasticGene.likelihoodfn(get_param(model),data,model)
    figure(data.gene)
    plot(h)
    plot(normalize_histogram(data.histRNA))
    return h
end

"""
residenceprob_G(file,G,header)
residenceprob_G(r,n)

Residence probability of G states
given by exact steady state solution
of master equation

"""
function residenceprob_G(file::String,G,header=false)
    r = get_all_rates(file,header)
    m = size(r)[1]
    p = Array{Any,2}(undef,m,G+1)
    n = G-1
    for i in 1:m
        p[i,1] = r[i,1]
        p[i,2:end] = residenceprob_G(r[i,2:2*n+1],n)
    end
    return p
end

function residenceprob_G(r::Vector,n::Int)
    Gss = Array{Float64,2}(undef,1,n+1)
    Gss[1,1] = 1.
    for k in 1:n
        Gss[1,k+1] = Gss[1,k]*r[2*k-1]/r[2*k]
    end
    Gss ./= sum(Gss)
end

"""
splicesiteusage()

splice site usage probability is the probabilty of ejection times
the probability of not ejection prior to that point
"""
splicesiteusage(model::GRSMmodel) = splicesiteusage(model.rates,model.G-1,model.R)
function splicesiteusage(r::Vector,n::Int,nr::Int)
    nu = get_nu(r,n,nr)
    eta = get_eta(r,n,nr)
	ssf = zeros(nr)
	survival = 1
	for j in 1:nr
		ssf[j] = eta[j]/(nu[j+1]+eta[j])*survival
		survival *= nu[j+1]/(nu[j+1]+eta[j])
	end
	return ssf
end

"""
burstsize(n,nr,r)
Burst size distribution of GRS  model
for total pre-mRNA occupancy and
unspliced (visible) occupancy
obtained by marginalizing over conditional steady state distribution
"""
burstsize(model::GRSMmodel) = burstsize(model.G-1,model.R,model.rates)
function burstsize(n::Int,nr::Int,r::Vector)
    T =  mat_GSR_T(r,n,nr)
    pss = normalized_nullspace(T)
	Rss = zeros(nr)
	Rssvisible = zeros(nr)
	ssf = zeros(nr)
	asum = 0
	for w=1:nr
		for i = 1:n+1, z = 1:3^nr
			zdigits = digits(z-1,base=3,pad=nr)
			a = i + (n+1)*(z-1)
			if sum(zdigits .== 2) == w
				Rssvisible[w] += pss[a]
			end
			if sum(zdigits .> 0) == w
				Rss[w] += pss[a]
			end
		end
	end
	Rss ./= sum(Rss), Rssvisible ./= sum(Rssvisible)
end

"""
teststeadystatemodel(model,nhist)

Compare chemical master solution to Gillespie simulation for steadystate mRNA distribution

"""
function teststeadystatemodel(model::AbstractGMmodel,nhist)
    G = model.G
    r = model.rates
    g1 = steady_state(r[1:2*G],G-1,nhist,model.nalleles)
    g2 = telegraph(G-1,r[1:2*G],10000000,1e-5,nhist,model.nalleles)
    return g1,g2
end

# function online_covariance(data1, data2)
#     meanx = meany = C = n = 0
#     for x in data1, y in data2
#         n += 1
#         dx = x - meanx
#         meanx += dx / n
#         meany += (y - meany) / n
#         C += dx * (y - meany)
#         population_covar = C / n
#         # Bessel's correction for sample variance
#         sample_covar = C / (n - 1)
#     end
# end

max_stars_repo_name: nih-niddk-mbs/StochasticGene

max_issues_repo_name: nih-niddk-mbs/StochasticGene

max_forks_repo_name: nih-niddk-mbs/StochasticGene


text:
#Note that the duration is software-timed so not precise
function ttl_pulse(;nsecs = 1.0, line_name = "Port0/line0")
    otsk = NIDAQ.digital_output(DEFAULT_DEVICE*line_name)
    try
        start(otsk)
        write(otsk, [UInt32(true)]) #TODO: element type here is device-dependent (depends on number of lines per port)
        sleep(nsecs)
        write(otsk, [UInt32(false)])
    finally
        stop(otsk)
        clear(otsk)
    end
    return 0
end

function sampleclockstring(sig1)
    if isoutput(sig1) && isdigital(sig1)
        return "do/SampleClock"
    elseif !isoutput(sig1) && isdigital(sig1)
        return "di/SampleClock"
    elseif isoutput(sig1) && !isdigital(sig1)
        return "ao/SampleClock"
    else
        return "ai/Sampleclock"
    end
end

function check_valid(sigs, sync_clocks::Bool, skip_ii_validation::Bool)
    if isempty(sigs)
        error("Empty signal list")
    end
    sig1 = first(sigs)
    if rig_name(sig1) == "dummy-6002" && sync_clocks
        error("The usb-6002 device does not support clock synchronization.  Please set the sync_clocks kwarg to false")
    end
    #Don't require a "sufficient" set of signals for an imaging experiment until this is a fully working alternative to Imagine (easier for testing)
    if !skip_ii_validation
        ImagineInterface.validate_all(sigs; check_is_sufficient = false)
    end
    outs = getoutputs(sigs)
    if isempty(outs)
        error("No output channels found, so unable to guess the number of samples to acquire.  Use the record_signals function instead.")
    end
    nothing
end

function run_imagine(base_name::T, sigs::Vector{S};
				   ai_trig_dest = "disabled",
				   ao_trig_dest = "disabled",
				   trigger_source = "Port0/Line0",
				   sync_clocks = true,
				   run_locally = false,
				   skip_validation = false) where {T<:AbstractString, S<:ImagineSignal}
    check_valid(sigs, sync_clocks, skip_validation)
    sig1 = first(sigs)
    ins = getinputs(sigs)
    outs = getoutputs(sigs)
    #determine shared clock source from first signal
    clk_src = sync_clocks ? sampleclockstring(sig1) : ""
    dev = DEFAULT_DEVICE
    if length(WORKERS) < NPERSISTENT_WORKERS #currently we always keep 4 workers ready.  Could instead ready them on demand but so far this seems like a win.
        to_add = NPERSISTENT_WORKERS - length(WORKERS)
        new_workers = add_workers(to_add, dev)
    end
    sigs_out = ImagineSignal[]
    rchans = RemoteChannel{Channel{Int}}[]
    rrs = []
    ids = Int[]
    nsamps = length(first(outs))
    (rchns_o, rr) = output_signals(outs; trigger_terminal = ao_trig_dest, clock = clk_src, run_locally = run_locally)
    append!(rchans, rchns_o)
    append!(rrs, rr)
    if !isempty(ins)
        #TODO: insert more logic for digital signals, di_trig_dest
        (rchns_i, rr) = record_signals(base_name, ins, nsamps; trigger_terminal = ai_trig_dest, clock = clk_src, run_locally = run_locally)
        append!(rchans, rchns_i)
        append!(rrs, rr)
    end
    print("Preparing all DAQ tasks...\n")
    finished = falses(length(rchans))
    #NOTE: this method of monitoring errors seems fragile.  The behavior of isready(::Future) may be changing in v0.7
    #see docs for isready()
    while !all(finished)
        for i in findall(.!(finished))
            c = rchans[i]
            if !isready(c)
                if isready(rrs[i]) #shouldn't be ready yet, must be error
                    print("Error while preparing experiment:\n")
                    fetch(rrs[i])
                end
            else
                push!(ids, take!(c))
                finished[i] = true
            end
        end
    end
    if rig_name(first(sigs)) == "dummy-6002"
        sleep(4.0) #this shouldn't be necessary, but it is (a digital trigger for AI was found inneffective immediately after the CfgDigEdgeStartTrig function returned with usb 6002)
    end
    print("Triggering tasks...\n")
    ttl_pulse(; line_name = trigger_source) #P0.0 is wired to PFI0 and PFI1 for testing with usb 6002
    rslts = Vector{Any}(undef, length(rrs))
    finished = falses(length(rrs))

    while !all(finished)
        for i = 1:length(rrs)
            if !finished[i] && isready(rrs[i])
                try
                    rslts[i] = fetch(rrs[i])
                    append!(sigs_out, rslts[i])
                catch err
                    rslts[i] = err
                    print("ERROR ")
                    showerror(Base.stderr, err) #rslts[i] #can also test later with isa(rslts[i], RemoteException)
                end
                finished[i] = true
            end
        end
        sleep(0.1)
    end

    if any(x->isa(x, RemoteException), rslts)
        error("One or more DAQ threads crashed.  See details in output above.")
    end

    free_workers(ids)
    return sigs_out
end

#ttl_pulse(; nsecs = 1.0, dev=DEFAULT_DEVICE, line_name="Port0/Line0") = _ttl_pulse(nsecs, line_name)

#This version spawns a worker
#function ttl_pulse(; nsecs = 1.0, dev=DEFAULT_DEVICE, line_name="Port0/Line0") = _ttl_pulse(nsecs, line_name)
#    id = get_worker(;dev=dev)
#    sendto([id;], nsecs=nsecs, line_name=line_name)
#    rslt = @spawnat id _ttl_pulse(nsecs, line_name)
#    fetch(rslt)
#    free_workers(id)
#end

#note: this shows whether start and stop triggering are supported on the device.
#NIDAQ.GetDevAOTrigUsage(r, rslt)
#NIDAQ.GetDevAITrigUsage(r, rslt)
#For USB-6002 devices that clock synchronization is not supported (see list of terminals).  bummer.
#strategy:  for the usb-6002 we can synchronize AO and AI starts by connecting the start-triggers with the PFI0 or PFI1 pins.  Note that the tasks will desynchronize over time since their
#clocks aren't synchronized
#note: with other boards we shouldn't need the PFI pins, maybe can just use the ao or ai start trigger terminals.  Then one needs to be triggered by software, the others can be triggered by its start sig.
#USB-6002 list of terminals (Ctr0Source is mysterious, probably useless) ["/Dev1/PFI0", "/Dev1/PFI1", "/Dev1/ai/StartTrigger", "/Dev1/ao/StartTrigger", "/Dev1/Ctr0Source"]
#will also implement another function that synchronizes the sample clocks of the two (and also of digital channels) for when this is used with a better DAQ
#   TODO: synchronize sample clocks for ai and ao tasks


text:
# This file was generated, do not modify it. # hide
logdensity(Normal(), Lebesgue(ℝ), 1.0)

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


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

@doc raw"""
    ObjsMessageIcons(;
        emoji=nothing,
        image_64=nothing,
    )

    - emoji::String
    - image_64::String
"""
mutable struct ObjsMessageIcons <: SwaggerModel
    emoji::Any # spec type: Union{ Nothing, String } # spec name: emoji
    image_64::Any # spec type: Union{ Nothing, String } # spec name: image_64

    function ObjsMessageIcons(;emoji=nothing, image_64=nothing)
        o = new()
        validate_property(ObjsMessageIcons, Symbol("emoji"), emoji)
        setfield!(o, Symbol("emoji"), emoji)
        validate_property(ObjsMessageIcons, Symbol("image_64"), image_64)
        setfield!(o, Symbol("image_64"), image_64)
        o
    end
end # type ObjsMessageIcons

const _property_map_ObjsMessageIcons = Dict{Symbol,Symbol}(Symbol("emoji")=>Symbol("emoji"), Symbol("image_64")=>Symbol("image_64"))
const _property_types_ObjsMessageIcons = Dict{Symbol,String}(Symbol("emoji")=>"String", Symbol("image_64")=>"String")
Base.propertynames(::Type{ ObjsMessageIcons }) = collect(keys(_property_map_ObjsMessageIcons))
Swagger.property_type(::Type{ ObjsMessageIcons }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsMessageIcons[name]))}
Swagger.field_name(::Type{ ObjsMessageIcons }, property_name::Symbol) =  _property_map_ObjsMessageIcons[property_name]

function check_required(o::ObjsMessageIcons)
    true
end

function validate_property(::Type{ ObjsMessageIcons }, name::Symbol, val)
end