text:
#entry point to acquire the parameters
function getparameters!(parameterfile::String=PARAMETER_FILE)

  #this will contain the parameters
  workbook = XLSX.readxlsx(parameterfile)
  (length(PARAM) > 0) && clearcollection(PARAM)

  #iterate across sheets
  for sheet::XLSX.Worksheet ∈ workbook.workbook.sheets
    startswith("$(sheet.name)", "dep_") && continue #allows for depreciaiton of worksheets
    getparameters(sheet, PARAM)
  end

  close(workbook)

  logparameters()

  #PARAM[] = param

  return nothing
end

#extracts and parses the parameters out of each workbook
function getparameters(sheet::XLSX.Worksheet, param::AbstractDict{Symbol, Any})

  #create a dataframe
  paramdf::DataFrame = DataFrame(XLSX.gettable(sheet)...)

  #println(describe(paramdf))

  #easier to access the vectors as a dictionary
  #paramdfindex = Dict(r.x2=>r.x1 for r ∈ eachrow(paramdf))

  #get the types
  typestrings::Vector{String} = paramdf.type

  #get the default values
  valuestrings::Vector{String} = paramdf.default

  #get the labels
  parameterlabels = (Symbol).(paramdf.parameter)

  #now add in the overrides
  for (i, s) ∈ enumerate(paramdf.override)
    if !ismissing(s)
      try
        valuestrings[i] = s
      catch err
        error("Failed to load override. parse failed.
          Label: $(parameterlabels[i]) Type: $(typestrings[i]) Val: $s
          error: $err
          trace:\n$(stacktrace(catch_backtrace()))")
      end
    end
  end

  #now parse the parameters to their correct type
  for (i, s) ∈ enumerate(valuestrings)
    try

      #each parameter needs a unique key
      haskey(param, parameterlabels[i]) && error(
        "Multiple parameters found with same key.")

      #core parsing logic
      T::Type = parseparam(Type, typestrings[i])
      param[parameterlabels[i]] = parseparam(T, s)

    catch err #throw a meaningful error following failure
      throw("Parameter parse failed.
        Label: $(parameterlabels[i]) Type: $(typestrings[i]) Val: $s
        error: $err")
    end
  end
end

#try to write out all the parameter info used
#WARNING: constants need to be actively added to this list, so try to avoid them
function formatparameters()
  b = IOBuffer()
  write(b, "Capacity log from time: $(Dates.format(now(),"yyyymmdd_HHMM"))")
  write(b, "\n*****************\n")
  write(b, formatdict(PARAM, title="PARAM dict:"))
  write(b, "\n*****************\n")
  write(b, formatdict(REGRESSION_METHOD, title="REGRESSION_METHOD dict:"))
  write(b, "\n*****************\nConstants: \n")
  write(b, "RANDOM_SEED_VALUE = $RANDOM_SEED_VALUE\n")
  write(b, "PARAMETER_FILE = $(PARAMETER_FILE)\n")
  write(b, "PARALLEL = $(PARALLEL)\n")
  write(b, "LARGE_VAL = $(LARGE_VAL)\n")
  write(b, "LARGE_VAL32 = $(LARGE_VAL32)\n")
  write(b, "SMALL_VAL = $(SMALL_VAL)\n")
  write(b, "SMALL_VAL32 = $(SMALL_VAL32)\n")

  return String(take!(b))
end

function logparameters(;logpath = "$(PARAM[:testpath])\\log",
    logname = "cap-$(Dates.format(now(),"yyyymmdd_HHMM"))")
  paramlog = formatparameters()
  open("$logpath\\$logname.txt", "w+") do f
    write(f, paramlog)
  end

  return nothing
end

#parsing scenarios
parseparam(::Type{Nothing}, args...) = nothing
parseparam(::Type{T}, s::String) where T<:Number = something(tryparse(T,s), parseparamalways(T,s))
parseparam(::Type{Type}, s::String) = eval(Meta.parse("""$(s)"""))
parseparam(::Type{String}, s::String) = s
parseparam(::Type{Union{String,Nothing}}, s::String) = s ≡"nothing" ? nothing : s

parseparam(::Type{Date}, i::Int, dateformat::DateFormat) = parseparam(Date, "$i", dateformat)
parseparam(::Type{DateFormat}, s::String) = DateFormat(s)
parseparam(::Type{T}, s::String) where T<:DatePeriod = eval(Meta.parse("""$T($(s))"""))
parseparam(::Type{Date}, s::String, dateformat::DateFormat) = Dates.Date(s, dateformat)

#fallback
parseparam(::Type{T}, s::String) where T = parseparamalways(T,s)
function parseparamalways(::Type{T}, s::String) where T
  x::T = eval(Meta.parse("""($(s))"""))
  return x
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module ThermalTraj

using JuMP
using NLopt
using JSON

export printTest
export heading2XY
export initTrajModel!
export solveTraj!

function printTest()
  println("Testing print from module ThermalTraj")
end

function heading2XY(h)
    e = -(h-90.0)#*pi/180.0
    if e <= -180.0
        e += 360.0
    end
    return e*pi/180.0
end

# function to setup the optimization model, but not optimize the traj
function initTrajModel!(m, x0, nodes, l_x, b_u)
    b_dx = 25.9 #80.0/3.6
    b_dy = 0
    dt = 0.6;
    #println("Constant forward body speed of [", b_dx, ", ", b_dy,"] [m/s]")

    umax = 40.0*pi/180.0;

    # Set the objective function
    @setNLObjective(m, Min,  sum{(l_x[1,i])^2+(l_x[2,i])^2, i = 1:nodes}
                    + 100.0*sum{(b_u[i])^2,i=1:floor(nodes*4.0/5.0)}
                    - 50.0 *sum{ l_x[4,i-3]*l_x[4,i] ,i=4:nodes})

    x = linspace(x0[1],0.,nodes);
    y = linspace(x0[2],0.,nodes);

    phi0 = -20.0*pi/180.0
    if x0[4] > 0
    	phi0 = 20.0*pi/180.0
    end

    for n in 1:nodes
        #initialize variables
        setValue(l_x[1,n],x[n])
        setValue(l_x[2,n],y[n])
        setValue(l_x[3,n],0.0)
        setValue(l_x[4,n],phi0)
        setValue(l_x[5,n],0.0)
        setValue(b_u[n],0.0)
    end

 	expdt = exp(dt)

	#a34 = 0.06039
	#a35 = 0.23181
	#a44 = 0.69768
	#a54 = 0.34478
	#a55 = 0.58275
	#b32 = 0.0091949
	#b42 = 0.42167
	#b52 = 0.090857

	#a34 = 0.092057
	#a35 = 0.19781
	#a44 = 0.74082
	#a54 = 0.50137
	#a55 = 0.40657
	#b32 = 0.011964
	#b42 = 0.36043
	#b52 = 0.11599

	a34 = 0.082852
	a35 = 0.17803
	a44 = 0.74082
	a54 = 0.50137
	a55 = 0.40657
	b32 = 0.010768
	b42 = 0.36043
	b52 = 0.11599

	# initial condition constraints on dynamics
    @addNLConstraint(m, b_u[1] <= umax)
    @addNLConstraint(m, -b_u[1] <= umax)
    @addNLConstraint(m, l_x[4,1] <= umax)
    @addNLConstraint(m, -l_x[4,1] <= umax)
    @addNLConstraint(m, l_x[1,1] == x0[1] + (expdt*cos(x0[3]) - cos(x0[3]))*b_dx  )
    @addNLConstraint(m, l_x[2,1] == x0[2] + (expdt*sin(x0[3]) - sin(x0[3]))*b_dx  )
    @addNLConstraint(m, l_x[3,1] == x0[3] + a34*x0[4] + a35*x0[5] + b32*b_u[1]  )
    @addNLConstraint(m, l_x[4,1] == x0[4]*a44 + b42*b_u[1]  )
    @addNLConstraint(m, l_x[5,1] == x0[5]*a55 + a54*x0[4] + b52*b_u[1]  )

    for n in 1:(nodes-1)
        # dynamics constraints
        @addNLConstraint(m, l_x[1,n+1] == l_x[1,n] + (expdt*cos(l_x[3,n]) - cos(l_x[3,n]))*b_dx  )
		@addNLConstraint(m, l_x[2,n+1] == l_x[2,n] + (expdt*sin(l_x[3,n]) - sin(l_x[3,n]))*b_dx  )
		@addNLConstraint(m, l_x[3,n+1] == l_x[3,n] + a34*l_x[4,n] + a35*l_x[5,n] + b32*b_u[n+1]  )
		@addNLConstraint(m, l_x[4,n+1] == l_x[4,n]*a44 + b42*b_u[n+1]  )
		@addNLConstraint(m, l_x[5,n+1] == l_x[5,n]*a55 + a54*l_x[4,n] + b52*b_u[n+1]  )
        @addNLConstraint(m, b_u[n+1] <= umax)
        @addNLConstraint(m, -b_u[n+1] <= umax)
        @addNLConstraint(m, l_x[4,n+1] <= umax)
        @addNLConstraint(m, -l_x[4,n+1] <= umax)
    end
end

function solveTraj!(x0,u,x)
    # Setup the optimization scheme
    m = Model(solver=NLoptSolver(algorithm=:LD_SLSQP))
    # state variables for for direct transcription
    nodes = size(u)[1]
    @defVar(m, l_x[1:5,1:nodes])
    # input decision variables -- u[1] is actually u[0]
    @defVar(m, b_u[1:nodes])

    initTrajModel!(m, x0, nodes, l_x, b_u)
    # Test the solving capability

    toc = @elapsed status = solve(m)

    # get optimal values in general numeric array
    for i in 1:nodes, j in 1:5
        x[j,i] = getValue(l_x[j,i])
        u[i] = getValue(b_u[i])
    end

    println("TrajOpt compute time: ", toc)
end

end

ext: jl

max_stars_repo_name: dehann/ThermalSoaring

max_issues_repo_name: dehann/ThermalSoaring

max_forks_repo_name: dehann/ThermalSoaring


text:
module FiniteDifferenceSchemes
    using LinearAlgebra

    export Weno5, Linear
    export linear_scheme, weno5_js
    export flux_positive_linear
    export flux_positive_weno5
    export flux_negative_function

    abstract type AbstractSchemes{T1<:Integer, T2<:AbstractFloat} end

    struct Weno5{T1<:Integer, T2<:AbstractFloat} <: AbstractSchemes{T1, T2}
        offset_min::T1
        offset_max::T1
        epsilon::T2
    end # weno5

    struct Linear{T1<:Integer, T2<:AbstractFloat} <: AbstractSchemes{T1, T2}
        offset_min::T1
        offset_max::T1
        flux_coefficients::Array{T2, 1}
    end # linear

    include("linear_scheme_coefficients.jl")
    include("inverted_vector.jl")
    include("flux_negative.jl")
    include("linear.jl")
    include("weno.jl")
end # module


text:
#=
Description

US presidents serve four year terms, with most presidents serving one or two terms.
Unless a president dies in office, they then live after leaving office.
This challenge, then, given a list of presidents and their dates of birth and dates
of death, is to figure out what year the most presidents - future, present, or
previous - were alive.

Challenge Input

In a CSV input of presidential birthdates and death dates. Find what year
in which the most number of people who would serve, are serving, or have served
as presidents. The answer might be multiple years, or only a partial year.

Challenge Output

Any of the following years is valid: 1822, 1823, 1824, 1825, 1826, 1831, 1833,
1834, 1835, 1836, 1837, 1838, 1839, 1840, 1841, 1843, 1844, 1845 (each year had
18 presidents, current, former, or to be, alive).
=#

using DataFrames

type President
    birth_year::Integer
    death_year::Integer
end

function main( )
    local data = readtable( "presidents.csv", header=true )
    local n_rows::Integer = n_cols::Integer = time::Integer = count::Integer = index::Integer = 0
    local all_presidents = President[]

    ( n_rows, n_cols ) = size( data )

    for i in 1:n_rows
        birth_year::Integer = parse( Int,( match( r"\d{4}", data[ i, 2 ] ) ).match )
        death_year::Integer = ( !is( data[ i, 4 ], NA ) ) ? parse( Int,( match( r"\d{4}", data[ i, 4 ] ) ).match ) : 0
        push!( all_presidents, President( birth_year, death_year ) )
    end

    isless( p1::President, p2::President ) = ( p1.birth_year < p2.birth_year )
    sort!( all_presidents, lt=isless )
    for i in 1:length( all_presidents )-1
        pivot::President = all_presidents[ i ]

        for j in i+1:length( all_presidents )
            if all_presidents[ j ].birth_year > pivot.death_year
                break
            else
                count += 1
            end
        end

        if count > time
            time = count
            index = i
        end

        count = 0
    end

    println( all_presidents[ index+time ].birth_year )
end

main( )

max_stars_repo_name: farmfirst/dailyprogrammer

max_issues_repo_name: farmfirst/dailyprogrammer

max_forks_repo_name: farmfirst/dailyprogrammer


text:
## AbstractArray wrapper for collection of variates
## Similar to https://github.com/JuliaLang/julia/pull/32310 - replace with EachSlice?
struct EachVariate{V,P,A,T,N} <: AbstractArray{T,N}
    parent::P
    axes::A
end

function EachVariate{V}(x::AbstractArray{<:Real,M}) where {V,M}
    ax = ntuple(i -> axes(x, i + V), Val(M - V))
    T = typeof(view(x, ntuple(i -> i <= V ? Colon() : firstindex(x, i), Val(M))...))
    return EachVariate{V,typeof(x),typeof(ax),T,M-V}(x, ax)
end

Base.IteratorSize(::Type{EachVariate{V,P,A,T,N}}) where {V,P,A,T,N} = Base.HasShape{N}()

Base.axes(x::EachVariate) = x.axes

Base.size(x::EachVariate) = map(length, x.axes)
Base.size(x::EachVariate, d::Int) = 1 <= ndims(x) ? length(axes(x)[d]) : 1

# We don't need `setindex!` (currently), therefore only `getindex` is implemented
Base.@propagate_inbounds function Base.getindex(
    x::EachVariate{V,P,A,T,N}, I::Vararg{Int,N},
) where {V,P,A,T,N}
    return view(x.parent, ntuple(_ -> Colon(), Val(V))..., I...)
end

# optimization for univariate distributions
eachvariate(x::AbstractArray{<:Real}, ::Type{Univariate}) = x
function eachvariate(x::AbstractArray{<:Real}, ::Type{ArrayLikeVariate{N}}) where {N}
    return EachVariate{N}(x)
end


text:
#=
  Set covering in Julia ConstraintSolver.jl

  Placing of firestations, from Winston "Operations Research", page 486

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK
const CS = ConstraintSolver
include("constraints_utils.jl")

function set_covering()

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   # "all_solutions"=> all_solutions,
                                                            "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS, # <-
                                                            # "traverse_strategy"=>:DBFS,

                                                            "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            # "branch_split"=>:InHalf, # <-

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>16,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                        ))
    min_distance = 15 # minimum distance
    # distances between the cities
    distance = [[0,10,20,30,30,20],
               [10,0,25,35,20,10],
               [20,25,0,15,30,20],
               [30,35,15,0,15,25],
               [30,20,30,15,0,14],
               [20,10,20,25,14,0]]
    num_cities = length(distance)

    # where to place the fire stations: 1 if placed in this city.
    @variable(model, x[1:num_cities], Bin)
    @variable(model, 0 <= z <= num_cities, Int)

    # calculate the number of covered fire stations
    for j in 1:num_cities
        @constraint(model,sum([x[i] for i in 1:num_cities if distance[i][j] <= min_distance]) >= 1)
    end

    # objective: minimize the number of fire stations
    @constraint(model, z == sum(x))

    @objective(model, Min, z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                xx = convert.(Integer,JuMP.value.(x; result=sol))
                zz = convert.(Integer,JuMP.value.(z; result=sol))
                println("x:$xx z:$zz")

            end
        end
    else
        println("status:$status")
    end

    return status
end

@time set_covering()

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
push!(LOAD_PATH, "../src/")
ENV["GKSwstype"] = "100"

using DifferentialEquations
using DelimitedFiles
using Documenter
using Literate
using Plots
using AnalogDataAssimilation

# generate examples
OUTPUT = joinpath(@__DIR__, "src", "generated")

# examples = filter(x -> occursin(r".jl", x), map(relpath, readdir(joinpath(@__DIR__, "..", "examples"))))

examples = String[]
push!(examples, "catalog.jl")
push!(examples, "model_forecasting.jl")
push!(examples, "analog_forecasting.jl")
push!(examples, "lorenz63.jl")

example_pages = Any[]

for example in examples
    EXAMPLE = joinpath(@__DIR__, "examples", example)
    @show page = string("generated/", example[1:end-3], ".md")

    open(`head -n1 $EXAMPLE`) do io
        title = join(readdlm(io)[3:end], " ")
        push!(example_pages, title => page)
    end

    Literate.markdown(EXAMPLE, OUTPUT)
    #Literate.notebook(EXAMPLE, OUTPUT)
    #Literate.script(EXAMPLE, OUTPUT)
end

@show example_pages

pages = Any[
    "Home"=>"index.md",
    "Catalog"=>"catalog.md",
    "Data Assimilation"=>"data_assimilation.md",
    "State Space"=>"state-space.md",
    "Models"=>"models.md",
    "Ensemble Kalman filters"=>"ensemble_kalman_filters.md",
    "Ensemble Kalman smoothers"=>"ensemble_kalman_smoothers.md",
    "Forecasting"=>"forecasting.md",
    "Particle filters"=>"particle_filters.md",
    "Time Series"=>"time-series.md",
    "Utilities"=>"utils.md",
    "Examples"=>example_pages,
]

@show pages

makedocs(
    modules = [AnalogDataAssimilation],
    sitename = "AnalogDataAssimilation.jl",
    doctest = true,
    authors = "Pierre Navaro",
    format = Documenter.HTML(),
    pages = pages,
)

deploydocs(repo = "github.com/pnavaro/AnalogDataAssimilation.jl.git", devbranch = "main")


text:
using ReferenceSequences
using Bio.Seq

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "Conversion" begin
    @test isa(convert(ReferenceSequence, dna""), ReferenceSequence)
    @test isa(convert(ReferenceSequence, dna"ACGTN"), ReferenceSequence)
    @test isa(convert(ReferenceSequence, dna"ACGTN"^50), ReferenceSequence)
    @test_throws ArgumentError convert(ReferenceSequence, dna"ACGRT")

    @test isa(convert(ReferenceSequence, ""), ReferenceSequence)
    @test isa(convert(ReferenceSequence, "ACGTN"), ReferenceSequence)

    refseq = ReferenceSequence(dna"ACGTN")
    @test isa(convert(DNASequence, refseq), DNASequence)
    @test isa(convert(ASCIIString, refseq), ASCIIString)
end

@testset "Basic Operations" begin
    seq = ReferenceSequence(dna"")
    @test length(seq) === endof(seq) === 0
    @test isempty(seq)
    @test seq == dna""
    @test cmp(seq, seq) == 0
    @test cmp(seq, dna"") == 0
    @test cmp(dna"", seq) == 0
    @test !isless(seq, seq)
    @test isless(seq, dna"A")
    @test_throws BoundsError seq[-1]
    @test_throws BoundsError seq[0]
    @test_throws BoundsError seq[1]
    @test_throws BoundsError seq[2:3]
    @test ReferenceSequence("") == seq

    seq = ReferenceSequence(dna"ACGTN")
    @test length(seq) === endof(seq) === 5
    @test !isempty(seq)
    @test seq[1] === DNA_A
    @test seq[2] === DNA_C
    @test seq[3] === DNA_G
    @test seq[4] === DNA_T
    @test seq[5] === DNA_N
    @test seq[2:3] == dna"CG"
    @test seq[3:5] == dna"GTN"
    @test seq[5:4] == dna""
    @test_throws BoundsError seq[0]
    @test_throws BoundsError seq[6]
    @test_throws BoundsError seq[0:2]
    @test_throws BoundsError seq[5:6]
    @test_throws MethodError seq[1] = DNA_A
    @test collect(seq) == [DNA_A, DNA_C, DNA_G, DNA_T, DNA_N]
    @test seq == seq
    @test seq == dna"ACGTN"
    @test dna"ACGTN" == seq
    @test cmp(seq, seq) == 0
    @test ReferenceSequence("ACGTN") == seq

    seq = ReferenceSequence(dna"ACGNNNANG")
    @test findnext(seq, DNA_A, 1) == 1
    @test findnext(seq, DNA_C, 1) == 2
    @test findnext(seq, DNA_G, 1) == 3
    @test findnext(seq, DNA_T, 1) == 0
    @test findnext(seq, DNA_N, 1) == 4
    @test findnext(seq, DNA_A, 2) == 7
    @test findnext(seq, DNA_N, 8) == 8

    @test_throws BoundsError findnext(seq, DNA_A, 0)

    @test findprev(seq, DNA_A, 9) == 7
    @test findprev(seq, DNA_C, 9) == 2
    @test findprev(seq, DNA_T, 9) == 0
    @test findprev(seq, DNA_N, 9) == 8
    @test findprev(seq, DNA_N, 5) == 5
    @test findprev(seq, DNA_N, 3) == 0

    @test_throws BoundsError findnext(seq, DNA_A, 10)

    @test findfirst(seq, DNA_A) == 1
    @test findfirst(seq, DNA_N) == 4
    @test findlast(seq, DNA_A) == 7
    @test findlast(seq, DNA_N) == 8
end

@testset "Long Sequences" begin
    seq = ReferenceSequence(dna"ACGT"^10000)
    @test length(seq) == 4 * 10000
    @test seq == dna"ACGT"^10000
    @test_throws BoundsError seq[0]
    @test_throws BoundsError seq[4 * 10000 + 1]
    @test ReferenceSequence("ACGT"^10000) == seq

    seq = ReferenceSequence(
        dna"NNNN"^1000 *
        dna"ACGT"^1000 *
        dna"NNNN"^1000 *
        dna"TGCA"^1000 *
        dna"NNNN"^1000
    )
    @test length(seq) == 4 * 1000 * 5
    @test seq == (
        dna"NNNN"^1000 *
        dna"ACGT"^1000 *
        dna"NNNN"^1000 *
        dna"TGCA"^1000 *
        dna"NNNN"^1000)
    @test ReferenceSequence(
        "NNNN"^1000 *
        "ACGT"^1000 *
        "NNNN"^1000 *
        "TGCA"^1000 *
        "NNNN"^1000) == seq
end


text:
module HyperSparseMatrices
using SparseArrays
using StorageOrders
export HyperSparseMatrix, HyperSparseCSC, HyperSparseCSR

struct HyperSparseMatrix{O, Tv, Tf, Ti<:Integer} <: AbstractSparseMatrix{Tv, Ti}
    # Comments here reflect column major ordering.
    # To understand as DCSR simply swap references to columns with rows and vice versa.
    vlen::Int # m in CSC, n in CSR. This is the length of the vectors. In DCSC this is thus the number of rows.
    vdim::Int # n in CSC, m in CSR. This is the number of vectors being stored. in DCSC this is the number of columns.

    p::Vector{Ti} # The pointers into i/nzval. The row indices found in the k'th stored column are
    # found in idx[p[k], p[k+1]-1]

    # If column (row) j has stored entries, then j = h[k] for some k.
    # j is the k'th stored column.
    h::Vector{Ti}
    idx::Vector{Ti} # the stored row indices.
    v::Vector{Tv} # the coefficients of stored indices in the matrix
    fill::Tf # the fill value. A[i,j] == fill if (i,j) is compressed out.
end

StorageOrders.storageorder(::HyperSparseMatrix{O}) where {O} = O

const HyperSparseCSC{Tv, Tf, Ti} = HyperSparseMatrix{ColMajor(), Tv, Tf, Ti}
const HyperSparseCSR{Tv, Tf, Ti} = HyperSparseMatrix{RowMajor(), Tv, Tf, Ti}

HyperSparseCSC(vlen, vdim, p::Ti, h::Ti, idx::Ti, v::Tv, fill::Tf) where {Ti, Tv, Tf} = HyperSparseCSC{Tv, Tf, Ti}(vlen, dvim, p, h, idx, v, fill)
HyperSparseCSR(vlen, vdim, p::Ti, h::Ti, idx::Ti, v::Tv, fill::Tf) where {Ti, Tv, Tf} = HyperSparseCSR{Tv, Tf, Ti}(vlen, dvim, p, h, idx, v, fill)

Base.size(A::HyperSparseCSC) = (A.vlen, A.vdim)
Base.size(A::HyperSparseCSR) = (A.vdim, A.vlen)

SparseArrays.nnz(A::HyperSparseMatrix) = length(A.v)
SparseArrays.nonzeros(A::HyperSparseMatrix) = A.v

nvec(A::HyperSparseCSC) = size(A, 2)
nvec(A::HyperSparseCSR) = size(A, 1)

# indexing adapted from SS:GrB and SparseArrays:
function Base.getindex(A::HyperSparseMatrix{O, Tv, Tf, Ti}, row::Integer, col::Integer) where {O, Tv, Tf, Ti}
    # A HyperSparseMatrix stores sparse vectors.
    # Since we support both ColMajor and RowMajor i and j don't mean what they usually do.
    # Instead j refers to the vector
    # and i refers a scalar within that vector.
    if O === ColMajor() # For column major this is as it is for SparseMatrixCSC
        i = row
        j = col
    elseif O === RowMajor() # For row major, j is a row vector (which in this case (DCSC), may be compressed out entirely).
        i = col
        j = row
    end
    @boundscheck checkbounds(A, row, col) # the overloaded size should take care of doing this correctly.
    nnz(A) == o && return A.fill # no values, return fill.

    # A.h contains a list of vectors that were not compressed out.
    # We determine the idx in A.h of j (or the next greater idx if j ∉ A.h)
    idx = searchsortedfirst(A.h, j, 1, nvec(A), Forward)

    # if this is true then we have a vector, and this just becomes the same as indexing in SparseMatrixCSC
    # If we don't, we just return the fill value.
    j == A.h[idx] || return A.fill

    # this is all adapted directly from sparsematrix.jl ∈ SparseArrays.jl
    r1 = A.p[idx]
    r2 = A.p[idx + 1] - 1
    (r1 > r2) && return A.fill
    r1 = searchsortedfirst(A.idx, i, r1, r2, Forward)
    return ((r1 > r2) || (A.idx[r1] != i)) ? A.fill : A.v[r1]
end

#TODO:
# indexing
# printing (braille?)
# implement a good subset of the AbstractArray and AbstractSparseMatrix interface.
end


text:
using Pkg
using Test
using StaticArrays
using AbstractPlotting
using AbstractPlotting.Observables
using AbstractPlotting.GeometryBasics
using AbstractPlotting.GeometryBasics: Pyramid
using AbstractPlotting.PlotUtils
using AbstractPlotting.FileIO
using AbstractPlotting.MakieLayout
using AbstractPlotting.IntervalSets

include("reference_tests.jl")
include("unit_tests/runtests.jl")


text:
using Images, TestImages, ColorTypes, FixedPointNumbers
using CuArrays, CUDAnative
using CuTextures

# Get the input image. Use RGBA to have 4 channels since CUDA textures can have only 1, 2 or 4 channels.
img = RGBA{N0f8}.(testimage("lighthouse"))

# Create a texture memory object (CUDA array) and initilaize it with the input image content (from host).
texturearray = CuTextureArray(img)

# Create a texture object and bind it to the texture memory created above
texture = CuTexture(texturearray)

# Define an image warping kernel
function warp(dst, texture)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    u = (Float32(i) - 1f0) / (Float32(size(dst, 1)) - 1f0)
    v = (Float32(j) - 1f0) / (Float32(size(dst, 2)) - 1f0)
    x = u + 0.02f0 * CUDAnative.sin(30v)
    y = v + 0.03f0 * CUDAnative.sin(20u)
    @inbounds dst[i,j] = texture(x,y)
    return nothing
end

# Create a 500x1000 CuArray for the output (warped) image
outimg_d = CuArray{eltype(img)}(undef, 500, 1000)

# Execute the kernel
@cuda threads = (size(outimg_d, 1), 1) blocks = (1, size(outimg_d, 2)) warp(outimg_d, texture)

# Get the output image into host memory and save it to a file
outimg = Array(outimg_d)
save("imgwarp.png", outimg)


text:
using StructuresKit

#Test shear and moment calculations for a three span continuous beam.

#Ix Iy Ixy J Cw
sectionProperties = [(100.0, 100.0, 0.0, 10.0, 10.0)]

#E  ν
materialProperties = [(29500, 0.30)]

#ax ay
loadLocation = [(0.0, 0.0)]

#kx kϕ
springStiffness = [(0.0, 0.0)]

#member information
#L dL SectionProperties MaterialProperties LoadLocation SpringStiffness
memberDefinitions = [(30*12,3.0,1,1,1,1)]

#end boundary conditions
#type=1 u''=v''=ϕ''=0 (simply supported), type=2 u'=v'=ϕ'=0  (fixed), type=3 u''=v''=ϕ''=u'''=v'''=ϕ'''=0 (free end, e.g., a cantilever)
endBoundaryConditions = [1 1]

#supports
#location where u=v=ϕ=0
supports = [0.0 30.0*12/3 30.0*12*(2/3) 30.0*12]

#load  qx   qy
uniformLoad = (0.0, 5.0/1000)

z, u, v, ϕ, beamProperties = Beam.solve(memberDefinitions, sectionProperties, materialProperties, loadLocation, springStiffness, endBoundaryConditions, supports, uniformLoad)

Mxx = InternalForces.moment(z, beamProperties.dm, -v, beamProperties.E, beamProperties.Ix)

Vyy = InternalForces.shear(z, beamProperties.dm, -v, beamProperties.E, beamProperties.Ix)

#check interior moment over a support

intIndex = findall(x-> x==30.0*12/3, z)

intMxx = Mxx[intIndex[1]]
intVyy = Vyy[intIndex[1]-1]

#http://faculty-legacy.arch.tamu.edu/anichols/index_files/courses/arch331/NS8-2beamdiagrams.pdf
#Figure 36
intMxx_AISC = -0.100*uniformLoad[2]*z[intIndex[1]].^2
intVyy_AISC = -0.600*uniformLoad[2]*z[intIndex[1]] + uniformLoad[2]*3.0  #shear is checked 3 in. away from support

MxxError= abs((intMxx - intMxx_AISC)/ intMxx_AISC)
VyyError= abs((intVyy - intVyy_AISC)/ intVyy_AISC)

#accept 1% error from numerical solution
@test MxxError <= 0.01
@test VyyError <= 0.01


text:
#=
Implement an LRU (Least Recently Used) cache. It should be able to be initialized with a cache size n, and contain the following methods:

    set(key, value): sets key to value. If there are already n items in the cache and we are adding a new item, then it should also remove the least recently used item.
    get(key): gets the value at key. If no such key exists, return null.

Each operation should run in O(1) time.
=#

using Test
include("Solutions/problem52_LRUcache.jl")

LRU = LRU_cache(3)
set("a", 1)
set("b", 2)
set("c", 3)
set("d", 4)

@test get("b") == 2
@test get("c") == 3
@test get("d") == 4
@test isnothing(get("a"))

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct VirtualInertia <: ActivePowerControl
        Ta::Float64
        kd::Float64
        kω::Float64
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters of a Virtual Inertia with SRF using VSM for active power controller

# Arguments
- `Ta::Float64`: VSM inertia constant, validation range: `(0, nothing)`
- `kd::Float64`: VSM damping constant, validation range: `(0, nothing)`
- `kω::Float64`: frequency droop gain, validation range: `(0, nothing)`
- `P_ref::Float64`: Reference Power Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the VirtualInertia model are:
	ω_oc: Speed of the rotating reference frame of the virtual synchronous generator model,
	θ_oc: Phase angle displacement of the virtual synchronous generator model
- `n_states::Int`: VirtualInertia has two states
"""
mutable struct VirtualInertia <: ActivePowerControl
    "VSM inertia constant"
    Ta::Float64
    "VSM damping constant"
    kd::Float64
    "frequency droop gain"
    kω::Float64
    "Reference Power Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    "The states of the VirtualInertia model are:
	ω_oc: Speed of the rotating reference frame of the virtual synchronous generator model,
	θ_oc: Phase angle displacement of the virtual synchronous generator model"
    states::Vector{Symbol}
    "VirtualInertia has two states"
    n_states::Int
end

function VirtualInertia(Ta, kd, kω, P_ref=1.0, ext=Dict{String, Any}(), )
    VirtualInertia(Ta, kd, kω, P_ref, ext, [:ω_oc, :θ_oc], 2, )
end

function VirtualInertia(; Ta, kd, kω, P_ref=1.0, ext=Dict{String, Any}(), states=[:ω_oc, :θ_oc], n_states=2, )
    VirtualInertia(Ta, kd, kω, P_ref, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function VirtualInertia(::Nothing)
    VirtualInertia(;
        Ta=0,
        kd=0,
        kω=0,
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`VirtualInertia`](@ref) `Ta`."""
get_Ta(value::VirtualInertia) = value.Ta
"""Get [`VirtualInertia`](@ref) `kd`."""
get_kd(value::VirtualInertia) = value.kd
"""Get [`VirtualInertia`](@ref) `kω`."""
get_kω(value::VirtualInertia) = value.kω
"""Get [`VirtualInertia`](@ref) `P_ref`."""
get_P_ref(value::VirtualInertia) = value.P_ref
"""Get [`VirtualInertia`](@ref) `ext`."""
get_ext(value::VirtualInertia) = value.ext
"""Get [`VirtualInertia`](@ref) `states`."""
get_states(value::VirtualInertia) = value.states
"""Get [`VirtualInertia`](@ref) `n_states`."""
get_n_states(value::VirtualInertia) = value.n_states

"""Set [`VirtualInertia`](@ref) `Ta`."""
set_Ta!(value::VirtualInertia, val) = value.Ta = val
"""Set [`VirtualInertia`](@ref) `kd`."""
set_kd!(value::VirtualInertia, val) = value.kd = val
"""Set [`VirtualInertia`](@ref) `kω`."""
set_kω!(value::VirtualInertia, val) = value.kω = val
"""Set [`VirtualInertia`](@ref) `P_ref`."""
set_P_ref!(value::VirtualInertia, val) = value.P_ref = val
"""Set [`VirtualInertia`](@ref) `ext`."""
set_ext!(value::VirtualInertia, val) = value.ext = val


text:
#do some pre-filtering to reduce the processing time
function prefiltercrsp!(crsp::DataFrame, validpermnos::Vector{Int})::Nothing
  #rename columns to lower case

  crsp.tokeep = trues(size(crsp,1))

  crspgrp::GroupedDataFrame = groupby(crsp, :permno)
  Threads.@threads for i ∈ 1:length(crspgrp)
    scrsp::SubDataFrame = crspgrp[i]

    if scrsp.permno[1] ∉ validpermnos
      scrsp[:,:tokeep] .= false
    end
  end

  filter!(r::DataFrameRow->r.tokeep, crsp)
  select!(crsp, Not([:tokeep]))

  return nothing

end

@inline normalizecrspnames!(crsp::DataFrame) = rename!(crsp, (s::String->Symbol(lowercase(s))).(names(crsp)))

#helper to handle default case
#prepcrsp(param::AbstractDict; args...) = prepcrsp(x->x, param; args...)::DataFrame

function prepcrsp(op::Function = (noop(x)=x);
  crsptype=error("crsptype is required"),
  refreshcrsp=true,
  incsvstream::Function = IN_CSV_STREAM,
  inbinstream::Function = IN_BIN_STREAM,
  outbinstream::Function = OUT_BIN_STREAM,
  csvextension::String = CSV_EXTENSION,
  binextension::String = BIN_EXTENSION,
  crspname::String = string(crsptype),
  crsppath::String = PARAM[:crsppath],
  workingpath::String = PARAM[:workingpath],
  yearrange::Union{AbstractRange, Nothing} = nothing,
  validpermnos::Union{Nothing, Vector{Int}} = nothing,
  binprefix::String = "",
  crspnameout::String = "$binprefix$crspname",
  crspcolumns::Union{Nothing, Vector{Symbol}} = nothing,
  crsprequiredcolumns::Union{Nothing, Vector{Symbol}} = nothing,
  #usestdsharesexchanges::Bool = PARAM[:crspusestdsharesexchanges],
  makebin::Bool = true)::DataFrame

  local binpath::String = "$workingpath\\$crspnameout.$binextension"
  local csvpath::String = "$crsppath\\$crspname.$csvextension"
  local crsp::DataFrame
  local crspdateformat::DateFormat = PARAM[:wrdsdateformat]

  #makes a serialization object so we don't have to keep parsing the CSV
  if refreshcrsp
    #do the below so we can pre-filter
    crsp = incsvstream(csvpath) |> CSV.File |> f->DataFrame(f,copycols=false)
    normalizecrspnames!(crsp)
    (!isnothing(validpermnos)) && (prefiltercrsp!(crsp, validpermnos))

    #filter by year first
    crsp.date = (s->parsecrsp(Date, s, crspdateformat)).(crsp.date)
    (!isnothing(yearrange)) && filter!(r::DataFrameRow->year(r.date) ∈ yearrange, crsp) #restrict to desired years

    crsp = preprocesscrsp(crsp, crsptype)

    if !isnothing(crsprequiredcolumns)
      crsp = crsp[completecases(crsp, crsprequiredcolumns),:]
    end

    for Fposcol ∈ [:price, :shares, :mc]
      if !(all(crsp[!, Fposcol] .≥ 0.0))
        display(crsp[Fposcol .≤ 0.0,[:permno, :date, :price, :shares, :mc]])
        throw("Found negative values for $Fposcol in crsp")
      end
    end

    crsp = op(crsp)
    #throw("crsp colummns: $(propertynames(crsp))")

    if !isnothing(crspcolumns)
      select!(crsp, crspcolumns)
    end

    makebin && outbinstream(binpath, crsp)

  else
    (!makebin) && (@warn "makebin = false and loading a bin! If this works, beware of stale data....")
    crsp = inbinstream(binpath)
  end

  @info "CRSP data loaded and/or saved into file $crspnameout.$BIN_EXTENSION"

  return crsp
end

#function
#creates crsp returns
function computecrspreturns(crsp::DataFrame)
  crspparsedmissing::Vector{Float64} = PARAM[:crspparsedmissing]
  retlisted::Vector{MFloat64} = (s::MString->parsecrsp(Float64, s, crspparsedmissing)).(crsp.ret)
  dlret::Vector{MFloat64} = (s::MString->parsecrsp(Float64, s, crspparsedmissing)).(crsp.dlret)

  retlisted = (r->
    !ismissing(r) && (r > -1.0) ? r : missing).(retlisted)
  dlret .= (r->
    !ismissing(r) && (r > -1.0) ? r : missing).(dlret)

  ret::Vector{MFloat64} = ((retlisted,dlret)->
    ismissing(retlisted) && ismissing(dlret) ? missing :
      ismissing(dlret) ? retlisted :
      ismissing(retlisted) ? dlret :
      (1+retlisted)*(1+dlret) - 1.
    ).(retlisted,dlret)

  return ret
end

#computes a cumulative index for each stock
function computecrsptotalreturnindex(crsp::AbstractDataFrame; Fret::Symbol=:ret)
  @assert (crsp[!, Fret] .≥ -1.0) |> all
  @assert (crsp[!, Fret] .!== missing) |> all
  @assert issorted(crsp, [:permno, :date])

  rets = crsp[!, [:permno, :date, Fret]] |> deepcopy
  rets.retP1 = rets[!, Fret] .+ 1

  transform!(groupby(rets, :permno), :retP1 => cumprod => :tridx)

  return rets.tridx
end

function computecrspprices(crsp::DataFrame, crsptype::Symbol)
  #Now work on the price
  #in this case, use the mean of the two values
  #price cannot be zero
  crsp.prc .= (
    prc->coalesce(iszero(prc),true) ? missing : prc).(crsp.prc)
  crsp.dlprc .= (
    dlprc->coalesce(iszero(dlprc),true) ? missing : dlprc).(crsp.dlprc)

  @inline ismissingorzero(p::Float64) = iszero(p) ? true : false
  @inline ismissingorzero(::Missing) = true

  price = Vector{MFloat64}(undef, size(crsp,1))
  Threads.@threads for i ∈ 1:length(price)
    if ismissingorzero(crsp.dlprc[i])
      price[i] = crsp.prc[i]
    elseif ismissingorzero(crsp.prc[i])
      price[i] = crsp.dlprc[i]
    else
      price[i] = (crsp.prc[i]+crsp.dlprc[i])/2
    end
  end

  #add in the alternative price if nothing else is available
  #(data available monthly only)
  if crsptype == :monthly
  Threads.@threads for i ∈ 1:length(price)
      if ismissing(price[i]) && (!ismissingorzero(crsp.altprc[i]))
        price[i] = crsp.altprc[i]
      end
    end
  end

  #negative values computed via bid-ask spread method
  price .= (abs).(price)

  return price
end

function filtercrsp!(crsp::DataFrame;)

    #NOTE: SHRCD==10,11 : US Common Shares
    #NOTE:  :EXCHCD ∈ [1,2,3] : NYSE, Amex, NASDAQ
    delete!(crsp, (!).(completecases(crsp[:,[:shrcd, :exchcd]])))
    local shrcdvalid::Set{Int} = Set{Int}((10,11))
    local exchcdvalid::Set{Int} = Set{Int}((1,2,3))
    filter!(r::DataFrameRow->
      (r.shrcd ∈ shrcdvalid) && (r[:exchcd] ∈ exchcdvalid), crsp)

    return crsp
end

computecrspshares(crsp) = crsp.shrout .* 1000.
function computecrspmarketcap(crsp::DataFrame)

  #adjust the units and compute the market cap
  mc::Vector{MInt} = Vector{MInt}((round).(crsp.price .* crsp.shares))
  mc .= (i::MInt->coalesce(iszero(i),true) ? missing : i).(mc)

  return mc
end

function computerangemeasures!(crsp::DataFrame; disttol::Float64=10^-4,
    periods2stale::Day = Day(30))

  issorted(crsp, [:permno, :date]) || error("crsp must be sorted on permno, date")
  posormissing(::Missing) = missing
  posormissing(x::Real) = x>0.0 ? x : missing

  #weird crsp thing that bidlo is a bid under certain circumstances, represented by a neg val
  crsp.lo = (posormissing).(crsp.bidlo)
  crsp.hi = (posormissing).(crsp.askhi)

  lagwithin2!(crsp, [:price], :permno, date=:date, maxnotstale = periods2stale,
    laggedvals = [:Ldistadjclose])

  #adjust the lagged close prices for distributions
  #note a distribuition here is defined as anything that causes the price return
  #to deveiate from the crsp ret, so it can be negative (consider splits)
  priceret::Vector{MFloat64} = crsp.price ./ crsp.Ldistadjclose .- 1.
  dist::Vector{MFloat64} = (crsp.ret .- priceret) .* crsp.Ldistadjclose
  dist[(d->((!ismissing(d)) && (abs(d)<disttol))).(dist)] .= 0.0 #assume such small changes are rounding error
  crsp.Ldistadjclose .-= dist
  crsp.Ldistadjclose .= (posormissing).(crsp.Ldistadjclose)

end

function computedollarvol(crsp::DataFrame)
  #use this to compute dollar volume- take the price midpoint if available
  #println(names(crsp))

  dprice::Vector{MFloat64} = ((open, close)->
    ifelse(open !== missing, (open + close)/2, close)).(crsp.open, crsp.price)

  parsevol(::Missing) = missing
  parsevol(f::Real) = ifelse(f > 0.0, f, missing)

  @assert  (v->(v===missing) || (v≥0)).(crsp.vol) |> all
  vol::Vector{MFloat64} = (parsevol).(crsp.vol)

  return dprice .* crsp.vol
end

#normalizes the dollar volume by the point in time total dollar volume
function computenvol!(crsp::DataFrame, Fdvol::Symbol)

  (:nvol ∈ names(crsp)) && (sum((!ismissing).(crsp.nvol)) > 0) && error(
    ":nvol already contains data")

  crsp.nvol =  Vector{MFloat64}(undef, size(crsp,1))

  #use this to compute dollar volume- take the price midpoint if available
  scrsps::GroupedDataFrame = groupby(crsp, :date)
  for i ∈ 1:length(scrsps)
    scrsp::SubDataFrame = scrsps[i]
    scrsp.nvol .= scrsp[!, Fdvol]./sum(skipmissing(scrsp[!, Fdvol]))
  end

  return nothing
end

#should probably make this based on time interval
function trailingsigma!(crspd::DataFrame; Fendpoint::Symbol = "Fendpoint is required")

  (:sigma ∉ names(crspd)) || (error(":sigma already in crspd"))
  crspd.sigma = Vector{MFloat64}(undef, size(crspd,1))

  trailingmeasure!(std, crspd, PARAM[:sigmacalendarinterval],PARAM[:sigmaminpoints],
    Ftarget=:ret,
    Fmeasure=PARAM[:sigmaname],
    Fconditional=Fendpoint)::Nothing
end

function preprocesscrsp(crsp::DataFrame, crsptype::Symbol,
  crspdateformat::DateFormat = PARAM[:wrdsdateformat])

  local scrsps::GroupedDataFrame

  #NOTE: SHRCD==10,11 : US Common Shares
  #NOTE:  :EXCHCD ∈ [1,2,3] : NYSE, Amex, NASDAQ
  local shrcdvalid::Set{Int} = Set{Int}((10,11))
  local exchcdvalid::Set{Int} = Set{Int}((1,2,3))

  #=crsp.tokeep = trues(size(crsp,1))
  Threads.@threads for r ∈ eachrow(crsp)
    r.tokeep = ((r.shrcd ∈ shrcdvalid) &&
      (r.exchcd ∈ exchcdvalid) &&
      (!ismissing(r.shrcd)) &&
      (!ismissing(r.exchcd)))
  end
  crsp = crsp[crsp.tokeep,:]=#
  delete!(crsp, (!).(completecases(crsp[:,[:shrcd, :exchcd]])))
  filter!(r::DataFrameRow->(r.shrcd ∈ shrcdvalid) && (r.exchcd ∈ exchcdvalid), crsp)

  #Doing this before the dedup eliminates rows without returns,
  #prices and shares outstanding first. Implicitly requires that
  #these fields exist
  crsp.ret = computecrspreturns(crsp)
  crsp.price = computecrspprices(crsp, crsptype)

  #require valid returns and prices
  crsp = crsp[completecases(crsp, [:ret,:price]),:]

  #following FF UMD, force price and return availability before computing TR
  crsp.tridx = computecrsptotalreturnindex(crsp)

  crsp.lret = (log).(1.0 .+ crsp.ret) #make log returns
  select!(crsp, Not([:prc, :dlprc, :dlret])) #cleanup

  #now dedup crsp
  dedupctr::Int = 0
  crsp.keep = trues(size(crsp,1))

  scrsps = groupby(crsp, [:permno, :date])
  Threads.@threads for i ∈ 1:length(scrsps)
    local scrsp::SubDataFrame = scrsps[i]
    local partialscrsp::SubDataFrame

    (size(scrsp,1) == 1) && continue
    scrsp.keep .= scrsp.distcd .== 1232 #sometimes dups on the distribution type- prefer to keep ord dividends
    (sum(scrsp.keep) == 1) && continue

    #if we have too many 1232s, we only will select from these
    partscrsp = (sum(scrsp.keep) > 1) ? view(scrsp, scrsp.keep, :) : scrsp

    #try to take the minimum distribution code
    nonmissingdistcds::Vector{Int} = collect(skipmissing(partscrsp.distcd))
    if length(nonmissingdistcds) ≠ 0
      partscrsp.keep .= (partscrsp.distcd .== minimum(nonmissingdistcds))
    end
    (sum(partscrsp.keep) == 1) && continue

    #again, see if we have too many  of the minimum distcd, or the distcd was missing
    partscrsp = (sum(partscrsp.keep) > 1) ? view(partscrsp, partscrsp.keep, :) : partscrsp

    #finally, check if the rows are identical- if so, keep 1, if not, drop the rows
    partscrsp.keep .= false
    if eachrowequal(partscrsp)
      partscrsp.keep[1] = true
    end

  end

  delete!(crsp, (!).(crsp.keep))

  #compute market cap here since it is used in selection criteria
  crsp.open = (abs).(crsp.openprc)
  crsp.shares = computecrspshares(crsp)
  crsp.mc = computecrspmarketcap(crsp)

  #=print("time assertion: ")
  @time =#@assert (!any(nonunique(crsp[!,[:permno, :date]])))

  return crsp
end

#parse a value where the type is known: Generic case
function parsecrsp(::Type{T}, s::String, crspparsedmissing::AbstractVector{Float64}) where T

  #if it doesn't parse to the right type, set to missing
  v::Union{T,Missing} = something(tryparse(T, s), missing)

  #check against the list of missing codes
  (!ismissing(v)) && (v ∈ crspparsedmissing) && (v=missing)

  return v
end

#the date case
function parsecrsp(::Type{Date}, s::String, crspdateformat::DateFormat)
  #try
    Dates.Date(s, crspdateformat)
  #catch err
  #  println("could not parse $s with date format $crspdateformat")
  #  error(err)
  #end
end

#helper methods and special cases
parsecrsp(::Type{<:Any}, ::Missing) = missing
parsecrsp(::Type{<:Any}, ::Missing, ::Any) = missing
parsecrsp(::Type{Date}, i::Int, crspdateformat::DateFormat) = parsecrsp(Date, "$i", crspdateformat)
parsecrsp(::Type{Int}, s::String) = something(tryparse(Int, s), missing)
parsecrsp(::Type{Int}, i::Int) = i


text:
# everything that should move to learnbase

# Eltype, Labelcount, Arraydimensions
abstract type LabelEncoding{T,K,M} end
# act as scalar in broadcast, see julia #18618
broadcastable(x::LabelEncoding) = Ref(x)

const BinaryLabelEncoding{T,M} = LabelEncoding{T,2,M}
const VectorLabelEncoding{T,K} = LabelEncoding{T,K,1}
const MatrixLabelEncoding{T,K} = LabelEncoding{T,K,2}

"""
    nlabel(obj) -> Int

Returns the number of labels represented in the given object `obj`.

    julia> nlabel([:yes,:no,:yes,:yes])
    2
"""
nlabel(::Type{<:BinaryLabelEncoding}) = 2
nlabel(::Type{LabelEncoding{T,K,M}}) where {T,K,M} = Int(K)
nlabel(::Type{Any}) = throw(ArgumentError("number of labels could not be inferred for the given type"))
nlabel(::Type{T}) where {T} = nlabel(supertype(T))
nlabel(::LabelEncoding{T,K}) where {T,K} = Int(K)
nlabel(itr) = length(label(itr))

"""
    label(obj) -> Vector

Returns the labels represented in the given object `obj`.
Note that the order of the labels matters.
In the case of two labels, the first element represents the positive
label and the second element the negative label.

    julia> label([:yes,:no,:yes,:yes])
    2-element Array{Symbol,1}:
     :yes
     :no

    julia> label(LabelEnc.ZeroOne())
    2-element Array{Float64,1}:
     1.0
     0.0
"""
label(itr) = _arrange_label(unique(itr))
label(A::AbstractVector) = _arrange_label(unique(A))
label(A::AbstractArray{T,N}) where {T,N} = throw(MethodError(label, (A,)))
label(A::AbstractMatrix{<:Union{Number,Bool}}; obsdim = LearnBase.default_obsdim(A)) = label(A, convert(LearnBase.ObsDimension,obsdim))
label(A::AbstractMatrix{<:Union{Number,Bool}}, ::Union{ObsDim.Constant{2},ObsDim.Last}) = collect(1:size(A,1))
label(A::AbstractMatrix{<:Union{Number,Bool}}, ::ObsDim.First) = collect(1:size(A,2))

# make sure pos label is first
_arrange_label(lbl::Vector) = lbl
_arrange_label(lbl::Vector{<:Bool}) = [true,false]
function _arrange_label(lbl::Vector{T}) where {T<:Number}
    if length(lbl) == 2
        if minimum(lbl) == 0 && maximum(lbl) == 1
            lbl[1] = T(1)
            lbl[2] = T(0)
        elseif minimum(lbl) == -1 && maximum(lbl) == 1
            lbl[1] = T(1)
            lbl[2] = T(-1)
        elseif minimum(lbl) == 1 && maximum(lbl) == 2
            lbl[1] = T(1)
            lbl[2] = T(2)
        end
    end
    lbl
end

"""
    labeltype(::Type{<:LabelEncoding}) -> DataType

Determine the type of the labels represented by the label encoding `T`
"""
labeltype(::Type{MatrixLabelEncoding{T}}) where {T} = T
labeltype(::Type{VectorLabelEncoding{T}}) where {T} = T
labeltype(::Type{LabelEncoding{T,K,M}}) where {T,K,M} = T
labeltype(::Type{Any}) = Any
labeltype(::Type{T}) where {T} = labeltype(supertype(T))
labeltype(lm::LabelEncoding{T}) where {T} = T

"""
    ind2label(index, encoding)

Converts the given `index` into the corresponding label defined
by the `encoding`. Note that in the binary case, `index == 1`
represents the positive label and `index == 2` the negative label.

    julia> ind2label(1, LabelEnc.ZeroOne(Float32))
    1.0f0

    julia> ind2label(2, LabelEnc.ZeroOne(Float32))
    0.0f0
"""
function ind2label end

"""
    label2ind(label, encoding) -> Int

Converts the given `label` into the corresponding index defined
by the encoding. Note that in the binary case, the positive label
will result in the index `1` and the negative label in the index `2`
respectively

    julia> label2ind(:no, LabelEnc.NativeLabels([:yes,:no]))
    2

    julia> label2ind(1, LabelEnc.MarginBased())
    1

    julia> label2ind(-1, LabelEnc.MarginBased())
    2
"""
function label2ind end

"""
    poslabel(encoding)

If the encoding is binary it will return the positive label of it.
The function will throw an error otherwise.

    julia> poslabel(LabelEnc.MarginBased(Float32))
    1.0f0
"""
function poslabel(values::AbstractArray)
    lbl = label(values)
    length(lbl) == 2 || throw(ArgumentError("The given object has more or less than two labels, thus poslabel is not defined."))
    lbl[1]
end

"""
    neglabel(encoding)

If the encoding is binary it will return the negative label of it.
The function will throw an error otherwise.

    julia> neglabel(LabelEnc.MarginBased(Float32))
    -1.0f0
"""
function neglabel(values::AbstractArray)
    lbl = label(values)
    length(lbl) == 2 || throw(ArgumentError("The given object has more or less than two labels, thus neglabel is not defined."))
    lbl[2]
end

"""
    labelenc(obj) -> LabelEncoding

Tries to determine the most approriate label-encoding to describe the
given object `obj` based on the result of `label(obj)`. Note that in
most cases this function is not typestable.

    julia> labelenc([:yes,:no,:no,:yes,:maybe])
    MLLabelUtils.LabelEnc.NativeLabels{Symbol,3}(Symbol[:yes,:no,:maybe],Dict(:yes=>1,:maybe=>3,:no=>2))

    julia> labelenc([1,0,0,1,0,1])
    MLLabelUtils.LabelEnc.ZeroOne{Int64,Float64}(0.5)

    julia> labelenc(Int8[1,-1,-1,1,-1,1])
    MLLabelUtils.LabelEnc.MarginBased{Int8}()
"""
function labelenc end

"""
    islabelenc(obj, encoding) -> Bool

Checks is the given object `obj` can be described as being produced
by the given `encoding` in which case the function returns true,
or false otherwise.

    julia> islabelenc([1,0,1], LabelEnc.ZeroOne)
    true

    julia> islabelenc([1,-1,1], LabelEnc.ZeroOne)
    false
"""
function islabelenc end

"""
    isposlabel(x, encoding) -> Bool

Checks if the given value `x` can be interpreted as the positive label
given the `encoding`. This function takes potential classification
rules into account.

    julia> isposlabel(0.6, LabelEnc.ZeroOne(0.5))
    true

    julia> isposlabel(0.4, LabelEnc.ZeroOne(0.5))
    false

    julia> isposlabel(:b, LabelEnc.NativeLabels([:a,:b]))
    false
"""
function isposlabel end

"""
    isneglabel(x, encoding) -> Bool

Checks if the given value `x` can be interpreted as the negative label
given the `encoding`. This function takes potential classification
rules into account.

    julia> isneglabel(0.6, LabelEnc.ZeroOne(0.5))
    false

    julia> isneglabel(0.4, LabelEnc.ZeroOne(0.5))
    true

    julia> isneglabel(:b, LabelEnc.NativeLabels([:a,:b]))
    true
"""
function isneglabel end

"""
    classify(x, encoding)

Returns the classified version of `x` given the `encoding`.
Which means that if `x` can be interpreted as a positive label,
the positive label of `encoding` is returned; the negative otherwise.

    julia> classify(0.6, LabelEnc.ZeroOne(UInt8))
    0x01

    julia> classify(0.4, LabelEnc.ZeroOne(UInt8))
    0x00

    julia> classify([0.1,0.2,0.6,0.1], LabelEnc.OneOfK)
    3
"""
function classify end

"""
    classify!(out, x, encoding)

Same as `classify`, but uses `out` to store the result.

    julia> buffer = zeros(2);
    julia> classify!(buffer, [0.4,0.6], LabelEnc.ZeroOne)
    2-element Array{Float64,1}:
     0.0
     1.0
"""
function classify! end

"""
    convertlabel(new_encoding, x, [old_encoding])

Converts the given value/array `x` from the `old_encoding` into the
`new_encoding`. Note that if `old_encoding` is not specified it will
be derived automaticaly using `labelenc`.

    julia> convertlabel(LabelEnc.MarginBased, [0, 1, 1, 0, 0])
    5-element Array{Int64,1}:
     -1
      1
      1
     -1
     -1

    julia> convertlabel([:yes,:no], [0, 1, 1, 0, 0])
    5-element Array{Symbol,1}:
     :no
     :yes
     :yes
     :no
     :no

For more information on the available encodings, see `?LabelEnc`.

    convertlabel(new_encoding, x, [old_encoding], [obsdim])

When working with `OneOfK` one can additionally specifify which
dimension of the array denotes the observations using `obsdim`

    julia> convertlabel(LabelEnc.OneOfK, [0, 1, 1, 0, 0], obsdim = 2)
    2×5 Array{Int64,2}:
     0  1  1  0  0
     1  0  0  1  1
"""
function convertlabel end
function convertlabel! end

"""
    convertlabel(new_encoding, vec::AbstractVector, [old_encoding]) -> (Readonly)MappedArray

Creates a lazy view into `vec` that makes it look like it is
in the encoding specified by `new_encoding`, while it is actually
preserved as being of `old_encoding`.

This method only works for label-encodings that are vector-based
(i.e. pretty much all but `OneOfK`). The resulting MappedArray
will be writeable unless `old_encoding` is of type `OneVsRest`,
in which case the result will be a `ReadonlyMappedArray`.
"""
function convertlabelview end

"""
    labelmap(obj) -> Dict

Computes a mapping from the labels in `obj` to all the individual
element-indices in `obj` that correspond to that label

    julia> labelmap([0, 1, 1, 0, 0])
    Dict{Int64,Array{Int64,1}} with 2 entries:
      0 => [1,4,5]
      1 => [2,3]
"""
function labelmap end

"""
    labelmap!(dict, idx, elem) -> Dict

Updates the given label-map `dict` with the new element `elem`,
which is assumed to be associated with the index `idx`.

    julia> lm = labelmap([0, 1, 1, 0, 0])
    Dict{Int64,Array{Int64,1}} with 2 entries:
      0 => [1,4,5]
      1 => [2,3]

    julia> labelmap!(lm, 6, 0)
    Dict{Int64,Array{Int64,1}} with 2 entries:
      0 => [1,4,5,6]
      1 => [2,3]

    julia> labelmap!(lm, 7:8, [1,0])
    Dict{Int64,Array{Int64,1}} with 2 entries:
      0 => [1,4,5,6,8]
      1 => [2,3,7]
"""
function labelmap! end

"""
    labelfreq(obj) -> Dict

Computes the absolute frequencies for each label in `obj`.

    julia> labelfreq([0, 1, 1, 0, 0])
    Dict{Int64,Int64} with 2 entries:
      0 => 3
      1 => 2
"""
function labelfreq end

"""
    labelfreq!(dict, obj) -> Dict

updates the given label-frequency-map `dict` with the absolute
frequencies for each label in `obj`

    julia> ld = labelfreq([0, 1, 1, 0, 0])
    Dict{Int64,Int64} with 2 entries:
      0 => 3
      1 => 2

    julia> labelfreq!(ld, [1,0,0])
    Dict{Int64,Int64} with 2 entries:
      0 => 5
      1 => 3
"""
function labelfreq! end

"""
    labelmap2vec(dict) -> Vector

Inverse function of labelmap.
Computes an `array` of labels by element-wise
traversal of the entries in `dict`.

    julia> labelvec = [:yes,:no,:no,:yes,:yes]

    julia> lm = labelmap(labelvec)
    Dict{Symbol,Array{Int64,1}} with 2 entries:
        :yes => [1, 4, 5]
        :no  => [2, 3]

    julia> labelmap2vec(lm)
    5-element Array{Symbol,1}:
        :yes
        :no
        :no
        :yes
        :yes
"""
function labelmap end


text:
include("molzahn.jl")
include("sliwak.jl")
include("treshold_functions.jl")

const heuristic_functions = Dict(
    "molzahn" => molzahn_heuristic,
    "sliwak" => sliwak_heuristic
)

const update_functions = Dict(
    "molzahn" => update_molzahn,
    "sliwak" => update_sliwak
)

function update_cliquetree!(cliquetree, nb_clique, edge)
    sorted_edge = sort(edge)
    for i in 1:length(cliquetree)
        for j in 1:2
            if cliquetree[i][j] in sorted_edge
                cliquetree[i][j] = nb_clique
            elseif cliquetree[i][j] >= sorted_edge[2]
                cliquetree[i][j] -= 2
            elseif cliquetree[i][j] >= sorted_edge[1]
                cliquetree[i][j] -= 1
            end
        end
    end
end

function merge_clique(clique, cliquetree, edge)
    c1 = clique[edge[1]]
    c2 = clique[edge[2]]
    new_clique = filter(c -> !(c in [c1, c2]), deepcopy(clique))
    push!(new_clique, union(c1, c2))
    new_cliquetree = filter(e -> e != edge, deepcopy(cliquetree))
    update_cliquetree!(new_cliquetree, length(new_clique), edge)
    return new_clique, new_cliquetree
end

function minimize(clique, cliquetree; heuristic_name, merge_kwargs)
    heuristic = heuristic_functions[heuristic_name]
    deltas = [heuristic(clique, edge; merge_kwargs=merge_kwargs) for edge in cliquetree]
    index = argmin(deltas)
    return cliquetree[index]
end

function merge_dec(clique::AbstractVector, cliquetree::AbstractVector,
                   heuristics=["sliwak", "molzahn"], heuristic_switch=[1];
                   treshold_name::AbstractString="clique_nv_up",
                   merge_kwargs::AbstractDict{String, <:Any})
    merge_kwargs = Dict{String, Any}(merge_kwargs)
    merge_kwargs["clique"] = clique
    merge_kwargs["cliquetree"] = cliquetree
    treshold = treshold_functions[treshold_name](merge_kwargs)
    iter = 1
    heuristic_index = 1
    heuristic_switch_index = 1
    while treshold(merge_kwargs)
        if iter > heuristic_switch[heuristic_switch_index]
            heuristic_switch_index = heuristic_switch_index == length(heuristic_switch) ? 1 : heuristic_switch_index + 1
            heuristic_index = heuristic_index == length(heuristics) ? 1 : heuristic_index + 1
        end
        heuristic_name = heuristics[heuristic_index]
        edge = minimize(merge_kwargs["clique"], merge_kwargs["cliquetree"]; heuristic_name=heuristic_name, merge_kwargs=merge_kwargs)
        update_functions[heuristic_name](clique=merge_kwargs["clique"], cliquetree=merge_kwargs["cliquetree"], edge=edge, merge_kwargs=merge_kwargs)
        merge_kwargs["clique"], merge_kwargs["cliquetree"] = merge_clique(merge_kwargs["clique"], merge_kwargs["cliquetree"], edge)
        iter += 1
    end
    clique, cliquetree = merge_kwargs["clique"], merge_kwargs["cliquetree"]
    return clique, map(c -> map(Int, c), cliquetree)
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using PlotlyJS

function house()
    trace1 = scatter()
    x0 = [2,   2,   5.5, 9,   9, 2,   5, 5, 6]
    y0 = [1,   5.5, 9.5, 5.5, 1, 5.5, 1, 4, 4]
    x1 = [2,   5.5, 9,   9,   2, 9,   5, 6, 6]
    y1 = [5.5, 9.5, 5.5, 1,   1, 5.5, 4, 4, 1]
    shapes = line(x0, x1, y0, y1; xref="x", yref="y")
    plot([trace1],
         Layout(;shapes=shapes, xaxis_range=(1, 10), yaxis_range=(0, 10)))
end

function house2()
    trace1 = scatter()
    _p = string("M 2 1 L 2 5.5 L 5.5 9.6 L 9 5.5 L 9 1 L 2 1 ",
                "M 2 5.5 L 9 5.5 ",
                "M 5 1 L 5 4 L 6 4 L 6 1 Z")
    plot([trace1],
         Layout(;shapes=[path(_p)], xaxis_range=(1, 10), yaxis_range=(0, 10)))
end

function clusters()
    @eval using Distributions
    x0 = rand(Normal(2, 0.45), 300)
    y0 = rand(Normal(2, 0.45), 300)
    x1 = rand(Normal(6, 0.4), 200)
    y1 = rand(Normal(6, 0.4), 200)
    x2 = rand(Normal(4, 0.3), 200)
    y2 = rand(Normal(4, 0.3), 200)

    data = [scatter(;x=x0, y=y0, mode="markers"),
              scatter(;x=x1, y=y1, mode="markers"),
              scatter(;x=x2, y=y2, mode="markers"),
              scatter(;x=x1, y=y0, mode="markers")]

    args = [(x0, y0, "blue"), (x1, y1, "orange"), (x2, y2, "green"),
            (x1, y0, "red")]
    shapes = [circle(x0=minimum(x), y0=minimum(y),
                     x1=maximum(x), y1=maximum(y);
                     opacity=0.2, fillcolor=c, line_color=c)
              for (x, y, c) in args]
    plot(data, Layout(;height=400, width=480, showlegend=false, shapes=shapes))
end

function temperature()
    x = ["2015-02-01", "2015-02-02", "2015-02-03", "2015-02-04", "2015-02-05",
         "2015-02-06", "2015-02-07", "2015-02-08", "2015-02-09", "2015-02-10",
         "2015-02-11", "2015-02-12", "2015-02-13", "2015-02-14", "2015-02-15",
         "2015-02-16", "2015-02-17", "2015-02-18", "2015-02-19", "2015-02-20",
         "2015-02-21", "2015-02-22", "2015-02-23", "2015-02-24", "2015-02-25",
         "2015-02-26", "2015-02-27", "2015-02-28"]
    y = rand(1:20, length(x))
    data = scatter(;x=x, y=y, name="temperature", mode="line")

    shapes = rect(["2015-02-04", "2015-02-20"], ["2015-02-06", "2015-02-22"],
                  0, 1; fillcolor="#d3d3d3", opacity=0.2, line_width=0,
                  xref="x", yref="paper")
    plot(data, Layout(shapes=shapes, width=500, height=500))
end

function vlines1()
    # one scalar argument produces one line. Need to wrap in an array because
    # layout.shapes should be an array
    trace1 = scatter(;x=1:10, y=(1:10).^2)
    shapes = [vline(2)]
    plot([trace1], Layout(;shapes=shapes))
end

function vlines2()
    # one argument draws a vertical line up the entire plot
    trace1 = scatter(;x=1:10, y=(1:10).^2)
    shapes = vline([2, 6])
    plot([trace1], Layout(;shapes=shapes))
end

function vlines3()
    # yref paper makes the 2nd and 3rd arguments on a (0, 1) scale vertically
    # so 0.5 is 1/2 through the plot regardless of the values on y-axis
    trace1 = scatter(;x=1:10, y=(1:10).^2)
    shapes = vline([2, 6], 0, 0.5; yref="paper")
    plot([trace1], Layout(;shapes=shapes))
end

function vlines4()
    # Whichever argument is a scalar is repeated
    trace1 = scatter(;x=1:10, y=(1:10).^2)
    shapes = vline([2, 6], 0, [0.5, 0.75]; yref="paper")
    plot([trace1], Layout(;shapes=shapes))
end

function vlines5()
    # we can also set arbitrary line attributes line color and dash
    trace1 = scatter(;x=1:10, y=(1:10).^2)
    shapes = vline([2, 6], 0, [0.5, 0.75]; yref="paper",
                   line_color="green", line_dash="dashdot")
    plot([trace1], Layout(;shapes=shapes))
end

function hlines1()
    # one scalar argument produces one line. Need to wrap in an array because
    # layout.shapes should be an array
    trace1 = scatter(;x=1:10, y=(1:10).^2)
    shapes = [hline(2)]
    plot([trace1], Layout(;shapes=shapes))
end

function hlines2()
    # one argument draws a horizontal line across the entire plot
    trace1 = scatter(;x=1:10, y=(1:10).^2)
    shapes = hline([25, 81])
    plot([trace1], Layout(;shapes=shapes))
end

function hlines3()
    # xref paper makes the 2nd and 3rd arguments on a (0, 1) scale horizontally
    # so 0.5 is 1/2 through the plot regardless of the values on x-axis
    trace1 = scatter(;x=1:10, y=(1:10).^2)
    shapes = hline([25, 81], 0, 0.5; xref="paper")
    plot([trace1], Layout(;shapes=shapes))
end

function hlines4()
    # Whichever argument is a scalar is repeated
    trace1 = scatter(;x=1:10, y=(1:10).^2)
    shapes = hline([25, 81], 0, [0.5, 0.75]; xref="paper")
    plot([trace1], Layout(;shapes=shapes))
end

function hlines5()
    # we can also set arbitrary line attributes line color and dash
    trace1 = scatter(;x=1:10, y=(1:10).^2)
    shapes = hline([25, 81], 0, [0.5, 0.75]; xref="paper",
                   line_color="green", line_dash="dashdot")
    plot([trace1], Layout(;shapes=shapes))
end


text:
import PAPI

const N = 2_000

@show LinAlg.peakflops(N)

function peakflops(n=2000)
    a = ones(Float64, 100,100)
    t = @elapsed a2 = a*a
    a = ones(Float64, n, n)
    res = @PAPI.flops a2 = a*a
    @assert a2[1,1] == n
    return res
end

flpops, mflops, rtime, ptime = peakflops(N)
println("flpops: $flpops, mflops: $mflops, rtime (s): $rtime, ptime (s): $ptime")
println("peakflops: $(2*Float64(N)^3 / rtime)")


text:
function ADMMetaData(arg0::ADMFile)
    return ADMMetaData((ADMFile,), arg0)
end

function get_adm_file(obj::ADMMetaData)
    return jcall(obj, "getADMFile", ADMFile, ())
end

function get_center_body(obj::ADMMetaData)
    return jcall(obj, "getCenterBody", CelestialBody, ())
end

function get_center_name(obj::ADMMetaData)
    return jcall(obj, "getCenterName", JString, ())
end

function get_comment(obj::ADMMetaData)
    return jcall(obj, "getComment", List, ())
end

function get_has_creatable_body(obj::ADMMetaData)
    return jcall(obj, "getHasCreatableBody", jboolean, ())
end

function get_launch_number(obj::ADMMetaData)
    return jcall(obj, "getLaunchNumber", jint, ())
end

function get_launch_piece(obj::ADMMetaData)
    return jcall(obj, "getLaunchPiece", JString, ())
end

function get_launch_year(obj::ADMMetaData)
    return jcall(obj, "getLaunchYear", jint, ())
end

function get_object_id(obj::ADMMetaData)
    return jcall(obj, "getObjectID", JString, ())
end

function get_object_name(obj::ADMMetaData)
    return jcall(obj, "getObjectName", JString, ())
end

function get_time_scale(obj::ADMMetaData)
    return jcall(obj, "getTimeScale", TimeScale, ())
end

function get_time_system(obj::ADMMetaData)
    return jcall(obj, "getTimeSystem", CcsdsTimeScale, ())
end

function set_center_body(obj::ADMMetaData, arg0::CelestialBody)
    return jcall(obj, "setCenterBody", void, (CelestialBody,), arg0)
end

function set_center_name(obj::ADMMetaData, arg0::JString)
    return jcall(obj, "setCenterName", void, (JString,), arg0)
end

function set_comment(obj::ADMMetaData, arg0::List)
    return jcall(obj, "setComment", void, (List,), arg0)
end

function set_has_creatable_body(obj::ADMMetaData, arg0::jboolean)
    return jcall(obj, "setHasCreatableBody", void, (jboolean,), arg0)
end

function set_launch_number(obj::ADMMetaData, arg0::jint)
    return jcall(obj, "setLaunchNumber", void, (jint,), arg0)
end

function set_launch_piece(obj::ADMMetaData, arg0::JString)
    return jcall(obj, "setLaunchPiece", void, (JString,), arg0)
end

function set_launch_year(obj::ADMMetaData, arg0::jint)
    return jcall(obj, "setLaunchYear", void, (jint,), arg0)
end

function set_object_id(obj::ADMMetaData, arg0::JString)
    return jcall(obj, "setObjectID", void, (JString,), arg0)
end

function set_object_name(obj::ADMMetaData, arg0::JString)
    return jcall(obj, "setObjectName", void, (JString,), arg0)
end

function set_time_system(obj::ADMMetaData, arg0::CcsdsTimeScale)
    return jcall(obj, "setTimeSystem", void, (CcsdsTimeScale,), arg0)
end


text:
using StitchDPI, ImageBase, FileIO, Test

img = rand(3,3,3,4)
multiproc_write("test.nhdr", img, Float32)

img2 = load("test.nhdr")

@test all(img2.==Float32.(img))


text:
mutable struct LMDBDict{K,V}
    env::LMDB.Environment
    dbi::LMDB.DBI
    function LMDBDict{K,V}(env::LMDB.Environment, dbi::LMDB.DBI) where {K,V}
        x = new{K,V}(env, dbi)
        finalizer(x) do d
            LMDB.close(d.env,d.dbi)
            LMDB.close(d.env)
        end
        x
    end
end
function LMDBDict{K,V}(path::String; readonly = false) where {K,V}
    flags = readonly ? MDB_RDONLY : zero(Cuint)
    env = LMDB.create()
    open(env, path)
    #A transaction just for getting a DBI handle
    dbi = LMDB.start(env,flags=flags) do txn
        LMDB.open(txn)
    end
    LMDBDict{K,V}(env, dbi)
end
LMDBDict(path::String; kwargs...) = LMDBDict{String, Vector{Uint8}}(path; kwargs...)

function cursor_do(f, d; readonly = false)
    txnflags = readonly ? Cuint(LMDB.MDB_RDONLY) : Cuint(0)
    LMDB.start(d.env, flags = txnflags) do txn
        LMDB.open(txn,d.dbi) do cur
            f(cur)
        end
    end
end

function txn_dbi_do(f, d; readonly = false)
    txnflags = readonly ? Cuint(LMDB.MDB_RDONLY) : Cuint(0)
    LMDB.start(d.env, flags = txnflags) do txn
        f(txn, d.dbi)
    end
end

function list_dirs(d::LMDBDict{String}; prefix = "", sep = '/')
    cursor_do(d, readonly = true) do cur
        bprefix = Vector{UInt8}(prefix)
        iter = LMDB.LMDBIterator(cur,LMDB.DirectoryLister(;sep = sep, lprefix = length(bprefix)),bprefix)
        collect(iter)
    end
end

function Base.keys(d::LMDBDict{K}; prefix=UInt8[]) where K
    cursor_do(d, readonly = true) do cur
        collect(keys(cur,K,prefix=prefix))
    end
end

function Base.values(d::LMDBDict{K,V}; prefix=UInt8[]) where {K,V}
    cursor_do(d, readonly = true) do cur
        collect(values(cur,V,prefix=prefix))
    end
end

function Base.collect(d::LMDBDict{K,V}; prefix=UInt8[]) where {K,V}
    cursor_do(d, readonly = true) do cur
        collect((LMDB.LMDBIterator(cur,LMDB.ReturnBoth{K,V}(),Vector{UInt8}(prefix))))
    end
end

function valuesize(d::LMDBDict; prefix = UInt8[])
    cursor_do(d, readonly = true) do cur
        iter = LMDB.LMDBIterator(cur,LMDB.ReturnValueSize(),Vector{UInt8}(prefix))
        sum(iter)
    end
end

function Base.getindex(d::LMDBDict{K,V},k) where {K,V}
    cursor_do(d, readonly = true) do cur
        LMDB.get(cur, k, V, LMDB.MDB_SET_KEY)
    end
end

function Base.haskey(d::LMDBDict{K}, key) where K
    txn_dbi_do(d, readonly = true) do txn, dbi
        mdb_key_ref = Ref(MDBValue(toref(key)))
        mdb_val_ref = Ref(MDBValue())
        @show key
        # Get value
        ret = _mdb_get(txn.handle, dbi.handle, mdb_key_ref, mdb_val_ref)
        if ret == MDB_NOTFOUND
            return false
        elseif ret == Cint(0)
            return true
        else
            throw(LMDB.LMDBError(ret))
        end
    end
end

function Base.setindex!(d::LMDBDict{K,V},v,k) where {K,V}
    txn_dbi_do(d) do txn, dbi
        LMDB.put!(txn,dbi,k,v)
    end
    v
end

function Base.delete!(d::LMDBDict{K},k) where K
    txn_dbi_do(d) do txn, dbi
        LMDB.delete!(txn, dbi, k)
    end
    d
end


text:
module Format

using ..NESSie
using ..NESSie: reverseindex, unpack, vertexnormals
using JSON: json
using LightXML: XMLDocument, create_root, new_child, add_text, set_attribute

include("format/hmo.jl")
export readhmo, writehmo

include("format/mcsf.jl")
export readmcsf

include("format/msms.jl")
export readmsms

include("format/off.jl")
export readoff

include("format/pqr.jl")
export readpqr

include("format/skel.jl")
export writeskel

include("format/stl.jl")
export readstl, writestl

include("format/vtk.jl")
export writevtk

include("format/xml3d.jl")
export writexml3d_json, writexml3d_xml

end # module


text:
using DataStructures
using ProgressMeter: Progress, finish!, next!

test = """
    1163751742
    1381373672
    2136511328
    3694931569
    7463417111
    1319128137
    1359912421
    3125421639
    1293138521
    2311944581
    """

function parse_input(input::AbstractString)
    lines = split(input, r"\s+"; keepempty = false)
    return hcat([parse.(Int, split(line, "")) for line in lines]...) |>
        permutedims
end

function neighbors(M::Matrix, ci::CartesianIndex)::Vector{CartesianIndex}
    return [ci - CartesianIndex(1, 0),
            ci - CartesianIndex(0, 1),
            ci + CartesianIndex(1, 0),
            ci + CartesianIndex(0, 1)] ∩ keys(M)
end

function dijkstra(cave::Matrix, start, stop)::Int
    p = Progress(length(cave))
    handles = similar(cave)
    costs = fill(-1, size(cave))
    open = BinaryMinHeap([(0, start)])
    sizehint!(open, length(cave))
    while !isempty(open)
        (ci, i) = pop!(open)
        costs[i] < 0 || continue
        costs[i] = ci
        if i == stop
            finish!(p)
            return ci
        end
        for j in neighbors(cave, i)
            push!(open, (ci + cave[j], j))
        end
        next!(p)
    end
    @assert false "Couldn't reach $(stop)!"
end

q1(cave)::Int = dijkstra(cave, first(keys(cave)), last(keys(cave)))

function q2(cave)::Int
    big_cave = [(cave ⊕ 0) (cave ⊕ 1) (cave ⊕ 2) (cave ⊕ 3) (cave ⊕ 4)
                (cave ⊕ 1) (cave ⊕ 2) (cave ⊕ 3) (cave ⊕ 4) (cave ⊕ 5)
                (cave ⊕ 2) (cave ⊕ 3) (cave ⊕ 4) (cave ⊕ 5) (cave ⊕ 6)
                (cave ⊕ 3) (cave ⊕ 4) (cave ⊕ 5) (cave ⊕ 6) (cave ⊕ 7)
                (cave ⊕ 4) (cave ⊕ 5) (cave ⊕ 6) (cave ⊕ 7) (cave ⊕ 8)]
    return q1(big_cave)
end

(⊕)(M::Matrix, i::Int) = @. (M + i - 1) % 9 + 1

let cave = parse_input(test)
    @assert q1(cave) == 40
    @assert q2(cave) == 315
end

cave = parse_input(read("day15.in", String))
println("Q1: ", q1(cave))
println("Q2: ", q2(cave))       # slow!

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
"""
    Permute

Subtype of `IndexTransform` that permutes each index prior to accessing.
"""
struct Permute{P,I,NI,NO} <: IndexTransform{NI,NO}
    perm::P
    iperm::I

    function Permute(
        perm::Tuple{Vararg{CanonicalInt,NI}},
        iperm::Tuple{Vararg{CanonicalInt,NO}}
    ) where {NI,NO}
        new{typeof(perm),typeof(iperm),NI,NO}(perm, iperm)
    end
end

Base.invperm(x::Permute) = getfield(x, :iperm)
perm(x::Permute) = getfield(x, :perm)

const VecPermute = Permute{Tuple{StaticInt{2},StaticInt{1}},Tuple{StaticInt{1}},2,1}
const MatPermute = Permute{Tuple{StaticInt{2},StaticInt{1}},Tuple{StaticInt{2},StaticInt{1}},2,2}

@inline function (p::Permute{P,I})(x::Tuple{Vararg{Any}}) where {P,I}
    Static.permute(x, getfield(p, :perm))
end
(t::Permute)(i::AbstractCartesianIndex) = t(Tuple(i))
@inline function (t1::Permute{P1,I1})(t2::Permute{P2,I2}) where {P1,I1,P2,I2}
    Permute(Static.permute(perm(t2), perm(t1)), Static.permute(invperm(t2), invperm(t1)))
end

function Base.show(io::IO, ::MIME"text/plain", @nospecialize(x::Permute))
    print(io, "Permute(x.perm, x.iperm)")
end


text:
function instantiate(model_type::DataType, data, pools, param::Dict{Symbol, <:Any})
    debug(LOGGER, "Instantiating model of type $(model_type) with params $(param)")
    m = model_type(data, pools)
    set_param!(m, param)
    return m
end

function set_param!(model, param::Dict{Symbol, <:Any})
    has_changed = false
    for (k, v) in param
        k in keys(get_model_params(model)) || throw(ArgumentError("Model $(typeof(model)) has no parameter $k. Valid parameters are $(get_model_params(model))."))
        fieldtype = typeof(getfield(model, k))
        typeof(v) == fieldtype || throw(ArgumentError("Cannot set parameter $k for $(typeof(model)), because the parameter type is $(typeof(v)) but the model expects $(fieldtype)."))
        is_valid_param_value(model, Val{k}, v) || throw(ArgumentError("Parameter $k has value $v which is outside the valid range for model $(typeof(model))."))
        if getfield(model,k) != v
            has_changed = true
            setfield!(model, k, v)
        end
    end
    has_changed && invalidate_solution!(model)
    return nothing
end

function set_pools!(model, pools::Vector{Symbol})
    islabelenc(pools, learning_pool_enc) || throw(ArgumentError("Pools must be a LabelEncoding of $(learning_pool_enc.label)."))
    set_pools!(model, labelmap(pools))
    return nothing
end

function set_pools!(model, pools::Dict{Symbol, Vector{Int}})
    if !isdefined(model, :pools) ||model.pools != pools
        model.pools = pools
        invalidate_solution!(model)
    end
    return nothing
end

function set_adjust_K!(model, adjust_K::Bool)
    if adjust_K && !model.adjust_K
        model.adjust_K = adjust_K
        update_K!(model)
        invalidate_solution!(model)
    elseif !adjust_K && model.adjust_K
        model.adjust_K = adjust_K
        invalidate_solution!(model)
    end
    return nothing
end

is_K_adjusted(model)::Bool = model.adjust_K

function update_K!(model)
    updated_K = MLKernels.kernelmatrix(Val(:col), model.kernel_fct, model.data)
    if !isdefined(model, :K) || updated_K != model.K || (!isdefined(model, :K_adjusted) && model.adjust_K)
        debug(LOGGER, "[UPDATE_K] Updating Kernel matrix.")
        invalidate_solution!(model)
        model.K = updated_K
        if model.adjust_K
            model.K_adjusted = adjust_kernel_matrix(updated_K)
        end
    else
        debug(LOGGER, "[UPDATE_K] Nothing to update, kernel matrix has not changed.")
    end
    return nothing
end

function set_data!(model, data::Array{T, 2}) where T <: Real
    if model.data != data
        model.data = data
        update_K!(model)
        invalidate_solution!(model)
    end
    return nothing
end

function set_kernel!(model::OCClassifier, kernel_fct::Kernel)
    if !isdefined(model, :kernel_fct) || model.kernel_fct != kernel_fct
        model.kernel_fct = kernel_fct
        update_K!(model)
        invalidate_solution!(model)
    end
    return nothing
end

function initialize!(model::OCClassifier, strategy::InitializationStrategy)
    debug(LOGGER, "[INITIALIZE] Initializing $(typeof(model)) with strategy: $(typeof(strategy))")
    C, kernel_fct = get_parameters(model, strategy)
    debug(LOGGER, "[INITIALIZE] Fitted parameters are: C = $(C) and kernel_fct = $(kernel_fct)")
    set_C!(model, C)
    set_kernel!(model, kernel_fct)
    model.state = model_initialized
    return nothing
end

get_kernel(model::OCClassifier) = model.kernel_fct


text:
## DateString Aliasing

yrmo2dir(date::TimeType) = Dates.format(date,dateformat"yyyy/mm")
yrmo2str(date::TimeType) = Dates.format(date,dateformat"yyyymm")
yr2str(date::TimeType)   = Dates.format(date,dateformat"yyyy")
ymd2str(date::TimeType)  = Dates.format(date,dateformat"yyyymmdd")
mo2str(date::TimeType)   = Dates.format(date,dateformat"mm")

yrmo2str(yr::Integer,mo::Integer) = @sprintf("%04d%02d",yr,mo)
mo2str(mo::Integer) = @sprintf("%02d",mo)
dy2str(dy::Integer) = @sprintf("%02d",dy)

function extractdate(startdate::TimeType,finish::TimeType);

    yrs = Dates.year(start);  mos = Dates.month(start);  dys = Dates.day(start);
    yrf = Dates.year(finish); mof = Dates.month(finish); dyf = Dates.day(finish);
    ndy = Dates.value((finish-start)/Dates.day(1));
    dvecs = Date(yrs,mos); dvecf = Date(yrf,mof);

    dvec = collect(dvecs:Month(1):dvecf);

    return dvec,dys,dyf,ndy

end

bold() = Crayon(bold=true)
reset() = Crayon(reset=true)

function checkegrid(step::Real)

    if rem(360,step) == 0
          return true
    else; return false
    end

end


text:
struct Box
    x1::Int
    x2::Int
    y1::Int
    y2::Int
    z1::Int
    z2::Int
end

function volume((;x1, x2, y1, y2, z1, z2)::Box)
    return (x2 - x1 + 1) * (y2 - y1 + 1) * (z2 - z1 + 1)
end

function volume50((;x1, x2, y1, y2, z1, z2)::Box)
    if x2 < -50 || x1 > 50 || y2 < -50 || y1 > 50 || z2 < -50 || z1 > 50
        return 0
    end
    x1 = max(-50, x1)
    x2 = min(50, x2)
    y1 = max(-50, y1)
    y2 = min(50, y2)
    z1 = max(-50, z1)
    z2 = min(50, z2)
    return (x2 - x1 + 1) * (y2 - y1 + 1) * (z2 - z1 + 1)
end

function step(boxes, (on, (X1, X2, Y1, Y2, Z1, Z2)))
    clipped_boxes = Box[]
    if on
        push!(clipped_boxes, Box(X1, X2, Y1, Y2, Z1, Z2))
    end
    for box in boxes
        (;x1, x2, y1, y2, z1, z2) = box
        if X1 > x2 || X2 < x1 || Y1 > y2 || Y2 < y1 || Z1 > z2 || Z2 < z1
            push!(clipped_boxes, box)
            continue
        end
        if x1 < X1 <= x2
            push!(clipped_boxes, Box(x1, X1 - 1, y1, y2, z1, z2))
            x1 = X1
        end
        if x1 <= X2 < x2
            push!(clipped_boxes, Box(X2 + 1, x2, y1, y2, z1, z2))
            x2 = X2
        end
        if y1 < Y1 <= y2
            push!(clipped_boxes, Box(x1, x2, y1, Y1 - 1, z1, z2))
            y1 = Y1
        end
        if y1 <= Y2 < y2
            push!(clipped_boxes, Box(x1, x2, Y2 + 1, y2, z1, z2))
            y2 = Y2
        end
        if z1 < Z1 <= z2
            push!(clipped_boxes, Box(x1, x2, y1, y2, z1, Z1 - 1))
            z1 = Z1
        end
        if z1 <= Z2 < z2
            push!(clipped_boxes, Box(x1, x2, y1, y2, Z2 + 1, z2))
        end
    end
    return clipped_boxes
end

function parse_line(line)
    regex = r"x=(-?\d+)..(-?\d+),y=(-?\d+)..(-?\d+),z=(-?\d+)..(-?\d+)"
    coords = parse.(Int, match(regex, last(split(line))).captures)
    return startswith(line, "on"), coords
end

function part1(input)
    instructions = parse_line.(readlines(input))
    boxes = foldl(step, instructions, init = Box[])
    return sum(volume50.(boxes))
end

function part2(input)
    instructions = parse_line.(readlines(input))
    boxes = foldl(step, instructions, init = Box[])
    return sum(volume.(boxes))
end


text:
# diffs.jl
# finite differences
# 2019-03-96 Jeff Fessler, University of Michigan

using LinearMaps

"""
`d = diff2d_forw(X)`

2D finite differences along both dimensions, for anisotropic TV regularization.
Performs the same operations as
d = ``[(I_N \\otimes D_M); (D_N \\otimes I_M)] X[:]``
where ``D_N`` denotes the ``N-1 \\times N`` 1D finite difference matrix
and ``\\otimes`` denotes the Kronecker product,
but does it efficiently
without using `spdiagm` (or any `SparseArrays` function).

in
* `X`		`M x N` array (typically a 2D image).
It cannot be a Vector!  (But it can be a `Mx1` or `1xN` 2D array.)

out
* `d`		vector of length `N*(M-1) + (N-1)*M`
"""
function diff2d_forw(x::AbstractMatrix{<:Number})

return [diff(x,dims=1)[:]; diff(x,dims=2)[:]]
end

"""
`z =  diff2d_adj(d, M, N; out2d=false)`

Adjoint of 2D finite differences along both dimensions.
Performs the same operations as
``z = [(I_N \\otimes D_M); (D_N \\otimes I_M)]' * d``
where D_N denotes the N-1 x N 1D finite difference matrix
and \\otimes denotes the Kronecker product,
but does it efficiently without using spdiagm (or any SparseArrays function).

in
* `d`		vector of length `N*(M-1) + (N-1)*M`
* `M,N`		desired output size

option
* `out2d`	if true then return `M x N` array, else `M*N` vector

out
* `z`		`M*N` vector or `M x N` array (typically a 2D image)

"""
function diff2d_adj(d::AbstractVector{<:Number}, M, N; out2d=false)

@assert length(d) == N*(M-1) + (N-1)*M

dx = d[1:(N*(M-1))]
dx = reshape(dx, M-1, N)
dy = d[1+(N*(M-1)):end]
dy = reshape(dy, M, N-1)
zx = [-transpose(dx[1,:]);
	(@views dx[1:end-1,:] - dx[2:end,:]);
	transpose(dx[end,:])]
zy = [-dy[:,1] (@views dy[:,1:end-1] - dy[:,2:end]) dy[:,end]]
z = zx + zy # M by N

if out2d # 2D array
	return z
else # vector
	return z[:]
end

end

"""
`T = diff2d_map(M,N)`
"""
function diff2d_map(M,N)
	return LinearMap(
        x -> diff2d_forw(reshape(x,M,N)),
        d -> diff2d_adj(d, M, N),
        N*(M-1)+M*(N-1), N*M)
end

"""
`T = diff_map(M,N)`

in
* `M,N` image size

out
* `T` a `LinearMap` object for regularizing via T*x
"""
function diff_map(M,N)
	return diff2d_map(M,N)
end

"""
`diff_map(:test)` self test
"""
function diff_map(test::Symbol)
	@assert test == :test
	M,N = 4,5
	T = diff_map(M,N)
	@assert Matrix(T)' == Matrix(T') # adjoint test
	true
