text:
####SAMPLING FUNCTIONS####
#function for a Distributed worker to produce a set of samples of given parameters from genomic sequences
function get_sample_set(input_sample_jobs::RemoteChannel, completed_sample_jobs::RemoteChannel, progress_updates::RemoteChannel)
    while isready(input_sample_jobs)
        genome_path, genome_index_path, partition_df, partitionid, sample_set_length, sample_window_min, sample_window_max, deterministic = take!(input_sample_jobs)

        stranded::Bool = get_strand_dict()[partitionid]
        scaffold_sequence_record_dict::Dict{String,LongSequence} = build_scaffold_seq_dict(genome_path, genome_index_path)

        sample_df = DataFrame(SampleScaffold=String[],SampleStart=Integer[],SampleEnd=Integer[],SampleSequence=LongSequence[],Strand=Char[])
        metacoordinate_bitarray = trues(partition_df.MetaEnd[end])
        sample_set_counter = 0

        while sample_set_counter < sample_set_length #while we don't yet have enough sample sequence
            sample_scaffold::String, sample_Start::Integer, sample_End::Integer, sample_metaStart::Integer, sample_metaEnd::Integer, sample_sequence::LongSequence, strand::Char = get_sample(metacoordinate_bitarray, sample_window_min, sample_window_max,  partition_df, scaffold_sequence_record_dict; stranded=stranded, deterministic=deterministic)
            push!(sample_df,[sample_scaffold, sample_Start, sample_End, sample_sequence, strand]) #push the sample to the df
            sample_length = sample_End - sample_Start + 1
            sample_set_counter += sample_length #increase the counter by the length of the sampled sequence
            metacoordinate_bitarray[sample_metaStart:sample_metaEnd] = [false for base in 1:sample_length] #mark these residues as sampled
            put!(progress_updates, (partitionid, min(sample_set_counter,sample_set_length)))
        end
        put!(completed_sample_jobs,(partitionid, sample_df))
    end
end
                #get_sample_set() SUBFUNCTIONS
                #function defining whether partitions respect stranding upon fetching sequence (ie is the sequence fetched in the feature strand orientation, or are we agnostic about the strand we sample?)
                function get_strand_dict()
                    return Dict("exon"=>true,"periexonic"=>true,"intergenic"=>false)
                end
                #function to obtain a dict of scaffold sequences from a FASTA reader
                function build_scaffold_seq_dict(genome_fa, genome_index)
                    genome_reader = open(FASTA.Reader, genome_fa, index=genome_index)
                    seq_dict::Dict{String, LongSequence} = Dict{String,FASTA.Record}()
                    @inbounds for record in genome_reader
                        id = FASTA.identifier(record)
                        seq_dict[rectify_identifier(id)]=FASTA.sequence(record)
                    end
                    close(genome_reader)
                    return seq_dict
                end

                # function to convert scaffold ID from that observed by the masked .fna to the more legible one observed by the GRCz11 GFF3
                function rectify_identifier(scaffold_id::String)
                    if length(scaffold_id) >= 4 && scaffold_id[1:4] == "CM00" #marks chromosome scaffold
                        chr_code = scaffold_id[5:10]
                        chr_no = "$(Int((parse(Float64,chr_code)) - 2884.2))"
                        return chr_no
                    else
                        return scaffold_id
                    end
                end

                #function to produce a single sample from a metacoordinate set and the feature df
                function get_sample(metacoordinate_bitarray::BitArray, sample_window_min::Integer, sample_window_max::Integer, partition_df::DataFrame,  scaffold_seq_dict::Dict{String,LongSequence}; stranded::Bool=false, deterministic::Bool=false)
                    proposal_acceptance = false
                    sample_metaStart = 0
                    sample_metaEnd = 0
                    sample_Start = 0
                    sample_End = 0
                    sample_sequence = LongSequence{DNAAlphabet{2}}("")
                    sample_scaffold = ""
                    strand = nothing

                    while proposal_acceptance == false
                        available_indices = findall(metacoordinate_bitarray) #find all unsampled indices
                        window = "FAIL"
                        start_index,feature_metaStart,feature_metaEnd, feature_length = 0,0,0,0
                        strand = '0'
                        while window == "FAIL"
                            start_index = rand(available_indices) #randomly choose from the unsampled indices
                            feature_metaStart, feature_metaEnd, strand = get_feature_params_from_metacoord(start_index, partition_df, stranded)
                            feature_length = length(feature_metaStart:feature_metaEnd) #find the metaboundaries of the feature the index occurs in
                            if feature_length >= sample_window_min #don't bother finding windows on features smaller than min
                                window = determine_sample_window(feature_metaStart, feature_metaEnd, start_index, metacoordinate_bitarray, sample_window_min, sample_window_max) #get an appropriate sampling window around the selected index, given the feature boundaries and params
                            end
                        end

                        sample_scaffoldid, sample_scaffold_start, sample_scaffold_end = meta_to_feature_coord(window[1],window[2],partition_df)

                        proposal_sequence = fetch_sequence(sample_scaffoldid, scaffold_seq_dict, sample_scaffold_start, sample_scaffold_end, strand; deterministic=deterministic) #get the sequence associated with the sample window

                        if mask_check(proposal_sequence)
                            proposal_acceptance = true #if the sequence passes the mask check, accept the proposed sample
                            sample_scaffold = sample_scaffoldid
                            sample_Start = sample_scaffold_start
                            sample_End = sample_scaffold_end
                            sample_metaStart = window[1]
                            sample_metaEnd = window[2]
                            sample_sequence=proposal_sequence
                        end
                    end
                    return sample_scaffold, sample_Start, sample_End, sample_metaStart, sample_metaEnd, sample_sequence, strand
                end

                #function to find a valid sampling window
                function determine_sample_window(feature_metaStart::Integer, feature_metaEnd::Integer, metacoord::Integer, metacoord_bitarray::BitArray, sample_window_min::Integer, sample_window_max::Integer)
                    window_start = 0
                    window_end = 0
                    feature_size = feature_metaEnd - feature_metaStart + 1
                    feature_sampled = any(!eval, metacoord_bitarray[feature_metaStart:feature_metaEnd])
                    #attempt to construct the sample window by taking the whole feature
                    if !feature_sampled && feature_size < sample_window_max && feature_size > sample_window_min
                        window_start = feature_metaStart
                        window_end = feature_metaEnd
                        return window_start, window_end
                    else #if this fails, build the biggest window we can from the sampling point, within the feature
                        featurepos = metacoord - feature_metaStart + 1
                        next_sampled_index = findnext(!eval, metacoord_bitarray[feature_metaStart:feature_metaEnd], featurepos)
                        prev_sampled_index = findprev(!eval, metacoord_bitarray[feature_metaStart:feature_metaEnd], featurepos)
                        next_sampled_index === nothing ? (window_end = feature_metaEnd) : (window_end = next_sampled_index + feature_metaStart - 1)
                        prev_sampled_index === nothing ? (window_start = feature_metaStart) : (window_start = prev_sampled_index + feature_metaStart - 1)
                        windowsize = window_end - window_start + 1
                        #check to see if this window is bigger than the min, if not, return a failure code
                        windowsize < sample_window_min && return "FAIL"
                        #check to see if this window is bigger than the max, if so, trim it before returning, removing as evenly as possible around the metacoordinate
                        if windowsize > sample_window_max
                            bases_to_trim = windowsize - sample_window_max
                            clearance_5P = metacoord - window_start + 1
                            clearance_3P = window_end - metacoord + 1
                            trimmed=0
                            while trimmed < bases_to_trim
                                if clearance_5P >= clearance_3P && clearance_5P > 0
                                    clearance_5P -= 1
                                elseif clearance_3P > clearance_5P && clearance_3P > 0
                                    clearance_3P -= 1
                                end
                                trimmed +=1
                            end
                            window_start = metacoord - clearance_5P + 1
                            window_end = metacoord + clearance_3P - 1
                        end
                        return window_start, window_end
                    end
                end

                # function to check for repetitive stretches or degenerate bases in proposal sequence
                function mask_check(proposal_sequence::LongSequence)
                    proposal_acceptance = true
                    if hasambiguity(proposal_sequence) || isrepetitive(proposal_sequence, (length(proposal_sequence) ÷ 10))
                        proposal_acceptance = false
                    end
                    return proposal_acceptance
                end

####SHARED SEQUENCE FETCHER####
# function to get proposal sequence from dict of scaffold sequences, given coords and scaffold id
function fetch_sequence(scaffold_id::String, scaffold_seq_dict::Dict{String, LongSequence}, proposal_start::Integer, proposal_end::Integer, strand::Char; deterministic=false)
    if strand == '0' #unstranded samples may be returned with no preference in either orientation
        deterministic ? strand = '+' : (rand(1)[1] <=.5 ? strand = '+' : strand = '-')
    end
    if strand == '+'
        proposal_sequence = scaffold_seq_dict[scaffold_id][proposal_start:proposal_end]
    elseif strand == '-'
        proposal_sequence = reverse_complement(scaffold_seq_dict[scaffold_id][proposal_start:proposal_end])
    else
        throw(ArgumentError("Invalid sample code! Must be '+', '-', or '0' (random strand)"))
    end

    return proposal_sequence
end

####SHARED BASIC COORDINATE SUBFUNCTIONS####
# function to push scaffold ID, start, and end points of given featuretype to supplied dataframe
function build_feature_df!(GFF3_path::String, feature_type::String, scaffold_exclusion::String, feature_df::DataFrame)
    reader = open(GFF3.Reader, GFF3_path) # access the GFF3
    @inbounds for record in reader # iterate over Gff3 records
        if GFF3.isfeature(record) # if the record is a feature
            if GFF3.featuretype(record) == feature_type #if the features is of the requested type, get the following info
                seqID = GFF3.seqid(record)
                if seqID != scaffold_exclusion
                    seq_start = GFF3.seqstart(record)
                    seq_end = GFF3.seqend(record)
                    if feature_type == "CDS" || feature_type == "gene"
                        seq_strand = convert(Char, GFF3.strand(record))
                        push!(feature_df, [seqID, seq_start, seq_end, seq_strand])
                    else
                        push!(feature_df, [seqID, seq_start, seq_end]) # push relevant info to the df
                    end
                end
            end
        end
    end
    close(reader)
end

#function to assemble dataframe of scaffold coords + metacoords given gff3
function build_scaffold_df(gff3_path)
    scaffold_df = DataFrame(SeqID = String[], Start = Integer[], End = Integer[])
    build_feature_df!(gff3_path, "supercontig", "MT", scaffold_df)
    build_feature_df!(gff3_path, "chromosome", "MT", scaffold_df)
    add_metacoordinates!(scaffold_df)
    return scaffold_df
end

#function to add pad to either side of some featureset
function add_pad_to_coordinates!(feature_df::DataFrame, scaffold_df::DataFrame, pad_size::Integer; col_symbols::Array{Symbol}=[:Start, :End])
    pad_start_array = zeros(Integer,size(feature_df,1))
    pad_end_array = zeros(Integer,size(feature_df,1))
    feature_df[!, col_symbols[1]] = [max(feature_df.Start[i]-pad_size,1) for i in 1:size(feature_df,1)] #truncate pads at beginning and end of scaffolds
    feature_df[!, col_symbols[2]] = [min(feature_df.End[i]+pad_size,scaffold_df.End[findfirst(isequal(feature_df.SeqID[i]),scaffold_df.SeqID)]) for i in 1:size(feature_df,1)]
end

####SHARED METACOORDINATE FUNCTIONS####
# function to add a metacoordinate column to a dataframe of scaffold positions, allowing sampling across all scaffolds
function add_metacoordinates!(feature_df::DataFrame)
    meta_start_array = Integer[]
    meta_end_array = Integer[]
    metaposition = 1
    @inbounds for feature in eachrow(feature_df)
        push!(meta_start_array, (metaposition))
        metaposition += feature.End - feature.Start
        push!(meta_end_array, (metaposition))
        metaposition += 1
    end
    feature_df.MetaStart = meta_start_array # metacoordinate contains 'start' metaposition across all genomic material
    feature_df.MetaEnd = meta_end_array # metacoordinate contains 'end' metaposition across all genomic material
end

# function to convert metacoordinate set to scaffold-relative coordinates
function meta_to_feature_coord(meta_start::Integer, meta_end::Integer, feature_df::DataFrame)
    feature_row = get_feature_row_index(feature_df, meta_start)
    seqid = feature_df.SeqID[feature_row]
    scaffold_start = feature_df.Start[feature_row] + (meta_start - feature_df.MetaStart[feature_row])
    scaffold_end = feature_df.End[feature_row] - (feature_df.MetaEnd[feature_row] - meta_end)
    return seqid, scaffold_start, scaffold_end
