text:
using AnimatedPlots

function circle(x)
	10cos(x), 10sin(x)
end

g = StaticGraph(circle)

plot(g)

waitfor(current_window())


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct DoubleValue <: ProtoType
    value::Float64
    DoubleValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct DoubleValue

mutable struct FloatValue <: ProtoType
    value::Float32
    FloatValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct FloatValue

mutable struct Int64Value <: ProtoType
    value::Int64
    Int64Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Int64Value

mutable struct UInt64Value <: ProtoType
    value::UInt64
    UInt64Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UInt64Value

mutable struct Int32Value <: ProtoType
    value::Int32
    Int32Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct Int32Value

mutable struct UInt32Value <: ProtoType
    value::UInt32
    UInt32Value(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UInt32Value

mutable struct BoolValue <: ProtoType
    value::Bool
    BoolValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct BoolValue

mutable struct StringValue <: ProtoType
    value::AbstractString
    StringValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct StringValue

mutable struct BytesValue <: ProtoType
    value::Array{UInt8,1}
    BytesValue(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct BytesValue

export DoubleValue, FloatValue, Int64Value, UInt64Value, Int32Value, UInt32Value, BoolValue, StringValue, BytesValue

ext: jl


text:
module SimpleGraphAlgorithms
using SimpleGraphs
using JuMP
using Cbc, ChooseOptimizer

_SOLVER = Cbc       # this is used by JuMP
_OPTS = Dict()        # to pass options to JuMP optimizers

export use_Cbc, use_Gurobi, use_optimizer

"""
`use_optimizer(choice::Module, options::Dict)` chooses the
optimization package and package options to be used by functions
in `SimpleGraphAlgorithms` (and `SimplePosetAlgorithms`).

See `use_Cbc` and `use_Gurobi`.
"""
function use_optimizer(choice::Module, options::Dict = Dict())
    global _SOLVER = choice
    global _OPTS = deepcopy(options)

    nothing
end

"""
`use_Gurobi(verbose=false)` sets the optimization software
to be the `Gurobi` solver. Requires the user to invoke
`using Gurobi` first.

See also `use_Cbc`.
"""
function use_Gurobi(verbose::Bool=false)
    set_solver(Main.Gurobi)
    set_solver_verbose(verbose)
    nothing
end

"""
`use_Cbc(verbose=false)` sets the optimization software
to be the `Cbc` solver. This is set by default when `SimpleGraphAlgorithms`
is loaded.

See also `use_Gurobi`.
"""
function use_Cbc(verbose::Bool=false)
    set_solver(Cbc)
    set_solver_verbose(verbose)
    nothing
end

# Use Cbc in nonverbose mode until told otherwise
use_Cbc()

export max_indep_set, max_clique, max_matching, min_dom_set
export min_vertex_cover, min_edge_cover

import SimpleGraphs.cache_save

"""
`max_indep_set(G)` returns a maximum size independent set of a
`SimpleGraph`.
"""
function max_indep_set(G::SimpleGraph)
    if cache_check(G,:max_indep_set)
        return cache_recall(G,:max_indep_set)
    end
    VV = vlist(G)
    EE = elist(G)
    n = NV(G)
    m = NE(G)

    MOD = Model(get_solver())
    @variable(MOD, x[VV],Bin)
    for e in EE
        u,v = e
        @constraint(MOD,x[u]+x[v]<=1)
    end
    @objective(MOD,Max,sum(x[v] for v in VV))
    optimize!(MOD)

    X = value.(x)
    A = Set([v for v in VV if X[v]>0.1])

    cache_save(G,:max_indep_set, A)
    return A
end

"""
`min_vertex_cover(G)` returns a smallest vertex cover `S` of `G`.
This is a smallest possible set of vertices such that every edge
of `G` has one or both end points in `S`.

`min_vertex_cover(G,k)` returns the smallest set of vertices `S`
such that at least `k` edges are indicent with a vertex in `S`.
"""
min_vertex_cover(G) = setdiff(G.V, max_indep_set(G))

function min_vertex_cover(G::SimpleGraph, k::Int)
    MOD = Model(get_solver())
    Vs = vlist(G)
    Es = elist(G)

    @variable(MOD, x[v in Vs], Bin)
    @variable(MOD, y[f in Es], Bin)

    # edge constraints
    for f in Es
        u,v = f
        @constraint(MOD, x[u]+x[v] >= y[f])
    end

    # assure proper size
    @constraint(MOD, sum(y[f] for f in Es) >= k)

    # we want to minimize the sum of the x[v]'s
    @objective(MOD, Min, sum(x[v] for v in Vs))

    optimize!(MOD)
    status = Int(termination_status(MOD))
    if status != 1
        error("Cannot find such a minimum vertex cover")
    end

    XX = value.(x)
    A = Set([v for v in Vs if XX[v]> 0.1])
    return A
end

"""
`max_clique(G)` returns a maximum size clique of a `SimpleGraph`.
"""
function max_clique(G::SimpleGraph)
  if cache_check(G,:max_clique)
    return cache_recall(G,:max_clique)
  end
  result = max_indep_set(G')
  cache_save(G,:max_clique,result)
  return result
end

"""
`max_matching(G)` returns a maximum matching of a `SimpleGraph`.
"""
function max_matching(G::SimpleGraph)
    if cache_check(G,:max_matching)
      return cache_recall(G,:max_matching)
    end
    VV = vlist(G)
    EE = elist(G)
    n = NV(G)
    m = NE(G)

    MOD = Model(get_solver())
    @variable(MOD, x[EE], Bin)

    for v in VV
        star = [e for e in EE if e[1]==v || e[2]==v]
        @constraint(MOD, sum(x[e] for e in star) <= 1)
    end

    @objective(MOD, Max, sum(x[e] for e in EE))

    optimize!(MOD)
    X = value.(x)
    M = Set([e for e in EE if X[e]>0.1])

    cache_save(G,:max_matching,M)
    return M
end

"""
`min_edge_cover(G)` finds a smallest subset `F` of the edges such that
every vertex of `G` is the end point of at least one edge in `F`. An
error is raised if `G` has a vertex of degree 0.
"""
function min_edge_cover(G::SimpleGraph)
    if cache_check(G,:min_edge_cover)
      return cache_recall(G,:min_edge_cover)
    end

    if NV(G) == 0
      T = vertex_type(G)
      S = Tuple{T,T}
      return Set{S}()
    end

    if minimum(deg(G))==0
        error("Graph has an isolated vertex; no minimum edge cover exists.")
    end

    VV = vlist(G)
    EE = elist(G)
    n = NV(G)
    m = NE(G)

    MOD = Model(get_solver())
    @variable(MOD, x[EE], Bin)

    for v in VV
        star = [e for e in EE if e[1]==v || e[2]==v]
        @constraint(MOD, sum(x[e] for e in star) >= 1)
    end

    @objective(MOD, Min, sum(x[e] for e in EE))

    optimize!(MOD)
    X = value.(x)
    A = Set([e for e in EE if X[e]>0.1])

    cache_save(G,:min_edge_cover,A)
    return A
end

"""
`min_dom_set(G)` returns a smallest dominating set of a
`SimpleGraph`. That is, a smallest set `S` with the property that
every vertex of `G` either is in `S` or is adjacent to a vertex of
`S`.
"""
function min_dom_set(G::SimpleGraph)
    if cache_check(G,:min_dom_set)
      return cache_recall(G,:min_dom_set)
    end
    n = NV(G)
    if n==0
      T = vertex_type(G)
      return Set{T}()
    end

    VV = vlist(G)

    MOD = Model(get_solver())
    @variable(MOD, x[VV], Bin)

    for v in VV
        Nv = G[v]
        append!(Nv,v)
        @constraint(MOD, sum(x[w] for w in Nv) >= 1)
    end

    @objective(MOD, Min, sum(x[v] for v in VV))
    optimize!(MOD)

    X = value.(x)

    S = Set(v for v in VV if X[v]>0.1)

    cache_save(G,:min_dom_set,S)
    return S
end

include("iso.jl")
include("kcolor.jl")
include("chrome_poly.jl")
include("edge-coloring.jl")
include("mad.jl")
include("kfactor.jl")
include("connectivity.jl")

end # of module SimpleGraphAlgorithms


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TestPkg

using Random

end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
# using Base.Test
# using MIToS.MSA

# Fields of MultipleSequenceAlignment
const msa_fields = Symbol[:id, :msa, :annotations]

# Fields of AlignedAlignedSequence
const seq_fields = Symbol[:id, :index, :sequence, :annotations]

print("""

Tests for Multiple Sequence Alignments
======================================
""")

print("""

Parse Pfam
----------
""")

print("""
Test pfam stockholm parser using the 4 sequence full MSA for PF09645
> Order: Tree
> Inserts lower case
> Gaps as "." or "-" (mixed)
> Pfam version 28.0, based on UniProt release 2014_07
""")

const pfam = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm)
const pfam_na = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm, MultipleSequenceAlignment)
const F112_SSV1 = collect(".....QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ.....")

let rawpfam = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm, Matrix{Residue})
  @test rawpfam == pfam.msa
  @test names(rawpfam) == IndexedArrays.IndexedArray(["1","2","3","4"])