end


text:
# test finite difference derivative with respect to the mesh vertex coordinates

"""
  Test complex stepping the derivative of the residual with repsect to the
  mesh coordinate nodes
"""
function test_coordsfd()

  println("testing coordinates finite difference")

  @testset "Testing Coordinate Field Derivative" begin
    fname = "input_vals_coordsfd.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    func = createFunctional(mesh, sbp, eqn, opts, "lift",  [1, 2, 3])

    testMeshDerivatives(mesh, sbp, opts)
    testResidualDerivatives(mesh, sbp, eqn, opts)
    testDJDx(mesh, sbp, eqn, opts, func)

    func = createFunctional(mesh, sbp, eqn, opts, "entropydissipation", Int[])
    testDJDx(mesh, sbp, eqn, opts, func)
  end

end

add_func1!(EulerTests, test_coordsfd, [TAG_SHORTTEST, TAG_COMPLEX, TAG_ADJOINT])


text:
using OrdinaryDiffEq
using Trixi

equations = CompressibleEulerEquations2D(1.4) # set gas gamma = 1.4

# initial condition
initial_condition = initial_condition_convergence_test
source_term       = source_terms_convergence_test

# boundary condition types
boundary_condition_dirichlet= BoundaryConditionDirichlet(initial_condition)

# boundary condition dictionary
boundary_conditions = Dict( :OuterCircle => boundary_condition_dirichlet,
                            :LeftSlant => boundary_condition_dirichlet,
                            :RightSlant => boundary_condition_dirichlet,
                            :IceCream => boundary_condition_dirichlet)