end

#function to obtain the feature boundaries and strand of the feature that a metacoordinate falls within
function get_feature_params_from_metacoord(metacoordinate::Integer, feature_df::DataFrame, stranded::Bool)
    feature_row = get_feature_row_index(feature_df, metacoordinate)
    feature_metaStart = feature_df.MetaStart[feature_row]
    feature_metaEnd = feature_df.MetaEnd[feature_row]
    stranded ? feature_strand = feature_df.Strand[feature_row] : feature_strand = '0'
    return feature_metaStart, feature_metaEnd, feature_strand
end

#function obtain the index of a feature given its metacoordinate
function get_feature_row_index(feature_df::DataFrame, metacoordinate::Integer)
    total_feature_bases = feature_df.MetaEnd[end]
    if metacoordinate == total_feature_bases #edge case of metacoord at end of range
        return size(feature_df,1) #last index in df
    else
        index = findfirst(end_coord->end_coord>metacoordinate,feature_df.MetaEnd) #find the index of the feature whose end metacoord is > the query metacoord
        if index > 1 && metacoordinate == feature_df.MetaEnd[index-1] #if the metacoordinate is the last base of the previous feature
            if metacoordinate >= feature_df.MetaStart[index-1] && metacoordinate <= feature_df.MetaEnd[index-1] #confirm that the metacoord is in the previous feature
                return index-1 #return the previous feature index
            end
        elseif metacoordinate >= feature_df.MetaStart[index] && metacoordinate <= feature_df.MetaEnd[index] #else confirm that the metacoordinate is in the found feature
            return index #return the found feature index
        else
            throw(DomainError("Unexpected metacoordinate $metacoordinate in partition of $total_feature_bases bases, with feature start $(feature_df.MetaStart[index]), end $(feature_df.MetaEnd[index])"))
        end
    end
end


text:
# DSL
# ===

# Types and methods that make up the Pseusoseq DSL API.

struct Amplifier{F<:Function}
    pred::F
    n::Int
end
amplify(n::Int) = Amplifier(x -> true, n)
amplify(pref::Function, n::Int) = Amplifier(pred, n)
(a::Amplifier)(p::Molecules) = amplify(a.pred, p, a.n)

struct Fragmenter
    meansize::Int
end

"""
    fragment(meansize::SequenceLength)

Construct a Fragmenter transform that fragments any `Molecules` passed to it, to
an average length of `meansize`.
"""
fragment(meansize::SequenceLength) = Fragmenter(meansize.val)
(f::Fragmenter)(p::Molecules) = fragment(p, f.meansize)

struct Tagger
    ntags::Int
end
tag(ntags::Int) = Tagger(ntags)
(t::Tagger)(p::Molecules) = tag(p, t.ntags)

struct NSubSampler
    nsamples::Int
end
subsample(n::Int) = NSubSampler(n)
(s::NSubSampler)(p::Molecules) = subsample(p, s.nsamples)

struct CovSubSampler{T<:Union{SequenceLength,Tuple{SequenceLength,SequenceLength}}}
    cov::ExpectedCoverage
    readlen::T
end
subsample(cov::ExpectedCoverage, rlens)  = CovSubSampler{typeof(rlens)}(cov, rlens)
(s::CovSubSampler{T})(p::Molecules) where {T} = subsample(p, s.cov, s.readlen)

struct Selector{F<:Function}
    f::F
end
select(f::Function) = Selector(f)
(s::Selector{F})(p::Molecules) where {F<:Function} = select(s.f, p)

makereads(len::SequenceLength) = unpaired_reads(len.val)
makereads() = unpaired_reads(nothing)
makereads(lena::SequenceLength, lenb::SequenceLength) = paired_reads(lena.val, lenb.val)
makereads(lens::Tuple{SequenceLength,SequenceLength}) = makereads(lens...)

struct PairedReads
    flen::Int
    rlen::Int
end
paired_reads(flen::Int, rlen::Int = flen) = PairedReads(flen, rlen)
(pr::PairedReads)(p::Molecules) = paired_reads(p, pr.flen, pr.rlen)

struct UnPairedReads{T<:Union{Nothing,Int}}
    len::T
end
unpaired_reads(len::T) where {T<:Union{Nothing,Int}} = UnPairedReads(len)
(ur::UnPairedReads)(p::Molecules) = unpaired_reads(p, ur.len)

struct SubstitutionMaker{F<:Function}
    fun::F
end
make_substitutions(f::F) where {F<:Function} = SubstitutionMaker{F}(f)
(sm::SubstitutionMaker{F})(p::Reads) where {F<:Function} = edit_substitutions(sm.fun, p)

struct FileGenerator
    filename::String
end
generate(filename::String) = FileGenerator(filename)
(fg::FileGenerator)(r::Reads) = generate(fg.filename, r)

struct FilePairGenerator
    F1::String
    F2::String
end
generate(F1::String, F2::String) = FilePairGenerator(F1, F2)
(fg::FilePairGenerator)(r::Reads) = generate(fg.F1, fg.F2, r)

ext: jl


text:
abstract type AbstractEnvironment end

"Move to the next step of an environment"
function step! end

"Reset an environment"
function reset! end

"Render the current environment"
function render end

"Close an environment"
function close! end

"Set the random seed of an environment"
function seed! end

"get the current state of an environment"
function get_state end

function is_end end

function get_reward end


text:
readintegers() = [parse(Int, x) for x in split(readline())]

n, m = readintegers()
votes = readintegers()

votes = sort!(votes, rev=true)
threshold = 1/(4m) * sum(votes)
sol = votes[m] >= threshold ? "Yes" : "No"

print(sol)

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using SafeTestsets

@safetestset "Estimating alpha diversity" begin include("test_alpha.jl") end

@safetestset "Estimating beta diversity" begin include("test_beta.jl") end

@safetestset "Generating LV model" begin include("test_lv.jl") end

@safetestset "Plotting time series" begin include("test_plot.jl") end

@safetestset "Transforming assays" begin include("test_transform.jl") end

@safetestset "Retrieving artifacts" begin include("test_artifacts.jl") end


text:
include("types.jl")

include("Header.jl")
include("Data.jl")

include("Rinex2Obs.jl")
include("Rinex3Obs.jl")

include("RinexObs.jl")
include("RinexNav.jl")

include("RinexNavCalculation.jl")


text:
module CalibrationAnalysis

using Reexport

@reexport using CalibrationErrors
@reexport using CalibrationTests
@reexport using ReliabilityDiagrams

end


text:
"""
    n_blocks(x,bs)

How many blocks of size `bs` are needed to process `x` items? Shortcut for
`div` with `RoundUp`.
"""
n_blocks(x, bs) = div(x, bs, RoundUp)

"""
    hasbits(flags, required)

Check that all `required` bits are present in `flags`.
"""
hasbits(x, y) = (x & y) == y


text:
include("../common/ortho_gram_schmidt.jl")
include("printMatrix.jl")