end

@test pfam_na == convert(MultipleSequenceAlignment, pfam)

@test size(pfam.msa, 1) == 4
@test size(pfam.msa, 2)  == length(F112_SSV1[ F112_SSV1 .!= '.' ]) # Without inserts
@test slice(pfam.msa,4,:) == convert(Vector{Residue}, F112_SSV1[ F112_SSV1 .!= '.' ])
@test pfam.id.items == ["C3N734_SULIY/1-95", "H2C869_9CREN/7-104", "Y070_ATV/2-70", "F112_SSV1/3-112"]
@test !isempty(pfam.annotations)
@test length(pfam.annotations.columns) == 2
@test length(pfam.annotations.residues) == 1
@test pfam.annotations.residues[("F112_SSV1/3-112","SS")] == "X---HHHHHHHHHHHHHHHSEE-HHHHHHHH---HHHHHHHHHHHHHHHHH-TTTEEEEE-SS-EEEEE--XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
@test pfam.annotations.columns["seq_cons"] == "...NshphAclhaKILppKtElolEDIlAQFEISsosAYsI.+sL+hICEpH.-ECpsppKsRKTlhh.hKpEphppptpEp..ppItKIhsAp................"
@test pfam.annotations.sequences[("F112_SSV1/3-112","DR")] == "PDB; 2VQC A; 4-73;"

print("""
Test parse for string inputs
""")

const pfam_string = """# STOCKHOLM 1.0
#=GS C3N734_SULIY/1-95   AC C3N734.1
#=GS H2C869_9CREN/7-104  AC H2C869.1
#=GS Y070_ATV/2-70       AC Q3V4T1.1
#=GS F112_SSV1/3-112     AC P20220.1
#=GS F112_SSV1/3-112     DR PDB; 2VQC A; 4-73;
C3N734_SULIY/1-95              ...mp---NSYQMAEIMYKILQQKKEISLEDILAQFEISASTAYNVQRTLRMICEKHPDECEVQTKNRRTIFKWIKNEETTEEGQEE--QEIEKILNAQPAE-------------k....
H2C869_9CREN/7-104             ...nk--LNDVQRAKLLVKILQAKGELDVYDIMLQFEISYTRAIPIMKLTRKICEAQ-EICTYDEKEHKLVSLNAKKEKVEQDEEENEREEIEKILDAH----------------trreq
Y070_ATV/2-70                  qsvne-------VAQQLFSKLREKKEITAEDIIAIYNVTPSVAYAIFTVLKVMCQQHQGECQAIKRGRKTVI-------------------------------------------vskq.
F112_SSV1/3-112                .....QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ.....
#=GR F112_SSV1/3-112     SS    .....X---HHHHHHHHHHHHHHHSEE-HHHHHHHH---HHHHHHHHHHHHHHHHH-TTTEEEEE-SS-EEEEE--XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX.....
#=GC SS_cons                   .....X---HHHHHHHHHHHHHHHSEE-HHHHHHHH---HHHHHHHHHHHHHHHHH-TTTEEEEE-SS-EEEEE--XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX.....
#=GC seq_cons                  ........NshphAclhaKILppKtElolEDIlAQFEISsosAYsI.+sL+hICEpH.-ECpsppKsRKTlhh.hKpEphppptpEp..ppItKIhsAp................h....
//
"""

@test parse(pfam_string, Stockholm) == pfam

print("""

Keep insert columns
-------------------
""")

let aln = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm, keepinserts=true)
    @test (collect(getannotcolumn(aln, "Aligned")) .== '1') == (F112_SSV1 .!= '.') # Aligned columns
    @test asciisequence(aln, 1) == replace(uppercase("...mp---NSYQMAEIMYKILQQKKEISLEDILAQFEISASTAYNVQRTLRMICEKHPDECEVQTKNRRTIFKWIKNEETTEEGQEE--QEIEKILNAQPAE-------------k...."), '.', '-')
end

print("""

Parse Fasta
-----------
""")

const fasta = read(joinpath(pwd(), "data", "PF09645_full.fasta.gz"), FASTA)
const small = read(joinpath(pwd(), "data", "Gaoetal2011.fasta"), FASTA)
const fasta_na = read(joinpath(pwd(), "data", "PF09645_full.fasta.gz"), FASTA, MultipleSequenceAlignment)
const small_na = read(joinpath(pwd(), "data", "Gaoetal2011.fasta"), FASTA, MultipleSequenceAlignment, deletefullgaps=false)

let rawfasta = read(joinpath(pwd(), "data", "PF09645_full.fasta.gz"), FASTA, Matrix{Residue})
  @test rawfasta == fasta.msa
end

@test fasta_na == convert(MultipleSequenceAlignment, fasta)
@test small_na == convert(MultipleSequenceAlignment, small)

print("""
Test download of read(URL, ...)
""")

@test read("https://raw.githubusercontent.com/diegozea/MIToS.jl/master/test/data/Gaoetal2011.fasta", FASTA) == small
@test read("https://raw.githubusercontent.com/diegozea/MIToS.jl/master/test/data/PF09645_full.fasta.gz", FASTA) == fasta

print("""
Test parse for string inputs
""")

const fasta_string = """>SEQ1
DAWAEE
>SEQ2
DAWAEF
>SEQ3
DAWAED
>SEQ4
DAYCMD
>SEQ5
DAYCMT
>SEQ6
DAYCMT
"""

@test parse(fasta_string, FASTA) == small

print("""
Test the FASTA parser usando Gao et.al. 2011 example
""")

const seq1 = res"DAWAEE"
const seq2 = res"DAWAEF"
const seq3 = res"DAWAED"
const seq4 = res"DAYCMD"
const seq5 = res"DAYCMT"
const seq6 = res"DAYCMT"
const gaores = [seq1'; seq2'; seq3'; seq4'; seq5'; seq6']

@test small.msa == gaores
@test small.id.items == ASCIIString[ "SEQ$i" for i in 1:6 ]

print("""
Sequence identifiers FASTA
""")

let file = joinpath(pwd(), "data", "ids.fasta")
    @test names(parse(readall(file), FASTA)) == names(read(file, FASTA))
end

print("""
Are the parsers (pfam/fasta) getting the same result?
""")

for field in msa_fields
  if field != :annotations
    @eval @test fasta.$field == pfam.$field
  end
end

print("""
Test parse with ambiguous or not standard residues (are gaps on MIToS)
""")

let default = read(joinpath(pwd(), "data", "alphabet.fasta"), FASTA, generatemapping=true),
    notused = read(joinpath(pwd(), "data", "alphabet.fasta"), FASTA, checkalphabet=true) # deletenotalphabetsequences! is called also for Raw & Stockholm

  @test Base.vec(default[1,:]) == res"ARNDCQEGHILKMFPSTWYV"
  for i in 2:nsequences(default)
    @test Base.vec(default[i,:]) == res"AR-DCQEGHILKMFPSTWYV"
    @test  getsequencemapping(default,1) == getsequencemapping(default,i)
  end

  @test nsequences(notused) == 1
  @test nsequences(read(joinpath(pwd(), "data", "alphabet.fasta"), FASTA, Matrix{Residue}, checkalphabet=true)) == 1 # _strings_to_msa is called also for Raw & Stockholm
end

print("""
Dots are part of the alphabet for input.
""")

let msa = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm, checkalphabet=true)
  @test size(pfam.msa, 1) == 4
end

print("""

Selection without Mappings
--------------------------
""")

@test getresidues(small) == gaores
@test getresidues(small_na) == gaores

print("""
Test getindex
""")

for index in [ 2, (2,:), (:,2), (2,2), (2:4, 2:4), (2,[3,4,5]),
              (:, [1,2,3,4,5,6] .< 4 ), ([1,2,3,4,5,6] .< 4, :),
              (:,:), ([1,2,3,4,5,6] .< 4, [1,2,3,4,5,6] .< 4) ]
  @test small[index...] == gaores[index...]
  @test small_na[index...] == gaores[index...]
end

print("""
Test setindex! with copy and deepcopy
""")
const copy_small = copy(small)
const deepcopy_small_na = deepcopy(small_na)

for (index, value) in [ ((1),Residue('H')), # ((1),Residue('D')),
                        ((:,1), res"HHHHHH"), # ((:,1), res"DDDDDD"),
                        ((1),Residue('H')) ] #, ((1),Residue('D')) ]
  copy_small[index...] = value
  @test copy_small[index...] == value
  @test small[index...] != value
  deepcopy_small_na[index...] = value
  @test deepcopy_small_na[index...] == value
  @test small_na[index...] != value
end

for (index, value) in [ (1,Residue('H')), (1,Residue('D')),
                        ( : , res"HHHHHH"), ( : , res"DAYCMD") ]
  tmp_sequence = copy(getsequence(small, 4))
  tmp_sequence[index] = value
  @test tmp_sequence[index] == value
end

print("""
Test getresidues, getsequence and getresiduesequence
""")

for i in 1:6
  println("Test sequence $i for getresidues and getsequence")
  @eval @test getresidues(getsequence(small,$i)) == $(Symbol("seq$i"))
  @eval @test getresidues(getsequence(small_na,$i)) == $(Symbol("seq$i"))
end

const rawseqs = getresiduesequences(small)
const rawseqs_na = getresiduesequences(small_na)

for i in 1:6
  println("Test sequence $i from getresiduesequences")
  @eval @test rawseqs[$i] == $(Symbol("seq$i"))
  @eval @test rawseqs_na[$i] == $(Symbol("seq$i"))
end

print("""

Size, sequences and columns
---------------------------
""")

for aln in (small, small_na)
  @test size(aln) == (6,6)
  @test length(aln) == 36
  @test ncolumns(aln) == 6
  @test nsequences(aln) == 6
end

@test ncolumns(getsequence(small, 4)) == 6

for aln in (fasta, fasta_na, pfam, pfam_na)
  @test size(aln) == (4,110)
  @test length(aln) == 440
  @test ncolumns(aln) == 110
  @test nsequences(aln) == 4
end

print("""

AnnotatedAlignedSequence and AlignedSequence
--------------------------------------------
""")

print("""
Test getsequence using Index and ID and convert to AlignedSequence
""")

for i in 1:4
  annseq = getsequence(pfam, i)
  @test isa(annseq, AnnotatedAlignedSequence)
  seq = getsequence(pfam_na, i)
  @test isa(seq, AlignedSequence)
  @test AlignedSequence(annseq) == seq
end

for id in ["C3N734_SULIY/1-95", "H2C869_9CREN/7-104", "Y070_ATV/2-70", "F112_SSV1/3-112"]
  annseq = getsequence(pfam, id)
  @test pfam[id] == annseq
  @test isa(annseq, AnnotatedAlignedSequence)
  seq = getsequence(pfam_na, id)
  @test pfam_na[id] == seq
  @test isa(seq, AlignedSequence)
  @test AlignedSequence(annseq) == seq
end

print("""

Test Filters!
-------------
""")

print("""
filtersequences! with and without annotations
""")
@test getsequence(filtersequences!(copy(pfam), [1,2,3,4] .== 3), 1) == getsequence(pfam,3)
@test getsequence(filtersequences!(copy(pfam), [1,2,3,4] .> 2), 2) == getsequence(pfam,4)

@test getsequence(filtersequences!(copy(pfam_na), [1,2,3,4] .== 3), 1) == getsequence(pfam_na,3)
@test getsequence(filtersequences!(copy(pfam_na), [1,2,3,4] .> 2), 2) == getsequence(pfam_na,4) # BitVector

@test getsequence(filtersequences!(copy(pfam_na), Bool[false,false,true,true] ), 2) == getsequence(pfam_na,4) # Vector{Bool}

print("""
filtercolumns! with and without annotations
""")

@test getresidues(getsequence(filtercolumns!(copy(pfam), collect(1:110) .<= 10), 4)) == res"QTLNSYKMAE"
@test getresidues(getsequence(filtercolumns!(copy(pfam_na), collect(1:110) .<= 10), 4)) == res"QTLNSYKMAE"

print("""
Test for BoundsError's on filters with bad masks
""")

@test_throws BoundsError filtersequences!(copy(pfam), [1,2,3,4,5,6,7,8,9,10] .> 2)
@test_throws BoundsError filtersequences!(copy(pfam_na), [1,2,3] .> 2)
@test_throws BoundsError filtercolumns!(copy(pfam), [1,2,3] .> 2)
@test_throws BoundsError filtercolumns!(copy(pfam_na), collect(1:200) .<= 10)

print("""
Test copy and deepcopy for sequences
""")

for (index, value) in [ (2,Residue('H')) , (2:3,res"HG") ]
  annseq = getsequence(pfam, 4)
  seq = getsequence(pfam_na, 4)
  copied_annseq = copy(annseq)
  deepcopied_seq = deepcopy(seq)
  copied_annseq[index] = value
  deepcopied_seq[index] = value
  @test copied_annseq[index] == value
  @test annseq[index] != value
  @test deepcopied_seq[index] == value
  @test seq[index] != value
end

print("""
Test filtercolumns! for sequences
""")

let i = 4
  annseq = getsequence(pfam, i)
  seq = getsequence(pfam_na, i)
  print("""
  Test Sequence .== Residue for mask generation
  """)
  @test annseq[annseq .== Residue('Q')] == res"QQQQQQQQQQQQQQ"
  @test seq[seq .== Residue('Q')] == res"QQQQQQQQQQQQQQ"
  print("""
  Test filtercolumns!
  """)
  filtered_annseq = filtercolumns!(copy(annseq), annseq .== Residue('Q'))
  filtered_seq = filtercolumns!(copy(seq), seq .== Residue('Q'))
  @test_throws BoundsError filtercolumns!(copy(annseq), collect(1:(length(seq)-10)) .> 2)
  @test_throws BoundsError filtercolumns!(copy(seq), collect(1:(length(seq)+10)) .<= 10)
  @test getresidues(filtered_annseq) == res"QQQQQQQQQQQQQQ"
  @test getresidues(filtered_seq) == res"QQQQQQQQQQQQQQ"
  @test AlignedSequence(filtered_annseq) == filtered_seq
  @test getannotcolumn(filtered_annseq, "SS_cons") == "XHHEXXXXXXXXXX"
  @test getannotresidue(filtered_annseq, "F112_SSV1/3-112", "SS") == "XHHEXXXXXXXXXX"
end

print("""

Test Annotation Getters
-----------------------
""")

@test getannotcolumn(pfam, "SS_cons") == getannotcolumn(getsequence(pfam, 4), "SS_cons")
@test getannotresidue(pfam, "F112_SSV1/3-112", "SS") == getannotresidue(getsequence(pfam, 4), "F112_SSV1/3-112", "SS")
@test getannotfile(pfam) == getannotfile(getsequence(pfam, 4))
@test getannotsequence(pfam,"F112_SSV1/3-112","DR") == getannotsequence(getsequence(pfam, 4), "F112_SSV1/3-112", "DR")

print("""

Test Gaps and Coverage
----------------------
""")

@test gapfraction(small) == 0.0
@test gapfraction(getsequence(small,1)) == 0.0
@test gapfraction(small[1,:]) == 0.0

@test gapfraction(getsequence(pfam,4)) == 0.0
@test gapfraction(getsequence(pfam_na,4)) == 0.0

@test gapfraction(convert(Vector{Residue}, F112_SSV1)) == mean(F112_SSV1 .== '.')

@test gapfraction(pfam[:,1]) == 0.75
@test gapfraction(pfam_na[:,1]) == 0.75

print("""
Test residuefraction
""")

@test residuefraction(small) == 1.0
@test residuefraction(getsequence(pfam, 4)) == 1.0
@test residuefraction(pfam[:,1]) == 0.25

print("""
Test coverage
""")

@test coverage(pfam)[4] == 1.0
@test sum(coverage(small)) == 6.0

print("""

Reference
---------
""")

print("""
Test setreference!
""")