# DGSEM solver.
#    1) polydeg must be >= the polynomial order set in the HOHQMesh control file to guarantee
#       freestream preservation. As a extra task try setting polydeg=3
#    2) VolumeIntegralFluxDifferencing with central volume flux is activated
#       for dealiasing
volume_flux = flux_ranocha
solver = DGSEM(polydeg=6, surface_flux=flux_hll,
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

# create the unstructured from your mesh file
mesh_file = joinpath(@__DIR__, "out/mesh_for_live_tutorial.mesh")
mesh = UnstructuredMesh2D(mesh_file)

# Create semidiscretization with all spatial discretization-related components
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions=boundary_conditions,
                                    source_terms=source_term)

# Create ODE problem from semidiscretization with time span from 0.0 to 2.0
tspan = (0.0, 5.0)
ode = semidiscretize(semi, tspan)

# Create the callbacks to output timing information, solution files, and adapt the time step
summary_callback = SummaryCallback()
save_solution = SaveSolutionCallback(interval=10,
                                     save_initial_solution=true,
                                     save_final_solution=true)
stepsize_callback = StepsizeCallback(cfl=1.0)

callbacks = CallbackSet(summary_callback, save_solution, stepsize_callback)

# Evolve ODE problem in time using `solve` from OrdinaryDiffEq
sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
# print the timer summary
summary_callback()


text:
@testset "Round2" begin
  @testset "Qt" begin
    qt, t = RationalFunctionField(QQ, "t")
    qtx, x = PolynomialRing(qt, "x")
    f = x^4 + t*x^3 - 6*x^2 - t*x + 1
    F, a = FunctionField(f, "a")
    O = integral_closure(Hecke.Localization(qt, degree), F)
    b = basis(O, F)
    mp = map(minpoly, b)
    @test all(i->iszero(mp[i](b[i])), 1:length(b))

    integral_closure(parent(numerator(t)), F)
    #= crashes julia in 1.6
    integral_closure(Hecke.Globals.Zx, F)
    factor(F.pol, F)
    =#
  end

  @testset "FldNum" begin
    R = Hecke.Localization(ZZ, 13)
    K, _ = wildanger_field(3, 9*13^2)
    O = integral_closure(R, K)
    b = basis(O, K)
    mp = map(minpoly, b)
    @test all(i->iszero(mp[i](b[i])), 1:length(b))
  end
end


text:
# SudokuService
# Sudoku-as-a-service
# By Iain Dunning
# https://github.com/IainNZ/SudokuService
# http://iaindunning.com/2013/sudoku-as-a-service.html
# License: MIT

# sudoku.jl
# A sudoku solver that uses a MIP to find solutions

