text:
using GSL

function expectation(M_B::Int16,u_B::Int16,p_m::Float64,p_h::Float64,betas::Vector{Float64})
	  log_beta = 0
    for k=0:M_B
        h_B = M_B - k
        m_B = k
    	  log_beta = sf_lnchoose(M_B, k) + m_B * log( p_m ) +  h_B * log ( p_h ) -  M_B * log( p_m + p_h )
        #println("k: $k, log_beta: $log_beta, p_m $p_m, p_h $p_h")
        betas[k+1] = exp(log_beta)
    end
    return betas
end

function maximisation(M_B::Int16, u_B::Int16, p_m::Float64, betas::Vector{Float64})
    (num,denom) = (0,0)
    #h_B = int16(0)
    #k=int16(0)
    for k=0:M_B
    	h_B = M_B - k
        #m_B = k
        num   +=  betas[k+1] * h_B * ( 1 - p_m)
        denom +=  betas[k+1] * ( u_B + h_B )
    end
    return num/denom
end

function update_pm(M_B::Int16, u_B::Int16, m_O::Int16, U_O::Int16, betas::Vector{Float64})
    num = m_O
    denom=0
    for k=0:M_B
    	m_B=k
    	num += betas[k+1]*m_B
    end
    denom = M_B + m_O + u_B + U_O
    return num/denom
end

 # Go for the estimate with the highest coverage??
function start_point(M_B,u_B,m_O,U_O)
        starting_estimate_p_m = M_B / (M_B + u_B)
        starting_estimate_p_h = 1 - starting_estimate_p_m
        return ( starting_estimate_p_m, starting_estimate_p_h )
end

# function returns the values
# p_m and p_h and the number of iterations taken to find the solution
function EM( M_B::Int16, u_B::Int16, m_O::Int16, U_O::Int16, p_m::Float64, p_h::Float64)
	tolerance = float64(1e-15)
	iter = int(0);
	verbose = 0
	betas = Array(Float64,M_B+1)
	while true
		betas = expectation(M_B, u_B, p_m, p_h,betas)
		p_m_old = float64(p_m)
    p_h_old = float64(p_h)
    if verbose == 1
       println("p_m: $p_m_old p_h: $p_h_old")
		end
		p_h   = maximisation(M_B, u_B, p_m, betas)
		p_m   = update_pm(M_B, u_B, m_O, U_O, betas)
		delta = max(abs(p_m_old - p_m), abs(p_h_old - p_h));
	      if iter > 500
            tolerance=delta
            if verbose == 1
                 println("Reached 500 iterations ( delta = $delta )")
            end
            return p_m,p_h, iter
        end

        if (delta <= tolerance)
        	if verbose == 1
                println( "delta <= tolerance ( $delta <= $tolerance ) - reached after $iter iterations ")
            end
            return p_m,p_h, iter
            #break
        end
        iter +=1
	end
end

#=
function EM(BS,OX,count)

   if count % 10000 == 0
           t = TmStruct( time() )
           t_str = join( [t.hour,t.min,t.sec], ":")
           println("On $count. current time: ", t_str)
   end

   if any( x -> x < 0, (  BS.T_count, BS.C_count, OX.C_count, OX.T_count )  )
      println( "error negative counts: ", join( [BS.pos.chr,BS.pos.start, BS.T_count, BS.C_count, OX.C_count, OX.T_count], "\t" ) )
      return CytosineMethylation( BS.pos, BS.strand, 2, 2 )
   end

   if ! ( BS.pos.chr == OX.pos.chr && BS.pos.start == OX.pos.start )
     println("Error: BS and OX positions do not match up:  ", join( [BS.pos.chr,BS.pos.start, BS.T_count, BS.C_count, OX.pos.chr, OX.pos.start, OX.C_count, OX.T_count], "\t" ) )
   end

   est_p_m = ( OX.C_count / (OX.C_count + OX.T_count ))
   est_p_h = ( BS.C_count / ( BS.T_count + BS.C_count ) ) -  est_p_m

   if( est_p_h >= 0  )
        return CytosineMethylation( BS.pos, BS.strand, est_p_m, est_p_h )
   else
       (starting_estimate_p_m, starting_estimate_p_h) = start_point(BS.C_count,BS.T_count,OX.C_count,OX.T_count)
       (em_p_m,em_p_h) = EM(BS.C_count,BS.T_count,OX.C_count,OX.T_count,starting_estimate_p_m, starting_estimate_p_h)
       return CytosineMethylation( BS.pos, BS.strand, em_p_m, em_p_h )
    end
end

function run_mlml(BS,OX)
   println("Start run_mlml")
   return pmap( EM, BS, OX, 1:length(BS) )
end

function EM( BS::Vector{CytosineCount}, OX::Vector{CytosineCount}, count )
    if count % 10000 == 0
           t = TmStruct( time() )
           t_str = join( [t.hour,t.min,t.sec], ":")
           println("On $count. current time: ", t_str)
   end
   if any( x -> x < 0, (  BS.T_count, BS.C_count, OX.C_count, OX.T_count )  )
      println( "error negative counts: ", join( [BS.pos.chr,BS.pos.start, BS.T_count, BS.C_count, OX.C_count, OX.T_count], "\t" ) )
      return CytosineMethylation( BS.pos, BS.strand, 2, 2 )
   end

   if ! ( BS.pos.start == OX.pos.start )
     println("Error: BS and OX positions do not match up:  ", join( [BS.pos.start, BS.T_count, BS.C_count, OX.pos.start, OX.C_count, OX.T_count], "\t" ) )
   end

   est_p_m = ( OX.C_count / (OX.C_count + OX.T_count ))
   est_p_h = ( BS.C_count / ( BS.T_count + BS.C_count ) ) -  est_p_m

   if( est_p_h >= 0  )
        return CytosineMethylation( BS.pos, BS.strand, est_p_m, est_p_h )
   else
       (starting_estimate_p_m, starting_estimate_p_h) = start_point(BS.C_count,BS.T_count,OX.C_count,OX.T_count)
       (em_p_m,em_p_h) = EM(BS.C_count,BS.T_count,OX.C_count,OX.T_count,starting_estimate_p_m, starting_estimate_p_h)
       return CytosineMethylation( BS.pos, BS.strand, em_p_m, em_p_h )
    end
end

function run_mlml(BS::SequenceFeature, OX::SequenceFeature)
    # get list of chr names
    chrs = keys(BS.features)
    MethRatios = SequenceFeatures(Dict(),Dict())
    # go through each chromosome sharing processing
    for chr in chrs
       BS_cytosine_counts = BS.features[chr]
       OX_cytosine_counts = OX.features[chr]
       cm = pmap(EM, BS_cytosine_counts, OX_cytosine_counts, 1:length(BS_cytosine_counts) )
       MethRatios.feature[chr] = cm
    end
    return MethRatios
end
=#


text:
@testset "ExpLogSinCosTan" begin
    @testset for T in (Posit8, Posit16, Posit32,
                    # Posit8_1, Posit16_1, Posit24_1,
                    Posit8_2, Posit16_2, Posit24_2)

        @test one(T) == sin(T(π/2))
        @test minusone(T) == cos(T(1π))

        if T != Posit8
            @test one(T) == tan(T(π/4))
        else # due to rounding errors for Posit8 this is slightly off
            @test one(T) == nextfloat(tan(T(π/4)))
        end

        @test zero(T) == tan(zero(T))

        x = T(0.5)  # doesn't work for all x due to rounding errors
        # Pythagorean Identity
        @test one(T) == sin(x)^2 + cos(x)^2

        # Parity
        @test -sin(x) == sin(-x)
        @test cos(x) == cos(-x)
        @test -tan(x) == tan(-x)

        # Euler's Identity
        @test one(T) == abs(exp(im*T(1π)))

        @test one(T) == exp(zero(T))
        @test T(16.0) == exp2(T(4.0))
        @test T(10.0) == exp10(T(1.0))

        @test zero(T) == log(one(T))
        @test T(4.0) == log2(T(16.0))
        @test one(T) == log10(T(10.0))
    end
end

ext: jl


text:
# ---
# title: 521. Longest Uncommon Subsequence I
# id: problem521
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String, Brainteaser
# link: <https://leetcode.com/problems/longest-uncommon-subsequence-i/description/>
# hidden: true
# ---
#
# Given two strings `a` and `b`, find the length of the  **longest uncommon
# subsequence**  between them.
#
# A  **subsequence**  of a string `s` is a string that can be obtained after
# deleting any number of characters from `s`. For example, `"abc"` is a
# subsequence of `"aebdc"` because you can delete the underlined characters in
# `"a _e_ b _d_ c"` to get `"abc"`. Other subsequences of `"aebdc"` include
# `"aebdc"`, `"aeb"`, and `""` (empty string).
#
# An  **uncommon subsequence**  between two strings is a string that is a
# **subsequence of one  but not the other**.
#
# Return _the length of the **longest uncommon subsequence**  between `a` and
# `b`_. If the longest uncommon subsequence doesn't exist, return `-1`.
#
#
#
# **Example 1:**
#
#
#
#     Input: a = "aba", b = "cdc"
#     Output: 3
#     Explanation: One longest uncommon subsequence is "aba" because "aba" is a subsequence of "aba" but not "cdc".
#     Note that "cdc" is also a longest uncommon subsequence.
#
#
# **Example 2:**
#
#
#
#     Input: a = "aaa", b = "bbb"
#     Output: 3
#     Explanation:  The longest uncommon subsequences are "aaa" and "bbb".
#
#
# **Example 3:**
#
#
#
#     Input: a = "aaa", b = "aaa"
#     Output: -1
#     Explanation:  Every subsequence of string a is also a subsequence of string b. Similarly, every subsequence of string b is also a subsequence of string a.
#
#
#
#
# **Constraints:**
#
#   * `1 <= a.length, b.length <= 100`
#   * `a` and `b` consist of lower-case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
const MONTHS_PER_YEAR = Int64(12)
const DAYS_PER_WEEK = Int64(7)
const HOURS_PER_DAY = Int64(24)
const HOURS_PER_WEEK = HOURS_PER_DAY * DAYS_PER_WEEK
const MINUTES_PER_HOUR = Int64(60)
const MINUTES_PER_DAY = MINUTES_PER_HOUR * HOURS_PER_DAY
const MINUTES_PER_WEEK = MINUTES_PER_DAY * DAYS_PER_WEEK
const SECONDS_PER_MINUTE = Int64(60)
const SECONDS_PER_HOUR = SECONDS_PER_MINUTE * MINUTES_PER_HOUR
const SECONDS_PER_DAY = SECONDS_PER_HOUR * HOURS_PER_DAY
const SECONDS_PER_WEEK = SECONDS_PER_DAY * DAYS_PER_WEEK
const MILLISECONDS_PER_SECOND = Int64(1_000)
const MICROSECONDS_PER_MILLISECOND = Int64(1_000)
const NANOSECONDS_PER_MICROSECOND  = Int64(1_000)
const MICROSECONDS_PER_SECOND = Int64(1_000_000)
const NANOSECONDS_PER_MILLISECOND  = Int64(1_000_000)
const NANOSECONDS_PER_SECOND  = Int64(1_000_000_000)
const MILLISECONDS_PER_DAY = MILLISECONDS_PER_SECOND * SECONDS_PER_DAY
const MICROSECONDS_PER_DAY = MICROSECONDS_PER_SECOND * SECONDS_PER_DAY
const NANOSECONDS_PER_DAY = NANOSECONDS_PER_SECOND * SECONDS_PER_DAY
const MILLISECONDS_PER_WEEK = MILLISECONDS_PER_DAY * DAYS_PER_WEEK
const MICROSECONDS_PER_WEEK = MICROSECONDS_PER_DAY * DAYS_PER_WEEK
const NANOSECONDS_PER_WEEK = NANOSECONDS_PER_DAY * DAYS_PER_WEEK

Day(x::Week) = Day(x.value * DAYS_PER_WEEK)
Hour(x::Week) = Hour(x.value * HOURS_PER_WEEK)
Minute(x::Week) = Minute(x.value * MINUTES_PER_WEEK)
Second(x::Week) = Second(x.value * SECONDS_PER_WEEK)
Millisecond(x::Week) = Millisecond(x.value * MILLISECONDS_PER_WEEK)
Microsecond(x::Week) = Microsecond(x.value * MICROSECONDS_PER_WEEK)
Nanosecond(x::Week) = Nanosecond(x.value * NANOSECONDS_PER_WEEK)

Hour(x::Day) = Hour(x.value * HOURS_PER_DAY)
Minute(x::Day) = Minute(x.value * MINUTES_PER_DAY)
Second(x::Day) = Second(x.value * SECONDS_PER_DAY)
Millisecond(x::Day) = Millisecond(x.value * MILLISECONDS_PER_DAY)
Microsecond(x::Day) = Microsecond(x.value * MICROSECONDS_PER_DAY)
Nanosecond(x::Day) = Nanosecond(x.value * NANOSECONDS_PER_DAY)

Minute(x::Hour) = Minute(x.value * MINUTES_PER_HOUR)
Second(x::Hour) = Second(x.value * SECONDS_PER_HOUR)
Millisecond(x::Hour) = Millisecond(x.value * (MILLISECONDS_PER_SECOND * SECONDS_PER_HOUR))
Microsecond(x::Hour) = Microsecond(x.value * (MICROSECONDS_PER_SECOND * SECONDS_PER_HOUR))
Nanosecond(x::Hour) = Nanosecond(x.value * (NANOSECONDS_PER_SECOND * SECONDS_PER_HOUR))

Second(x::Minute) = Second(x.value * SECONDS_PER_MINUTE)
Millisecond(x::Minute) = Millisecond(x.value * (MILLISECONDS_PER_SECOND * SECONDS_PER_MINUTE))
Microsecond(x::Minute) = Microsecond(x.value * (MICROSECONDS_PER_SECOND * SECONDS_PER_MINUTE))
Nanosecond(x::Minute) = Nanosecond(x.value * (NANOSECONDS_PER_SECOND * SECONDS_PER_MINUTE))

Millisecond(x::Second) = Millisecond(x.value * MILLISECONDS_PER_SECOND)
Microsecond(x::Second) = Microsecond(x.value * MICROSECONDS_PER_SECOND)
Nanosecond(x::Second) = Nanosecond(x.value * NANOSECONDS_PER_SECOND)

Microsecond(x::Millisecond) = Microsecond(x.value * MICROSECONDS_PER_MILLISECOND)
Nanosecond(x::Millisecond) = Nanosecond(x.value * NANOSECONDS_PER_MILLISECOND)

Nanosecond(x::Microsecond) = Nanosecond(x.value * NANOSECONDS_PER_MICROSECOND)

#=
for P in (:Hour, :Minute, :Second, :Millisecond, :Microsecond, :Nanosecond)
  @eval Day(x::$P) = Day(canonical(x))
end
for P in (:Minute, :Second, :Millisecond, :Microsecond, :Nanosecond)
  @eval Hour(x::$P) = Hour(canonical(x))
end
=#


text:
#!/usr/bin/env julia6
using ReadGlobal, Decomp, StaticArrays

function TensorDirec_part!(t11::Vector{T},t22::Vector{T},t33::Vector{T},t12::Vector{T},t13::Vector{T},t23::Vector{T},dir::Vector{T}) where T<:AbstractFloat
    eigg  = [MVector{3,T}() for i in 1:Threads.nthreads()]
    eigv1 = [MVector{3,T}() for i in 1:Threads.nthreads()]
    eigv2 = [MVector{3,T}() for i in 1:Threads.nthreads()]
    eigv3 = [MVector{3,T}() for i in 1:Threads.nthreads()]
    Threads.@threads for i in 1:length(t11)
        j = Threads.threadid()
        @inbounds eigen!(t11[i],t22[i],t33[i],t12[i],t13[i],t23[i],eigg[j],eigv1[j],eigv2[j],eigv3[j])
        @inbounds dir[i] = abs(eigv1[j][3])
    end
end

function TensorDirec(t11::String,t22::String,t33::String,t12::String,t13::String,t23::String,out::String)
    doinchunks(TensorDirec_part!,input=(t11,t22,t33,t12,t13,t23),output=(out,))
    return 0
end

TensorDirec(ARGS[1],ARGS[4],ARGS[6],ARGS[2],ARGS[3],ARGS[5],ARGS[7])


text:
######################################################################
# NOTE THIS IS NOT OBSOLETE AND NEEDS TO BE REWRITTEN
######################################################################

error("THIS IS NOT OBSOLETE AND NEEDS TO BE REWRITTEN")

######################################################################
# count the AM column status for AMDB by year. large runtime, save
# results at the end, processed by another file (count_AM_tabulate.jl)
######################################################################

using AMDB
using TranscodingStreams
using CodecZlib
using JLD
using DocStringExtensions

"""
    $(SIGNATURES)

Read lines from `io`, split them on `delim`, and call `f` on the resulting vector of strings.

`limit` is passed to `split`, and sets the maximum number of fields.

Output a `.` every `progress` lines when the latter is positive.

When `maxlines` > 0, stop after that many lines.
"""
function dofields(f, io;
                  progress = 1000000, limit = 0, delim = ';', maxlines = 0)
    line = 0
    while !eof(io)
        if isa(progress, Integer) && line % progress == 0
            print(".")
        end
        line += 1
        if maxlines > 0 && line ≥ maxlines
            break
        end
        f(split(readline(io), delim; limit = limit))
    end
end

"""
    $(SIGNATURES)

Open `filename`, and count the unique strings in column `col_index`.

`options` are passed to `dofields` as keyword parameters.
"""
function count_field(filename, col_index; options...)
    io = open(filename, "r")
    c = Dict{String,Int}()
    dofields(GzipDecompressionStream(io);
             limit = col_index + 1, delim = ';', options...) do fields
                 kind = fields[col_index]
                 c[kind] = get(c, kind, 0) + 1
             end
    close(io)
    c
end

"""
    $(SIGNATURES)

Count all unique occurrences of strings in the given column `colname` for each year. Return a dictionary with elements `year => counts`, where counts is a `string => count` pair.
"""
function count_field_by_year(colname)
    counts = Dict{Int,Dict{String,Int}}()
    colnames = AMDB.colnames()
    colindex = findfirst(colnames, colname)
    @assert colindex > 0 "column $(colname) not found"
    for year in AMDB.all_years()
        println("processing $(year)")
        counts[year] = count_field(AMDB.data_file(year), colindex)
    end
    counts
end

status_counts = count_field_by_year("AM")

save(joinpath(AMDB.files_directory(), "status_counts.jld"), "status_counts",
     status_counts)


text:
## LEARNING NETWORK MACHINES

surrogate(::Type{<:Deterministic})  = Deterministic()
surrogate(::Type{<:Probabilistic})  = Probabilistic()
surrogate(::Type{<:Unsupervised}) = Unsupervised()
surrogate(::Type{<:Static}) = Static()

"""
    model_supertype(signature)

Return, if this can be deduced, which of `Deterministic`,
`Probabilistic` and `Unsupervised` is the appropriate supertype for a
composite model obtained by exporting a learning network with the
specified `signature`.

A learning network *signature* is a named tuple, such as
`(predict=yhat, transform=W)`, specifying what nodes of the network
are called to produce output of each operation represented by the
keys, in an exported version of the network.

If a supertype cannot be deduced, `nothing` is returned.

If the network with given `signature` is not exportable, this method
will not error but it will not a give meaningful return value either.

**Private method.**

"""
function model_supertype(signature)

    operations = keys(signature)

    length(intersect(operations, (:predict_mean, :predict_median))) == 1 &&
        return Deterministic

    if :predict in operations
        node = signature.predict
        if node isa Source
            return Deterministic
        end
        if node.machine !== nothing
            model = node.machine.model
            model isa Deterministic && return Deterministic
            model isa Probabilistic && return Probabilistic
        end
    end

    return nothing

end

caches_data_by_default(::Type{<:Surrogate}) = false

const ERR_MUST_PREDICT = ArgumentError(
    "You must specify at least `predict=<some node>`. ")
const ERR_MUST_TRANSFORM = ArgumentError(
    "You must specify at least `transform=<some node>`. ")
const ERR_MUST_OPERATE = ArgumentError(
    "You must specify at least one operation, as in `predict=<some node>`. ")
const ERR_MUST_SPECIFY_SOURCES = ArgumentError(
    "You must specify at least one source `Xs`, as in "*
    "`machine(surrogate_model, Xs, ...; kwargs...)`. ")

function check_surrogate_machine(::Surrogate, signature, _sources)
    isempty(signature) && throw(ERR_MUST_OPERATE)
    isempty(_sources) && throw(ERR_MUST_SPECIFY_SOURCES)
    return nothing
end

function check_surrogate_machine(::Union{Supervised,SupervisedAnnotator},
                                 signature,
                                 _sources)
    isempty(signature) && throw(ERR_MUST_PREDICT)
    length(_sources) > 1 || throw(err_supervised_nargs())
    return nothing
end

function check_surrogate_machine(::Union{Unsupervised},
                                 signature,
                                 _sources)
    isempty(signature) && throw(ERR_MUST_TRANSFORM)
    length(_sources) < 2 || throw(err_unsupervised_nargs())
    return nothing
end

function machine(model::Surrogate, _sources::Source...; pair_itr...)

    # named tuple, such as `(predict=yhat, transform=W)`:
    signature = (; pair_itr...)
    for op in keys(signature)
        op in OPERATIONS || throw(ArgumentError(
            "`$op` is not an admissible operation. "))
    end

    check_surrogate_machine(model, signature, _sources)

    mach = Machine(model, _sources...)

    mach.fitresult = signature

    return mach

end

function machine(_sources::Source...; pair_itr...)

    signature = (; pair_itr...)

    isempty(signature) && throw(ERR_MUST_OPERATE)

    T = model_supertype(signature)
    if T == nothing
        @warn "Unable to infer surrogate model type. \n"*
            "Using Deterministic(). To override specify "*
            "surrogate model, as in "*
        "`machine(Probabilistic(), ...)` or `machine(Interval(), ...)`"
        model = Deterministic()
    else
        model = surrogate(T)
    end

    return machine(model, _sources...; pair_itr...)

end

"""
    N = glb(mach::Machine{<:Surrogate})

A greatest lower bound for the nodes appearing in the signature of
`mach`.

**Private method.**

"""
glb(mach::Machine{<:Union{Composite,Surrogate}}) =
    glb(values(mach.fitresult)...)

"""
    fit!(mach::Machine{<:Surrogate};
         rows=nothing,
         acceleration=CPU1(),
         verbosity=1,
         force=false))

Train the complete learning network wrapped by the machine
`mach`.

More precisely, if `s` is the learning network signature used to
construct `mach`, then call `fit!(N)`, where `N = glb(values(s)...)`
is a greatest lower bound on the nodes appearing in the signature. For
example, if `s = (predict=yhat, transform=W)`, then call
`fit!(glb(yhat, W))`. Here `glb` is `tuple` overloaded for nodes.

See also [`machine`](@ref)

"""
function fit!(mach::Machine{<:Surrogate}; kwargs...)

    glb_node = glb(mach)
    fit!(glb_node; kwargs...)

    mach.state += 1
    mach.report = report(glb_node)
    return mach

end

MLJModelInterface.fitted_params(mach::Machine{<:Surrogate}) =
    fitted_params(glb(mach))

## CONSTRUCTING THE RETURN VALUE FOR A COMPOSITE FIT METHOD

# Identify which properties of `model` have, as values, a model in the
# learning network wrapped by `mach`, and check that no two such
# properties have have identical values (#377). Return the property name
# associated with each model in the network (in the order appearing in
# `models(glb(mach))`) using `nothing` when the model is not
# associated with any property.
function network_model_names(model::M,
                             mach::Machine{<:Surrogate}) where M<:Model

    signature = mach.fitresult
    network_model_ids = objectid.(MLJBase.models(glb(mach)))

    names = propertynames(model)

    # intialize dict to detect duplicity a la #377:
    name_given_id = Dict{UInt64,Vector{Symbol}}()

    # identify location of properties whose values are models in the
    # learning network, and build name_given_id:
    for name in names
        id = objectid(getproperty(model, name))
        if id in network_model_ids
            if haskey(name_given_id, id)
                push!(name_given_id[id], name)
            else
                name_given_id[id] = [name,]
            end
        end
    end

    # perform #377 check:
    no_duplicates = all(values(name_given_id)) do name
        length(name) == 1
    end
    if !no_duplicates
        for (id, name) in name_given_id
            if length(name) > 1
                @error "The hyperparameters $name of "*
                    "$model have identical model "*
                "instances as values. "
            end
        end
        throw(ArgumentError(
        "Two distinct hyper-parameters of a "*
            "composite model that are both "*
            "associated with models in the underlying learning "*
            "network (eg, any two components of a `@pipeline` model) "*
            "cannot have identical values, although they can be `==` "*
            "(corresponding nested properties are `==`). "*
            "Consider constructing instances "*
            "separately or use `deepcopy`. "))
    end

    return map(network_model_ids) do id
        if id in keys(name_given_id)
            return name_given_id[id] |> first
        else
            return nothing
        end
    end

end

"""

    return!(mach::Machine{<:Surrogate}, model, verbosity)

The last call in custom code defining the `MLJBase.fit` method for a
new composite model type. Here `model` is the instance of the new type
appearing in the `MLJBase.fit` signature, while `mach` is a learning
network machine constructed using `model`. Not relevant when defining
composite models using `@pipeline` or `@from_network`.

For usage, see the example given below. Specificlly, the call does the
following:

- Determines which hyper-parameters of `model` point to model
  instances in the learning network wrapped by `mach`, for recording
  in an object called `cache`, for passing onto the MLJ logic that
  handles smart updating (namely, an `MLJBase.update` fallback for
  composite models).

- Calls `fit!(mach, verbosity=verbosity)`.

- Moves any data in source nodes of the learning network into `cache`
  (for data-anonymization purposes).

- Records a copy of `model` in `cache`.

- Returns `cache` and outcomes of training in an appropriate form
  (specifically, `(mach.fitresult, cache, mach.report)`; see [Adding
  Models for General
  Use](https://alan-turing-institute.github.io/MLJ.jl/dev/adding_models_for_general_use/)
  for technical details.)

### Example

The following code defines, "by hand", a new model type `MyComposite`
for composing standardization (whitening) with a deterministic
regressor:

```
mutable struct MyComposite <: DeterministicComposite
    regressor
end

function MLJBase.fit(model::MyComposite, verbosity, X, y)
    Xs = source(X)
    ys = source(y)

    mach1 = machine(Standardizer(), Xs)
    Xwhite = transform(mach1, Xs)

    mach2 = machine(model.regressor, Xwhite, ys)
    yhat = predict(mach2, Xwhite)

    mach = machine(Deterministic(), Xs, ys; predict=yhat)
    return!(mach, model, verbosity)
end
```

"""
function return!(mach::Machine{<:Surrogate},
                 model::Union{Model,Nothing},
                 verbosity)

    _network_model_names = network_model_names(model, mach)

    verbosity isa Nothing || fit!(mach, verbosity=verbosity)

    # anonymize the data:
    sources = mach.args
    data = Tuple(s.data for s in sources)
    [MLJBase.rebind!(s, nothing) for s in sources]

    # record the current hyper-parameter values:
    old_model = deepcopy(model)

    cache = (sources = sources,
             data=data,
             network_model_names=_network_model_names,
             old_model=old_model)

    return mach.fitresult, cache, mach.report

end

#legacy code:
function (mach::Machine{<:Surrogate})()
    Base.depwarn("Calling a learning network machine `mach` "*
                 "with no arguments, as in"*
                 "`mach()`, is "*
                 "deprecated and could lead "*
                 "to unexpected behaviour for `Composite` models "*
                 "with hyper-parameters that are not models. "*
                 "Instead of `fit!(mach, verbosity=verbosity); return mach()` "*
                 "use `return!(mach, model, verbosity)`, "*
                 "where `model` is the `Model` instance appearing in your "*
                 "`MLJBase.fit` signature. Query the `return!` doc-string "*
                 "for details. ",
                 nothing)

    return!(mach, nothing, nothing)
end
network_model_names(model::Nothing, mach::Machine{<:Surrogate}) =
    nothing

## DUPLICATING AND REPLACING PARTS OF A LEARNING NETWORK MACHINE

"""
    replace(mach, a1=>b1, a2=>b2, ...; empty_unspecified_sources=false)

Create a deep copy of a learning network machine `mach` but replacing
any specified sources and models `a1, a2, ...` of the original
underlying network with `b1, b2, ...`.

If `empty_unspecified_sources=true` then any source nodes not
specified are replaced with empty source nodes.

"""
function Base.replace(mach::Machine{<:Surrogate},
                      pairs::Pair...; empty_unspecified_sources=false)

    signature = mach.fitresult
    interface_nodes = values(signature)

    W = glb(interface_nodes...)

    # Note: We construct nodes of the new network as values of a
    # dictionary keyed on the nodes of the old network. Additionally,
    # there are dictionaries of models keyed on old models and
    # machines keyed on old machines. The node and machine
    # dictionaries must be built simultaneously.

    # build model dict:
    model_pairs = filter(collect(pairs)) do pair
        first(pair) isa Model
    end
    models_ = models(W)
    models_to_copy = setdiff(models_, first.(model_pairs))
    model_copy_pairs = [model=>deepcopy(model) for model in models_to_copy]
    newmodel_given_old = IdDict(vcat(model_pairs, model_copy_pairs))

    # build complete source replacement pairs:
    sources_ = mach.args # sources(W)
    specified_source_pairs = filter(collect(pairs)) do pair
        first(pair) isa Source
    end
    unspecified_sources = setdiff(sources_, first.(specified_source_pairs))
    unspecified_sources_wrapping_something =
        filter(s -> !isempty(s), unspecified_sources)
    if !isempty(unspecified_sources_wrapping_something) &&
        !empty_unspecified_sources
        @warn "No replacement specified for one or more non-empty source "*
        "nodes. Contents will be duplicated. "
    end
    if empty_unspecified_sources
        unspecified_source_pairs = [s => source() for
                                    s in unspecified_sources]
    else
        unspecified_source_pairs = [s => deepcopy(s) for
                                    s in unspecified_sources]
    end

    all_source_pairs = vcat(specified_source_pairs, unspecified_source_pairs)

    # drop source nodes from all nodes of network terminating at W:
    nodes_ = filter(nodes(W)) do N
        !(N isa Source)
    end
    isempty(nodes_) && error("All nodes in network are source nodes. ")
    # instantiate node and machine dictionaries:
    newnode_given_old =
        IdDict{AbstractNode,AbstractNode}(all_source_pairs)
    newsources = [newnode_given_old[s] for s in sources_]
    newinterface_node_given_old =
        IdDict{AbstractNode,AbstractNode}()
    newmach_given_old = IdDict{Machine,Machine}()

    # build the new network:
    for N in nodes_
       args = [newnode_given_old[arg] for arg in N.args]
         if N.machine === nothing
             newnode_given_old[N] = node(N.operation, args...)
         else
             if N.machine in keys(newmach_given_old)
                 m = newmach_given_old[N.machine]
             else
                 train_args = [newnode_given_old[arg] for arg in N.machine.args]
                 m = Machine(newmodel_given_old[N.machine.model],
                                train_args...)
                 newmach_given_old[N.machine] = m
             end
             newnode_given_old[N] = N.operation(m, args...)
         end
        if N in interface_nodes
            newinterface_node_given_old[N] = newnode_given_old[N]
        end
    end

    newinterface_nodes = Tuple(newinterface_node_given_old[N] for N in
                          interface_nodes)
    newsignature = NamedTuple{keys(signature)}(newinterface_nodes)

    return machine(mach.model, newsources...; newsignature...)

end


text:
##### Data configuration file #####
# Copy this file to 'config/data.jl' and modify as needed

# source names by region
push!(dataconfig.sources,
    "Germany" => "JHU",
    "Italy" => "JHU",
    "France" => "JHU",
    "China" => "JHU",
    "US" => "JHU",
    "King" => "JHU",
    "New York" => "JHU",
    "United Kingdom" => "JHU",
    "Miami" => "JHU"
)

# data loaders by source
push!(dataconfig.loaders,
    "JHU" => "JHU",
)

# data dirs by source
push!(dataconfig.data_dirs,
    "native" => "native",
    "JHU" => "JHU",
)

# filename builder configuration
push!(dataconfig.args_filename_builder,
    "King" => ("King","Washington","US"),
    "New York" => ("New York","New York","US"),
    "Miami" => ("Miami-Dade","Florida","US"),
    "Laputa" => ("Laputa","Balnibarbi"),
)

# data loader arguments by region
push!(dataconfig.args_data_loader,
    "France" => (r"^$","France"),
    "King" => ("King","Washington","US"),
    "New York" => ("New York","New York","US"),
    "Miami" => ("Miami-Dade","Florida","US"),
)

# population loader arguments by region
push!(dataconfig.args_population_loader,
    "France" => (r"^$","France"),
    "King" => ("King","Washington"),
    "New York" => ("New York","New York"),
    "Miami" => ("Miami-Dade","Florida"),
    "Laputa" => ("Laputa","Balnibarbi"),
)

# data maps by loader
push!(dataconfig.data_maps,
    "JHU" => [0 1 1 1; 0 0 0 1],
    "native" => [0 1 1 1; 0 0 0 1],
)

# outbreak conditions by loader
push!(dataconfig.outbreak_conditions,
    "JHU" => cases -> findfirst(values(cases)[:,1].>0),
    "native" => cases -> findfirst(values(cases)[:,1].>0)
)


text:
#=  ================================================================================================
    Read partitioned dataframes from local or s3 paths.

    TODO:       dataset_root currently works with AbstractPath and S3Path types
                partition labels are detected by a simple regex (not timestamp formats)
                kwargs better
                add column type defintion for partition columns
                prune partions from the read stage via kwarg

    author:     mahiki@users.noreply.github.com
=   ===============================================================================================#
using DataFrames, Pipe, Dates
using FilePathsBase, AWS, AWSS3, CSV

"""
    read_partitioned_dataframe_csv(dataset_root, date_format::String, column_types::Dict)

Read and return a dataframe from the root of a local path (AbstractPath) or s3 path (S3Path).
Partition columns are typed as String.
Dataset partitions are detected from regex matches to the path using:

    r"/([\\w]+?)=([\\w-]+?)/"
"""
function read_partitioned_dataframe_csv(dataset_root, date_format::String, column_types::Dict)

    # word/digit column, word/digit value and hyphenated dates
    partition_regex = r"/([\w]+?)=([\w-]+?)/"

    datafiles = @pipe walkpath(dataset_root) |>
        collect |>
        filter(x -> extension(x) == "csv", _)

    partitions = @pipe datafiles |>
        eachmatch.(partition_regex, convert.(String, _), overlap = true) |>
        collect.(_)

    file_index = [
        [ datafiles[i], [ convert.(String, j.captures) for j in partitions[i] ]
        ] for i in keys(datafiles)
        ]

    # assemble the df with all paritions fro
    df = DataFrame()

    for i in keys(file_index)
        df_raw = @pipe file_index[i][1] |>
            read(_) |>
            CSV.read(_
                , DataFrame
                , dateformat = date_format
                , types = column_types
            )

        # add partition columns
        for parts in file_index[i][2]
            part_col = parts[1]
            part_val = parts[2]

            df_raw[:, "$part_col"] .= part_val
        end

        append!(df, df_raw, cols = :orderequal)
    end

    return df
end

max_stars_repo_name: mahiki/DesertIslandDisk

max_issues_repo_name: mahiki/DesertIslandDisk

max_forks_repo_name: mahiki/DesertIslandDisk


text:
function add_res_source!(repo::ResRepoTree, sources::Vector)
    for src in sources
        dst = joinpath(repo.source_dir, basename(src));
        cp(src, dst, follow_symlinks = true, remove_destination = true);
    end
end

function clear_res_source!(repo::ResRepoTree)
    for src in readdir(repo.source_dir)
        rm(joinpath(repo.source_dir, basename(src)), force = true, recursive = true)
    end
end

function set_res_source!(repo::ResRepoTree, sources::Vector)
    clear_res_source!(repo);
    add_res_source!(repo, sources);
end


text:
module conv_code

struct ConvolutionalCode
	quant_mem::Int64
	outs::Array{Array{Int64, 1}, 1}
	ConvolutionalCode(quant_mem, outs) = new(quant_mem, outs)
end

export ConvolutionalCode

end


text:
"""
    encoding!(actr, chunk, args...; kwargs...)

Sets imaginal module as busy and registers a new event to create a new `chunk`

# Arguments

- `actr`: an ACT-R model object
- `chunk`: a memory chunk
"""
function encoding!(actr, chunk)
    actr.imaginal.state.busy = true
    description = "Create New Chunk"
    type = "model"
    id = get_name(actr)
    tΔ = rnd_time(.200)
    register!(actr, encode!, after, tΔ , actr, chunk; id, description, type)
    return tΔ
end

"""
    encode!(actr, chunk, args...; kwargs...)

Completes the creation of a chunk and adds resulting `chunk` to the imaginal buffer. The buffer
states are set to busy = false and empty = false.

# Arguments

- `actr`: an ACT-R model object
- `chunk`: a memory chunk
"""
function encode!(actr, chunk)
    actr.imaginal.state.busy = false
    actr.imaginal.state.empty = false
    add_to_buffer!(actr.imaginal, chunk)
    return nothing
end


text:
cd(@__DIR__)
using Pkg
Pkg.activate(@__DIR__); Pkg.instantiate()

include("synapse_setup.jl")

indexof(sym,syms) = findfirst(isequal(sym), syms)

simplified_system = structural_simplify(network)
# since the index of states is non-deterministic, we look it up at runtime
idx = indexof(neuron1.sys.Vₘ, states(simplified_system))

byhand_prob = ODEProblem{true}(simple_synapse!, u0, (0.,t), p)
byhand_sol = solve(byhand_prob, Rosenbrock23(), reltol=1e-9, abstol=1e-9, saveat=0.025);
current_mtk_sol = solve(simul, Rosenbrock23(), reltol=1e-9, abstol=1e-9, saveat=0.025);

tsteps = 0.0:0.025:t
byhand_out = Array(byhand_sol(tsteps, idxs=3))
current_mtk_out = Array(current_mtk_sol(tsteps, idxs=idx))

using Test

@test isapprox(byhand_out, current_mtk_out, rtol=0.0001)


text:
using Ipopt

const IPO_OPTIMIZER = Ipopt.Optimizer(print_level=0)
MOI.set(IPO_OPTIMIZER, MOI.Silent(), true)
const IPO_CACHE = MOIU.UniversalFallback(MOIU.Model{Float64}())
const IPO_CACHED = MOIU.CachingOptimizer(IPO_CACHE, IPO_OPTIMIZER)
const IPO_BRIDGED = MOI.Bridges.full_bridge_optimizer(IPO_CACHED, Float64)
# MOI.Bridges.add_bridge(IPO_BRIDGED, MOI.Bridges.Constraint.SOCtoNonConvexQuadBridge{Float64})

push!(solvers_nlp_lowtol, (opt = IPO_BRIDGED, mode = BilevelJuMP.ProductMode{Float64}(1e-5)))
push!(solvers_nlp, (opt = IPO_BRIDGED, mode = BilevelJuMP.ProductMode{Float64}(1e-9)))
push!(solvers_nlp, (opt = IPO_BRIDGED, mode = BilevelJuMP.ProductWithSlackMode{Float64}(1e-9)))

push!(solvers, (opt = IPO_BRIDGED, mode = BilevelJuMP.ProductMode{Float64}(1e-9)))
push!(solvers_quad, (opt = IPO_BRIDGED, mode = BilevelJuMP.ProductMode{Float64}(1e-9)))

#=
push!(solvers, (opt = IPO_BRIDGED, mode = BilevelJuMP.StrongDualityInequalityMode{Float64}(1e-9)))
push!(solvers_quad, (opt = IPO_BRIDGED, mode = BilevelJuMP.StrongDualityInequalityMode{Float64}(1e-9)))
#push!(solvers_nlp, (opt = IPO_BRIDGED, mode = BilevelJuMP.StrongDualityInequalityMode{Float64}(1e-9)))
=#


text:
function DiscreteValueIteration.transition_matrix_a_s_sp(mdp::MDP; tranfile::String="transition.jld2")
    if isfile(tranfile)
        JLD2.@load tranfile trans_mat
        @info "Loading transition from $tranfile"
        return trans_mat
    else
        # Thanks to zach
        na = n_actions(mdp)
        ns = n_states(mdp)
        transmat_row_A = [Float64[] for _ in 1:n_actions(mdp)]
        transmat_col_A = [Float64[] for _ in 1:n_actions(mdp)]
        transmat_data_A = [Float64[] for _ in 1:n_actions(mdp)]

        @showprogress for s in states(mdp)
            si = stateindex(mdp, s)
            for a in actions(mdp, s)
                ai = actionindex(mdp, a)
                if !isterminal(mdp, s) # if terminal, the transition probabilities are all just zero
                    td = transition(mdp, s, a)
                    for (sp, p) in weighted_iterator(td)
                        if p > 0.0
                            spi = stateindex(mdp, sp)
                            push!(transmat_row_A[ai], si)
                            push!(transmat_col_A[ai], spi)
                            push!(transmat_data_A[ai], p)
                        end
                    end
                end
            end
        end
        transmats_A_S_S2 = [sparse(transmat_row_A[a], transmat_col_A[a], transmat_data_A[a], n_states(mdp), n_states(mdp)) for a in 1:n_actions(mdp)]
        # Note: not valid for terminal states
        # @assert all(all(sum(transmats_A_S_S2[a], dims=2) .≈ ones(n_states(mdp))) for a in 1:n_actions(mdp)) "Transition probabilities must sum to 1"
        return transmats_A_S_S2
    end
end

function DiscreteValueIteration.reward_s_a(mdp::MDP, rewfile::String="reward.jld2")
    if isfile(rewfile)
        JLD2.@load rewfile reward_mat
        @info "Loading reward from $rewfile"
        return reward_mat
    else
        reward_S_A = fill(-Inf, (n_states(mdp), n_actions(mdp))) # set reward for all actions to -Inf unless they are in actions(mdp, s)
        @showprogress for s in states(mdp)
            if isterminal(mdp, s)
                reward_S_A[stateindex(mdp, s), :] .= 0.0
            else
                for a in actions(mdp, s)
                    td = transition(mdp, s, a)
                    r = 0.0
                    for (sp, p) in weighted_iterator(td)
                        if p > 0.0
                            r += p*reward(mdp, s, a, sp)
                        end
                    end
                    reward_S_A[stateindex(mdp, s), actionindex(mdp, a)] = r
                end
            end
        end
        return reward_S_A
    end
end


text:
# Autogenerated wrapper script for assimp_win_jll for x86_64-w64-mingw32-cxx03
export assimp

using Zlib_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "PATH"
LIBPATH_default = ""

# Relative path to `assimp`
const assimp_splitpath = ["bin", "libassimp.dll"]

# This will be filled out by __init__() for all products, as it must be done at runtime
assimp_path = ""

# assimp-specific global declaration
# This will be filled out by __init__()
assimp_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const assimp = "libassimp.dll"

# Initialize PATH and LIBPATH environment variable listings.
# From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
# then append them to our own.
foreach(p -> append!(PATH_list, p), (Zlib_jll.PATH_list,))
foreach(p -> append!(LIBPATH_list, p), (Zlib_jll.LIBPATH_list,))
# Inform that the wrapper is available for this platform
is_available() = true

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"assimp_win")

    global PATH_list, LIBPATH_list
    global assimp_path = normpath(joinpath(artifact_dir, assimp_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global assimp_handle = dlopen(assimp_path)
    push!(LIBPATH_list, dirname(assimp_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ';')
    global LIBPATH = join(vcat(LIBPATH_list, [Sys.BINDIR, joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ';')

end  # __init__()


text:
using DistributedArrays, DistributedOperations, Documenter

makedocs(sitename="DistributedOperations", modules=[DistributedOperations])

deploydocs(
    repo = "git@github.com:ChevronETC/DistributedOperations.jl.git"
)


text:
cd(joinpath(@__DIR__, "."))
Pkg.activate(".")
using Chordal
using LinearAlgebra, SparseArrays, LightGraphs
using Plots: spy
using GraphPlot: gplot

function plot_clique_graph(cg::Chordal.CliqueGraph)
        nc = length(cg.active_cliques)
        ind = sort(collect(cg.active_cliques))
        reduced_edge_mat = cg.edge_mat[ind,ind]
        clique_graph = SimpleGraph(reduced_edge_mat)
        node_labels = ["$c" for (i, c) in enumerate(get_cliques(cg))]
        edge_labels = [string(reduced_edge_mat[src(e),dst(e)]) for e in edges(clique_graph)]
        plt = gplot(
                clique_graph,
                nodelabel=node_labels,
                edgelabel=edge_labels,
        )
end

## Example from A Clique Graph Based Merging Strategy for Decomposable SDPs
# by Michael Garstka, Mark Cannon, Paul Goulart
#       See Figures 1 and 3
#       Note: This does not include chordal extension (already chordal)
# Setup sparsity pattern
n = 9
nonzero_inds = vcat([
        (3,1), (3,2), (5,4), (6,1), (6,3), (7,3), (7,6), (8,3), (8,4), (8,5),
        (8,6), (8,7), (9,6), (9,7), (9,8)
    ],
    [(i,i) for i in 1:n]
)
sp = sparse(Chordal.unzip(nonzero_inds)..., [10i+j for (i, j) in nonzero_inds])
spy(sp, ms=10)
# c.f. Figure 1

## Find cliques
sp = sp + tril(sp, -1)'
Chordal.preprocess!(sp)

sparsity_graph = SimpleGraph(sp)
cliques = maximal_cliques(sparsity_graph)
@info join(vcat(["Cliques are:"], ["\n\t\tC$i: $(cliques[i])" for i in 1:length(cliques)]))
gplot(sparsity_graph, nodelabel=1:nv(sparsity_graph))

## Generate clique graph
cg = generate_clique_graph(cliques, n)
@info ("Initial cliques: $([string(a) for a in get_cliques(cg)])")
plot_clique_graph(cg)
# c.f. Figure 3

## Merge cliques
merge_cliques!(cg; verbose=true)
@info ("Final cliques: $([string(a) for a in get_cliques(cg)])")
plot_clique_graph(cg)
# c.f. Figure 3

## Selector Matrices
mat = Matrix(sp)
Cℓs = get_cliques(cg)
Tℓs = make_selectors_from_cliques(Cℓs, n)
println("Full Matrix:")
display(sp)
for i in 1:length(Cℓs)
        println("\nClique $i: $(string(Cℓs[i]))")
        display((Tℓs[i]*mat*Tℓs[i]'))
end


text:
module SchumakerSpline

include("SchumakerFunctions.jl")

export schumaker

end


text:
using OneHotArrays
using Test

@testset "OneHotArray" begin
  include("array.jl")
end

@testset "Constructors" begin
  include("onehot.jl")
end

@testset "Linear Algebra" begin
  include("linalg.jl")
end

using Zygote
import CUDA
if CUDA.functional()
  using CUDA  # exports CuArray, etc
  @info "starting CUDA tests"
else
  @info "CUDA not functional, testing via GPUArrays"
  using GPUArrays
  GPUArrays.allowscalar(false)

  # GPUArrays provides a fake GPU array, for testing
  jl_file = normpath(joinpath(pathof(GPUArrays), "..", "..", "test", "jlarray.jl"))
  using Random  # loaded within jl_file
  include(jl_file)
  using .JLArrays
  cu = jl
  CuArray{T,N} = JLArray{T,N}
end

@test cu(rand(3)) .+ 1 isa CuArray

@testset "GPUArrays" begin
  include("gpu.jl")
end


text:
using HelloWorld
using Test

@testset "HelloWorld.jl" begin
    @test "hello" == hello("hello")
end


text:
using aoc2019.computer: load_program, io, run
using aoc2019.robot: Point, Black, White, deploy, color, paint!, move!
using Match

input = joinpath(@__DIR__, "input")
p = load_program(input)

function activate_robot(program, initial_color)
    surface = Dict{Point, Int}()
    robot = deploy(surface)
    paint!(robot, initial_color)
    stdin, stdout = io()
    task = @async run(p, stdin, stdout)
    while !istaskdone(task)
        put!(stdin, color(robot))
        current_color = take!(stdout)
        turn = take!(stdout)
        paint!(robot, current_color)
        move!(robot, turn)
        yield()
    end
    return surface
end

ship = activate_robot(copy(p), Black)
p1 = length(collect(keys(ship)))
@assert p1 == 1951

ship = activate_robot(copy(p), White)
p2 = length(collect(keys(ship)))
@assert p2 == 249

xmin = minimum([p.x for p in collect(keys(ship))])
ymin = minimum([p.y for p in collect(keys(ship))])
xmax = maximum([p.x + abs(xmin) + 1 for p in collect(keys(ship))])
ymax = maximum([p.y + abs(ymin) + 1 for p in collect(keys(ship))])

panel = zeros(Int, (xmax, ymax))
for p in ship
    loc = p.first
    color = p.second
    panel[loc.x + abs(xmin) + 1, loc.y + abs(ymin) + 1] = color
end

print("-----------------------------------------------------------------------\n")
print("space police -- part one\n    panels painted      : $p1\n")
print("space police -- part two\n    identification code :\n")
for y in 1:ymax
    print("        ")
    for x in 1:xmax
        # print bottom to top to reflect the image and make it readable
        glyph = panel[x, ymax - y + 1] == 1 ? "#" : " "
        print("$glyph")
    end
    println()
end
print("-----------------------------------------------------------------------\n")


text:
using SparseIR
using Test

@testset "sampling.moveaxis" begin
    shape = (2, 3, 4)
    x = reshape(collect(1:prod(shape)), shape)
    ref = zeros(Int64, 3, 4, 2)
    for k in 1:4
        for j in 1:3
            for i in 1:2
                ref[j, k, i] = x[i, j, k]
            end
        end
    end
    res = SparseIR._move_axis(x, 1, 3)
    @test all(res .== ref)
end

@testset "sampling._matop_along_axis" begin
    op = reshape(collect(1:6), (3, 2))
    x = reshape(collect(1:4), (2, 2))
    res = SparseIR._matop_along_axis(op, x, 2)
    ref = zeros(Int64, (2, 3))
    for j in 1:2
        for i in 1:3
            for k in 1:2
                ref[j, i] += op[i, k] * x[j, k]
            end
        end
    end
    @test all(res .== ref)
end

@testset "sampling.XSampling" begin
    beta = 1e+4
    wmax = 1.0
    eps = 2e-8

    poles = wmax * Float64[-0.999, -0.01, 0.5]
    coeff = Float64[0.8, -0.2, 0.5]
    for XSampling in [TauSampling, MatsubaraSampling]
        for stat in [fermion, boson]
            basis = FiniteTempBasis(stat, beta, wmax, eps)
            rhol = basis.v(poles) * coeff
            gl = -basis.s .* rhol

            smpl = XSampling(basis)
            gtau = evaluate(smpl, gl)
            gl_reconst = fit(smpl, gtau)

            @test isapprox(gl, gl_reconst; rtol=0, atol=smpl.cond * eps * maximum(abs.(gl)))
        end
    end
end

@testset "sampling.XSamplingHigherDim" begin
    beta = 1e+4
    wmax = 1.0
    eps = 2e-8
    newaxis = [CartesianIndex()]

    npoles = 3
    poles = wmax * transpose(Float64[-0.999 -0.01 0.5; 0.5 -0.01 -0.999])
    coeff = transpose(Float64[0.8 -0.2 0.5; 0.5 -0.2 0.8])
    @test size(poles) == (npoles,  2)
    @test size(coeff) == (npoles,  2)
    for XSampling in [TauSampling, MatsubaraSampling]
        for stat in [fermion, boson]
            basis = FiniteTempBasis(stat, beta, wmax, eps)
            vval = basis.v(vec(poles))
            @assert size(vval) == (basis.size, npoles*2)
            vval = reshape(vval, (basis.size, npoles, 2))
            rhol = vval .* coeff[newaxis, :, :]
            @assert size(rhol) == (basis.size, npoles, 2)
            gl = -basis.s[:, newaxis, newaxis] .* rhol

            smpl = XSampling(basis)
            gtau = evaluate(smpl, gl, axis=1)
            gl_reconst = fit(smpl, gtau, axis=1)

            @test isapprox(gl, gl_reconst; rtol=0, atol=smpl.cond * eps * maximum(abs.(gl)))
        end
    end
end

@testset "sampling.return_types" begin
    @test Base.return_types(fit, (MatsubaraSampling, Array{Float64,3})) ==
          [Array{ComplexF64,3}]
    @test Base.return_types(fit, (MatsubaraSampling, Array{ComplexF64,3})) ==
          [Array{ComplexF64,3}]
    @test Base.return_types(fit, (TauSampling, Array{Float64,3})) == [Array{ComplexF64,3}]
    @test Base.return_types(fit, (TauSampling, Array{ComplexF64,3})) ==
          [Array{ComplexF64,3}]
end


text:
# Copyright 2021 The CLVR Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

# Script for executing algorithms on selected datasets on the DRO problem with Wasserstein metric
# based ambiguity sets. Command line usage:
# julia scripts/run_algo.jl <dataset> <gamma> <algo1>...<algo4>

using Arpack
using CSV
using Dates
using LinearAlgebra
using Logging
using SparseArrays

BLAS.set_num_threads(1)

include("../src/problems/standardLP.jl")
include("../src/algorithms/utils/exitcriterion.jl")
include("../src/algorithms/utils/results.jl")
include("../src/algorithms/utils/helper.jl")
include("../src/problems/dro/utils/libsvm_parser.jl")
include("../src/problems/dro/wasserstein.jl")

include("../src/algorithms/clvr_lazy_restart.jl")
include("../src/algorithms/pdhg_restart.jl")
include("../src/algorithms/purecd_restart.jl")
include("../src/algorithms/spdhg_restart.jl")

DATASET_INFO = Dict([
    ("a1a", (123, 1605)),
    ("a9a", (123, 32561)),
    ("gisette", (5000, 6000)),
    ("news20", (1355191, 19996)),
    ("rcv1", (47236, 20242)),
])

# Dataset parameters
outputdir = "./run_results/"
dataset = ARGS[1]
if !haskey(DATASET_INFO, dataset)
    throw(ArgumentError("Invalid dataset name supplied."))
end
dim_dataset, num_dataset = DATASET_INFO[dataset]
filepath = "./data/$(dataset).txt"

# Problem instance parameters
κ = 0.1
ρ = 10.

# Problem instance instantiation
yX_T = read_libsvm_into_yXT_sparse(filepath, dim_dataset, num_dataset)
A_T, b, c = droreformuation_wmetric_hinge_standardformnormalized(yX_T, κ, ρ)
problem = StandardLinearProgram(A_T, b, c)
L = svds(A_T, nsv = 1)[1].S[1]

# Exit criterion
maxiter = 1e12
maxtime = 3600 * 12
targetaccuracy = 1e-7
loggingfreq = 5
exitcriterion = ExitCriterion(maxiter, maxtime, targetaccuracy, loggingfreq)

# Common algo parameters
blocksize = 50
R = sqrt(blocksize)
γ = parse(Float64, ARGS[2])  # TODO: Use ArgParse
restartfreq = Inf  # For restart when metric halves, set restartfreq=Inf

timestamp = Dates.format(Dates.now(), "yyyy-mm-dd_HH-MM-SS-sss")
loggingfilename = "$(outputdir)/$(dataset)-$(join(ARGS[3:end], "_"))-execution_log-$(timestamp).txt"
io = open(loggingfilename, "w+")
logger = SimpleLogger(io)

println("timestamp = $(timestamp)")
println("Completed initialization.")

with_logger(logger) do

    @info "Running on $(dataset) dataset."
    @info "--------------------------------------------------"
    @info "κ = $(κ)"
    @info "ρ = $(ρ)"
    @info "--------------------------------------------------"
    @info "A_T has size: $(size(A_T))"
    @info "A_T has nnz: $(size(findnz(A_T)[1])[1]))"
    @info "nnz ratio: $(size(findnz(A_T)[1])[1] / (size(A_T)[1] * size(A_T)[2])))"
    @info "L = $(L)"
    @info "--------------------------------------------------"
    @info "maxiter = $(maxiter)"
    @info "maxtime = $(maxtime)"
    @info "targetaccuracy = $(targetaccuracy)"
    @info "loggingfreq = $(loggingfreq)"
    @info "--------------------------------------------------"
    @info "blocksize = $(blocksize)"
    @info "R = $(R)"
    @info "γ = $(γ)"
    @info "restartfreq = $(restartfreq)"

    if "0" in ARGS[3:end]  # TODO: Use algo names instead
        println("========================================")
        println("Running clvr_lazy_restart_x_y with blocksize=1.")

        println("clvr_blocksize = 1")
        println("clvr_R = 1.")

        r_clvr_lazy_restart = clvr_lazy_restart_x_y(
            problem,
            exitcriterion;
            blocksize=1,
            R=1.,
            γ=γ,
            restartfreq=restartfreq,
            io=io,
        )

        export_filename = "$(outputdir)/$(dataset)-clvr_lazy_restart_x_y-blocksize=1-R=1-$(timestamp).csv"
        exportresultstoCSV(r_clvr_lazy_restart, export_filename)

        println("========================================")
    end

    if "1" in ARGS[3:end]  # TODO: Use algo names instead
        println("========================================")
        println("Running clvr_lazy_restart_x_y.")

        clvr_R_multiplier = 1.0
        println("clvr_R_multiplier = $(clvr_R_multiplier)")

        r_clvr_lazy_restart = clvr_lazy_restart_x_y(
            problem,
            exitcriterion;
            blocksize=blocksize,
            R=R * clvr_R_multiplier,
            γ=γ,
            restartfreq=restartfreq,
            io=io,
        )

        export_filename = "$(outputdir)/$(dataset)-clvr_lazy_restart_x_y-$(timestamp).csv"
        exportresultstoCSV(r_clvr_lazy_restart, export_filename)

        println("========================================")
    end

    if "2" in ARGS[3:end]  # TODO: Use algo names instead
        println("========================================")
        println("Running pdhg_restart_x_y.")

        pdhg_L_multiplier = 1.0
        println("pdhg_L_multiplier = $(pdhg_L_multiplier)")

        r_pdhg_restart = pdhg_restart_x_y(
            problem,
            exitcriterion;
            L=L * pdhg_L_multiplier,
            γ=γ,
            restartfreq=restartfreq,
            io=io,
        )

        export_filename = "$(outputdir)/$(dataset)-pdhg_restart_x_y-$(timestamp).csv"
        exportresultstoCSV(r_pdhg_restart, export_filename)

        println("========================================")
    end

    if "3" in ARGS[3:end]  # TODO: Use algo names instead
        println("========================================")
        println("Running spdhg_restart_x_y.")

        spdhg_R_multiplier = 1.0
        println("spdhg_R_multiplier = $(spdhg_R_multiplier)")

        r_spdhg_restart = spdhg_restart_x_y(
            problem,
            exitcriterion;
            blocksize=blocksize,
            R=R * spdhg_R_multiplier,
            γ=γ,
            restartfreq=restartfreq,
            io=io,
        )

        export_filename = "$(outputdir)/$(dataset)-spdhg_restart_x_y-$(timestamp).csv"
        exportresultstoCSV(r_spdhg_restart, export_filename)

        println("========================================")
    end

    if "4" in ARGS[3:end]  # TODO: Use algo names instead
        println("========================================")
        println("Running purecd_restart_x_y.")

        println("purecd_blocksize = 1")
        println("purecd_R = 1")

        r_purecd_restart = purecd_restart_x_y(
            problem,
            exitcriterion;
            blocksize=1,
            R=1,
            γ=γ,
            restartfreq=restartfreq,
            io=io,
        )

        export_filename = "$(outputdir)/$(dataset)-purecd_restart_x_y-$(timestamp).csv"
        exportresultstoCSV(r_purecd_restart, export_filename)

        println("========================================")
    end
end

max_stars_repo_name: ericlincc/Efficient-GLP

max_issues_repo_name: ericlincc/Efficient-GLP

max_forks_repo_name: ericlincc/Efficient-GLP


text:
function createpropertydf(properties)

  prop_types = [Float64, Float64, Float64]
  prop_names = [:ea, :er, :gj]

  propertydf = DataFrame(prop_types, prop_names, 0)
  for i in 1:size(properties, 1)
      append!(propertydf,
        DataFrame(
          ea=properties[i, 1],
          er=properties[i, 2],
          gj=properties[i, 3]
        )
      )
  end

  propertydf

end


text:
#### Turbulence closures
using DocStringExtensions
using CLIMA.PlanetParameters: grav
using CLIMA.SubgridScaleParameters
export ConstantViscosityWithDivergence, SmagorinskyLilly, Vreman, AnisoMinDiss
export turbulence_tensors

abstract type TurbulenceClosure end

vars_state(::TurbulenceClosure, FT) = @vars()
vars_aux(::TurbulenceClosure, FT) = @vars()

function atmos_init_aux!(
    ::TurbulenceClosure,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
) end
function atmos_nodal_update_aux!(
    ::TurbulenceClosure,
    ::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function gradvariables!(
    ::TurbulenceClosure,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function diffusive!(
    ::TurbulenceClosure,
    ::Orientation,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
) end

"""
    ν, τ = turbulence_tensors(::TurbulenceClosure, state::Vars, diffusive::Vars, aux::Vars, t::Real)

Compute the kinematic viscosity tensor (`ν`) and SGS momentum flux tensor (`τ`).
"""
function turbulence_tensors end

"""
    principal_invariants(X)

Calculates principal invariants of a tensor `X`. Returns 3 element tuple containing the invariants.
"""
function principal_invariants(X)
    first = tr(X)
    second = (first^2 - tr(X .^ 2)) / 2
    third = det(X)
    return (first, second, third)
end

"""
    symmetrize(X)

Compute `(X + X')/2`, returning a `SHermitianCompact` object.
"""
function symmetrize(X::StaticArray{Tuple{3, 3}})
    SHermitianCompact(SVector(
        X[1, 1],
        (X[2, 1] + X[1, 2]) / 2,
        (X[3, 1] + X[1, 3]) / 2,
        X[2, 2],
        (X[3, 2] + X[2, 3]) / 2,
        X[3, 3],
    ))
end

"""
    norm2(X)

Compute
```math
\\sum_{i,j} S_{ij}^2
```
"""
function norm2(X::SMatrix{3, 3, FT}) where {FT}
    abs2(X[1, 1]) +
    abs2(X[2, 1]) +
    abs2(X[3, 1]) +
    abs2(X[1, 2]) +
    abs2(X[2, 2]) +
    abs2(X[3, 2]) +
    abs2(X[1, 3]) +
    abs2(X[2, 3]) +
    abs2(X[3, 3])
end
function norm2(X::SHermitianCompact{3, FT, 6}) where {FT}
    abs2(X[1, 1]) +
    2 * abs2(X[2, 1]) +
    2 * abs2(X[3, 1]) +
    abs2(X[2, 2]) +
    2 * abs2(X[3, 2]) +
    abs2(X[3, 3])
end

"""
    strain_rate_magnitude(S)

Compute
```math
|S| = \\sqrt{2\\sum_{i,j} S_{ij}^2}
```
"""
function strain_rate_magnitude(S::SHermitianCompact{3, FT, 6}) where {FT}
    return sqrt(2 * norm2(S))
end

"""
    ConstantViscosityWithDivergence <: TurbulenceClosure

Turbulence with constant dynamic viscosity (`ρν`).
Divergence terms are included in the momentum flux tensor.

# Fields

$(DocStringExtensions.FIELDS)
"""
struct ConstantViscosityWithDivergence{FT} <: TurbulenceClosure
    "Dynamic Viscosity [kg/m/s]"
    ρν::FT
end

vars_gradient(::ConstantViscosityWithDivergence, FT) = @vars()
vars_diffusive(::ConstantViscosityWithDivergence, FT) =
    @vars(S::SHermitianCompact{3, FT, 6})

function diffusive!(
    ::ConstantViscosityWithDivergence,
    ::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    diffusive.turbulence.S = symmetrize(∇transform.u)
end

function turbulence_tensors(
    m::ConstantViscosityWithDivergence,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)

    S = diffusive.turbulence.S
    ν = m.ρν / state.ρ
    τ = (-2 * ν) * S + (2 * ν / 3) * tr(S) * I
    return ν, τ
end

"""
    SmagorinskyLilly <: TurbulenceClosure

See § 1.3.2 in CliMA documentation

article{doi:10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2,
  author = {Smagorinksy, J.},
  title = {General circulation experiments with the primitive equations},
  journal = {Monthly Weather Review},
  volume = {91},
  number = {3},
  pages = {99-164},
  year = {1963},
  doi = {10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2},
  URL = {https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2},
  eprint = {https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2}
  }

article{doi:10.1111/j.2153-3490.1962.tb00128.x,
  author = {LILLY, D. K.},
  title = {On the numerical simulation of buoyant convection},
  journal = {Tellus},
  volume = {14},
  number = {2},
  pages = {148-172},
  doi = {10.1111/j.2153-3490.1962.tb00128.x},
  url = {https://onlinelibrary.wiley.com/doi/abs/10.1111/j.2153-3490.1962.tb00128.x},
  eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.2153-3490.1962.tb00128.x},
  year = {1962}
  }

Brunt-Vaisala frequency N² defined as in equation (1b) in
  Durran, D.R. and J.B. Klemp, 1982:
  On the Effects of Moisture on the Brunt-Väisälä Frequency.
  J. Atmos. Sci., 39, 2152–2158,
  https://doi.org/10.1175/1520-0469(1982)039<2152:OTEOMO>2.0.CO;2

# Fields

$(DocStringExtensions.FIELDS)
"""
struct SmagorinskyLilly{FT} <: TurbulenceClosure
    "Smagorinsky Coefficient [dimensionless]"
    C_smag::FT
end

vars_aux(::SmagorinskyLilly, FT) = @vars(Δ::FT)
vars_gradient(::SmagorinskyLilly, FT) = @vars(θ_v::FT)
vars_diffusive(::SmagorinskyLilly, FT) =
    @vars(S::SHermitianCompact{3, FT, 6}, N²::FT)

function atmos_init_aux!(
    ::SmagorinskyLilly,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.turbulence.Δ = lengthscale(geom)
end

function gradvariables!(
    m::SmagorinskyLilly,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.turbulence.θ_v = aux.moisture.θ_v
end

function diffusive!(
    ::SmagorinskyLilly,
    orientation::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    diffusive.turbulence.S = symmetrize(∇transform.u)
    ∇Φ = ∇gravitational_potential(orientation, aux)
    diffusive.turbulence.N² =
        dot(∇transform.turbulence.θ_v, ∇Φ) / aux.moisture.θ_v
end

function turbulence_tensors(
    m::SmagorinskyLilly,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)

    FT = eltype(state)
    S = diffusive.turbulence.S
    normS = strain_rate_magnitude(S)

    # squared buoyancy correction
    Richardson = diffusive.turbulence.N² / (normS^2 + eps(normS))
    f_b² = sqrt(clamp(1 - Richardson * inv_Pr_turb, 0, 1))
    ν = normS * f_b² * FT(m.C_smag * aux.turbulence.Δ)^2
    τ = (-2 * ν) * S
    return ν, τ
end

"""
    Vreman{FT} <: TurbulenceClosure

  §1.3.2 in CLIMA documentation
Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Smagorinsky coefficient
is specified and used to compute the equivalent Vreman coefficient.

1) ν_e = √(Bᵦ/(αᵢⱼαᵢⱼ)) where αᵢⱼ = ∂uⱼ∂uᵢ with uᵢ the resolved scale velocity component.
2) βij = Δ²αₘᵢαₘⱼ
3) Bᵦ = β₁₁β₂₂ + β₂₂β₃₃ + β₁₁β₃₃ - β₁₂² - β₁₃² - β₂₃²
βᵢⱼ is symmetric, positive-definite.
If Δᵢ = Δ, then β = Δ²αᵀα

@article{Vreman2004,
  title={An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications},
  author={Vreman, AW},
  journal={Physics of fluids},
  volume={16},
  number={10},
  pages={3670--3681},
  year={2004},
  publisher={AIP}
}

# Fields

$(DocStringExtensions.FIELDS)
"""
struct Vreman{FT} <: TurbulenceClosure
    "Smagorinsky Coefficient [dimensionless]"
    C_smag::FT
end
vars_aux(::Vreman, FT) = @vars(Δ::FT)
vars_gradient(::Vreman, FT) = @vars(θ_v::FT)
vars_diffusive(::Vreman, FT) = @vars(∇u::SMatrix{3, 3, FT, 9}, N²::FT)

function atmos_init_aux!(::Vreman, ::AtmosModel, aux::Vars, geom::LocalGeometry)
    aux.turbulence.Δ = lengthscale(geom)
end
function gradvariables!(
    m::Vreman,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.turbulence.θ_v = aux.moisture.θ_v
end
function diffusive!(
    ::Vreman,
    orientation::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.turbulence.∇u = ∇transform.u
    ∇Φ = ∇gravitational_potential(orientation, aux)
    diffusive.turbulence.N² =
        dot(∇transform.turbulence.θ_v, ∇Φ) / aux.moisture.θ_v
end

function turbulence_tensors(
    m::Vreman,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    α = diffusive.turbulence.∇u
    S = symmetrize(α)

    normS = strain_rate_magnitude(S)
    Richardson = diffusive.turbulence.N² / (normS^2 + eps(normS))
    f_b² = sqrt(clamp(1 - Richardson * inv_Pr_turb, 0, 1))

    β = f_b² * (aux.turbulence.Δ)^2 * (α' * α)
    Bβ = principal_invariants(β)[2]

    ν = max(0, m.C_smag^2 * FT(2.5) * sqrt(abs(Bβ / (norm2(α) + eps(FT)))))
    τ = (-2 * ν) * S

    return ν, τ
end

"""
    AnisoMinDiss{FT} <: TurbulenceClosure

  §1.3.2 in CLIMA documentation
Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Poincare coefficient
is specified and used to compute the equivalent AnisoMinDiss coefficient (computed as the solution to the
eigenvalue problem for the Laplacian operator).

@article{doi:10.1063/1.4928700,
author = {Rozema,Wybe  and Bae,Hyun J.  and Moin,Parviz  and Verstappen,Roel },
title = {Minimum-dissipation models for large-eddy simulation},
journal = {Physics of Fluids},
volume = {27},
number = {8},
pages = {085107},
year = {2015},
doi = {10.1063/1.4928700},
URL = {https://aip.scitation.org/doi/abs/10.1063/1.4928700},
eprint = {https://aip.scitation.org/doi/pdf/10.1063/1.4928700}
}
-------------------------------------------------------------------------------------
# TODO: Future versions will include modifications of Abkar(2016), Verstappen(2018)
@article{PhysRevFluids.1.041701,
title = {Minimum-dissipation scalar transport model for large-eddy simulation of turbulent flows},
author = {Abkar, Mahdi and Bae, Hyun J. and Moin, Parviz},
journal = {Phys. Rev. Fluids},
volume = {1},
issue = {4},
pages = {041701},
numpages = {10},
year = {2016},
month = {Aug},
publisher = {American Physical Society},
doi = {10.1103/PhysRevFluids.1.041701},
url = {https://link.aps.org/doi/10.1103/PhysRevFluids.1.041701}
}

"""
struct AnisoMinDiss{FT} <: TurbulenceClosure
    C_poincare::FT
end
vars_aux(::AnisoMinDiss, FT) = @vars(Δ::FT)
vars_gradient(::AnisoMinDiss, FT) = @vars()
vars_diffusive(::AnisoMinDiss, FT) = @vars(∇u::SMatrix{3, 3, FT, 9})

function atmos_init_aux!(
    ::AnisoMinDiss,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.turbulence.Δ = lengthscale(geom)
end
function gradvariables!(
    m::AnisoMinDiss,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function diffusive!(
    ::AnisoMinDiss,
    ::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.turbulence.∇u = ∇transform.u
end

function turbulence_tensors(
    m::AnisoMinDiss,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    α = diffusive.turbulence.∇u
    S = symmetrize(α)

    coeff = (aux.turbulence.Δ * m.C_poincare)^2
    βij = -(α' * α)
    ν = max(0, coeff * (dot(βij, S) / (norm2(α) + eps(FT))))
    τ = (-2 * ν) * S

    return ν, τ
end

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
module wind_module
wind_profs=Dict()
function save_wind4_module(wnd,name)
        push!(wind_profs,(name=>wnd))
end
end


text:
module LoggedDicts

using Logging

import Base: get, delete!, haskey, push!, pop!

export LoggedDict, set!, get, delete!,         # Create, Read, Delete
       pop!, push!,                            # Update
       write_logged_dict, read_logged_dict,
       haskey

type LoggedDict
    name::AbstractString    # LoggedDict is identifiable if it shares the log file with another source of logs
    dct::Dict
    logger::Logging.Logger
    write_counter::Int64    # Counts the writes

    function LoggedDict(name::AbstractString, dct::Dict, logfile::AbstractString, logging_off = false)
	lgr = Logger(name)
	if logging_off
	    Logging.configure(lgr, output = open(logfile, "a"), level = OFF)
	else
	    Logging.configure(lgr, output = open(logfile, "a"), level = INFO)
	end
	ld = new(name, dct, lgr, 0)
        log(ld, "INIT")
	ld
    end
    LoggedDict(name::AbstractString, logfile::AbstractString, logging_off=false) = LoggedDict(name, Dict(), logfile, logging_off)
end

"Write LoggedDict to disk."
function write_logged_dict(filename::AbstractString, ld::LoggedDict)
    f = open(filename, "r+")
    serialize(f, ld)
    close(f)
end

"Read LoggedDict from disk."
function read_logged_dict(filename::AbstractString)
    f  = open(filename, "r")
    ld = deserialize(f)
    close(f)
    ld
end

"Returns the value of ld::LoggedDict located at the path defined by keys..."
function get(ld::LoggedDict, keys...)
    dct = get_containing_dict(ld::LoggedDict, keys...)
    dct[keys[length(keys)]]
end

"Deletes the key-value pair located at the path defined by keys..."
function delete!(ld::LoggedDict, keys...)
    dct = get_containing_dict(ld::LoggedDict, keys...)
    log(ld, "DEL $keys")
    delete!(dct, keys[length(keys)])
end

"""
Sets the value of ld::LoggedDict located at the path defined by keys... equal to value.

Notes:
1. This function assumes that the final argument is the value.
2. If the specified path doesn't exist, it is created.
"""
function set!(ld::LoggedDict, keys_value...)
    nkeys = length(keys_value) - 1    # Final argument is the value (not a key)
    assert(nkeys > 0)
    dct = ld.dct
    for i in 1:(nkeys - 1)
        k = keys_value[i]
	if haskey(dct, k)
            if !(typeof(dct[k]) <: Dict)
                dct[k] = Dict()
	    end
	else
            dct[k] = Dict()
        end
        dct = dct[k]
    end
    k = keys_value[nkeys]
    log(ld, "SET $keys_value")
    dct[k] = keys_value[nkeys + 1]
end

"""
Pushes value onto the value of ld::LoggedDict located at the path defined by keys...

Notes:
1. This function assumes that the final argument is the value.
2. The specified path must exist, otherwise an error is raised.
"""
function push!(ld::LoggedDict, keys_value...)
    val_idx = length(keys_value)    # value is the last element of keys_value
    nkeys   = val_idx - 1           # Final argument is the value (not a key)
    keys    = keys_value[1:nkeys]
    dct     = get_containing_dict(ld::LoggedDict, keys...)
    log(ld, "PUSH $keys_value")
    push!(dct[keys_value[nkeys]], keys_value[val_idx])
end

"""
Pops value from the value of ld::LoggedDict located at the path defined by keys...

Notes:
1. This function assumes that the final argument is the value.
2. The specified path must exist, otherwise an error is raised.
"""
function pop!(ld::LoggedDict, keys_value...)
    val_idx = length(keys_value)    # value is the last element of keys_value
    nkeys   = val_idx - 1           # Final argument is the value (not a key)
    keys    = keys_value[1:nkeys]
    dct     = get_containing_dict(ld::LoggedDict, keys...)
    log(ld, "POP $keys_value")
    pop!(dct[keys_value[nkeys]], keys_value[val_idx])
end

################################################################################
### Utils
################################################################################

"Increments ld.write_counter, prepends write_counter to msg and writes msg to ld.logger."
function log(ld::LoggedDict, msg::AbstractString)
    ld.write_counter += 1
    info(ld.logger, "$(ld.write_counter): " * msg)
end

"""
Returns: Dict that contains the path defined by keys.
Raises an error if the path doesn't exist.
"""
function get_containing_dict(ld::LoggedDict, keys...)
    nkeys = length(keys)
    assert(nkeys > 0)
    dct = ld.dct
    for i in 1:(nkeys - 1)
        k = keys[i]
        if !haskey(dct, k)
            error("LoggedDict has no path $(keys[1:i]).")
        end
        dct = dct[k]
    end
    dct
end

function haskey(ld::LoggedDict, keys...)
    result = true
    try
	v = get(ld, keys...)
    catch
	result = false
    end
    result
end

end # module


text:
struct rotsymCTMIterable{T,
        TA <: AbstractTensor{T,4},
         TC <: AbstractTensor{T,2},
         TT <: AbstractTensor{T,3}} <: AbstractCTMIterable
    A::TA
    χ::Int
    Cinit::Union{TC,Nothing}
    Tinit::Union{TT,Nothing}
end

function rotsymctmiterable(A::DTensor{T,4}, χ::Int, Cinit = nothing, Tinit = nothing) where T
    rotsymCTMIterable{T,DTensor{T,4},DTensor{T,2},DTensor{T,3}}(A, χ, Cinit, Tinit)
end

function rotsymctmiterable(A::DASTensor{T,4,SYM,CHS,SS,CH}, χ::Int,
        Cinit = nothing, Tinit = nothing) where {T,N,SYM,CHS,SS,CH}
    TC = DASTensor{T,2,SYM,CHS,SS,CH}
    TT = DASTensor{T,3,SYM,CHS,SS,CH}
    TA = DASTensor{T,4,SYM,CHS,SS,CH}
    rotsymCTMIterable{T,TA,TC,TT}(A, χ, Cinit, Tinit)
end

struct rotsymCTMState{S,
        TA <: AbstractTensor,
        TC <: AbstractTensor,
        TT <: AbstractTensor} <: AbstractCTMState
    C::TC
    T::TT
    oldsvdvals::Vector{S}
    diffs::Vector{S}
    n_it::Ref{Int}
end

function iterate(iter::rotsymCTMIterable{S,TA,TC,TT}) where {S,TA,TC,TT}
    @unpack A, χ, Cinit, Tinit = iter
    if isnothing(Cinit)
        C = initializeC(A,χ)
        @tensor C[1,2] := C[1,2] + C[2,1]
    else
        C = Cinit
    end
    if isnothing(Tinit)
        T = initializeT(A,χ)
        @tensor T[1,2,3] := T[1,2,3] + T[3,2,1]
    else
        T = Tinit
    end

    l = ifelse(C isa DASTensor, 2χ, χ)
    oldsvdvals = zeros(real(S),l)
    state = rotsymCTMState{real(S),TA,TC,TT}(C, T, oldsvdvals, [], Ref(0))
    return state, state
end

function iterate(iter::rotsymCTMIterable, state::rotsymCTMState)
    @unpack A, χ = iter
    @unpack C, T, oldsvdvals, diffs, n_it = state
    #grow
    @tensor begin
        Cp[1,2,3,4]   := C[-1,-2]   * T[1,-3,-1] *
                         T[-2,-4,4] * A[2,3,-4,-3]
        Tp[1,2,3,4,5] := T[1,-1,5]  * A[3,4,-1 ,2]
    end

    #renormalize
    Z = tensorsvd(Cp, ((1,2),(3,4)), svdtrunc = svdtrunc_maxχ(χ))[1]
    @tensor begin
        Ctmp[1,2]   := Cp[-1,-2,-3,-4]   * Z[-1,-2,1] * Z'[-4,-3,2]
        Ttmp[1,2,3] := Tp[-1,-2,2,-3,-4] * Z'[-1,-2,1] * Z[-4,-3,3]
    end

    #symmetrize
    @tensor begin
        C[1,2]   = Ctmp[1,2] + Ctmp[2,1]
        T[1,2,3] = Ttmp[1,2,3] + Ttmp[3,2,1]
    end
    _, s, = tensorsvd(C)
    vals = diag(s)
    maxval = maximum(vals)
    apply!(C, x -> x .= x ./ maxval)
    apply!(T, x -> x .= x ./ maxval)
    normalize!(vals,1)

    #compare
    push!(diffs, sum(abs, oldsvdvals - vals))
    oldsvdvals[:] = vals
    n_it[] += 1

    return state, state
end


text:
import StatsBase: BaseRadixSortSafeTypes

"""
    sorttwo!(vs, index)

Sort both the `vs` and reorder `index` at the same. This allows for faster sortperm
for radix sort.
"""
function sorttwo!(vs::Vector{T}, index, lo::Int = 1, hi::Int=length(vs), RADIX_SIZE = 16, RADIX_MASK = 0xffff) where T # <:Union{BaseRadixSortSafeTypes}
    # Input checking
    if lo >= hi;  return (vs, index);  end
    #println(vs)
    o = Forward

    # Init
    iters = ceil(Integer, sizeof(T)*8/RADIX_SIZE)
    # number of buckets in the counting steuint_histp
    nbuckets = 2^RADIX_SIZE

    # Histogram for each element, radix
    bin = uint_hist(vs, RADIX_SIZE, RADIX_MASK)

    # bin = zeros(UInt32, nbuckets, iters)
    # if lo > 1;  bin[1,:] = lo-1;  end

    # Sort!
    swaps = 0
    len = hi-lo+1

    index1 = similar(index)
    ts=similar(vs)
    for j = 1:iters
        # Unroll first data iteration, check for degenerate case
        v = uint_mapping(o, vs[hi])
        idx = Int((v >> ((j-1)*RADIX_SIZE)) & RADIX_MASK) + 1

        # are all values the same at this radix?
        if bin[idx,j] == len;  continue;  end

        # cbin = cumsum(bin[:,j])
        # tries to achieve the above one-liner with more efficiency
        cbin = zeros(UInt32, nbuckets)
        cbin[1] = bin[1,j]
        @inbounds for i in 2:nbuckets
            cbin[i] = cbin[i-1] + bin[i,j]
        end

        ci = cbin[idx]
        #println((ci, hi))
        ts[ci] = vs[hi]
        index1[ci] = index[hi]
        # println(cbin[idx])
        cbin[idx] -= 1

        # Finish the loop...
        @inbounds for i in hi-1:-1:lo
            v = uint_mapping(o, vs[i])
            idx = Int((v >> ((j-1)*RADIX_SIZE)) & RADIX_MASK) + 1
            ci = cbin[idx]
            #println(ci)
            ts[ci] = vs[i]
            index1[ci] = index[i]
            cbin[idx] -= 1
        end
        vs,ts = ts,vs
        index, index1 = index1, index
        swaps += 1
    end

    if isodd(swaps)
        vs,ts = ts,vs
        index, index1 = index1, index
        for i = lo:hi
            @inbounds vs[i] = ts[i]
            @inbounds index[i] = index1[i]
        end
    end
    (vs, index)
end


text:
module MemoryHunter

function print_sizes(threshold=5#=mb=#)
    @show Base.format_bytes(Sys.free_memory())
    biggest_vars = Pair{String, Float64}[]
    visited = IdDict()
    # just make sure the numbers just look different because
    # of gc running inconistantly
    GC.gc(true)
    for m in values(Base.loaded_modules)
        if !haskey(visited, m)
            visited[m] = true
            for vs in names(m, all = true)
                if isdefined(m, vs)
                    v = getfield(m, vs)
                    if !isa(v, Function) && !haskey(visited, v)
                        visited[v] = true
                        x = Base.summarysize(v)
                        if x > threshold * 10^6
                            push!(biggest_vars, string(m, ".", vs) => x)
                        end
                    end
                end
            end
        end
    end
    sort!(biggest_vars, by=last, rev=true)
    for (var, s) in biggest_vars
        println(var, ": ", Base.format_bytes(s))
    end
end

end


text:
include("eggify.jl")
using Metatheory.Library
using Metatheory.EGraphs.Schedulers

mult_t = commutative_monoid(:(*), 1)
plus_t = commutative_monoid(:(+), 0)

minus_t = @theory begin
  a - a    => 0
  a + (-b) => a - b
end

mulplus_t = @theory begin
  0 * a => 0
  a * 0 => 0
  a * (b + c) == ((a * b) + (a * c))
  a + (b * a) => ((b + 1) * a)
end

pow_t = @theory begin
  (y^n) * y => y^(n + 1)
  x^n * x^m == x^(n + m)
  (x * y)^z == x^z * y^z
  (x^p)^q == x^(p * q)
  x^0 => 1
  0^x => 0
  1^x => 1
  x^1 => x
  inv(x) == x^(-1)
end

function customlt(x, y)
  if typeof(x) == Expr && Expr == typeof(y)
    false
  elseif typeof(x) == typeof(y)
    isless(x, y)
  elseif x isa Symbol && y isa Number
    false
  else
    true
  end
end

canonical_t = @theory begin
  # restore n-arity
  (x + (+)(ys...)) => +(x, ys...)
  ((+)(xs...) + y) => +(xs..., y)
  (x * (*)(ys...)) => *(x, ys...)
  ((*)(xs...) * y) => *(xs..., y)

  (*)(xs...) |> Expr(:call, :*, sort!(xs; lt = customlt)...)
  (+)(xs...) |> Expr(:call, :+, sort!(xs; lt = customlt)...)
end

cas = mult_t ∪ plus_t ∪ minus_t ∪ mulplus_t ∪ pow_t
theory = cas

query = cleanast(:(a + b + (0 * c) + d))

function simplify(ex)
  g = EGraph(ex)
  params = SaturationParams(
    scheduler = BackoffScheduler,
    timeout = 20,
    schedulerparams = (1000, 5), # fuel and bantime
  )
  report = saturate!(g, cas, params)
  println(report)
  res = extract!(g, astsize)
  res = rewrite(res, canonical_t; clean = false, m = @__MODULE__) # this just orders symbols and restores n-ary plus and mult
  res
end

###########################################

params = SaturationParams(timeout = 20, schedulerparams = (1000, 5))

for i in 1:2
  ex = simplify(:(a + b + (0 * c) + d))
  println("Best found: $ex")
end

open("src/main.rs", "w") do f
  write(f, rust_code(theory, query))
end


text:
module BitemporalPostgres
import SearchLight: AbstractModel, DbId, query, save!
import SearchLightPostgreSQL
import Base: @kwdef
import Intervals, Dates, TimeZones
using Intervals,
    Dates,
    SearchLight,
    SearchLight.Transactions,
    SearchLightPostgreSQL,
    TimeZones,
    AbstractTrees
include("DDL.jl")
using .DDL
export History,
    Version,
    Component,
    get_typeof_revision,
    SubComponent,
    ComponentRevision,
    ValidityInterval,
    TestDummyComponent,
    TestDummyComponentRevision,
    TestDummySubComponent,
    TestDummySubComponentRevision,
    Testtstzrange,
    Workflow,
    findversion,
    findcomponentrevision,
    create_entity!,
    create_component!,
    create_subcomponent!,
    update_entity!,
    update_component!,
    delete_component!,
    commit_workflow!,
    rollback_workflow!,
    MaxVersion,
    MaxDate,
    InfinityKey,
    up,
    down,
    Node,
    print_tree,
    mkforest

"""
  Workflow
  Holds application based transaction data.
  * current version
  * world validFrom date
  * db validfrom date
  * is_committed means: bitemporal transaction is
      * pending: 0
      * committed: 1
      * rolled back: 2

"""
@kwdef mutable struct Workflow <: AbstractModel
    id::DbId = DbId()
    ref_history::DbId = 0
    ref_version::DbId = 0
    tsdb_validfrom::ZonedDateTime = MaxDate
    tsw_validfrom::ZonedDateTime = MaxDate
    is_committed::Integer = 0
end

@kwdef mutable struct Testtstzrange <: AbstractModel
    id::DbId = DbId()
    i::Interval{ZonedDateTime,Closed,Open} = Interval{ZonedDateTime,Closed,Open}(
        ZonedDateTime(2013, 2, 13, 0, 30, 0, 1, tz"America/Winnipeg"),
        ZonedDateTime(2030, 2, 13, 0, 30, 0, 1, tz"America/Winnipeg"),
    )
end

"""
  History

  The root of an aggregate entity, which groups versions and entities

"""
@kwdef mutable struct History <: AbstractModel
    id::DbId = DbId()
    dummy::Integer = 0
end
# One ms less than the maximum date, as
const MaxDate = ZonedDateTime(DateTime(2038, 1, 19, 3, 14, 6, 999), tz"UTC")
const MaxDateSQL = SQLInput(MaxDate)
const InfinityKey = 9223372036854775807::Integer
const MaxVersion = DbId(InfinityKey)
"""
  Version

  version number of the state of aggregate entity

"""
@kwdef mutable struct Version <: AbstractModel
    id::DbId = DbId()
    ref_history::DbId = InfinityKey
end

"""
ValidityInterval

  a 2-dimensional interval of validity of a version
  with respect to the state

  * of the data base storage and
  * of the world
  See the theory page for background

"""
@kwdef mutable struct ValidityInterval <: AbstractModel
    id::DbId = DbId()
    ref_history::DbId = InfinityKey
    ref_version::DbId = InfinityKey
    tsdb_validfrom::ZonedDateTime = now(tz"UTC")
    tsworld_validfrom::ZonedDateTime = now(tz"UTC")
    tsdb_invalidfrom::ZonedDateTime = MaxDate
    tsworld_invalidfrom::ZonedDateTime = MaxDate
    is_committed::Integer = 0
end

"""
Component

  an abstract component of a bitemporal entity
  i.e. a contract, a reference to a partner
  a component's history of states is represented by its revisions

"""
abstract type Component <: AbstractModel end

function getid(component::Component)::DbId
    component.id
end

function getref_history(component::Component)::DbId
    component.ref_history
end

function setref_history(component::Component, ref::DbId)
    component.ref_history = ref
end

function getref_version(component::Component)::DbId
    component.ref_version
end

function setref_version(component::Component, ref::DbId)
    component.ref_version = ref
end

"""
SubComponent

  an abstract subcomponent of a bitemporal entity
  i.e. product item of a contract

"""
abstract type SubComponent <: Component end

function getref_super(component::Component)::DbId
    component.ref_super
end

function setref_super(component::SubComponent, ref::DbId)
    component.ref_super = ref
end

"""
ComponentRevision

  a component of a bitemporal entity
  i.e. a contract, a reference to a partner
  a component's history of states is represented by its revisions

"""
abstract type ComponentRevision <: AbstractModel end

"""
get_typeof_revision(component::T) :: R where {T<: Component, R <: ComponentRevision}
    returns the actual subtype of ComponentRevision that fits the actual type of component
"""
function get_typeof_revision(
    component::T,
)::Type{R} where {T<:Component,R<:ComponentRevision}
    Type{R}
end

function getid(rev::ComponentRevision)::DbId
    rev.id
end

function getref_component(rev::ComponentRevision)::DbId
    rev.ref_component
end

function setref_component(rev::ComponentRevision, ref::DbId)
    rev.ref_component = ref
end

function getref_validfrom(ref::DbId)::ComponentRevision
    rev.ref_validfrom = ref
end

function setref_validfrom(rev::ComponentRevision, ref::DbId)
    rev.ref_validfrom = ref
end

function getref_invalidfrom(ref::DbId)::ComponentRevision
    rev.ref_invalidfrom = ref
end

function setref_invalidfrom(rev::ComponentRevision, ref::DbId)
    rev.ref_invalidfrom = ref
end

"""
TestDummyComponent <: Component

  a component of a bitemporal entity

"""
@kwdef mutable struct TestDummyComponent <: Component
    id::DbId = DbId()
    ref_history::DbId = InfinityKey
    ref_version::DbId = InfinityKey
end

"""
get_typeof_revision(component::TestDummyComponent) :: T where {T<:ComponentRevision}
"""
function get_typeof_revision(component::TestDummyComponent)::T where {T<:ComponentRevision}
    TestDummyComponentRevision
end

"""
TestDummyComponentRevision <: ComponentRevision

  a revision of a contract component of a bitemporal entity

"""
@kwdef mutable struct TestDummyComponentRevision <: ComponentRevision
    id::DbId = DbId()
    ref_component::DbId = InfinityKey
    ref_validfrom::DbId = InfinityKey
    ref_invalidfrom::DbId = InfinityKey
    description::String = ""
end

Base.copy(src::TestDummyComponentRevision) = TestDummyComponentRevision(
    ref_component = src.ref_component,
    description = src.description,
)

"""
TestDummySubComponent <: Component
  a component of a bitemporal entity

"""
@kwdef mutable struct TestDummySubComponent <: SubComponent
    id::DbId = DbId()
    ref_history::DbId = InfinityKey
    ref_super::DbId = InfinityKey
    ref_version::DbId = InfinityKey
end

"""
get_typeof_revision(component::TestDummySubComponent) :: T where {T<:ComponentRevision}
"""
function get_typeof_revision(
    component::TestDummySubComponent,
)::Type{TestDummySubComponentRevision}
    TestDummySubComponentRevision
end

"""
TestDummySubComponentRevision <: ComponentRevision

  a revision of a contract component of a bitemporal entity

"""
@kwdef mutable struct TestDummySubComponentRevision <: ComponentRevision
    id::DbId = DbId()
    ref_component::DbId = InfinityKey
    ref_validfrom::DbId = InfinityKey
    ref_invalidfrom::DbId = InfinityKey
    description::String = ""
end

Base.copy(src::TestDummySubComponentRevision) = TestDummySubComponentRevision(
    ref_component = src.ref_component,
    description = src.description,
)
"""
findversion(ref_history::DbId, tsdb::ZonedDateTime, tsw::ZonedDateTime)::DbId
       retrieves the version_id of a bitemporal history asof tsdb as per tsw
"""
function findversion(ref_history::DbId, tsdb::ZonedDateTime, tsw::ZonedDateTime)::DbId
    find(
        ValidityInterval,
        SQLWhereExpression(
            "ref_history=? and tsrworld @> TIMESTAMPTZ ? AND tsrdb @> TIMESTAMPTZ ?",
            ref_history,
            SQLInput(tsw),
            SQLInput(tsdb),
        ),
    )[1].ref_version
end

"""
findcomponentrevision(t::Type{T},ref_component::DbId,ref_version::DbId,)::Vector{T} where {T<:ComponentRevision}
    retrieves the version_id of a bitemporal history asof tsdb as per tsw
"""
function findcomponentrevision(
    t::Type{T},
    ref_component::DbId,
    ref_version::DbId,
)::Vector{T} where {T<:ComponentRevision}
    find(
        t,
        SQLWhereExpression(
            "ref_component = ? AND ref_valid @> BIGINT ? ",
            ref_component,
            ref_version,
        ),
    )
end

"""
create_entity!(w::Workflow)
* opens a bitemporal transaction identified by ref_version
* persists a history, version, validityInterval and a Workflow
* requires: w.tsw_validfrom is a valid date

"""
function create_entity!(w::Workflow)
    transaction() do
        h = History()
        save!(h)
        v = Version(ref_history = h.id)

        save!(v)

        w.ref_history = h.id
        w.ref_version = v.id
        save!(w)

        i = ValidityInterval(
            ref_history = h.id,
            ref_version = v.id,
            tsworld_validfrom = w.tsw_validfrom,
            tsworld_invalidfrom = MaxDate,
            tsdb_validfrom = now(tz"Africa/Porto-Novo"),
            tsdb_invalidfrom = MaxDate,
        )
        save!(i)
    end
end

"""
create_component!(c::Component, cr :: ComponentRevision, w:: Workflow)

creates a component and a componentRevision

"""
function create_component!(c::Component, cr::ComponentRevision, w::Workflow)
    transaction() do
        hid = w.ref_history
        vid = w.ref_version
        setref_history(c, hid)
        setref_version(c, vid)
        save!(c)

        setref_component(cr, getid(c))
        setref_validfrom(cr, vid)
        save!(cr)
    end
end

"""
create_subcomponent!(c::Component, sc::SubComponent, cr :: ComponentRevision, w:: Workflow)
creates a subcomponent and a componentRevision
"""
function create_subcomponent!(
    c::Component,
    sc::SubComponent,
    cr::ComponentRevision,
    w::Workflow,
)
    transaction() do
        hid = w.ref_history
        vid = w.ref_version
        setref_super(sc, c.id)
        setref_history(sc, hid)
        save!(sc)

        setref_component(cr, getid(sc))
        setref_validfrom(cr, vid)
        save!(cr)
    end
end

"""
update_entity!(w::Workflow)
* opens a bitemporal transaction identified by ref_version
* persists a version, a validityInterval and a Workflow
* requires:
    * w.tsw_validfrom is a valid date
    * w.ref_history is a valid history id
    * w.ref_version is a valid version of w.ref_history

"""
function update_entity!(w::Workflow)
    transaction() do
        hid = w.ref_history
        v = Version(ref_history = hid)
        save!(v)
        w.ref_version = v.id
        save!(w)

        i = ValidityInterval(
            ref_history = hid,
            ref_version = v.id,
            tsworld_validfrom = w.tsw_validfrom,
            tsworld_invalidfrom = MaxDate,
            tsdb_validfrom = now(tz"Africa/Porto-Novo"),
            tsdb_invalidfrom = MaxDate,
        )
        save!(i)
    end
end

"""
update_component!(cr :: ComponentRevision, crNew :: ComponentRevision, w:: Workflow)

terminates a componentRevision and persists its successor

"""
function update_component!(cr::ComponentRevision, crNew::ComponentRevision, w::Workflow)
    transaction() do
        vid = w.ref_version
        setref_invalidfrom(cr, vid)
        save!(cr)

        setref_validfrom(crNew, vid)
        save!(crNew)
    end
end

"""
delete_component!(c::T, w::Workflow)  where {T<:Component}

 * deletes a component if it was created for the current version or
 * mark its latest component revision as invalid

"""
function delete_component!(c::T, w::Workflow) where {T<:Component}
    transaction() do
        vid = w.ref_version
        revs = find(get_typeof_revision(c), SQLWhereExpression("ref_component=?", c.id))
        if length(revs) == 1 && revs[1].ref_validfrom == vid
            """
            c was created for current version, foreign key constraint cascades  deletion to revision
            """
            delete(c)
        else
            for r in revs
                if vid.value >= r.ref_validfrom.value && vid.value < r.ref_invalidfrom.value
                    """
                    the current revision gets terminated
                    """
                    r.ref_invalidfrom = vid
                    save!(r)
                end
            end
        end
    end
end

"""
commit_workflow!(w::Workflow)
  commits the bitemporal transaction identified by ref_version
  begun with the creation of w
  * sets w's version's validityInterval to is_committed
  * ends all overlapping intervals to tsdb_invalidfrom w.tsdb_validfrom
  * creates a new interval for the non_overlapping part of a partially interval if any
"""
function commit_workflow!(w::Workflow)
    transaction() do
        w.tsdb_validfrom = now(tz"Africa/Porto-Novo")

        uncommitted = find(
            ValidityInterval,
            SQLWhereExpression("ref_version = ?  AND is_committed = 0", w.ref_version),
        )

        shadowed = find(
            ValidityInterval,
            SQLWhereExpression(
                "ref_history = ?  AND tsrdb @> TIMESTAMPTZ ? AND tsrworld <@ tstzrange(?,?) AND is_committed=1",
                w.ref_history,
                MaxDate - Dates.Day(1),
                uncommitted[1].tsworld_validfrom + Dates.Day(1),
                uncommitted[1].tsworld_invalidfrom,
            ),
        )

        for i in shadowed
            i.tsdb_invalidfrom = w.tsdb_validfrom
            save!(i)
        end

        overlapped = find(
            ValidityInterval,
            SQLWhereExpression(
                "ref_history = ?  AND tsrdb @> TIMESTAMPTZ ? AND tsrworld && tstzrange(?,?) AND is_committed=1",
                w.ref_history,
                MaxDate - Dates.Day(1),
                uncommitted[1].tsworld_validfrom + Dates.Day(1),
                uncommitted[1].tsworld_invalidfrom,
            ),
        )
        for i in overlapped
            i.tsdb_invalidfrom = w.tsdb_validfrom
            j = ValidityInterval(
                ref_history = i.ref_history,
                ref_version = i.ref_version,
                tsworld_validfrom = i.tsworld_validfrom,
                tsworld_invalidfrom = uncommitted[1].tsworld_validfrom,
                tsdb_validfrom = w.tsdb_validfrom,
                tsdb_invalidfrom = MaxDate,
                is_committed = 1,
            )
            save!(i)
            save!(j)
        end
        for i in uncommitted
            i.is_committed = 1
            i.tsdb_validfrom = w.tsdb_validfrom
            save!(i)
        end
        w.is_committed = 1
        save!(w)
    end
end

"""
rollback_workflow!(w::Workflow)
    rolls back the bitemporal transaction identified by ref_version
    begun with the creation of w
"""
function rollback_workflow!(w::Workflow)
    transaction() do
        v = find(Version, SQLWhereExpression("id=?", w.ref_version))[1]
        delete(v)
        w.is_committed = 2
        save(w)
    end
end

"""
Node

node of tree of mutations
"""
struct Node
    interval::ValidityInterval
    shadowed::Vector{Node}
end

function children(n::Node)
    return n.shadowed
end

printnode(io::IO, n::Node) =
    print(io, n.interval.ref_history.value * n.interval.ref_version.value)

"""
mkforest(hid::DbId,tsdb_invalidfrom::ZonedDateTime,tsworld_validfrom::ZonedDateTime,tsworld_invalidfrom::ZonedDateTime,)::Vector{Node}

    builds a tree of vectors of version nodes where

    * nodes in one vector denote conecutive mutations and
    * child node vectors denote mutations which have been retrospectively corrected by their predecessor

    see: Theory: Textual representation of mutation histories

"""
function mkforest(
    hid::DbId,
    tsdb_invalidfrom::ZonedDateTime,
    tsworld_validfrom::ZonedDateTime,
    tsworld_invalidfrom::ZonedDateTime,
)::Vector{Node}
    map(
        i::ValidityInterval -> Node(
            i,
            mkforest(hid, i.tsdb_validfrom, i.tsworld_validfrom, i.tsworld_invalidfrom),
        ),
        find(
            ValidityInterval,
            SQLWhereExpression(
                "ref_history=? AND  upper(tsrdb)=? AND tstzrange(?,?) * tsrworld = tsrworld",
                hid,
                tsdb_invalidfrom,
                tsworld_validfrom,
                tsworld_invalidfrom,
            ),
        ),
    )
end

end # module


text:
using SimulacoesTemplate
using Test

@testset "SimulacoesTemplate.jl" begin
   @test f(2) == 4
end


text:
include("constants.jl")

#Intensity
function intensity_fq(ν,T)
    I_λ= (2 * ν^2 * k_b * T) / (c^2)
    return I_λ
end

function intensity_wl(λ,T)
    I_λ= (2 * c * k_b * T) ./ (λ.^4)
    return I_λ
end

#Flux received
function flux_ll(R,r,ν,T)
    f_ll = (R / r)^2 * π * intensity_fq(ν,T)
    return f_ll
end

#Observed flux
function flux_obs(A,R,r,ν,T)
    f_obs = A * flux_ll(R,r,ν,T)
    return f_obs
end

#Δν
function flux_obs_dnu(Δν,R,r,ν,T)
    f_obs_nu = flux_ll(R,r,ν,T) * Δν
    return f_obs_nu
end

export intensity_fq, intensity_wl,
       flux_ll, flux_obs, flux_obs_dnu