let copy_pfam = copy(pfam), copy_pfam_na = copy(pfam_na)
  setreference!(copy_pfam, 4)
  @test getresidues(getsequence(copy_pfam,1)) == res"QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ"
  setreference!(copy_pfam, "C3N734_SULIY/1-95")
  @test getresidues(getsequence(copy_pfam,4)) == res"QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ"
  @test copy_pfam == pfam

  setreference!(copy_pfam_na, 4)
  @test getresidues(getsequence(copy_pfam_na,1)) == res"QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ"
  setreference!(copy_pfam_na, "C3N734_SULIY/1-95")
  @test getresidues(getsequence(copy_pfam_na,4)) == res"QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ"
  @test copy_pfam_na == pfam_na
end

print("""
Test gapstrip!
""")

let copy_pfam = copy(pfam), copy_pfam_na = copy(pfam_na)
  setreference!(copy_pfam, 4)
  setreference!(copy_pfam_na, 4)
  gapstrip!(copy_pfam, gaplimit=1.0, coveragelimit=0.0)
  gapstrip!(copy_pfam_na, gaplimit=1.0, coveragelimit=0.0)
  @test size(copy_pfam) == (4, 110)
  @test copy_pfam == copy_pfam_na

  setreference!(copy_pfam, 1)
  setreference!(copy_pfam_na, 1)
  gapstrip!(copy_pfam)
  gapstrip!(copy_pfam_na)
  residuefraction(copy_pfam[1,:]) == 1.0
  residuefraction(copy_pfam_na[1,:]) == 1.0
end

print("""
Test adjustreference!
""")

let copy_pfam = copy(pfam), copy_pfam_na = copy(pfam_na)
  adjustreference!(copy_pfam)
  adjustreference!(copy_pfam_na)
  residuefraction(copy_pfam[1,:]) == 1.0
  residuefraction(copy_pfam_na[1,:]) == 1.0
end

print("""

Printers
--------
""")

print("""
Test asciisequence
""")

@test asciisequence(small, 4) == "DAYCMD"
@test asciisequence(pfam, 1) == asciisequence(pfam, "C3N734_SULIY/1-95")

print("""
Test printfasta
""")

let io = IOBuffer()
  print(io, small, FASTA)
  @test takebuf_string(io) == fasta_string
end

print("""
Test printpfam
""")

let io = IOBuffer()
  print(io, pfam, Stockholm)
  @test parse(takebuf_string(io), Stockholm) == pfam
end

let io = IOBuffer()
  print(io, small_na, Stockholm)
  @test parse(takebuf_string(io), Stockholm) == small_na
end


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

JLLWrappers.@generate_main_file_header("SCS")
JLLWrappers.@generate_main_file("SCS", UUID("f4f2fc5b-1d94-523c-97ea-2ab488bedf4b"))
end  # module SCS_jll


text:
using Distributions: Dirichlet
using SpecialFunctions: lgamma, digamma, lbeta, beta
using Base.MathConstants: γ

struct BatchKumaraswamy{T}
    a::T
    b::T
end

function _u2kumaraswamysample(u, kuma::BatchKumaraswamy)
    _one = one(eltype(u))
    return (_one .- u.^(_one ./ kuma.b)).^(_one ./ kuma.a)
end

function _u2logkumaraswamysample(u, kuma::BatchKumaraswamy)
    _one = one(eltype(u))
    _eps = eps(eltype(u))
    return log1p.(-exp.(log.(u .+ _eps) ./ kuma.b)) ./ kuma.a
end

"""
    rand(kuma::BatchKumaraswamy) where

Sample from Kumaraswamy distribution.

NOTE: `k.a` and `k.b` are assumed to be in batch

Ref: https://arxiv.org/abs/1605.06197
"""
function rand(kuma::BatchKumaraswamy)
    u = randsimilar(kuma.a)
    x = _u2kumaraswamysample(u, kuma)
    return x
end

function rand(kuma::BatchKumaraswamy{T}, dims::Int...) where {T<:Real}
    @assert length(kuma.a) == 1 "`rand` for multiple samples only supports for univariate case"
    @assert length(kuma.b) == 1 "`rand` for multiple samples only supports for univariate case"
    u = length(dims) == 0 ? rand(T) : randarr(dims)
    x = _u2kumaraswamysample(u, kuma)
    return x
end

function logrand(kuma::BatchKumaraswamy)
    u = randsimilar(kuma.a)
    logx = _u2logkumaraswamysample(u, kuma)
    return logx
end

function logrand(kuma::BatchKumaraswamy{T}, dims::Int...) where {T<:Real}
    @assert length(kuma.a) == 1 "`rand` for multiple samples only supports for univariate case"
    @assert length(kuma.b) == 1 "`rand` for multiple samples only supports for univariate case"
    u = length(dims) == 0 ? rand(T) : randarr(dims)
    logx = _u2logkumaraswamysample(u, kuma)
    return logx
end

"""
    logpdf(kuma::BatchKumaraswamy, x)

Compute ``Kumaraswamy(x; a, b)``.

WARN: this function is not tested.
"""
function logpdf(kuma::BatchKumaraswamy, x)
    _one = one(eltype(x))
    _eps = eps(eltype(x))
    lp = log.(kuma.a) .+
         log.(kuma.b) .+
         (kuma.a .- _one) .* log.(x .+ _eps) .+
         (kuma.b .- _one) .* log.(_one + _eps .- x.^kuma.a)
    return lp
end

"""
    kldiv(d1::Dirichlet, d2::Dirichlet)

Compute ``KL(Dir(α)||Dir(β))`` where ``α = [α_1, \\dots, α_K]`` and ``β = [β_1, \\dots, β_K]``.

Ref: http://bariskurt.com/kullback-leibler-divergence-between-two-dirichlet-and-beta-distributions/
"""
function kldiv(d1::Dirichlet, d2::Dirichlet)
    α = d1.alpha
    β = d2.alpha

    α0 = sum(α)
    β0 = sum(β)

    kl = (lgamma.(α0) .- sum(lgamma.(α))) .-
         (lgamma.(β0) - sum(lgamma.(β))) .+
         sum((α .- β) .* (digamma.(α) .- digamma(α0)))
    return kl
end

struct BatchBeta{T}
    α::T
    β::T
end

function logpdf(bb::BatchBeta, x)
    _one = one(eltype(x))
    _eps = eps(eltype(x))
    lp = (bb.α .- _one) .* log.(x .+ _eps) .+
         (bb.β .- _one) .* log.(_one + _eps .- x) .-
         lbeta(bb.α, bb.β)
    return lp
end

"""
    kldiv(bb1::BatchBeta, bb2::BatchBeta)

Compute ``KL(Beta(α1, β1)||Beta(α2, β2))``.
"""
function kldiv(bb1::BatchBeta, bb2::BatchBeta)
    α1 = bb1.α
    β1 = bb1.β
    α2 = bb2.α
    β2 = bb2.β

    α0 = α1 .+ β1
    β0 = α2 .+ β2

    kl = (lgamma.(α0) .- (lgamma.(α1) .+ lgamma.(β1))) .-
         (lgamma.(β0) .- (lgamma.(α2) .+ lgamma.(β2))) .+
         (α1 - α2) .* (digamma.(α1) - digamma.(α0)) .+
         (β1 - β2) .* (digamma.(β1) - digamma.(α0))
    return kl
end

"""
    kldiv(kuma::BatchKumaraswamy, bb::BatchBeta; M=10)

Compute ``KL(Kumaraswamy(a, b)||Beta(α, β))``.

NOTE: only `a` and `b` are assumed to be in batch
"""
function kldiv(kuma::BatchKumaraswamy, bb::BatchBeta; M::Int=11)
    a = kuma.a
    b = kuma.b
    α = bb.α
    β = bb.β
    T = eltype(a)
    _one = one(eltype(a))

    @assert M > 0
    a_times_b = a .* b
    acc = _one ./ (_one .+ a_times_b) .* beta(_one ./ a, b)
    for m = 2:M
        acc = acc .+ _one ./ (T(m) .+ a_times_b) .* beta(T(m) ./ a, b)
    end

    kl = (a .- α) ./ a .* (-T(γ) .- digamma.(b) .- _one ./ b) .+ log.(a_times_b) .+
         lbeta(α, β) .- (b .- _one) ./ b .+ (β - _one) .* b .* acc
    return kl
end


text:
__precompile__()

module EvolvingCurves

using FastGaussQuadrature
using FixedSizeArrays
using PyPlot
using FEFunctions

##########
# Plotting

export update

function PyPlot.plot{T}(
    x::LinearFEFunction{Vec{2,T}},
    c::LinearFEFunction{T};
    size = 6,
    vmin = nothing,
    vmax = nothing
)
    scat = scatter(
        vec(x)[1:2:end],
        vec(x)[2:2:end];
        c = vec(c),
        s = size,
        edgecolors = "face",
        vmin = vmin,
        vmax = vmax
    )
    gca()[:set_aspect]("equal")
    return scat
end
function update{T}(
    scat,
    x::LinearFEFunction{Vec{2,T}},
    c::LinearFEFunction{T}
)
    scat[:set_offsets](reshape(vec(x), (2,length(x)))')
    scat[:set_color](scat[:to_rgba](vec(c)))
end
function PyPlot.plot{T}(
    x::FEFunction{Vec{2,T}},
    c::FEFunction{T};
    major_size = 10,
    minor_size = 0.5,
    nsamples = 10,
    vmin = nothing,
    vmax = nothing
)
    scat_minor = plot(sample(x,nsamples),sample(c,nsamples), size=minor_size, vmin=vmin, vmax=vmax)
    scat_major = plot(sample(x,1),sample(c,1), size=major_size, vmin=vmin, vmax=vmax)
    return nsamples,scat_major,scat_minor
end
function update{T}(
    plotobj,
    x::FEFunction{Vec{2,T}},
    c::FEFunction{T}
)
    nsamples,scat_major,scat_minor = plotobj
    update(scat_major, sample(x,1), sample(c,1))
    update(scat_minor, sample(x,nsamples), sample(c,nsamples))
end

function PyPlot.plot{T}(
    x::LinearFEFunction{Vec{2,T}};
    size = 6,
    color = "k"
)
    scat = scatter(
        vec(x)[1:2:end],
        vec(x)[2:2:end];
        s = size,
        c = color,
        edgecolors = "face"
    )
    gca()[:set_aspect]("equal")
    return scat
end
function update{T}(
    scat,
    x::LinearFEFunction{Vec{2,T}}
)
    scat[:set_offsets](reshape(vec(x), (2,length(x)))')
end
function PyPlot.plot{T}(
    x::FEFunction{Vec{2,T}},
    major_size = 10,
    minor_size = 0.5,
    color = "k",
    nsamples = 10
)
    scat_minor = plot(sample(x,nsamples), size=minor_size)
    scat_major = plot(sample(x,1), size=major_size)
    return nsamples,scat_major,scat_minor
end
function update{T}(
    plotobj,
    x::FEFunction{Vec{2,T}},
)
    nsamples,scat_major,scat_minor = plotobj
    update(scat_major, sample(x,1))
    update(scat_minor, sample(x,nsamples))
end

############################
# Geometry utility functions

export perp
perp{T}(v::Vec{2,T}) = Vec{2,T}(-v[2],v[1])

export integrate
function integrate(x,c, nq = 5)
    (P,W) = gausslegendre(nq)
    P = 0.5*(P+1); W /= 2

    I = 0
    for i = 1:length(x)
        for (p,w) in zip(P,W)
            I += w*c[i,p]*norm(x[i,p,Val{1}])
        end
    end
    return I
end

##################
# Comparing curves

using Optim
export project
function project{T}(v::Vec{2,T}, x::FEFunction{Vec{2,T}})
    imin = 0
    pmin = zero(T)
    dmin = convert(T,Inf)
    for i = 1:length(x)
        res = optimize(p -> norm(v - x[i,p]), 0,1, abs_tol=1e-6)
        if dmin > res.f_minimum
            imin = i
            pmin = res.minimum
            dmin = res.f_minimum
        end
    end
    return imin,pmin,dmin
end

export linf_errors
function linf_errors(x,c, xref,cref)
    xlinf = zero(scalartype(x))
    clinf = zero(scalartype(c))
    for iref = 1:length(xref)
        i,p,d = project(xref[iref],x)
        xlinf = max(xlinf,d)
        clinf = max(clinf,abs(c[i,p] - cref[iref]))
    end
    return xlinf, clinf
end

########
# Evolve

export evolve
function evolve(
    x,T;
    α = 0.1,
    ks = 1.0,
    kb = 1.0,
    conserve_volume = true,
    nq = 5
)
    for it = 2:length(T)
        dt = T[it] - T[it-1]

        Mx = @assemble_mat(x, nq, begin
            xp = x[i,p,Val{1}]
        end,begin
            Mat{2,2,Float64}((
                -perp(xp)*shape(x,s,p)*shape(x,ss,p),
                  α * xp *shape(x,s,p)*shape(x,ss,p)
            ))'
        end)
        Ax = @assemble_mat(x, nq, begin
            xp = x[i,p,Val{1}]
            xpp = x[i,p,Val{2}]
        end, begin
            dt*Mat{2,2,Float64}((
                (
                    # Mean curvature
                    -ks*perp(xp)/norm(xp)^2*shape(x,s,p,Val{2})*shape(x,ss,p,Val{0})
                    # Willmore
                    +kb*(
                        perp(xp)/norm(xp)^2*shape(x,s,p,Val{2})*(
                            + shape(x,ss,p,Val{2})/norm(xp)^2
                            - dot(xp,xpp)/norm(xp)^4*shape(x,ss,p,Val{1})
                        )
                        +
                        0.5*dot(perp(xp),xpp)^2/norm(xp)^8*perp(xp)*shape(x,s,p,Val{2})*shape(x,ss,p,Val{0})
                    )
                ),
                (xp/norm(xp)^2*shape(x,s,p,Val{2})*shape(x,ss,p,Val{0}))
            ))'
        end)

        # Volume conservation
        if conserve_volume
            ν = @assemble_vec(x,nq, begin
                xp = x[i,p,Val{1}]
            end,begin
                norm(xp)*Vec{2}(shape(x,s,p),0)
            end)
            I = @assemble_vec(x,nq, begin
                xp = x[i,p,Val{1}]
            end,begin
                -perp(xp)*shape(x,s,p)
            end)
            A = [
                Mx - Ax          sparsevec(ν)
                sparsevec(I)'    spzeros(1,1)
            ]
            b = [
                Mx*vec(x);
                dot(I,vec(x))
            ]
            xnew = func(x,(A\b)[1:end-1])
        else
            xnew = func(x,(Mx - Ax)\(Mx*vec(x)))
        end
        x = xnew
    end
    return x
end
function evolve(
    x,c,f,T;
    α = 0.1,
    ks = 1.0,
    kb = 1.0,
    kc = 1.0,
    conserve_volume = true,
    nq = 5
)
    Mc = @assemble_mat(c, nq, begin
        xp = x[i,p,Val{1}]
    end,begin
         shape(c,s,p)*shape(c,ss,p)*norm(xp)
    end)
    for it = 2:length(T)
        dt = T[it] - T[it-1]

        Mx = @assemble_mat(x, nq, begin
            xp = x[i,p,Val{1}]
        end,begin
            Mat{2,2,Float64}((
                -perp(xp)*shape(x,s,p)*shape(x,ss,p),
                  α * xp *shape(x,s,p)*shape(x,ss,p)
            ))'
        end)
        Ax = @assemble_mat(x, nq, begin
            xp = x[i,p,Val{1}]
            xpp = x[i,p,Val{2}]
        end, begin
            dt*Mat{2,2,Float64}((
                (
                    # Mean curvature
                    -ks*perp(xp)/norm(xp)^2*shape(x,s,p,Val{2})*shape(x,ss,p,Val{0})
                    # Willmore
                    +kb*(
                        perp(xp)/norm(xp)^2*shape(x,s,p,Val{2})*(
                            + shape(x,ss,p,Val{2})/norm(xp)^2
                            - dot(xp,xpp)/norm(xp)^4*shape(x,ss,p,Val{1})
                        )
                        +
                        0.5*dot(perp(xp),xpp)^2/norm(xp)^8*perp(xp)*shape(x,s,p,Val{2})*shape(x,ss,p,Val{0})
                    )
                ),
                (xp/norm(xp)^2*shape(x,s,p,Val{2})*shape(x,ss,p,Val{0}))
            ))'
        end)
        fx = @assemble_vec(x,nq, begin
            xp = x[i,p,Val{1}]
            fci = f(c[i,p])
        end,begin
            Vec{2}(dt*fci*shape(x,s,p)*norm(xp),0.0)
        end)

        # Volume conservation
        if conserve_volume
            ν = @assemble_vec(x,nq, begin
                xp = x[i,p,Val{1}]
            end,begin
                norm(xp)*Vec{2}(shape(x,s,p),0)
            end)
            I = @assemble_vec(x,nq, begin
                xp = x[i,p,Val{1}]
            end,begin
                -perp(xp)*shape(x,s,p)
            end)
            A = [
                Mx - Ax          sparsevec(ν)
                sparsevec(I)'    spzeros(1,1)
            ]
            b = [
                Mx*vec(x) + fx;
                dot(I,vec(x))
            ]
            xnew = func(x,(A\b)[1:end-1])
        else
            xnew = func(x,(Mx - Ax)\(Mx*vec(x) + fx))
        end

        Mcnew = @assemble_mat(c, nq, begin
            xnewp = xnew[i,p,Val{1}]
        end,begin
            shape(c,s,p)*shape(c,ss,p)*norm(xnewp)
        end)
        Ac = @assemble_mat(c, nq, begin
            xnewp = xnew[i,p,Val{1}]
            xnewi = xnew[i,p]
            xi = x[i,p]
        end,begin
            (
                -dt*kc*shape(c,s,p,Val{1})*shape(c,ss,p,Val{1})/norm(xnewp)
                -dot(xnewp/norm(xnewp), xnewi - xi)*shape(c,s,p)*shape(c,ss,p,Val{1})
            )
        end)
        c = func(c,(Mcnew - Ac)\(Mc*vec(c)))

        x = xnew
        Mc = Mcnew
    end
    return x,c
end

end

max_stars_repo_name: ettersi/CellMotility

max_issues_repo_name: ettersi/CellMotility

max_forks_repo_name: ettersi/CellMotility


text:
using SystemImageBuilder

buildimage()
resetimage()


text:
"""
A handle for an individual database in the DB environment.
"""
type DBI
    handle::Cuint
    name::String
    DBI(dbi::Cuint, name::String) = new(dbi, name)
end

"Check if database is open"
isopen(dbi::DBI) = dbi.handle != zero(Cuint)

"Open a database in the environment"
function open(txn::Transaction, dbname::String = ""; flags::Cuint=zero(Cuint))
    cdbname = length(dbname) > 0 ? bytestring(dbname) : convert(Cstring, Ptr{UInt8}(C_NULL))
    handle = Cuint[0]
    ret = ccall((:mdb_dbi_open, liblmdb), Cint,
                (Ptr{Void}, Cstring, Cuint, Ptr{Cuint}),
                 txn.handle, cdbname, flags, handle)
    (ret != 0) && throw(LMDBError(ret))
    return DBI(handle[1], dbname)
end

"Wrapper of DBI `open` for `do` construct"
function open(f::Function, txn::Transaction, dbname::String = ""; flags::Cuint=zero(Cuint))
    dbi = open(txn, dbname, flags=flags)
    tenv = env(txn)
    try
        f(dbi)
    finally
        close(tenv, dbi)
    end
end

"Close a database handle"
function close(env::Environment, dbi::DBI)
    if !isopen(env)
        warn("Environment is closed")
    end
    ccall((:mdb_dbi_close, liblmdb), Void, (Ptr{Void}, Cuint), env.handle, dbi.handle)
    dbi.handle = zero(Cuint)
    return
end

"Retrieve the DB flags for a database handle"
function flags(txn::Transaction, dbi::DBI)
    flags = Cuint[0]
    ret = ccall((:mdb_dbi_flags, liblmdb), Cint,
                (Ptr{Void}, Cuint, Ptr{Cuint}),
                 txn.handle, dbi.handle, flags)
    (ret != 0) && throw(LMDBError(ret))
    return flags[1]
end

"""Empty or delete+close a database.

If parameter `delete` is `false` DB will be emptied, otherwise
DB will be deleted from the environment and DB handle will be closed
"""
function drop(txn::Transaction, dbi::DBI; delete=false)
    del = delete ? int32(1) : int32(0)
    ret = ccall((:mdb_drop, liblmdb), Cint,
                (Ptr{Void}, Cuint, Cint),
                 txn.handle, dbi.handle, del)
    (ret != 0) && throw(LMDBError(ret))
    return ret
end

"Store items into a database"
function put!(txn::Transaction, dbi::DBI, key, val; flags::Cuint=zero(Cuint))
    mdb_key_ref = Ref(MDBValue(key))
    mdb_val_ref = Ref(MDBValue(val))

    ret = ccall((:mdb_put, liblmdb), Cint,
                (Ptr{Void}, Cuint, Ptr{MDBValue}, Ptr{MDBValue}, Cuint),
                txn.handle, dbi.handle, mdb_key_ref, mdb_val_ref, flags)

    (ret != 0) && throw(LMDBError(ret))
    return ret
end

"Delete items from a database"
function delete!(txn::Transaction, dbi::DBI, key, val)
    mdb_key_ref = Ref(MDBValue(key))
    mdb_val_ref = Ref(MDBValue(val))

    ret = ccall((:mdb_del, liblmdb), Cint,
                (Ptr{Void}, Cuint, Ptr{MDBValue}, Ptr{MDBValue}),
                txn.handle, dbi.handle, mdb_key_ref, mdb_val_ref)

    (ret != 0) && throw(LMDBError(ret))
    return ret
end

"Get items from a database"
function get{T}(txn::Transaction, dbi::DBI, key, ::Type{T})
    # Setup parameters
    mdb_key_ref = Ref(MDBValue(key))
    mdb_val_ref = Ref(MDBValue())

    # Get value
    ret = ccall((:mdb_get, liblmdb), Cint,
                 (Ptr{Void}, Cuint, Ptr{MDBValue}, Ptr{MDBValue}),
                 txn.handle, dbi.handle, mdb_key_ref, mdb_val_ref)
    (ret != 0) && throw(LMDBError(ret))

    # Convert to proper type
    mdb_val = mdb_val_ref[]
    if T <: AbstractString
        return bytestring(convert(Ptr{UInt8}, mdb_val.data), mdb_val.size)
    else
        nvals = floor(Int, mdb_val.size/sizeof(T))
        value = pointer_to_array(convert(Ptr{T}, mdb_val.data), nvals)
        return length(value) == 1 ? value[1] : value
    end
end


text:
module adapters_postgresql_window_frames_test

using Test # @test
using Octo.Adapters.PostgreSQL # to_sql window SELECT FROM RANK OVER PARTITION BY ORDER DESC

u = from([PARTITION BY :depname ORDER BY :salary DESC])
@test to_sql([SELECT (:depname, :empno, :salary, [RANK() OVER u]) FROM :empsalary]) ==
             "SELECT depname, empno, salary, RANK() OVER (PARTITION BY depname ORDER BY salary DESC) FROM empsalary"
@test to_sql([WINDOW :w AS u]) == "WINDOW w AS (PARTITION BY depname ORDER BY salary DESC)"

w = from([PARTITION BY :depname ORDER BY :salary DESC], :w)
@test to_sql([SELECT (:depname, :empno, :salary, [RANK() OVER w]) FROM :empsalary]) ==
             "SELECT depname, empno, salary, RANK() OVER (PARTITION BY depname ORDER BY salary DESC) AS w FROM empsalary"
@test to_sql([WINDOW :w AS w]) == "WINDOW w AS (PARTITION BY depname ORDER BY salary DESC)"
@test to_sql([w.salary]) == "w.salary"

Repo.debug_sql()

include("options.jl")

Repo.connect(;
    adapter = Octo.Adapters.PostgreSQL,
    Options.for_postgresql...
)

# https://robots.thoughtbot.com/postgres-window-functions
struct Post
end
Schema.model(Post, table_name="posts")

struct Comment
end
Schema.model(Comment, table_name="comments")

Repo.execute([DROP TABLE IF EXISTS Post])
Repo.execute([DROP TABLE IF EXISTS Comment])

Repo.execute(Raw("""
CREATE TABLE posts (
  id integer PRIMARY KEY,
  body varchar,
  created_at timestamp DEFAULT current_timestamp
);

CREATE TABLE comments (
 id INTEGER PRIMARY KEY,
 post_id integer NOT NULL,
 body varchar,
 created_at timestamp DEFAULT current_timestamp
);

/* make two posts */
INSERT INTO posts VALUES (1, 'foo');
INSERT INTO posts VALUES (2, 'bar');

/* make 4 comments for the first post */
INSERT INTO comments VALUES (1, 1, 'foo old');
INSERT INTO comments VALUES (2, 1, 'foo new');
INSERT INTO comments VALUES (3, 1, 'foo newer');
INSERT INTO comments VALUES (4, 1, 'foo newest');

/* make 4 comments for the second post */
INSERT INTO comments VALUES (5, 2, 'bar old');
INSERT INTO comments VALUES (6, 2, 'bar new');
INSERT INTO comments VALUES (7, 2, 'bar newer');
INSERT INTO comments VALUES (8, 2, 'bar newest');
"""))

posts = from(Post, :posts)
comments = from(Comment, :comments)

q = [SELECT (as(posts.id, :post_id), as(comments.id, :comment_ids), as(comments.body, :body)) FROM posts LEFT OUTER JOIN comments ON posts.id == comments.post_id]
@test to_sql(q) == "SELECT posts.id AS post_id, comments.id AS comment_ids, comments.body AS body FROM posts LEFT OUTER JOIN comments ON posts.id = comments.post_id"
df = Repo.query(q)
@test size(df) == (8,)

comment_rank = from([PARTITION BY :post_id ORDER BY comments.created_at DESC], :comment_rank)
@test to_sql([DENSE_RANK() OVER comment_rank]) == "DENSE_RANK() OVER (PARTITION BY post_id ORDER BY comments.created_at DESC) AS comment_rank"
q = [SELECT (as(posts.id, :post_id), as(comments.id, :comment_id), as(comments.body, :body), [DENSE_RANK() OVER comment_rank]) FROM posts LEFT OUTER JOIN comments ON posts.id == comments.post_id]
@test to_sql(q) == "SELECT posts.id AS post_id, comments.id AS comment_id, comments.body AS body, DENSE_RANK() OVER (PARTITION BY post_id ORDER BY comments.created_at DESC) AS comment_rank FROM posts LEFT OUTER JOIN comments ON posts.id = comments.post_id"
df = Repo.query(q)
@test size(df) == (8,)

comment_id = as(comments.id, :comment_id)
post_id = as(posts.id, :post_id)
body = as(comments.body, :body)
ranked_comments = from([SELECT (post_id, comment_id, body, [DENSE_RANK() OVER comment_rank]) FROM posts LEFT OUTER JOIN comments ON posts.id == comments.post_id], :ranked_comments)
@test to_sql(ranked_comments) == "(SELECT posts.id AS post_id, comments.id AS comment_id, comments.body AS body, DENSE_RANK() OVER (PARTITION BY post_id ORDER BY comments.created_at DESC) AS comment_rank FROM posts LEFT OUTER JOIN comments ON posts.id = comments.post_id) AS ranked_comments"
q = [SELECT (:comment_id, :post_id, :body) FROM ranked_comments WHERE comment_rank < 4]
@test to_sql(q) == "SELECT comment_id, post_id, body FROM (SELECT posts.id AS post_id, comments.id AS comment_id, comments.body AS body, DENSE_RANK() OVER (PARTITION BY post_id ORDER BY comments.created_at DESC) AS comment_rank FROM posts LEFT OUTER JOIN comments ON posts.id = comments.post_id) AS ranked_comments WHERE comment_rank < 4"
df = Repo.query(q)
@test size(df) == (8,)

ranked_comments = from([SELECT (post_id, comment_id, body, [DENSE_RANK() OVER comment_rank]) FROM posts LEFT OUTER JOIN comments ON posts.id == comments.post_id])
with = [WITH :ranked_comments AS ranked_comments]
@test to_sql(with) == "WITH ranked_comments AS (SELECT posts.id AS post_id, comments.id AS comment_id, comments.body AS body, DENSE_RANK() OVER (PARTITION BY post_id ORDER BY comments.created_at DESC) AS comment_rank FROM posts LEFT OUTER JOIN comments ON posts.id = comments.post_id)"

q = [with... SELECT (:post_id, :comment_id, :body) FROM :ranked_comments WHERE comment_rank < 4]
@test to_sql(q) == string(to_sql(with), ' ', "SELECT post_id, comment_id, body FROM ranked_comments WHERE comment_rank < 4")
df = Repo.query(q)
@test size(df) == (8,)

# https://stackoverflow.com/questions/12410791/difficult-for-me-postgres-sql-query
struct C
end
Schema.model(C, table_name="c")

struct D
end
Schema.model(D, table_name="d")

Repo.execute([DROP TABLE IF EXISTS C])
Repo.execute([DROP TABLE IF EXISTS D])
Repo.execute(Raw("""
CREATE TABLE c
  (
    id int NOT NULL PRIMARY KEY,
    aid int NOT NULL,
    bid int NOT NULL,
    units int NOT NULL
  );
CREATE TABLE d
  (
    id int NOT NULL PRIMARY KEY,
    open boolean NOT NULL,
    cid int NOT NULL
  );
"""))
Repo.insert!(C, [(2,3,5,4),(3,3,5,6),(4,4,6,8),(5,4,6,10),(6,7,8,9)])
Repo.insert!(D, [(1,true,2),(2,true,3),(3,true,3),(4,true,4)])

c = from(C, :c)
d = from(D, :d)
n = from([SELECT (:aid, as(COUNT(*), :cnt)) FROM c JOIN d ON (d.cid == c.id) GROUP BY :aid], :n)
@test to_sql(n)              ==           "(SELECT aid, COUNT(*) AS cnt FROM c JOIN d ON d.cid = c.id GROUP BY aid) AS n"
@test to_sql([WITH :n AS n]) == "WITH n AS (SELECT aid, COUNT(*) AS cnt FROM c JOIN d ON d.cid = c.id GROUP BY aid)"
@test to_sql([JOIN n])       == "JOIN n"
q = [WITH :n AS n SELECT (:aid, as(SUM(c.units) - COALESCE(n.cnt, 0), :difference)) FROM c LEFT JOIN n USING (:aid) GROUP BY (:aid, n.cnt) ORDER BY :aid]
@test to_sql(q)              == "WITH n AS (SELECT aid, COUNT(*) AS cnt FROM c JOIN d ON d.cid = c.id GROUP BY aid) SELECT aid, (SUM(c.units) - COALESCE(n.cnt, 0)) AS difference FROM c LEFT JOIN n USING (aid) GROUP BY aid, n.cnt ORDER BY aid"

df = Repo.query(q)
@test Pretty.table(df) == """
|   aid |   difference |
| ----- | ------------ |
|     3 |            7 |
|     4 |           17 |
|     7 |            9 |
3 rows."""

Repo.disconnect()

end # module adapters_postgresql_window_frames_test


text:
function multi_startswith(str::AbstractString, pats...)
    for pat in pats
        startswith(str, pat) && return true
    end
    false
end

function constantly(value)
    function(_...)
        value
    end
end

function tolist(a::AbstractArray)
  l = list()
  for item in a
    l = cons(item, l)
  end
  reverse(l)
end

date_formats = ("yyyy-mm-ddTHH:MM:SS.sss", "yyyy-mm-ddTHH:MM:SS")

function parsedatetime(s::AbstractString)
    s = replace(s, r"Z$", "")
    s = replace(s, r"-\d\d:\d\d$", "")
    for fmt in date_formats
        try
            return DateTime(s, fmt)
	catch ex
            println(ex)
	end
    end
    throw(ArgumentError("Dont know how to parse date/time: $s"))
end


text:
using Test
using Repos

include("mechanics.jl")
include("usecase.jl")

@testset "Internal Mechanics" begin
    swaptest()
    classtest()
end

@testset "Usecase Tests     " begin
    divisortest()
end


text:
#====================================================

    generate_SFC_models.jl

    generates definitions for models
    with same formula for each cycle (SFC models)

====================================================#

import DataStructures.OrderedDict
using JuMP

## function definitions

"Generate prediction function given SFC model definition."
function add_funcs_pred!(
    md          ::SFCModelDef,
    verbose     ::Bool = false
)
    _x_args_str = join(map(_x_str -> "$_x_str::Real", md._x_strs), ", ")
    coefs_str = join(map(str -> "$str::Real", md.coef_strs), ", ")
    #
    for func_key in MD_func_keys
        func_name = "$(md.name)_$func_key"
        func_str = join([
            "function $func_name($_x_args_str, $coefs_str)", ## In v0.4.5, `;` can't be parsed correctly for anonymous function lambda, because the signature will be parsed as a tuple not allowing `;`; default values are not allowed either. Heard it'll be fixed in v0.5
                md.pred_strs[func_key],
            "end",
            "function $func_name(X_cases::AbstractVector, $coefs_str)", ## X_cases is a vector (indexed by cases) of vectors (indexed by features)
                "map(X_case -> $func_name(X_case..., $coefs_str), X_cases)",
            "end"
        ], "; ")
        md.func_pred_strs[func_key] = func_str
        func_expr = Base.parse(func_str) ## not JSON.parse
        md.funcs_pred[func_key] = @eval $func_expr
    end
    #
    return nothing ## side effects only
end

"Generate fit function given SFC model definition."
function add_func_fit!( ## vco = variable constraints objective
    md          ::SFCModelDef;
    Y_str       ::String = "Y",
    obj_algrt   ::Symbol = :RSS,
    sense       ::Symbol = :Min, ## :Min, :Max
)
    ## X_strs
    X_strs = md.X_strs
    num_fts = length(X_strs)
    X_args_str = join(map(X_str -> "$X_str::AbstractVector", X_strs), ", ")
    #
    ## function signature
    sig_str = "function $(md.name)_func_fit(
        $X_args_str,
        $Y_str::AbstractVector,
        wts::AbstractVector=ones(length($Y_str));
        kwargs_Model...
    )"
    #
    ## initiate model
    mod_init_str = "jmp_model = Model(;kwargs_Model...)"
    #
    ## define coefficients as variables
    coef_init_str = "$X_args_str, $Y_str = map(abs_vec -> Array(abs_vec), ($X_args_str, $Y_str)); init_coefs = $(md.name)_func_init_coefs($X_args_str, $Y_str)" ## `Array` because `linreg` doesn't work on `DataArray`
    var_str = join(
        map(md.coef_strs) do coef_str
            "@variable(
                jmp_model,
                $coef_str,
                start=init_coefs[\"$coef_str\"])"
        end,
        "; ")
    #
    ## add constraints for the coefficents
    cnstrnts_str = join(
        map(md.coef_cnstrnts) do coef_constrnt
            "@constraint(jmp_model, $coef_constrnt)"
        end,
        "; ")
    #
    ## set objective
    obj_macro_str = md.linear ? "@objective" : "@NLobjective" ## `a1 = :(@some_macro); :($a1(arg1, arg2))` is equivalent to :((@some_macro()(arg1,arg2))), both of which raises "syntax: invalid macro use \"@($a1)\"". This why obj_macro need to be string, and other expressions are started as strings too for convenience.
    #
    func_str_replaced = md.pred_strs[:f]
    for j in 1:num_fts
        func_str_replaced = replace(
            func_str_replaced,
            md._x_strs[j],
            "$(X_strs[j])[i]")
    end
    #
    residual_str = "($func_str_replaced - $Y_str[i])"
    iter_str = "for i in eachindex($(X_strs[1]))" ## assuming X and Y have the same length
    #
    if obj_algrt == :RSS
        obj_expr_str = "sum(wts[i] * $residual_str^2 $iter_str)"
        # sumabs2(map(eval(x_symbol) -> eval(func_expr), X) .- Y)
    elseif obj_algrt == :l1_norm
        obj_expr_str = "sum(abs($residual_str) $iter_str)" ## `abs()` may cause "Ipopt finished with status Restoration_Faild"
    end
    #
    obj_str = "$obj_macro_str(jmp_model, $sense, $obj_expr_str)"
    #
    ## return
    return_str = join([
        "status = solve(jmp_model)",
        "coef_strs = [\"$(join(md.coef_strs, "\", \""))\"]",
        "coefs = map(getvalue, [$(join(md.coef_strs, ", "))])",
        "obj_val = getobjectivevalue(jmp_model)",
        "return SFCFit(coef_strs, coefs, status, obj_val, jmp_model, init_coefs); end"
    ], "; ")
    #
    ## add definition of func_fit
    func_str = join(
        [sig_str, mod_init_str, coef_init_str, var_str, cnstrnts_str, obj_str, return_str],
        "; ")
    md.func_fit_str = func_str
    #
    func_expr = Base.parse(func_str) ## not JSON.parse
    md.func_fit = @eval $func_expr
    #
    return nothing ## side effects only
end

## constants (MD_ = model definition)

## when `num_fts > 1`, "d*" are partial derivatives in vector of length `num_fts`
const MD_func_keys = [:f, :inv, :bl, :dr1, :dr2]

# function empty_func() end
# function empty_func(arg1::Any=0, args...; kwargs...) end
function empty_func(args...; kwargs...) end

## `EMPTY_fp` for `func_pred_strs` and `funcs_pred`
const EMPTY_fp = map(("", empty_func)) do empty_val
    # OrderedDict(map(MD_func_keys) do func_key ## v0.4, `supertype` not defined, `typeof(some_function) == Function`
    OrderedDict{Symbol,supertype(typeof(empty_val))}(map(MD_func_keys) do func_key ## v0.5, `super` becomes `supertype`, `typeof(some_function) == some_function AND supertype(typeof(some_function)) == Function`
        func_key => empty_val
    end) ## do func_key
end ## do empty_val

## generate generic model definition objects
const MD_EMPTY_vals = (
    EMPTY_fp...,    ## :func_pred_strs, :funcs_pred
    "",             ## func_fit_str
    empty_func)     ## func_fit
const SFC_MODEL_DEFS = OrderedDict(map(SFC_MODEL_BASES) do SFC_model_base
    SFC_model_base[1] => SFCModelDef(
        SFC_model_base...,
        deepcopy(MD_EMPTY_vals)...)
end)
for md_ in values(SFC_MODEL_DEFS)
    add_funcs_pred!(md_)
    add_func_fit!(md_)
end

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
# # SEIR
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/SEIR.ipynb)
#
#md # !!! note "Overview"
#md #     System type: Nonlinear system\
#md #     State dimension: 7\
#md #     Application domain: Epidemiology
#
# ## Model description
#
# The SEIR Model is an Compartmental model, these try to predict things such as
# how a disease spreads, or the total number infected, or the duration of an
# epidemic, and to estimate various epidemiological parameters such as the
# reproductive number.
# The dynamics are described as follows:
#
# ```math
#   \left\{ \begin{array}{lcl}
#   \frac{dS}{dt} & = & \beta I S \\
#   \frac{dE}{dt} & = & \beta I S - \alpha E \\
#   \frac{dI}{dt} & = & -\gamma I - \alpha E \\
#   \frac{dR}{dt} & = & \gamma I
#   \end{array} \right.
# ```
#
# where ``S`` is the stock of susceptible population, ``E`` is the stock of
# exposed population``I`` is the stock of infected, ``R`` is the stock of
# removed population (either by death or recovery), with ``S + E + I + R = N``

using ReachabilityAnalysis, Plots

@taylorize function seir2!(du,u,p,t)
  S, E, I, R, α, β, γ = u

  βIS = β * (I * S)
  αE = α*E
  γI = γ*I

  du[1] = -βIS      # dS
  du[2] = βIS - αE  # dE
  du[3] = -γI + αE  # dI
  du[4] = γI        # dR

  #uncertain parameters
  local zerou = zero(u[1])
  du[5] = zerou
  du[6] = zerou
  du[7] = zerou
end

E₀ = 1e-4
u₀ = [1-E₀, E₀, 0, 0]
param_error = 0.01
α = 0.2 ± param_error
β = 1.0 ± 0.0
γ = 0.5 ± param_error
p = [α, β, γ]
X0 = IntervalBox(vcat(u₀, p));
prob = @ivp(x' = seir2!(x), dim=7, x(0) ∈ X0);

# ## Reachability settings
#
# The initial values used were ``E₀ = 1e-4, u₀ = [1-E₀, E₀, 0, 0]``,
# ``α = 0.2 ± 0.01, β = 1.0 ± 0.0, γ = 0.5 ± 0.01`` and ``p = [α, β, γ]``, for a
# time span of 200.
# The algorithm used was `TMJets` with ``n_T=7`` and ``n_Q=1``

# ## Results

sol = solve(prob, tspan=(0.0, 200.0), alg=TMJets21a(orderT=7, orderQ=1));
solz = overapproximate(sol, Zonotope);

LazySets.set_ztol(Float64, 1e-13)
fig = plot(legend=:outerright)
plot!(fig, solz, vars=(0, 1), color=:blue, lw=0.0, lab="S")
plot!(fig, solz, vars=(0, 2), color=:green, lw=0.0, lab="E")
plot!(fig, solz, vars=(0, 3), color=:red, lw=0.0, lab="I")
plot!(solz, vars=(0, 4), color=:grey, lw=0.0, lab="R")
fig