# We have binary variables x[i,j,k] which, if = 1, say that cell (i,j)
# contains the number k
# The constraints are:
# 1 - Each row contains each number exactly once
# 2 - Each column contains each number exactly once
# 3 - Each 3x3 subgrid contains each number exactly once
# and the obvious one that
# 4 - Each cell has one value only
# We take as input a 9x9 matrix, where 0s are "blanks"

using JuMP

function SolveModel(initgrid)
  m = Model(:Max)  # Feasibility problem, so sense not important

  @defVar(m, 0 <= x[1:9, 1:9, 1:9] <= 1, Int)

  for val in 1:9
    for dim1 in 1:9
      # Constraint 1 - Each row...
      @addConstraint(m, sum{x[dim1, dim2, val], dim2=1:9} == 1)
      # Constraint 2 - Each column...
      @addConstraint(m, sum{x[dim2, dim1, val], dim2=1:9} == 1)
    end
  end

  # Constraint 3 - Each sub-grid...
  for i in 0:3:6
    for j in 0:3:6
      for val in 1:9
        @addConstraint(m, sum{x[i+off1, j+off2, val], off1=1:3, off2=1:3} == 1)
      end
    end
  end

  # Constraint 4 - Cells...
  for row in 1:9
    for col in 1:9
      @addConstraint(m, sum{x[row, col, val], val=1:9} == 1)
    end
  end

  # Initial solution
  for row in 1:9
    for col in 1:9
      if initgrid[row,col] != 0
        @addConstraint(m, x[row, col, initgrid[row, col]] == 1)
      end
    end
  end

  # Solve it
  status = solve(m)

  # Check solution
  if status == :Infeasible
    error("No solution found!")
  else
    mipSol = getValue(x)
    sol = zeros(Int,9,9)
    for row in 1:9
      for col in 1:9
        for val in 1:9
          if mipSol[row, col, val] >= 0.9
            sol[row, col] = val
          end
        end
      end
    end
    return sol
  end

end

max_stars_repo_name: IainNZ/SudokuService

max_issues_repo_name: IainNZ/SudokuService

max_forks_repo_name: IainNZ/SudokuService


text:
abstract type GeneralizedKroneckerProduct{T} <: AbstractMatrix{T} end

Base.eltype(K::GeneralizedKroneckerProduct{T}) where {T} = T

abstract type AbstractKroneckerProduct{T} <: GeneralizedKroneckerProduct{T} end

Base.IndexStyle(::Type{<:GeneralizedKroneckerProduct}) = IndexCartesian()

"""
    KroneckerProduct{T,TA<:AbstractMatrix, TB<:AbstractMatrix} <: AbstractKroneckerProduct{T}

Concrete Kronecker product between two matrices `A` and `B`.
"""
struct KroneckerProduct{T<:Any,TA<:AbstractMatrix, TB<:AbstractMatrix} <: AbstractKroneckerProduct{T}
    A::TA
    B::TB
    function KroneckerProduct(A::AbstractMatrix{T}, B::AbstractMatrix{V}) where {T, V}
        return new{promote_type(T, V), typeof(A), typeof(B)}(A, B)
    end
end

"""
    kronecker(A::AbstractMatrix, B::AbstractMatrix)

Construct a Kronecker product object between two arrays. Does not evaluate the
Kronecker product explictly.
"""
kronecker(A::AbstractMatrix, B::AbstractMatrix) = KroneckerProduct(A, B)

"""
    kronecker(A::AbstractMatrix, B::AbstractMatrix)

Higher-order Kronecker lazy kronecker product, e.g.
```
kronecker(A, B, C, D)
```
"""
kronecker(A::AbstractMatrix, B::AbstractMatrix...) = kronecker(A,
                                                            kronecker(B...))

"""
    ⊗(A::AbstractMatrix, B::AbstractMatrix)

Binary operator for `kronecker`, computes as Lazy Kronecker product. See
`kronecker` for documentation.
"""
⊗(A::AbstractMatrix, B::AbstractMatrix) = kronecker(A, B)
⊗(A::AbstractMatrix...) = kronecker(A...)

"""
    getindex(K::AbstractKroneckerProduct, i1::Integer, i2::Integer)

Computes and returns the (i,j)-th element of an `AbstractKroneckerProduct` K.
Uses recursion if `K` is of an order greater than two.
"""
function getindex(K::AbstractKroneckerProduct, i1::Integer, i2::Integer)
    A, B = getmatrices(K)
    m, n = size(A)
    k, l = size(B)
    return A[cld(i1, k), cld(i2, l)] * B[(i1 - 1) % k + 1, (i2 - 1) % l + 1]
end

"""
    getmatrices(K::AbstractKroneckerProduct)

Obtain the two matrices of an `AbstractKroneckerProduct` object.
"""
getmatrices(K::AbstractKroneckerProduct) = (K.A, K.B)

"""
    getmatrices(A::AbstractArray)

Returns a matrix itself. Needed for recursion.
"""
getmatrices(A::AbstractArray) = (A,)

"""
    size(K::AbstractKroneckerProduct)

Returns a the size of an `AbstractKroneckerProduct` instance.
"""
function size(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    (m, n) = size(A)
    (k, l) = size(B)
    return m * k, n * l
end

"""
    size(K::GeneralizedKroneckerProduct)

Returns a the size of an `GeneralizedKroneckerProduct` instance.
"""
size(K::GeneralizedKroneckerProduct, dim::Integer) = size(K)[dim]

# CHECKS

"""
    issquare(A::AbstractMatrix)

Checks if an array is a square matrix.
"""
function issquare(A::AbstractMatrix)
    m, n = size(A)
    return m == n
end

"""
    issquare(K::AbstractKroneckerProduct)

Checks if all matrices of a Kronecker product are square.
"""
issquare(K::AbstractKroneckerProduct) = issquare(K.A) && issquare(K.B)

squarecheck(K::AbstractKroneckerProduct) = issquare(K) || throw(
            DimensionMismatch(
                "kronecker system is not composed of two square matrices: " *
                                        "$size(K.A) and $size(K.B)"))

"""
    issymmetric(K::AbstractKroneckerProduct)

Checks if a Kronecker product is symmetric.
"""
function LinearAlgebra.issymmetric(K::AbstractKroneckerProduct)
    return squarecheck(K) && issymmetric(K.A) && issymmetric(K.B)
end

"""
    isposdef(K::AbstractKroneckerProduct)

Test whether a Kronecker product is positive definite (and Hermitian) by trying
to perform a Cholesky factorization of K.
"""
function LinearAlgebra.isposdef(K::AbstractKroneckerProduct)
    squarecheck(K)
    return isposdef(K.A) && isposdef(K.B)
end

# SCALAR PROPERTIES

"""
    order(M::AbstractMatrix)

Returns the order of a matrix, i.e. how many matrices are involved in the
Kronecker product (default to 1 for general matrices).
"""
order(M::AbstractMatrix) = 1
order(M::AbstractKroneckerProduct) = order(M.A) + order(M.B)

"""
    det(K::AbstractKroneckerProduct)

Compute the determinant of a Kronecker product.
"""
function LinearAlgebra.det(K::AbstractKroneckerProduct)
    squarecheck(K)
    A, B = getmatrices(K)
    m = size(A)[1]
    n = size(B)[1]
    return det(K.A)^n * det(K.B)^m
end

"""
    logdet(K::AbstractKroneckerProduct)

Compute the logarithm of the determinant of a Kronecker product.
"""
function LinearAlgebra.logdet(K::AbstractKroneckerProduct)
    squarecheck(K)
    A, B = getmatrices(K)
    m = size(A)[1]
    n = size(B)[1]
    return n * logdet(A) + m * logdet(B)
end

"""
    tr(K::AbstractKroneckerProduct)

Compute the trace of a Kronecker product.
"""
function LinearAlgebra.tr(K::AbstractKroneckerProduct)
    squarecheck(K)
    return tr(K.A) * tr(K.B)
end

"""
    inv(K::AbstractKroneckerProduct)

Compute the inverse of a Kronecker product.
"""
function inv(K::AbstractKroneckerProduct)
    squarecheck(K)
    A, B = getmatrices(K)
    return KroneckerProduct(inv(A), inv(B))
end

"""
    adjoint(K::AbstractKroneckerProduct)

Compute the adjoint of a Kronecker product.
"""
function adjoint(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return kronecker(A', B')
end

"""
    transpose(K::AbstractKroneckerProduct)

Compute the transpose of a Kronecker product.
"""
function Base.transpose(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return kronecker(transpose(A), transpose(B))
end

function Base.conj(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return kronecker(conj(A), conj(B))
end

# COLLECTING

"""
    collect(K::AbstractKroneckerProduct)

Collects a lazy instance of the `AbstractKroneckerProduct` type into a full,
native matrix. Equivalent with `Matrix(K::AbstractKroneckerProduct)`.
"""
function collect(K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    return kron(A, B)
end

"""
    Matrix(K::GeneralizedKroneckerProduct)

Converts a `GeneralizedKroneckerProduct` instance to a Matrix type.
"""
Base.Matrix(K::GeneralizedKroneckerProduct) = collect(K)

function Base.:+(A::AbstractKroneckerProduct, B::StridedMatrix)
    C = Matrix(A)
    C .+= B
    return C
end
Base.:+(A::StridedMatrix, B::AbstractKroneckerProduct) = B + A

function Base.kron(K::AbstractKroneckerProduct, C::AbstractMatrix)
    A, B = getmatrices(K)
    return kron(kron(A, B), C)
end

function Base.kron(A::AbstractMatrix, K::AbstractKroneckerProduct)
    B, C = getmatrices(K)
    return kron(A, kron(B, C))
end

Base.kron(K1::AbstractKroneckerProduct,
            K2::AbstractKroneckerProduct) = kron(collect(K1), collect(K2))

# mixed-product property
function Base.:*(K1::AbstractKroneckerProduct, K2::AbstractKroneckerProduct)
    A, B = getmatrices(K1)
    C, D = getmatrices(K2)
    # check for size
    if size(A, 2) != size(C, 1)
        throw(DimensionMismatch("Mismatch between A and C in (A ⊗ B)(C ⊗ D)"))
    end
    if size(B, 2) != size(D, 1)
        throw(DimensionMismatch("Mismatch between B and D in (A ⊗ B)(C ⊗ D)"))
    end
    return (A * C) ⊗ (B * D)
end

"""
    lmul!(a::Number, K::AbstractKroneckerProduct)

Scale an `AbstractKroneckerProduct` `K` inplace by a factor `a` by rescaling the
**left** matrix.
"""
function LinearAlgebra.lmul!(a::Number, K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    lmul!(a, A)
end

"""
    rmul!(K::AbstractKroneckerProduct, a::Number)

Scale an `AbstractKroneckerProduct` `K` inplace by a factor `a` by rescaling the
**right** matrix.
"""
function LinearAlgebra.rmul!(K::AbstractKroneckerProduct, a::Number)
    A, B = getmatrices(K)
    rmul!(B, a)
end

function Base.:*(a::Number, K::AbstractKroneckerProduct)
    A, B = getmatrices(K)
    kronecker(a * A, B)
end

function Base.:*(K::AbstractKroneckerProduct, a::Number)
    A, B = getmatrices(K)
    kronecker(A, B * a)
end

# SOLVING
function LinearAlgebra.:\(K::AbstractKroneckerProduct{T}, c::AbstractVector{T}) where {T}
    size(K, 1) != length(c) && throw(DimensionMismatch("size(K, 1) != length(c)"))
    C = reshape(c, size(K.B, 1), size(K.A, 1)) # matricify
    return vec((K.B \ C) / K.A') #(A ⊗ B)vec(X) = vec(C) <=> BXA' = C => X = B^{-1} C A'^{-1}
end