function test_main()
    srand(2222)
    Ngwx = 1000
    Nstates = 4
    psi = randn(Ngwx,Nstates) + im*randn(Ngwx,Nstates)
    psi_ortho1 = ortho_gram_schmidt(psi)

    U = inv(sqrtm(psi'*psi))
    printMatrix(U)

    U2 = (psi*U)'*(psi*U)
    printMatrix(U2)

    U1 = psi_ortho1'*psi_ortho1
    printMatrix(U1)

    Haux = zeros(ComplexF64,Nstates,Nstates)

end

test_main()


text:
const XC_LDA_X = 1   # Exchange
const XC_LDA_C_WIGNER = 2  # Wigner parametrization
const XC_LDA_C_RPA = 3  # Random Phase Approximation
const XC_LDA_C_HL = 4  # Hedin & Lundqvist
const XC_LDA_C_GL = 5  # Gunnarson & Lundqvist
const XC_LDA_C_XALPHA = 6  # Slater Xalpha
const XC_LDA_C_VWN = 7  # Vosko, Wilk, & Nusair (5)
const XC_LDA_C_VWN_RPA = 8  # Vosko, Wilk, & Nusair (RPA)
const XC_LDA_C_PZ =  9  # Perdew & Zunger
const XC_LDA_C_PZ_MOD = 10  # Perdew & Zunger (Modified)
const XC_LDA_C_OB_PZ = 11  # Ortiz & Ballone (PZ)
const XC_LDA_C_PW = 12  # Perdew & Wang
const XC_LDA_C_PW_MOD                 13  # Perdew & Wang (Modified)
const XC_LDA_C_OB_PW                  14  # Ortiz & Ballone (PW)
const XC_LDA_C_2D_AMGB                15  # Attaccalite et al
const XC_LDA_C_2D_PRM                 16  # Pittalis, Rasanen & Marques correlation in 2D
const XC_LDA_C_vBH                    17  # von Barth & Hedin
const XC_LDA_C_1D_CSC                 18  # Casula, Sorella, and Senatore 1D correlation
const XC_LDA_X_2D                     19  # Exchange in 2D
const XC_LDA_XC_TETER93               20  # Teter 93 parametrization
const XC_LDA_X_1D                     21  # Exchange in 1D
const XC_LDA_C_ML1                    22  # Modified LSD (version 1) of Proynov and Salahub
const XC_LDA_C_ML2                    23  # Modified LSD (version 2) of Proynov and Salahub
const XC_LDA_C_GOMBAS                 24  # Gombas parametrization
const XC_LDA_C_PW_RPA                 25  # Perdew & Wang fit of the RPA
const XC_LDA_C_1D_LOOS                26  # P-F Loos correlation LDA
const XC_LDA_C_RC04                   27  # Ragot-Cortona
const XC_LDA_C_VWN_1                  28  # Vosko, Wilk, & Nusair (1)
const XC_LDA_C_VWN_2                  29  # Vosko, Wilk, & Nusair (2)
const XC_LDA_C_VWN_3                  30  # Vosko, Wilk, & Nusair (3)
const XC_LDA_C_VWN_4                  31  # Vosko, Wilk, & Nusair (4)
const XC_LDA_XC_ZLP                   43  # Zhao, Levy & Parr, Eq. (20)
const XC_LDA_K_TF                     50  # Thomas-Fermi kinetic energy functional
const XC_LDA_K_LP                     51  # Lee and Parr Gaussian ansatz
const XC_LDA_XC_KSDT                 259  # Karasiev et al. parametrization
const XC_GGA_X_GAM                    32  # GAM functional from Minnesota
const XC_GGA_C_GAM                    33  # GAM functional from Minnesota
const XC_GGA_X_HCTH_A                 34  # HCTH-A
const XC_GGA_X_EV93                   35  # Engel and Vosko
const XC_GGA_X_BGCP                   38  # Burke, Cancio, Gould, and Pittalis
const XC_GGA_C_BGCP                   39  # Burke, Cancio, Gould, and Pittalis
const XC_GGA_X_LAMBDA_OC2_N           40  # lambda_OC2(N) version of PBE
const XC_GGA_X_B86_R  = 41  # Revised Becke 86 Xalpha,beta,gamma (with mod. grad. correction)
const XC_GGA_X_LAMBDA_CH_N            44  # lambda_CH(N) version of PBE
const XC_GGA_X_LAMBDA_LO_N            45  # lambda_LO(N) version of PBE
const XC_GGA_X_HJS_B88_V2             46  # HJS screened exchange corrected B88 version
const XC_GGA_C_Q2D                    47  # Chiodo et al
const XC_GGA_X_Q2D                    48  # Chiodo et al
const XC_GGA_X_PBE_MOL                49  # Del Campo, Gazquez, Trickey and Vela (PBE-like)
const XC_GGA_K_TFVW                   52  # Thomas-Fermi plus von Weiszaecker correction
const XC_GGA_K_REVAPBEINT             53  # interpolated version of REVAPBE
const XC_GGA_K_APBEINT                54  # interpolated version of APBE
const XC_GGA_K_REVAPBE                55  # revised APBE
const XC_GGA_X_AK13                   56  # Armiento & Kuemmel 2013
const XC_GGA_K_MEYER                  57  # Meyer,  Wang, and Young
const XC_GGA_X_LV_RPW86               58  # Berland and Hyldgaard
const XC_GGA_X_PBE_TCA                59  # PBE revised by Tognetti et al
const XC_GGA_X_PBEINT                 60  # PBE for hybrid interfaces
const XC_GGA_C_ZPBEINT                61  # spin-dependent gradient correction to PBEint
const XC_GGA_C_PBEINT                 62  # PBE for hybrid interfaces
const XC_GGA_C_ZPBESOL                63  # spin-dependent gradient correction to PBEsol
const XC_GGA_XC_OPBE_D                65  # oPBE_D functional of Goerigk and Grimme
const XC_GGA_XC_OPWLYP_D              66  # oPWLYP-D functional of Goerigk and Grimme
const XC_GGA_XC_OBLYP_D               67  # oBLYP-D functional of Goerigk and Grimme
const XC_GGA_X_VMT84_GE               68  # VMT{8,4} with constraint satisfaction with mu = mu_GE
const XC_GGA_X_VMT84_PBE              69  # VMT{8,4} with constraint satisfaction with mu = mu_PBE
const XC_GGA_X_VMT_GE                 70  # Vela, Medel, and Trickey with mu = mu_GE
const XC_GGA_X_VMT_PBE                71  # Vela, Medel, and Trickey with mu = mu_PBE
const XC_GGA_C_N12_SX                 79  # N12-SX functional from Minnesota
const XC_GGA_C_N12                    80  # N12 functional from Minnesota
const XC_GGA_X_N12                    82  # N12 functional from Minnesota
const XC_GGA_C_REGTPSS                83  # Regularized TPSS correlation (ex-VPBE)
const XC_GGA_C_OP_XALPHA              84  # one-parameter progressive functional (XALPHA version)
const XC_GGA_C_OP_G96                 85  # one-parameter progressive functional (G96 version)
const XC_GGA_C_OP_PBE                 86  # one-parameter progressive functional (PBE version)
const XC_GGA_C_OP_B88                 87  # one-parameter progressive functional (B88 version)
const XC_GGA_C_FT97                   88  # Filatov & Thiel correlation
const XC_GGA_C_SPBE                   89  # PBE correlation to be used with the SSB exchange
const XC_GGA_X_SSB_SW                 90  # Swarta, Sola and Bickelhaupt correction to PBE
const XC_GGA_X_SSB                    91  # Swarta, Sola and Bickelhaupt
const XC_GGA_X_SSB_D                  92  # Swarta, Sola and Bickelhaupt dispersion
const XC_GGA_XC_HCTH_407P             93  # HCTH/407+
const XC_GGA_XC_HCTH_P76              94  # HCTH p=7/6
const XC_GGA_XC_HCTH_P14              95  # HCTH p=1/4
const XC_GGA_XC_B97_GGA1              96  # Becke 97 GGA-1
const XC_GGA_C_HCTH_A                 97  # HCTH-A
const XC_GGA_X_BPCCAC                 98  # BPCCAC (GRAC for the energy)
const XC_GGA_C_REVTCA                 99  # Tognetti, Cortona, Adamo (revised)
const XC_GGA_C_TCA                   100  # Tognetti, Cortona, Adamo
const XC_GGA_X_PBE                   101  # Perdew, Burke & Ernzerhof exchange
const XC_GGA_X_PBE_R                 102  # Perdew, Burke & Ernzerhof exchange (revised)
const XC_GGA_X_B86                   103  # Becke 86 Xalpha,beta,gamma
const XC_GGA_X_HERMAN                104  # Herman et al original GGA
const XC_GGA_X_B86_MGC               105  # Becke 86 Xalpha,beta,gamma (with mod. grad. correction)
const XC_GGA_X_B88                   106  # Becke 88
const XC_GGA_X_G96                   107  # Gill 96
const XC_GGA_X_PW86                  108  # Perdew & Wang 86
const XC_GGA_X_PW91                  109  # Perdew & Wang 91
const XC_GGA_X_OPTX                  110  # Handy & Cohen OPTX 01
const XC_GGA_X_DK87_R1               111  # dePristo & Kress 87 (version R1)
const XC_GGA_X_DK87_R2               112  # dePristo & Kress 87 (version R2)
const XC_GGA_X_LG93                  113  # Lacks & Gordon 93
const XC_GGA_X_FT97_A                114  # Filatov & Thiel 97 (version A)
const XC_GGA_X_FT97_B                115  # Filatov & Thiel 97 (version B)
const XC_GGA_X_PBE_SOL               116  # Perdew, Burke & Ernzerhof exchange (solids)
const XC_GGA_X_RPBE                  117  # Hammer, Hansen & Norskov (PBE-like)
const XC_GGA_X_WC                    118  # Wu & Cohen
const XC_GGA_X_MPW91                 119  # Modified form of PW91 by Adamo & Barone
const XC_GGA_X_AM05                  120  # Armiento & Mattsson 05 exchange
const XC_GGA_X_PBEA                  121  # Madsen (PBE-like)
const XC_GGA_X_MPBE                  122  # Adamo & Barone modification to PBE
const XC_GGA_X_XPBE                  123  # xPBE reparametrization by Xu & Goddard
const XC_GGA_X_2D_B86_MGC            124  # Becke 86 MGC for 2D systems
const XC_GGA_X_BAYESIAN              125  # Bayesian best fit for the enhancement factor
const XC_GGA_X_PBE_JSJR              126  # JSJR reparametrization by Pedroza, Silva & Capelle
const XC_GGA_X_2D_B88                127  # Becke 88 in 2D
const XC_GGA_X_2D_B86                128  # Becke 86 Xalpha,beta,gamma
const XC_GGA_X_2D_PBE                129  # Perdew, Burke & Ernzerhof exchange in 2D
const XC_GGA_C_PBE                   130  # Perdew, Burke & Ernzerhof correlation
const XC_GGA_C_LYP                   131  # Lee, Yang & Parr
const XC_GGA_C_P86                   132  # Perdew 86
const XC_GGA_C_PBE_SOL               133  # Perdew, Burke & Ernzerhof correlation SOL
const XC_GGA_C_PW91                  134  # Perdew & Wang 91
const XC_GGA_C_AM05                  135  # Armiento & Mattsson 05 correlation
const XC_GGA_C_XPBE                  136  # xPBE reparametrization by Xu & Goddard
const XC_GGA_C_LM                    137  # Langreth and Mehl correlation
const XC_GGA_C_PBE_JRGX              138  # JRGX reparametrization by Pedroza, Silva & Capelle
const XC_GGA_X_OPTB88_VDW            139  # Becke 88 reoptimized to be used with vdW functional of const Dion et al
const XC_GGA_X_PBEK1_VDW             140  # PBE reparametrization for vdW

const XC_GGA_X_OPTPBE_VDW            141  # PBE reparametrization for vdW
const XC_GGA_X_RGE2                  142  # Regularized PBE
const XC_GGA_C_RGE2                  143  # Regularized PBE
const XC_GGA_X_RPW86                 144  # refitted Perdew & Wang 86
const XC_GGA_X_KT1                   145  # Keal and Tozer version 1
const XC_GGA_XC_KT2                  146  # Keal and Tozer version 2
const XC_GGA_C_WL                    147  # Wilson & Levy
const XC_GGA_C_WI                    148  # Wilson & Ivanov
const XC_GGA_X_MB88                  149  # Modified Becke 88 for proton transfer
const XC_GGA_X_SOGGA                 150  # Second-order generalized gradient approximation
const XC_GGA_X_SOGGA11               151  # Second-order generalized gradient approximation 2011
const XC_GGA_C_SOGGA11               152  # Second-order generalized gradient approximation 2011
const XC_GGA_C_WI0                   153  # Wilson & Ivanov initial version
const XC_GGA_XC_TH1                  154  # Tozer and Handy v. 1
const XC_GGA_XC_TH2                  155  # Tozer and Handy v. 2
const XC_GGA_XC_TH3                  156  # Tozer and Handy v. 3
const XC_GGA_XC_TH4                  157  # Tozer and Handy v. 4
const XC_GGA_X_C09X                  158  # C09x to be used with the VdW of Rutgers-Chalmers
const XC_GGA_C_SOGGA11_X             159  # To be used with HYB_GGA_X_SOGGA11_X
const XC_GGA_X_LB                    160  # van Leeuwen & Baerends
const XC_GGA_XC_HCTH_93              161  # HCTH functional fitted to  93 molecules
const XC_GGA_XC_HCTH_120             162  # HCTH functional fitted to 120 molecules
const XC_GGA_XC_HCTH_147             163  # HCTH functional fitted to 147 molecules
const XC_GGA_XC_HCTH_407             164  # HCTH functional fitted to 407 molecules
const XC_GGA_XC_EDF1                 165  # Empirical functionals from Adamson, Gill, and Pople
const XC_GGA_XC_XLYP                 166  # XLYP functional
const XC_GGA_XC_B97_D                170  # Grimme functional to be used with C6 vdW term
const XC_GGA_XC_PBE1W                173  # Functionals fitted for water
const XC_GGA_XC_MPWLYP1W             174  # Functionals fitted for water
const XC_GGA_XC_PBELYP1W             175  # Functionals fitted for water
const XC_GGA_X_LBM                   182  # van Leeuwen & Baerends modified
const XC_GGA_X_OL2                   183  # Exchange form based on Ou-Yang and Levy v.2
const XC_GGA_X_APBE                  184  # mu fixed from the semiclassical neutral atom
const XC_GGA_K_APBE                  185  # mu fixed from the semiclassical neutral atom
const XC_GGA_C_APBE                  186  # mu fixed from the semiclassical neutral atom
const XC_GGA_K_TW1                   187  # Tran and Wesolowski set 1 (Table II)
const XC_GGA_K_TW2                   188  # Tran and Wesolowski set 2 (Table II)
const XC_GGA_K_TW3                   189  # Tran and Wesolowski set 3 (Table II)
const XC_GGA_K_TW4                   190  # Tran and Wesolowski set 4 (Table II)
const XC_GGA_X_HTBS                  191  # Haas, Tran, Blaha, and Schwarz
const XC_GGA_X_AIRY                  192  # Constantin et al based on the Airy gas
const XC_GGA_X_LAG                   193  # Local Airy Gas
const XC_GGA_XC_MOHLYP               194  # Functional for organometallic chemistry
const XC_GGA_XC_MOHLYP2              195  # Functional for barrier heights
const XC_GGA_XC_TH_FL                196  # Tozer and Handy v. FL
const XC_GGA_XC_TH_FC                197  # Tozer and Handy v. FC
const XC_GGA_XC_TH_FCFO              198  # Tozer and Handy v. FCFO
const XC_GGA_XC_TH_FCO               199  # Tozer and Handy v. FCO
const XC_GGA_C_OPTC                  200  # Optimized correlation functional of Cohen and Handy
const XC_GGA_C_PBELOC                246  # Semilocal dynamical correlation
const XC_GGA_XC_VV10                 255  # Vydrov and Van Voorhis
const XC_GGA_C_PBEFE                 258  # PBE for formation energies
const XC_GGA_C_OP_PW91               262  # one-parameter progressive functional (PW91 version)
const XC_GGA_X_PBEFE                 265  # PBE for formation energies
const XC_GGA_X_CAP                   270  # Correct Asymptotic Potential
const XC_GGA_K_VW                    500  # von Weiszaecker functional
const XC_GGA_K_GE2                   501  # Second-order gradient expansion (l = 1/9)
const XC_GGA_K_GOLDEN                502  # TF-lambda-vW form by Golden (l = 13/45)
const XC_GGA_K_YT65                  503  # TF-lambda-vW form by Yonei and Tomishima (l = 1/5)
const XC_GGA_K_BALTIN                504  # TF-lambda-vW form by Baltin (l = 5/9)
const XC_GGA_K_LIEB                  505  # TF-lambda-vW form by Lieb (l = 0.185909191)
const XC_GGA_K_ABSP1                 506  # gamma-TFvW form by Acharya et al [g = 1 - 1.412/N^(1/3)]
const XC_GGA_K_ABSP2                 507  # gamma-TFvW form by Acharya et al [g = 1 - 1.332/N^(1/3)]
const XC_GGA_K_GR                    508  # gamma-TFvW form by Gazquez and Robles
const XC_GGA_K_LUDENA                509  # gamma-TFvW form by Ludena
const XC_GGA_K_GP85                  510  # gamma-TFvW form by Ghosh and Parr
const XC_GGA_K_PEARSON               511  # Pearson
const XC_GGA_K_OL1                   512  # Ou-Yang and Levy v.1
const XC_GGA_K_OL2                   513  # Ou-Yang and Levy v.2
const XC_GGA_K_FR_B88                514  # Fuentealba & Reyes (B88 version)
const XC_GGA_K_FR_PW86               515  # Fuentealba & Reyes (PW86 version)
const XC_GGA_K_DK                    516  # DePristo and Kress
const XC_GGA_K_PERDEW                517  # Perdew
const XC_GGA_K_VSK                   518  # Vitos, Skriver, and Kollar
const XC_GGA_K_VJKS                  519  # Vitos, Johansson, Kollar, and Skriver
const XC_GGA_K_ERNZERHOF             520  # Ernzerhof
const XC_GGA_K_LC94                  521  # Lembarki & Chermette
const XC_GGA_K_LLP                   522  # Lee, Lee & Parr
const XC_GGA_K_THAKKAR               523  # Thakkar 1992
const XC_GGA_X_WPBEH                 524  # short-range version of the PBE
const XC_GGA_X_HJS_PBE               525  # HJS screened exchange PBE version
const XC_GGA_X_HJS_PBE_SOL           526  # HJS screened exchange PBE_SOL version
const XC_GGA_X_HJS_B88               527  # HJS screened exchange B88 version
const XC_GGA_X_HJS_B97X              528  # HJS screened exchange B97x version
const XC_GGA_X_ITYH                  529  # short-range recipe for exchange GGA functionals
const XC_GGA_X_SFAT                  530  # short-range recipe for exchange GGA functionals
const XC_HYB_GGA_X_N12_SX             81  # N12-SX functional from Minnesota
const XC_HYB_GGA_XC_B97_1p           266  # version of B97 by Cohen and Handy
const XC_HYB_GGA_XC_B3PW91           401  # The original (ACM) hybrid of Becke
const XC_HYB_GGA_XC_B3LYP            402  # The (in)famous B3LYP
const XC_HYB_GGA_XC_B3P86            403  # Perdew 86 hybrid similar to B3PW91
const XC_HYB_GGA_XC_O3LYP            404  # hybrid using the optx functional
const XC_HYB_GGA_XC_mPW1K            405  # mixture of mPW91 and PW91 optimized for kinetics
const XC_HYB_GGA_XC_PBEH             406  # aka PBE0 or PBE1PBE
const XC_HYB_GGA_XC_B97              407  # Becke 97
const XC_HYB_GGA_XC_B97_1            408  # Becke 97-1
const XC_HYB_GGA_XC_B97_2            410  # Becke 97-2
const XC_HYB_GGA_XC_X3LYP            411  # hybrid by Xu and Goddard
const XC_HYB_GGA_XC_B1WC             412  # Becke 1-parameter mixture of WC and PBE
const XC_HYB_GGA_XC_B97_K            413  # Boese-Martin for Kinetics
const XC_HYB_GGA_XC_B97_3            414  # Becke 97-3
const XC_HYB_GGA_XC_MPW3PW           415  # mixture with the mPW functional
const XC_HYB_GGA_XC_B1LYP            416  # Becke 1-parameter mixture of B88 and LYP
const XC_HYB_GGA_XC_B1PW91           417  # Becke 1-parameter mixture of B88 and PW91
const XC_HYB_GGA_XC_mPW1PW           418  # Becke 1-parameter mixture of mPW91 and PW91
const XC_HYB_GGA_XC_MPW3LYP          419  # mixture of mPW and LYP
const XC_HYB_GGA_XC_SB98_1a          420  # Schmider-Becke 98 parameterization 1a
const XC_HYB_GGA_XC_SB98_1b          421  # Schmider-Becke 98 parameterization 1b
const XC_HYB_GGA_XC_SB98_1c          422  # Schmider-Becke 98 parameterization 1c
const XC_HYB_GGA_XC_SB98_2a          423  # Schmider-Becke 98 parameterization 2a
const XC_HYB_GGA_XC_SB98_2b          424  # Schmider-Becke 98 parameterization 2b
const XC_HYB_GGA_XC_SB98_2c          425  # Schmider-Becke 98 parameterization 2c
const XC_HYB_GGA_X_SOGGA11_X         426  # Hybrid based on SOGGA11 form
const XC_HYB_GGA_XC_HSE03            427  # the 2003 version of the screened hybrid HSE
const XC_HYB_GGA_XC_HSE06            428  # the 2006 version of the screened hybrid HSE
const XC_HYB_GGA_XC_HJS_PBE          429  # HJS hybrid screened exchange PBE version
const XC_HYB_GGA_XC_HJS_PBE_SOL      430  # HJS hybrid screened exchange PBE_SOL version
const XC_HYB_GGA_XC_HJS_B88          431  # HJS hybrid screened exchange B88 version
const XC_HYB_GGA_XC_HJS_B97X         432  # HJS hybrid screened exchange B97x version
const XC_HYB_GGA_XC_CAM_B3LYP        433  # CAM version of B3LYP
const XC_HYB_GGA_XC_TUNED_CAM_B3LYP  434  # CAM version of B3LYP tuned for excitations
const XC_HYB_GGA_XC_BHANDH           435  # Becke half-and-half
const XC_HYB_GGA_XC_BHANDHLYP        436  # Becke half-and-half with B88 exchange
const XC_HYB_GGA_XC_MB3LYP_RC04      437  # B3LYP with RC04 LDA
const XC_HYB_GGA_XC_MPWLYP1M         453  # MPW with 1 par. for metals/LYP
const XC_HYB_GGA_XC_REVB3LYP         454  # Revised B3LYP
const XC_HYB_GGA_XC_CAMY_BLYP        455  # BLYP with yukawa screening
const XC_HYB_GGA_XC_PBE0_13          456  # PBE0-1/3
const XC_HYB_GGA_XC_B3LYPs           459  # B3LYP* functional
const XC_HYB_GGA_XC_WB97             463  # Chai and Head-Gordon
const XC_HYB_GGA_XC_WB97X            464  # Chai and Head-Gordon
const XC_HYB_GGA_XC_LRC_WPBEH        465  # Long-range corrected functional by Rorhdanz et al
const XC_HYB_GGA_XC_WB97X_V          466  # Mardirossian and Head-Gordon
const XC_HYB_GGA_XC_LCY_PBE          467  # PBE with yukawa screening
const XC_HYB_GGA_XC_LCY_BLYP         468  # BLYP with yukawa screening
const XC_HYB_GGA_XC_LC_VV10          469  # Vydrov and Van Voorhis
const XC_HYB_GGA_XC_CAMY_B3LYP       470  # B3LYP with Yukawa screening
const XC_HYB_GGA_XC_WB97X_D          471  # Chai and Head-Gordon
const XC_HYB_GGA_XC_HPBEINT          472  # hPBEint
const XC_HYB_GGA_XC_LRC_WPBE         473  # Long-range corrected functional by Rorhdanz et al
const XC_HYB_GGA_XC_B3LYP5           475  # B3LYP with VWN functional 5 instead of RPA
const XC_HYB_GGA_XC_EDF2             476  # Empirical functional from Lin, George and Gill
const XC_HYB_GGA_XC_CAP0             477  # Correct Asymptotic Potential hybrid
const XC_MGGA_C_DLDF                  37  # Dispersionless Density Functional
const XC_MGGA_XC_ZLP                  42  # Zhao, Levy & Parr, Eq. (21)
const XC_MGGA_XC_OTPSS_D              64  # oTPSS_D functional of Goerigk and Grimme
const XC_MGGA_C_CS                    72  # Colle and Salvetti
const XC_MGGA_C_MN12_SX               73  # Worker for MN12-SX functional
const XC_MGGA_C_MN12_L                74  # MN12-L functional from Minnesota
const XC_MGGA_C_M11_L                 75  # M11-L functional from Minnesota
const XC_MGGA_C_M11                   76  # Worker for M11 functional
const XC_MGGA_C_M08_SO                77  # Worker for M08-SO functional
const XC_MGGA_C_M08_HX                78  # Worker for M08-HX functional
const XC_MGGA_X_LTA                  201  # Local tau approximation of Ernzerhof & Scuseria
const XC_MGGA_X_TPSS                 202  # Perdew, Tao, Staroverov & Scuseria exchange
const XC_MGGA_X_M06_L                203  # M06-Local functional of Minnesota
const XC_MGGA_X_GVT4                 204  # GVT4 from Van Voorhis and Scuseria
const XC_MGGA_X_TAU_HCTH             205  # tau-HCTH from Boese and Handy
const XC_MGGA_X_BR89                 206  # Becke-Roussel 89
const XC_MGGA_X_BJ06                 207  # Becke & Johnson correction to Becke-Roussel 89
const XC_MGGA_X_TB09                 208  # Tran & Blaha correction to Becke & Johnson
const XC_MGGA_X_RPP09 = 209  # Rasanen, Pittalis, and Proetto correction to Becke & Johnson

const XC_MGGA_X_2D_PRHG07            210  # Pittalis, Rasanen, Helbig, Gross Exchange Functional
const XC_MGGA_X_2D_PRHG07_PRP10      211  # PRGH07 with PRP10 correction
const XC_MGGA_X_REVTPSS              212  # revised Perdew, Tao, Staroverov & Scuseria exchange
const XC_MGGA_X_PKZB                 213  # Perdew, Kurth, Zupan, and Blaha
const XC_MGGA_X_M05                  214  # Worker for M05 functional
const XC_MGGA_X_M05_2X               215  # Worker for M05-2X functional
const XC_MGGA_X_M06_HF               216  # Worker for M06-HF functional
const XC_MGGA_X_M06                  217  # Worker for M06 functional
const XC_MGGA_X_M06_2X               218  # Worker for M06-2X functional
const XC_MGGA_X_M08_HX               219  # Worker for M08-HX functional
const XC_MGGA_X_M08_SO               220  # Worker for M08-SO functional
const XC_MGGA_X_MS0                  221  # MS exchange of Sun, Xiao, and Ruzsinszky
const XC_MGGA_X_MS1                  222  # MS1 exchange of Sun, et al
const XC_MGGA_X_MS2                  223  # MS2 exchange of Sun, et al
const XC_MGGA_X_M11                  225  # Worker for M11 functional
const XC_MGGA_X_M11_L                226  # M11-L functional from Minnesota
const XC_MGGA_X_MN12_L               227  # MN12-L functional from Minnesota
const XC_MGGA_C_CC06                 229  # Cancio and Chou 2006
const XC_MGGA_X_MK00                 230  # Exchange for accurate virtual orbital energies
const XC_MGGA_C_TPSS                 231  # Perdew, Tao, Staroverov & Scuseria correlation
const XC_MGGA_C_VSXC                 232  # VSxc from Van Voorhis and Scuseria (correlation part)
const XC_MGGA_C_M06_L                233  # M06-Local functional from Minnesota
const XC_MGGA_C_M06_HF               234  # Worker for M06-HF functional
const XC_MGGA_C_M06                  235  # Worker for M06 functional
const XC_MGGA_C_M06_2X               236  # Worker for M06-2X functional
const XC_MGGA_C_M05                  237  # Worker for M05 functional
const XC_MGGA_C_M05_2X               238  # Worker for M05-2X functional
const XC_MGGA_C_PKZB                 239  # Perdew, Kurth, Zupan, and Blaha
const XC_MGGA_C_BC95                 240  # Becke correlation 95
const XC_MGGA_C_REVTPSS              241  # revised TPSS correlation
const XC_MGGA_XC_TPSSLYP1W           242  # Functionals fitted for water
const XC_MGGA_X_MK00B                243  # Exchange for accurate virtual orbital energies (v. B)
const XC_MGGA_X_BLOC                 244  # functional with balanced localization
const XC_MGGA_X_MODTPSS              245  # Modified Perdew, Tao, Staroverov & Scuseria exchange
const XC_MGGA_C_TPSSLOC              247  # Semilocal dynamical correlation
const XC_MGGA_X_MBEEF                249  # mBEEF exchange
const XC_MGGA_X_MBEEFVDW             250  # mBEEF-vdW exchange
const XC_MGGA_XC_B97M_V              254  # Mardirossian and Head-Gordon
const XC_MGGA_X_MVS                  257  # MVS exchange of Sun, Perdew, and Ruzsinszky
const XC_MGGA_X_MN15_L               260  # MN15-L functional from Minnesota
const XC_MGGA_C_MN15_L               261  # MN15-L functional from Minnesota
const XC_MGGA_X_SCAN                 263  # SCAN exchange of Sun, Ruzsinszky, and Perdew
const XC_MGGA_C_SCAN                 267  # SCAN correlation
const XC_MGGA_C_MN15                 269  # MN15 functional from Minnesota
const XC_HYB_MGGA_X_DLDF              36  # Dispersionless Density Functional
const XC_HYB_MGGA_X_MS2H             224  # MS2 hybrid exchange of Sun, et al
const XC_HYB_MGGA_X_MN12_SX          248  # MN12-SX hybrid functional from Minnesota
const XC_HYB_MGGA_X_SCAN0            264  # SCAN hybrid
const XC_HYB_MGGA_X_MN15             268  # MN15 functional from Minnesota
const XC_HYB_MGGA_XC_M05             438  # M05 functional from Minnesota
const XC_HYB_MGGA_XC_M05_2X          439  # M05-2X functional from Minnesota
const XC_HYB_MGGA_XC_B88B95          440  # Mixture of B88 with BC95 (B1B95)
const XC_HYB_MGGA_XC_B86B95          441  # Mixture of B86 with BC95
const XC_HYB_MGGA_XC_PW86B95         442  # Mixture of PW86 with BC95
const XC_HYB_MGGA_XC_BB1K            443  # Mixture of B88 with BC95 from Zhao and Truhlar
const XC_HYB_MGGA_XC_M06_HF          444  # M06-HF functional from Minnesota
const XC_HYB_MGGA_XC_MPW1B95         445  # Mixture of mPW91 with BC95 from Zhao and Truhlar
const XC_HYB_MGGA_XC_MPWB1K          446  # Mixture of mPW91 with BC95 for kinetics
const XC_HYB_MGGA_XC_X1B95           447  # Mixture of X with BC95
const XC_HYB_MGGA_XC_XB1K            448  # Mixture of X with BC95 for kinetics
const XC_HYB_MGGA_XC_M06             449  # M06 functional from Minnesota

const XC_HYB_MGGA_XC_M06_2X          450  # M06-2X functional from Minnesota
const XC_HYB_MGGA_XC_PW6B95          451  # Mixture of PW91 with BC95 from Zhao and Truhlar
const XC_HYB_MGGA_XC_PWB6K           452  # Mixture of PW91 with BC95 from Zhao and Truhlar for kinetics
const XC_HYB_MGGA_XC_TPSSH           457  # TPSS hybrid
const XC_HYB_MGGA_XC_REVTPSSH        458  # revTPSS hybrid
const XC_HYB_MGGA_XC_M08_HX          460  # M08-HX functional from Minnesota
const XC_HYB_MGGA_XC_M08_SO          461  # M08-SO functional from Minnesota
const XC_HYB_MGGA_XC_M11             462  # M11    functional from Minnesota
const XC_HYB_MGGA_X_MVSH             474  # MVS hybrid
const XC_HYB_MGGA_XC_WB97M_V         531  # Mardirossian and Head-Gordon


text:
using Distributed, Distributions, Serialization

test_e_pth = "/bench/PhD/NGS_binaries/BMI/test_e"

@info "Assembling worker pool..."

#DISTRIBUTED CLUSTERs CONSTANTS
remote_machine = "10.0.0.3"
no_local_processes = 1

@info "Spawning local cluster workers..."
worker_pool=addprocs(no_local_processes, topology=:master_worker)

@info "Loading worker libraries everywhere..."
@everywhere using BioMotifInference,Random
@everywhere Random.seed!(myid()*100000)

#JOB CONSTANTS
models_to_permute=5000
func_limit=30
clamp=.02
funcvec=full_perm_funcvec
push!(funcvec, BioMotifInference.permute_source)
push!(funcvec, BioMotifInference.permute_source)
args=[Vector{Tuple{Symbol,Any}}() for i in 1:length(funcvec)]
args[end-1]=[(:weight_shift_freq,0.),(:length_change_freq,1.),(:length_perm_range,1:1)]
args[end]=[(:weight_shift_freq,.1),(:length_change_freq,0.),(:weight_shift_dist,Uniform(.00001,.01))]

instruct = Permute_Instruct(funcvec, ones(length(funcvec))./length(funcvec), models_to_permute, func_limit, clamp; args=args)

display_rotation=[true,10,1,[[:tuning_disp,:lh_disp,:src_disp],[:conv_plot,:liwi_disp,:ens_disp]]]

@info "Beginning nested sampling for test ensemble..."
test_e=deserialize(test_e_pth*"/ens")
test_logZ = converge_ensemble!(test_e, instruct, worker_pool, .001, backup=(true,1),tuning_disp=true,lh_disp=true,src_disp=true, disp_rotate_inst=display_rotation)

rmprocs(worker_pool)

@info "Job done!"


text:
using PkgBenchmark
results = benchmarkpkg("TaylorSeries")
show(results)

#=
# specify tag and uncommit to benchmark versus prior tagged version
tag =
results = judge("TaylorSeries", tag)
show(results)
=#

export_markdown("results.md", results)


text:
# # Test Registration Functions
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/example_placeholder.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/example_placeholder.ipynb)

# In this example page I will show how the registration algorithms ([`lap`](@ref), [`pflap`](@ref), [`sparse_lap`](@ref) and [`sparse_pflap`](@ref))
# perform. I will will show their speed and the quality of their outputs.

# First I some testing images:

## using the package and TimerOutputs
using LAP_julia, TimerOutputs

## default arguments:
img, imgw, flow = gen_init();

# See the [`Basic Interaction`](@ref basic_interaction) section or the [`Public Documentation`](@ref public_api) for more advanced input generation

# __These are the differences between the `target` (`img`) and `source` (`imgw`) images.__

imgoverlay(img, imgw, figtitle="Target vs Source")

# ## [`test_registration_alg`](@ref)

# This function times the registration algorithm on the inside and then test the output flow - `flow_est` and aligned source - `source_reg`.

## choose a method
method = sparse_pflap
## start a timer
timer = TimerOutput("ALG: " * string(method))
## set the keyword arguments of the method
method_kwargs = Dict(:display => false, :timer => timer, :match_source_histogram => false)
## set the arguments of the method if any
method_args = []
## start the test
flow_est, source_reg, timer, results = test_registration_alg(method,
                                                             img,
                                                             imgw,
                                                             flow,
                                                             timer=timer,
                                                             method_args=method_args,
                                                             only_flow_compare=false, # this adds source reg tests.
                                                             method_kwargs=method_kwargs);

# #### Check the generated flow and the original flow:
# __Ground Truth:__
showflow(flow.*(-1), figtitle="Ground Truth")
# __Estimated Flow:__
showflow(flow_est, figtitle="Estimated Flow")

# ### Compare it with the [`pflap`](@ref) function

method = pflap
## same inputs, new timer
timer = TimerOutput("ALG: " * string(method))
flow_est, source_reg, timer, results = test_registration_alg(method,
                                                             img,
                                                             imgw,
                                                             flow,
                                                             timer=timer,
                                                             method_args=method_args,
                                                             only_flow_compare=false, # this adds source reg tests.
                                                             method_kwargs=method_kwargs);

# #### Check the generated flow and the original flow:
# __Ground Truth:__
showflow(flow.*(-1), figtitle="Ground Truth")
# __Estimated Flow:__
showflow(flow_est, figtitle="Estimated Flow")

# The [`time_reg_alg`](@ref) is used to time the registration algorithm, when the original flow isn't known.

# ## using the `display` keyword
# The registration algorithms have a display keyword argument, that can be set to true to output figures
# showing the `flow_est` at each iteration and print some extra debug info.
method = sparse_pflap
timer = TimerOutput("ALG: " * string(method))
method_kwargs = Dict(:display => true, :timer => timer, :match_source_histogram => false)
method_args = []
flow_est, source_reg, timer, results, (figs,) = test_registration_alg(method,
                                                                      img,
                                                                      imgw,
                                                                      flow,
                                                                      timer=timer,
                                                                      method_args=method_args,
                                                                      only_flow_compare=false, # this adds source reg tests.
                                                                      method_kwargs=method_kwargs);


text:
using DynamicalSystemsBase
using DelayEmbeddings
using Test
import Peaks

println("\nTesting garcia_almeida.jl...")
@testset "Garcia & Almeida method" begin

lo = Systems.lorenz()

tr = trajectory(lo, 80; Δt = 0.01, Ttr = 10)

x = tr[:, 1]
Y = Dataset(x)
Y = standardize(Y)

@testset "N-statistic" begin
    τs = 0:100
    N , NN_distances = n_statistic(Y, x, w=0, T=1, τs = τs)

    T = 17
    N2 , NN_distances2 = n_statistic(Y, x, w=0, T=T, τs = τs)

    # check whether the `d_E1`-statistic is the same
    @test NN_distances[1][1:50] == NN_distances2[1][1:50]
    @test NN_distances[5][1:100] == NN_distances2[5][1:100]

    min_dist = 12
    (max_1_idx,_) = Peaks.findmaxima(N,min_dist)
    (max_2_idx,_) = Peaks.findmaxima(N2,min_dist)

    @test 30 ≤ max_1_idx[1] ≤ 35
    @test 30 ≤ max_2_idx[1] ≤ 35
    @test 69 ≤ max_1_idx[2] ≤ 75
    @test 69 ≤ max_2_idx[2] ≤ 75

    # # plot N-Statistic for the Lorenz system as in Fig. 2(a) in [^Garcia2005a]
    # using Plots
    # plot(N,linewidth = 2, label = "T=1 (as in the Paper)", xaxis=:log)
    # plot!(N2,linewidth = 2, label = "T=$T", xaxis=:log)
    # plot!(N,seriestype = :scatter, label = "",xaxis=:log)
    # plot!(N2,seriestype = :scatter, label = "",xaxis=:log)
    # plot!(title = "Lorenz System as in Fig. 2(a) in the Paper")
    # xlabel!("τ")
    # ylabel!("N-Statistic")
end

@testset "garcia_embed univariate" begin
    Y_act, τ_vals, ts_vals, FNNs, NS = garcia_almeida_embedding(x; τs=0:100,  w = 17, T = 17)
    Y_act2, τ_vals2, ts_vals2, FNNs2, NS2 = garcia_almeida_embedding(x; τs=0:100,  w = 1, T = 1)

    @test size(Y_act,2) == 2
    @test size(Y_act2,2) == 3
    @test 10 ≤ τ_vals[2] ≤ 18
    @test 1 ≤ τ_vals2[2] ≤ 3
    @test FNNs[end]<=0.05
    @test FNNs2[end]>FNNs2[end-1]

    # using Plots
    # plot3d(Y_act[:,1],Y_act[:,2],Y_act[:,3], marker=2, camera = (6, 4))
    # plot!(title = "Reconstruction of Lorenz attractor using Garcia&Almeida univariate embedding")
    #
    # plot3d(Y_act2[:,1],Y_act2[:,2],Y_act2[:,3], marker=2, camera = (6, 4))
    # plot!(title = "Reconstruction of Lorenz attractor using Garcia&Almeida univariate embedding as in the paper")
end

@testset "garcia_embed multivariate" begin
    Y_act, τ_vals, ts_vals, FNNs, NS = garcia_almeida_embedding(tr; τs=0:100,  w = 17, T = 17)
    Y_act2, τ_vals2, ts_vals2, FNNs2, NS2 = garcia_almeida_embedding(tr; τs=0:100,  w = 1, T = 1)

    @test size(Y_act,2) == size(Y_act2,2) == 3
    @test ts_vals[1] == ts_vals2[1] == ts_vals2[2] == ts_vals2[3] ==1
    @test ts_vals[2] == 2
    @test ts_vals[3] == 3
    @test 10 ≤ τ_vals[3] ≤ 18
    @test 1 ≤ τ_vals2[3] ≤ 3

    # try to reproduce Fig.2a in [^Garcia2005b]
    tra = Dataset(hcat(tr[:,1], tr[:,3]))
    tra = standardize(tra)
    taus = 0:100

    N , _ = n_statistic(tra, tra[:,1], w=17, T=17, τs = taus)
    N2 , _ = n_statistic(tra, tra[:,2], w=17, T=17, τs = taus)

    _, minis2 = findmin(N2)

    @test 20 ≤ taus[minis2] ≤ 30

    # using Plots
    # NN = hcat(N,N2)
    # plot(NN, xaxis=:log)
    # plot(NN)

end

end


text:
@testset "flatten/expand" begin

    @testset "fllatten" begin
        @test flatten(ones(3, 4, 4)) == ones(3, 16)

        @test_throws MethodError flatten(ones(3, 4))
    end

    @testset "expand" begin
        @test expand(ones(3, 16)) == ones(3, 4, 4)

        @test_throws ErrorException expand(ones(3, 15))
    end

    # Simple regression
    X = rand(10, 512, 512)
    @test expand(flatten(X)) == X
end

@testset "rotation" begin
    @testset "rotate" begin
        X = rand(10, 512, 512)

        # test no-op
        @test rotate(X, zeros(10)) === X
    end
    @testset "derotate" begin
        X = rand(10, 512, 512)

        # test no-op
        @test derotate(X, zeros(10)) === X
    end

    # Simple regression (but be careful to only test some inner part that won't be clipped)
    X = rand(10, 512, 512)
    θ = 90 .* rand(10)
    @test filter(!isnan, derotate(X, θ)) ≈ filter(!isnan, rotate(X, -θ))
    @test filter(!isnan, rotate(X, θ)) ≈ filter(!isnan, derotate(X, -θ))

    targ = deepcopy(X)
    derotate!(X, θ)
    @test X !== targ
    X = rand(10, 512, 512)
    rotate!(X, θ)
    @test X !== targ

end

@testset "combining" begin
    X = rand(10, 512, 512)
    angles = sort!(90 .* rand(10))

    @test combine(X) == combine(X, zeros(10))
    @test combine(X) == combine(X; method = median)
    @test combine(X; method = mean) == mean(X, dims = 1)[1, :, :]

    @test filter(!isnan, combine(X, angles)) ≈ filter(!isnan, combine(derotate(X, angles)))
    @test any(combine(X, angles) .=== NaN)

end


text:
using PublicSuffix
using Base.Test

const NULL = ""

function checkPublicSuffix(domain::String, expected_ps::String)
    ps = NULL
    try
        ps = public_suffix(domain)
    catch ex
        #Base.error_show(STDERR, ex, catch_backtrace())
    end
    #println("domain:[$domain]\t\tps:[$ps]\t\texpected:[$expected_ps]")
    @test ps == expected_ps
end

# generated from
# Any copyright is dedicated to the Public Domain.
# http://creativecommons.org/publicdomain/zero/1.0/

# NULL input.
#checkPublicSuffix(NULL, NULL)
# Mixed case.
checkPublicSuffix("COM", NULL)
checkPublicSuffix("example.COM", "example.com")
checkPublicSuffix("WwW.example.COM", "example.com")
# Leading dot.
checkPublicSuffix(".com", NULL)
checkPublicSuffix(".example", NULL)
checkPublicSuffix(".example.com", NULL)
checkPublicSuffix(".example.example", NULL)
# Unlisted TLD.
checkPublicSuffix("example", NULL)
checkPublicSuffix("example.example", "example.example")
checkPublicSuffix("b.example.example", "example.example")
checkPublicSuffix("a.b.example.example", "example.example")
# Listed, but non-Internet, TLD.
#checkPublicSuffix("local", NULL)
#checkPublicSuffix("example.local", NULL)
#checkPublicSuffix("b.example.local", NULL)
#checkPublicSuffix("a.b.example.local", NULL)
# TLD with only 1 rule.
checkPublicSuffix("biz", NULL)
checkPublicSuffix("domain.biz", "domain.biz")
checkPublicSuffix("b.domain.biz", "domain.biz")
checkPublicSuffix("a.b.domain.biz", "domain.biz")
# TLD with some 2-level rules.
checkPublicSuffix("com", NULL)
checkPublicSuffix("example.com", "example.com")
checkPublicSuffix("b.example.com", "example.com")
checkPublicSuffix("a.b.example.com", "example.com")
checkPublicSuffix("uk.com", NULL)
checkPublicSuffix("example.uk.com", "example.uk.com")
checkPublicSuffix("b.example.uk.com", "example.uk.com")
checkPublicSuffix("a.b.example.uk.com", "example.uk.com")
checkPublicSuffix("test.ac", "test.ac")
# TLD with only 1 (wildcard) rule.
checkPublicSuffix("cy", NULL)
checkPublicSuffix("c.cy", NULL)
checkPublicSuffix("b.c.cy", "b.c.cy")
checkPublicSuffix("a.b.c.cy", "b.c.cy")
# More complex TLD.
checkPublicSuffix("jp", NULL)
checkPublicSuffix("test.jp", "test.jp")
checkPublicSuffix("www.test.jp", "test.jp")
checkPublicSuffix("ac.jp", NULL)
checkPublicSuffix("test.ac.jp", "test.ac.jp")
checkPublicSuffix("www.test.ac.jp", "test.ac.jp")
checkPublicSuffix("kyoto.jp", NULL)
checkPublicSuffix("test.kyoto.jp", "test.kyoto.jp")
checkPublicSuffix("ide.kyoto.jp", NULL)
checkPublicSuffix("b.ide.kyoto.jp", "b.ide.kyoto.jp")
checkPublicSuffix("a.b.ide.kyoto.jp", "b.ide.kyoto.jp")
checkPublicSuffix("c.kobe.jp", NULL)
checkPublicSuffix("b.c.kobe.jp", "b.c.kobe.jp")
checkPublicSuffix("a.b.c.kobe.jp", "b.c.kobe.jp")
checkPublicSuffix("city.kobe.jp", "city.kobe.jp")
checkPublicSuffix("www.city.kobe.jp", "city.kobe.jp")
# TLD with a wildcard rule and exceptions.
checkPublicSuffix("ck", NULL)
checkPublicSuffix("test.ck", NULL)
checkPublicSuffix("b.test.ck", "b.test.ck")
checkPublicSuffix("a.b.test.ck", "b.test.ck")
checkPublicSuffix("www.ck", "www.ck")
checkPublicSuffix("www.www.ck", "www.ck")
# US K12.
checkPublicSuffix("us", NULL)
checkPublicSuffix("test.us", "test.us")
checkPublicSuffix("www.test.us", "test.us")
checkPublicSuffix("ak.us", NULL)
checkPublicSuffix("test.ak.us", "test.ak.us")
checkPublicSuffix("www.test.ak.us", "test.ak.us")
checkPublicSuffix("k12.ak.us", NULL)
checkPublicSuffix("test.k12.ak.us", "test.k12.ak.us")
checkPublicSuffix("www.test.k12.ak.us", "test.k12.ak.us")
# IDN labels.
checkPublicSuffix("食狮.com.cn", "食狮.com.cn")
checkPublicSuffix("食狮.公司.cn", "食狮.公司.cn")
checkPublicSuffix("www.食狮.公司.cn", "食狮.公司.cn")
checkPublicSuffix("shishi.公司.cn", "shishi.公司.cn")
checkPublicSuffix("公司.cn", NULL)
checkPublicSuffix("食狮.中国", "食狮.中国")
checkPublicSuffix("www.食狮.中国", "食狮.中国")
checkPublicSuffix("shishi.中国", "shishi.中国")
checkPublicSuffix("中国", NULL)
# Same as above, but punycoded.
# TODO: puny code support to be added
# checkPublicSuffix("xn--85x722f.com.cn", "xn--85x722f.com.cn")
# checkPublicSuffix("xn--85x722f.xn--55qx5d.cn", "xn--85x722f.xn--55qx5d.cn")
# checkPublicSuffix("www.xn--85x722f.xn--55qx5d.cn", "xn--85x722f.xn--55qx5d.cn")
# checkPublicSuffix("shishi.xn--55qx5d.cn", "shishi.xn--55qx5d.cn")
# checkPublicSuffix("xn--55qx5d.cn", NULL)
# checkPublicSuffix("xn--85x722f.xn--fiqs8s", "xn--85x722f.xn--fiqs8s")
# checkPublicSuffix("www.xn--85x722f.xn--fiqs8s", "xn--85x722f.xn--fiqs8s")
# checkPublicSuffix("shishi.xn--fiqs8s", "shishi.xn--fiqs8s")
# checkPublicSuffix("xn--fiqs8s", NULL)


text:
using Plots

function bracket_minimum(f, x=0; s=1e-2, k=2.0)
    a, ya = x, f(x)
    b, yb = a + s, f(a + s)
    if yb > ya
        a, b = b, a
        ya, yb = yb, ya
        s = -s
    end
    while true
        c, yc = b + s, f(b + s)
        if yc > yb
            return a < c ? (a,b, c) : (c, b, a)
        end
        a, ya, b, yb = b, yb, c, yc
        s *= k
    end
end

function quadratic_fit_search(f, a, b, c, n)
    ya, yb, yc = f(a), f(b), f(c)
    for i in 1:n-3
        x = 0.5*(ya*(b^2-c^2)+yb*(c^2-a^2)+yc*(a^2-b^2)) /
            (ya*(b-c) +yb*(c-a) +yc*(a-b))
        yx = f(x)
        if x > b
            if yx > yb
                c, yc = x, yx
            else
                a, ya, b, yb = b, yb, x, yx
            end
        elseif x < b
            if yx > yb
                a, ya = x, yx
            else
                c, yc, b, yb = b, yb, x, yx
            end
        end
    end
    return (a, b, c)
end

f = x->sin(10x)+cos(3x)
x0 = 3.4
#busca do intervalo
a,b,c = bracket_minimum(f,x0)

#numero de iteracoes
φ = MathConstants.φ
n=(c-a)/(10^-3*log(φ))

am,bm,cm = quadratic_fit_search(f,a,b,c,n)
(am,f(am))
am,bm,cm

plot(f,0.,4.,lw=2,draw_arrow="true")
scatter!([a,b,c],[f(a),f(b),f(c)])
scatter!([am,bm,cm],[f(am),f(bm),f(cm)])


text:
@doc """`care(A, B, Q, R)`

Compute 'X', the solution to the continuous-time algebraic Riccati equation,
defined as A'X + XA - (XB)R^-1(B'X) + Q = 0, where R is non-singular.

Algorithm taken from:
Laub, "A Schur Method for Solving Algebraic Riccati Equations."
http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf
""" ->
function care(A, B, Q, R)
    G = try
        B*inv(R)*B'
    catch
        error("R must be non-singular.")
    end

    Z = [A  -G;
        -Q  -A']

    S = schurfact(Z)
    S = ordschur(S, real(S.values).<0)
    U = S.Z

    (m, n) = size(U)
    U11 = U[1:div(m, 2), 1:div(n,2)]
    U21 = U[div(m,2)+1:m, 1:div(n,2)]
    return U21/U11
end

@doc """`dare(A, B, Q, R)`

Compute `X`, the solution to the discrete-time algebraic Riccati equation,
defined as A'XA - X - (A'XB)(B'XB + R)^-1(B'XA) + Q = 0, where A and R
are non-singular.

Algorithm taken from:
Laub, "A Schur Method for Solving Algebraic Riccati Equations."
http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf
""" ->
function dare(A, B, Q, R)
    G = try
        B*inv(R)*B'
    catch
        error("R must be non-singular.")
    end

    Ait = try
        inv(A)'
    catch
        error("A must be non-singular.")
    end

    Z = [A + G*Ait*Q   -G*Ait;
         -Ait*Q        Ait]

    S = schurfact(Z)
    S = ordschur(S, abs(S.values).<=1)
    U = S.Z

    (m, n) = size(U)
    U11 = U[1:div(m, 2), 1:div(n,2)]
    U21 = U[div(m,2)+1:m, 1:div(n,2)]
    return U21/U11
end

@doc """`dlyap(A, Q)`

Compute the solution "X" to the discrete Lyapunov equation
"AXA' - X + Q = 0".
""" ->
function dlyap(A, Q)
    lhs = kron(A, conj(A))
    lhs = eye(size(lhs, 1)) - lhs
    x = lhs\reshape(Q, prod(size(Q)), 1)
    return reshape(x, size(Q))
end

@doc """`gram(sys, opt)`

Compute the grammian of system `sys`. If `opt` is `:c`, computes the
controllability grammian. If `opt` is `:o`, computes the observability
grammian.""" ->
function gram(sys::StateSpace, opt::Symbol)
    if !isstable(sys)
        error("gram only valid for stable A")
    end
    func = iscontinuous(sys) ? lyap : dlyap
    if opt == :c
        return func(sys.A, sys.B*sys.B')
    elseif opt == :o
        return func(sys.A', sys.C'*sys.C)
    else
        error("opt must be either :c for controllability grammian, or :o for
                observability grammian")
    end
end

@doc """`obsv(A, C)` or `obsv(sys)`

Compute the observability matrix for the system described by `(A, C)` or `sys`.

Note that checking for observability by computing the rank from `obsv` is
not the most numerically accurate way, a better method is checking if
`gram(sys, :o)` is positive definite.""" ->
function obsv(A, C)
    n = size(A, 1)
    ny = size(C, 1)
    if n != size(C, 2)
        error("C must have the same number of columns as A")
    end
    res = zeros(n*ny, n)
    res[1:ny, :] = C
    for i=1:n-1
        res[(1 + i*ny):(1 + i)*ny, :] = res[((i - 1)*ny + 1):i*ny, :] * A
    end
    return res
end
obsv(sys::StateSpace) = obsv(sys.A, sys.C)

@doc """`ctrb(A, B)` or `ctrb(sys)`

Compute the controllability matrix for the system described by `(A, B)` or
`sys`.

Note that checking for controllability by computing the rank from
`obsv` is not the most numerically accurate way, a better method is
checking if `gram(sys, :c)` is positive definite.""" ->
function ctrb(A, B)
    n = size(A, 1)
    nu = size(B, 2)
    if n != size(B, 1)
        error("B must have the same number of rows as A")
    end
    res = zeros(n, n*nu)
    res[:, 1:nu] = B
    for i=1:n-1
        res[:, (1 + i*nu):(1 + i)*nu] = A * res[:, ((i - 1)*nu + 1):i*nu]
    end
    return res
end
ctrb(sys::StateSpace) = ctrb(sys.A, sys.B)

@doc """`P = covar(sys, W)`

Calculate the stationary covariance `P = E[y(t)y(t)']` of an lti-model `sys`, driven by gaussian
white noise 'w' of covariance `E[w(t)w(τ)]=W*δ(t-τ)` where δ is the dirac delta.

The ouput is if Inf if the system is unstable. Passing white noise directly to
the output will result in infinite covariance in the corresponding outputs
(D*W*D.' .!= 0) for contunuous systems.""" ->
function covar(sys::StateSpace, W::StridedMatrix)
    (A, B, C, D) = (sys.A, sys.B, sys.C, sys.D)
    if size(B,2) != size(W, 1) || size(W, 1) != size(W, 2)
        error("W must be a square matrix the same size as `sys.B` columns")
    end
    if !isstable(sys)
        return fill(Inf,(size(C,1),size(C,1)))
    end
    func = iscontinuous(sys) ? lyap : dlyap
    Q = try
        func(A, B*W*B')
    catch
        error("No solution to the Lyapunov equation was found in covar")
    end
    P = C*Q*C'
    if iscontinuous(sys)
        #Variance and covariance infinite for direct terms
        directNoise = D*W*D'
        for i in 1:size(C,1)
            if directNoise[i,i] != 0
                P[i,:] = Inf
                P[:,i] = Inf
            end
        end
    else
        P += D*W*D'
    end
    return P
end

covar(sys::TransferFunction, W::StridedMatrix) = covar(ss(sys), W)

# Note: the following function returns a Float64 when p=2, and a Tuple{Float64,Float64} when p=Inf.
# This varargout behavior might not be a good idea.
# The H∞ norm computation is probably not as accurate as with SLICOT, but this seems to be still
# reasonably ok as a first step
@doc """
`..  norm(sys[, p; tol=...])`

Compute the H`p`-norm of the LTI system `sys`. `p` can be either `2` or `Inf`
(default is 2). For the H∞ norm, the function returns the pair `(val,fpeak)`,
where `fpeak` (in rad/TimeUnit) is where the gain achieves its peak value `val`.

sys is first converted to a state space model if needed.

`tol` is an optional keyword argument, used currently only for the computation of H∞ norms.
It represents the desired relative accuracy for the computed H∞ norm
(default 1e-12 - not an absolute certificate however).

The H∞ norm computation implements the 'two-step algorithm' in:
N.A. Bruinsma and M. Steinbuch, 'A fast algorithm to compute the H∞-norm
of a transfer function matrix', Systems and Control Letters 14 (1990), pp. 287-293.
""" ->
function Base.norm(sys::StateSpace, p::Real=2; tol=1e-12)
    if p == 2
        return sqrt(trace(covar(sys, eye(size(sys.B, 2)))))
    elseif p == Inf
        if sys.Ts == 0
            return normHinf_twoSteps_ct(sys,tol)
        else
            error("H∞ norm computation not yet implemented for discrete-time systems!")
        end
    else
        error("`p` must be either `2` or `Inf`")
    end
end

function Base.norm(sys::TransferFunction, p::Real=2; tol=1e-12)
    return Base.norm(ss(sys),p)
end

function normHinf_twoSteps_ct(sys::StateSpace,tol=1e-12,maxIters=1000,approximag=1e-10)
    # `maxIters`: the maximum  number of iterations allowed in the algorithm (default 1000)
    # approximag is a tuning parameter: what does it mean for a number to be on the imaginary axis
    # Because of this tuning for example, the relative precision that we provide on the norm computation
    # is not a true guarantee, more an order of magnitude
    # outputs: pair of Float64, namely H∞ norm approximation and frequency fpeak at which it is achieved
    p = pole(sys)
    if any(real(p).>=0)
        # N.B.: more efficient than calling isstable, which would recompute the poles
        return Inf
    else
        # Initialization: computation of a lower bound from 3 terms
        lb = maximum(svdvals(sys.D)); fpeak = Inf
        (lb, idx) = findmax([lb, maximum(svdvals(evalfr(sys,0)))])
        if idx == 2
            fpeak = 0
        end
        p = pole(sys)
        if isreal(p)  # only real poles
            omegap = minimum(abs(p))
        else  # at least one pair of complex poles
            tmp = maximum(abs(imag(p)./(real(p).*abs(p))))
            omegap = abs(p[indmax(tmp)])
        end
        (lb, idx) = findmax([lb, maximum(svdvals(evalfr(sys,omegap*1im)))])
        if idx == 2
            fpeak = omegap
        end

        # Iterations
        iter = 1;
        while iter <= maxIters
            res = (1+2tol)*lb
            R = sys.D'*sys.D - res^2*eye(sys.nu)
            S = sys.D*sys.D' - res^2*eye(sys.ny)
            M = sys.A-sys.B*(R\sys.D')*sys.C
            H = [         M              -res*sys.B*(R\sys.B') ;
                   res*sys.C'*(S\sys.C)            -M'            ]
            omegas = eigvals(H)
            omegaps = imag(omegas[ (abs(real(omegas)).<=approximag) & (imag(omegas).>=0) ])
            if isempty(omegaps)
                return (1+tol)*lb, fpeak
            else  # if not empty, omegaps contains at least two values
                ms = [(x+y)/2 for x=omegaps[1:end-1], y=omegaps[2:end]]
                for mval in ms
                    (lb, idx) = findmax([lb, maximum(svdvals(evalfr(sys,mval*1im)))])
                            if idx == 2
                                fpeak = mval
                            end
                end
            end
            iter += 1
        end
        println("The computation of the H-infinity norm did not converge in $maxIters iterations")
    end
end

@doc """`T, B = balance(A[, perm=true])`

Compute a similarity transform `T` resulting in `B = T\\A*T` such that the row
and column norms of `B` are approximately equivalent. If `perm=false`, the
transformation will only scale, and not permute `A`.""" ->
function balance(A, perm::Bool=true)
    n = Base.LinAlg.checksquare(A)
    B = copy(A)
    job = perm ? 'B' : 'S'
    ilo, ihi, scaling = LAPACK.gebal!(job, B)

    S = diagm(scaling)
    for j = 1:(ilo-1)   S[j,j] = 1 end
    for j = (ihi+1):n   S[j,j] = 1 end

    P = eye(Int, n)
    if perm
        if ilo > 1
            for j = (ilo-1):-1:1 cswap!(j, round(Int, scaling[j]), P) end
        end
        if ihi < n
            for j = (ihi+1):n    cswap!(j, round(Int, scaling[j]), P) end
        end
    end
    return S, P, B
end

function cswap!{T<:Number}(i::Integer, j::Integer, X::StridedMatrix{T})
    for k = 1:size(X,1)
        X[i, k], X[j, k] = X[j, k], X[i, k]
    end
end

"""
`sysr, G = balreal(sys::StateSpace)`

Calculates a balance realization of the system sys, such that the observability and reachability gramians of the balanced system are equal and diagonal `G`

See also `gram`

Glad, Ljung, Reglerteori: Flervariabla och Olinjära metoder
"""
function balreal(sys::StateSpace)
P = gram(sys, :c)
Q = gram(sys, :o)
Q1 = chol(Q)
U,Σ,V = svd(Q1*P*Q1')
Σ = sqrt(Σ)
Σ1 = diagm((sqrt(Σ)))
T = Σ1\(U'Q1)

Pz = T*P*T'
Qz = inv(T')*Q*inv(T)
if vecnorm(Pz-Qz) > sqrt(eps())
    warn("balreal: Result may be inaccurate")
    println("Controllability gramian before transform")
    display(P)
    println("Controllability gramian after transform")
    display(Pz)
    println("Observability gramian before transform")
    display(Q)
    println("Observability gramian after transform")
    display(Qz)
    println("Singular values of PQ")
    display(Σ)
end

sysr = ss(T*sys.A/T, T*sys.B, sys.C/T, sys.D), diagm(Σ)
end


text:
import TableTraits
import NamedTuples

struct K_Table  <: AbstractDataFrame
    a::Array{K_,0}
    function K_Table(x::K_)
        a = asarray(x)
        t = xt(x)
        if(t != XT)
            throw(ArgumentError("type mismatch: t=$t ≠ $XT"))
        end
        return new(a)
    end
    function K_Table(columns::Vector, colnames::Vector{Symbol})
        ncols = length(columns)
        x = K_new(colnames)
        y = K_new(columns)
        a = asarray(xT(xD(x, y)))
        new(a)
    end
    function K_Table(; kwargs...)
        x = ktn(KS, 0)
        y = ktn(KK, 0)
        rx, ry = map(Ref{K_}, [x, y])
        for (k, v) in kwargs
            x = js(rx, ss(k))
            y = jk(ry, K_new(v))
        end
        a = asarray(xT(xD(x, y)))
        new(a)
    end
end
K_Table(df::AbstractDataFrame) = K_Table(K_new(df))
function K_Table(::Type{T}, n::Integer) where T <: NamedTuples.NamedTuple
    cols = fieldnames(T)
    x = K_new(cols)
    y = knk(length(cols), (ktn(ktypecode(S), n) for S in T.types)...)
    K_Table(xT(xD(x, y)))
end
kpointer(x::K_Table) = K_(pointer(x.a)-8)
valptr(x::K_Table, i) = unsafe_load(Ptr{K_}(xy(x.a[])+16), i)
colnames(x::K_Table) = K(r1(xx(x.a[])))
coldata(x::K_Table) = K(r1(xy(x.a[])))

Base.serialize(io::AbstractSerializer, x::K_Table) =
    _serialize(io, x, typeof(x))

DataFrames.ncol(x::K_Table) = Int(xn(xx(x.a[])))
DataFrames.nrow(x::K_Table) = Int(xn(valptr(x, 1)))
DataFrames.index(x::K_Table) = DataFrames.Index(Array(colnames(x)))
function DataFrames.names!(x::K_Table, vals; allow_duplicates=true)
    u = DataFrames.make_unique(vals, allow_duplicates=allow_duplicates)
    kS(kK(x.a[])[1])[:] = map(ss, u)
    x
end

Base.getindex(x::K_Table, i::Integer) = K(r1(valptr(x, i)))
Base.getindex(x::K_Table, i::Integer, j::Integer) = x[j][i]
Base.getindex(x::K_Table, i::Symbol) = x[DataFrames.index(x)[i]]

## IterableTable protocol

# T is the type of the elements produced
# TS is a tuple type that stores the columns of the DataFrame
struct K_Table_Iter{T, TS}
    x::K_Table
    # This field hodls a tuple with the columns of the DataFrame.
    # Having a tuple of the columns here allows the iterator
    # functions to access the columns in a type stable way.
    columns::TS
end

TableTraits.isiterable(x::K_Table) = true
TableTraits.isiterabletable(x::K_Table) = true

function TableTraits.getiterator(x::K_Table)
    names = colnames(x)
    columns = coldata(x)
    col_exprs = Expr[]
    columns_tuple_type = :(Tuple{})
    for i in 1:length(columns)
        push!(col_exprs, Expr(:(::), Symbol(names[i]), eltype(columns[i])))
        push!(columns_tuple_type.args, typeof(columns[i]))
    end
    named_tuple_type = NamedTuples.make_tuple(col_exprs)
    iter_type_expr = :(K_Table_Iter{$named_tuple_type,$columns_tuple_type})
    iter_type = eval(iter_type_expr)
    iter_type(x, (columns...))
end

Base.length{T,TS}(iter::K_Table_Iter{T,TS}) = size(iter.x,1)
Base.eltype{T,TS}(iter::K_Table_Iter{T,TS}) = T
Base.start{T,TS}(iter::K_Table_Iter{T,TS}) = 1

@generated function Base.next{T,TS}(iter::K_Table_Iter{T,TS}, state)
    constructor_call = :($T())
    for i in 1:length(iter.types[2].types)
        push!(constructor_call.args, :(iter.columns[$i][state]))
    end
    :($constructor_call, state + 1)
end

Base.done{T,TS}(iter::K_Table_Iter{T,TS}, state) = state > length(iter)

function K_Table(source)
    iter = TableTraits.getiterator(source)
    _table(Base.iteratorsize(iter), iter)
end

_table(::Base.IsInfinite, iter) = error("infinite source")
_table(::Base.HasShape, iter) = _table_prealloc(iter)
_table(::Base.HasLength, iter) = _table_prealloc(iter)

function _table_prealloc(iter)
    x = K_Table(eltype(iter), length(iter))
    for (i, row) in enumerate(iter)
        for (j, v) in enumerate(row)
            x[j][i] = v
        end
    end
    x
end

function _raw_eltype(T)
    exprs = [Expr(:(::), n, eltype(t))
        for (n, t) in zip(fieldnames(T), T.types)]
    NamedTuples.make_tuple(exprs)
end

function _table(::Base.SizeUnknown, iter)
    x = K_Table(eltype(iter), 0)
    for row in iter
        for (i, v) in enumerate(row)
            Q.coldata(x)[i] = push!(x[i], row[i])
        end
    end
    x
end


text:
Pkg.instantiate()

using Distributed, StanSample, CSV, DataFrames
Pkg.instantiate()

ProjDir = @__DIR__
cd(ProjDir)

rcd = CSV.read(joinpath(ProjDir, "RedCardData.csv"), DataFrame; missingstring = "NA")
rcd = dropmissing(rcd, :rater1)
rcd[!, :rater1] = Float64.(rcd.rater1)

stan_logistic_0 ="
data {
  int<lower=0> N;
  int<lower=0> n_redcards[N];
  int<lower=0> n_games[N];
  vector[N] rating;
}
parameters {
  vector[2] beta;
}
model {
  beta[1] ~ normal(0, 10);
  beta[2] ~ normal(0, 1);
  n_redcards ~ binomial_logit(n_games, beta[1] + beta[2] * rating);
}
";

stan_logistic_1 = "
functions {
  real partial_sum_lpmf(int[] slice_n_redcards,
                        int start, int end,
                        int[] n_games,
                        vector rating,
                        vector beta) {
    return binomial_logit_lupmf(slice_n_redcards |
                               n_games[start:end],
                               beta[1] + beta[2] * rating[start:end]);
  }
}
data {
  int<lower=0> N;
  int<lower=0> n_redcards[N];
  int<lower=0> n_games[N];
  vector[N] rating;
  int<lower=1> grainsize;
}
parameters {
  vector[2] beta;
}
model {

  beta[1] ~ normal(0, 10);
  beta[2] ~ normal(0, 1);
  target += reduce_sum(partial_sum_lupmf, n_redcards, grainsize,
                       n_games, rating, beta);
}
";

#isdir("tmp0") && rm("tmp0", recursive=true);
if !isdir("tmp0")
  tmp0 = mkdir(joinpath(ProjDir, "tmp0"))
else
  tmp0 = joinpath(ProjDir, "tmp0")
end;

logistic_0 = SampleModel("logistic0", stan_logistic_0, [1]; tmpdir=tmp0)

#isdir("tmp1") && rm("tmp1", recursive=true);
if !isdir("tmp1")
  tmp1 = mkdir(joinpath(ProjDir, "tmp1"))
else
  tmp1 = joinpath(ProjDir, "tmp1")
end;

logistic_1 = SampleModel("logistic1", stan_logistic_1; tmpdir=tmp1)

data = Dict(
    :N => size(rcd, 1),
    :grainsize => 1,
    :n_redcards => rcd.redCards,
    :n_games => rcd.games,
    :rating => rcd.rater1
)

for i in 1:1
  #addprocs(1)
  println("\nUsing $(nprocs()) cores.\n")
  @time rc_0 = stan_sample(logistic_0; data);

  if success(rc_0)
      dfs_0 = read_summary(logistic_0)
      dfs_0[8:9, [1,2,4,8,9,10]] |> display
      println()
  end

  @time rc_1 = stan_sample(logistic_1; data);

  if success(rc_1)
      dfs_1 = read_summary(logistic_1)
      dfs_1[8:9, [1,2,4,8,9,10]] |> display
      println()
  end
end


text:
#= using UnicodePlots =#

# Set up model
function build_model()
    F = diagm([1.0])
    V = diagm([2.0])
    G = reshape([1, 2, -0.5], 3, 1)
    W = diagm([8.0, 2.5, 4.0])
    x0 = randn(1)
    P0 = diagm([1e7])
    StateSpaceModel(F, V, G, W, x0, P0)
end

# Time varying model
function sinusoid_model(omega::Real; fs::Int=256, x0=[0.5, -0.5], W::AbstractFloat=0.1)
    F(n)  = [1.0 0; 0 1.0]
    V(n)  = diagm([1e-10, 1e-10])
    G(n)  = [cos(2*pi*omega*(1/fs)*n) -sin(2*pi*omega*(1/fs)*n)]
    w(n)  = diagm([W])
    P0    = diagm([1e-1, 1e-1])
    StateSpaceModel(F, V, G, w, x0, P0)
end

@testset "Kalman Filter" begin

    srand(1)

    @testset "Time invariant models" begin

        @testset "Building model" begin
            build_model()
        end

        @testset "Simulations" begin
            mod1 = build_model()
            x, y = TimeModels.simulate(mod1, 100)
        end

        @testset "Filtering" begin
            mod1 = build_model()
            x, y = TimeModels.simulate(mod1, 100)
            filt = kalman_filter(y, mod1)
        end

        @testset "Smoothing" begin
            mod1 = build_model()
            x, y = TimeModels.simulate(mod1, 100)

            smooth = kalman_filter(y, mod1) |> kalman_smooth
            smooth2 = kalman_smooth(y, mod1)

            @test smooth.smoothed ≈ smooth2.smoothed atol=1e-2
        end

        @testset "Missing data" begin
            mod1 = build_model()
            x, y = simulate(mod1, 100)
            y[1:9:end] = NaN
            y[100] = NaN
            filt = kalman_filter(y, mod1)
            smooth = kalman_smooth(filt)
            smooth2 = kalman_smooth(y, mod1)
            @test !any(isnan, filt.filtered)
            @test !any(isnan, smooth.filtered)
            @test !any(isnan, smooth2.smoothed)
        end

        @testset "Return sizes" begin
            mod1 = build_model()
            x, y = TimeModels.simulate(mod1, 100)
            filt = kalman_filter(y, mod1)
            smooth = kalman_smooth(filt)
            smooth2 = kalman_smooth(y, mod1)

            @test size(filt.filtered) == size(filt.predicted)
            @test size(filt.filtered) == size(smooth.filtered)
            @test size(filt.filtered) == size(smooth.smoothed)
            @test size(filt.error_cov) == size(smooth.error_cov)
            @test size(filt.filtered) == size(smooth2.smoothed)
            @test size(filt.error_cov[:,:,1]) == size(smooth2.error_cov[1])
            @test size(filt.error_cov, 3) == length(smooth2.error_cov)
        end

        @testset "Linear regression test" begin
            m, b, s, dt = 5, 2, 2, .1
            t = 0:dt:10
            y_true = m*t + b
            input = 100*[sin.(t./2) sin.(t./4) cos.(t./2) cos.(t./4)] .+ 10
            y_noisy = [y_true zeros(length(t)) -y_true] +
                        100*[sin.(t./2) + sin.(t./4) sin.(t./2) + cos.(t./2) cos.(t./2) + cos.(t./4)] .+ 10 .+ randn(length(t), 3)
            lm = StateSpaceModel([1 dt; 0 1], zeros(2,2),
                                  [1. 0; 0 0; -1 0], s*eye(3),
                                  zeros(2), 100*eye(2),
                                  B=zeros(2, 4), D=[1. 1 0 0; 1 0 1 0; 0 0 1 1])
            lm_filt = kalman_filter(y_noisy, lm, u=input)
            @test lm_filt.filtered[end,1] ≈ y_true[end, 1] atol=3*sqrt(lm_filt.error_cov[1,1,end])

            lm_smooth = kalman_smooth(lm_filt)
            stderr = sqrt.(lm_smooth.error_cov[1:1,1:1,:][:])
            @test lm_filt.filtered[end:end,:] == lm_smooth.smoothed[end:end,:]
            @test all(abs.(y_true - lm_smooth.smoothed[:,1]) .< 3*stderr)
            @test ones(t) * lm_smooth.smoothed[1,2] ≈ lm_smooth.smoothed[:, 2] atol=1e-12

            # Repeat with DK smoother
            lm_smooth = kalman_smooth(y_noisy, lm, u=input)
            stderr = sqrt.(Float64[P[1,1] for P in lm_smooth.error_cov])
            @test all(abs.(y_true - lm_smooth.smoothed[:,1]) .< 3*stderr)
            @test ones(t) * lm_smooth.smoothed[1,2] ≈ lm_smooth.smoothed[:, 2] atol=1e-12

        end

    end

    @testset "Time varying models" begin

        @testset "Building Model" begin
            sinusoid_model(40)
        end

        @testset "Simulations" begin
            srand(1)
            fs = 256
            mod2 = sinusoid_model(4, fs = 256)
            x, y = TimeModels.simulate(mod2, fs*2)

            #= display(lineplot(collect(1:size(x, 1)) / fs, vec(y), width = 120, title = "Original Data")) =#
        end

        @testset "Filtering" begin
            srand(1)
            fs = 256
            mod2 = sinusoid_model(4, fs = 256)
            x, y = TimeModels.simulate(mod2, fs*2)

            @testset "Correct initial guess" begin
                filt = kalman_filter(y, mod2)
                @test filt.predicted[end, :] ≈ [0.5, -0.5] atol=0.3
            end

            @testset "Incorrect initial guess" begin
                mod3 = sinusoid_model(4, fs = 256, x0=[1.7, -0.2])
                filt = kalman_filter(y, mod3)
                @test filt.predicted[end, :] ≈ [0.5, -0.5] atol= 0.3

            end

            @testset "Model error" begin
                mod4 = sinusoid_model(4, fs = 256, x0=[1.7, -0.2], W=3.0)
                filt = kalman_filter(y, mod4)
                @test filt.predicted[end, :] ≈ [0.5, -0.5] atol= 0.3
            end

            @testset "Standalone log-likelihood function works" begin
                loglik1 = kalman_filter(y, mod2).loglik
                loglik2 = loglikelihood(y, mod2)
                @test loglik1 ≈ loglik2
            end

        end

        @testset "Smoothing" begin
            srand(1)
            fs = 256
            mod2 = sinusoid_model(4, fs = 8192)
            x, y = TimeModels.simulate(mod2, fs*10)
            smooth = kalman_smooth(y, sinusoid_model(4, fs = 8192, x0=[1.7, -0.2]) )
            @test mean(smooth.smoothed, 1) ≈ [0.5 -0.5] atol= 0.1

            #= x_est = round(smooth.smoothed[end:end, :], 3) =#
            #= display(lineplot(collect(1:size(x, 1)) / fs, vec(smooth.smoothed[1:end, 1]), width = 120, title="Smoothed State 1: $(x_est[1])")) =#
            #= display(lineplot(collect(1:size(x, 1)) / fs, vec(smooth.smoothed[1:end, 2]), width = 120, title="Smoothed State 2: $(x_est[2])")) =#
        end

        @testset "Sparse regression test" begin
            # To be fair, not an exemplary example of sparse matrices in action
            # since the state covariance off-diagonals are small but nonzero
            n, m, s = 5000, 10, 5.
            x = 100*rand(n, m)
            coeffs = randn(m)
            y = x * coeffs + s*randn(n)

            I_m   = speye(m)
            I0_m  = spzeros(m,m)
            S     = diagm(s)
            mlm   = StateSpaceModel(
                _->I_m, # A
                _->I0_m, # V
                t->x[t:t,:], # C
                _->S, # W
                zeros(m), # x1
                speye(m) # P1
            )
            mlm_smooth = kalman_smooth(y, mlm)
            @test ones(n,m) .* mlm_smooth.smoothed[1:1,:] ≈ mlm_smooth.smoothed atol=1e-12
            @test vec(mlm_smooth.smoothed[1:1,:]) ≈ coeffs atol=1e-2

            @testset "Correct failure" begin
                @test_throws MethodError StateSpaceModel([1 0.1; 0 1], zeros(2,3), zeros(2,2),
                    [1. 0; 0 0; -1 0], [1. 1 0 0; 1 0 1 0; 0 0 1 1], eye(3), zeros(2), 100*eye(2))
            end
        end

    end

end


text:
#= Load these if running these tests in terminal
import Pkg; Pkg.activate("../")  # CytofDensityEstimation
using CytofDensityEstimation; const CDE=CytofDensityEstimation
using Distributions

include("runtests.jl")
=#

@testset "metropolis adaptive" begin
  easy = false
  if easy
    K = 3
    S = CDE.MCMC.eye(K) * 1.0
    S[1,2] = S[2,1] = -0.8
    S[1,3] = S[3,1] = 0.5
  else
    K = 32
    S = let
      M = randn(K,K)
      M*M'
    end
  end
  inv_S = inv(S)
  log_prob(x) = -x' * inv_S * x / 2.0
  mvn = MvNormal(S)

  tuner = CDE.MCMC.MvTuner(K)
  xnew = randn(K)
  out = Vector{Float64}[]
  @time for i in 1:100000
    xnew = CDE.MCMC.metropolisAdaptive(xnew, log_prob, tuner)
    i > 50000 && append!(out, [xnew])
  end
  println("$((cov(out)[1], S[1]))")
  println("$((cov(out)[2], S[2]))")
  # plot([getindex.(out, 1) getindex.(out, 2)])
  # plot(getindex.(out, 1), getindex.(out, 4))
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation