text:
function test_block_structure(t, o, l, u)
    B = FEDVR(t, o)
    M = Matrix(undef, B)

    @test M isa BlockSkylineMatrix
    bs = M.block_sizes
    @test length(bs.l) == length(l)+1
    @test length(bs.u) == length(u)+1

    # TODO Test block sizes

    @test all(bs.l[1:end-2] .== l[1:end-1])
    length(l) > 0 && @test bs.l[end-1] ≥ l[end]
    @test all(bs.u[3:end] .== u[2:end])
    length(u) > 0 && @test bs.u[2] ≥ u[1]
end

@testset "Block structure" begin
    @testset "Simple" begin
    test_block_structure(1.0:2, [4], [], [])
    test_block_structure(1.0:7, [2,2,3,4,2,4], [1,1,2,1,2,1,1,1], [1,1,1,2,1,2,1,1])
    test_block_structure(1.0:8, [2,2,3,4,2,4,2], [1,1,2,1,2,1,1,2,1,1], [1,1,1,2,1,2,1,1,2,1])
    test_block_structure(range(0,1,length=7), 1 .+ (1:6), [1,2,1,2,1,2,1,2,1,1], [1,1,2,1,2,1,2,1,2,1])
    test_block_structure(range(0,1,length=7), reverse(1 .+ (1:6)), [1,2,1,2,1,2,1,2,1,1], [1,1,2,1,2,1,2,1,2,1])
    test_block_structure(range(0,1,length=7), 2 .+ (1:6), [1,2,1,2,1,2,1,2,1,1], [1,1,2,1,2,1,2,1,2,1])
        @test Matrix(undef, FEDVR(1.0:7, 2)) isa Tridiagonal
    end
    @testset "Restrictions" begin
        L = 20.0
        n = 10
        t = range(-L/2, stop=L/2, length=n)
        B = FEDVR(t, [4,3,4,2,5,4,2,4,8])
        @test CompactBases.elements(B) == 1:9
        @test CompactBases.block_structure(B) == [3, 1, 1, 1, 2, 1, 1, 3, 1, 2, 1, 1, 2, 1, 7]
        # @test CompactBases.block_bandwidths(B) == ([1,2,1,2,1,1,2,1,2,1,1,2,1,1],
        #                                          [1,2,1,2,1,1,2,1,2,1,1,2,1,1])

        N = size(B,2)
        a = round(Int, 0.4N)
        b = round(Int, 0.6N)+1

        @test CompactBases.elements(B[:,1:b]) == 1:8
        @test CompactBases.block_structure(B[:,1:b]) == [3, 1, 1, 1, 2, 1, 1, 3, 1, 2, 1, 1]
        @test CompactBases.elements(B[:,2:b]) == 1:8
        @test CompactBases.block_structure(B[:,2:b]) == [2, 1, 1, 1, 2, 1, 1, 3, 1, 2, 1, 1]
        @test CompactBases.elements(B[:,3:b]) == 1:8
        @test CompactBases.block_structure(B[:,3:b]) == [1, 1, 1, 1, 2, 1, 1, 3, 1, 2, 1, 1]
        @test CompactBases.elements(B[:,4:b]) == 1:8
        @test CompactBases.block_structure(B[:,4:b]) == [2, 1, 2, 1, 1, 3, 1, 2, 1, 1]

        @test CompactBases.elements(B[:,1:b-1]) == 1:7
        @test CompactBases.block_structure(B[:,1:b-1]) == [3, 1, 1, 1, 2, 1, 1, 3, 1, 3]
        @test CompactBases.elements(B[:,1:b-2]) == 1:6
        @test CompactBases.block_structure(B[:,1:b-2]) == [3, 1, 1, 1, 2, 1, 1, 3, 1, 2]
        @test CompactBases.elements(B[:,1:b-3]) == 1:6
        @test CompactBases.block_structure(B[:,1:b-3]) == [3, 1, 1, 1, 2, 1, 1, 3, 1, 1]
        @test CompactBases.elements(B[:,1:b-4]) == 1:6
        @test CompactBases.block_structure(B[:,1:b-4]) == [3, 1, 1, 1, 2, 1, 1, 4]

        @test CompactBases.elements(B[:,a:end]) == 5:9
        @test CompactBases.block_structure(B[:,a:end]) == [3, 1, 2, 1, 1, 2, 1, 7]
        @test CompactBases.elements(B[:,a:end-1]) == 5:9
        @test CompactBases.block_structure(B[:,a:end-1]) == [3, 1, 2, 1, 1, 2, 1, 6]
        @test CompactBases.elements(B[:,a:end-7]) == 5:9
        @test CompactBases.block_structure(B[:,a:end-7]) == [3, 1, 2, 1, 1, 3]
        @test CompactBases.elements(B[:,a:end-8]) == 5:8
        @test CompactBases.block_structure(B[:,a:end-8]) == [3, 1, 2, 1, 1, 2]

        # @test CompactBases.block_bandwidths(B[:,a:end-1]) == ([1,2,1,1,2,1,1],[1,1,2,1,1,2,1])
        # @test CompactBases.block_bandwidths(B[:,a:end-7]) == ([1,2,1,1,1],[1,1,2,1,1])

        @testset "Single element" begin
            B = FEDVR(range(-L/2, stop=L/2, length=2), 8)
            @test CompactBases.block_structure(B) == [8]
            @test CompactBases.block_structure(B[:,1:end-1]) == [7]
            @test CompactBases.block_structure(B[:,2:end]) == [7]
            @test CompactBases.block_structure(B[:,2:end-1]) == [6]
            @test CompactBases.block_bandwidths(B) == ([0],[0])
            @test CompactBases.block_bandwidths(B[:,2:end-1]) == ([0],[0])
        end
    end
end

function test_blocks(f::Function, t, o)
    B = FEDVR(t, o)
    A = Matrix(undef, B)

    blocks = map(f, 1:CompactBases.nel(B))
    coords = vcat(1,1 .+ cumsum(B.order[1:end-1] .- 1))
    CompactBases.set_elements!(f, A, B)

    for (i,(b,c)) in enumerate(zip(blocks,coords))
        sel = c .+ (0:size(b,1)-1)
        b′ = copy(A[sel,sel])
        i > 1 && (b′[1,1] -= blocks[i-1][end,end])
        i < length(blocks) && (b′[end,end] -= blocks[i+1][1,1])
        @test all(b′ .== b)
    end
end

@testset "Set blocks" begin
    test_blocks(1.0:3,[2,3]) do i
        i*ones(i+1,i+1)
    end
    begin
        o = [2,2,3,4,2,4]
        test_blocks(1.0:7,o) do i
            i*ones(o[i],o[i])
        end
    end
    begin
        o = 4 .+ (1:6)
        test_blocks(range(0,1,length=7), o) do i
            i*ones(o[i],o[i])
        end
    end
end

@testset "Set blocks 2" begin
    L = 20.0
    n = 10
    t = range(-L/2, stop=L/2, length=n)
    B = FEDVR(t, [4,3,4,2,5,4,2,4,8])

    N = size(B,2)

    Bl = B[:,1:round(Int, 0.6N)+1]
    M = Matrix(undef, Bl, Int)
    @test M isa BlockSkylineMatrix
    @test eltype(M) == Int
    @test size(M) == (18,18)

    CompactBases.set_elements!(M, Bl) do i
        o = CompactBases.order(Bl, i)
        i*ones(Int,o,o)
    end

    for (v,ijs) in [1 => [(1,1), (1,2), (2,1)],
                    3 => [(2,2)],
                    2 => [(2,3), (2,4), (3,2), (4,2), (3,3), (3,4), (4,3)],
                    5 => [(4,4)],
                    3 => [(5,4), (6,4), (4,5), (4,6), (5,5), (6,5), (5,6)],
                    7 => [(6,6)],
                    4 => [(7,6), (6,7)],
                    9 => [(7,7)],
                    5 => [(8,7), (9,7), (7,8), (7,9), (8,8), (9,8), (8,9)],
                    11 => [(9,9)],
                    6 => [(10,9), (11,9), (9,10), (9,11), (10,10), (11,10), (10,11)],
                    13 => [(11,11)],
                    7 => [(12,11), (11,12)],
                    15 => [(12,12)]]
        for (i,j) in ijs
            @test all(M[Block(i,j)] .== v)
        end
    end

    Bl2 = B[:,4:round(Int, 0.6N)]
    M2 = Matrix(undef, Bl2)
    @time CompactBases.set_elements!(M2, Bl2) do i
        o = CompactBases.order(Bl2, i)
        i*ones(Int,o,o)
    end
    @test M2 == M[4:end-1,4:end-1]
end


text:
# linear-map/getindex.jl
#
# Provides getindex() capabilities like A[:,j] for LinearMap objects.
#
# Currently this provides only a partial set of the possible ways
# one can use indexing for a matrix, because currently
# LinearMap is not a subtype of an AbstractMatrix.
#
# If LinearMap were a subtype of an AbstractMatrix, then all possible
# indexing will be supported by Base.getindex, albeit very likely
# by quite inefficient iterators.
#
# These capabilities are provided as a user convenience.
# The user must recognize that LinearMap objects are not designed
# for efficient element-wise indexing.
#
# 2018-01-19, Jeff Fessler, University of Michigan

using LinearMaps
using Test

# A[end]
function Base.lastindex(A::LinearMap)
    return prod(size(A))
end

# A[?,end] and A[end,?]
function Base.lastindex(A::LinearMap, d::Integer)
    return size(A, d)
end

# A[i,j]
function Base.getindex(A::LinearMap, i::Integer, j::Integer)
	e = zeros(size(A,2)); e[j] = 1
	tmp = A * e
	return tmp[i]
end

# A[k]
function Base.getindex(A::LinearMap, k::Integer)
	c = CartesianIndices(size(A))[k] # is there a more elegant way?
	return A[c[1], c[2]]
#	return A[c[:]] # fails
#	return A[c...] # fails
#	return A[Tuple(c)...] # works
#	return getindex(A, c[1], c[2]) # works
end

# A[:,j]
# it is crucial to provide this function rather than to inherit from
# Base.getindex(A::AbstractArray, ::Colon, ::Integer)
# because Base.getindex does this by iterating (I think).
function Base.getindex(A::LinearMap, ::Colon, j::Integer)
	e = zeros(size(A,2)); e[j] = 1
	return A * e
end

# A[i,:]
function Base.getindex(A::LinearMap, i::Integer, ::Colon)
	# in Julia: A[i,:] = A'[:,i] for real matrix A else need conjugate
	return eltype(A) <: Complex ? conj.(A'[:,i]) : A'[:,i]
end

# A[:,j:k]
# this one is also important for efficiency
function Base.getindex(A::LinearMap, ::Colon, ur::UnitRange)
	return hcat([A[:,j] for j in ur]...)
end

# A[i:k,:]
Base.getindex(A::LinearMap, ur::UnitRange, ::Colon) = A'[:,ur]'

# A[:,:] = Matrix(A)
Base.getindex(A::LinearMap, ::Colon, ::Colon) = Matrix(A)

# A[i:k,j:l]
# this one is inefficient so could be inherited from Base?
Base.getindex(A::LinearMap, r1::UnitRange, r2::UnitRange) = A[:,r2][r1,:]

# A[???]
# informative error message in case we have overlooked any types
#=
function Base.getindex(A::LinearMap, kw...)
	@show kw
	for arg in kw
		@show typeof(arg)
	end
	throw("unsupported indexing type")
end
=#

# tests for getindex for LinearMaps

function ir_LinearMap_test_getindex(A::LinearMap)
	B = Matrix(A)
	@test all(size(A) .>= (4,4)) # required by tests
	@test B[1] == A[1]
	@test B[7] == A[7]
	@test B[:,5] == A[:,5]
	@test B[3,:] == A[3,:]
	@test B[1,3] == A[1,3]
	@test B[:,1:3] == A[:,1:3]
	@test B[1:3,:] == A[1:3,:]
	@test B[1:3,2:4] == A[1:3,2:4]
	@test B == A[:,:]
	@test B'[3] == A'[3]
	@test B'[:,4] == A'[:,4]
	@test B'[2,:] == A'[2,:]
	@test B[end] == A[end] # lastindex
	@test B[3,end-1] == A[3,end-1]

	# The following do not work because currently LinearMap is not a
	# subtype of AbstractMatrix.  If it were such a subtype, then LinearMap
	# would inherit general Base.getindex abilities
#=
	# todo later, if/when LinearMap is a subtype of AbstractMatrix
	@test B[[1, 3, 4]] == A[[1, 3, 4]]
	@test B[:, [1, 3, 4]] == A[:, [1, 3, 4]]
	@test B[[1, 3, 4], :] == A[[1, 3, 4], :]
	@test B[4:7] == A[4:7]
=#

	true
end

# tests for cumsum()
# note: the adjoint of cumsum is reverse(cumsum(reverse(y)))
function ir_lm_test_getindex_cumsum( ; N::Integer=5)
	A = LinearMap(cumsum, y -> reverse(cumsum(reverse(y))), N)
	ir_LinearMap_test_getindex(A)
	@test Matrix(A)' == Matrix(A')
#	@show A[:,4]
	true
end

ext: jl


text:
using Documenter, TraitSimulation

ENV["DOCUMENTER_DEBUG"] = "true"

makedocs(
    format = Documenter.HTML(),
    sitename = "TraitSimulation.jl",
    authors = "Sarah Ji, Chris German, Kenneth Lange, Janet Sinsheimer, Jin Zhou, Hua Zhou, Eric Sobel",
    clean = true,
    debug = true,
    pages = Any[
        "Home" => "index.md",
	"Testing Methods" => "examples/testing_MendelIHT_glm.md",
        "Real Data Examples" => Any[
            "examples/ukbiobank_vcm_power.md",
	    "examples/ukbiobank_ordered_multinomial_power.md"],
	"Exploring GPU Potential" => "examples/GPU.md",
     ]
)

deploydocs(
    repo   = "github.com/OpenMendel/TraitSimulation.jl.git",
    target = "build"
)


text:
module VTKExportMod

using FEM

using Zlib
using LightXML
using Codecs

import Base.push!
import FEM: write_data, get_coord, get_geotype, get_ncomponents, get_cell_data,
            get_displacement, get_field, get_galphas, DofVals
export VTKExporter, set_binary!, set_compress!

type VTKExporter <: AbstractDataExporter
    fields::Vector{DataType}
    binary::Bool
    compress::Bool
end
fields(vtkexp::VTKExporter) = vtkexp.fields
VTKExporter() = VTKExporter(Vector{DataType}[], true, false)

function set_binary!(vtkexp::VTKExporter, val::Bool)
    vtkexp.binary = val
    if !val
        vtkexp.compress = false
    end
end

function set_compress!(vtkexp::VTKExporter, val::Bool)
    if !vtkexp.binary && val
        warning("Can only compress when VTKEXporter is in binary mode")
        return
    else
        vtkexp.compress = val
    end
end
push!(vtkexp::VTKExporter, field::DataType) = push!(vtkexp.fields, field)

#Long ass names in spirit of the VTK library :)
abstract AbstractVTKXMLWriter

abstract AbstractVTKXMLBinaryWriter <: AbstractVTKXMLWriter

function add_data!(vtkw::AbstractVTKXMLBinaryWriter, data)
    write(vtkw.buffer, data)
end

type VTKXMLBinaryCompressedWriter <: AbstractVTKXMLBinaryWriter
    buffer::IOBuffer
end
VTKXMLBinaryCompressedWriter() = VTKXMLBinaryCompressedWriter(IOBuffer())

function write_data!(vtkw::VTKXMLBinaryCompressedWriter, xmlele::XMLElement)
    uncompressed_size = vtkw.buffer.size
    buff = takebuf_array(vtkw.buffer)
    #compressed_data = encode(Zlib, buff, 5)
    compressed_data = compress(buff, 5)
    compressed_size = length(compressed_data)
    header = UInt32[1, uncompressed_size, uncompressed_size, compressed_size]
    header_binary = bytestring(encode(Base64, reinterpret(UInt8, header)))
    data_binary = bytestring(encode(Base64, compressed_data))
    add_text(xmlele, header_binary)
    add_text(xmlele, data_binary)
end

type VTKXMLBinaryUncompressedWriter <: AbstractVTKXMLBinaryWriter
    buffer::IOBuffer
end
VTKXMLBinaryUncompressedWriter() = VTKXMLBinaryUncompressedWriter(IOBuffer())

function write_data!(vtkw::VTKXMLBinaryUncompressedWriter, xmlele::XMLElement)
    uncompressed_size = vtkw.buffer.size
    uncompressed_data = takebuf_array(vtkw.buffer)
    header = UInt32[uncompressed_size]
    header_binary = bytestring(encode(Base64, reinterpret(UInt8, header)))
    data_binary = bytestring(encode(Base64, uncompressed_data))
    add_text(xmlele, header_binary)
    add_text(xmlele, data_binary)
end

type VTKXMLASCIIWriter <: AbstractVTKXMLWriter
    buffer::IOBuffer
end
VTKXMLASCIIWriter() = VTKXMLASCIIWriter(IOBuffer())

function add_data!(vtkw::VTKXMLASCIIWriter, data)
    print(vtkw.buffer, data, " ")
end

function add_data!{T <: AbstractArray}(vtkw::VTKXMLASCIIWriter, data::T)
    for comp in data
        print(vtkw.buffer, comp, " ")
    end
end

function write_data!(vtkw::VTKXMLASCIIWriter, xmlele::XMLElement)
    add_text(xmlele, takebuf_string(vtkw.buffer))
end

get_vtk_num(::GeoTrig) = 5
get_vtk_num(::Type{GeoTrig}) = 5
get_vtk_num(::GeoQTrig) = 22
get_vtk_num(::Type{GeoQTrig}) = 22
get_vtk_num(::GeoQuad) = 9
get_vtk_num(::Type{GeoQuad}) = 9
get_vtk_num(::GeoTetra) = 10
get_vtk_num(::Type{GeoTetra}) = 10

function write_data(fp::FEProblem, vtkexp::VTKExporter, tstep::Int)

    if vtkexp.binary
        if vtkexp.compress
            writer = VTKXMLBinaryCompressedWriter()
        else
            writer = VTKXMLBinaryUncompressedWriter()
        end
    else
        writer = VTKXMLASCIIWriter()
    end

    for section in fp.sections
        _write_VTKXML_section(fp.name, fp.nodes, fp.dof_vals, section, vtkexp, writer, tstep)
    end
end

function _write_VTKXML_section(filename::ASCIIString, nodes::Vector{FENode2}, dof_vals::DofVals,
                               section::FESection, vtkexp:: VTKExporter, vtkw::AbstractVTKXMLWriter, tstep::Int)

    if vtkexp.binary
        const VTK_FORMAT = "binary"
    else
        const VTK_FORMAT = "ascii"
    end

    xdoc = XMLDocument()
    xroot = create_root(xdoc, "VTKFile")
    set_attribute(xroot, "type", "UnstructuredGrid")
    set_attribute(xroot, "version", "0.1")
    set_attribute(xroot, "byte_order", "LittleEndian")

    if vtkexp.compress
        set_attribute(xroot, "compressor", "vtkZLibDataCompressor")
    end

    xgrid = new_child(xroot, "UnstructuredGrid")

    xpiece = new_child(xgrid, "Piece")
    set_attribute(xpiece, "NumberOfPoints", length(nodes))

    ncells = length(section.elements)
    set_attribute(xpiece, "NumberOfCells", length(section.elements))

    # Points
    xpoints = new_child(xpiece, "Points")

    # Coordinates for points
    xcoords = new_child(xpoints, "DataArray")
    set_attribute(xcoords, "type", "Float64")
    set_attribute(xcoords, "name", "Points")
    set_attribute(xcoords, "format", VTK_FORMAT)
    set_attribute(xcoords, "NumberOfComponents", "3")

    for node in nodes
        coords = get_coord(node)
        for c in coords
            add_data!(vtkw, c)
        end
    end
    write_data!(vtkw, xcoords)

    # Cells
    xcells = new_child(xpiece, "Cells")

    # Cell connectivity
    xcellconn = new_child(xcells, "DataArray")
    set_attribute(xcellconn, "type", "Int64")
    set_attribute(xcellconn, "Name", "connectivity")
    set_attribute(xcellconn, "format", VTK_FORMAT)

    for element in section.elements
        for vertex in element.vertices
            add_data!(vtkw, vertex-1)
        end
    end
    write_data!(vtkw, xcellconn)

    # Cell location data
    xcell_offsets = new_child(xcells, "DataArray")
    set_attribute(xcell_offsets, "type", "Int64")
    set_attribute(xcell_offsets, "Name", "offsets")
    set_attribute(xcell_offsets, "format", VTK_FORMAT)
    nverts = length(section.elements[1].vertices)
    offsets = collect(nverts:nverts:nverts*length(section.elements))
    add_data!(vtkw, offsets)
    write_data!(vtkw, xcell_offsets)

    xcell_types = new_child(xcells, "DataArray")
    set_attribute(xcell_types, "type", "UInt8")
    set_attribute(xcell_types, "Name", "types")
    set_attribute(xcell_types, "format", VTK_FORMAT)
    ele_vtknum = get_vtk_num(get_geotype(section.elements[1]))
    cell_types = UInt8[ele_vtknum for _ in 1:length(section.elements)]
    add_data!(vtkw, cell_types)
    write_data!(vtkw, xcell_types)

    # Cell data
    xcell_data = new_child(xpiece, "CellData")
    for field in fields(vtkexp)
        # Separated this into its own function for type stability w.r.t. fields
        xcellfield = new_child(xcell_data, "DataArray")
        write_celldata_field!(xcellfield, field, section, vtkw, VTK_FORMAT)
    end

    # Point data
    xpoint_data = new_child(xpiece, "PointData")

    # Points
    xdisp = new_child(xpoint_data, "DataArray")
    set_attribute(xdisp, "Name", "Displacement")
    set_attribute(xdisp, "NumberOfComponents", "3")
    set_attribute(xdisp, "type", "Float64")
    set_attribute(xdisp, "format", VTK_FORMAT)
    for node in nodes
        disp = get_displacement(node, dof_vals)
        for i in disp
            add_data!(vtkw, i)
        end
    end
    write_data!(vtkw, xdisp)

    println("Finished writing data")

    save_file(xdoc, string(filename, "_$(tstep).vtu"))

    #free(xdoc) # Does this free everything?
end

function write_celldata_field!{T <: AbstractField}(xcellfield::XMLElement, field::Type{T},
                               section::FESection, vtkw::AbstractVTKXMLWriter, VTK_FORMAT::ASCIIString)
    set_attribute(xcellfield, "NumberOfComponents", get_ncomponents(field))
    set_attribute(xcellfield, "Name", string(field))
    set_attribute(xcellfield, "type", "Float64")
    set_attribute(xcellfield, "format", VTK_FORMAT)
    for elem in section.elements
        add_data!(vtkw, get_cell_data(elem, field))
    end
    write_data!(vtkw, xcellfield)
end

end # module


text:
import PredictMD
import Test

Test.@test( Base.VERSION >= VersionNumber("1.0") )

Test.@test( PredictMD.version() > VersionNumber(0) )

Test.@test(
    PredictMD.version() ==
        PredictMD.version(PredictMD)
    )

Test.@test(
    PredictMD.version() ==
        PredictMD.version(first(methods(PredictMD.eval)))
    )

Test.@test(
    PredictMD.version() ==
        PredictMD.version(PredictMD.eval)
    )

Test.@test(
    PredictMD.version() ==
        PredictMD.version(PredictMD.eval, (Any,))
    )

Test.@test( PredictMD.version(TestModuleA) == VersionNumber("1.2.3") )

Test.@test( PredictMD.version(TestModuleB) == VersionNumber("4.5.6") )

Test.@test_throws(
    ErrorException,
    PredictMD.TomlFile(joinpath(PredictMD.maketempdir(),"1","2","3","4")),
    )


text:
module Issue689

using Gridap
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.CellData
using Test

function main(transf,T)

  𝒯 = CartesianDiscreteModel((0,1,0,1),(10,10)) |> transf
  labels = get_face_labeling(𝒯)
  add_tag_from_tags!(labels,"surface",6)
  Ω = Interior(𝒯)
  Γ = Boundary(𝒯,tags="surface")

  k = 2
  dΩ = Measure(Ω,2*k)
  dΓ = Measure(Γ,2*k)

  reffe = ReferenceFE(lagrangian,T,k)
  V = FESpace(Ω,reffe)
  Q = FESpace(Γ,reffe)

  u((x,y)) = x^k + y^k
  uh = interpolate(u,V)
  qh = interpolate(u,Q)
  Δ_Γ(f) = x->∇∇(f)(x)[1,1]

  tol = 1.0e-9
  @test sqrt(sum( ∫( abs2( Δ(u) - Δ(uh) )  )dΩ )) < tol
  @test sqrt(sum( ∫( abs2( Δ(u) - Δ(uh) )  )dΓ )) < tol
  @test sqrt(sum( ∫( abs2( Δ_Γ(u) - Δ(qh) )  )dΓ )) < tol

end

main(simplexify,Float64)
main(identity,Float64)

end


text:
to_subscript(x::Union{Int, AbstractString}) = sprint(io->to_subscript(io, x))
to_subscript(io::IO, x::Char) = print(io, to_subscript(x))
function to_subscript(io::IO, x::AbstractString)
    for char in x
        print(io, to_subscript(char))
    end
end

to_superscript(x::Union{Int, AbstractString}) = sprint(io->to_superscript(io, x))
to_superscript(io::IO, x::Char) = print(io, to_superscript(x))
function to_superscript(io::IO, x::AbstractString)
    for char in x
        print(io, to_superscript(char))
    end
end

function to_subscript(io::IO, x::Int)
    if x in 0:9
        print(io, Char(0x2080+x))
    elseif x < 0
        print(io, Char(0x208B)); to_subscript(io, abs(x))
    else # positive + more than one digit
        for d in reverse(digits(x))
            to_subscript(io, d)
        end
    end
end

let subscript_map = Dict(
        '0' => '₀',
        '1' => '₁',
        '2' => '₂',
        '3' => '₃',
        '4' => '₄',
        '5' => '₅',
        '6' => '₆',
        '7' => '₇',
        '8' => '₈',
        '9' => '₉',
        '+' => '₊',
        '-' => '₋',
        '=' => '₌',
        '(' => '₍',
        ')' => '₎',
        'a' => 'ₐ',
        'e' => 'ₑ',
        'i' => 'ᵢ',
        'o' => 'ₒ',
        'r' => 'ᵣ',
        'u' => 'ᵤ',
        'v' => 'ᵥ',
        'x' => 'ₓ',
        'β' => 'ᵦ',
        'γ' => 'ᵧ',
        'ρ' => 'ᵨ',
        'φ' => 'ᵩ',
        'χ' => 'ᵪ',
    )
    global to_subscript
    function to_subscript(x::Char)
        haskey(subscript_map, x) || error("Char $x doesn't have a unicode superscript")
        subscript_map[x]
    end
end

function to_superscript(io::IO, x::Int)
    if x == 1
        print(io, Char(0x00B9))
    elseif x in 2:3
        print(io, Char(0x00B0+x))
    elseif x in 0:9
        print(io, Char(0x2070+x))
    elseif x < 0
        print(io, Char(0x207B)); to_superscript(io, abs(x))
    else # positive + more than one digit
        for d in reverse(digits(x))
            to_superscript(io, d)
        end
    end
end

let superscript_map = Dict(
        '.' => '⋅',
        '0' => '⁰',
        '1' => '¹',
        '2' => '²',
        '3' => '³',
        '4' => '⁴',
        '5' => '⁵',
        '6' => '⁶',
        '7' => '⁷',
        '8' => '⁸',
        '9' => '⁹',
        '+' => '⁺',
        '-' => '⁻',
        '=' => '⁼',
        '(' => '⁽',
        ')' => '⁾',
        'a' => 'ᵃ',
        'b' => 'ᵇ',
        'c' => 'ᶜ',
        'd' => 'ᵈ',
        'e' => 'ᵉ',
        'f' => 'ᶠ',
        'g' => 'ᵍ',
        'h' => 'ʰ',
        'i' => 'ⁱ',
        'j' => 'ʲ',
        'k' => 'ᵏ',
        'l' => 'ˡ',
        'm' => 'ᵐ',
        'n' => 'ⁿ',
        'o' => 'ᵒ',
        'p' => 'ᵖ',
        'r' => 'ʳ',
        's' => 'ˢ',
        't' => 'ᵗ',
        'u' => 'ᵘ',
        'v' => 'ᵛ',
        'w' => 'ʷ',
        'x' => 'ˣ',
        'y' => 'ʸ',
        'z' => 'ᶻ',
        'A' => 'ᴬ',
        'B' => 'ᴮ',
        'D' => 'ᴰ',
        'E' => 'ᴱ',
        'G' => 'ᴳ',
        'H' => 'ᴴ',
        'I' => 'ᴵ',
        'J' => 'ᴶ',
        'K' => 'ᴷ',
        'L' => 'ᴸ',
        'M' => 'ᴹ',
        'N' => 'ᴺ',
        'O' => 'ᴼ',
        'P' => 'ᴾ',
        'R' => 'ᴿ',
        'T' => 'ᵀ',
        'U' => 'ᵁ',
        'V' => 'ⱽ',
        'W' => 'ᵂ',
        'α' => 'ᵅ',
        'β' => 'ᵝ',
        'γ' => 'ᵞ',
        'δ' => 'ᵟ',
        '∊' => 'ᵋ',
        'θ' => 'ᶿ',
        'ι' => 'ᶥ',
        'Φ' => 'ᶲ',
        'φ' => 'ᵠ',
        'χ' => 'ᵡ'
    )
    global to_superscript
    function to_superscript(x::Char)
        haskey(superscript_map, x) || error("Char $x doesn't have a unicode superscript")
        superscript_map[x]
    end
end

"""
Turns given `numerator` and `denominator` into a fraction:
```
to_fraction("a-123", 392) -->
ᵃ⁻¹²³⁄₃₉₂
 ```
 Restricted to characters that can be turned into superscript and subscript.
 For a more general translation, see to_fraction_nl (newline)
"""
function to_fraction(numerator, denominator)
    sprint() do io
        to_fraction(io, numerator, denominator)
    end
end

function to_fraction(io::IO, numerator, denominator)
    to_superscript(io, numerator)
    print(io, Char(0x2044))
    to_subscript(io, denominator)
end

"""
Turns given `numerator` and `denominator` into a fraction with a newline:
```
to_fraction("α² ⋅ α²⁺³ ≡ α⁷", " ℝ: 𝐴𝐯 = λᵢ𝐯") -->

α̲²̲ ̲⋅̲ ̲α̲²̲⁺̲³̲ ̲≡̲ ̲α̲⁷̲
 ℝ: 𝐴𝐯 = λᵢ𝐯
 ```
"""
function to_fraction_nl(numerator, denominator)
    sprint() do io
        to_fraction_nl(io, numerator, denominator)
    end
end
function to_fraction_nl(io::IO, numerator, denominator)
    to_underline(io, numerator)
    println(io)
    print(io, denominator)
end


text:
## Alpha conversion

# α-conversion renames all bound variables to fresh ones
# This is necessary for transformations such as inlining

function αconvert(ex)
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
######################
# Functions
######################

# Basic functions
function func(x, y)
    x + y
end;

println(func(2+3im, 3)) # 5+3im

# Compact "assignment form"
f(x) = 5x + 3

println("$(f(0))") # 3

# Copying functions

function func2(x)
    x*(1+1im)
end

g = f # g(x) = 5x + 3

println("g(0) = $(g(0)), f(0) = $(f(0)), func2(0) = $(func2(0))")

f(x) = func2(x) # g = f = func2 = x*(1+1im)

println("g(0) = $(g(0)), f(0) = $(f(0)), func2(0) = $(func2(0))")

# Argument-type declarations
function exp(x::Integer, n::Integer) # exp(x, n) = x^n
    if n==1
        return x
    end

    temp = exp(x, div(n, 2))

    if n%2==1
        return temp*temp*x
    else
        return temp*temp
    end
end

println(exp(2, 5))

# Anonymous Functions
anonym = x -> (x-5)*x
println(anonym(10))

# Do-block syntax
# do keyword create a function that is passed as the first argument
function need_cb(callback::Function, x, y)
    callback(x^2 + 2*x*y + y^2)
end

need_cb(3, 5) do result
                println(result)
            end

# Function composition and piping
# data |> f1 |> f2 |> ...

function arr_to_string(arr)
    join(arr, " ")
end

function print_arr(arr)
    arr |> arr_to_string |> println
end

print_arr(1:10)

# combination of the chaining/piping and dot vectorization syntax
vec = [1, -4, 9, 16, -25]
vec .|> abs .|> sqrt .|> x->print(x," ")

max_stars_repo_name: Nikola-Mircic/julia-cheatsheet

max_issues_repo_name: Nikola-Mircic/julia-cheatsheet

max_forks_repo_name: Nikola-Mircic/julia-cheatsheet


text:
module ObjectDetector
export YOLO
export prepareImage, prepareImage!, resizekern, sizethatfits, emptybatch, drawBoxes

import Flux.gpu
import Flux.cpu
export gpu, cpu

using Pkg.Artifacts
getArtifact(name::String) = joinpath(@artifact_str(name), "$(name).weights")

using ImageFiltering
using ImageTransformations
using ImageCore

using BenchmarkTools
using PrettyTables
using ImageDraw

abstract type Model end
function getModelInputSize end

include("prepareimage.jl")

## YOLO models
include(joinpath(@__DIR__,"yolo","yolo.jl"))
import .YOLO

include("utils.jl")

end #module


text:
# In `solver_tests.jl` we want to load the preamble only once but if we load
# several time, say `csdp_tests.jl`, from the REPL it may be because we modify
# the tests in-between so we want to reload it
if !(@isdefined shared_preamble) || !shared_preamble
    include("Tests/Tests.jl")
    using Test, JuMP
end


text:
export resampleα, randomindex

"""
	resampleα(α::Float64, N::Int, K::Int)

Resample α given the approach by Escobar and West. (page 585)
"""
function resampleα(α::Float64, N::Int, K::Int; k = 1.0, Θ = 1.0, maxiter = 1)

    w = zeros(2)
    w[2] = N

    f0 = k + K
    f1 = Θ + K - 1

    newα = α

    for i in 1:maxiter

        η = rand(Beta(newα + 1, N))
        f3 = Θ - log(η)

        w[1] = f1 / f3

        z = randomindex(w)

        if z == 1
            newα = rand(Gamma(f0, 1/f3))
        else
            newα = rand(Gamma(f1, 1/f3))
        end
    end

    return newα
end

function resampleα(α::Float64, N::Array{Int}, K::Array{Int}; k = 1.0, Θ = 1.0, maxiter = 1)
    # implementation according to escobar - west page
    totalK = sum(K)
    num = length(N)

    η = zeros(num)

		newα = α
    for i in 1:maxiter

        for j in 1:num
            η[j] = rand(Beta(newα + 1, N[j]))
        end

        z = rand(num) .* (newα + N) .< N

        g_a = k + totalK - sum(z)
        g_b = Θ - sum(log(η))

        newα = rand(Gamma(g_a, 1/g_b))
    end

    newα
end

"""
	randomindex(p::Vector{Float64})

Randomly select an index propotional to its probability.
"""
function randomindex(p::Vector{Float64})

    max = sum(p)
    csum = 0.0
    thres = rand()

    @inbounds for i in 1:length(p)
        csum += p[i]
        if csum >= thres * max
            return convert(Int, i)
        end
    end

    return length(p)
end

function randnumtable(weights::Array{Float64}, table::Array{Int})

    numtable = zeros(Int, size(table))

    B = unique( table )
    w = log( weights )

    for i in 1:length(B)

        max = B[i]
        if max > 0

            m = 1:max

            stirnums = map( x -> abs(stirlings1(max, x)), m)
            stirnums /= maximum(stirnums)

            for (idx, j) in enumerate(findall(table .== max))

                clike = m .* w[idx]
                clike = cumsum(stirnums .* exp(clike - maximum(clike)))

                numtable[j] = 1+sum(rand() * clike[max] .> clike)
            end

        end

    end

    return numtable
end


text:
module BasePlus

include("os.jl")
include("text.jl")

end # module


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

@doc raw"""VolumeNodeResources is a set of resource limits for scheduling of volumes.

    IoK8sApiStorageV1VolumeNodeResources(;
        count=nothing,
    )

    - count::Int32 : Maximum number of unique volumes managed by the CSI driver that can be used on a node. A volume that is both attached and mounted on a node is considered to be used once, not twice. The same rule applies for a unique volume that is shared among multiple pods on the same node. If this field is not specified, then the supported number of volumes on this node is unbounded.
"""
mutable struct IoK8sApiStorageV1VolumeNodeResources <: SwaggerModel
    count::Any # spec type: Union{ Nothing, Int32 } # spec name: count

    function IoK8sApiStorageV1VolumeNodeResources(;count=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1VolumeNodeResources, Symbol("count"), count)
        setfield!(o, Symbol("count"), count)
        o
    end
end # type IoK8sApiStorageV1VolumeNodeResources

const _property_map_IoK8sApiStorageV1VolumeNodeResources = Dict{Symbol,Symbol}(Symbol("count")=>Symbol("count"))
const _property_types_IoK8sApiStorageV1VolumeNodeResources = Dict{Symbol,String}(Symbol("count")=>"Int32")
Base.propertynames(::Type{ IoK8sApiStorageV1VolumeNodeResources }) = collect(keys(_property_map_IoK8sApiStorageV1VolumeNodeResources))
Swagger.property_type(::Type{ IoK8sApiStorageV1VolumeNodeResources }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiStorageV1VolumeNodeResources[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1VolumeNodeResources }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1VolumeNodeResources[property_name]

function check_required(o::IoK8sApiStorageV1VolumeNodeResources)
    true
end

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


text:
module AccurateArithmetic

include("SIMDops.jl")

include("EFT.jl")
using .EFT
export two_sum, two_prod

include("Summation.jl")
using .Summation
export sum_naive, sum_kbn, sum_oro
export dot_naive, dot_oro

include("Test.jl")

end # module


text:
using Statistics
using PyCall
using JSON

datasets = pyimport("sklearn.datasets")

function get_data()
    url = "https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/"
    train_file = "a9a"
    test_file = train_file * ".t"
    !isfile(train_file) && download(url * train_file, train_file)
    !isfile(test_file) && download(url * test_file, test_file)
    X_train, y_train = datasets.load_svmlight_file(train_file)
    X_train = X_train.todense()
    X_test, y_test = datasets.load_svmlight_file(test_file)
    X_test = hcat(X_test.todense(), zeros(eltype(X_train), size(y_test)[1]))

    y_train = ((y_train .+ 1)/2)
    y_test = ((y_test .+ 1)/2)
    Array{Float32}(X_train), Array{Float32}(y_train), Array{Float32}(X_test), Array{Float32}(y_test)
end

function hmc(U, ∇U, ϵ, L, current_q)
    q = current_q
    p = randn(Float32, length(q))
    current_p = p
    # half step for momentum
    p = p - ϵ*∇U(q)/2
    @simd for i in 1:L
        # full step for the pos
        q = q + ϵ*p
        # full step for the momentum
        if i != L
            p = p - ϵ*∇U(q)
        end
    end
    # half step for momentum
    p = p - ϵ*∇U(q)/2
    # Negate momentum for symmetry
    p = -p
    current_U = U(current_q)
    current_K = (current_p'*current_p)/2
    proposed_U = U(q)
    proposed_K = (p'*p)/2
    # Accept or reject
    rand() < exp.(current_U .- proposed_U .+ current_K .- proposed_K)[1] ? q : current_q
end

function lr_hmc(X, y, ϵ, L, α, n_iter)
    U(β) = sum(log.(1 .+ exp.(X*β))) .- y'*(X*β) .+ (β'*β)/2α
    ∇U(β) = X'*(exp.(X*β)./(1 .+ exp.(X*β))-y).+ β/α
    D = size(X)[2]
    let q = zeros(Float32, size(X)[2], 1)
    hcat([q = hmc(U, ∇U, ϵ, L, q) for _ in 1:n_iter]...)'
    end
end

function benchmark()
    sigmoid(x) = 1 ./(1 .+ exp.(-x))
    X_train, y_train, X_test, y_test = get_data()
    params = JSON.parsefile("params.json")
    α = Float32(params["alpha"])
    n_iter = Int(params["n_iter"])
    ϵ = Float32(params["epsilon"])
    n_leaps = Int(params["n_leaps"])
    burn_in =  Int(params["burn_in"])
    lr_hmc(X_train, y_train, ϵ, n_leaps, α, 1)
    t = @elapsed z = lr_hmc(X_train, y_train, ϵ, n_leaps, α, n_iter)
    coef = mean(z[burn_in:end, :], dims=1)
    acc = mean(float(sigmoid((X_test*coef')) .> 0.5) .== y_test)
    @assert acc > 0.8
    t
end

println(JSON.json(Dict("julia" => benchmark())))

max_stars_repo_name: AndreyKolev/MLLibMark

max_issues_repo_name: AndreyKolev/MLLibMark

max_forks_repo_name: AndreyKolev/MLLibMark


text:
module LaneReversal

# package code goes here
using LightGraphs, Optim, BinDeps
using Distributed, Printf, LinearAlgebra, SparseArrays

include("load_network.jl")
include("frank_wolfe.jl")

export
        load_ta_network, download_tntp, read_ta_network,
        ta_frank_wolfe,
        TA_Data, solve_TAP

end # module

max_stars_repo_name: salomonw/tnet

max_issues_repo_name: salomonw/tnet

max_forks_repo_name: salomonw/tnet


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct PLL <: FrequencyEstimator
        ω_lp::Float64
        kp_pll::Float64
        ki_pll::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of a Phase-Locked Loop (PLL) for VSM

# Arguments
- `ω_lp::Float64`: PLL low-pass filter frequency (rad/sec), validation range: (0, nothing)
- `kp_pll::Float64`: PLL proportional gain, validation range: (0, nothing)
- `ki_pll::Float64`: PLL integral gain, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct PLL <: FrequencyEstimator
    "PLL low-pass filter frequency (rad/sec)"
    ω_lp::Float64
    "PLL proportional gain"
    kp_pll::Float64
    "PLL integral gain"
    ki_pll::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function PLL(ω_lp, kp_pll, ki_pll, ext=Dict{String, Any}(), )
    PLL(ω_lp, kp_pll, ki_pll, ext, [:vpll_d, :vpll_q, :ε_pll, :δθ_pll], 4, InfrastructureSystemsInternal(), )
end

function PLL(; ω_lp, kp_pll, ki_pll, ext=Dict{String, Any}(), )
    PLL(ω_lp, kp_pll, ki_pll, ext, )
end

# Constructor for demo purposes; non-functional.
function PLL(::Nothing)
    PLL(;
        ω_lp=0,
        kp_pll=0,
        ki_pll=0,
        ext=Dict{String, Any}(),
    )
end

"""Get PLL ω_lp."""
get_ω_lp(value::PLL) = value.ω_lp
"""Get PLL kp_pll."""
get_kp_pll(value::PLL) = value.kp_pll
"""Get PLL ki_pll."""
get_ki_pll(value::PLL) = value.ki_pll
"""Get PLL ext."""
get_ext(value::PLL) = value.ext
"""Get PLL states."""
get_states(value::PLL) = value.states
"""Get PLL n_states."""
get_n_states(value::PLL) = value.n_states
"""Get PLL internal."""
get_internal(value::PLL) = value.internal


text:
# Read all entries in PCREO_DATABASE_DIRECTORY and verify that they are
# syntactically valid PCREO files. This performs some basic sanity checks
# (e.g., floating-point numbers and integers can be correctly parsed, the
# advertised dimension and number of points in each configuration are correct)
# but does not recompute points configurations, energies, or nearest neighbors.

using ProgressMeter

push!(LOAD_PATH, @__DIR__)
using PCREO

function main()
    @assert isdir(ENV["PCREO_DATABASE_DIRECTORY"])
    for num_dir in lsdir(ENV["PCREO_DATABASE_DIRECTORY"]; join=true)
        println(num_dir)
        flush(stdout)
        @showprogress for entry_dir in lsdir(num_dir)
            for record_path in lsdir(num_dir, entry_dir; join=true)
                old = read(record_path, String)
                new = string(PCREORecord(record_path))
                if length(old) != length(new)
                    error(record_path)
                end
                for i = 1 : length(new)
                    if !(old[i] == new[i] ||
                         (old[i:i+3] == "-0.0" && new[i:i+3] == "+0.0"))
                        error(record_path)
                    end
                end
                write(record_path, new)
            end
        end
    end
end

main()


text:
using Base.Test
import Optics

# write your own tests here
@testset "SphericalLens Tests" begin
    println("Biconcave spherical lens")
    @show R1 = 2
    @show R2 = -2
    @show d = 1
    @show t = 1
    @show nl = 1.5

    (xR1, xR2, xV1, xV2) = Optics.sphericalLensVertexDiam(R1, R2, d)
    @test xR1 ≈ 2
    @test xR2 ≈ -1.872983346207417
    @test xV1 ≈ 0
    @test xV2 ≈ 0.12701665379258298

    (xe1, xe2) = Optics.sphericalLensxEdge(R1, R2, d)
    @test xe1 ≈ 0.06350832689629149
    @test xe2 ≈ 0.06350832689629149

    l = Optics.SphericalLens(R1, R2, d, nl)
    @test l.R1 ≈ R1
    @test l.R2 ≈ R2
    @test l.t ≈ Optics.sphericalLensMinThick(R1, R2, d)
    @test l.d ≈ d
    @test l.nl ≈ nl
    @test l.xV1 ≈ xV1
    @test l.xV2 ≈ xV2
    @test l.xR1 ≈ xR1
    @test l.xR2 ≈ xR2
    @test l.xe1 ≈ xe1
    @test l.xe2 ≈ xe2
    @test l.xmin ≈ xV1
    @test l.xmax ≈ xV2
    @test l.xmid ≈ (xV1 + xV2)/2
    @test l.yR ≈ 0
    @test l.ymin ≈ -d/2
    @test l.ymax ≈ d/2

    (xR1, xR2, xV1, xV2) = Optics.sphericalLensVertexThick(R1, R2, t)
    @test xR1 ≈ 2
    @test xR2 ≈ -1
    @test xV1 ≈ 0
    @test xV2 ≈ 1

    (xe1, xe2) = Optics.sphericalLensxEdge(R1, R2, d, t)
    @test xe1 ≈ 0.06350832689629149
    @test xe2 ≈ 0.9364916731037085

    l = Optics.SphericalLens(R1, R2, t, d, nl)
    @test l.R1 ≈ R1
    @test l.R2 ≈ R2
    @test l.t ≈ t
    @test l.d ≈ d
    @test l.nl ≈ nl
    @test l.xV1 ≈ xV1
    @test l.xV2 ≈ xV2
    @test l.xR1 ≈ xR1
    @test l.xR2 ≈ xR2
    @test l.xe1 ≈ xe1
    @test l.xe2 ≈ xe2
    @test l.xmin ≈ xV1
    @test l.xmax ≈ xV2
    @test l.xmid ≈ (xV1 + xV2)/2
    @test l.yR ≈ 0
    @test l.ymin ≈ -d/2
    @test l.ymax ≈ d/2

    println("R1 > 0, R2 > 0, spherical lens")
    @show R1 = 2
    @show R2 = 2
    @show d = 1
    @show t = 1
    @show nl = 1.5

    (xR1, xR2, xV1, xV2) = Optics.sphericalLensVertexDiam(R1, R2, d)
    @test xR1 ≈ 2.0
    @test xR2 ≈ 2.0
    @test xV1 ≈ 0.0
    @test xV2 ≈ 0.0

    (xe1, xe2) = Optics.sphericalLensxEdge(R1, R2, d)
    @test xe1 ≈ 0.06350832689629149
    @test xe2 ≈ 0.06350832689629149

    l = Optics.SphericalLens(R1, R2, d, nl)
    @test l.R1 ≈ R1
    @test l.R2 ≈ R2
    @test l.t ≈ Optics.sphericalLensMinThick(R1, R2, d)
    @test l.d ≈ d
    @test l.nl ≈ nl
    @test l.xV1 ≈ xV1
    @test l.xV2 ≈ xV2
    @test l.xR1 ≈ xR1
    @test l.xR2 ≈ xR2
    @test l.xe1 ≈ xe1
    @test l.xe2 ≈ xe2
    @test l.xmin ≈ xV1
    @test l.xmax ≈ xe2
    @test l.xmid ≈ (xV1 + xV2)/2
    @test l.yR ≈ 0
    @test l.ymin ≈ -d/2
    @test l.ymax ≈ d/2

    (xR1, xR2, xV1, xV2) = Optics.sphericalLensVertexThick(R1, R2, t)
    @test xR1 ≈ 2.0
    @test xR2 ≈ 3.0
    @test xV1 ≈ 0.0
    @test xV2 ≈ 1.0

    (xe1, xe2) = Optics.sphericalLensxEdge(R1, R2, d, t)
    @test xe1 ≈ 0.06350832689629149
    @test xe2 ≈ 1.06350832689629149

    l = Optics.SphericalLens(R1, R2, t, d, nl)
    @test l.R1 ≈ R1
    @test l.R2 ≈ R2
    @test l.t ≈ t
    @test l.d ≈ d
    @test l.nl ≈ nl
    @test l.xV1 ≈ xV1
    @test l.xV2 ≈ xV2
    @test l.xR1 ≈ xR1
    @test l.xR2 ≈ xR2
    @test l.xe1 ≈ xe1
    @test l.xe2 ≈ xe2
    @test l.xmin ≈ xV1
    @test l.xmax ≈ xe2
    @test l.xmid ≈ (xV1 + xV2)/2
    @test l.yR ≈ 0
    @test l.ymin ≈ -d/2
    @test l.ymax ≈ d/2

    println("R1 < 0, R2 > 0, spherical lens")
    @show R1 = -2
    @show R2 = 2
    @show d = 1
    @show t = 1
    @show nl = 1.5

    (xR1, xR2, xV1, xV2) = Optics.sphericalLensVertexDiam(R1, R2, d)
    @test xR1 ≈ -2.0
    @test xR2 ≈ 2.0
    @test xV1 ≈ 0.0
    @test xV2 ≈ 0.0

    (xe1, xe2) = Optics.sphericalLensxEdge(R1, R2, d)
    @test xe1 ≈ -0.06350832689629149
    @test xe2 ≈ 0.06350832689629149

    l = Optics.SphericalLens(R1, R2, d, nl)
    @test l.R1 ≈ R1
    @test l.R2 ≈ R2
    @test l.t ≈ Optics.sphericalLensMinThick(R1, R2, d)
    @test l.d ≈ d
    @test l.nl ≈ nl
    @test l.xV1 ≈ xV1
    @test l.xV2 ≈ xV2
    @test l.xR1 ≈ xR1
    @test l.xR2 ≈ xR2
    @test l.xe1 ≈ xe1
    @test l.xe2 ≈ xe2
    @test l.xmin ≈ xe1
    @test l.xmax ≈ xe2
    @test l.xmid ≈ (xV1 + xV2)/2
    @test l.yR ≈ 0
    @test l.ymin ≈ -d/2
    @test l.ymax ≈ d/2

    (xR1, xR2, xV1, xV2) = Optics.sphericalLensVertexThick(R1, R2, t)
    @test xR1 ≈ -2.0
    @test xR2 ≈ 3.0
    @test xV1 ≈ 0.0
    @test xV2 ≈ 1.0

    (xe1, xe2) = Optics.sphericalLensxEdge(R1, R2, d, t)
    @test xe1 ≈ -0.06350832689629149
    @test xe2 ≈ 1.06350832689629149

    l = Optics.SphericalLens(R1, R2, t, d, nl)
    @test l.R1 ≈ R1
    @test l.R2 ≈ R2
    @test l.t ≈ t
    @test l.d ≈ d
    @test l.nl ≈ nl
    @test l.xV1 ≈ xV1
    @test l.xV2 ≈ xV2
    @test l.xR1 ≈ xR1
    @test l.xR2 ≈ xR2
    @test l.xe1 ≈ xe1
    @test l.xe2 ≈ xe2
    @test l.xmin ≈ xe1
    @test l.xmax ≈ xe2
    @test l.xmid ≈ (xV1 + xV2)/2
    @test l.yR ≈ 0
    @test l.ymin ≈ -d/2
    @test l.ymax ≈ d/2

    println("R1 < 0, R2 < 0, spherical lens")
    @show R1 = -2
    @show R2 = -2
    @show d = 1
    @show t = 1
    @show nl = 1.5

    (xR1, xR2, xV1, xV2) = Optics.sphericalLensVertexDiam(R1, R2, d)
    @test xR1 ≈ -2.0
    @test xR2 ≈ -2.0
    @test xV1 ≈ 0.0
    @test xV2 ≈ 0.0

    (xe1, xe2) = Optics.sphericalLensxEdge(R1, R2, d)
    @test xe1 ≈ -0.06350832689629149
    @test xe2 ≈ -0.06350832689629149

    l = Optics.SphericalLens(R1, R2, d, nl)
    @test l.R1 ≈ R1
    @test l.R2 ≈ R2
    @test l.t ≈ Optics.sphericalLensMinThick(R1, R2, d)
    @test l.d ≈ d
    @test l.nl ≈ nl
    @test l.xV1 ≈ xV1
    @test l.xV2 ≈ xV2
    @test l.xR1 ≈ xR1
    @test l.xR2 ≈ xR2
    @test l.xe1 ≈ xe1
    @test l.xe2 ≈ xe2
    @test l.xmin ≈ xe1
    @test l.xmax ≈ xV2
    @test l.xmid ≈ (xV1 + xV2)/2
    @test l.yR ≈ 0
    @test l.ymin ≈ -d/2
    @test l.ymax ≈ d/2

    (xR1, xR2, xV1, xV2) = Optics.sphericalLensVertexThick(R1, R2, t)
    @test xR1 ≈ -2.0
    @test xR2 ≈ -1.0
    @test xV1 ≈ 0.0
    @test xV2 ≈ 1.0

    (xe1, xe2) = Optics.sphericalLensxEdge(R1, R2, d, t)
    @test xe1 ≈ -0.06350832689629149
    @test xe2 ≈ 0.9364916731037085

    l = Optics.SphericalLens(R1, R2, t, d, nl)
    @test l.R1 ≈ R1
    @test l.R2 ≈ R2
    @test l.t ≈ t
    @test l.d ≈ d
    @test l.nl ≈ nl
    @test l.xV1 ≈ xV1
    @test l.xV2 ≈ xV2
    @test l.xR1 ≈ xR1
    @test l.xR2 ≈ xR2
    @test l.xe1 ≈ xe1
    @test l.xe2 ≈ xe2
    @test l.xmin ≈ xe1
    @test l.xmax ≈ xV2
    @test l.xmid ≈ (xV1 + xV2)/2
    @test l.yR ≈ 0
    @test l.ymin ≈ -d/2
    @test l.ymax ≈ d/2
end

max_stars_repo_name: larschrj/Optics


text:
export Horocycle, RandomHorocycle

"""
`Horocycle(P::HPoint, theta::Real)` create the horocycle
containing the point `P` and the ideal point at `exp(im*theta)`.
"""
struct Horocycle <: HObject
    pt::HPoint
    t::Float64
    attr::Dict{Symbol,Any}
    function Horocycle(P::HPoint, theta::Real)
        theta = mod(theta,2pi)
        f = rotation(-theta)
        g = rotation(theta)
        # pt = P # put something

        z = f(getz(P))  # map to +x axis
        x,y = reim(z)

        a = (x^2 + y^2 - x)/(x - 1)
        w = a + 0im
        PP = g(HPoint(w))

        HC = new(PP,theta,Dict{Symbol,Any}())
        set_color(HC)
        set_thickness(HC)
        set_line_style(HC)
        return HC
    end
end
Horocycle() = Horocycle(HPoint(), 0)

Horocycle(t::Real,p::HPoint) = Horocycle(p,t)

# Copy constructor
Horocycle(HC::Horocycle) = Horocycle(HC.pt, HC.t)

"""
`RandomHorocycle()` creates a random horocycle by choosing
a point at random by `RandomHPoint` and a random ideal
point (uniformly between 0 and 2π).
"""
RandomHorocycle() = Horocycle(RandomHPoint(), 2*pi*rand())

function show(io::IO, HC::Horocycle)
    print(io,"Horocycle($(HC.pt),$(HC.t))")
end

function euclidean_center(HC::Horocycle)::Complex
    z = getz(HC.pt)
    w = exp(im*HC.t)
    return (z+w)/2
end

function(==)(H1::Horocycle, H2::Horocycle)
    if abs(H1.t-H2.t) > THRESHOLD * eps(1.0)
        return false
    end
    return H1.pt == H2.pt
end


text:
"""
    add_noise(S,0.05)
    add_noise(-||-, k=0.15)

A function that adds noise to a signal in a consistent way. The function adds both background noise as well as
noise to the signal itself. The level of the background noise and the signal noise is determined by the b
and k variables respectively. The levels are input as a fraction of the mean of the original signal S strength.
"""
function add_noise(s, b; k=0.1)
    sn = max.(s,0.0) .+ k.*(mean(sqrt.(abs.(s)))).*rand.(Normal.(0.0, max.(sqrt.(max.(s,0)), sqrt.(b))))
    err = k.*mean(sqrt.(abs.(s))).*max.(sqrt.(max.(s,0)), sqrt.(b))
    return sn, err
end

"""
    get_orbel_volume(myOrbitGrid, os_equidistant)

Get the orbit element volume of this specific orbit-grid.
If os_equidistant=false, then return a 3D array with all the orbit volumes.
Otherwise, just a scalar.
"""
function get_orbel_volume(og::OrbitGrid, os_equidistant::Bool)

    if os_equidistant
        dRm = abs(og.r[2]-og.r[1]) # The difference between the first and second element should be representative for the whole grid, due to equidistancy
        dE = abs(og.energy[2]-og.energy[1])
        dpm = abs(og.pitch[2]-og.pitch[1])
        dO = dE*dpm*dRm
        return dO
    else
        dO = zeros(length(og.energy), length(og.pitch), length(og.r)) # If not equidistant, create a 3D array with zeros
        for Ei=1:length(og.energy) # For all energies...
            if Ei==length(og.energy) # If at the edge of the orbit grid...
                # Assume edge orbit-element volume to be same as next-to-edge
                dEi = abs(og.energy[end]-og.energy[end-1])
            else
                dEi = abs(og.energy[Ei+1]-og.energy[Ei])
            end
            for pmi=1:length(og.pitch)
                if pmi==length(og.pitch)
                    dpmi = abs(og.pitch[end]-og.pitch[end-1])
                else
                    dpmi = abs(og.pitch[pmi+1]-og.pitch[pmi])
                end
                for Rmi=1:length(og.r)
                    if Rmi==length(og.r)
                        dRmi = abs(og.r[end]-og.r[end-1])
                    else
                        dRmi = abs(og.r[Rmi+1]-og.r[Rmi])
                    end

                    dO[Ei, pmi, Rmi] = dEi*dpmi*dRmi # Replace the zero-element in the 3D array with the resulting orbit-element volume
                end
            end
        end

        return dO
    end
end


text:
function Euler4()
	ThreeDigits = [100:999]
	Products = unique(reshape(ThreeDigits*ThreeDigits',1,length(ThreeDigits)^2))
	return maximum(filter(x -> dec(x)==dec(x)[end:-1:1], Products))
end

@time println(Euler4())

max_stars_repo_name: gribeill/ProjectEuler

max_issues_repo_name: gribeill/ProjectEuler

max_forks_repo_name: gribeill/ProjectEuler


text:
module UTF16

import Base: UTF16String, utf16
export UTF16String, utf16

end


text:
using NeXLParticle
using Test

@testset "Blob" begin
    spec=loadspectrum(ASPEXTIFF, joinpath(@__DIR__, "images", "test1.tif"), withImgs=true)
    bse = spec[:Image1]

    bs = blob(bse, i->i>0.7)
    @test area.(bs) == [ 2882, 43, 1, 1, 1, 1]
    @test perimeter(bs[1])[24]==CartesianIndex(66,12)
    @test length(perimeter(bs[1]))==240
    @test size(CartesianIndices(bs[1]))==(69, 66)
    @test isapprox(perimeterlength(bs[1]),292.19,atol=0.01)

    @test isapprox(ecd(bs[1],true), 60.81, atol=0.01)
    @test isapprox(ecd(bs[1]), 60.81, atol=0.01)
    @test isapprox(ecd(bs[1],false), 60.58, atol=0.01)

    @test length(curvature(bs[1], 8))==length(perimeter(bs[1]))
    bs2=separate(bs[1])
    @test length(bs2)==10
    @test area.(bs2) == [ 1415, 836, 183, 227, 57, 2, 3, 3, 1, 1 ]
    @test_nowarn colorizedimage(bs, bse)
    @test_nowarn colorizedimage(bs2, bse)

    ms=multiseparate(bse, 0.6:0.05:1.0)

    ir=interiorregions(bs[1])
    @test area.(ir)== [15,2,1,1,1,1,1,1]
    @test NeXLParticle.filledarea(bs[1]) > area(bs[1])
    @test NeXLParticle.filledarea(bs[1]) == area(bs[1]) + sum(area.(ir))
    @test isapprox(sum(maskedimage(bs[1], bse), init=1.0), 2798.8, atol=0.1)
end

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
import Base.convert

export mat, vecs
export SVec, SMat, STen, JVec, JVecs, JVecF, JVecsF
export JMat, JMatF
export zerovecs, maxdist, maxnorm

"typealias fora static vector type; currently `StaticArrays.SVector`"
const SVec = SVector
"typealias fora static matrix type; currently `StaticArrays.SMatrix`"
const SMat = SMatrix
"typealias fora static array type; currently `StaticArrays.SArray`"
const STen = SArray

"`JVec{T}` : 3-dimensional immutable vector"
const JVec{T} = SVec{3,T}
const JVecF = JVec{Float64}
const JVecI = JVec{Int}

# Base.zero{T}(::Type{JVec{T}}) = JVec([zero(T) for i=1:3])

# which of these types shall we deprecate?

"`JVecs{T}` : List of 3-dimensional immutable vectors"
const JVecs{T} = Vector{JVec{T}}
const JVecsF = JVecs{Float64}
const JVecsI = JVecs{Int}

"""
`vecs(V::Matrix)` : convert (as reference) a 3 x N matrix representing
N vectors (e.g. forces) in R³ to a list (vector) of fixed-size-array vectors.

`vecs(V::Vector)` : assumes that V is morally a 3 x N matrix, stored in a long
vector

`vecs(V::Array{T,N})` : If `V` has dimensions 3 x n2 x ... x nN then
it gets converted to an n2 x ... x nN array with JVec{T} entries.
"""
vecs{T}(V::Matrix{T}) = reinterpret(JVec{T}, V, (size(V,2),))
vecs{T}(V::Vector{T}) = reinterpret(JVec{T}, V, (length(V) ÷ 3,))
vecs{T,N}(V::Array{T,N}) = reinterpret(JVec{T}, V, tuple(size(V)[2:end]...))

"`JMat{T}` : 3 × 3 immutable marix"

const JMat{T,N} = SMatrix{3,3,T,N}
const JMatF = SMatrix{3,3,Float64,9}
const JMatI = JMat{Int}

# Base.zero{T}(::Type{JMat{T}}) = JMat([zero(T) for i = 1:9])
# Base.eye{T}(::Type{JMat{T}}) = JMat(T, eye(3))

"`JMats{T}` : (2-dimensional) Array of 3 × 3 immutable matrices"
const JMats{T} = Array{JMat{T}}
const JMatsF = JMats{Float64}
const JMatsI = JMats{Int}
#
# """
# `mats(V::Matrix)` : convert (as reference) a 3 x 3 x N x N tensor representing
# N × N matrices (e.g. local hess) in R³ to a list (vector) of fixed-size-array vectors.
#
# `mats(V::Vector)` : assumes that V is morally a list of  x N matrix, stored in a long
# vector
#
# `mats(V::Array{T,N})` : If `V` has dimensions 3 x n2 x ... x nN then
# it gets converted to an n2 x ... x nN array with JVec{T} entries.
# """
# mats{T}(V::Array{T}) = reinterpret(JMat{T}, V, tuple(size(V)[3:end]...))
# mats{T}(V::Vector{T}) = reinterpret(JVec{T}, V, (length(V) ÷ 3,)
# mats{T,N}(V::Array{T,N}) = reinterpret(JMat{T}, V, tuple(size(V)[3:end]...))
# mats{T}(V::Array{T}) = vecs(vecs(V))
# mats{T}(V::Matrix{T}) = reshape(permutedims(V, [1 3 2 4]), 6, 6)

export mats
function mats{T,N}(V::Array{T, N})
  @assert size(V,1) == size(V,2) == 3
  @assert ndims(V) > 2
  return reinterpret(JMat{T}, V, tuple(size(V)[3:end]...))
end

export xyz

"""
convert a Vector{JVec{T}} or a 3 x N Matrix{T} into a 3-tuple
(x, y, z). E.g.,
```
x, y, z = positions(at) |> xyz
```
a short-cut is
```
x, y, z = xyz(at)
```
Conversely, `set_positions!(at, x, y, z)` is also allowed.
"""
xyz{T}(V::Vector{JVec{T}}) = xyz(mat(V))

function xyz(V::Matrix)
   @assert size(V, 1) == 3
   return (view(V, 1, :), view(V, 2, :), view(V, 3, :))
end

"""
`mat(X::JVecs)`: convert (as reference) a list (Vector) of
fixed-size-array points or vecs to a 3 x N matrix representing those.

### Usage:
```
X = positions(at)          # returns a Vector{JVec}
X = positions(at) |> mat   # returns a Matrix
```
"""
mat{N,T}(V::Vector{SVec{N,T}}) = reinterpret(T, V, (N, length(V)))
mat{N,T}(X::AbstractVector{SVec{N,T}}) = mat(collect(X))

# rewrite all of this in terms of `convert` (TODO: is this needed?)
convert{T}(::Type{Matrix{T}}, V::JVecs{T}) = mat(V)
convert{T}(::Type{JVecs{T}}, V::Matrix{T}) = vec(V)

# initialise a vector of vecs or points
zerovecs(n::Integer) = zerovecs(Float64, n)
zerovecs(T::Type, n::Integer) = zeros(T, 3, n) |> vecs

zeromats(T::Type, n::Integer) = zeros()
# initialise a matrix of

"""
maximum of distances between two sets of JVec's, usually positions;
`maximum(a - b for (a,b) in zip(x,y))`
"""
function maxdist(x::AbstractArray{JVec{T}}, y::AbstractArray{JVec{T}}) where T
   @assert length(x) == length(y)
   return maximum( norm(a-b)  for (a,b) in zip(x,y) )
end

"`maximum(norm(y) for y in x);` typically, x is a vector of forces"
maxnorm{T}(x::JVecs{T}) = maximum( norm.(x) )


text:
"""
    compute_mf_diffusion_tensor(mesh, M, lap_eig, gradient::ScalarGradient)

Compute effective diffusion tensors using the matrix formalism.
"""
function compute_mf_diffusion_tensor(mesh, M, lap_eig, gradient::ScalarGradient)
    volumes = get_cmpt_volumes(mesh)
    p = reduce(hcat, mesh.points)
    a = p * M * lap_eig.funcs
    j = j_integral.(gradient.profile, lap_eig.values)
    D = a * Diagonal(j) * a' / sum(volumes)
end


text:
module random

export randstate, randoperator

using ..bases, ..states, ..operators_dense

"""
    randstate(basis)

Calculate a random normalized ket state.
"""
function randstate(b::Basis)
    psi = Ket(b, rand(Complex128, length(b)))
    normalize!(psi)
    psi
end

"""
    randoperator(b1[, b2])

Calculate a random unnormalized dense operator.
"""
randoperator(b1::Basis, b2::Basis) = DenseOperator(b1, b2, rand(Complex128, length(b1), length(b2)))
randoperator(b::Basis) = randoperator(b, b)

end #module


text:
abstract type BrainFlowLib end
struct BoardControllerLib <: BrainFlowLib end
struct DataHandlerLib <: BrainFlowLib end
struct MlModuleLib <: BrainFlowLib end
library_path(::BoardControllerLib) = BOARD_CONTROLLER_INTERFACE
library_path(::DataHandlerLib) = DATA_HANDLER_INTERFACE
library_path(::MlModuleLib) = ML_MODULE_INTERFACE

const BOARD_CONTROLLER = BoardControllerLib()
const DATA_HANDLER = DataHandlerLib()
const ML_MODULE = MlModuleLib()

@enum LogLevels begin

    LEVEL_TRACE = 0
    LEVEL_DEBUG = 1
    LEVEL_INFO = 2
    LEVEL_WARN = 3
    LEVEL_ERROR = 4
    LEVEL_CRITICAL = 5
    LEVEL_OFF = 6

end

enable_brainflow_logger(log_lib::BrainFlowLib) = set_log_level(Integer(LEVEL_INFO), log_lib)
enable_dev_brainflow_logger(log_lib::BrainFlowLib) = set_log_level(Integer(LEVEL_TRACE), log_lib)
disable_brainflow_logger(log_lib::BrainFlowLib) = set_log_level(Integer(LEVEL_OFF), log_lib)

# available only for BoardController, no need to provide log_lib
@brainflow_rethrow function log_message(log_level::Integer, message::String)
    ccall((:log_message, BOARD_CONTROLLER_INTERFACE), Cint, (Cint, Ptr{UInt8}), Int32(log_level), message)
end

@brainflow_rethrow function set_log_file(log_file, log_lib::BrainFlowLib)
    lib_cglobal = log_file_cglobal(log_lib)
    ccall(lib_cglobal, Cint, (Ptr{UInt8},), log_file)
end

# need to hardcode the cglobal input symbols, related to https://github.com/JuliaLang/julia/issues/29602
log_file_cglobal(::BoardControllerLib) = cglobal((:set_log_file, BOARD_CONTROLLER_INTERFACE))
log_file_cglobal(::DataHandlerLib) = cglobal((:set_log_file, DATA_HANDLER_INTERFACE))
log_file_cglobal(::MlModuleLib) = cglobal((:set_log_file, ML_MODULE_INTERFACE))

@brainflow_rethrow function set_log_level(log_level::Integer, log_lib::BrainFlowLib)
    lib_cglobal = log_level_cglobal(log_lib)
    ccall(lib_cglobal, Cint, (Cint,), Int32(log_level))
end

# need to hardcode the cglobal input symbols, related to https://github.com/JuliaLang/julia/issues/29602
log_level_cglobal(::BoardControllerLib) = cglobal((:set_log_level, BOARD_CONTROLLER_INTERFACE))
log_level_cglobal(::DataHandlerLib) = cglobal((:set_log_level, DATA_HANDLER_INTERFACE))
log_level_cglobal(::MlModuleLib) = cglobal((:set_log_level, ML_MODULE_INTERFACE))

max_stars_repo_name: matthijscox/brainflow

max_issues_repo_name: matthijscox/brainflow

max_forks_repo_name: matthijscox/brainflow


text:
module DataLoaders

include("sbloader.jl")
using .Statsbomb

end ##end of module


text:
@testset "Rotations" begin
    v₁ = [3, 2, 1]
    v₂ = [0.5, 0.8, 0.1]

    # Datamine convention
    datamine_angles = DatamineAngles(70, 30, -20)
    datamine_dcm = angle_to_dcm(datamine_angles)

    @test datamine_dcm * v₁ ≈ [2.533789309500169, 1.6971296470206085, 2.1678705441668704]
    @test datamine_dcm * v₂ ≈ [0.593857346855506, 0.7024802506855807, 0.2320666908077572]

    # GSLIB convention
    gslib_angles = GslibAngles(70, -30, 20)
    gslib_dcm = angle_to_dcm(gslib_angles)

    @test gslib_dcm * v₁ ≈ [2.533789309500169, 1.6971296470206085, 2.1678705441668704]
    @test gslib_dcm * v₂ ≈ [0.593857346855506, 0.7024802506855807, 0.2320666908077572]

    # Vulcan convention
    vulcan_angles = VulcanAngles(70, -30, -20)
    vulcan_dcm = angle_to_dcm(vulcan_angles)

    @test vulcan_dcm * v₁ ≈ [2.533789309500169, 1.6971296470206085, 2.1678705441668704]
    @test vulcan_dcm * v₂ ≈ [0.593857346855506, 0.7024802506855807, 0.2320666908077572]

    # comparing conventions
    @test datamine_dcm ≈ gslib_dcm
    @test datamine_dcm ≈ vulcan_dcm
    @test gslib_dcm ≈ vulcan_dcm
end


text:
"""
    shiftRight!(matrixIn,nn)

  Shift the rows of matrixIn to the right by nn columns, leaving zeros in the
  first nn columns.
"""
function shiftRight!(matrixIn::Array{Float64,2}, nn::Int64)

	 numColumns = size(matrixIn, 2)

	 # left side of input matrix that holds data from column 1 to nn
	 left = matrixIn[:, 1:(numColumns - nn)]

	 # initialize new array same size as input matrix to all zeros
	 result = zero(matrixIn)

	 # replace right size of input matrix with the first nn columns
	 result[:, (nn + 1):numColumns] = left

	 return result
end # shiftRight


text:
# just use the real part for LS.
#yr, yi assumed to have LS_inds applied already.
# This does A\b but with bounds.
function solveBLScL( W::Matrix{Complex{T}},
    b,
    x_lower::Vector{T},
    x_upper::Vector{T}) where T <: Real

    # This is the formula.
    # BoundedLeastSquares.Quadratic(A'*A, A'*b)

    A = reinterpret(T, W)

    AtA = forcesymmetric(A'*A)
    Aty = A'*b

    Q = BoundedLeastSquares.Quadratic(AtA, Aty)

    x_opt, out_flag = BoundedLeastSquares.min_bound_constrained_quadratic(Q,
                        x_lower, x_upper)

    return x_opt, out_flag
end

####### w.

function setupupdateLS(N_positions::Int, N_compounds::Int, y::Vector{Complex{T}}) where T <: Real

    A = Matrix{Complex{T}}(undef, N_positions, N_compounds)
    w = zeros(T, N_compounds)
    b = collect(reinterpret(T, y))

    return A, w, b
end

"""
A is an M x L complex-valued matrix, b is 2*M 1D array.
Updates solution to w (mutates it), and mutates intermediate variables A, b.
"""
function updatew!(  A::Matrix{Complex{T}},
    b,
    w::Vector{T},
    U_rad_LS,
    Es,
    As,
    w_lower::Vector{T},
    w_upper::Vector{T}) where T <: Real

    @assert length(b) == 2*length(U_rad_LS)

    evaldesignmatrixw2!(A, U_rad_LS, Es, As)

    # ## check.
    # A2 = copy(A)
    # evaldesignmatrixw!(A, U_rad_LS, Es, As)
    #
    # if norm(A-A2) > 1e-10
    #     println("problem)")
    # end

    ### LS solve.
    w[:], status_flag = solveBLScL(A, b, w_lower, w_upper)

    return status_flag
end

function evaldesignmatrixw!(B::Matrix{Complex{T}},
    U_rad, Es, As) where T <: Real

    M = length(U_rad)
    N = length(Es)

    @assert size(B) == (M,N)

    fill!(B, NaN) # debug.

    for n = 1:N
        for m = 1:M
            # speed up later.
            B[m,n] = NMRSpectraSimulator.evalitpproxycompound(U_rad[m], As[n], Es[n])
        end
    end

    return nothing
end

function evaldesignmatrixw2!(B::Matrix{Complex{T}},
    U_rad, Es, As) where T <: Real

    M = length(U_rad)
    N = length(Es)

    @assert size(B) == (M,N)

    fill!(B, NaN) # debug.

    for n = 1:N
        E = Es[n]
        A = As[n]

        for m = 1:M

            s = zero(Complex{T})
            for i in eachindex(E.κs_α)
                for k in eachindex(E.κs_α[i])
                    s += E.core.qs[i][k](U_rad[m] - E.core.ss_params.d[i], E.core.ss_params.κs_λ[i])
                end
            end
            B[m,n] = s

            #B[m,n] = sum( sum( E.core.qs[i][k](U_rad[m] - E.core.ss_params.d[i], E.core.ss_params.κs_λ[i]) for k in eachindex(E.κs_α[i])) for i in eachindex(E.κs_α))

            # singlets.
            for k = 1:length(E.κs_α_singlets)

                B[m,n] += NMRSpectraSimulator.evalsinglets(U_rad[m], E.core.d_singlets,
                A.αs_singlets, A.Ωs_singlets,
                E.core.β_singlets, E.core.λ0, E.core.κs_λ_singlets, E.κs_α_singlets)
            end

        end
    end

    # if !all(isfinite.(B))
    #     println("problem")
    # end

    return nothing
end

######### κ, compensation for α.

function updateκ2!(  A::Matrix{Complex{T}},
    #r_buffer::Vector{T},
    b,
    κ::Vector{T},
    U_rad_LS,
    Es,
    As,
    w::Vector{T},
    κ_lower::Vector{T},
    κ_upper::Vector{T}) where T <: Real

    #@assert size(A) == (length(b), length(αS))
    @assert length(b) == 2*length(U_rad_LS)

    evaldesignmatrixκ!(A, U_rad_LS, Es, As, w)

    ### LS solve.
    # w[:] = NonNegLeastSquares.nonneg_lsq(A, b; alg = :fnnls)
    # status_flag = true

    κ[:], status_flag = solveBLScL(A, b, κ_lower, κ_upper)
    # κ[:] = A\b
    # clamp!(κ, κ_lower[1], κ_upper[1])
    # status_flag = true

    return status_flag
end

function evaldesignmatrixκ!(B::Matrix{Complex{T}},
    U_rad,
    Es,
    As,
    w::Vector{T}) where T <: Real

    #
    M = length(U_rad)
    N = length(Es)

    N_κ, N_κ_singlets = countκ(Es)
    @assert size(B) == (M, N_κ + N_κ_singlets)
    fill!(B, Inf) # debug.

    resetκ!(Es)
    j = 0

    # loop over each κ partition element in Es.
    for n = 1:N
        E = Es[n]
        A = As[n]
        @assert length(E.κs_α) == length(E.core.qs)

        #spin system.
        #for i = 1:length(E.κs_α)
        for i in eachindex(E.κs_α)

            # partition
            #for k = 1:length(E.κs_α[i])
            for k in eachindex(E.κs_α[i])

                j += 1

                # # loop over each fit position.
                # for m = 1:M
                #
                #     # taken from evalitproxysys()
                #     r = U_rad[m] - E.core.ss_params.d[i]
                #
                #     #B[m,j] = w[n]*E.core.qs[i][k](r, E.core.ss_params.κs_λ[i])
                #     #w[n]*E.core.qs[i][k](r, E.core.ss_params.κs_λ[i])
                #     #B[m,j] = 1.9
                #
                #     out = w[n]*E.core.qs[i][k](r, E.core.ss_params.κs_λ[i])
                #
                #     #Q[j,m] = w[n]*E.core.qs[i][k](r, E.core.ss_params.κs_λ[i])
                # end

                #out = collect( w[n]*E.core.qs[i][k](U_rad[m] - E.core.ss_params.d[i], E.core.ss_params.κs_λ[i]) for m = 1:M )
                B[:,j] = collect( w[n]*E.core.qs[i][k](U_rad[m] - E.core.ss_params.d[i], E.core.ss_params.κs_λ[i]) for m = 1:M )
            end
        end

        # singlets.
        for k = 1:length(E.κs_α_singlets)
            j += 1

            # for m = 1:M
            #
            #     B[m,j] = w[n]*NMRSpectraSimulator.evalsinglets(U_rad[m], E.core.d_singlets,
            #     A.αs_singlets, A.Ωs_singlets,
            #     E.core.β_singlets, E.core.λ0, E.core.κs_λ_singlets, E.κs_α_singlets)
            # end

            B[:,j] = collect( w[n]*NMRSpectraSimulator.evalsinglets(U_rad[m], E.core.d_singlets,
            A.αs_singlets, A.Ωs_singlets,
            E.core.β_singlets, E.core.λ0, E.core.κs_λ_singlets, E.κs_α_singlets) for m = 1:M )

        end
    end

    #println("all(isfinite.(B)) = ", all(isfinite.(B)))
    @assert all(isfinite.(B))

    return nothing
end


text:
"""
Simple flat basis (all symmetry blocks) intended for iteration.
"""
struct FlatBasis{N}
    "Local basis truncation."
    l_max::Int
    "Many-body basis truncation."
    l_total_max::Int
end

mutable struct FlatBasisIterState
    v::Vector{Int}
    l_total::Int
end

function FlatBasisIterState(::FlatBasis{N}) where {N}
    v = zeros(Int, 2N)
    # m = -1
    v[2] -= 1

    l_total = 0

    FlatBasisIterState(v, l_total)
end

function Base.iterate(fb::FlatBasis{N}, state::FlatBasisIterState=FlatBasisIterState(fb)) where {N}
    for i in 1:N
        # m += 1
        state.v[2i] += 1

        # m > l
        if state.v[2i] > state.v[2i-1]
            state.l_total += 1
            # l += 1
            state.v[2i-1] += 1

            # l > l_max || l_total > l_total_max
            if state.v[2i-1] > fb.l_max || state.l_total > fb.l_total_max
                state.l_total -= state.v[2i-1]
                # l = 0
                state.v[2i-1] = 0
                # m = 0
                state.v[2i] = 0

                i == N && return nothing

                continue
            else
                # m = -l
                state.v[2i] = -state.v[2i-1]
            end
        end

        break
    end

    state.v, state
end

Base.eltype(::Type{FlatBasis}) = Vector{Int}

"""
Symmetry block label.
"""
struct BlockLabel
    "Total l parity."
    lp::Int
    "Total m value."
    m::Int
end

"""
Symmetry block of a basis for `N` rotors.
"""
struct Block{N}
    "Basis vectors, each of the form [l_1, m_1, ..., l_N, m_N]."
    vectors::Vector{Vector{Int}}
    "Mapping from basis vectors to their indices."
    lookup::Dict{Vector{Int},Int}

    "Number of basis vectors."
    size::Int
end

"""
    Block(N::Int, vectors::Vector{Vector{Int}}, lookup::Dict{Vector{Int},Int})

Create a block of `vectors` with corresponding `lookup` table for `N` rotors.
"""
function Block(N::Int, vectors::Vector{Vector{Int}}, lookup::Dict{Vector{Int},Int})
    # For non-empty blocks, make sure N is correct.
    length(vectors) == 0 || length(vectors[1]) == 2N || throw(DomainError(N, "Invalid number of rotors."))

    Block{N}(vectors, lookup, length(vectors))
end

abstract type AbstractBasis{N} end

"""
Truncated basis for `N` rotors, containing a single symmetry block.
"""
struct SingleBlockBasis{N} <: AbstractBasis{N}
    "Local basis truncation."
    l_max::Int
    "Many-body basis truncation."
    l_total_max::Int

    "Symmetry label for the block."
    label::BlockLabel
    "The block of basis vectors."
    block::Block{N}

    "Total number of basis vectors."
    size::Int
end

"""
    SingleBlockBasis(N::Int, l_max::Int, l_total_max::Int, sym_lp::Int, sym_m::Int)

Generate a truncated `N`-body basis of spherical harmonics.

The local bases are truncated at `l_max`. The many-body basis is truncated so
that the sum of `l` values does not exceed `l_total_max`.

Only states with total l parity `sym_lp` and total m value `sym_m` are
included.
"""
function SingleBlockBasis(N::Int, l_max::Int, l_total_max::Int, sym_lp::Int, sym_m::Int)
    N >= 1 || throw(DomainError(N, "At least one rotor."))
    l_max >= 0 || throw(DomainError(l_max, "Non-positive l_max."))
    l_total_max >= 0 || throw(DomainError(l_total_max, "Non-positive l_total_max."))

    vectors = Vector{Int}[]
    lookup = Dict{Vector{Int},Int}()

    for v in FlatBasis{N}(l_max, l_total_max)
        lp_total = sum(v[2i-1] for i in 1:N) % 2
        m_total = sum(v[2i] for i in 1:N)

        lp_total == sym_lp || continue
        m_total == sym_m || continue

        vv = copy(v)
        push!(vectors, vv)
        lookup[vv] = length(vectors)
    end

    label = BlockLabel(sym_lp, sym_m)
    block = Block(N, vectors, lookup)

    SingleBlockBasis{N}(l_max, l_total_max, label, block, block.size)
end

"""
Truncated basis for `N` rotors, containing multiple symmetry blocks.
"""
struct MultiBlockBasis{N} <: AbstractBasis{N}
    "Local basis truncation."
    l_max::Int
    "Many-body basis truncation."
    l_total_max::Int

    "Blocks of basis vectors."
    blocks::Dict{BlockLabel,Block{N}}

    "Total number of basis vectors."
    size::Int
end

"""
    MultiBlockBasis(N::Int, l_max::Int, l_total_max::Int)

Generate a truncated `N`-body basis of spherical harmonics.

The local bases are truncated at `l_max`. The many-body basis is truncated so
that the sum of `l` values does not exceed `l_total_max`.
"""
function MultiBlockBasis(N::Int, l_max::Int, l_total_max::Int)
    N >= 1 || throw(DomainError(N, "At least one rotor."))
    l_max >= 0 || throw(DomainError(l_max, "Non-positive l_max."))
    l_total_max >= 0 || throw(DomainError(l_total_max, "Non-positive l_total_max."))

    # For each site, the magnitude of m is bounded by l, so (by the triangle
    # inequality) the magnitude of the sum of the m values can't exceed the sum
    # of the l values, which is in turn bounded by l_total_max.
    vectors = Dict(((lp, m), Vector{Int}[]) for lp in 0:1 for m in -l_total_max:l_total_max)
    lookups = Dict(((lp, m), Dict{Vector{Int},Int}()) for lp in 0:1 for m in -l_total_max:l_total_max)

    for v in FlatBasis{N}(l_max, l_total_max)
        lp_total = sum(v[2i-1] for i in 1:N) % 2
        m_total = sum(v[2i] for i in 1:N)

        vv = copy(v)
        push!(vectors[(lp_total, m_total)], vv)
        lookups[(lp_total, m_total)][vv] = length(vectors[(lp_total, m_total)])
    end

    blocks = Dict{BlockLabel,Block{N}}()
    size = 0
    for (lp, m) in keys(vectors)
        label = BlockLabel(lp, m)
        blocks[label] = Block(N, vectors[(lp, m)], lookups[(lp, m)])
        size += blocks[label].size
    end

    MultiBlockBasis{N}(l_max, l_total_max, blocks, size)
end


text:
module TestBouncingBall

using ModiaLang
using DifferentialEquations
@usingModiaPlot

# Does not work: edge(..) appears in a linear equation system

BouncingBall = Model(
    e = 0.7,
    g = 9.81,
    h = Var(init = 1.0),
    v = Var(init = 0.0),
    flying = Var(start = true),   # due to pre(flying) -> SimulationModel(....; pre_startValues = [true], ....)
    equations = :[
        # desired:
        #   flying = edge(-h) ? reinit(v, -e*v) : pre(flying)
        flying = edge(instantiatedModel, 1, -h, "-h", _leq_mode) ?
                    !reinit(instantiatedModel, _x, 2, -e*v, _leq_mode) : pre(flying),
        der(h) = v,
        der(v) = flying ? -g : 0.0
    ]
)

model = @instantiateModel(BouncingBall, logCode=true)
@show model.equationInfo.x_info

# Temporary code, until edge, reinit and pre are supported in the Modia language
    eh = model.eventHandler
    eh.nz = 1
    eh.z  = ones(eh.nz)
    eh.zPositive = fill(false, eh.nz)
    model.pre = [true]

simulate!(model, Tsit5(), stopTime = 3.0, log=true, logEvents=true)   # requiredFinalStates = [-0.3617373025974107]

plot(model, ["h", "v", "flying"])

end

max_stars_repo_name: ModiaSim/ModiaLang

max_issues_repo_name: ModiaSim/ModiaLang

max_forks_repo_name: ModiaSim/ModiaLang


text:
export get_TD_operator

"""
   input:
# Arguments
  - TD_type   : string indicates which transform-domain operator to return
  - comp_grid : structure with computational grid information; comp_grid.n = [n1 n2] number of gridpoints in each direction; comp_grid.d = [d1 d2] spacing between gridpoints in each direction.

   output
   - TD_OP   : transform domain operator as a sparse matrix or matrix-free object
"""
function get_TD_operator(comp_grid,TD_type,TF)

# if TF==Float64
#   TI=Int64
# else
#   TI=Int32
# end
TI=Int64

h1 = TF(comp_grid.d[1])
h2 = TF(comp_grid.d[2])
n1 = comp_grid.n[1]
n2 = comp_grid.n[2]

if length(comp_grid.n)==3 && comp_grid.n[3]>1 #use 3d version
  h3 = TF(comp_grid.d[3])
  n3 = comp_grid.n[3]
  #(D3D, D3x, D3y, D3z) = get_discrete_Grad(n1,n2,n3,h1,h2,h3);
  if TD_type=="TV" || TD_type=="D3D"# anisotropic TV operator
      TD_OP = get_discrete_Grad(n1,n2,n3,h1,h2,h3,TD_type)
      AtA_diag = false ; dense = false ; TD_n = (n1-1+n1+n1 , n2-1+n2+n2 , n3-1+n3+n3) ; banded=true
  elseif TD_type=="D_z" # vertical derivative operator
      TD_OP = get_discrete_Grad(n1,n2,n3,h1,h2,h3,TD_type)
      AtA_diag = false ; dense = false ; TD_n = (n1 , n2 , n3-1) ; banded=true
  elseif TD_type=="D_x" # lateral derivative operator
      TD_OP = get_discrete_Grad(n1,n2,n3,h1,h2,h3,TD_type)
      AtA_diag = false ; dense = false ; TD_n = (n1-1 , n2 , n3) ; banded=true
  elseif TD_type=="D_y" # lateral derivative operator
      TD_OP = get_discrete_Grad(n1,n2,n3,h1,h2,h3,TD_type)
      AtA_diag = false ; dense = false ; TD_n = (n1 , n2-1 , n3) ; banded=true
  elseif TD_type == "identity"
      TD_OP = SparseMatrixCSC{TF}(LinearAlgebra.I,n1*n2*n3,n1*n2*n3)
      AtA_diag = true ; dense = false ; TD_n = (n1,n2,n3) ; banded=true
  elseif TD_type == "DFT"
      TD_OP=joDFT(convert(Int64,n1),convert(Int64,n2),convert(Int64,n3);planned=false,DDT=TF,RDT=Complex{TF})
      AtA_diag = true ; dense = true ; TD_n = (n1,n2,n3) ; banded=false
    elseif TD_type == "DCT"
      TD_OP=joDCT(convert(Int64,n1),convert(Int64,n2),convert(Int64,n3);planned=false,DDT=TF,RDT=TF)
      AtA_diag = true ; dense = true ; TD_n = (n1,n2,n3) ; banded=false
  elseif TD_type == "curvelet"
      error("currently no 3D curvelet implemented, needs to be done soon")
    elseif TD_type == "wavelet"
      error("currently no 3D wavelet implemented, needs to be done soon")
  else
      error("provided an unknown transform domain operator. check function get_TD_operator(comp_grid,TD_type,TF) for options")
  end
else #use 2d version
  #(D2D, D2x, D2z) = get_discrete_Grad(n1,n2,h1,h2);
  if TD_type=="TV" || TD_type=="D2D"# anisotropic TV operator
      TD_OP = get_discrete_Grad(n1,n2,h1,h2,TD_type)
      AtA_diag = false ; dense = false ; TD_n = ( (n1-1)+n1 , n2+(n2-1) ) ; banded=true
  elseif TD_type=="D_z" # vertical derivative operator
      TD_OP = get_discrete_Grad(n1,n2,h1,h2,TD_type)
      AtA_diag = false ; dense = false ; TD_n = (n1 , n2-1) ; banded=true
  elseif TD_type=="D_x" # lateral derivative operator
      TD_OP = get_discrete_Grad(n1,n2,h1,h2,TD_type)
      AtA_diag = false ; dense = false ; TD_n = (n1-1 , n2) ; banded=true
  elseif TD_type=="D_xz" # lateral derivative operator composed with vertical derivative TD_OP=D_z*D_x
      D_x = get_discrete_Grad(n1,n2,h1,h2,"D_x")
      D_z = get_discrete_Grad(n1-1,n2,h1,h2,"D_z") #slightly different dimensions corresponding to the size of D_x*v (v=vector)
      TD_OP = D_z*D_x
      AtA_diag = false ; dense = false ; TD_n = (n1-1 , n2-1) ; banded=true
  elseif TD_type == "identity"
    TD_OP = SparseMatrixCSC{TF}(LinearAlgebra.I,n1*n2,n1*n2)
    AtA_diag = true ; dense = false ; TD_n = (n1,n2) ; banded=true
  elseif TD_type == "DCT"
    TD_OP=joDCT(convert(Int64,n1),convert(Int64,n2);planned=false,DDT=TF,RDT=TF)
    AtA_diag = true ; dense = true ; TD_n = (n1,n2) ; banded=false
  elseif TD_type == "DFT"
    TD_OP=joDFT(convert(Int64,n1),convert(Int64,n2);planned=false,DDT=TF,RDT=Complex{TF})
    AtA_diag = true ; dense = true ; TD_n = (n1,n2) ; banded=false
  elseif TD_type == "curvelet"
    TD_OP = joCurvelet2D(convert(Int64,n1),convert(Int64,n2);DDT=TF,RDT=TF,real_crvlts=true,all_crvlts=true)
    AtA_diag = true ; dense = true ; TD_n = (0,0) ; banded=false
  elseif TD_type == "wavelet"
     TD_OP = joDWT(n1,n2,wavelet(WT.db4);L=maxtransformlevels(min(n1,n2)),DDT=TF,RDT=TF)
     AtA_diag = true ; dense = true ; TD_n = (n1,n2) ; banded=false
  else
    error("provided an unknown transform domain operator. check function get_TD_operator(comp_grid,TD_type,TF) for options")
  end
end

return TD_OP, AtA_diag, dense, TD_n, banded
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
module NISTStatisticalReferenceData

using Downloads
using OrderedCollections: OrderedDict
using CSV
using DataFrames
using Parsers

depsdir(args...) = abspath(joinpath(@__DIR__, "..", "deps", args...))

function datasets()
    out = OrderedDict{String,Vector{String}}()
    for dir in ["anova", "mcmc", "nonlinear", "regression", "univariate"]
        dsets = map(x -> replace(x, ".dat" => ""), readdir(depsdir(dir)))
        append!(get!(out, dir, String[]), dsets)
    end
    return out
end

#-----------------------------------------------------------------------------# TestCase
struct TestCase
    name::String
    localfile::String
    url::String
    format::String
    procedure::String
    reference::String
    data_description::String
    model::String
    data::DataFrame
    certified_values::Vector
end
function TestCase(name::String)
    name = endswith(name, ".dat") ? name : name * ".dat"
    if isfile(depsdir("anova", name))
        url = "https://www.itl.nist.gov/div898/strd/anova/$name"
        localfile = depsdir("anova", name)
    elseif isfile(depsdir("mcmc", name))
        url = "https://www.itl.nist.gov/div898/strd/mcmc/$name"
        localfile = depsdir("mcmc", name)
    elseif isfile(depsdir("nonlinear", name))
        url = "https://www.itl.nist.gov/div898/strd/nls/data/LINKS/DATA/$name"
        localfile = depsdir("nonlinear", name)
    elseif isfile(depsdir("regression", name))
        url = "https://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA/$name"
        localfile = depsdir("regression", name)
    elseif isfile(depsdir("univariate", name))
        url = "https://www.itl.nist.gov/div898/strd/univ/data/$name"
        localfile = depsdir("univariate", name)
    else
        error("No dataset with the name \"$name\" was found.")
    end
    lines = readlines(localfile)
    TestCase(name, localfile, url, get_format(lines), get_procedure(lines), get_reference(lines),
        get_data_description(lines), get_model(lines), get_data(lines), get_certified_values(lines))
end

function extract_section(lines, heading::String; leftalign=true, op=findfirst)
    i = op(startswith(heading * ":"), lines)
    isnothing(i) && return "Univariate"
    line = replace(lines[i], "$heading:" => "")
    j = findnext(!startswith(r" |-"), lines, i + 1)
    out = isnothing(j) ? vcat(line, lines[i+1:end]) : vcat(line, lines[i+1:j-1])
    out = leftalign ? lstrip.(out) : out
    replace(join(out, '\n'), "$heading:"=>"")
end

function get_format(lines)
    extract_section(lines, "File Format")
end

function get_procedure(lines)
   extract_section(lines, "Procedure")
end

function get_reference(lines)
    extract_section(lines, "Reference")
end

function get_data_description(lines)
    extract_section(lines, "Data")
end

function get_model(lines)
    extract_section(lines, "Model")
end

function get_data(lines)
    # Convert from fixed width to CSV
    s = replace(extract_section(lines, "Data"; op=findlast), r"[ \t]+" => ',')
    CSV.read(IOBuffer(s), DataFrame)
end

function get_certified_values(lines)
    out = []
    for T in [AnovaTable, Coefficients, RSquared, ResidualStandardError, MCMCStats, UnivariateStats]
        x = extract(T, lines)
        isnothing(x) || push!(out, x)
    end
    return out
end

function Base.show(io::IO, t::TestCase)
    printstyled(io, "Test Case\n", color=:light_cyan, bold=true)
    for (k,v) in [
            "LocalFile" => t.localfile,
            "URL      " => t.url,
            "Procedure" => t.procedure,
        ]
        printstyled(io, "    $k: ")
        printstyled(io, v, color=:light_green)
        println(io)
    end
    for (k,v) in [
            # "File Format" => t.format,
            "Data Description" => t.data_description,
            "Model" => t.model,
            "Data" => summary(t.data)
        ]
        v2 = split(string(v), '\n')
        printstyled(io, "    $k: \n")
        for val in v2
            printstyled(io, "        $val\n", color=:light_green)
        end
    end
    printstyled(io, "    Certified Values:\n")
    for val in t.certified_values
        for row in split(string(val), '\n')
            printstyled(io, "        $row\n", color=:light_green)
        end
    end
end

#----------------------------------------------------# AnovaTable (anova, regression)
struct AnovaTable
    model::NamedTuple{(:df, :ss, :ms), Tuple{Int, Float64, Float64}}
    residual::NamedTuple{(:df, :ss, :ms), Tuple{Int, Float64, Float64}}
    f_statistic::Float64
end
function Base.show(io::IO, t::AnovaTable)
    m, r = t.model, t.residual
    d = DataFrame(source=["Model", "Residual"], df=[m.df, r.df], ss=[m.ss, r.ss], ms=[m.ms, r.ms])
    printstyled(io, "AnovaTable                           F = $(t.f_statistic)\n", bold=true)
    io2 = IOBuffer()
    show(io2, d; summary=false, eltypes=false, show_row_number=false)
    for row in split(String(take!(io2)), '\n')
        println(io, "    $row")
    end
end

function extract(::Type{AnovaTable}, lines)
    # anova
    i = findfirst(startswith("Between"), lines)
    if !isnothing(i)
        _, _, df, ss, ms, f = split(lines[i])
        _, _, df2, ss2, ms2 = split(lines[i+1])
        return AnovaTable(
            (df = Parsers.parse(Int, df), ss = Parsers.parse(Float64, ss), ms=Parsers.parse(Float64, ms)),
            (df = Parsers.parse(Int, df2), ss = Parsers.parse(Float64, ss2), ms=Parsers.parse(Float64, ms2)),
            Parsers.parse(Float64, f)
        )
    end
    # regression
    i = findfirst(x -> occursin("Certified Analysis of Variance Table", x), lines)
    if !isnothing(i)
        _, df, ss, ms, f = split(lines[i+5])
        _, df2, ss2, ms2 = split(lines[i+6])
        return AnovaTable(
            (df = Parsers.parse(Int, df), ss = Parsers.parse(Float64, ss), ms=Parsers.parse(Float64, ms)),
            (df = Parsers.parse(Int, df2), ss = Parsers.parse(Float64, ss2), ms=Parsers.parse(Float64, ms2)),
            Parsers.parse(Float64, f)
        )
    end
    return nothing
end

#----------------------------------------------------# Coefficients (nonlinear, regression)
struct Coefficients
    β::Vector{Float64}
    β_std_err::Vector{Float64}
end
function Base.show(io::IO, c::Coefficients)
    println(io, "Coefficients")
    d = DataFrame(β = c.β, se = c.β_std_err)
    io2 = IOBuffer()
    show(io2, d; summary=false, eltypes=false, show_row_number=false)
    for row in split(String(take!(io2)), '\n')
        println(io, "    $row")
    end
end

function extract(::Type{Coefficients}, lines)
    # regression
    i = findfirst(x -> occursin("Certified Regression Statistics", x), lines)
    if !isnothing(i)
        β, βse = Float64[], Float64[]
        j = 5
        line = lines[i + j]
        while startswith(lstrip(line), 'B')
            _, βi, βsei = split(line)
            push!(β, Parsers.parse(Float64, βi))
            push!(βse, Parsers.parse(Float64, βsei))
            j += 1
            line = lines[i + j]
        end
        return Coefficients(β, βse)
    end
    # TODO: MCMC
    # TODO: Nonlinear
    return nothing
end

#----------------------------------------------------# RSquared (anova, regression)
struct RSquared
    r2::Float64
end
Base.show(io::IO, r::RSquared) = print(io, "R-Squared: ", r.r2)
function extract(::Type{RSquared}, lines)
    i = findfirst(x -> occursin("R-Squared", x), lines)
    if !isnothing(i)
        line = lines[i]
        j = occursin("Certified", line) ? 3 : 2
        return RSquared(Parsers.parse(Float64, split(lines[i])[j]))
    end
    return nothing
end

#----------------------------------------------------# ResidualStandardError (anova, regression, nonlinear)
struct ResidualStandardError
    resid_std_err::Float64
end
Base.show(io::IO, r::ResidualStandardError) = print(io, "Residual Standard Error: ", r.resid_std_err)
function extract(::Type{ResidualStandardError}, lines)
    i = findfirst(x -> occursin("Residual", x), lines)
    if !isnothing(i)
        return ResidualStandardError(Parsers.parse(Float64, split(lines[i+1])[end]))
    end
    return nothing
end

#-----------------------------------------------------------------------------# MCMCStats
struct MCMCStats
    μ_posterior::NamedTuple{(:m, :s, :q), Tuple{Float64, Float64, Vector{Float64}}}
    σ_posterior::NamedTuple{(:m, :s, :q), Tuple{Float64, Float64, Vector{Float64}}}
end
function Base.show(io::IO, m::MCMCStats)
    io2 = IOContext(io, :compact=>true)
    println(io2, "MCMCStats")
    println(io2, "    μ Posterior: |", m.μ_posterior)
    println(io2, "    σ Posterior: |", m.σ_posterior)
end
function extract(::Type{MCMCStats}, lines)
    i = findfirst(x -> occursin("Posterior mean of mu", x), lines)
    if !isnothing(i)
        return MCMCStats(
            (
                m = Parsers.parse(Float64, split(lines[i])[end]),
                s = Parsers.parse(Float64, split(lines[i+1])[end]),
                q = [
                    Parsers.parse(Float64, split(lines[i+2])[end]),
                    Parsers.parse(Float64, split(lines[i+3])[end]),
                    Parsers.parse(Float64, split(lines[i+4])[end]),
                ]
            ),
            (
                m = Parsers.parse(Float64, split(lines[i+6])[end]),
                s = Parsers.parse(Float64, split(lines[i+7])[end]),
                q = [
                    Parsers.parse(Float64, split(lines[i+8])[end]),
                    Parsers.parse(Float64, split(lines[i+9])[end]),
                    Parsers.parse(Float64, split(lines[i+10])[end]),
                ]
            )
        )
    end
    return nothing
end

#----------------------------------------------------# UnivariateStats
struct UnivariateStats
    mean::Float64
    std::Float64
    autocor::Float64
    is_exact::Bool
    nobs::Int
end
function Base.show(io::IO, s::UnivariateStats)
    print(io, "UnivariateStats: m=$(s.mean), s=$(s.std), a=$(s.autocor), n=$(s.nobs) | exact=$(s.is_exact)")
end
function extract(::Type{UnivariateStats}, lines)
    i = findfirst(x -> occursin("Sample Mean", x), lines)
    if !isnothing(i)
        s1 = split(lines[i])
        j = s1[end] == "(exact)" ? 1 : 0
        m = Parsers.parse(Float64, s1[end-j])
        s = Parsers.parse(Float64, split(lines[i+1])[end-j])
        a = Parsers.parse(Float64, split(lines[i+2])[end-j])
        nobs = Parsers.parse(Int, split(lines[i+4])[end])
        return UnivariateStats(m, s, a, j == 1, nobs)
    end
    return nothing
end

end


text:
# Use
#
#     DOCUMENTER_DEBUG=true julia --color=yes make.jl local [nonstrict] [fixdoctests]
#
# for local builds.

using Documenter
using Literate
using BAT

import HDF5

function fix_literate_output(content)
    content = replace(content, "EditURL = \"@__REPO_ROOT_URL__/\"" => "")
    return content
end

gen_content_dir = joinpath(@__DIR__, "src")
tutorial_src = joinpath(@__DIR__, "src", "tutorial_lit.jl")
Literate.markdown(tutorial_src, gen_content_dir, name = "tutorial", documenter = true, credit = true, postprocess = fix_literate_output)
#Literate.markdown(tutorial_src, gen_content_dir, name = "tutorial", codefence = "```@repl tutorial" => "```", documenter = true, credit = true)
Literate.notebook(tutorial_src, gen_content_dir, execute = false, name = "bat_tutorial", documenter = true, credit = true)
Literate.script(tutorial_src, gen_content_dir, keep_comments = false, name = "bat_tutorial", documenter = true, credit = false)

makedocs(
    sitename = "BAT",
    modules = [BAT],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://bat.github.io/BAT.jl/stable/"
    ),
    authors = "The BAT development team",
    pages=[
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Tutorial" => "tutorial.md",
        "API Documentation" => "stable_api.md",
        "Experimental Features" => "experimental_api.md",
        "Internal API" => "internal_api.md",
        "Developer instructions" => "developing.md",
        "License" => "license.md",
    ],
    doctest = ("fixdoctests" in ARGS) ? :fix : true,
    linkcheck = ("linkcheck" in ARGS),
    strict = !("nonstrict" in ARGS),
)

deploydocs(
    repo = "github.com/bat/BAT.jl.git",
    forcepush = true
)


text:
# ## Data

y = [-1, 1]

# ## Model

import Random

using Turing

Random.seed!(1)

@model function m8_2(y)
    α ~ Flat() ## improper prior with pobability one everywhere
    σ ~ FlatPos(0.0) ## improper prior with probability one everywhere above 0.0

    y ~ Normal(α, σ)
end;

# ## Output

chns = sample(m8_2(y), NUTS(), 1000)

# \defaultoutput{}


text:
# abstract types -
abstract type AbstractLatticeModel end
abstract type AbstractAssetModel end
abstract type AbstractDerivativeContractModel <: AbstractAssetModel end

# option contract types -
mutable struct PutContractModel <: AbstractDerivativeContractModel

    # data -
    ticker::String
    strike_price::Float64
    expiration_date::Date
    premium::Float64
    current_price::Float64
    direction::Int64
    number_of_contracts::Int64

    # constructor -
    PutContractModel() = new()
end

mutable struct CallContractModel <: AbstractDerivativeContractModel

    # data -
    ticker::String
    strike_price::Float64
    expiration_date::Date
    premium::Float64
    current_price::Float64
    direction::Int64
    number_of_contracts::Int64

    # constructor -
    CallContractModel() = new()
end

mutable struct EquityModel <: AbstractAssetModel

    # data -
    ticker::String
    purchase_price::Float64
    current_price::Float64
    direction::Int64
    number_of_shares::Int64

    # constructor -
    EquityModel() = new()
end

mutable struct CRRLatticeModel <: AbstractLatticeModel

    # data -
    number_of_levels::Int64
    branch_factor::Int64
    connectivity::Array{Int64,2}

    σ::Float64
    p::Float64
    μ::Float64
    u::Float64
    d::Float64
    ΔT::Float64
    T::Float64
    Sₒ::Float64

    data::Array{Float64,2}

    # constructor -
    CRRLatticeModel() = new()
end

mutable struct PQContractPremiumLatticePoint

    # data -
    s::Int64
    d::Int64
    i::Int64
    j::Int64

    # constructor -
    PQContractPremiumLatticePoint() = new()
end

mutable struct CRRJITContractPremiumLatticeModel

    # data -
    underlying::Array{Float64,1}
    strike::Union{Array{Float64,1}, Array{Float64,2}}
    iv::Array{Float64,1}
    dte::Array{Float64,1}
    number_of_levels::Int64
    risk_free_rate::Float64
    contractType::Type{T} where {T <: AbstractDerivativeContractModel}

    # constructor -
    CRRJITContractPremiumLatticeModel() = new()
end

mutable struct CRRContractPremiumLatticeModel

    # data -
    grid::Dict{NamedTuple,Union{Float64, Array{Float64,1}}}
    underlying::Array{Float64,1}
    strike::Union{Array{Float64,1}, Array{Float64,2}}
    iv::Array{Float64,1}
    dte::Array{Float64,1}
    number_of_levels::Int64
    risk_free_rate::Float64
    contractType::Type{T} where {T <: AbstractDerivativeContractModel}

    # constructor -
    CRRContractPremiumLatticeModel() = new()
end


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using DataFrames
using CSV

function plant_outputs_report(solution)::DataFrame
    df = DataFrame()
    df."plant type" = String[]
    df."location name" = String[]
    df."year" = Int[]
    df."product name" = String[]
    df."amount produced (tonne)" = Float64[]
    df."amount sent (tonne)" = Float64[]
    df."amount disposed (tonne)" = Float64[]
    df."disposal cost (\$)" = Float64[]
    T = length(solution["Energy"]["Plants (GJ)"])
    for (plant_name, plant_dict) in solution["Plants"]
        for (location_name, location_dict) in plant_dict
            for (product_name, amount_produced) in location_dict["Total output"]
                send_dict = location_dict["Output"]["Send"]
                disposal_dict = location_dict["Output"]["Dispose"]

                sent = zeros(T)
                if product_name in keys(send_dict)
                    for (dst_plant_name, dst_plant_dict) in send_dict[product_name]
                        for (dst_location_name, dst_location_dict) in dst_plant_dict
                            sent += dst_location_dict["Amount (tonne)"]
                        end
                    end
                end
                sent = round.(sent, digits = 2)

                disposal_amount = zeros(T)
                disposal_cost = zeros(T)
                if product_name in keys(disposal_dict)
                    disposal_amount += disposal_dict[product_name]["Amount (tonne)"]
                    disposal_cost += disposal_dict[product_name]["Cost (\$)"]
                end
                disposal_amount = round.(disposal_amount, digits = 2)
                disposal_cost = round.(disposal_cost, digits = 2)

                for year = 1:T
                    push!(
                        df,
                        [
                            plant_name,
                            location_name,
                            year,
                            product_name,
                            round(amount_produced[year], digits = 2),
                            sent[year],
                            disposal_amount[year],
                            disposal_cost[year],
                        ],
                    )
                end
            end
        end
    end
    return df
end

write_plant_outputs_report(solution, filename) =
    CSV.write(filename, plant_outputs_report(solution))

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
"""
NOAA provides nighttime lights images as tif files. They can be opened as 2D arrays using the ```load_img``` function.
```julia
load_img("example.tif")
```
"""
function load_img(filepath)
    img = ArchGDAL.read((ArchGDAL.read(filepath)),1)
    img_trans = img' #tif loaded using ArchGDAL needs to be transposed
    return Array{Union{Missing, Float16}, 2}(img_trans)
end

"""
NOAA provides nighttime lights images as tif files. They can be opened as 2D arrays using the ```load_img``` function. The top-left and bottom-right parameters can be used to crop the image.
```julia
load_img("example.tif")
```
"""
function load_img(filepath, top_left, bottom_right)
    img = ArchGDAL.readraster(filepath)
    img = img[top_left[2]+1:bottom_right[2], top_left[1]+1:bottom_right[1], 1]
    GC.gc()
    return Array{Union{Missing, Float16}, 2}(img')
end

"""
NOAA provides nighttime lights images as tif files. They can be opened as 2D arrays using the load_img function. The top-left and bottom-right parameters can be used to crop the image.
```julia
load_img("example.tif")
```
"""
function load_img(filepath, top_left::Coordinate, bottom_right::Coordinate, geometry::CoordinateSystem)
    top_left = coordinate_to_image(geometry, top_left)
    bottom_right = coordinate_to_image(geometry, bottom_right)
    img = load_img(filepath, top_left, bottom_right)
    return img
end

"""
Images in the form of 2D arrays can be saved as tif files.
```julia
save_img("example.tif", img)
```
"""
function save_img(filepath,img)
    img = Float64.(img)
    GeoArrays.write!(filepath,GeoArray(img'))
end

"""
NOAA provides images for each month since April 2012. Images of the same place taken at different times can be stacked together to make a 3D array representating a panel data. In this package, we refer to such arrays as datacubes. JLD files containing 3D arrays can be loaded using the ```load_datacube``` function.
```julia
load_datacube("example.jld")
```
"""
function load_datacube(datacube_path)
    datacube = load(datacube_path)["data"]
    return datacube
end

"""
3D arrays can be saved as JLD files.
```julia
save_datacube("example.jld", datacube)
```
"""
function save_datacube(filepath, datacube)
    save(filepath, "data", datacube)
end

"""
Loads all images (tif files) in a folder and generates a datacube. The function prints the file names to you the order in which they are loaded.
```julia
make_datacube("~/Downloads/ntl_images")
```
"""
function make_datacube(folder_path, display_names = false)
        files = readdir(folder_path)
        if display_names == true
            print(files)
        end
        paths = folder_path .* "/".* files
        datacube = []
        for path in paths
            push!(datacube, load_img(path))
        end
        return Array{Union{Missing, Float16}, 3}(cat(datacube...,dims = 3))
end

"""
Loads all images (tif files) in a folder and generates a datacube. The function prints the file names to you the order in which they are loaded. The ```top_left``` and the ```bottom_right``` parameters can be used to crop the datacube.
```julia
make_datacube("~/Downloads/ntl_images")
```
"""
function make_datacube(folder_path, top_left, bottom_right, display_names = false)
        files = readdir(folder_path)
        if display_names == true
            display(files)
        end
        paths = folder_path .* "/".* files
        datacube = []
        for path in paths
            push!(datacube, load_img(path, top_left, bottom_right))
        end
        return Array{Union{Missing, Float16}, 3}(cat(datacube...,dims = 3))
end

"""
Loads all images (tif files) in a folder and generates a datacube. The function prints the file names to you the order in which they are loaded. The ```top_left``` and the ```bottom_right``` parameters can be used to crop the datacube.
```julia
make_datacube("~/Downloads/ntl_images")
```
"""
function make_datacube(folder_path, top_left::Coordinate, bottom_right::Coordinate, geometry::CoordinateSystem, display_names = false)
        files = readdir(folder_path)
        if display_names == true
            display(files)
        end
        paths = folder_path .* "/".* files
        datacube = []
        for path in paths
            push!(datacube, load_img(path, top_left, bottom_right, geometry))
        end
        return Array{Union{Missing, Float16}, 3}(cat(datacube...,dims = 3))
end