text:
using FieldDistributionNonuniformMedium
using Documenter

DocMeta.setdocmeta!(FieldDistributionNonuniformMedium, :DocTestSetup, :(using FieldDistributionNonuniformMedium); recursive=true)

makedocs(;
    modules=[FieldDistributionNonuniformMedium],
    authors="JingYu Ning <foldfelis@gmail.com> and contributors",
    repo="https://github.com/foldfelis/FieldDistributionNonuniformMedium.jl/blob/{commit}{path}#{line}",
    sitename="FieldDistributionNonuniformMedium.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://foldfelis.github.io/FieldDistributionNonuniformMedium.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Simulation" => "simulation.md",
        "APIs" => "api.md",
    ],
)

deploydocs(;
    repo="github.com/foldfelis/FieldDistributionNonuniformMedium.jl",
    devbranch="master",
)


text:
# this file is useful for interactive REPL testing

# include("arima.jl")
# include("diagnostic_tests.jl")
# include("garch.jl")
include("kalman_filter.jl")

ext: jl


text:
pkgs = ["NDFiles", "TrackRoots"]
for pkg in pkgs
    if Pkg.cd(!isdir, pkg)
        Pkg.clone("https://github.com/yakir12/$pkg.jl")
        Pkg.build(pkg)
    end
end


text:
function main(::Robot)
  i = 1
  while isborder(r, Nord)
    move_to_the_side!(r, i, West)
    move_to_the_side!(r, i, Ost)
    move_to_the_side!(r, i, Ost)
    move_to_the_side!(r, i, West)
    i += 1
  end
end

function move_to_the_side!(::Robot, j::Int, side::HorizonSide)
  for i in 1:j
    if isborder(r,Nord)
      move!(r, side)
    end
  end
end

max_stars_repo_name: VladislavBukin/mirea-progs

max_issues_repo_name: VladislavBukin/mirea-progs

max_forks_repo_name: VladislavBukin/mirea-progs


text:
module NetDecOPF

# Import necessary packages
import JuMP
import PowerModels
import DualDecomposition
import BundleMethod
import LinearAlgebra: I
const PM = PowerModels
const DD = DualDecomposition
const BM = BundleMethod

using Metis
using LightGraphs, SimpleWeightedGraphs

include("model.jl")
include("data.jl")
include("prob.jl")
include("solution.jl")
include("utils.jl")

export NetDecModel, get_partition, get_models, get_split_vars, set_subnet_optimizer!

export find_neighbor_buses, generate_subnet_data, ref_add_cut_bus!, ref_add_cut_branch!

export decompose, build_acopf_with_free_lines, build_socbf_with_free_lines, build_acots_with_free_lines, collect_split_vars

export init_DD_algo, add_split_vars_to_algo!

export metis_cluster

include("prob_ext.jl")

export W_ACRModel, WU_ACRModel

end # module


text:
export LowRankMatrix, OuterProduct, projection, outerprod
abstract type LowRankMatrix{T} <: AbstractMatrix{T} end

"""
    OuterProduct{T,AT<:AbstractArray{T}} <: LowRankMatrix{T}

If `AT <: AbstractVector`, it represents an outer product `x.left*transpose(x.right)`.
Else if `AT <: AbstractMatrix`, then it is an outer product `x.left*transpose(x.right)`.
"""
struct OuterProduct{T,AT<:AbstractArray{T}} <: LowRankMatrix{T}
    left::AT
    right::AT
end

function OuterProduct(left::AT, right::AT) where {T,AT<:AbstractMatrix{T}}
    size(left,2) != size(right,2) && throw(DimensionMismatch("The seconds dimension of left ($(size(left,2))) and right $(size(right,2)) does not match."))
    return OuterProduct{T,AT}(left, right)
end

const BatchedOuterProduct{T,MT} = OuterProduct{T,MT} where {MT<:AbstractMatrix}
const SimpleOuterProduct{T,VT} = OuterProduct{T,VT} where {VT<:AbstractVector}

LinearAlgebra.rank(op::OuterProduct) = size(op.left, 2)
Base.getindex(op::OuterProduct{T,<:AbstractVector}, i::Int, j::Int) where {T} =
    op.left[i] * op.right[j]
Base.getindex(op::BatchedOuterProduct, i::Int, j::Int) =
    sum(k -> op.left[i, k] * op.right[j, k], 1:rank(op))
Base.size(op::OuterProduct) = (size(op.left, 1), size(op.right, 1))
Base.size(op::OuterProduct, i::Int) =
    i == 1 ? size(op.left, 1) : (i == 2 ? size(op.right, 1) : throw(DimensionMismatch("")))

LinearAlgebra.mul!(a::OuterProduct, b) = a.left * (transpose(a.right) * b)
LinearAlgebra.mul!(a::OuterProduct, b::AbstractMatrix) =
    OuterProduct(a.left, transpose(transpose(a.right) * b))
LinearAlgebra.mul!(a, b::OuterProduct) = (a * b.left) * transpose(b.right)
LinearAlgebra.mul!(a::AbstractMatrix, b::OuterProduct) = OuterProduct(a * b.left, b.right)
LinearAlgebra.mul!(a::OuterProduct, b::OuterProduct) =
    OuterProduct(a.left * (transpose(a.right) * b.left), b.right)

Base.conj!(op::OuterProduct) = OuterProduct(conj!(op.left), conj!(op.right))
Base.transpose(op::OuterProduct) = OuterProduct(op.right, op.left)
Base.adjoint(op::OuterProduct) = OuterProduct(conj(op.right), conj(op.left))

outerprod(left::AbstractVector, right::AbstractVector) = OuterProduct(left, right)
outerprod(left::AbstractMatrix, right::AbstractMatrix) = OuterProduct(left, right)
outerprod(outδ::ArrayReg{1}, in::ArrayReg{1}) = outerprod(statevec(outδ), conj(statevec(in)))
outerprod(outδ::ArrayReg{B}, in::ArrayReg{B}) where {B} =
    outerprod(statevec(outδ), conj(statevec(in)))

"""
    projection(y::AbstractMatrix, op::AbstractMatrix) -> typeof(y)

Project `op` to sparse matrix with same sparsity as `y`.
"""
function projection(y::AbstractMatrix, op::AbstractMatrix)
    size(y) == size(op) || throw(DimensionMismatch("can not project a matrix of size $(size(op)) to target size $(size(y))"))
    out = zero(y)
    unsafe_projection!(out, op)
end

unsafe_projection!(y::Diagonal, m::AbstractMatrix) = (y.diag .= diag(m); y)
unsafe_projection!(y::Diagonal, op::OuterProduct) = (y.diag .+= op.left .* op.right; y)

unsafe_projection!(y::Matrix, adjy, v) = y .+= adjy .* v

@inline function unsafe_projection!(y::AbstractSparseMatrix, m::AbstractMatrix)
    is, js, vs = findnz(y)
    for (k, (i, j)) in enumerate(zip(is, js))
        @inbounds y.nzval[k] += m[i, j]
    end
    y
end

@inline function unsafe_projection!(y::PermMatrix, m::AbstractMatrix)
    for i in 1:size(y, 1)
        @inbounds y.vals[i] = m[i, y.perm[i]]
    end
    y
end

projection(x::RT, adjx::Complex) where {RT<:Real} = RT(real(adjx))
projection(x::T, y::T) where {T} = y
projection(x::T1, y::T2) where {T1,T2} = convert(T1, y)


text:
"""
    abstract type SetMapBridge{T,S2,S1,F,G} <: AbstractBridge end

Consider two type of sets `S1`, `S2` and a linear mapping `A` that
the image of a set of type `S1` under `A` is a set of type `S2`.
A `SetMapBridge{T,S2,S1,F,G}` is a bridge that maps `G`-in-`S2` constraints
into `F`-in-`S1` by mapping the function through `A`.

The linear map `A` is described by
[`MathOptInterface.Bridges.map_set`](@ref),
[`MathOptInterface.Bridges.map_function`](@ref).
Implementing a method for these two functions is sufficient to bridge
constraints. In order for the getters and setters of dual solutions,
starting values, etc...  to work as well a method for the following
functions should be implemented as well:
[`MathOptInterface.Bridges.inverse_map_set`](@ref),
[`MathOptInterface.Bridges.inverse_map_function`](@ref),
[`MathOptInterface.Bridges.adjoint_map_function`](@ref) and
[`MathOptInterface.Bridges.inverse_adjoint_map_function`](@ref).
See the docstrings of the function
to see which feature would be missing it it was not implemented for a given
bridge.
"""
abstract type SetMapBridge{T,S2,S1,F,G} <: AbstractBridge end

function bridge_constraint(
    BT::Type{<:SetMapBridge{T,S2,S1,F,G}},
    model::MOI.ModelLike,
    func::G,
    set::S1,
) where {T,S2,S1,F,G}
    mapped_func = MOIB.map_function(BT, func)
    constraint = MOI.add_constraint(model, mapped_func, MOIB.map_set(BT, set))
    return BT(constraint)
end

function MOI.supports_constraint(
    ::Type{<:SetMapBridge{T,S2,S1}},
    ::Type{<:MOI.AbstractScalarFunction},
    ::Type{S1},
) where {T,S2,S1<:MOI.AbstractScalarSet}
    return true
end

function MOI.supports_constraint(
    ::Type{<:SetMapBridge{T,S2,S1}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{S1},
) where {T,S2,S1<:MOI.AbstractVectorSet}
    return true
end

function MOIB.added_constrained_variable_types(::Type{<:SetMapBridge})
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(
    ::Type{<:SetMapBridge{T,S2,S1,F}},
) where {T,S2,S1,F}
    return Tuple{Type,Type}[(F, S2)]
end

# Attributes, Bridge acting as a model
function MOI.get(
    ::SetMapBridge{T,S2,S1,F},
    ::MOI.NumberOfConstraints{F,S2},
)::Int64 where {T,S2,S1,F}
    return 1
end

function MOI.get(
    bridge::SetMapBridge{T,S2,S1,F},
    ::MOI.ListOfConstraintIndices{F,S2},
) where {T,S2,S1,F}
    return [bridge.constraint]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::SetMapBridge)
    MOI.delete(model, bridge.constraint)
    return
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::SetMapBridge{T,S2,S1,F,G},
) where {T,S2,S1,F,G}
    mapped_func = MOI.get(model, attr, bridge.constraint)
    func = MOIB.inverse_map_function(typeof(bridge), mapped_func)
    return MOIU.convert_approx(G, func)
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::SetMapBridge{T,S2,S1,F,G},
    func::G,
) where {T,S2,S1,F,G}
    MOI.set(
        model,
        attr,
        bridge.constraint,
        MOIB.map_function(typeof(bridge), func),
    )
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::SetMapBridge,
)
    set = MOI.get(model, attr, bridge.constraint)
    return MOIB.inverse_map_set(typeof(bridge), set)
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::SetMapBridge{T,S2,S1},
    new_set::S1,
) where {T,S2,S1}
    MOI.set(
        model,
        attr,
        bridge.constraint,
        MOIB.map_set(typeof(bridge), new_set),
    )
    return
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{<:SetMapBridge{T,S2,S1,F}},
) where {T,S2,S1,F}
    return MOI.supports(model, attr, MOI.ConstraintIndex{F,S2})
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::SetMapBridge,
)
    value = MOI.get(model, attr, bridge.constraint)
    return MOIB.inverse_map_function(typeof(bridge), value)
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::SetMapBridge,
    value,
)
    mapped_value = MOIB.map_function(typeof(bridge), value)
    MOI.set(model, attr, bridge.constraint, mapped_value)
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::SetMapBridge,
)
    value = MOI.get(model, attr, bridge.constraint)
    return MOIB.adjoint_map_function(typeof(bridge), value)
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDualStart,
    bridge::SetMapBridge,
    value,
)
    mapped_value = MOIB.inverse_adjoint_map_function(typeof(bridge), value)
    MOI.set(model, attr, bridge.constraint, mapped_value)
    return
end

# By linearity of the map, we can just change the constant/coefficient
function _map_change(::Type{BT}, change::MOI.ScalarConstantChange) where {BT}
    constant = MOIB.map_function(BT, change.new_constant)
    return MOI.ScalarConstantChange(constant)
end
function _map_change(::Type{BT}, change::MOI.VectorConstantChange) where {BT}
    constant = MOIB.map_function(BT, change.new_constant)
    return MOI.VectorConstantChange(constant)
end
function _map_change(::Type{BT}, change::MOI.ScalarCoefficientChange) where {BT}
    coefficient = MOIB.map_function(BT, change.new_coefficient)
    return MOI.ScalarCoefficientChange(change.variable, coefficient)
end
function _map_change(::Type{BT}, change::MOI.MultirowChange) where {BT}
    # It is important here that `change.new_coefficients` contains
    # the complete new sparse column associated to the variable.
    # Calling modify twice with part of the column won't work since
    # the linear map might reset all the column each time.
    coefficients = MOIB.map_function(BT, change.new_coefficients)
    return MOI.MultirowChange(change.variable, coefficients)
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::SetMapBridge,
    change::MOI.AbstractFunctionModification,
)
    MOI.modify(model, bridge.constraint, _map_change(typeof(bridge), change))
    return
end

include("flip_sign.jl")
const GreaterToLess{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{GreaterToLessBridge{T},OT}
const LessToGreater{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{LessToGreaterBridge{T},OT}
const NonnegToNonpos{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NonnegToNonposBridge{T},OT}
const NonposToNonneg{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{NonposToNonnegBridge{T},OT}

include("soc_rsoc.jl")
const RSOC{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{RSOCtoSOCBridge{T},OT}
const SOCR{T,OT<:MOI.ModelLike} = SingleBridgeOptimizer{SOCtoRSOCBridge{T},OT}

include("ltgt_to_interval.jl")
const GreaterToInterval{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{GreaterToIntervalBridge{T},OT}
const LessToInterval{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{LessToIntervalBridge{T},OT}

include("soc_to_psd.jl")
const SOCtoPSD{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{SOCtoPSDBridge{T},OT}
const RSOCtoPSD{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{RSOCtoPSDBridge{T},OT}


text:
include("data.jl")
using Flux, Statistics
using Flux: onehot, onehotbatch, throttle, crossentropy, reset!, onecold

const epochs = 20

train = gendata(100, 2)
val = gendata(10, 2)

scanner = LSTM(length(alphabet), 20)
encoder = Dense(20, length(alphabet))

function model(x)
    state = scanner.(x.data)[end]
    reset!(scanner)
    softmax(encoder(state))
end

loss(x, y) = crossentropy(model(x), y)
batch_loss(data) = mean(loss(d...) for d in data)

opt = ADAM(params(scanner, encoder))
evalcb = () -> @show batch_loss(val)

for i=1:epochs
    Flux.train!(loss, train, opt, cb=throttle(evalcb, 10))
end

# sanity test
tx = map(c -> onehotbatch(c, alphabet), [
    [false, true], # 01 -> 1
    [true, false], # 10 -> 1
    [false, false], # 00 -> 0
    [true, true]]) # 11 -> 0
[onecold(model(x)) - 1 for x in tx] |> println

max_stars_repo_name: lcy-seso/model-zoo

max_issues_repo_name: lcy-seso/model-zoo

max_forks_repo_name: lcy-seso/model-zoo


text:
using MathOptInterface
using Test

const MOI = MathOptInterface

@testset "MOI.FileFormats tests" begin
    @testset "$(file)" for file in ["CBF", "LP", "MOF", "MPS", "SDPA"]
        include(joinpath(@__DIR__, file, "$(file).jl"))
    end

    @testset "Copying options" begin
        models = [
            MOI.FileFormats.CBF.Model,
            MOI.FileFormats.LP.Model,
            MOI.FileFormats.MOF.Model,
            MOI.FileFormats.MPS.Model,
            MOI.FileFormats.SDPA.Model,
        ]
        for src in models
            model_src = src()
            for dest in models
                model_dest = dest()
                MOI.copy_to(model_dest, model_src)
                @test !isempty(sprint(write, model_dest))
            end
            model_dest = MOI.Utilities.MockOptimizer(
                MOI.Utilities.Model{Float64}()
            )
            MOI.copy_to(model_dest, model_src)
        end
    end

    @testset "Calling MOF.compressed_open" begin
        for cs in [MOI.FileFormats.Bzip2(), MOI.FileFormats.Gzip()]
            for open_type in ["a", "r+", "w+", "a+"]
                @test_throws ArgumentError MOI.FileFormats.compressed_open(
                    (x) -> nothing, "dummy.gz", open_type, cs
                )
            end
        end
    end

    @testset "Provided compression schemes" begin
        filename = joinpath(@__DIR__, "MPS", "free_integer.mps")
        model = MOI.FileFormats.Model(filename = filename)
        MOI.read_from_file(model, filename)
        filename = joinpath(@__DIR__, "free_integer.mps")
        MOI.write_to_file(model, filename * ".garbage")
        for ext in ["", ".bz2", ".gz"]
            MOI.write_to_file(model, filename * ext)
            model2 = MOI.FileFormats.Model(filename = filename * ext)
            MOI.read_from_file(model2, filename)
        end

        sleep(1.0)  # Allow time for unlink to happen.
        for ext in ["", ".garbage", ".bz2", ".gz"]
            rm(filename * ext, force = true)
        end
    end

    @testset "Model" begin
        for (format, model) in [
            (MOI.FileFormats.FORMAT_CBF, MOI.FileFormats.CBF.Model()),
            (MOI.FileFormats.FORMAT_LP, MOI.FileFormats.LP.Model()),
            (MOI.FileFormats.FORMAT_MOF, MOI.FileFormats.MOF.Model()),
            (MOI.FileFormats.FORMAT_MPS, MOI.FileFormats.MPS.Model()),
            (MOI.FileFormats.FORMAT_SDPA, MOI.FileFormats.SDPA.Model()),
        ]
            @test typeof(
                MOI.FileFormats.Model(format = format, filename = "foo.bar")
            ) == typeof(model)
        end
        @test_throws(
            ErrorException(
                "When `format==FORMAT_AUTOMATIC` you must pass a `filename`."
            ),
            MOI.FileFormats.Model(format = MOI.FileFormats.FORMAT_AUTOMATIC)
        )
        for (ext, model) in [
            (".cbf", MOI.FileFormats.CBF.Model()),
            (".lp", MOI.FileFormats.LP.Model()),
            (".mof.json", MOI.FileFormats.MOF.Model()),
            (".mps", MOI.FileFormats.MPS.Model()),
            (".sdpa", MOI.FileFormats.SDPA.Model()),
        ]
            @test typeof(MOI.FileFormats.Model(filename = "a$(ext)")) ==
                typeof(model)
            @test typeof(MOI.FileFormats.Model(filename = "a$(ext).gz")) ==
                typeof(model)
        end
        @test_throws(
            ErrorException("Unable to automatically detect format of a.b."),
            MOI.FileFormats.Model(filename = "a.b")
        )
    end
end


text:
using Test
using MiniTB

include( "../src/taxcalcs.jl" )
include( "../src/piecewise_linear_generator.jl" )

# @testset begin
 p1 = Point2D( 1.011,2.011)
 p2 = Point2D( 1.011, 2.011 )
 p3 = Point2D( 1.011, 3.011 )
 p4 = Point2D( -1.011, -3.011 )
 p5 = Point2D( 10.011, 13.011 )

 bc = BudgetConstraint( [p1,p2,p3,p4,p5])

sort!( bc )

@test size( bc )[1] == 5

ps = PointsSet([p1,p2,p3,p4,p5] )
push!( ps, p1 )
push!( ps, p2 )
push!( ps, p3 )
push!( ps, p4 )
push!( ps, p5 )
println( "ps = $ps ")
push!( ps, Point2D( 1.0011, 1.999999 ))

bc2 = censor( ps )

println( "bc2 = $bc2 ")
@test size( bc2 )[1] == 5
pers = DEFAULT_PERSON
pars = deepcopy(DEFAULT_PARAMS)
res = calculate( pers, pars )
println( "res=$res" )

# @test size( ps ) == 4
function getnet( data :: Dict, gross :: Float64 ) :: Float64
     person = data[:person]
     person.wage = gross
     person.hours = gross/DEFAULT_WAGE
     rc = calculate( person, data[:params] )
     return rc[:netincome]
end

function makebc( person :: Person, params :: TBParameters ) :: BudgetConstraint
    data = Dict(
        :person=>person,
        :params=>params )
    bc = makebc( data, getnet )
    return bc
end

bc = makebc( DEFAULT_PERSON, DEFAULT_PARAMS )
println( "\nDefault Case\n")
println(pointstoarray( bc))
annotations = annotate_bc( bc )
println( annotations )

bc = makebc( DEFAULT_PERSON, ZERO_PARAMS )
println( "\nZero Params\n")
println( pointstoarray( bc))
annotations = annotate_bc( bc )
println( annotations )

person = DEFAULT_PERSON
# rc = calculate( DEFAULT_PERSON, pars )
pars999 = deepcopy(DEFAULT_PARAMS)
pars999.it_allow = 999.0

println( "\n999 allowance case\n")
bc = makebc( person, pars999 )
println( pointstoarray( bc ))
annotations = annotate_bc( bc )
println( annotations )


text:
export defaults

"Dictionary of default values for parameters."
const defaults = Dict{String, Any}(
    "out_dir"                => "exprose_out",
    "n_strucs"               => 250,
    "tolerance_weight"       => 0.3,
    "other_ratio"            => 20.0,
    "out_prefix"             => "out",
    "iters_per_atom"         => 60000,
    "discard_ratio"          => 1.5,
    "box_size"               => 100.0,
    "cyc_step_ratio"         => 50.0,
    "ens_align_cutoff"       => 0.001,
    "align_cycles"           => 5,
    "align_cutoff"           => 0.5,
    "mod_n_points"           => 120,
    "mod_intra_cutoff"       => 0.0,
    "mod_intra_tolerance"    => 0.1,
    "mod_inter_cutoff"       => 7.0,
    "mod_inter_tolerance"    => 0.1,
    "mod_min_constraint_dist"=> 0.001,
    "out_dir_param"          => "exprose_param",
    "n_strucs_param"         => 10,
    "frac_between"           => 0.9,
    "tmscore_path"           => "TMscore",
    "param_tw_start"         => 1.0,
    "param_tw_increment"     => 0.1,
    "cluster_dist"           => 3.0,
    "exprose_repo_url"       => "https://github.com/jgreener64/ProteinEnsembles.jl",
)


text:
function sequencesToMatrix(samples::Int64,l::Int64,sequences::Array{Tuple{String,String},1})

	matrix = Array{Char}(undef,samples + 2,l)

	deleteIndex = Array{Int8}[]

	@simd for i in 1:samples
		# Extract each sample sequence
		# tmp = multiFasta[samples[i]][:].seq
		tmp = sequences[i][2]
		@inbounds matrix[i+1,:] = collect(tmp)
		# if(len(tmp) != seqLen)
		# 	print('errorAlign')
		# 	break
		# if('N' in tmp)
		# 	deleteIndex.append(i)
		# else
	end
	# matrix = np.delete(matrix,deleteIndex,0)
	return matrix
end

function degeneracy(data::String,codonDict::String)
	#degeneracy DICTIONARIES
	standardDict = Dict{String,String}(
		"TTT"=> "002", "TTC"=> "002", "TTA"=> "202", "TTG"=> "202",
		"TCT"=> "004", "TCC"=> "004", "TCA"=> "004", "TCG"=> "004",
		"TAT"=> "002", "TAC"=> "002", "TAA"=> "022", "TAG"=> "002",
		"TGT"=> "002", "TGC"=> "002", "TGA"=> "020", "TGG"=> "000",
		"CTT"=> "004", "CTC"=> "004", "CTA"=> "204", "CTG"=> "204",
		"CCT"=> "004", "CCC"=> "004", "CCA"=> "004", "CCG"=> "004",
		"CAT"=> "002", "CAC"=> "002", "CAA"=> "002", "CAG"=> "002",
		"CGT"=> "004", "CGC"=> "004", "CGA"=> "204", "CGG"=> "204",
		"ATT"=> "003", "ATC"=> "003", "ATA"=> "003", "ATG"=> "000",
		"ACT"=> "004", "ACC"=> "004", "ACA"=> "004", "ACG"=> "004",
		"AAT"=> "002", "AAC"=> "002", "AAA"=> "002", "AAG"=> "002",
		"AGT"=> "002", "AGC"=> "002", "AGA"=> "202", "AGG"=> "202",
		"GTT"=> "004", "GTC"=> "004", "GTA"=> "004", "GTG"=> "004",
		"GCT"=> "004", "GCC"=> "004", "GCA"=> "004", "GCG"=> "004",
		"GAT"=> "002", "GAC"=> "002", "GAA"=> "002", "GAG"=> "002",
		"GGT"=> "004", "GGC"=> "004", "GGA"=> "004", "GGG"=> "004")

	vectorIter = collect(1:3:length(data));
	degen = Array{String}(undef,size(vectorIter,1))
	for i in 1:size(vectorIter,1)
		codon = data[vectorIter[i]:vectorIter[i]+2]
		if ('N' in codon) || ('-' in codon)
			degen = string(degen,codon)
		else
		degen[i] =  standardDict[codon]
		end
	end
	degen = join(degen)
	return degen
end

function iterFastaMatrix(sequenceMatrix::Array{Char,2})

	# output = DataFrame([Float64,Int64,String],[:daf,:div,:degen])
	output = Array{Any,2}[]

	for n in eachcol(sequenceMatrix)

		degen   = n[1]
		aa      = n[end]
		pol     = n[2:end-1]

		# Undefined Ancestra Allele. Try to clean out of the loop
		if (aa == 'N' || aa == '-')
			continue
		# Monomorphic sites. Try to clean out of the loop
		elseif (size(unique(n[2:end][n[2:end] .!= 'N']),1) == 1)
			continue
		else

			if degen == '4'
				functionalClass = "4fold"
			else
				functionalClass = "0fold"
			end

			# Check if pol != AA and monomorphic
			if (size(unique(pol),1) == 1 && unique(pol)!= aa)
					div = 1.0; af = 0.0
					tmp = [af div functionalClass]
					push!(output,tmp);
			else
					div    = 0.0
					an     = size(pol,1)
					ac     = countmap(pol)

					if (aa ∉ keys(ac))
						continue
					else
						da = collect(keys(ac))
						da = da[da .!= aa ]

						if (size(da,1) < 2)
							af = ac[da[1]]/an
						else
							continue
						end
					end
					tmp = [af div functionalClass]
					# println(tmp);
					push!(output,tmp);
			end
		end
	end

	return reduce(vcat,output)
end

function formatSfs(rawSfsOutput::Array{Any,2},samples::Int64,bins::Int64)

	freq = OrderedDict(round.(collect(1:samples)/samples,digits=5) .=> 0)

	sfs = rawSfsOutput[rawSfsOutput[:,2].!=1,[1,3]]
	pn = sort!(OrderedDict(StatsBase.countmap(sfs[sfs[:,2] .!= "4fold",1])))
	ps = sort!(OrderedDict(StatsBase.countmap(sfs[sfs[:,2] .== "4fold",1])))

	sfsPn        = reduce(vcat,values(merge(+,freq,pn)))
	sfsPs        = reduce(vcat,values(merge(+,freq,ps)))
	sfs::DataFrame = DataFrame(reduceSfs(hcat(freq.keys,sfsPn,sfsPs),bins),[:f,:pi,:p0])

	divergence = rawSfsOutput[rawSfsOutput[:,2].==1,2:3]

	d::DataFrame = DataFrame(countmap(divergence[:,2]))
	return(sfs,d)
end

function uSfsFromFasta(;file::String,reference::String,outgroup::String,samples::Int64,bins::Int64,codonTable::String)

	multiFasta = readfasta(file);
	ref        = readfasta(reference);
	out        = readfasta(outgroup);

    samples    = size(multiFasta)[1]
    seqLen     = length(ref[1][2])

	multiFastaMatrix = sequencesToMatrix(samples,seqLen,multiFasta);

	degenCode = collect(degeneracy(ref[1][2],codonTable));
	multiFastaMatrix[1,:] = degenCode
	multiFastaMatrix[end,:] = collect(out[1][2])

	m = multiFastaMatrix[:,(multiFastaMatrix[end,:].!= 'N') .& (multiFastaMatrix[end,:].!= '-')]
	m = m[:,(m[1,:] .== '4') .| (m[1,:] .=='0')]

	joinedSfsDiv = iterFastaMatrix(m)

	sfs, d = formatSfs(joinedSfsDiv,samples,bins)
	return sfs,d
end

function uSfsFromVcf(file::String)

	reader = VCF.Reader(GzipDecompressorStream(open(file, "r")))
	output = []

	for record in reader
		ref =  VCF.ref(record)
		alt = VCF.alt(record)
		info = Dict(VCF.info(record))
		aa = info["AA"]
		ac = parse(Int64,info["AC"]);
		an = parse(Int64,info["AN"]);
		af = ac/an

		if (aa == ref)
			daf = af
			tmp = [parse(Int64,VCF.chrom(record)),VCF.pos(record),daf]
			push!(output,tmp)
		elseif (aa == alt)
			daf = 1 - af
			tmp = [parse(Int64,VCF.chrom(record)),VCF.pos(record),daf]
			push!(output,tmp)
		else
			continue
		end
	end
	close(reader)
	return output
end


text:
include("../utility.jl")
using Turing, Test

N = 10
setchunksize(N)
# alg = HMCDA(2000, 0.65, 1.5)
alg = HMC(1000, 0.2, 4)

@model vdemo() = begin
  x = Vector{Real}(undef, N)
  for i = 1:N
    x[i] ~ Normal(0, sqrt(4))
  end
end

t_loop = @elapsed res = sample(vdemo(), alg)

# Test for vectorize UnivariateDistribution
@model vdemo() = begin
  x = Vector{Real}(undef, N)
  x ~ [Normal(0, 2)]
end

t_vec = @elapsed res = sample(vdemo(), alg)

@model vdemo() = begin
  x ~ MvNormal(zeros(N), 2 * ones(N))
end

t_mv = @elapsed res = sample(vdemo(), alg)

println("Time for")
println("  Loop : $t_loop")
println("  Vec  : $t_vec")
println("  Mv   : $t_mv")

# Transformed test
@model vdemo() = begin
  x = Vector{Real}(undef, N)
  x ~ [InverseGamma(2, 3)]
end

sample(vdemo(), alg)


text:
using Dates
# using GitCommand
using GitHub
using JSON
using Pkg
using Printf
using RegistryCI
using Test
using TimeZones

@testset "Public interface" begin
    @testset "RegistryCI.test" begin
        path = joinpath(DEPOT_PATH[1], "registries", "General")
        RegistryCI.test(path)
    end
end

@testset "Internal functions (private)" begin
    @testset "RegistryCI.load_registry_dep_uuids" begin
        all_registry_deps_names = [
            ["General"],
            ["https://github.com/JuliaRegistries/General"],
            ["https://github.com/JuliaRegistries/General.git"],
        ]
        for registry_deps_names in all_registry_deps_names
            extrauuids = RegistryCI.load_registry_dep_uuids(registry_deps_names)
            @test extrauuids isa Set{Base.UUID}
            @test length(extrauuids) > 1_000
        end
    end
end


text:
import Base
import REPL

have_color = false
is_interactive = false

"exec_options(opts): Adapted from base/client.jl"
function exec_options(opts)
    quiet                 = (opts.quiet != 0)
    startup               = (opts.startupfile != 2)
    history_file          = (opts.historyfile != 0)
    color_set             = (opts.color != 0) # --color!=auto
    global have_color     = (opts.color == 1) # --color=on
    global is_interactive = (opts.isinteractive != 0)
    interactiveinput = isa(stdin, Base.TTY)
    if interactiveinput
        global is_interactive = true
        banner = (opts.banner != 0) # --banner!=no
    else
        banner = (opts.banner == 1) # --banner=yes
    end
    run_main_repl(interactiveinput, quiet, banner, history_file, color_set)
    nothing
end
"exec_options(): Use Base.JLOptions())"
exec_options() = exec_options(Base.JLOptions())

"run_main_repl(interactive, quiet, banner, history_file, color_set): See base/client.jl"
function run_main_repl(interactive::Bool, quiet::Bool, banner::Bool, history_file::Bool, color_set::Bool)
    # load interactive-only libraries
    if !isdefined(Main, :InteractiveUtils)
        try
            let InteractiveUtils = Base.require(Base.PkgId(Base.UUID(0xb77e0a4c_d291_57a0_90e8_8db25a27a240), "InteractiveUtils"))
                Core.eval(Main, :(const InteractiveUtils = $InteractiveUtils))
                Core.eval(Main, :(using .InteractiveUtils))
            end
        catch ex
            @warn "Failed to import InteractiveUtils into module Main" exception=(ex, catch_backtrace())
        end
    end

    if interactive && isassigned(Base.REPL_MODULE_REF)
        Base.invokelatest(Base.REPL_MODULE_REF[]) do REPL
            term_env = get(ENV, "TERM", @static Sys.iswindows() ? "" : "dumb")
            term = REPL.Terminals.TTYTerminal(term_env, stdin, stdout, stderr)
            color_set || (global have_color = REPL.Terminals.hascolor(term))
            set_have_color()
            banner && Base.banner(term)
            if term.term_type == "dumb"
                active_repl = REPL.BasicREPL(term)
                quiet || @warn "Terminal not fully functional"
            else
                active_repl = REPL.LineEditREPL(term, have_color, true)
                active_repl.history_file = history_file
            end
            set_active_repl(active_repl)
            # Make sure any displays pushed in .julia/config/startup.jl ends up above the
            # REPLDisplay
            Base.pushdisplay(REPL.REPLDisplay(active_repl))
            Base._atreplinit(active_repl)
            if Base.roottask === Base.current_task()
                active_repl.interface.modes[1].prompt = "julia-root> "
            end
            set_active_repl(active_repl)
            run_repl(active_repl, set_active_repl_backend)
        end
    end
end
run_main_repl() = exec_options()

function set_active_repl(repl::REPL.AbstractREPL)
    @eval Base active_repl = $repl
    @eval Base have_color = $have_color
    @eval Base is_interactive = $is_interactive
end

function set_have_color()
    @eval Base have_color = $have_color
end

function set_active_repl_backend(backend::REPL.REPLBackend)
    @eval Base active_repl_backend = $backend
end

"run_repl(repl, consumer): Run frontend on another Task, run backend on current Task. See stdlib/REPL/src/REPL.jl"
function run_repl(repl::REPL.AbstractREPL, @nospecialize(consumer = x -> nothing))
    repl_channel = Channel(1)
    response_channel = Channel(1)
    backend = REPL.REPLBackend(repl_channel, response_channel, false)
    backend.backend_task = Base.current_task()
    consumer(backend)
    @debug "RootTaskREPL: Starting frontend_task"
    frontend_task = @async begin
            REPL.run_frontend(repl, REPL.REPLBackendRef(repl_channel, response_channel))
            @debug "RootTaskREPL: Frontend exit"
            # Explicitly stop the backend REPL
            put!(repl_channel, (nothing, -1))
    end
    @debug "RootTaskREPL: Starting backend"
    repl_backend_loop(backend)
    @debug "RootTaskREPL: Backend exit"
    # Make non-interactive so normal REPL does not start
    @eval Base is_interactive = false
    return backend
end

"repl_backend_loop(backend): Run backend loop on current Task rather than using @async. See stdlib/REPL/src/REPL.jl"
function repl_backend_loop(backend::REPL.REPLBackend)
    # include looks at this to determine the relative include path
    # nothing means cwd
    while true
        tls = task_local_storage()
        tls[:SOURCE_PATH] = nothing
        ast, show_value = take!(backend.repl_channel)
        if show_value == -1
            # exit flag
            break
        end
        REPL.eval_user_input(ast, backend)
    end
    nothing
end


text:
struct Difference{PA, PB}
    a::PA
    b::PB
end

Base.:*(n::Number, d::Difference) = Difference(*(n, d.a), *(n, d.b))
Base.:+(d1::Difference, d2::Difference) = Difference(d1.a + d2.a, d1.b + d2.b)

struct CollisionCache{GeomA, GeomB, M, D <: Difference}
    bodyA::GeomA
    bodyB::GeomB
    simplex_points::MVector{M, D}
end

function CollisionCache(geomA, geomB)
    N = dimension(typeof(geomA))
    @assert dimension(typeof(geomA)) == dimension(typeof(geomB))
    CollisionCache(N, geomA, geomB)
end

@generated function initial_simplex(dim::Val{N}, geomA, geomB) where {N}
    return quote
        interior_point = Difference(any_inside(geomA), any_inside(geomB))
        simplex = $(Expr(:call, :(MVector), [:(interior_point) for i in 1:(N + 1)]...))
    end
end

function CollisionCache(::Val{N}, geomA, geomB) where N
    simplex = initial_simplex(Val(N), geomA, geomB)
    CollisionCache(geomA, geomB, simplex)
end

dimension(::Type{CollisionCache{G1, G2, M, D}}) where {G1,G2,M,D} = dimension(G1)

function support_vector_max(geometry::gt.GeometryPrimitive, direction, initial_guess::Tagged)
    best_pt, score = gt.support_vector_max(geometry, gt.Vec(direction))
    Tagged(SVector(best_pt))
end

function support_vector_max(pt::gt.Vec{N, T}, direction, initial_guess::Tagged) where {N,T}
    Tagged(SVector(pt))
end

function support_vector_max(simplex::Union{gt.AbstractSimplex, gt.AbstractFlexibleGeometry}, direction, initial_guess::Tagged)
    best_pt, score = gt.support_vector_max(simplex, gt.Vec(direction))
    Tagged(SVector(best_pt))
end

function support_vector_max(mesh::gt.HomogenousMesh{gt.Point{N, T}}, direction, initial_guess::Tagged) where {N,T}
    best_arg, best_value = gt.argmax(x-> dot(SVector(x), direction), gt.vertices(mesh))
    best_vec = SVector(best_arg)
    Tagged(best_vec)
end

any_inside(pt::SVector) = Tagged(pt)
support_vector_max(pt::SVector, direction, initial_guess::Tagged) = Tagged(pt)

function transform_simplex(cache::CollisionCache, poseA, poseB)
    transform_simplex(dimension(typeof(cache)), cache, poseA, poseB)
end

@generated function transform_simplex(::Val{N}, cache::CollisionCache, poseA, poseB) where {N}
    transform_simplex_impl(N, cache, poseA, poseB)
end

function transform_simplex_impl(N, cache, poseA, poseB)
    Expr(:call, :(SVector),
        [:((poseA(value(cache.simplex_points[$i].a)) -
            poseB(value(cache.simplex_points[$i].b)))) for i in 1:(N + 1)]...)
end

struct GJKResult{M, N, T}
    simplex::SVector{M, SVector{N, T}}
    closest_point_in_body::Difference{SVector{N, T}, SVector{N, T}}
    signed_distance::T
end

function gjk!(cache::CollisionCache,
              poseA::Transformation,
              poseB::Transformation,
              max_iter=100,
              atol=1e-6)
    rotAinv = transform_deriv(inv(poseA), 0)
    rotBinv = transform_deriv(inv(poseB), 0)
    simplex = transform_simplex(cache, poseA, poseB)
    iter = 1

    while true
        weights = projection_weights(simplex)
        min_weight, index_to_replace = findmin(weights)
        if min_weight > 0
            # in collision
            return GJKResult(
                simplex,
                dot(weights, cache.simplex_points),
                penetration_distance(simplex)
            )
        end
        best_point = dot(weights, simplex)

        direction = -best_point
        direction_in_A = rotAinv * direction
        direction_in_B = rotBinv * direction

        starting_vertex_index = 1
        starting_vertex = cache.simplex_points[starting_vertex_index]
        starting_vertex_score =
            dot(value(starting_vertex.a), direction_in_A) +
            dot(value(starting_vertex.b), direction_in_B)
        for j in 2:length(cache.simplex_points)
            candidate = cache.simplex_points[j]
            candidate_score =
                dot(value(candidate.a), direction_in_A) +
                dot(value(candidate.b), direction_in_B)
            if candidate_score > starting_vertex_score
                starting_vertex_score = candidate_score
                starting_vertex = candidate
            end
        end

        improved_vertex = Difference(
            support_vector_max(cache.bodyA, direction_in_A, starting_vertex.a),
            support_vector_max(cache.bodyB, -direction_in_B, starting_vertex.b))
        improved_point = poseA(value(improved_vertex.a)) - poseB(value(improved_vertex.b))
        score = dot(improved_point, direction)
        if score <= dot(best_point, direction) + atol || iter >= max_iter
            return GJKResult(
                simplex,
                dot(weights, cache.simplex_points),
                norm(best_point)
            )
        else
            cache.simplex_points[index_to_replace] = improved_vertex
            simplex = setindex(simplex, improved_point, index_to_replace)
            # simplex[index_to_replace] = improved_point
        end
        iter += 1
    end
end

function penetration_distance(simplex)
    _, penetration_distance = gt.argmax(1:length(simplex)) do i
        face = simplex_face(simplex, i)
        weights = projection_weights(face)
        closest_point = dot(weights, face)
        distance_to_face = norm(closest_point)
        -distance_to_face
    end
    return penetration_distance
end

function gjk(geomA, geomB,
             poseA::Transformation=IdentityTransformation(),
             poseB::Transformation=IdentityTransformation())
    cache = CollisionCache(geomA, geomB)
    gjk!(cache, poseA, poseB)
end


text:
using Potree
using FileManager
using Visualization
using BenchmarkTools
using Profile

# potree source folder
# big files preferable
#potree = "C:/Users/Alessio/Documents/potreeDirectory/pointclouds/CAVA" # replace this path with local potree directory
potree = "C:/Users/Alessio/.julia/dev/Potree/potree/stairsLAS(v1.7)"
potreeBin = "C:/Users/Alessio/.julia/dev/Potree/potree/stairsBIN(v1.7)"

println("@btime for bin2pointcloud WITHOUT Multithreading")
@btime bin2pointcloudNoMultithreading($potreeBin)
println("@btime for bin2pointcloud WITH Multithreading")
@btime bin2pointcloud($potreeBin)
println("@btime for bin2pointcloud WITH NEW function")
@btime bin2pointcloudNEW($potreeBin)


text:
using Documenter, GRAYTR

push!(LOAD_PATH,"../src/")

makedocs(
    sitename = "GRAYTR"
)

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using CSV, Distributions, HypothesisTests

data1 = CSV.read("../data/machine1.csv", header=false)[:,1]
data2 = CSV.read("../data/machine2.csv", header=false)[:,1]
xBar1, xBar2 = mean(data1), mean(data2)
n1, n2 = length(data1), length(data2)
alpha = 0.05
t = quantile(TDist(n1+n2-2),1-alpha/2)

s1, s2 = std(data1), std(data2)
sP = sqrt(((n1-1)*s1^2 + (n2-1)*s2^2) / (n1+n2-2))

println("Calculating formula: ", (xBar1 - xBar2 - t*sP* sqrt(1/n1 + 1/n2),
				  xBar1 - xBar2 + t*sP* sqrt(1/n1 + 1/n2)))
println("Using confint(): ", confint(EqualVarianceTTest(data1,data2),alpha))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
"""
    AbstractTaylorModelReachSet{N}

Abstract type for all reach sets types that represent a Taylor model.

### Notes

The parameter `N` refers to the numerical type of the representation.
"""
abstract type AbstractTaylorModelReachSet{N} <: AbstractReachSet{N} end


text:
#
# hierarchicalsolve
#

# Continuous analogues for low-rank operators case

\{U<:Operator,V<:AbstractLowRankOperator}(H::HierarchicalOperator{U,V},f::Fun) = hierarchicalsolve(H,f)
\{U<:Operator,V<:AbstractLowRankOperator,F<:Fun}(H::HierarchicalOperator{U,V},f::Vector{F}) = hierarchicalsolve(H,f)

hierarchicalsolve{U<:Operator,V<:AbstractLowRankOperator}(H::HierarchicalOperator{U,V},f::Fun) = hierarchicalsolve(H,[f])[1]

hierarchicalsolve(H::Operator,f::Fun) = H\f
hierarchicalsolve{F<:Fun}(H::Operator,f::Vector{F}) = map(g->H\g,f)#vec(H\transpose(f))

function hierarchicalsolve{U<:Operator,V<:AbstractLowRankOperator,F<:Fun}(H::HierarchicalOperator{U,V},f::Vector{F})
    # Pre-compute Factorization

    !isfactored(H) && factorize!(H)

    # Partition HierarchicalOperator

    (H11,H22),(H21,H12) = partition(H)

    # Off-diagonal low-rank matrix assembly

    U12,V12 = H12.U,H12.V
    U21,V21 = H21.U,H21.V

    # Partition the right-hand side

    f1,f2 = partition(f,space(first(U12)),space(first(U21)))

    # Solve recursively

    H22f2,H11f1 = hierarchicalsolve(H22,f2),hierarchicalsolve(H11,f1)

    # Compute pivots

    v12,v21 = computepivots(V12,V21,H11f1,H22f2,H.factorization)

    # Solve again with updated right-hand sides

    RHS1 = f1 - v12.'*U12
    RHS2 = f2 - v21.'*U21

    sol1,sol2 = hierarchicalsolve(H11,RHS1),hierarchicalsolve(H22,RHS2)

    return assemble(sol1,sol2)
end

function factorize!{U<:Operator,V<:AbstractLowRankOperator}(H::HierarchicalOperator{U,V})
    # Partition HierarchicalOperator

    (H11,H22),(H21,H12) = partition(H)

    # Off-diagonal low-rank matrix assembly

    U12,V12 = H12.U,H12.V
    U21,V21 = H21.U,H21.V

    # Solve recursively

    H22U21,H11U12 = hierarchicalsolve(H22,U21),hierarchicalsolve(H11,U12)

    # Compute A

    fillpivotmatrix!(H.A,V12,V21,H22U21,H11U12)

    # Compute factorization

    H.factorization = pivotldufact(H.A,length(V12),length(V21))#lufact(H.A)
    H.factored = true
end

function fillpivotmatrix!{A1,A2,T}(A::Matrix{T},V12::Vector{Operator{T}},V21::Vector{Operator{T}},H22U21::Vector{Fun{A1,T}},H11U12::Vector{Fun{A2,T}})
    r1,r2 = length(V12),length(V21)
    for i=1:r1,j=1:r2
        @inbounds A[i,j+r1] += V12[i]*H22U21[j]
        @inbounds A[j+r2,i] += V21[j]*H11U12[i]
    end
end

function fillpivotmatrix!{V1,V2,A1,A2,T}(A::Matrix{T},V12::Vector{Fun{V1,T}},V21::Vector{Fun{V2,T}},H22U21::Vector{Fun{A1,T}},H11U12::Vector{Fun{A2,T}})
    r1,r2 = length(V12),length(V21)
    for i=1:r1,j=1:r2
        @inbounds A[i,j+r1] += linebilinearform(V12[i],H22U21[j])
        @inbounds A[j+r2,i] += linebilinearform(V21[j],H11U12[i])
    end
end

function computepivots{A1,A2,T}(V12::Vector{Operator{T}},V21::Vector{Operator{T}},H11f1::Vector{Fun{A1,T}},H22f2::Vector{Fun{A2,T}},A::PivotLDU{T})
    @assert length(H11f1) == length(H22f2)
    nf,r1,r2 = length(H11f1),length(V12),length(V21)
    b1,b2 = zeros(T,r1,nf),zeros(T,r2,nf)
    for i=1:nf
        for j=1:r1
            @inbounds b1[j,i] += V12[j]*H22f2[i]
        end
        for j=1:r2
            @inbounds b2[j,i] += V21[j]*H11f1[i]
        end
    end
    A_ldiv_B1B2!(A,b1,b2)
end

function computepivots{V1,V2,A1,A2,T}(V12::Vector{Fun{V1,T}},V21::Vector{Fun{V2,T}},H11f1::Vector{Fun{A1,T}},H22f2::Vector{Fun{A2,T}},A::PivotLDU{T})
    @assert length(H11f1) == length(H22f2)
    nf,r1,r2 = length(H11f1),length(V12),length(V21)
    b1,b2 = zeros(T,r1,nf),zeros(T,r2,nf)
    for i=1:nf
        for j=1:r1
            @inbounds b1[j,i] += linebilinearform(V12[j],H22f2[i])
        end
        for j=1:r2
            @inbounds b2[j,i] += linebilinearform(V21[j],H11f1[i])
        end
    end
    A_ldiv_B1B2!(A,b1,b2)
end

# Utilities

##
# partition(f::Fun,sp1,sp2) uses the spaces sp1 and sp2
# to achieve type-stability by specializing on the four cases:
#
# [x₁ | ⋯ ⋯ xₙ]
# [x₁ ⋯ | ⋯ xₙ]
# [x₁ ⋯ ⋯ | xₙ]
# [x₁ | x₂]
##

function partition{PWS<:PiecewiseSpace,S1<:PiecewiseSpace,S2<:PiecewiseSpace,T}(f::Fun{PWS,T},sp1::S1,sp2::S2)
    p,N1,N2 = pieces(f),length(sp1),length(sp2)
    return (depiece(p[1:N1])),(depiece(p[1+N1:N1+N2]))
end

function partition{PWS<:PiecewiseSpace,S1<:PiecewiseSpace,S2,T}(f::Fun{PWS,T},sp1::S1,sp2::S2)
    p,N1 = pieces(f),length(sp1)
    return (depiece(p[1:N1])),(p[1+N1])
end

function partition{PWS<:PiecewiseSpace,S1,S2<:PiecewiseSpace,T}(f::Fun{PWS,T},sp1::S1,sp2::S2)
    p,N2 = pieces(f),length(sp2)
    return (p[1]),(depiece(p[2:1+N2]))
end

function partition{PWS<:PiecewiseSpace,S1,S2,T}(f::Fun{PWS,T},sp1::S1,sp2::S2)
    p = pieces(f)
    return (p[1]),(p[2])
end

function partition{PWS<:PiecewiseSpace,S1,S2,T}(f::Vector{Fun{PWS,T}},sp1::S1,sp2::S2)
    p1 = partition(f[1],sp1,sp2)
    ret1,ret2 = Vector{typeof(p1[1])}(length(f)),Vector{typeof(p1[2])}(length(f))
    ret1[1],ret2[1] = p1[1],p1[2]
    for k=2:length(f)
        @inbounds ret1[k],ret2[k] = partition(f[k],sp1,sp2)
    end
    ret1,ret2
end

##
# assemble(sol1::Vector{Fun},sol2::Vector{Fun}) uses the spaces of sol1 and sol2
# to achieve type-stability by specializing on the four cases:
#
# [s₁ | ⋯ ⋯ sₙ]
# [s₁ ⋯ | ⋯ sₙ]
# [s₁ ⋯ ⋯ | sₙ]
# [s₁ | s₂]
##

function assemble{S1,S2,T1,T2}(sol1::Vector{Fun{S1,T1}},sol2::Vector{Fun{S2,T2}})
    @assert length(sol1) == length(sol2)
    p = depiece((sol1[1],sol2[1]))
    ret = Vector{typeof(p)}(length(sol1))
    ret[1] = p
    for k=2:length(sol1)
        @inbounds ret[k] = depiece((sol1[k],sol2[k]))
    end
    ret
end

function assemble{S1<:PiecewiseSpace,S2<:PiecewiseSpace,T1,T2}(sol1::Vector{Fun{S1,T1}},sol2::Vector{Fun{S2,T2}})
    @assert length(sol1) == length(sol2)
    p = depiece(vcat(pieces(sol1[1]),pieces(sol2[1])))
    ret = Vector{typeof(p)}(length(sol1))
    ret[1] = p
    for k=2:length(sol1)
        @inbounds ret[k] = depiece(vcat(pieces(sol1[k]),pieces(sol2[k])))
    end
    ret
end

function assemble{S1<:PiecewiseSpace,S2,T1,T2}(sol1::Vector{Fun{S1,T1}},sol2::Vector{Fun{S2,T2}})
    @assert length(sol1) == length(sol2)
    p = depiece(vcat(pieces(sol1[1]),sol2[1]))
    ret = Vector{typeof(p)}(length(sol1))
    ret[1] = p
    for k=2:length(sol1)
        @inbounds ret[k] = depiece(vcat(pieces(sol1[k]),sol2[k]))
    end
    ret
end

function assemble{S1,S2<:PiecewiseSpace,T1,T2}(sol1::Vector{Fun{S1,T1}},sol2::Vector{Fun{S2,T2}})
    @assert length(sol1) == length(sol2)
    p = depiece(vcat(sol1[1],pieces(sol2[1])))
    ret = Vector{typeof(p)}(length(sol1))
    ret[1] = p
    for k=2:length(sol1)
        @inbounds ret[k] = depiece(vcat(sol1[k],pieces(sol2[k])))
    end
    ret
end


text:
module VersusMatch

using Benchmarkplotting
using Statistics
using Gadfly
using MLStyle
using Match
import Base.getindex
getindex(asoc_lst :: Vector{Pair{Symbol, T}}, key ::Symbol) where T =
    for (search_key, value) in asoc_lst
        if search_key === key
            return value
        end
    end

data = [
    :node_fn1 => :(function f(a, b) a + b end),
    :node_fn2 => :(function f(a, b, c...) c end),
    :node_chain => :(subject.method(arg1, arg2)),
    :node_struct => :(
            struct name <: base
                field1 :: Int
                field2 :: Float32
            end
    ),
    :node_const => :(
            const a = value
    ),
    :node_assign => :(a = b + c)
]
implementations = [
    Symbol(:MLStyle, " Expr-pattern") => let
        function extract_name(e)
            MLStyle.@match e begin
                ::Symbol                           => e
                Expr(:<:, a, _)                    => extract_name(a)
                Expr(:struct, _, name, _)          => extract_name(name)
                Expr(:call, f, _...)               => extract_name(f)
                Expr(:., subject, attr, _...)      => extract_name(subject)
                Expr(:function, sig, _...)         => extract_name(sig)
                Expr(:const, assn, _...)           => extract_name(assn)
                Expr(:(=), fn, body, _...)         => extract_name(fn)
                Expr(expr_type,  _...)             => error("Can't extract name from ",
                                                            expr_type, " expression:\n",
                                                            "    $e\n")
            end
        end
        @assert extract_name(data[:node_fn1]) == :f
        @assert extract_name(data[:node_fn2]) == :f
        @assert extract_name(data[:node_chain]) == :subject
        @assert extract_name(data[:node_struct]) == :name
        @assert extract_name(data[:node_const]) == :a
        @assert extract_name(data[:node_assign]) == :a
        extract_name
    end,
    Symbol(:MLStyle, " AST-pattern") => let
        function extract_name_homoiconic(e)
            MLStyle.@match e begin
                ::Symbol                           => e
                :($a <: $_)                        => extract_name_homoiconic(a)
                :(struct $name <: $_
                    $(_...)
                end)                             => name
                :($f($(_...)))                     => extract_name_homoiconic(f)
                :($subject.$_)                     => extract_name_homoiconic(subject)
                :(function $name($(_...))
                    $(_...)
                end)                            => extract_name_homoiconic(name)
                :(const $assn = $_)                => extract_name_homoiconic(assn)
                :($fn = $_)                        => extract_name_homoiconic(fn)
                Expr(expr_type,  _...)             => error("Can't extract name from ",
                                                            expr_type, " expression:\n",
                                                            "    $e\n")
            end
        end
        @assert extract_name_homoiconic(data[:node_fn1]) == :f
        @assert extract_name_homoiconic(data[:node_fn2]) == :f
        @assert extract_name_homoiconic(data[:node_chain]) == :subject
        @assert extract_name_homoiconic(data[:node_struct]) == :name
        @assert extract_name_homoiconic(data[:node_const]) == :a
        @assert extract_name_homoiconic(data[:node_assign]) == :a
        extract_name_homoiconic
    end,
    Symbol("Match.jl") =>  let
        extract_name(e :: Symbol) = e
        function extract_name(e::Expr)
            Match.@match e begin
                Expr(:<:, [a, b])                  => extract_name(a)
                Expr(:struct,      [_, name, _])   => extract_name(name)
                Expr(:call,      [f, _...])        => extract_name(f)
                Expr(:., [subject, attr, _...])    => extract_name(subject)
                Expr(:function,  [sig, _...])      => extract_name(sig)
                Expr(:const,     [assn, _...])     => extract_name(assn)
                Expr(:(=),       [fn, body, _...]) => extract_name(fn)
                Expr(expr_type,  _...)             => error("Can't extract name from ",
                                                            expr_type, " expression:\n",
                                                            "    $e\n")
            end
        end
        @assert extract_name(data[:node_fn1]) == :f
        @assert extract_name(data[:node_fn2]) == :f
        @assert extract_name(data[:node_chain]) == :subject
        @assert extract_name(data[:node_struct]) == :name
        @assert extract_name(data[:node_const]) == :a
        @assert extract_name(data[:node_assign]) == :a
        extract_name
    end
]

criterion(x) = (meantime = mean(x.times), allocs = 1 + x.allocs)
df = bcompare(criterion, data, implementations)

theme = Theme(
    guide_title_position = :left,
    colorkey_swatch_shape = :circle,
    minor_label_font = "Consolas",
    major_label_font = "Consolas",
    point_size=6px
)
report_meantime, df_time = report(:meantime, df, Scale.y_log10, theme)
report_allocs, df_allocs = report(:allocs, df, theme)

open("stats/vs-match(expr).txt", "w") do f
    write(f, string(df))
end

draw(SVG("stats/vs-match(expr)-on-time.svg", 10inch, 4inch), report_meantime);
draw(SVG("stats/vs-match(expr)-on-allocs.svg", 10inch, 4inch), report_allocs);

end


text:
"""

The structure containing the set of neaerest-neighbor gates when
applied on the initial configuration `initconf` generates the
approximate MPS representation for a gaussian fermionic model. `lx` is
the size of the system. The left acting site of each gate is given in
the `positions` vector and each gate can be constructed using the
`θs` or possibly `ϕs`.

"""
struct GaussianMPS
    lx :: Int
    initconf :: Vector{Int}
    xs :: Vector{Int}
    θs :: Vector{Float64}

    function GaussianMPS(lx::Int,
                         initconf::Vector{Int},
                         xs::Vector{Int},
                         θs::Vector{Float64})
        @assert length(initconf) == lx
        @assert length(xs) == length(θs)

        new(lx, initconf, xs, θs)
    end
end

"""
    corrmat2gmps(corrmat, threshold, verbose)

For a given two-point correlation matrix `corrmat` called lambda
perform the approach explained in arxiv:1504.07701 and generate the
local (nearest-neighbor) gates.

"""

function corrmat2gmps(corrmat::Matrix{T},
                      threshold::Float64=1.e-8;
                      alternate::Bool=false) where {T}
    T <: Complex && error("correlation matrix is complex!")
    lx = size(corrmat)[1]

    #initconf = Vector{Int}(lx)
    initconf = zeros(Int, lx)
    positions = Int[]
    θs = Float64[]

    expected_next_evalue = 1
    for site=1:lx
        block_end = site
        evalue = corrmat[site, site]
        corr_block = [evalue]
        vals = corr_block

        if alternate
            delta = abs(expeted_next_evalue - evalue)
        else
            delta = min(abs(evalue), abs(1-evalue))
        end

        while (delta > threshold && block_end < lx)
            block_end += 1

            corr_block = corrmat[site:block_end, site:block_end]
            vals = eigvals(corr_block)

            if alternate
                delta = (expected_next_evalue == 0) ? abs.(minimum(vals)) : abs.(1-maximum(vals))
            else
                delta = min(abs.(minimum(vals)), abs.(1-maximum(vals)))
            end
        end

        if block_end > site
            block_size = block_end - site + 1
            block_range = site:block_end

            if alternate
                v_index = (expected_next_evalue == 0) ? indmin(abs.(vals)) : argmin(abs.(1-vals))
            else
                v_index = argmax(abs.(vals .- 0.5))
            end
            evalue = vals[v_index]
            v = eigvecs(corr_block)[:,v_index]

            # Find the set of (block_size-1) unitary gate that
            # diagonalize the current correlation block
            ugate_block = Matrix{T}(I, block_size, block_size)

            for p = block_size-1:-1:1 # p here stands for "pivot"

                θ = atan(v[p + 1]/ v[p])
                ugate = Matrix{Float64}(1.0I, block_size, block_size)
                ugate[p:p+1, p:p+1] =
                    [
                        cos(θ) -sin(θ);
                        sin(θ)  cos(θ)
                    ]
                push!(positions, site+p-1)
                push!(θs, θ)
                ## QQQ? acting on the vectors from left! Explain, why
                ## is this a better choice?
                v = transpose(transpose(v) * ugate)

                ugate_block = ugate_block * ugate
            end

            ugate_extended = Matrix{T}(I, lx, lx)
            ugate_extended[block_range, block_range] = ugate_block
            corrmat= Symmetric(ugate_extended' * corrmat * ugate_extended)
            # display(corrmat)
        end
        initconf[site] = round(evalue)
        expected_next_evalue = (round(evalue) == 1) ? 0 : 1
    end
    # display(corrmat)
    GaussianMPS(lx, initconf, positions, θs)
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
# Power flow module. The implementation is a modification of
# MATPOWER's code. We attach the original MATPOWER's license in makeYbus.m:
#
# MATPOWER
# Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC)
# by Ray Zimmerman, PSERC Cornell
#
# Covered by the 3-clause BSD License.
module ExaPF

using CUDA
using CUDA.CUSPARSE
using CUDA.CUSOLVER
using ForwardDiff
using IterativeSolvers
using KernelAbstractions
using Krylov
using LinearAlgebra
using MathOptInterface
using Printf
using SparseArrays
using SparseDiffTools
using TimerOutputs

import Base: show, get

const MOI = MathOptInterface
const TIMER = TimerOutput()

const VERBOSE_LEVEL_HIGH = 3
const VERBOSE_LEVEL_MEDIUM = 2
const VERBOSE_LEVEL_LOW = 1
const VERBOSE_LEVEL_NONE = 0

include("utils.jl")
# Import submodules
include("autodiff.jl")
using .AutoDiff
include("indexes.jl")
using .IndexSet
include("LinearSolvers/LinearSolvers.jl")
using .LinearSolvers
include("parsers/parse_mat.jl")
using .ParseMAT
include("parsers/parse_psse.jl")
using .ParsePSSE
include("PowerSystem/PowerSystem.jl")
using .PowerSystem

const PS = PowerSystem

# Modeling
include("models/models.jl")
# Evaluators
include("Evaluators/Evaluators.jl")

end


text:
#=
###############################################################################
#
# Leaf photosynthesis-related parameter set
# This struct passed the FT test
# This struct is documented in the Leaf page
# More documentation required
#
###############################################################################
"""
    LeafParams{FT<:AbstractFloat}

A structure that contains all parameters for photosynthesis calculations of a single leaf.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef mutable struct LeafParams{FT<:AbstractFloat}
    # TODO Group these variables

    # broadband albedo and emissivity
    "broadband shortwave albedo"
    α::FT = -999
    "longwave emissivity"
    ε::FT = -999

    # thermal-related characteristics
    "Leaf dry mass per area `[kg m⁻²]`, not WET mass"
    LMA::FT = 0.1
    "leaf relative water content"
    RWC::FT = 0.8
    "leaf specific heat capacity `[J kg⁻¹ K⁻¹]`"
    Cleaf::FT = 1000.0

    # TODO remove? pore density in pores/m^2
    #pore_density::FT = 1e8

    # chloroplast volume airspace (m^3 per pore)
    # http://kirschner.med.harvard.edu/files/bionumbers/A%20semi%20quantitative%20analysis%20of%20the%20photosynthetic%20system%20in%20an%20'average'%20C3%20plant%20leaf.pdf
    "Chloroplast volume airspace `[m³ per pore]`"
    Chloroplast_rel_volume::FT = 2.0/100.0

    # Photosynthesis
    "Kn in steady state"
    Kn_ss::FT = 0

    # Fluorescence-related
    "Dynamic Steady state `[Bool]`"
    dynamic_state::Bool = false

    # number of CO₂ per electrons - typically 1/5 for C3 and 1/6 for C4, i.e. about 1/10 for both PSI and PSII!
    "Efficiency of photosynthesis light reactions`[mol CO₂/mol e-]`"
    effcon::FT = 1/5

    "Derivative of saturation vapor pressure at leaf temperature wrt T `[Pa/K]`"
    desat::FT = 0.0

    "Steady state Stomatal conductance `[mol H₂O m⁻² s⁻¹]`"
    gs_ss::FT = 0.1

    # TODO Update this with Plant Hydraulics
    "Sap flow `[?]`"
    Sap::FT = 0;
    "relative humidity at the surface of the leaf `[-]`"
    RH::FT = 100.0

    # VPD gradient ACROSS THE LEAF interface - NOT the weather station one 8-)
    "Vapor pressure difference across leaf interface `[Pa]`"
    VPD::FT = 0.1

    # TODO add elastic modulus? and P-V curve?
    #ε_modulus::FT = 20e6

    # TODO link to Plant Hydraulics
    "tree height `[m]`"
    height = 20
    # -Scholz et al. 2011 Book, Hydraulic Capacitance: Biophysics and Functional Significance,
    # can also be related to P50 see same book
    "tree capacitance `[kg m⁻³ Pa⁻¹]`"
    Ctree::FT = (79 + 8*height)*1e-6;

    "Tree roughness `[m]`"
    z0m::FT= -999.
    # TODO should be changed later
    "tree roughness `[m]`"
    z0h::FT= -999.
    "tree displacement height `[m]`"
    d::FT = -999.

    "leaf thickness (m)"
    dleaf::FT = 2e-3
    "Leaf area index"
    LAI::FT= 1.0
    "turbulent transfer coefficient `[m s``^2``]`"
    Cd::FT= 0.01
end

###############################################################################
#
# Meteo data set
# This struct passed the FT test
# More documentation required
#
###############################################################################
"""
    struct MeteoParams{FT}

A Placeholder for all Meteo data for the moment

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef mutable struct MeteoParams{FT<:Number}
    # TODO Make this struct work first, then clean it up
    "?"
    S_down::FT = 0.0
    "?"
    L_down::FT = 0.0

    "Measurement height `[m]`"
    zscreen::FT= 10.0
    "atmospheric Obukhov length"
    L::FT      = 1e6

    # parameter to define stability function for stable case
    "2 Webb correction tends to be unstable at night - suggest not using"
    stab_type_stable::FT = 1
    "?"
    ustar::FT = 1e-6

    "Air resistance `[Pa]`"
    ra::FT    = 1e6
end
=#

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
"""
	grd2cpt(cmd0::String="", arg1=[], kwargs...)

Make linear or histogram-equalized color palette table from grid

Full option list at [`grd2cpt`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html)

Parameters
----------

- **A** : **alpha** : **transparency** : -- Str --

    Sets a constant level of transparency (0-100) for all color slices.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#a)
- $(GMT.opt_C)
- **D** : -- Str or [] --			Flags = [i|o]

    Select the back- and foreground colors to match the colors for lowest and highest
    z-values in the output CPT.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#d)
- **E** : **nlevels** : -- Int or [] --		Flags = [nlevels]

    Create a linear color table by using the grid z-range as the new limits in the CPT.
    Alternatively, append nlevels and we will resample the color table into nlevels equidistant slices.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#e)
- **F** : **force_rgb** : -- Str or [] --		Flags = [R|r|h|c][+c]]

    Force output CPT to written with r/g/b codes, gray-scale values or color name.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#f)
- **G** : **truncate** : -- Str --              Flags = zlo/zhi

    Truncate the incoming CPT so that the lowest and highest z-levels are to zlo and zhi.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#g)
- **I** : **inverse** : **reverse** : -- Str --	Flags = [c][z]

    Reverse the sense of color progression in the master CPT.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#i)
- **L** : **range** : -- Str --			Flags = minlimit/maxlimit

    Limit range of CPT to minlimit/maxlimit, and don’t count data outside this range when estimating CDF(Z).
    [Default uses min and max of data.]
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#l)
- **M** : **overrule_bg** -- Bool or [] --

    Overrule background, foreground, and NaN colors specified in the master CPT with the values of
    the parameters COLOR_BACKGROUND, COLOR_FOREGROUND, and COLOR_NAN.
    [`-M`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#m)
- **N** : **no_bg** -- Bool or [] --

    Do not write out the background, foreground, and NaN-color fields.
- **Q** : **log** : -- Bool or [] --

    Selects a logarithmic interpolation scheme [Default is linear].
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#q)
- **C** : **row_col** : -- Bool --

    Replace the x- and y-coordinates on output with the corresponding column and row numbers.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#c)
- $(GMT.opt_R)
- **S** : **steps** : -- Bool or [] or Str --			Flags = zstart/zstop/zinc or n

    Set steps in CPT. Calculate entries in CPT from zstart to zstop in steps of (zinc). Default
    chooses arbitrary values by a crazy scheme based on equidistant values for a Gaussian CDF.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#s)
- **T** : **symetric** : -- Str --			Flags = -|+|_|=

    Force the color table to be symmetric about zero (from -R to +R).
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#t)
- $(GMT.opt_V)
- **W** : **wrap** : **categorical** : -- Bool or Str or [] --      Flags = [w]

    Do not interpolate the input color table but pick the output colors starting at the
    beginning of the color table, until colors for all intervals are assigned.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#w)
- **Z** : **continuous** : -- Bool or [] --

    Creates a continuous CPT [Default is discontinuous, i.e., constant colors for each interval].
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/grd2cpt.html#z)
- $(GMT.opt_V)
- $(GMT.opt_write)
"""
function grd2cpt(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grd2cpt", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd, = parse_R("", d)
	cmd = parse_V_params(cmd, d)

	cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', 0, arg1)

	cmd = add_opt(cmd, 'A', d, [:A :alpha :transparency])
	cmd = add_opt(cmd, 'D', d, [:D :bg :background])
	cmd = add_opt(cmd, 'E', d, [:E :nlevels])
	cmd = add_opt(cmd, 'F', d, [:F :force_rgb])
	cmd = add_opt(cmd, 'G', d, [:G :truncate])
	cmd = add_opt(cmd, 'I', d, [:I :inverse :reverse])
	cmd = add_opt(cmd, 'L', d, [:L :limit])
	cmd = add_opt(cmd, 'M', d, [:M :overrule_bg])
	cmd = add_opt(cmd, 'N', d, [:N :no_bg])
	cmd = add_opt(cmd, 'Q', d, [:Q :log])
	cmd = add_opt(cmd, 'S', d, [:S :steps])
	cmd = add_opt(cmd, 'T', d, [:T :symetric])
	cmd = add_opt(cmd, 'W', d, [:W :no_interp])
	cmd = add_opt(cmd, 'Z', d, [:Z :continuous])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	if (isa(arg1, Array{<:Number}))		arg1 = mat2grid(arg1)	end
	return common_grd(d, cmd, got_fname, 1, "grd2cpt", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grd2cpt(arg1=[], cmd0::String=""; kw...) = grd2cpt(cmd0, arg1; kw...)


text:
import WaveletScattering: kthrange, PathKey, Literal

r = kthrange([:time], rand(2), 1)
@test r == (PathKey(Literal[Literal(:time,1)])=>1:1:2)


text:
function DBInterface.transaction(f, con::Connection)
    begin_transaction(con)
    try
        f()
    catch
        rollback(con)
        rethrow()
    finally
        commit(con)
    end
end

function DBInterface.transaction(f, db::DB)
    return DBInterface.transaction(f, db.main_connection)
end

"""
    DuckDB.begin(db)

begin a transaction
"""
function begin_transaction end

begin_transaction(con::Connection) = execute(con, "BEGIN TRANSACTION;")
begin_transaction(db::DB) = begin_transaction(db.main_connection)
transaction(con::Connection) = begin_transaction(con)
transaction(db::DB) = begin_transaction(db)

"""
    DuckDB.commit(db)

commit a transaction
"""
function commit end

commit(con::Connection) = execute(con, "COMMIT TRANSACTION;")
commit(db::DB) = commit(db.main_connection)

"""
    DuckDB.rollback(db)

rollback transaction
"""
function rollback end

rollback(con::Connection) = execute(con, "ROLLBACK TRANSACTION;")
rollback(db::DB) = rollback(db.main_connection)

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
module TestComprehensions

include("preamble.jl")

@testset "eduction(comprehension)" begin
    iter = (x for x in 1:10 if x % 2 == 0)
    ed = eduction(iter)
    @test collect(iter) == collect(ed)

    iter = (y for x in 1:10 if x % 2 == 0 for y in (x, x + 1))
    ed = eduction(iter)
    @test collect(iter) == collect(ed)

    iter = ((x, y) for x in 1:3, y in 1:2)
    ed = eduction(iter)
    @test vec(collect(iter)) == collect(ed)
end

@testset "Transducer(comprehension)" begin
    xf1 = Transducer(2x for x in nothing if x % 2 == 0)
    xf2 = Filter(x -> x % 2 == 0) |> Map(x -> 2x)
    xs = 1:10
    @test collect(xf1, xs) == collect(xf2, xs)
end

end  # module


text:
using Test, Bridge

# test the tricky bridge sampling when endpoint is or is not in grid
#srand(10)
G =  GammaProcess(10.0, 1.5)
GB = GammaBridge(1.0, 2.0, G)
n = 1000
@test abs(mean(sample([0.0, 1.0, 3.0],  G).yy[end] for i in 1:1000) -  mean(Bridge.increment(3, G))) < 4*std(Bridge.increment(0.5, G))/sqrt(n)

X1 = sample([0.0, 0.5, 1.0], GB, 0.2)
X2 = sample([0.0, 0.5, 2.0], GB, 0.2)
X3 = sample([0.0, 0.5], GB, 0.2)

@test  X1.tt ≈ [0,1,2]/2
@test  X2.tt ≈ [0,1,4]/2
@test  X3.tt ≈ [0,1]/2

@test  X1.yy[1] == X1.yy[1] == X1.yy[1] == 0.2

@test  X1.yy[3] ≈ 2.0
@test  X2.yy[3] >= 2.0
@test  X3.yy[2] <= 2.0

γᵒ, γ = 0.75, 1.2
t = 10000.
α = 16.0
n = 10000
tt = range(0, stop=t, length=n)
P = GammaProcess(γ, α)
X = sample(tt, P)

Y = copy(X)
Bridge.uniform_thinning!(Y, P, γᵒ)

@test abs(γ*t/n/α - mean(diff(X.yy))) < 0.2/sqrt(n)
@test abs(γᵒ*t/n/α - mean(diff(Y.yy))) < 0.2/sqrt(n)

γ = 10.
α = 5.
P = GammaProcess(γ, α)
T = 5.0
tt = 0.0:T/n:T
m = 1000
sigma = (t, x = 0.0) -> exp(-cos(t))
L = SamplePath(tt, sample(tt, P).yy .- sample(tt, P).yy)
X = solve(EulerMaruyama(), 0.0, L, ((t,x)-> -0*x, sigma))


text:
module OptionalArgChecks

using IRTools: @dynamo, IR, recurse!, block, branches, branch!
using MacroTools: postwalk

export @mark, @skip, @unsafe_skipargcheck

# reexport @argcheck and @check
using ArgCheck: @argcheck, @check, LABEL_ARGCHECK
export @argcheck, @check

"""
    @mark label ex

Marks `ex` as an optional argument check, so when a function is called via
[`@skip`](@ref) with label `label`, `ex` will be omitted.

```jldoctest
julia> function half(x::Integer)
           @mark check_even iseven(x) || throw(DomainError(x, "x has to be an even number"))
           return x ÷ 2
       end
half (generic function with 1 method)

julia> half(4)
2

julia> half(3)
ERROR: DomainError with 3:
x has to be an even number
[...]

julia> @skip check_even half(3)
1
```
"""
macro mark(label, ex)
    label isa Symbol || error("label has to be a Symbol")
    return Expr(
        :block,
        Expr(:meta, :begin_optional, label),
        esc(ex),
        Expr(:meta, :end_optional, label),
    )
end

struct Skip{labels,recursive}
    Skip(labels...; recursive=true) = new{labels,recursive}()
end

@dynamo function (::Skip{labels,recursive})(x...) where {labels,recursive}
    ir = IR(x...)
    ir === nothing && return

    tape = ()
    local orig

    for (x, st) in ir
        is_layer_begin = Meta.isexpr(st.expr, :meta) &&
            st.expr.args[1] === :begin_optional &&
            st.expr.args[2] in labels
        is_layer_end =  Meta.isexpr(st.expr, :meta) &&
            st.expr.args[1] === :end_optional &&
            !isempty(tape) && st.expr.args[2] === last(tape)

        if is_layer_begin
            tape = (tape..., st.expr.args[2])
        elseif is_layer_end
            tape = Base.front(tape)
        end

        is_begin = is_layer_begin && length(tape) == 1
        is_end   = is_layer_end   && isempty(tape)

        if is_begin
            orig = block(ir, x)
        elseif is_end
            dest = block(ir, x)
            if orig != dest
                empty!(branches(orig))
                branch!(orig, dest)
            end
        end

        if is_layer_begin || is_layer_end || !isempty(tape)
            delete!(ir, x)
        end
    end

    @assert isempty(tape)

    recursive && recurse!(ir)
    return ir
end

function _skip(l, ex, recursive=true)
    labels::Vector{Symbol} = if l isa Symbol
        [l]
    elseif Meta.isexpr(l, :vect)
        labels = l.args
    else
        error("label has to be a name or array of names")
    end
    return _skip(labels, ex, recursive)
end

function _skip(labels::Vector, ex, recursive)
    ex = postwalk(ex) do x
        if Meta.isexpr(x, :call)
            pushfirst!(x.args, Expr(
                :call,
                GlobalRef(@__MODULE__, :Skip),
                Expr(:parameters, Expr(:kw, :recursive, recursive)),
                Expr.(:quote, labels)...
            ))
        end
        return x
    end
    return esc(ex)
end

"""
    @skip label ex[[ recursive=true]]
    @skip [label1, label2, ...] ex[[ recursive=true]]

For every function call in `ex`, expressions marked with label `label` (or any of the labels
`label*` respectively) using the macro [`@mark`](@ref) get omitted recursively.
"""
macro skip end

macro skip(l, ex)
    return _skip(l, ex, true)
end

macro skip(l, ex, r)
    r = parse_recursive(r)
    return _skip(l, ex, r)
end

function parse_recursive(r)::Bool
    Meta.isexpr(r, :(=)) || error("expected keyword argument instead of `$r`")
    argname = r.args[1]
    argname == :recursive || error("unknown kewyword argument `$argname`")
    recursive = r.args[2]
    recursive isa Bool || error("keyword argument `recursive` has to be a `Bool` literal")
    return recursive
end

"""
    @unsafe_skipargcheck ex[[ recursive=true]]

Elides argument checks created with [`@argcheck`](@ref) or [`@check`](@ref),
provided by the package `ArgCheck.jl`.
"""
macro unsafe_skipargcheck(ex)
    return _skip([LABEL_ARGCHECK], ex, true)
end

macro unsafe_skipargcheck(ex, r)
    r = parse_recursive(r)
    return _skip([LABEL_ARGCHECK], ex, r)
end

end


text:
module Config

    export add_config, get_config, set_config,
           get_config_dictionary

    if !(@isdefined GICONFIG)
        global GICONFIG = Dict()
    end

    function add_config(name, value)
        if !haskey(GICONFIG, name)
            GICONFIG[name] = value
        end
    end

    function set_config(name, value)
        if haskey(GICONFIG, name)
            GICONFIG[name] = value
        else
            println("  WARNING: Unknown parameter name.")
        end
    end

    function get_config(name)
        if haskey(GICONFIG, name)
            return GICONFIG[name]
        else
            println("  WARNING: Unknown parameter name.")
            return nothing
        end
    end

    function get_config_dictionary()
        return Dict(GICONFIG)
    end

end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMIImport.FMICore: fmi2VariableNamingConventionStructured

tool = ENV["EXPORTINGTOOL"]
pathToFMU = "https://github.com/ThummeTo/FMI.jl/raw/main/model/" * ENV["EXPORTINGTOOL"] * "/SpringFrictionPendulum1D.fmu"

myFMU = fmi2Load(pathToFMU)

@test fmi2GetVersion(myFMU) == "2.0"
@test fmi2GetTypesPlatform(myFMU) == "default"

@test fmi2GetModelName(myFMU) == "SpringFrictionPendulum1D"
@test fmi2GetVariableNamingConvention(myFMU) == fmi2VariableNamingConventionStructured
@test fmi2IsCoSimulation(myFMU) == true
@test fmi2IsModelExchange(myFMU) == true

if tool == "Dymola/2020x"
    @test fmi2GetGUID(myFMU) == "{b02421b8-652a-4d48-9ffc-c2b223aa1b94}"
    @test fmi2GetGenerationTool(myFMU) == "Dymola Version 2020x (64-bit), 2019-10-10"
    @test fmi2GetGenerationDateAndTime(myFMU) == "2021-11-23T13:36:30Z"
    @test fmi2GetNumberOfEventIndicators(myFMU) == 24
    @test fmi2CanGetSetState(myFMU) == true
    @test fmi2CanSerializeFMUstate(myFMU) == true
    @test fmi2ProvidesDirectionalDerivative(myFMU) == true

    @test fmi2GetDefaultStartTime(myFMU.modelDescription) ≈ 0.0
    @test fmi2GetDefaultStopTime(myFMU.modelDescription) ≈ 1.0
    @test fmi2GetDefaultTolerance(myFMU.modelDescription) ≈ 1e-4
    @test fmi2GetDefaultStepSize(myFMU.modelDescription) === nothing

elseif tool == "OpenModelica/v1.17.0"
    @test fmi2GetGUID(myFMU) == "{8584aa5b-179e-44ed-9ba6-d557ed34541e}"
    @test fmi2GetGenerationTool(myFMU) == "OpenModelica Compiler OMCompiler v1.17.0"
    @test fmi2GetGenerationDateAndTime(myFMU) == "2021-06-21T11:48:49Z"
    @test fmi2GetNumberOfEventIndicators(myFMU) == 14
    @test fmi2CanGetSetState(myFMU) == false
    @test fmi2CanSerializeFMUstate(myFMU) == false
    @test fmi2ProvidesDirectionalDerivative(myFMU) == false

    @test fmi2GetDefaultStartTime(myFMU.modelDescription) ≈ 0.0
    @test fmi2GetDefaultStopTime(myFMU.modelDescription) ≈ 1.0
    @test fmi2GetDefaultTolerance(myFMU.modelDescription) ≈ 1e-6
    @test fmi2GetDefaultStepSize(myFMU.modelDescription) === nothing
else
    @warn "Unknown exporting tool `$tool`. Skipping model description tests."
end

fmi2Unload(myFMU)


text:
# TODO: what about @objective(m, Max, count(x .== 1))? Automatically add a constraint (i.e. bridge)? And/or support the constraint as a function?
# TODO: -> rather all constraints, more consistent with how other solvers work.

function _update_objective(model::Optimizer)
    # If the sense is feasibility and there is an internal Concert objective, remove it.
    # Otherwise, this is an optimisation problem.
    if model.objective_sense == MOI.FEASIBILITY_SENSE &&
       model.objective_cp !== nothing
        cpo_java_remove(model.inner, model.objective_cp)
        model.objective_cp = nothing
    end

    # If only no function is available, don't do anything.
    if model.objective_function_cp === nothing
        return
    end

    # Set the new objective.
    obj = if model.objective_sense == MOI.MIN_SENSE
        cpo_java_minimize(model.inner, model.objective_function_cp)
    else
        cpo_java_maximize(model.inner, model.objective_function_cp)
    end
    cpo_java_add(model.inner, obj)
    return
end

function MOI.set(
    model::Optimizer,
    ::MOI.ObjectiveSense,
    sense::MOI.OptimizationSense,
)
    model.objective_sense = sense
    _update_objective(model)
    return
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveSense)
    return model.objective_sense
end

function MOI.set(
    model::Optimizer,
    ::MOI.ObjectiveFunction{F},
    f::F,
) where {F <: MOI.VariableIndex}
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        convert(MOI.ScalarAffineFunction{Float64}, f),
    )
    model.objective_type = MOI.SINGLE_VARIABLE
    return
end

function MOI.get(
    model::Optimizer,
    ::MOI.ObjectiveFunction{F},
) where {F <: MOI.AbstractScalarFunction}
    if typeof(model.objective_function) <: F
        return model.objective_function
    else
        error(
            "Unable to get objective function. Current objective: $(model.objective_function).",
        )
    end
end

function MOI.set(
    model::Optimizer,
    ::MOI.ObjectiveFunction{F},
    f::F,
) where {F <: MOI.AbstractScalarFunction}
    model.objective_function = f
    model.objective_function_cp = _parse(model, f)
    _update_objective(model)
    return
end

# TODO: modifications. Easy to do, as we have a pointer on the Concert expression! Hard to do, as the MOI function must be rebuilt.
# function MOI.modify(
#     model::Optimizer,
#     ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
#     chg::MOI.ScalarConstantChange{Float64}
# )
#     CPLEX.c_api_chgobjoffset(model.inner, chg.new_constant)
#     return
# end


text:
using AbstractTrees
using AbstractTrees: TreeCharSet

# Node type that wraps an integer value `n` and has two children with value `n + 1`.
# Resulting tree has infinite height.
abstract type AbstractNumTree; end
struct Num <: AbstractNumTree
    n::Int
end
Base.show(io::IO, x::AbstractNumTree) = print(io, x.n)
AbstractTrees.children(x::Num) = (Num(x.n+1), Num(x.n+1))

struct SingleChildInfiniteDepth end
AbstractTrees.children(::SingleChildInfiniteDepth) = (SingleChildInfiniteDepth(),)

# Wrapper around a collection which does not define getindex() or keys()
struct Unindexable{C}
    col::C

    Unindexable(col) = new{typeof(col)}(col)
end

Base.eltype(::Type{Unindexable{C}}) where C = eltype(C)
Base.length(u::Unindexable) = length(u.col)
Base.iterate(u::Unindexable, args...) = iterate(u.col, args...)

# Wrapper around a node which wraps children in UnIndexable
struct UnindexableChildren{N}
    node::N
end

AbstractTrees.children(u::UnindexableChildren) = Unindexable(children(u.node))
AbstractTrees.printnode(io::IO, u::UnindexableChildren) = AbstractTrees.printnode(io, u.node)

@testset "Truncation" begin

    # test that `print_tree(headnode, maxdepth)` truncates the output at right depth
    # julia > print_tree(Num(0), 3)
    # 0
    # ├─ 1
    # │  ├─ 2
    # │  │  ├─ 3
    # │  │  └─ 3
    # │  └─ 2
    # │     ├─ 3
    # │     └─ 3
    # └─ 1
    #    ├─ 2
    #    │  ├─ 3
    #    │  └─ 3
    #    └─ 2
    #       ├─ 3
    #       └─ 3
    #

    truncation_str = TreeCharSet().trunc

    for maxdepth in [3,5,8]
        ptxt = repr_tree(Num(0), maxdepth=maxdepth)

        # Check that we see depth #s the expected number of times
        @test length(findall(string(maxdepth-1), ptxt)) == 2^(maxdepth-1)
        @test length(findall(string(maxdepth), ptxt)) == 2^maxdepth
        @test length(findall(string(maxdepth+1), ptxt)) == 0

        # test that `print_tree(headnode)` prints truncation characters under each
        # node at the maximum depth
        lines = split(ptxt, '\n')
        for i in 1:length(lines)
            if occursin(string(maxdepth), lines[i])
                @test endswith(lines[i+1], truncation_str)
            end
        end
    end

    # test correct number of lines printed 1
    ptxt = repr_tree(SingleChildInfiniteDepth())
    numlines = sum([1 for c in split(ptxt, '\n') if ~isempty(strip(c))])
    @test numlines == 7 # 1 (head node) + 5 (default depth) + 1 (truncation char)

    # test correct number of lines printed 2
    ptxt = repr_tree(SingleChildInfiniteDepth(), maxdepth=3)
    numlines = sum([1 for c in split(ptxt, '\n') if ~isempty(strip(c))])
    @test numlines == 5 # 1 (head node) + 3 (depth) + 1 (truncation char)

    # test correct number of lines printed 3
    ptxt = repr_tree(SingleChildInfiniteDepth(), maxdepth=3, indicate_truncation=false)
    numlines = sum([1 for c in split(ptxt, '\n') if ~isempty(strip(c))])
    @test numlines == 4 # 1 (head node) + 3 (depth)
end

@testset "Child keys" begin
    @testset "AbstractVector" begin
        tree = 1:2

        # let's not bother testing the exact string of the type itself, it breaks too easily

        @test endswith(repr_tree(tree), """
                       ├─ 1
                       └─ 2
                       """)

        @test endswith(repr_tree(tree, printkeys=true), """
                       ├─ 1 ⇒ 1
                       └─ 2 ⇒ 2
                       """)
    end

    @testset "Tuple" begin
        tree = (1, 2)

        @test endswith(repr_tree(tree), """
                       ├─ 1
                       └─ 2
                       """)

        @test endswith(repr_tree(tree, printkeys=true), """
                       ├─ 1 ⇒ 1
                       └─ 2 ⇒ 2
                       """)
    end

    @testset "Matrix" begin
        tree = [1 2; 3 4]

        @test endswith(repr_tree(tree), """
        ├─ (1, 1) ⇒ 1
        ├─ (2, 1) ⇒ 3
        ├─ (1, 2) ⇒ 2
        └─ (2, 2) ⇒ 4
        """)

        @test endswith(repr_tree(tree, printkeys=false), """
        ├─ 1
        ├─ 3
        ├─ 2
        └─ 4
        """)
    end

    @testset "No keys" begin
        tree = UnindexableChildren(1:2)  # Has no method for Base.keys()

        @test repr_tree(tree) == repr_tree(tree.node)
        @test repr_tree(tree, printkeys=true) == repr_tree(tree.node)
    end
end

@testset "print_tree with non-indexable children" begin
    tree = Num(0)
    @test repr_tree(UnindexableChildren(tree), maxdepth=4) == repr_tree(tree, maxdepth=4)
end

# Prints node as cool message box
struct BoxNode
    s::String
    children::Vector
end

AbstractTrees.children(n::BoxNode) = n.children
function AbstractTrees.printnode(io::IO, n::BoxNode)
    println(io, "┌", "─" ^ textwidth(n.s), "┐")
    println(io, "│", n.s, "│")
    print(io, "└", "─" ^ textwidth(n.s), "┘")
end

@testset "printnode multiline" begin
    tree = ["foo", BoxNode("bar", [1, 2:4, 5]), "baz"]

    @test endswith(repr_tree(tree), """
                   ├─ "foo"
                   ├─ ┌───┐
                   │  │bar│
                   │  └───┘
                   │  ├─ 1
                   │  ├─ 2:4
                   │  │  ├─ 2
                   │  │  ├─ 3
                   │  │  └─ 4
                   │  └─ 5
                   └─ "baz"
                   """)
end


text:
# This is the module file. It allows communication between different files within the src folder
# Also, determines which functions are provided by AutomotiveInteraction whenever someone says
# using AutomotiveInteraction
module AutomotiveInteraction

using AutomotiveSimulator
using AutomotiveVisualization
using Reel

using DataFrames
using DelimitedFiles
using Parameters
using Random # For AbstractRNG in Base.rand for c-idm

using Distributions # For pdf(Normal) in weight_and_resample
using StatsBase # For weight_and_resample
using LinearAlgebra # Form norm calculation for iterwise particle dist

using Combinatorics # make param pairs for filtering progress pairwise plots
using PGFPlots

export
    MergingEnvironment,
    MergingEnvironmentLower,
    FilteringEnvironment
include("Driving/environment.jl")

export
    keep_vehicle_subset!,
    make_def_models,
    make_IDM_models,
    make_cidm_models,
    make_lmidm_models,
    make_TimLaneChanger_models,
    get_hallucination_scenes,
    run_vehicles,
    compare2truth,
    run_vehicles_curvept_overlay,
    test_barrier_vehicle,
    test_jumpy_vehicle
include("Driving/driving_simulation.jl")

export
    centerlines_txt2tracks,
    make_roadway_interaction,
    make_roadway_ngsim,
    make_roadway_interaction_with_extensions,
    MergingRoadway,
    add_renderable!
include("Data_Processing/roadway_building.jl")

export
    INTERACTIONTrajdata,
    carsinframe,
    car_df_index,
    read_veh_tracks,
    extract_timespan,
    sample_simultaneous_vehs
include("Data_Processing/veh_track_reading.jl")

export
    video_trajdata_replay,
    scenelist2video,
    video_overlay_scenelists,
    scenelist2video_curvepts
include("Driving/visualization.jl")

export
    CooperativeIDM,
    find_merge_vehicle
include("Driving/cooperative_IDM.jl")

export
    MergeOverlay
include("Driving/overlays.jl")

export
    cidm_from_particle,
    hallucinate_a_step,
    weight_and_resample,
    multistep_update,
    obtain_driver_models
include("Filtering/particle_filtering.jl")

export
    initial_pmat,
    plot_pairwise_particles,
    gen_imitation_traj,
    imitation_with_replay,
    compute_rmse,
    rmse_dict2mean,
    test_collision,
    frac_colliding_timesteps,
    pgfplot_vector,
    truncate_vecs,
    pad_zeros,
    reality_metrics,
    replay_scenelist
include("Filtering/utils.jl")

end # End module

max_stars_repo_name: sisl/AutomotiveInteraction

max_issues_repo_name: sisl/AutomotiveInteraction

max_forks_repo_name: sisl/AutomotiveInteraction


text:
using HDF5

runtest(filename) = (println(filename); include(filename))

runtest("plain.jl")
runtest("readremote.jl")
runtest("extend_test.jl")
runtest("gc.jl")

nothing


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you 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.

# Autograd for NDArray
# this is a port of Python's autograd module
# https://github.com/apache/incubator-mxnet/blob/master/python/mxnet/autograd.py

###############################################################################
#  Private util functions
###############################################################################

"""
    _set_recording(state::Bool)::Bool

Set status to recording/not recording. When recording, graph will be constructed
for gradient computation.

## Parameters

* `state::Bool`

## Returns

Previous state before this set
"""
function _set_recording(state::Bool)::Bool
  prev = Ref{Cint}(C_NULL)
  @mxcall(:MXAutogradSetIsRecording, (Cint, Ref{Cint}), state, prev)
  prev[]
end

_set_recording(::Void) = nothing

"""
Set status to training/predicting.
For example, Dropout will drop inputs randomly when
`train_mode = true` while simply passing through if `train_mode = false`.

## Parameters
* `train_mode::Bool`

## Returns

Previous state before this set.
"""
function _set_training(train_mode::Bool)::Bool
  prev = Ref{Cint}(C_NULL)
  @mxcall(:MXAutogradSetIsTraining, (Cint, Ref{Cint}), train_mode, prev)
  prev[]
end

_set_training(::Void) = nothing

###############################################################################
#  Public API
###############################################################################

"""
    is_recording()::Bool

Get status on recording/not recording.
"""
function is_recording()::Bool
  state = Ref{Cint}(C_NULL)
  @mxcall(:MXAutogradIsRecording, (Ref{Cint},), state)
  state[]
end

"""
    is_training()::Bool

Get status on recording/not recording.
"""
function is_training()::Bool
  state = Ref{Cint}(C_NULL)
  @mxcall(:MXAutogradIsTraining, (Ref{Cint},), state)
  state[]
end

@inline function _record(f, is_record::Union{Void,Bool}, train_mode::Union{Void,Bool})
  # Port from Python's `_RecordingStateScope` context manager
  # __enter__
  prev_is_record = _set_recording(is_record)
  prev_train_mode = _set_training(train_mode)

  try
    f()
  finally
    # __exit__
    if is_record != nothing && prev_is_record != is_record
      _set_recording(prev_is_record)
    end
    if train_mode != nothing && prev_train_mode != train_mode
      _set_recording(prev_train_mode)
    end
  end
end

"""
    record(f, train_mode = true)
    record(translates = true) do
      ...
    end

Returns an autograd recording scope context to be used in `do` block
and captures code that needs gradients to be calculated.

Parameter `train_mode::Bool` controls whether the forward pass is in training
or predicting mode.
This controls the behavior of some layers such as `Dropout`, `BatchNorm`.

!!! note
    When forwarding with `train_mode = false`, the corresponding backward
    should also use `train_mode = false`, otherwise gradient is undefined.

```julia
x = mx.NDArray([1 2; 3 4])
∇ = mx.attach_grad!(x)
y = mx.record() do
  2x
end
mx.backward!(y)

julia> ∇
2×2 mx.NDArray{Int64,2} @ CPU0:
 2  2
 2  2
```
"""
record(f, train_mode::Bool = true) = _record(f, true, train_mode)

"""
    pause(f, train_mode = false)
    pause(train_mode = false) do
      ...
    end

Create a scope context for codes that do not need gradients to be calculated.

```julia
record() do
  ...
  pause() do
    # testing, IO, gradient updates...
  end
end
```
"""
pause(f, train_mode::Bool = false) = _record(f, false, train_mode)

"""
    train_mode(f)
    train_mode() do
      ...
    end

Create a scope context in which forward pass behavior is set to training mode,
without changing the recording states.

```julia
y = model(x)
train_mode() do
  z = mx.Dropout(y)
  ...
end
```
"""
train_mode(f) = _record(f, nothing, true)

"""
    predict_mode(f)
    predict_mode() do
      ...
    end

Create a scope context in which forward pass behavior is set to inference mode,
without changing the recording states.

```julia
record() do
  y = model(x)
  predict_mode() do
    y = sampling(y)
  end
end
```
"""
predict_mode(f) = _record(f, nothing, false)

"""
    backward!(head,  head_grad;  retain_graph = false, train_mode = true)
    backward!(heads, head_grads; retain_graph = false, train_mode = true)

Compute the gradients of heads w.r.t previously marked variables.

## Parameters

- `head::NDArray`: output NDArray

- `head_grad::NDArray` or `Void`: gradient coefficient with respect to head.

- `heads::Vector{NDArray}`: a list of output NDArray

- `head_grads::Vector`: a list of gradient coefficient with respect ot heads.
  the element should be `NDArray` or `Void`

- `retain_graph::Bool`: whether to keep the graph after backward. e.g:
  If you want to differentiate the same graph twice,
  you need to pass `retain_graph=true`.

- `train_mode::Bool`: whether to do backward for training or predicting.
"""
backward!(head::NDArray, head_grad::NDArray; kws...) =
  backward!([head], [head_grad]; kws...)

backward!(head::NDArray, head_grad::Void = nothing; kws...) =
  backward!([head], head_grad; kws...)

function backward!(heads::VecOfNDArray, head_grad::Void;
                   retain_graph::Bool = false, train_mode::Bool = true)
  @mxcall(
    :MXAutogradBackwardEx,
    (MX_uint,
     Ptr{MX_handle},
     Ptr{MX_handle},
     MX_uint,
     Ptr{MX_handle},
     Cint,
     Cint,
     Cint,
     Ptr{MX_handle},
     Ptr{MX_handle}),
    length(heads),
    map(x -> x.handle, heads),
    C_NULL,
    0,
    C_NULL,
    retain_graph,
    false,  # create_graph
    train_mode,
    C_NULL,
    C_NULL)
end

function backward!(heads::VecOfNDArray, head_grads::Vector;
                   retain_graph::Bool = false, train_mode::Bool = true)
  output_handles = map(x -> x.handle, heads)
  ograd_handles  = map(head_grads) do x
    if x isa NDArray
      x.handle
    elseif x isa Void
      MX_handle(C_NULL)
    else
      throw(ArgumentError("element of head_grads should be NDArray or Void"))
    end
  end
  @assert length(output_handles) == length(ograd_handles)
  @mxcall(
    :MXAutogradBackwardEx,
    (MX_uint,
     Ptr{MX_handle},
     Ptr{MX_handle},
     MX_uint,
     Ptr{MX_handle},
     Cint,
     Cint,
     Cint,
     Ptr{MX_handle},
     Ptr{MX_handle}),
    length(output_handles),
    output_handles,
    ograd_handles,
    0,
    C_NULL,
    retain_graph,
    false,  # create_graph
    train_mode,
    C_NULL,
    C_NULL)
end

"""
    getgrad(arr::NDArray)

Returns the gradient buffer attached to this `NDArray`.
If the gradient buffer isn't attached yet, return `nothing`.
"""
function getgrad(arr::NDArray)
  out = Ref{MX_handle}(C_NULL)
  @mxcall(:MXNDArrayGetGrad, (MX_handle, Ref{MX_handle}), arr.handle, out)
  (out[] == C_NULL) ? nothing : NDArray(MX_NDArrayHandle(out[]))
end

"""
    attach_grad!(x::NDArray, grad_req::Symbol = :write)

Attach a gradient buffer to this `NDArray`,
so that [`backward!`](@ref) can compute gradient with respect to it.

## Parameters

- `x::NDArray`
- `grad_req::Symbol` (default is `:write`)

## Return

The attached gradient buffer

## See also

- [`getgrad`](@ref)
"""
function attach_grad!(x::NDArray, grad_req::Symbol = :write)
  # TODO: support storage type (stype in Python)
  # TODO: make sure it works with gpu array
  grad = zeros_like(x)
  _mark_variables!([x], [grad], grad_req)
  grad
end

"""
    mark_variables!(var,  grad,  grad_req)
    mark_variables!(vars, grads, grad_reqs)

Mark `NDArrays` as variables to compute gradient for autograd.

## Parameters

- `var::NDArray`
- `grad::NDArray`
- `grad_req::Symbol`: `:nop`, `:write`, `:inplace` or `:add`
- `vars::Vector{NDArray}`
- `grads::Vector{NDArray}`
- `grad_req::Vector{Symbol}`
"""
mark_variables!(var::NDArray, grad::NDArray, grad_reqs::Symbol = :write) =
  _mark_variables!([var], [grad], grad_reqs)

mark_variables!(var::VecOfNDArray, grads::VecOfNDArray, grad_reqs = :write) =
  _mark_variables!(var, grads, grad_reqs)

@inline function _getgrad_req(x::Symbol)::GRAD_REQ
  val = get(grad_req_map, x, false)
  if val == false
    throw(ArgumentError("invalid grad_reqs $x"))
  end
  val
end

@inline _getgrad_reqs(x::Symbol, n::Int) =
  map((_) -> MX_uint(_getgrad_req(x)), Base.OneTo(n))

@inline function _getgrad_reqs(xs::Vector{Symbol}, n::Int)
  if length(xs) != n
    throw(ArgumentError("number of variables and grad_reqs not matched"))
  end
  map(MX_uint ∘ _getgrad_req, xs)
end

@inline function _mark_variables!(vars::VecOfNDArray, grads::VecOfNDArray,
                                  grad_reqs = :write)
  n = length(vars)
  if n != length(grads)
    throw(ArgumentError("number of variables and gradients not matched"))
  end

  var_hdls  = map(x -> x.handle, vars)
  grad_hdls = map(x -> x.handle, grads)
  grad_reqs = _getgrad_reqs(grad_reqs, n)

  @mxcall(:MXAutogradMarkVariables,
          (MX_uint, Ref{MX_handle}, Ptr{MX_uint}, Ref{MX_handle}),
          length(vars), var_hdls, grad_reqs, grad_hdls)
end

"""
    symbol(x::NDArray)

Retrieve recorded computation history as `SymbolicNode`,
 where `x` is a `NDArray` representing the head of computation graph.
 """
function symbol(x::NDArray)
  ref = Ref{MX_handle}(C_NULL)
  @mxcall(:MXAutogradGetSymbol, (MX_handle, Ref{MX_handle}), x, ref)
  SymbolicNode(MX_SymbolHandle(ref[]))
end

###############################################################################
#  TODO: User-defined differentiable function
###############################################################################

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet