text:
module ConstantCellValuesTests

using Test
using Gridap
using TensorValues
using ..MapsMocks

# Constructors

l = 10
sv = 1.0
scv = ConstantCellNumber(sv,l)

sv2 = 1.1
scv2 = ConstantCellNumber(sv2,l)

test_index_cell_value( scv, fill(sv,l) )

sa = [sv, sv, sv]
sca = ConstantCellArray(sa,l)

test_index_cell_array( sca, fill(sa,l) )

sca = ConstantCellVector(sa,l)

test_index_cell_array( sca, fill(sa,l) )

sa2 = [sv sv; sv sv; sv sv]
sca2 = ConstantCellArray(sa2,l)

test_index_cell_array( sca2, fill(sa2,l) )

sca2 = ConstantCellMatrix(sa2,l)

test_index_cell_array( sca2, fill(sa2,l) )

a = VectorValue(10,10)
b = VectorValue(15,20)
p1 = VectorValue(1,1)
p2 = VectorValue(2,2)
p3 = VectorValue(3,3)
p = [p1,p2,p3]
m = MockMap(a)
r = evaluate(m,p)

cm = ConstantCellMap(m,l)
cp = ConstantCellVector(p,l)
rm = [ r for i in 1:l]
test_index_cell_map_with_index_arg(cm,cp,rm)
cr = evaluate(cm,cp)
@test isa(cr,ConstantCellArray)

scv3 = apply(-,scv,scv2)
@test isa(scv3,ConstantCellValue{Float64})
test_index_cell_value( scv3, fill(-(sv,sv2),l) )

for op in (:+,:-)
  @eval begin
    scv3 = $op(scv,scv2)
    @test isa(scv3,ConstantCellValue{Float64})
    test_index_cell_value( scv3, fill($op(sv,sv2),l) )
  end
end

sca3 = apply(-,sca2,broadcast=true)
@test isa(sca3,ConstantCellArray{Float64,2})
test_index_cell_array( sca3, fill(broadcast(-,sa2),l) )

sca3 = apply(-,sca,sca2,broadcast=true)
@test isa(sca3,ConstantCellArray{Float64,2})
test_index_cell_array( sca3, fill(broadcast(-,sa,sa2),l) )

cm = ConstantCellMap(m,l)
cp = ConstantCellVector(p,l)
rm = [ -r for i in 1:l]
cm2 = -cm
test_index_cell_map_with_index_arg(cm2,cp,rm)
cr = evaluate(cm2,cp)
@test isa(cr,ConstantCellArray)

cm = ConstantCellMap(m,l)
cp = ConstantCellVector(p,l)
rm = [ r.-p for i in 1:l]
cm2 = apply(-,cm,cp,broadcast=true)
test_index_cell_map_with_index_arg(cm2,cp,rm)
cr = evaluate(cm2,cp)
@test isa(cr,ConstantCellArray)

@test scv == scv
@test scv ≈ scv

@test cp == cp
@test cp ≈ cp

ca = ConstantCellNumber(3,5)
sv = 1.0
scv = ConstantCellNumber(sv,l)
cc = reindex(scv,ca)
r = fill(sv,5)
@assert isa(cc,ConstantCellValue)
test_iter_cell_value(cc,r)

f = ( x -> 3*x )
fca = broadcast(f,ca)
@test all(fca .== ca*3)

end # module ConstantCellValuesTests


text:
using DifferentialEquations
using DataStructures
#here t denotes sensitivity, 0.3 is a constant which can be changed later after testing
eqn(i,p,s) =  i * 0.3
i0 = 0.5  #intial brightness
#=this would just generate a sol object from which any value of intensity corresponding to
a certain sensitivity value can be generated
eg :  sol(115)
=#
s_span = (0.0,10.0)
#I've taken a value of 1/3.0; might have to change this after testing
prob = ODEProblem(eqn,i0,s_span)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, dt=0.1, save_everystep=true)
print(sol)
#=
On printing sol, you'll find 2 lists, one for t, one for u
here t is s
u is intensity
=#
arsol = Array(sol)
println(arsol)        #solution array

#############
#This part should logically be in  a different script. Will break down while integration
#############
#I've assumed my current intensity to be 1.0, hence i computed its index
searchsortedfirst(arsol,1.0) #result 12
#Assumed required intensity to be 7.0 (Note : intensity can't be more than 10); will have to discuss solution handling
searchsortedfirst(arsol,7.0) # result 34
 #now iterate through 12 to 34 and change brightness accordingly

ext: jl


text:
using Makie
using FileIO, Colors

"""
mp_init()

makie plotの初期化

"""
function mp_init()
    scene = Scene()
end

"""
mp_earth(earth_radius=1)

地球を表示する関数

# Arguments
- （earth_radius:地球の半径）
"""
function mp_earth(;earth_radius=1)
    earth = try
        load(download("https://svs.gsfc.nasa.gov/vis/a000000/a002900/a002915/bluemarble-2048.png"))
    catch e
        @warn("Downloadinging the earth failed. Using random image, so this test will fail! (error: $e)")
        rand(RGBAf0, 100, 100) # don't error test when e.g. offline
    end
    m = GLNormalUVMesh(Sphere(Point3f0(0), earth_radius), 60)
    mesh(m, color = earth, shading = false)
end

"""
mp_basic_coordinate(scene,r=[0,0,0])

基底を示す矢印を表示・重ねて表示

# Arguments
- myscene：上書き対象となるscene
- （r:座標系の原点）
- （arrow_size:矢印の長さ）

"""
function mp_basic_coordinate(myscene; r=[0,0,0],arrow_size=1)
    arrows!(myscene,[r[1],r[1],r[1]],[r[2],r[2],r[2]],[r[3],r[3],r[3]],[arrow_size,0,0],[0,arrow_size,0],[0,0,arrow_size],normalize=false,linewidth=5,arrowsize=arrow_size/10)
end

function mp_r(myscene, r)
    arrows!(myscene, [0],[0],[0],r[1],r[2],r[3])
end

function mp_sat(myscene,r)
    mp_basic_coordinate(myscene)
    arrows!(myscene,[0],[0],[0],[r[1]],[r[2]],[r[3]],normalize=true,linewidth=5,arrowsize=0.1)
end

"""
mp_orbit(r)

軌道の表示
# Arguments
- myscene：上書き対象となるscene
- r_log：軌道の記録

"""
function mp_orbit(myscene, r_log)
    mp_basic_coordinate(myscene,arrow_size=10000)
    for i = 1:size(r_log)[1]
        arrows!(myscene,[0],[0],[0],[r_log[i,1]],[r_log[i,2]],[r_log[i,3]],normalize=false,linewidth=5,arrowsize=0.1)
    end

    Base.prompt("終了します")
end

"""
mp_orbit_normal(r)

軌道の表示
# Arguments
- myscene：上書き対象となるscene
- r_log：軌道の記録

"""
function mp_orbit_normal(myscene, r_log)
    mp_basic_coordinate(myscene,arrow_size=1)
    for i = 1:size(r_log)[1]
        arrows!(myscene,[0],[0],[0],[r_log[i,1]],[r_log[i,2]],[r_log[i,3]],normalize=true,linewidth=1,arrowsize=0.01)
    end

    Base.prompt("終了します")
end

"""
mp_sat_env(r)

軌道面座標系上での各パラメータの表示
# Arguments
- myscene：上書き対象となるscene
- vecs：表示したいベクトル（正規化される）

"""
function mp_sat_env(myscene, vecs,vec1...)
    mp_basic_coordinate(myscene,arrow_size=1)
    for i = 1:size(vecs)[1]
        arrows!(myscene,[0],[0],[0],[vecs[i,1]],[vecs[i,2]],[vecs[i,3]],normalize=true,linewidth=1,arrowsize=0.05)
    end

    Base.prompt("終了します")
end

function mp_delete(myscene,nums)
    for i = 1:size(nums)
        delete!(myscene, myscene[i])
    end
end

#Makie.save("/Users/takuya/Desktop/julia_plots/fig.png")

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
plot(mod::PhyloNetworks.TraitSubstitutionModel) = error("plot not defined for $(typeof(mod)).")

"""
    plot(model::TwoBinaryTraitSubstitutionModel)

Display substitution rates for a trait evolution model for two
possibly dependent binary traits; using `R` via `RCall`.
Adapted from fitPagel functions found in the `R` package `phytools`.

# Examples

```
julia> using PhyloNetworks

julia> m = TwoBinaryTraitSubstitutionModel([2.0,1.2,1.1,2.2,1.0,3.1,2.0,1.1],
                ["carnivory", "noncarnivory", "wet", "dry"])
Substitution model for 2 binary traits, with rate matrix:
                     carnivory-wet    carnivory-dry noncarnivory-wet noncarnivory-dry
    carnivory-wet                *           1.0000           2.0000           0.0000
    carnivory-dry           3.1000                *           0.0000           1.1000
 noncarnivory-wet           1.2000           0.0000                *           2.0000
 noncarnivory-dry           0.0000           2.2000           1.1000                *

julia> plot(m);
```
"""
function plot(object::PhyloNetworks.TwoBinaryTraitSubstitutionModel)
    R"""
    signif<-3
    plot.new()
    par(mar=c(1.1,2.1,3.1,2.1))
    plot.window(xlim=c(0,2),ylim=c(0,1),asp=1)
    """
    R"""
    mtext("Two Binary Trait Substitution Model",side=3,adj=0,line=1.2,cex=1.2)
    arrows(x0=0.15,y0=0.15,y1=0.85,lwd=2,length=0.1)
    arrows(x0=0.2,y0=0.85,y1=0.15,lwd=2,length=0.1)
    arrows(x0=1.6,y0=0.05,x1=0.4,lwd=2,length=0.1)
    arrows(x0=0.4,y0=0.1,x1=1.6,lwd=2,length=0.1)
    arrows(x0=1.8,y0=0.15,y1=0.85,lwd=2,length=0.1)
    arrows(x0=1.85,y0=0.85,y1=0.15,lwd=2,length=0.1)
    arrows(x0=1.6,y0=0.9,x1=0.4,lwd=2,length=0.1)
    arrows(x0=0.4,y0=0.95,x1=1.6,lwd=2,length=0.1)
    text(x=0.175,y=0.95,$(object.label[1]))
    text(x=1.825,y=0.95,$(object.label[2]))
    text(x=1.825,y=0.05,$(object.label[4]))
    text(x=0.175,y=0.05,$(object.label[3]))
    """
    R"""
    text(x=1,y=1,round($(object.rate[5]),signif),cex=0.8)
    """
    R"""
    text(x=1,y=0.85,round($(object.rate[6]),signif),cex=0.8)
    text(x=1.9,y=0.5,round($(object.rate[3]),signif),cex=0.8,srt=90)
    text(x=1.75,y=0.5,round($(object.rate[4]),signif),cex=0.8,srt=90)
    """
    R"""
    text(x=1,y=0,round($(object.rate[8]),signif),cex=0.8)
    text(x=1,y=0.15,round($(object.rate[7]),signif),cex=0.8)
    text(x=0.1,y=0.5,round($(object.rate[2]),signif),cex=0.8,srt=90)
    text(x=0.25,y=0.5,round($(object.rate[1]),signif),cex=0.8,srt=90)
    """
    return nothing
end


text:
immutable Hypergeometric <: DiscreteUnivariateDistribution
    ns::Int # number of successes in population
    nf::Int # number of failures in population
    n::Int  # sample size

    function Hypergeometric(ns::Int, nf::Int, n::Int)
        ns >= 0 || error("ns must be non-negative.")
        nf >= 0 || error("nf must be non-negative.")
        0 < n < ns + nf || error("n must have 0 < n < ns + nf")
        new(ns, nf, n)
    end

    # Note: `convert` ensures that the inputs are integers, otherwise it throws InexactError
    Hypergeometric(ns::Real, nf::Real, n::Real) =
        Hypergeometric(convert(Int, ns), convert(Int, nf), convert(Int, n))
end

@_jl_dist_3p Hypergeometric hyper

# handling support
minimum(d::Hypergeometric) = max(0, d.n - d.nf)
maximum(d::Hypergeometric) = min(d.ns, d.n)
support(d::Hypergeometric) = minimum(d):maximum(d)

islowerbounded(d::Hypergeometric) = true
isupperbounded(d::Hypergeometric) = true

function _probs(d::Hypergeometric, f::Int, l::Int)
    ns = d.ns
    nf = d.nf
    n = d.n
    r = Array(Float64, l - f + 1)
    r[1] = v = pdf(d, f)
    b = f - 1
    if l > f
        for x = f+1:l
            c = ((ns - x + 1) / x) * ((n - x + 1) / (nf - n + x))
            r[x-b] = (v *= c)
        end
    end
    return r
end

probs(d::Hypergeometric) = _probs(d, minimum(d), maximum(d))

function probs(d::Hypergeometric, rgn::UnitRange)
    f, l = rgn[1], rgn[end]
    minimum(d) <= f <= l <= maximum(d) || throw(BoundsError())
    _probs(d, f, l)
end

# properties
mean(d::Hypergeometric) = d.n * d.ns / (d.ns + d.nf)

function var(d::Hypergeometric)
    N = d.ns + d.nf
    p = d.ns / N
    d.n * p * (1.0 - p) * (N - d.n) / (N - 1.0)
end
mode(d::Hypergeometric) = int(floor((d.n+1) * (d.ns+1) / (d.ns+d.nf+2)))

function modes(d::Hypergeometric)
    if (d.ns == d.nf) && mod(d.n, 2) == 1
        [(d.n-1)/2, (d.n+1)/2]
    else
        [mode(d)]
    end
end

skewness(d::Hypergeometric) = (d.nf-d.ns)*sqrt(d.ns+d.nf-1)*(d.ns+d.nf-2*d.n)/sqrt(d.n*d.ns*d.nf*(d.ns+d.nf-d.n))/(d.ns+d.nf-2)
function kurtosis(d::Hypergeometric)
    N = d.ns + d.nf
    a = (N-1) * N^2 * (N * (N+1) - 6*d.ns * (N-d.ns) - 6*d.n*(N-d.n)) + 6*d.n*d.ns*(d.nf)*(N-d.n)*(5*N-6)
    b = (d.n*d.ns*(N-d.ns) * (N-d.n)*(N-2)*(N-3))
    a/b
end


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(13)
    n = 3000
    data = randn(n)
    a = fill(1, n)
    fig = Figure(resolution = (700, 450), font = "sans")
    ax1 = Axis(fig, xlabel = "variable", ylabel = "values",
                xticks = ([1],["normal Distribution"]) )
    violin!(ax1, a, data; width = 0.35,  color = (:yellow, 0.45),
                show_median= true, mediancolor = :navy, strokecolor = :black,
                strokewidth = 1, label = "vertical")
    limits!(ax1, 0,2,-5,5)
    axislegend(ax1, position = :rt)
    fig[1,1] = ax1
    fig
    save(joinpath(@__DIR__, "output", "ViolinPlot.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
module FuzzySets

name() = return "FuzzySets.jl"

using GLMakie
using LaTeXStrings
import Plots
using FLoops

include("Intervals.jl")
include("FuzzySet.jl")
include("FuzzyNumbers.jl")
include("FuzzyVectors.jl")
include("FuzzyMatrices.jl")
include("FuzzyWeightedAverages.jl")
include("FCM.jl")
include("PCM.jl")
include("IT2FCM.jl")
include("Dampening.jl")
include("LFCM.jl")
include("LPCM.jl")

export
    Interval, .., mid, rad,
    FuzzySet, support, core, height, cut, centroid,
    FuzzyNumber, SingletonFuzzyNumber, TriangularFuzzyNumber, TrapezoidalFuzzyNumber, GaussianFuzzyNumber,
    FuzzyVector,
    peak, draw,
    fuzzy_weighted_average,
    euclidean_distance,
    u_uncertainty
end


text:
@testset "descent tests" begin
    @testset "basic usage tests" begin
        data = [rand(5) for _ in 1:20]
        k = 4
        metric = SqEuclidean()
        # vector of vectors
        @inferred nndescent(data, k, metric; max_iters=1)
        @test nndescent(data, k, metric; max_iters=1) isa HeapKNNGraph{Int, Float64}

        # matrix
        data = rand(5, 20)
        @inferred nndescent(data, k, metric; max_iters=1)
        @test nndescent(data, k, metric; max_iters=1) isa HeapKNNGraph{Int, Float64}
    end
    @testset "GraphT method tests" begin
        # test passing graph types to nndescent
        for GraphT in [HeapKNNGraph, LockHeapKNNGraph]
            data = [rand(5) for _ in 1:20]
            k = 4
            metric = Euclidean()
            @test nndescent(GraphT, data, k, metric; max_iters=1) isa GraphT
            # matrix
            data = rand(5, 20)
            @test nndescent(GraphT, data, k, metric; max_iters=1) isa GraphT
        end
        @test true
    end
    test_data = [rand(2) for _ in 1:5]
    test_metric = Euclidean()
    test_inds = Int[2 4 2 1 1;
                    5 3 4 3 4]
    test_dsts = Float64[1. 2. 3. 4. 5.;
                        2. 3. 4. 5. 6.]
    @testset "get_neighbors! tests" begin
        g = HeapKNNGraph(test_data, test_metric, test_inds, test_dsts)
        olds, news = NearestNeighborDescent.get_neighbors!(g)
        @test all(length.(olds) .== 0) # no old edges
        @test news[1] == [2, 4, 5]
        @test news[2] == [1, 3, 4]
        @test news[3] == [2, 4]
        @test news[4] == [1, 2, 3, 5]
        @test news[5] == [1, 4]
        olds, news = NearestNeighborDescent.get_neighbors!(g)
        @test all(length.(news) .== 0)
        @test olds[1] == [2, 4, 5]
        @test olds[2] == [1, 3, 4]
        @test olds[3] == [2, 4]
        @test olds[4] == [1, 2, 3, 5]
        @test olds[5] == [1, 4]
    end

    cube_points = [[0., 0., 0.],
            [0., 0., 1.],
            [0., 1., 0.],
            [0., 1., 1.],
            [1., 0., 0.],
            [1., 0., 1.],
            [1., 1., 0.],
            [1., 1., 1.]]
    cube_true_nn = transpose([2 3 5;
                         1 4 6;
                         1 4 7;
                         2 3 8;
                         1 6 7;
                         2 5 8;
                         3 5 8;
                         4 6 7])

    @testset "local_join tests" begin
        mean_knn_diameter(g) = mean([knn_diameter(g, v) for v in vertices(g)])
        Random.seed!(0)
        data = [rand(5) for _ in 1:8]
        g1 = HeapKNNGraph(data, 3, Euclidean())
        init_inds, init_dsts = knn_matrices(g1)
        g2 = LockHeapKNNGraph(data, Euclidean(), init_inds, init_dsts)
        g1_diam = mean_knn_diameter(g1)

        NearestNeighborDescent.local_join!(g1)
        NearestNeighborDescent.local_join!(g2)
        @test mean_knn_diameter(g1) ≤ g1_diam # test the knn diameter is nonincreasing
        inds1, dsts1 = knn_matrices(g1)
        inds2, dsts2 = knn_matrices(g2)
        @test all(inds1 .== inds2)
        @test all(dsts1 .== dsts2)

    end

end


text:
using GeneticJulia
include("fitnessFunctions.jl")
import Random
println("hola")
println("hola")
println("hola")
"""
rng = Random.MersenneTwister(6)

setIndividualType(CGPGenotype)
setRNG(rng)
setGPExperimentInfo()
setGenerator(rampedHalfHalfGenerator, popSize = 10, generateOneByOne=false)
genPopulation!(GenJ)
gen1, gen2 = getGenotype(GenJ._population[1]), getGenotype(GenJ._population[6])
phen1 = getPhenotype(gen1)
phen2 = getPhenotype(gen2)
_,_, pointDepth = selectCrossPoint(gen1, rng)
selectCrossPoint(gen2, rng, minDepth=pointDepth, maxDepth=3)
#children = subtreeCross(gen1, gen2, GenJ._experimentInfo._GPExperimentInfo, rng)
#mutatedParent = subtreeMutation(parent1, GenJ._experimentInfo._GPExperimentInfo._functionSet, GenJ._experimentInfo._GPExperimentInfo._terminalSet, 3, GenJ._experimentInfo._rng, canonicalFullGenerator)
#fenotipo = getPhenotype(parent1, GenJ._experimentInfo._GPExperimentInfo)

f(x, y) = (y/5.5)+x^2

println("hola")
println("hola")
x = collect(-5:0.5:5)
y = collect(0:0.5:10)
nvalues = length(x)
objs = Array{Number}(undef, nvalues)
for i=1:nvalues
    objs[i] = f(x[i], y[i])
end

GeneticJulia.setIndividualType(GeneticJulia.CGPGenotype)
GeneticJulia.setRandomSeed(2)
GeneticJulia.setCGPInfo(nodesFile="GeneticJulia/src/utils/GeneticProgramming/Canonical/exampleNodesCGP.json")
GeneticJulia.setAlgorithm(GeneticJulia.SPEA, 10)
GeneticJulia.setStopCondition(maxIterations=50)
GeneticJulia.setEvaluator([GeneticJulia.FitnessFunction(GeneticJulia.compareFunctions, objs, x, y, weight=-1), GeneticJulia.FitnessFunction(nNodes, weight=-0.2)],
            globalFitnessFunction=GeneticJulia.FitnessFunction(GeneticJulia.pareto), compareFunctionArgs="global")
GeneticJulia.setGenerator(GeneticJulia.rampedHalfHalfGenerator, popSize = 50, generateOneByOne = false)
GeneticJulia.setSelector(GeneticJulia.tournamentSelector, 2)
GeneticJulia.setCrossoverOperator(GeneticJulia.subtreeCross, probability=0.9)
GeneticJulia.setMutationOperator(GeneticJulia.subtreeMutation, probability=0.2)
GeneticJulia.setReplacementOperator(GeneticJulia.replaceAllPopulation, eliteSize=5)
GeneticJulia.setExperimentSummary(batchSize=20, printDuringExperiment=true)

@time GeneticJulia.runGenJ(verbose=false)
"""

f(x, y) = (y/5.5)+x^2

println("hola")
println("hola")
x = collect(-5:0.5:5)
y = collect(0:0.5:10)
nvalues = length(x)
objs = Array{Number}(undef, nvalues)
for i=1:nvalues
    objs[i] = f(x[i], y[i])
end

GeneticJulia.setIndividualType(GeneticJulia.CGPGenotype)
GeneticJulia.setRandomSeed(54)
GeneticJulia.setCGPInfo(nodesFile="exampleNodesCGP.json")
GeneticJulia.setStopCondition(maxIterations=50)
GeneticJulia.setEvaluator([GeneticJulia.FitnessFunction(GeneticJulia.compareFunctions, objs, x, y, weight=-1)])
GeneticJulia.setGenerator(GeneticJulia.rampedHalfHalfGenerator, popSize = 50, generateOneByOne = false)
GeneticJulia.setSelector(GeneticJulia.tournamentSelector, 4)
GeneticJulia.setCrossoverOperator(GeneticJulia.subtreeCross, probability=0.9)
GeneticJulia.setMutationOperator(GeneticJulia.pointMutation, 0.2, probability=0.1)
GeneticJulia.setReplacementOperator(GeneticJulia.replaceAllPopulation, eliteSize=5)
GeneticJulia.setExperimentSummary(batchSize=-1)

@time GeneticJulia.runGenJ(verbose=false)

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
export TraceFlag,
    TraceState,
    SpanContext,
    Attributes,
    SPAN_KIND_CLIENT,
    SPAN_KIND_SERVER,
    SPAN_KIND_PRODUCER,
    SPAN_KIND_CONSUMER,
    SPAN_KIND_INTERNAL,
    Link,
    Event,
    SpanStatus,
    SPAN_STATUS_UNSET,
    SPAN_STATUS_ERROR,
    SPAN_STATUS_OK,
    Span,
    span_context,
    is_recording,
    add_event!,
    end!,
    set_status!,
    update_name!,
    record_exception!,
    NonRecordingSpan,
    global_tracer_provider,
    INVALID_SPAN_CONTEXT,
    current_span,
    create_span,
    with_span,
    get_tracer

using Dates: time

const TraceIdType = UInt128
const SpanIdType = UInt64

const SPAN_KEY = create_key("current-span")

Base.@kwdef struct TraceFlag
    sampled::Bool = false
end

const MAXIMUM_TRACESTATE_KEYS = 32
const TRACESTATE_KEY_PATTERN = r"^[a-z][_0-9a-z\\-\\*\\/]{0,255}|[a-z0-9][_0-9a-z\\-\\*\\/]{0,240}@[a-z][_0-9a-z\\-\\*\\/]{0,13}$"
const TRACESTATE_VAL_PATTERN = r"^[\x20-\x2b\x2d-\x3c\x3e-\x7e]{0,255}[\x21-\x2b\x2d-\x3c\x3e-\x7e]$"

function is_valid_trace_state_kv(k,v)
    if isnothing(match(TRACESTATE_KEY_PATTERN, k))
        @warn "Invalid key $k found"
        false
    elseif isnothing(match(TRACESTATE_VAL_PATTERN, v))
        @warn "Invalid val $v found"
        false
    else
        true
    end
end

"""
    TraceState(entries::Pair{String,String}...)

`TraceState` carries vendor-specific trace identification data, represented as a list of key-value pairs. `TraceState` allows multiple tracing systems to participate in the same trace. It is fully described in the [W3C Trace Context specification](https://www.w3.org/TR/trace-context/#tracestate-header).
"""
struct TraceState
    kv::NamedTuple
    function TraceState(entries::Pair{String,String}...)
        if length(entries) > MAXIMUM_TRACESTATE_KEYS
            @warn "There can't be more than $MAXIMUM_TRACESTATE_KEYS key/value pairs"
        else
            kv = NamedTuple(Symbol(k) => v for (k,v) in entries if is_valid_trace_state_kv(k,v))
            new(kv)
        end
    end
end

Base.getindex(s::TraceState, key::String) = s.kv[Symbol(key)]
Base.haskey(s::TraceState, key::String) = haskey(s.kv, Symbol(key))

"""
    SpanContext(;span_id, trace_id, is_remote, trace_flag=TraceFlag(), trace_state=TraceState())

A `SpanContext` represents the portion of a `Span` which must be serialized and
propagated along side of a distributed context. `SpanContext`s are immutable.

The OpenTelemetry `SpanContext` representation conforms to the [W3C TraceContext
specification](https://www.w3.org/TR/trace-context/). It contains two
identifiers - a `TraceId` and a `SpanId` - along with a set of common
`TraceFlags` and system-specific `TraceState` values.

`TraceId` A valid trace identifier is a 16-byte array with at least one
non-zero byte.

`SpanId` A valid span identifier is an 8-byte array with at least one non-zero
byte.

`TraceFlags` contain details about the trace. Unlike TraceState values,
TraceFlags are present in all traces. The current version of the specification
only supports a single flag called [sampled](https://www.w3.org/TR/trace-context/#sampled-flag).

`TraceState` carries vendor-specific trace identification data, represented as a list
of key-value pairs. TraceState allows multiple tracing
systems to participate in the same trace. It is fully described in the [W3C Trace Context
specification](https://www.w3.org/TR/trace-context/#tracestate-header).
"""
Base.@kwdef struct SpanContext
    span_id::SpanIdType
    trace_id::TraceIdType
    is_remote::Bool
    trace_flag::TraceFlag = TraceFlag()
    trace_state::TraceState = TraceState()
end

const INVALID_SPAN_ID = zero(SpanIdType)
const INVALID_TRACE_ID = zero(TraceIdType)

const INVALID_SPAN_CONTEXT = SpanContext(;
    trace_id = INVALID_TRACE_ID,
    span_id = INVALID_SPAN_ID,
    is_remote=false,
)

const TAttrVal = Union{
    String,
    Bool,
    Int64,
    Float64,
    Vector{String},
    Vector{Bool},
    Vector{Int64},
    Vector{Float64}
}

struct Attributes{T}
    count_limit::Int
    value_length_limit::Int
    kv::T
    n_dropped::Ref{Int}
end

_truncate(x, limit) = x

function _truncate(s::String, limit)
    if length(s) <= limit
        s
    else
        @warn "value is truncated to length $limit"
        s[1:min(end,limit)]
    end
end

function _truncate(xs::Vector{String}, limit)
    is_truncated = false
    for (i,s) in enumerate(xs)
        if length(s) > limit
            xs[i] = s[1:limit]
            is_truncated = true
        end
    end
    is_truncated && @warn "some elements in the value are truncated."
    xs
end

"""
    Attributes(kv::Pair{String,<:TAttrVal}...;count_limit=128, value_length_limit=typemax(Int), is_mutable=false)

Tha value type must be either `String`, `Bool`, `Int`, `Float64` or a `Vector` of above types. By default, we use a `NamedTuple` to represent the key-value pairs. If `is_mutable` is set to `true`, we'll use a `Dict` instead internally. If the value is of type `String` or `Vector{String}`, then each entry with length bigger than `value_length_limit` will be truncated.
"""
function Attributes(
    kv::Pair{String,<:TAttrVal}...
    ;count_limit=128,
    value_length_limit=typemax(Int),
    is_mutable=false
)
    if length(kv) > count_limit
        kv = kv[1:count_limit]
        n_dropped = Ref(length(kv) - count_limit)
    else
        n_dropped = Ref(0)
    end
    kv_truncated = (k=>_truncate(v, value_length_limit) for (k, v) in kv)
    if is_mutable
        data = Dict{String,TAttrVal}(kv)
    else
        data = NamedTuple(Symbol(k) => v for (k,v) in kv)
    end
    Attributes(count_limit, value_length_limit, data, n_dropped)
end

Base.getindex(d::Attributes, k::String) = getindex(d.kv, k)
Base.getindex(d::Attributes{<:NamedTuple}, k::String) = getindex(d.kv, Symbol(k))

function Base.setindex!(d::Attributes{<:Dict}, k::String, v::TAttrVal)
    if length(d.kv) == d.count_limit
        @warn "Count limit exceeded. Discarded."
        d.n_dropped[] += 1
    else
        d.kv[k] = _truncate(v, d.value_length_limit)
    end
end

@enum SpanKind begin
    SPAN_KIND_CLIENT
    SPAN_KIND_SERVER
    SPAN_KIND_PRODUCER
    SPAN_KIND_CONSUMER
    SPAN_KIND_INTERNAL
end

"""
    Link(span_context, attributes)

See more details at [links between spans](https://github.com/open-telemetry/opentelemetry-specification/blob/main/specification/overview.md#links-between-spans).
"""
struct Link
    context::SpanContext
    attributes::Attributes
end

"""
    Event(name, timestamp=time(), attributes=Attributes())
"""
Base.@kwdef struct Event
    name::String
    timestamp::Float64 = time()
    attributes::Attributes = Attributes()
end

@enum SpanStatusCode begin
    SPAN_STATUS_UNSET
    SPAN_STATUS_ERROR
    SPAN_STATUS_OK
end

"""
    SpanStatus(code, description=nothing)

Possible codes are:

- `SPAN_STATUS_UNSET`
- `SPAN_STATUS_ERROR`
- `SPAN_STATUS_OK`

`description` is required when `code` is `SPAN_STATUS_ERROR`.
"""
struct SpanStatus
    code::SpanStatusCode
    description::Union{String,Nothing}
    function SpanStatus(code::SpanStatusCode, description=nothing)
        if code === SPAN_STATUS_ERROR
            isnothing(description) && @error "description not provided"
            new(code, description)
        else
            # description is ignored
            new(code, nothing)
        end
    end
end

abstract type AbstractSpan end

mutable struct Span <: AbstractSpan
    name::String
    span_context::SpanContext
    parent::Union{Nothing,Span,SpanContext}
    kind::SpanKind
    start_time::Float64
    end_time::Union{Nothing,Float64}
    attributes::Attributes
    links::Vector{Link}
    events::Vector{Event}
    status::SpanStatus
end

"""
    Span(;kw...)

# Keyword Arguments

- `name`
- `span_context::SpanContext`
- `parent_context`
- `kind=SPAN_KIND_INTERNAL`
- `attributes=Attributes()`
- `links=[]`
- `start_time=time()`
"""
function Span(
    ;name,
    span_context,
    parent_context::Context=current_context(),
    kind=SPAN_KIND_INTERNAL,
    attributes=Attributes(;is_mutable=true),
    links=[],
    start_time=time()
)
    parent = get(parent_context, SPAN_KEY, nothing)
    Span(name, span_context, parent, kind, start_time, nothing, attributes, links, [], SpanStatus(SPAN_STATUS_UNSET))
end

is_end(s::Span) = !isnothing(s.end_time)

"""
    span_context(s::AbstractSpan)

Get the `SpanContext` from a span `s`.
"""
span_context(s::Span) = s.span_context

"""
    is_recording(s::AbstractSpan)

Returns `true` if this span `s` is recording information like events with the [`add_event!`](@ref) operation, attributes using [`setindex!`](@ref), status with [`set_status!`](@ref), etc.
"""
is_recording(s::Span) = !is_end(s)

"""
    s::AbstractSpan[key] = val

Set the attributes in span `s`. Only valid when the span is not ended yet.
"""
function Base.setindex!(s::Span, key, val)
    if is_end(s)
        @warn "the span is ended."
    else
        s.attributes[key] = val
    end
end

"""
    add_event!(s::AbstractSpan, event::Event)

Record a new `event`. Only valid when the span is not ended yet.
"""
function add_event!(s::Span, event::Event)
    if is_end(s)
        @warn "the span is ended."
    else
        push!(s.events, event)
    end
end

"""
    set_status!(s::Span, code::SpanStatusCode, description=nothing)

Update the status of span `s` by following the original specification. `description` is only considered when the `code` is `SPAN_STATUS_ERROR`. Only valid when the span is not ended yet.
"""
function set_status!(s::Span, code::SpanStatusCode, description=nothing)
    if is_end(s)
        @warn "the span is ended."
    else
        if s.status.code === OK
            # no further updates
        else
            if code === SPAN_STATUS_UNSET
                # ignore
            else
                s.status = StateStatus(code, description)
            end
        end
    end
end

"""
    update_name!(s::Span, name::String)

Update the name of span `s` with `name`. Only valid when the span is not ended yet.
"""
function update_name!(s::Span, name::String)
    if is_end(s)
        @warn "the span is ended."
    else
        s.name = name
    end
end

"""
    end!(s::Span, t=time())

Set the end time of the span.
"""
function end!(s::Span, t=time())
    if is_end(s)
        @warn "the span is ended"
    else
        s.end_time = t
    end
end

"""
    record_exception!(s::Span, ex::Exception, is_rethrow_followed=false)

A specialized variant of `add_event!` to record exceptions. Usually used in a `try... catch...end` to capture the backtrace. If the `ex` is `rethrow`ed in the `catch...end`, `is_rethrow_followed` should be set to `true`.
"""
function record_exception!(s::Span, ex::Exception, is_rethrow_followed=false)
    attrs = Attributes()

    attrs["exception.type"] = string(typeof(ex))

    msg_io = IOBuffer()
    showerror(msg_io, ex)
    attrs["exception.message"] = String(take!(msg_io))

    st_io = IOBuffer()
    showerror(st_io, CapturedException(ex, catch_backtrace()))
    attrs["exception.stacktrace"] = String(take!(st_io))
    attrs["exception.escaped"] = is_rethrow_followed

    add_event!(s, Event(name="exception", attributes=attrs))
end

struct NonRecordingSpan <: AbstractSpan
    span_context::SpanContext
end

for f in (:record_exception!, :end!, :update_name!, :set_status!, :add_event!)
    @eval $f(::NonRecordingSpan, args...) = nothing
end

Base.setindex!(::NonRecordingSpan, args...) = nothing
is_recording(::NonRecordingSpan) = false
span_context(s::NonRecordingSpan) = s.span_context

const INVALID_SPAN = NonRecordingSpan(INVALID_SPAN_CONTEXT)

"""
    current_span([current_context])

Get the span in the current context.
"""
current_span() = current_span(current_context())
current_span(ctx::Context) = get(ctx, SPAN_KEY, INVALID_SPAN)

#####
# Tracer
#####

abstract type AbstractTracer end

struct DefaultTracer <: AbstractTracer
end

"""
    create_span(name, t::AbstractTracer)

Create a new span based on tracer `t`.
"""
create_span(name, ::DefaultTracer) = INVALID_SPAN

"""
    with_span(f, s::AbstractSpan;kw...)

Call function `f` with the current span set to `s`.

# Keyword arguments

- `end_on_exit`, controls whether to call [`end!`](@ref) after `f` or not.
- `record_exception`, controls whether to record the exception.
- `set_status_on_exception`, decides whether to set status to [`SPAN_STATUS_ERROR`](@ref) automatically when an exception is caught.
"""
function with_span(
    f,
    s::AbstractSpan
    ;end_on_exit=true,
    record_exception=true,
    set_status_on_exception=true
)
    with_context(SPAN_KEY => s) do
        try
            f()
        catch ex
            if is_recording(s)
                if record_exception
                    record_exception!(s, ex, true)
                end
                if set_status_on_exception
                    set_status!(s, SPAN_STATUS_ERROR, string(ex))
                end
            end
            rethrow(ex)
        finally
            if end_on_exit
                end!(s)
            end
        end
    end
end

"""
    with_span(f, name, t::AbstractTracer; kw...)

The same with `with_span(f, create_span(name, t; kw...); kw...)`.
"""
function with_span(
    f,
    name::String,
    t::AbstractTracer
    ;end_on_exit=true,
    record_exception=true,
    set_status_on_exception=true,
    kw...
)
    with_span(
        f,
        create_span(name, t; kw...)
        ;end_on_exit=end_on_exit,
        record_exception=record_exception,
        set_status_on_exception=set_status_on_exception
    )
end

#####
# TracerProvider
#####

abstract type AbstractTracerProvider end

struct DefaultTracerProvider <: AbstractTracerProvider
end

"""
    get_tracer([global_tracer_provider], name="", version="", schema_url="")

Generate a new tracer based on tracer provider.
"""
get_tracer(args...) = get_tracer(global_tracer_provider(), args...)
get_tracer(::DefaultTracerProvider, args...) = DefaultTracer()

const GLOBAL_TRACER_PROVIDER = Ref{AbstractTracerProvider}(DefaultTracerProvider())

"""
    global_tracer_provider(p::AbstractTracerProvider)

Set the global tracer provider to `p`.
"""
global_tracer_provider(p::AbstractTracerProvider) = GLOBAL_TRACER_PROVIDER[] = p

"""
    global_tracer_provider()

Get the global tracer provider.
"""
global_tracer_provider() = GLOBAL_TRACER_PROVIDER[]


text:
suite["hamiltonian"] = BenchmarkGroup(["intracolgauge", "vertical_terms", "HIs", "field_terms", "left_terms", "right_terms", "davidson", "norm", "optimize"])
Ls   = buildLs(A, H; env_maxdim=χ, maxdim=D, cutoff=0.0)
Rs   = buildRs(A, H; env_maxdim=χ, maxdim=D, cutoff=0.0)
col  = 2
row  = 1

suite["hamiltonian"]["intracolgauge"] = @benchmarkable PEPS.intraColumnGauge(B, $col; maxdim=D, cutoff=0.0) setup=(B = copy($A))
A = PEPS.intraColumnGauge(A, col; maxdim=D, cutoff=0.0)
L = Ls[col - 1]
R = Rs[col + 1]
AncEnvs = PEPS.buildAncs(A, L, R, H, col)
ϕ = A[row, col]

field_H_terms  = PEPS.getDirectional(vcat(H[:, col]...), PEPS.Field)
vert_H_terms   = PEPS.getDirectional(vcat(H[:, col]...), PEPS.Vertical)
left_H_terms   = PEPS.getDirectional(vcat(H[:, col - 1]...), PEPS.Horizontal)
right_H_terms  = PEPS.getDirectional(vcat(H[:, col]...), PEPS.Horizontal)
suite["hamiltonian"]["norm"] = @benchmarkable PEPS.buildN($A, $L, $R, $AncEnvs[:I], $row, $col, $ϕ)
suite["hamiltonian"]["HIs"]  = @benchmarkable PEPS.buildHIs(B, L_, R_, $row, $col, φ) setup=(B = copy($A); φ=copy($ϕ); L_=copy(L); R_=copy(R))
suite["hamiltonian"]["vertical_terms"] = @benchmarkable PEPS.verticalTerms(B, L_, R_, $AncEnvs[:I], $AncEnvs[:V], $vert_H_terms, $row, $col, φ)  setup=(B = copy($A); φ=copy($ϕ); L_=copy(L); R_=copy(R))
suite["hamiltonian"]["field_terms"] = @benchmarkable PEPS.fieldTerms(B, L_, R_, $AncEnvs[:I], $AncEnvs[:F], $field_H_terms, $row, $col, φ)  setup=(B = copy($A); φ=copy($ϕ); L_=copy(L); R_=copy(R))
suite["hamiltonian"]["left_terms"] = @benchmarkable PEPS.connectLeftTerms(B, L_, R_, $AncEnvs[:I], $AncEnvs[:L], $left_H_terms, $row, $col, φ) setup=(B = copy($A); φ=copy($ϕ); L_=copy(L); R_=copy(R))
suite["hamiltonian"]["right_terms"] = @benchmarkable PEPS.connectRightTerms(B, L_, R_, $AncEnvs[:I], $AncEnvs[:R], $right_H_terms, $row, $col, φ) setup=(B = copy($A); φ=copy($ϕ); L_=copy(L); R_=copy(R))
mapper   = PEPS.ITensorMap(A, H, L, R, AncEnvs, row, col)
suite["hamiltonian"]["davidson"] = @benchmarkable PEPS.davidson($mapper, φ; maxiter=1, maxdim=D, cutoff=0.0) setup=(φ=copy($A[$row, $col]))
suite["hamiltonian"]["optimize"] = @benchmarkable PEPS.optimizeLocalH(B, $L, $R, AEs, $H, $row, $col; maxiter=1, maxdim=D, cutoff=0.0) setup=(B = copy($A); AEs=PEPS.buildAncs($A, $L, $R, $H, $col))


text:
export JKFIT, RIFIT, Chonky, AbstractDFERI, AbstractERI, SparseERI

"""
    Fermi.Integrals.AbstractERI

Abstract type common to all ERI types
"""
abstract type AbstractERI end

"""
    Fermi.Integrals.AbstractDFERI

Abstract type common to all density-fitted ERI types

# Struct tree

**AbstractDFERI** <: AbstractERI
"""
abstract type AbstractDFERI <: AbstractERI end

"""
    Fermi.Integrals.JKFIT

Concrete type representing a density-fitted ERI using a JKFIT auxiliar basis set.

# Fields

| Name       | Type       | Description |
|:-----------|:-----------|:----------------------------------------------------|
|`basisset`  | `BasisSet` | BasisSet object for the auxiliar JK basis |

# Examples

The JKFIT structure is used to build a IntegralHelper
```julia
julia> jk = Fermi.Integrals.JKFIT()
julia> ints =
```

# Struct tree

JKFIT <: AbstractDFERI <: AbstractERI
"""
struct JKFIT <: AbstractDFERI
    basisset::BasisSet
end

function JKFIT(mol::Molecule = Molecule())

    auxjk = Options.get("jkfit")
    # If aux is auto, determine the aux basis from the basis
    if auxjk == "auto"
        basis = Options.get("basis")
        dunning_name = Regex("cc-pv.z")
        auxjk = occursin(dunning_name, basis) ? basis*"-jkfit" : "cc-pvqz-jkfit"
    end

    return JKFIT(mol, auxjk)
end

function JKFIT(mol::Molecule, basis::String)
    return JKFIT(BasisSet(basis, mol.atoms))
end

"""
    Fermi.Integrals.RIFIT

Concrete type representing a density-fitted ERI using a RIFIT auxiliar basis set.

# Fields

| Name       | Type       | Description |
|:-----------|:-----------|:----------------------------------------------------|
|`basisset`  | `BasisSet` | BasisSet object for the auxiliar RI basis |

# Struct tree

RIFIT <: AbstractDFERI <: AbstractERI
"""
struct RIFIT <: AbstractDFERI
    basisset::BasisSet
end

function RIFIT(mol::Molecule = Molecule())

    auxri = Options.get("rifit")
    # If aux is auto, determine the aux basis from the basis
    if auxri == "auto"
        basis = Options.get("basis")
        dunning_name = Regex("cc-pv.z")
        auxri = occursin(dunning_name, basis) ? basis*"-rifit" : "cc-pvqz-rifit"
    end

    return RIFIT(mol, auxri)
end

function RIFIT(mol::Molecule, basis::String)
    return RIFIT(BasisSet(basis, mol.atoms))
end

struct Chonky <:AbstractERI end
struct SparseERI <: AbstractERI end


text:
const initialObjValue = -876.0447622921811
const objGap = -0.005*initialObjValue

const riskLowerBound = 1E-04

user_initial_guess = Float64[9.39814141926694 9.421255515714048 11.756567771905383 9.03105985090215 6.495199822103873 7.893513952143964 9.6428376293843 7.796408894631191 8.335034648446358 8.657465991368657 7.115048695747427 8.452274343408492 6.998361092184627 8.27908968589544 6.908705226052372 8.54305207251165 11.315525149090394 7.998047905866601 8.465636933461479 9.127946402350714 9.026701690459463 10.832136091852085 8.3641860057966 8.864957681128347 10.235575417962949 5.529769780666923 9.85918456151803 7.360198460447348 6.729160844842408 11.570312295489012 9.901915376888278 8.91009103195549 8.336630478924612 10.290425603128611 9.17860620357558 10.490179712858469 7.663934401259311 10.399145154402284 9.020943332984954 10.197797434727285 6.559779487299241 8.780536452378188 10.64847123161823 7.4750622615785 11.659781874493973 8.858396105736306 9.191913534866927 4.23897010607676 8.679633187607054 10.06638469811706 9.053463054014724 6.236216730628251 8.836919445418262 8.491230259785949 9.561506252924456 9.188621248249477 10.822455031590241 8.9416089784855 8.818711541022592 7.325298334614627 7.441309070055175 7.524325212452145 8.308213434537889 10.291516371249418 7.446636388832638 12.29687453286503 11.537077437784527 6.96186877808385 8.775252752454156 9.933468634937586 6.324603154919467 7.352431151406648 9.920476881329401 7.765810545828663 8.816917292091745 9.2271680967317 6.26568492081681 9.70235798324651 7.02846053279833 8.503747188803647 8.263642814765703 8.300435395377038 10.025236219254408 8.512723556042687 9.935048912509274 9.799580788378226 8.105558642569312 9.67748408825378 8.084039385733206 7.534874912384693 5.196313787027405 10.58781970474596 9.00010108726941 10.069119041911852 7.955708909302817 7.780073920344252 9.998997145069973 10.529530777752436 8.727780905531258 8.708508004514204]
user_initial_guess = user_initial_guess'

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
using ScalarSolvers
using Test

using JuMP, NLPModels, NLPModelsJuMP
using State, Stopping
using OptimizationProblems
using Polynomials

for solver in scalar_solvers
    printstyled("we solve AMPGO02 using $solver \n", color = :green)
    nlp = MathProgNLPModel(AMPGO02());
    nlp_at_x = NLPAtX([3.7])
    nlp_stop = NLPStopping(nlp, Stopping.unconstrained, nlp_at_x)
    optimal, stop = eval(solver)(nlp, nlp_stop, verbose = true)
    @show optimal
    @show typeof(stop)
    status(stop) == :Optimal ? color_status = :light_green : color_status = :light_yellow
    printstyled("status = $(String(status(stop))) \n", color = color_status)
    @test status(stop) == :Optimal
end


text:
VERSION >= v"0.4.0-dev+6641" && __precompile__()

module TerminalUI

using Compat

import Base.Terminals: width, height,
    cmove, cmove_col, cmove_line_up, cmove_line_down,
    cmove_up, cmove_down, UnixTerminal, CSI

import Base: UnitRange, size
import Base.Multimedia: display

# Libuv signal support to support SIGWINCH
import Base: wait, close, _uv_hook_close

const _sizeof_uv_signal = Base.uv_sizeof_handle(Base.UV_SIGNAL)

function __init__()
    global uv_jl_signalcb = cfunction(uv_signalcb,Void,(Ptr{Void},Cint))
end

type SignalListener
    handle::Ptr{Void}
    c::Condition
    let _check = function(err)
            err < 0 ? throw(Base.UVError("signal",err)) : nothing
        end
        function SignalListener(signum)
            global uv_jl_signalcb
            this = new(Libc.malloc(_sizeof_uv_signal), Condition())
            Base.associate_julia_struct(this.handle, this)
            Base.preserve_handle(this)
            _check(ccall(:uv_signal_init,Cint,(Ptr{Void},Ptr{Void}),
                Base.eventloop(), this.handle))
            _check(ccall(:uv_signal_start,Cint,(Ptr{Void},Ptr{Void},Cint),
                this.handle, uv_jl_signalcb::Ptr{Void}, signum))
            this
        end
    end
end

function uv_signalcb(handle::Ptr{Void},signum)
    signal = @Base.handle_as handle SignalListener
    notify(signal.c)
    nothing
end

wait(s::SignalListener) = wait(s.c)
close(s::SignalListener) = if s.handle != C_NULL
    ccall(:jl_close_uv, Void, (Ptr{Void},), s.handle)
end

function _uv_hook_close(s::SignalListener)
    Base.unpreserve_handle(s)
    Base.disassociate_julia_struct(s)
    s.handle = C_NULL
    notify(s.c)
    nothing
end

# Things missing from Base.Terminals
cmove(t::UnixTerminal,row,col) = write(t.out_stream,"$(CSI)$(row);$(col)H")

using Reactive

abstract Widget
abstract Screen

function draw(::Screen,::Void)
end

function draw(::Void,::Any)
end

# subsreen

abstract Subscreen <: Screen
istracking(s::Subscreen) = false

immutable SimpleSubscreen <: Subscreen
    parent::Screen
    rows::UnitRange
    cols::UnitRange
end

immutable TrackingSubscreen <: Subscreen
    _s::SimpleSubscreen
    children::Vector{TrackingSubscreen}
    scroll::Bool
    topwidget::Widget
end
simple(s::TrackingSubscreen) = s._s
istracking(s::TrackingSubscreen) = true
topwidget(s::TrackingSubscreen) = s.topwidget
subscreens(s::TrackingSubscreen) = s.children

TrackingSubscreen(_s::SimpleSubscreen,scroll::Bool,w::Widget) =
    TrackingSubscreen(_s, Vector{Screen}(), scroll, w::Widget)

function add_child(s::TrackingSubscreen, ss::TrackingSubscreen)
    push!(s.children, ss)
end

width(s::SimpleSubscreen) = length(s.cols)
height(s::SimpleSubscreen) = length(s.rows)
size(s::SimpleSubscreen) = (height(s),width(s))

topscreen(s::SimpleSubscreen) = topscreen(s.parent)
topscreen(s) = s

for f in (:width,:height,:size,:topscreen)
    @eval ($f)(s::TrackingSubscreen) = ($f)(simple(s))
end

swrite(s::TrackingSubscreen,args...;kwargs...) = swrite(simple(s),args...; kwargs...)
swrite_image(s::TrackingSubscreen,args...;kwargs...) = swrite_image(simple(s),args...; kwargs...)
clear(s::TrackingSubscreen,args...) = clear(simple(s),args...)

function translate(s::SimpleSubscreen, rows, cols)
    rows = (first(s.rows)+first(rows)-1):min(last(s.rows), first(s.rows)+last(rows)-1)
    cols = (first(s.cols)+first(cols)-1):min(last(s.cols), first(s.cols)+last(cols)-1)
    (rows,cols)
end

function translate⁻¹(s::SimpleSubscreen, rows, cols)
    rows = (first(rows)-first(s.rows)+1):(last(rows)-last(s.rows)+1)
    cols = (first(cols)-first(s.cols)+1):(last(cols)-last(s.cols)+1)
    (rows,cols)
end

function translate⁻¹(s::SimpleSubscreen, p::Tuple{Int,Int})
    r = translate⁻¹(s,p[1]:p[1],p[2]:p[2])
    (first(r[1]),first(r[2]))
end

translate⁻¹(s::TrackingSubscreen, p::Tuple{Int,Int}) = translate⁻¹(simple(s),p)

function clear(s::SimpleSubscreen, rows=1:length(s.rows), cols=1:length(s.cols))
    clear(s.parent,translate(s,rows,cols)...)
end

function swrite(s::SimpleSubscreen, rows, cols, char; args...)
    swrite(s.parent,translate(s,rows,cols)...,char; args...)
end

function swrite_image(s::SimpleSubscreen, rows, cols, img; args...)
    swrite_image(s.parent,translate(s,rows,cols)..., img; args...)
end

function subscreen(s::Screen, rows, cols, widget = nothing; scroll = false)
    if !(first(rows) >= 1 && last(rows) <= height(s))
        error("Tryed to create subscreen with rows $(rows) of a $(1:height(s)) screen")
    end
    if !(first(cols) >= 1 && last(cols) <= width(s))
        error("Tryed to create subscreen with columns $(cols) of a $(1:width(s)) screen")
    end
    ss = SimpleSubscreen(s, rows, cols)
    if scroll && istracking(s)
        @assert widget !== nothing
        ss = TrackingSubscreen(ss, true, widget)
        add_child(s,ss)
    end
    ss
end

# Unbuffered write directly to the terminal
# Quite inefficient but useful for debugging
type DirectTerminalScreen <: Screen
    tty::Base.Terminals.TTYTerminal
    # Buffer the size to make sure it doesn't change midway through the draw operation
    size::Tuple{Int,Int}
end
width(s::DirectTerminalScreen) = s.size[2]
height(s::DirectTerminalScreen) = s.size[1]
size(s::DirectTerminalScreen) = s.size
resize!(s, size) = s.size = size

function swrite(s::DirectTerminalScreen, rows, cols, char)
    cmove(s.tty,1,1)
    cur_row = 1
    for r in rows
        cmove(s.tty,r,1)
        for c in cols
            cmove_col(s.tty,c)
            write(s.tty,char)
        end
    end
end

# Screensize widget

immutable Screensize <: Widget
end

function draw(s::Screen,::Screensize)
    screensize = size(topscreen(s))
    towrite = string(typeof(topscreen(s)),':',screensize[2],'x',screensize[1])
    swrite(s,div(height(s),2),div(width(s),2)-div(sizeof(towrite),2),towrite)
end

# Trees
using AbstractTrees
immutable Tree
    t::Tree
end

# Most terminals use $(CSI)2J to indicate that the current contents
# should be scrolled into the history buffer, even though this is not
# super evident from the spec.
function fresh_screen(s::DirectTerminalScreen)
    w = width(s)
    h = height(s)
    write(s.tty.out_stream,"$(CSI)H$(CSI)2J")
end

function clear_screen(s::DirectTerminalScreen)
    w = width(s)
    h = height(s)
    cmove(s.tty,h,w)
    write(s.tty.out_stream,"$(CSI)1J")
end

const SIGWINCH = 28

function monitor_resize(tty)
    last_size = Base.displaysize(tty)
    s = Signal(Tuple{Int,Int}, last_size)
    signal = TerminalUI.SignalListener(TerminalUI.SIGWINCH)
    @schedule while true
        # First check if the terminal size changed on us while
        # we were processing the request. This is possible if we
        # have to re-enter the event loop while drawing, in which
        # case the signal would have been delivered with nobody
        # waiting for it. If this condition is true, we are guaranteed
        # that any resize will at least ping the event loop again.
        if Base.displaysize(tty) == last_size
            wait(signal)
        end
        last_size = Base.displaysize(tty)
        push!(s,last_size)
    end
    s
end

import VT100: readdec

global curspos = (1,1)
global curspos_condition = Condition()

function parse_terminal_reply(stream)
    (c1,f1) = readdec(stream)
    if c1 == ';'
        (c2,f2) = readdec(stream)
        if c2 == 'R'
            global curspos
            curspos = (f1,f2)
            notify(curspos_condition)
        end
    else
        error("Unrecognized reply message")
    end
end

invalidated = Signal(Any, nothing)
invalidated2 = Signal(Any, nothing)

include("focus.jl")
include("render.jl")
include("widgets.jl")
include("dialog.jl")
include("inputs.jl")

end # module


text:
"""
    AppVeyor(; config_file::Union{AbstractString, Nothing}="") -> AppVeyor

Add `AppVeyor` to a template's plugins to add a `.appveyor.yml` configuration file to
generated repositories, and an appropriate badge to the README.

# Keyword Arguments
* `config_file::Union{AbstractString, Nothing}=""`: Path to a custom `.appveyor.yml`.
  If `nothing` is supplied, no file will be generated.
"""
struct AppVeyor <: GenericPlugin
    gitignore::Vector{String}
    src::Union{String, Nothing}
    dest::String
    badges::Vector{Badge}
    view::Dict{String, Any}

    function AppVeyor(; config_file::Union{AbstractString, Nothing}="")
        if config_file != nothing
            config_file = if isempty(config_file)
                config_file = joinpath(DEFAULTS_DIR, "appveyor.yml")
            elseif isfile(config_file)
                abspath(config_file)
            else
                throw(ArgumentError("File $(abspath(config_file)) does not exist"))
            end
        end
        new(
            [],
            config_file,
            ".appveyor.yml",
            [
                Badge(
                    "Build Status",
                    "https://ci.appveyor.com/api/projects/status/github/{{USER}}/{{PKGNAME}}.jl?svg=true",
                    "https://ci.appveyor.com/project/{{USER}}/{{PKGNAME}}-jl",
                )
            ],
            Dict{String, Any}(),
        )
    end
end


text:
# a few utility functions copied/adapted from Images.jl

# github.com/JuliaGraphics/ColorVectorSpace.jl#abs-and-abs2
_abs(c::CT) where CT<:Color = mapreducec(v->abs(float(v)), +, zero(eltype(CT)), c)
_abs(c) = mapreducec(v->abs(float(v)), +, 0, c)

"blur images `A` and `B` with `sigma` and then compute a difference"
function blurdiff(A::AbstractArray, B::AbstractArray, sigma)
    # make sure arrays have the same size
    if size(A) != size(B)
        newsize = map(max, size(A), size(B))
        if size(A) != newsize
            A = copyto!(zeros(eltype(A), newsize...), A)
        end
        if size(B) != newsize
            B = copyto!(zeros(eltype(B), newsize...), B)
        end
    end

    kern = KernelFactors.IIRGaussian(sigma)
    Af = imfilter(A, kern, NA())
    Bf = imfilter(B, kern, NA())
    d = sad(Af, Bf)
    diffscale = max(_abs(maxabsfinite(A)), _abs(maxabsfinite(B)))
    diffpct = d / (length(Af) * diffscale)

    diffpct
end


text:
module Day8
export day8

import ..Console

function readInput(inputfile="inputs/input08.txt")
    Console.parseInstructions(eachline(inputfile))
end

function day8(inputfile="inputs/input08.txt")
    instructions = readInput(inputfile)
    (day8Part1(instructions), day8Part2(instructions))
end

function day8Part1(instr)
    gotoLoop!(Console.Program(instr)).acc
end

function gotoLoop!(prog)
    instructionsRun = Set{Int}()
    while !(prog.ip in instructionsRun || Console.hasFinished(prog))
        push!(instructionsRun, prog.ip)
        Console.step!(prog)
    end
    prog
end

function day8Part2(instructions)
    for (i, instr) in enumerate(instructions)
        if instr.name == "nop"
            instructions[i] = Console.Instruction("jmp", instr.arg)
            prog = gotoLoop!(Console.Program(instructions))
            if Console.hasFinished(prog)
                return prog.acc
            end
            instructions[i] = instr
        elseif instr.name == "jmp"
            instructions[i] = Console.Instruction("nop", instr.arg)
            prog = gotoLoop!(Console.Program(instructions))
            if Console.hasFinished(prog)
                return prog.acc
            end
            instructions[i] = instr
        end
    end
end

end #module


text: largez06() = largez(6)


text:
list_ex_str = "acedgfb cdfbe gcdfa fbcad dab cefabd cdfgeb eafb cagedb ab | cdfeb fcadb cdfeb cdbaf"

list_test_str = split(
"be cfbegad cbdgef fgaecd cgeb fdcge agebfd fecdb fabcd edb | fdgacbe cefdb cefbgd gcbe
edbfga begcd cbg gc gcadebf fbgde acbgfd abcde gfcbed gfec | fcgedb cgb dgebacf gc
fgaebd cg bdaec gdafb agbcfd gdcbef bgcad gfac gcb cdgabef | cg cg fdcagb cbg
fbegcd cbd adcefb dageb afcb bc aefdc ecdab fgdeca fcdbega | efabcd cedba gadfec cb
aecbfdg fbg gf bafeg dbefa fcge gcbea fcaegb dgceab fcbdga | gecf egdcabf bgf bfgea
fgeab ca afcebg bdacfeg cfaedg gcfdb baec bfadeg bafgc acf | gebdcfa ecba ca fadegcb
dbcfg fgd bdegcaf fgec aegbdf ecdfab fbedc dacgb gdcebf gf | cefg dcbef fcge gbcadfe
bdfegc cbegaf gecbf dfcage bdacg ed bedf ced adcbefg gebcd | ed bcgafe cdgba cbgef
egadfb cdbfeg cegd fecab cgb gbdefca cg fgcdab egfdb bfceg | gbdfcae bgc cg cgb
gcafb gcf dcaebfg ecagb gf abcdeg gaef cafbge fdbac fegbdc | fgae cfgab fg bagce",
'\n'
)

list_str = let input = "day8_input.txt"
    isfile(input) ? readlines(input) : list_test_str
end

# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 1

encoding = [
    "abcef"
    "cf"
    "acdeg"
    "acdfg"
    "bcdf"
    "abdfg"
    "abdefg"
    "abcdefg"
    "abcdfg"
]

function part1(list)
    output = []
    for line in list
        u, v = map(split, split(line, " | "))
        ulen = length.(u)

        for uleni in ulen
            if count(ulen .== uleni) == 1
                for j in 1:4
                    if length(v[j]) == uleni
                        push!(output, v[j])
                    end
                end
            end
        end
    end
    return output, length(output)
end

testit("Test Part 1", part1(list_test_str)[2], ==(26))
testit("Challenge Part 1", part1(list_str)[2], ==(255))

# Part 2

function part2(list)
    digs = Dict(
        "abcefg" => 0,
        "cf" => 1,
        "acdeg" => 2,
        "acdfg" => 3,
        "bcdf" => 4,
        "abdfg" => 5,
        "abdefg" => 6,
        "acf" => 7,
        "abcdefg" => 8,
        "abcdfg" => 9
        )

    total = 0
    for line in list
        u, v = map(split, split(line, " | "))
        ulen = length.(u)

        d = Dict{Char, Char}()
        push!(d, first(reduce(intersect, Set.(u[broadcast(∈((3, 5)), ulen)]))) => 'a')
        push!(d, first(reduce(intersect, Set.(u[broadcast(∈((2, 6)), ulen)]))) => 'f')
        push!(d, first(setdiff(Set(u[ulen .== 2][1]), reduce(intersect, Set.(u[broadcast(∈((2, 6)), ulen)])))) => 'c')
        push!(d, first(reduce(intersect, Set.(u[broadcast(∈((4, 5)), ulen)]))) => 'd')
        push!(d, first(setdiff(reduce(symdiff, Set.(u[broadcast(∈((3, 4)), ulen)])), keys(d))) => 'b')
        push!(d, first(setdiff(reduce(intersect, Set.(u[broadcast(∈((6)), ulen)])), Set(keys(d)))) => 'g')
        push!(d, first(setdiff('a':'g', keys(d))) => 'e')

        total += sum([10^(4-j) * digs[prod(sort(collect(Set(d[k] for k in v[j]))))] for j in eachindex(v)])
    end
    return total
end

testit("Test Part 2", part2(list_test_str), ==(61229))
testit("Challenge Part 2", part2(list_str), ==(982158))


text:
#=
plot/setup
2019-06-24, Jeff Fessler, University of Michigan
=#

# module MIRT_plot_setup # an experiment; would need Reexport

export mirt_setup_plot

using Plots: plot

"""
`function mirt_setup_plot()`

MIRT tests use Plots (with default GR backend).
When running on a remote linux server over X11 I get error messages from GKS.
The following kludge (?) tries to overcome this problem.
"""
function mirt_setup_plot()

	try
		plot(ones(3))
	catch
		if Sys.islinux
			ENV["GKS_WSTYPE"] = "x11"
		else
			@warn "Plot not working for Sys.KERNEL=$(Sys.KERNEL); test may fail"
		end
	end

	try
		plot(ones(3))
	catch
		@warn "Plot still not working; test may fail"
	end

	nothing
end

# end


text:
using Flux, DiffEqFlux, DiffEqSensitivity, OrdinaryDiffEq, LinearAlgebra, Zygote, Test

GDP = [11394358246872.6, 11886411296037.9, 12547852149499.6, 13201781525927, 14081902622923.3, 14866223429278.3, 15728198883149.2, 16421593575529.9, 17437921118338, 18504710349537.1, 19191754995907.1, 20025063402734.2, 21171619915190.4, 22549236163304.4, 22999815176366.2, 23138798276196.2, 24359046058098.6, 25317009721600.9, 26301669369287.8, 27386035164588.8, 27907493159394.4, 28445139283067.1, 28565588996657.6, 29255060755937.6, 30574152048605.8, 31710451102539.4, 32786657119472.8, 34001004119223.5, 35570841010027.7, 36878317437617.5, 37952345258555.4, 38490918890678.7, 39171116855465.5, 39772082901255.8, 40969517920094.4, 42210614326789.4, 43638265675924.6, 45254805649669.6, 46411399944618.2, 47929948653387.3, 50036361141742.2, 51009550274808.6, 52127765360545.5, 53644090247696.9, 55995239099025.6, 58161311618934.2, 60681422072544.7, 63240595965946.1, 64413060738139.7, 63326658023605.9, 66036918504601.7, 68100669928597.9, 69811348331640.1, 71662400667935.7, 73698404958519.1, 75802901433146, 77752106717302.4, 80209237761564.8, 82643194654568.3]
function monomial(cGDP, parameters, t)
    α1, β1, nu1, nu2, δ, δ2 = parameters

    [α1 * ((cGDP[1]))^β1]
end
GDP0 = GDP[1]

tspan = (1.0, 59.0)
p = [474.8501513113645, 0.7036417845990167, 0.0, 1e-10, 1e-10, 1e-10]
u0 = [GDP0]
if false
    prob = ODEProblem(monomial,[GDP0],tspan,p)
else ## false crashes. that is when i am tracking the initial conditions
    prob = ODEProblem(monomial,u0,tspan,p)
end
function predict_rd() # Our 1-layer neural network
  Array(solve(prob,Tsit5(),p=p,saveat=1.0:1.0:59.0,reltol=1e-4,sensealg=TrackerAdjoint()))
end

function loss_rd() ##L2 norm biases the newer times unfairly
    ##Graph looks better if we minimize relative error squared
    c = 0.0
    a = predict_rd()
    d = 0.0
    for i=1:59
        c += (a[i][1]/GDP[i]-1)^2 ## L2 of relative error
    end
    c + 3 * d
end

data = Iterators.repeated((), 100)
opt = ADAM(0.01)

peek = function () #callback function to observe training
    #reduces training speed by a lot
    println("Loss: ",loss_rd())
end

peek()
Flux.train!(loss_rd, Flux.params(p,u0), data, opt, cb=peek)
peek()

@test loss_rd() < 0.2

function monomial(dcGDP, cGDP, parameters, t)
    α1, β1, nu1, nu2, δ, δ2 = parameters
    dcGDP[1] = α1 * ((cGDP[1]))^β1
end

GDP0 = GDP[1]

tspan = (1.0, 59.0)
p = [474.8501513113645, 0.7036417845990167, 0.0, 1e-10, 1e-10, 1e-10]
u0 = [GDP0]
if false
    prob = ODEProblem(monomial,[GDP0],tspan,p)
else ## false crashes. that is when i am tracking the initial conditions
    prob = ODEProblem(monomial,u0,tspan,p)
end
function predict_adjoint() # Our 1-layer neural network
  Array(solve(prob,Tsit5(),p=p,saveat=1.0,reltol=1e-4,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true))))
end

function loss_adjoint() ##L2 norm biases the newer times unfairly
    ##Graph looks better if we minimize relative error squared
    c = 0.0
    a = predict_adjoint()
    d = 0.0
    for i=1:59
        c += (a[i][1]/GDP[i]-1)^2 ## L2 of relative error
    end
    c + 3 * d
end

data = Iterators.repeated((), 100)
opt = ADAM(0.01)

peek = function () #callback function to observe training
    #reduces training speed by a lot
    println("Loss: ",loss_adjoint())
end

peek()
Flux.train!(loss_adjoint, Flux.params(p,u0), data, opt, cb=peek)
peek()
@test loss_adjoint() < 0.2


text:
using Test, AdvancedHMC

include("common.jl")

θ_init = randn(D)
h = Hamiltonian(UnitEuclideanMetric(D), logπ, ∂logπ∂θ)
r_init = ones(D)

@test AdvancedHMC.kinetic_energy(h, r_init, θ_init) == D / 2


text:
"""
    ObsView(data, [indices])

Used to represent a subset of some `data` of arbitrary type by
storing which observation-indices the subset spans. Furthermore,
subsequent subsettings are accumulated without needing to access
actual data.

The main purpose for the existence of `ObsView` is to delay
data access and movement until an actual batch of data (or single
observation) is needed for some computation. This is particularily
useful when the data is not located in memory, but on the hard
drive or some remote location. In such a scenario one wants to
load the required data only when needed.

Any data access is delayed until `getindex` is called,
and even `getindex` returns the result of
[`obsview`](@ref) which in general avoids data movement until
[`getobs`](@ref) is called.
If used as an iterator, the view will iterate over the dataset
once, effectively denoting an epoch. Each iteration will return a
lazy subset to the current observation.

# Arguments

- **`data`** : The object describing the dataset. Can be of any
    type as long as it implements [`getobs`](@ref) and
    [`numobs`](@ref) (see Details for more information).

- **`indices`** : Optional. The index or indices of the
    observation(s) in `data` that the subset should represent.
    Can be of type `Int` or some subtype of `AbstractVector`.

# Methods

- **`getindex`** : Returns the observation(s) of the given
    index/indices. No data is copied aside
    from the required indices.

- **`numobs`** : Returns the total number observations in the subset.

- **`getobs`** : Returns the underlying data that the
    `ObsView` represents at the given relative indices. Note
    that these indices are in "subset space", and in general will
    not directly correspond to the same indices in the underlying
    data set.

# Details

For `ObsView` to work on some data structure, the desired type
`MyType` must implement the following interface:

- `getobs(data::MyType, idx)` :
    Should return the observation(s) indexed by `idx`.
    In what form is up to the user.
    Note that `idx` can be of type `Int` or `AbstractVector`.

- `numobs(data::MyType)` :
    Should return the total number of observations in `data`

The following methods can also be provided and are optional:

- `getobs(data::MyType)` :
    By default this function is the identity function.
    If that is not the behaviour that you want for your type,
    you need to provide this method as well.

- `obsview(data::MyType, idx)` :
    If your custom type has its own kind of subset type, you can
    return it here. An example for such a case are `SubArray` for
    representing a subset of some `AbstractArray`.

- `getobs!(buffer, data::MyType, [idx])` :
    Inplace version of `getobs(data, idx)`. If this method
    is provided for `MyType`, then `eachobs` can preallocate a buffer that is then reused
    every iteration. Note: `buffer` should be equivalent to the
    return value of `getobs(::MyType, ...)`, since this is how
    `buffer` is preallocated by default.

# Examples

```julia
X, Y = MLUtils.load_iris()

# The iris set has 150 observations and 4 features
@assert size(X) == (4,150)

# Represents the 80 observations as a ObsView
v = ObsView(X, 21:100)
@assert numobs(v) == 80
@assert typeof(v) <: ObsView
# getobs indexes into v
@assert getobs(v, 1:10) == X[:, 21:30]

# Use `obsview` to avoid boxing into ObsView
# for types that provide a custom "subset", such as arrays.
# Here it instead creates a native SubArray.
v = obsview(X, 1:100)
@assert numobs(v) == 100
@assert typeof(v) <: SubArray

# Also works for tuples of arbitrary length
subset = obsview((X, Y), 1:100)
@assert numobs(subset) == 100
@assert typeof(subset) <: Tuple # tuple of SubArray

# Use as iterator
for x in ObsView(X)
    @assert typeof(x) <: SubArray{Float64,1}
end

# iterate over each individual labeled observation
for (x, y) in ObsView((X, Y))
    @assert typeof(x) <: SubArray{Float64,1}
    @assert typeof(y) <: String
end

# same but in random order
for (x, y) in ObsView(shuffleobs((X, Y)))
    @assert typeof(x) <: SubArray{Float64,1}
    @assert typeof(y) <: String
end

# Indexing: take first 10 observations
x, y = ObsView((X, Y))[1:10]
```

# See also

[`obsview`](@ref),  [`getobs`](@ref), [`numobs`](@ref),
[`splitobs`](@ref), [`shuffleobs`](@ref),
[`kfolds`](@ref).
"""
struct ObsView{Tdata, I<:Union{Int,AbstractVector}} <: AbstractDataContainer
    data::Tdata
    indices::I

    function ObsView(data::T, indices::I) where {T,I}
        1 <= minimum(indices) || throw(BoundsError(data, indices))
        maximum(indices) <= numobs(data) || throw(BoundsError(data, indices))
        new{T,I}(data, indices)
    end
end

ObsView(data) = ObsView(data, 1:numobs(data))

# # don't nest subsets
ObsView(subset::ObsView) = subset

function ObsView(subset::ObsView, indices::Union{Int,AbstractVector})
    ObsView(subset.data, subset.indices[indices])
end

function Base.show(io::IO, subset::ObsView)
    if get(io, :compact, false)
        print(io, "ObsView{", typeof(subset.data), "} with " , numobs(subset), " observations")
    else
        print(io, summary(subset), "\n ", numobs(subset), " observations")
    end
end

function Base.summary(subset::ObsView)
    io = IOBuffer()
    print(io, typeof(subset).name.name, "(")
    Base.showarg(io, subset.data, false)
    print(io, ", ")
    Base.showarg(io, subset.indices, false)
    print(io, ')')
    first(readlines(seek(io,0)))
end

# compare if both subsets cover the same observations of the same data
# we don't care how the indices are stored, just that they match
# in order and values
function Base.:(==)(s1::ObsView, s2::ObsView)
    s1.data == s2.data && s1.indices == s2.indices
end

Base.IteratorEltype(::Type{<:ObsView}) = Base.EltypeUnknown()

Base.getindex(subset::ObsView, idx) =
    obsview(subset.data, subset.indices[idx])

Base.length(subset::ObsView) = length(subset.indices)

getobs(subset::ObsView) = getobs(subset.data, subset.indices)
getobs(subset::ObsView, idx) = getobs(subset.data, subset.indices[idx])

getobs!(buffer, subset::ObsView) = getobs!(buffer, subset.data, subset.indices)
getobs!(buffer, subset::ObsView, idx) = getobs!(buffer, subset.data, subset.indices[idx])

Base.parent(x::ObsView) = x.data

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

"""
    obsview(data, [indices])

Returns a lazy view of the observations in `data` that
correspond to the given `indices`. No data will be copied except
of the indices. It is similar to constructing an [`ObsView`](@ref),
but returns a `SubArray` if the type of
`data` is `Array` or `SubArray`. Furthermore, this function may
be extended for custom types of `data` that also want to provide
their own subset-type.

In case `data` is a tuple, the constructor will be mapped
over its elements. That means that the constructor returns a
tuple of `ObsView` instead of a `ObsView` of tuples.

If instead you want to get the subset of observations
corresponding to the given `indices` in their native type, use
`getobs`.

See [`ObsView`](@ref) for more information.
"""
obsview(data, indices=1:numobs(data)) = ObsView(data, indices)

##### Arrays / SubArrays

obsview(A::SubArray) = A

function obsview(A::AbstractArray{T,N}, idx) where {T,N}
    I = ntuple(_ -> :, N-1)
    return view(A, I..., idx)
end

getobs(a::SubArray) = getobs(a.parent, last(a.indices))

##### Tuples / NamedTuples
function obsview(tup::Union{Tuple, NamedTuple}, indices)
    map(data -> obsview(data, indices), tup)
end


text:
# general array for testing input to SpecialFunctions.jl
x = 1e-6:0.01:100.0

### Tests for besselj0
j0_SpecialFunctions = SpecialFunctions.besselj0.(big.(x)) # array to be tested against computed in BigFloats
@assert j0_SpecialFunctions[1] isa BigFloat               # just double check the higher precision

j0_64 = besselj0.(Float64.(x))
j0_32 = besselj0.(Float32.(x))
j0_big = besselj0.(big.(x))

# make sure output types match input types
@test j0_64[1] isa Float64
@test j0_32[1] isa Float32
@test j0_big[1] isa BigFloat

# test against SpecialFunctions.jl
@test j0_32 ≈ j0_SpecialFunctions

# BigFloat precision only computed to 128 bits
@test isapprox(j0_big, j0_SpecialFunctions, atol=1.5e-34)

# NaN should return NaN
@test isnan(besselj0(NaN))

# zero should return one
@test isone(besselj0(zero(Float32)))
@test isone(besselj0(zero(Float64)))
@test isone(besselj0(zero(BigFloat)))

# test that Inf inputs go to zero
@test besselj0(Inf32) == zero(Float32)
@test besselj0(Inf64) == zero(Float64)

### Tests for besselj1
j1_SpecialFunctions = SpecialFunctions.besselj1.(big.(x)) # array to be tested against computed in BigFloats
@assert j1_SpecialFunctions[1] isa BigFloat               # just double check the higher precision

j1_64 = besselj1.(Float64.(x))
j1_32 = besselj1.(Float32.(x))

# make sure output types match input types
@test j1_64[1] isa Float64
@test j1_32[1] isa Float32

# test against SpecialFunctions.jl
@test j1_64 ≈ j1_SpecialFunctions
@test j1_32 ≈ j1_SpecialFunctions

# NaN should return NaN
@test isnan(besselj1(NaN))

# zero should return zero
@test iszero(besselj1(zero(Float32)))
@test iszero(besselj1(zero(Float64)))

# test that Inf inputs go to zero
@test besselj1(Inf32) == zero(Float32)
@test besselj1(Inf64) == zero(Float64)

## Tests for besselj
# note this is not complete just a simple test
# this needs work and removing for now

#@test besselj(3, 1.0) ≈ SpecialFunctions.besselj(3, 1.0)
#@test besselj(-5, 6.1) ≈ SpecialFunctions.besselj(-5, 6.1)


text:
module Pitchjx

export
    pitchjx

using EzXML
using DataFrames
using Dates

include(joinpath(dirname(@__FILE__), "extractor.jl"))

"""
Scrape MLBAM pitchfx data.

# How to use

`pitchjx("2018-10-20")`
"""
function pitchjx(start, fin=start)
    @info "Initialize: Start"
    result = DataFrame()
    date = Date(start)
    findate = Date(fin)
    @info "Initialize: Finish!"
    @info "Extract dataset: Start"
    while date <= findate
        df = extract(date)
        if size(result) == (0, 0)
            result = df
        else
            result = vcat(result, df)
        end
        date += Dates.Day(1)
    end
    @info "Extract dataset: Finish!"
    return result
end

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ghr_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ghr")
JLLWrappers.@generate_main_file("ghr", UUID("07c12ed4-43bc-5495-8a2a-d5838ef8d533"))
end  # module ghr_jll


text:
using CSTParser
using Test

import CSTParser: parse, remlineinfo!, span, flisp_parse, typof, kindof, valof

@testset "CSTParser" begin

    include("parser.jl")
    include("interface.jl")
    include("display.jl")
    CSTParser.check_base()

end


text:
using LanguageServer, SymbolServer
using LanguageServer: process
using LanguageServer.JSONRPC: Request

function start_server()
    server = LanguageServerInstance(IOBuffer(), IOBuffer(), true, "", "", Dict())
    @async run(server)

    initstr = """{"jsonrpc":"2.0","id":0,"method":"initialize","params":{"processId":0,"rootPath":"","rootUri":null,"capabilities":{"workspace":{"applyEdit":true,"workspaceEdit":{"documentChanges":true},"didChangeConfiguration":{"dynamicRegistration":true},"didChangeWatchedFiles":{"dynamicRegistration":true},"symbol":{"dynamicRegistration":true,"symbolKind":{"valueSet":[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]}},"executeCommand":{"dynamicRegistration":true},"configuration":true,"workspaceFolders":true},"textDocument":{"publishDiagnostics":{"relatedInformation":true},"synchronization":{"dynamicRegistration":true,"willSave":true,"willSaveWaitUntil":true,"didSave":true},"completion":{"dynamicRegistration":true,"contextSupport":true,"completionItem":{"snippetSupport":true,"commitCharactersSupport":true,"documentationFormat":["markdown","plaintext"],"deprecatedSupport":true},"completionItemKind":{"valueSet":[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]}},"hover":{"dynamicRegistration":true,"contentFormat":["markdown","plaintext"]},"signatureHelp":{"dynamicRegistration":true,"signatureInformation":{"documentationFormat":["markdown","plaintext"]}},"definition":{"dynamicRegistration":true},"references":{"dynamicRegistration":true},"documentHighlight":{"dynamicRegistration":true},"documentSymbol":{"dynamicRegistration":true,"symbolKind":{"valueSet":[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]}},"codeAction":{"dynamicRegistration":true,"codeActionLiteralSupport":{"codeActionKind":{"valueSet":["","quickfix","refactor","refactor.extract","refactor.inline","refactor.rewrite","source","source.organizeImports"]}}},"codeLens":{"dynamicRegistration":true},"formatting":{"dynamicRegistration":true},"rangeFormatting":{"dynamicRegistration":true},"onTypeFormatting":{"dynamicRegistration":true},"rename":{"dynamicRegistration":true},"documentLink":{"dynamicRegistration":true},"typeDefinition":{"dynamicRegistration":true},"implementation":{"dynamicRegistration":true},"colorProvider":{"dynamicRegistration":true}}},"trace":"off","workspaceFolders":[]}}"""
    process(parse(Request, LanguageServer.JSON.parse(initstr)), server)
    process(Request{Val{Symbol("initialized")},Any}(0, nothing), server)
    return server
end

function load_file(text, server)
    process(Request{Val{Symbol("textDocument/didOpen")},LanguageServer.DidOpenTextDocumentParams}(0, LanguageServer.DidOpenTextDocumentParams(LanguageServer.TextDocumentItem("none", "julia", 0, text))), server)
    doc = server.documents[LanguageServer.URI2("none")]
end
function hover(text, line, col, server)
    doc = load_file(text, server)
    r = Request{Val{Symbol("textDocument/hover")},LanguageServer.TextDocumentPositionParams}(0, LanguageServer.TextDocumentPositionParams(LanguageServer.TextDocumentIdentifier("none"), LanguageServer.Position(line, col)))
    process(r, server)
end

function completion(text, line, col, server)
    doc = load_file(text, server)
    r = Request{Val{Symbol("textDocument/completion")},LanguageServer.CompletionParams}(0, LanguageServer.CompletionParams(LanguageServer.TextDocumentIdentifier("none"), LanguageServer.Position(line, col), nothing))
    process(r, server)
end

load_file("using CSTParser", server)
completion("ra", 0, 2, server)


text:
"""
Graphically represents a MixedDistribution structure, use the keyword func = :cdf
to represent the CDF instead of PDF
"""
function RecipesBase.plot(md::MixedDistribution; func = :pdf, fillalpha = 0.4, xlow = -Inf64, xhigh = Inf64)
    if func == :cdf
        return plot_cdf(md, xlow, xhigh)
    end
    p = RecipesBase.plot()
    for (wi,di) in zip(md.cont_weights, md.cont_dists)
        wf = x -> wi * Dst.pdf(di,x)
        m = mean(di)
        s = std(di)
        xaxis = max(m-5.0*s,minimum(di), xlow):0.01:min(m+5.0*s,maximum(di), xhigh) |> collect
        RecipesBase.plot!(p, xaxis, wf.(xaxis), fillalpha = fillalpha, fillcolor = :match, fillrange = 0)
    end

    for (xi,pi) in zip(md.mass_points,md.mass_probs)
        RecipesBase.plot!(p, [xi, xi], [0.0, pi], line = :arrow, color = :blue)
    end
    return p
end

function plot_cdf(md::MixedDistribution, xlow, xhigh; shallow = false)
    xlow_trunc  = max(xlow,  mean(md)-5.0*std(md))
    xhigh_trunc = min(xhigh, mean(md)+5.0*std(md))
    xaxis = collect(minimum(md.mass_points):0.01:maximum(md.mass_points))
    for itr in (max(minimum(md),xlow_trunc):0.01:first(md.mass_points),last(md.mass_points):0.01:min(maximum(md),xhigh_trunc))
        append!(xaxis, itr)
    end
    sort!(xaxis)
    cdf = x -> Dst.cdf(md, x)
    return RecipesBase.plot(xaxis, cdf.(xaxis))
end


text:
# Adapted from https://github.com/tpapp/TransformVariables.jl/blob/master/src/special_arrays.jl

export CorrCholesky

"""
    CorrCholesky(n)
Cholesky factor of a correlation matrix of size `n`.
Transforms ``n(n-1)/2`` real numbers to an ``n×n`` lower-triangular matrix `L`, such that
`L*L'` is a correlation matrix (positive definite, with unit diagonal).

# Notes
If
- `z` is a vector of `n` IID standard normal variates,
- `σ` is an `n`-element vector of standard deviations,
- `C` is obtained from `CorrCholesky(n)`,
then `Diagonal(σ) * C.L * z` is a zero-centered multivariate normal variate with the standard deviations `σ` and
correlation matrix `C.L * C.U`.
"""
struct CorrCholesky{N} <: TV.VectorTransform
    n::N
end

TV.dimension(t::CorrCholesky{Int}) = TV.dimension(CorrCholeskyUpper(t.n))

TV.dimension(t::CorrCholesky{StaticInt{n}}) where {n} = TV.dimension(CorrCholeskyUpper(t.n))

# TODO: For now we just transpose the CorrCholeskyUpper result. We should
# consider whether it can help performance to implement this directly for the
# lower triangular case
function TV.transform_with(flag::TV.LogJacFlag, t::CorrCholesky{StaticInt{n}}, x::AbstractVector, index) where {n}
    @nospecialize t
    U, ℓ, index = TV.transform_with(flag, CorrCholeskyUpper(n), x, index)
    U = UpperTriangular(SizedMatrix{n,n}(U.data))
    return Cholesky(U,'U', 0), ℓ, index
end

function TV.transform_with(flag::TV.LogJacFlag, t::CorrCholesky, x::AbstractVector, index)
    n = t.n
    U, ℓ, index = TV.transform_with(flag, CorrCholeskyUpper(n), x, index)
    return Cholesky(U,'U', 0), ℓ, index
end

TV.inverse_eltype(t::CorrCholesky, x::AbstractMatrix) = TV.extended_eltype(x)

function TV.inverse_at!(x::AbstractVector, index, t::CorrCholesky, L::LowerTriangular)
    return TV.inverse_at!(x, index, CorrCholeskyUpper(t.n), L')
end

function TV.inverse_at!(x::AbstractVector, index, t::CorrCholesky, U::UpperTriangular)
    return TV.inverse_at!(x, index, CorrCholeskyUpper(t.n), U)
end

function TV.inverse_at!(x::AbstractVector, index, t::CorrCholesky, C::Cholesky)
    return TV.inverse_at!(x, index, t, C.UL)
end


text:
using LinearAlgebra, ForwardDiff, StaticArrays, Plots

# rod
struct Rod
      m # mass
      J # inertia
      l # length
      g # gravity

      nq # configuration dimension
      nλ # force-input dimension
end

# kinematics
function kinematics(model, q)
      y, z, θ = q
      l = model.l

      # (p1y, p1z, p2y, p2z)
      @SVector [y + 0.5 * l * sin(θ),
                z - 0.5 * l * cos(θ),
                y - 0.5 * l * sin(θ),
                z + 0.5 * l * cos(θ)]
end

# jacobian
function jacobian(model, q)
      y, z, θ = q
      l = model.l

      @SMatrix [1.0 0.0 0.5 * l * cos(θ);
                0.0 1.0 0.5 * l * sin(θ);
                1.0 0.0 -0.5 * l * cos(θ);
                0.0 1.0 -0.5 * l * sin(θ)]
end

# mass matrix
function mass_matrix(model)
      m = model.m
      J = model.J

      Diagonal(@SVector [m, m, J])
end

# gravity
function gravity(model, q)
      y, z, θ = q
      m = model.m
      l = model.l
      g = model.g

      @SVector [0.0, m * g, 0.0]
end

# dynamics
function dynamics(model, q1, q2, q3, λ, h)
      nq = model.nq
      SVector{nq}(mass_matrix(model) * (2.0 * q2 - q1 - q3) / h
            - h * gravity(model, q2)
            + h * jacobian(model, q3)' * λ)
end

# create model
m = 1.0                    # mass
l = 1.0                    # length
J = 1.0 / 12.0 * m * l^2.0 # inertia
g = 9.18                   # gravity
nq = 3                     # configuration dimension
nλ = 4                     # force input dimension
h = 0.01                    # time step

model = Rod(m, J, l, g, nq, nλ)
q = zeros(model.nq)
λ = zeros(model.nλ)

kinematics(model, q)
jacobian(model, q)
mass_matrix(model)
gravity(model, q)
dynamics(model, q, q, q, λ, h)

# 1-step simulate
include(joinpath(pwd(), "src/solvers/newton.jl"))
dynamics(x) = dynamics(model, q, q, x, λ, h)
dynamics(q)
q⁺ = newton(dynamics, q)
dynamics(x) = dynamics(model, q, q⁺, x, λ, h)
q⁺ = newton(dynamics, q⁺)

# pinned dynamics
function simulated_pinned_dynamics()
      # conditions
      q1 = @SVector [0.5 * model.l, 0.0, 0.5 * π]
      q2 = copy(q1)

      # initial force guess
      λ1 = zeros(2)

      # trajectories
      q_hist = [q1, q2]
      λ_hist = []

      for t = 1:100
            # pinned dynamics
            function pinned_dynamics(x)
                  nq = model.nq
                  nλ = model.nλ

                  q = view(x, 1:nq)
                  λ = [zeros(2); view(x, nq .+ (1:2))]

                  SVector{nq + 2}([dynamics(model, q1, q2, q, λ, h);
                                   kinematics(model, q)[3:4]])
            end

            # step
            x3 = newton(pinned_dynamics, [q2; λ1], tol_r = 1.0e-12)

            # cache results
            push!(q_hist, view(x3, 1:nq))
            push!(λ_hist, view(x3, nq .+ (1:2)))

            q2 = q_hist[end]
            q1 = q_hist[end-1]
            λ1 = λ_hist[end]
      end

      return q_hist, λ_hist
end

q_hist, λ_hist = simulated_pinned_dynamics()

plot(hcat(q_hist...)', labels = ["y" "z" "t"])
plot(hcat(λ_hist...)')

# pinned dynamics (double pendulum)
function simulated_pinned_dynamics()
      # conditions
      q1 = @SVector [0.5 * model.l, 1.5, 0.5 * π, model.l + 0.5 * model.l, 1.5, 0.5 * π]
      q2 = copy(q1)

      # initial force guess
      λ1 = zeros(4)

      # trajectories
      q_hist = [q1, q2]
      λ_hist = []

      for t = 1:500
            println("t = $t")
            # pinned dynamics
            function pinned_dynamics(x)
                  nq = model.nq
                  nλ = model.nλ

                  q = view(x, 1:2 * nq)
                  λ = view(x, 2 * nq .+ (1:4))

                  SVector{2 * nq + 4}([dynamics(model,
                                    view(q1, 1:nq),
                                    view(q2, 1:nq),
                                    view(q, 1:nq),
                                    [λ[3:4]; λ[1:2]],
                                    h);
                                   dynamics(model,
                                    view(q1, nq .+ (1:nq)),
                                    view(q2, nq .+ (1:nq)),
                                    view(q, nq .+ (1:nq)),
                                    [zeros(2); -1.0 * λ[3:4]],
                                    h);
                                   kinematics(model, view(q, 1:nq))[3:4] - [0.0; 1.5];
                                   kinematics(model, view(q, 1:nq))[1:2] -
                                    kinematics(model, view(q, nq .+ (1:nq)))[3:4]])
            end

            # step
            x3 = newton(pinned_dynamics, [q2; zeros(4)], tol_r = 1.0e-12)

            # cache results
            push!(q_hist, view(x3, 1:2 * nq))
            push!(λ_hist, view(x3, 2 * nq .+ (1:4)))

            q2 = q_hist[end]
            q1 = q_hist[end-1]
            λ1 = λ_hist[end]
      end

      return q_hist, λ_hist
end

q_hist, λ_hist = simulated_pinned_dynamics()

plot(hcat([q[1:nq] for q in q_hist]...)', labels = ["y1" "z1" "t1"])
plot(hcat([q[nq .+ (1:nq)] for q in q_hist]...)', labels = ["y2" "z2" "t2"])
plot(hcat(λ_hist...)')

# Visualize
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

function visualize!(vis, model, x;
        color=RGBA(0.0, 0.0, 0.0, 1.0),
        r = 0.1, Δt = 0.1)

    default_background!(vis)

    l1 = Cylinder(Point3f0(0.0, 0.0, 0.0), Point3f0(0.0, 0.0, model.l),
        convert(Float32, 0.025))
    setobject!(vis["l1"], l1, MeshPhongMaterial(color = color))
    l2 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l),
        convert(Float32, 0.025))
    setobject!(vis["l2"], l2, MeshPhongMaterial(color = color))

    setobject!(vis["elbow"], Sphere(Point3f0(0.0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = RGBA(0.0, 1.0, 0.0, 1.0)))
    setobject!(vis["ee"], Sphere(Point3f0(0.0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = RGBA(0.0, 1.0, 0.0, 1.0)))

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

    T = length(x)
    for t = 1:T
        k1 = kinematics(model, view(x[t], 1:model.nq))
        k2 = kinematics(model, view(x[t], model.nq .+ (1:model.nq)))

        MeshCat.atframe(anim,t) do
            settransform!(vis["l1"], cable_transform([k1[3]; 0.0; k1[4]], [k1[1]; 0.0; k1[2]]))
            settransform!(vis["l2"], cable_transform([k2[3]; 0.0; k2[4]], [k2[1]; 0.0; k2[2]]))

            settransform!(vis["elbow"], Translation([k1[1]; 0.0; k1[2]]))
            settransform!(vis["ee"], Translation([k2[1]; 0.0; k2[2]]))
        end
    end

    settransform!(vis["/Cameras/default"],
       compose(Translation(0.0 , 1.0 , -1.0), LinearMap(RotZ(pi / 2.0))))

    MeshCat.setanimation!(vis, anim)
end

visualize!(vis, model, q_hist, Δt = h)

# pinned dynamics (double pendulum w/ impact)
"""
	Newton
"""
function newton_slack(res::Function, x;
        idx = (1:0),
		tol_r = 1.0e-8, tol_d = 1.0e-6)
	y = copy(x)
	Δy = copy(x)

    r = res(y)

    iter = 0

    while norm(r, 2) > tol_r && iter < 25
        ∇r = ForwardDiff.jacobian(res, y)

		try
        	Δy = -1.0 * ∇r \ r
		catch
			@warn "implicit-function failure"
			return y
		end

        α = 1.0

		iter_ls = 0

        while any((y + α * Δy)[idx] .< 0.0)
            α *= 0.5
            iter_ls += 1
            iter_ls == 50 && (@warn "line search (slack) failed")
        end

        iter_ls = 0

        while α > 1.0e-8 && iter_ls < 10
            ŷ = y + α * Δy
            r̂ = res(ŷ)

            if norm(r̂) < norm(r)
                y = ŷ
                r = r̂
                break
            else
                α *= 0.5
				iter_ls += 1
            end

			iter_ls == 10 && (@warn "line search failed")
        end

        iter += 1
    end

	iter == 10 && (@warn "Newton failure")

    return y
end

function simulated_pinned_impact_dynamics()
      # conditions
      q1 = @SVector [0.5 * model.l, 1.5, 0.5 * π, model.l + 0.5 * model.l, 1.5, 0.5 * π]
      q2 = copy(q1)

      function pinned_dynamics_nom(x)
            nq = model.nq
            nλ = model.nλ

            q = view(x, 1:2 * nq)
            λ = view(x, 2 * nq .+ (1:4))

            SVector{2 * nq + 4}([dynamics(model,
                              view(q1, 1:nq),
                              view(q2, 1:nq),
                              view(q, 1:nq),
                              [λ[3:4]; λ[1:2]],
                              h);
                             dynamics(model,
                              view(q1, nq .+ (1:nq)),
                              view(q2, nq .+ (1:nq)),
                              view(q, nq .+ (1:nq)),
                              [zeros(2); -1.0 * λ[3:4]],
                              h);
                             kinematics(model, view(q, 1:nq))[3:4] - [0.0; 1.5];
                             kinematics(model, view(q, 1:nq))[1:2] -
                              kinematics(model, view(q, nq .+ (1:nq)))[3:4]])
      end

      x1 = newton(pinned_dynamics_nom, [q2; zeros(4)],
          tol_r = 1.0e-8)

      # initial force guess
      λ1 = [x1[2nq .+ (1:4)]; 0.0]
      # λ1 = zeros(5)

      # trajectories
      q_hist = [q1, q2]
      λ_hist = []

      for t = 1:1000
            println("t = $t")

            ρ = 1.0
            # pinned dynamics
            function pinned_dynamics(x)
                  nq = model.nq
                  nλ = model.nλ

                  q = view(x, 1:2 * nq)
                  λ = view(x, 2 * nq .+ (1:5))
                  s = x[2 * nq + 5 + 1]

                  SVector{2 * nq + 6}([dynamics(model,
                                    view(q1, 1:nq),
                                    view(q2, 1:nq),
                                    view(q, 1:nq),
                                    [λ[3:4]; λ[1:2]],
                                    h);
                                   dynamics(model,
                                    view(q1, nq .+ (1:nq)),
                                    view(q2, nq .+ (1:nq)),
                                    view(q, nq .+ (1:nq)),
                                    [0.0; λ[5]; -1.0 * λ[3:4]],
                                    h);
                                   kinematics(model, view(q, 1:nq))[3:4] - [0.0; 1.5];
                                   kinematics(model, view(q, 1:nq))[1:2] -
                                    kinematics(model, view(q, nq .+ (1:nq)))[3:4];
                                   s - kinematics(model, view(q, nq .+ (1:nq)))[2];
                                   1.0 / ρ + λ[5] * s])
            end

            # step
            x3 = newton_slack(pinned_dynamics, [q2; λ1[1:4]; 0.1; 0.1],
                idx = (2 * nq + 4 .+ (1:2)),
                tol_r = 1.0e-8)

            for i = 1:5
                println("   i = $i")

                ρ *= 10.0
                # pinned dynamics
                function _pinned_dynamics(x)
                      nq = model.nq
                      nλ = model.nλ

                      q = view(x, 1:2 * nq)
                      λ = view(x, 2 * nq .+ (1:5))
                      s = x[2 * nq + 6]

                      SVector{2 * nq + 6}([dynamics(model,
                                        view(q1, 1:nq),
                                        view(q2, 1:nq),
                                        view(q, 1:nq),
                                        [λ[3:4]; λ[1:2]],
                                        h);
                                       dynamics(model,
                                        view(q1, nq .+ (1:nq)),
                                        view(q2, nq .+ (1:nq)),
                                        view(q, nq .+ (1:nq)),
                                        [0.0; λ[5]; -1.0 * λ[3:4]],
                                        h);
                                       kinematics(model, view(q, 1:nq))[3:4] - [0.0; 1.5];
                                       kinematics(model, view(q, 1:nq))[1:2] -
                                        kinematics(model, view(q, nq .+ (1:nq)))[3:4];
                                       s - kinematics(model, view(q, nq .+ (1:nq)))[2];
                                       1.0 / ρ - λ[5] * s])
                end

                # step
                x3 = newton_slack(_pinned_dynamics, x3,
                    idx = (2 * nq + 4 .+ (1:2)),
                    tol_r = 1.0e-8)

                println("COMP: $(x3[2nq + 5] * x3[2nq + 6])")
                println("ρ: $ρ")
                println("error: $(norm(_pinned_dynamics(x3)))")
            end

            # cache results
            push!(q_hist, view(x3, 1:2 * nq))
            push!(λ_hist, view(x3, 2 * nq .+ (1:5)))

            q2 = q_hist[end]
            q1 = q_hist[end-1]
            λ1 = λ_hist[end]
      end

      return q_hist, λ_hist
end

@time q_hist, λ_hist = simulated_pinned_impact_dynamics()
vis = Visualizer()
render(vis)
visualize!(vis, model, q_hist, Δt = h)


text:
@testset "permutations" begin
    for N = 1:4
        ps = collect(permutations(1:N))
        perminversion = [findfirst(isequal(invperm(p)), ps) for p in ps]
        permcomposition = [findfirst(isequal(p[q]), ps) for p in ps, q in ps]

        for i = 1:length(ps), j = 1:length(ps)
            @test inv(Bennu.Permutation{N}(i)) ==
                      Bennu.Permutation{N}(perminversion[i])
            @test Bennu.Permutation{N}(i) ∘ Bennu.Permutation{N}(j) ==
                      Bennu.Permutation{N}(permcomposition[i, j])
        end
    end

    for js in ([1,2,3,4], [1,3,2,4], [2,1,4,3], [2,4,1,3], [3,1,4,2], [3,4,1,2],
           [4,2,3,1], [4,3,2,1])
        p = inv(Bennu.tuplesortpermutation(js...))
        @test js == Bennu.getpermutedindex([1 3; 2 4], p, 1:4)
    end

    is = collect(LinearIndices((1:3,1:2)))
    @test [1,2,3,4,5,6] == Bennu.getpermutedindex(is, Bennu.Permutation{4}(1), 1:6)
    @test [1,4,2,5,3,6] == Bennu.getpermutedindex(is, Bennu.Permutation{4}(3), 1:6)
    @test [3,2,1,6,5,4] == Bennu.getpermutedindex(is, Bennu.Permutation{4}(8), 1:6)
    @test [3,6,2,5,1,4] == Bennu.getpermutedindex(is, Bennu.Permutation{4}(11), 1:6)
    @test [4,1,5,2,6,3] == Bennu.getpermutedindex(is, Bennu.Permutation{4}(14), 1:6)
    @test [4,5,6,1,2,3] == Bennu.getpermutedindex(is, Bennu.Permutation{4}(17), 1:6)
    @test [6,3,5,2,4,1] == Bennu.getpermutedindex(is, Bennu.Permutation{4}(22), 1:6)
    @test [6,5,4,3,2,1] == Bennu.getpermutedindex(is, Bennu.Permutation{4}(24), 1:6)

    is = collect(LinearIndices((1:3,)))
    @test [1,2,3] == Bennu.getpermutedindex(is, Bennu.Permutation{2}(1), 1:3)
    @test [3,2,1] == Bennu.getpermutedindex(is, Bennu.Permutation{2}(2), 1:3)
end


text:
# common functions (initialization, destruction etsc)

function __init__()

  # initialize MPI first so Petsc won't finalize it
  if !MPI.Initialized()
    MPI.Init()
  end

  for i=1:3
    if have_petsc[i]
      PetscInitialize(petsc_type[i])
    end
  end

  # we want Petsc to return errors to us, rather than using its own
  # error handlers, so that we can catch error codes and throw exceptions
  # need to do this for all Petsc versions

  for i=1:3
    if have_petsc[i]
      libname = C.petsc_libs[i]
      #val = @eval(cglobal((:PetscIgnoreErrorHandler, C.$libname)))
      #C.PetscPushErrorHandler(C.petsc_type[i], val, C_NULL)
    end
  end

  # register atexit finalizers for the Petsc libraries
  atexit() do
    for i=1:3
      if have_petsc[i]
        C.PetscFinalize(petsc_type[i])
      end
    end
  end
end

PetscInitialize(lib::DataType) = PetscInitialize(lib, [""], "", "")
PetscInitialize(lib::DataType, args) = PetscInitialize(lib, args, "", "")

function PetscInitialize(lib::DataType, args, filename, help)
  # argument list starts with program name
  args = filter(a->!isempty(a), args)
  args = map(string, vcat("julia", args))

  nargs = Int32[length(args)]

  # args are passed in as an Ptr{Ptr{Ptr{UInt8}}}
  args_ptr = map(pointer, args)
  args_ptr_ptr = [pointer(args_ptr)]
  args_ptr_ptr_ptr = pointer(args_ptr_ptr)

  err = C.PetscInitialize(lib, nargs, args_ptr_ptr_ptr, String(filename), String(help))
  if err != 0
    error("PETSc initialization error, code: $err")
  end
  return err
end

function petsc_sizeof(T::C.PetscDataType)
  sz = Array{Csize_t,1}(undef,1)
  chk(C.PetscDataTypeGetSize(C.petsc_type[1], T, sz))
  sz[1]
end

function PetscFinalized(T::Type)
  ret = Array{PetscBool,1}(undef,1)
  C.PetscFinalized(T, ret)
  ret[1] != 0
end