text:
module API

using CEnum

using BlueZ_jll

struct bdaddr_t
    b::NTuple{6, UInt8}
end

function ntoh64(n)
    ccall((:ntoh64, libbluetooth), UInt64, (UInt64,), n)
end

struct uint128_t
    data::NTuple{16, UInt8}
end

function ntoh128(src, dst)
    ccall((:ntoh128, libbluetooth), Cvoid, (Ptr{uint128_t}, Ptr{uint128_t}), src, dst)
end

function btoh128(src, dst)
    ccall((:btoh128, libbluetooth), Cvoid, (Ptr{uint128_t}, Ptr{uint128_t}), src, dst)
end

struct bt_security
    level::UInt8
    key_size::UInt8
end

struct bt_power
    force_active::UInt8
end

struct bt_voice
    setting::UInt16
end

@cenum __JL_Ctag_21::UInt32 begin
    BT_CONNECTED = 1
    BT_OPEN = 2
    BT_BOUND = 3
    BT_LISTEN = 4
    BT_CONNECT = 5
    BT_CONNECT2 = 6
    BT_CONFIG = 7
    BT_DISCONN = 8
    BT_CLOSED = 9
end

function bt_get_le64(ptr)
    ccall((:bt_get_le64, libbluetooth), UInt64, (Ptr{Cvoid},), ptr)
end

function bt_get_be64(ptr)
    ccall((:bt_get_be64, libbluetooth), UInt64, (Ptr{Cvoid},), ptr)
end

function bt_get_le32(ptr)
    ccall((:bt_get_le32, libbluetooth), UInt32, (Ptr{Cvoid},), ptr)
end

function bt_get_be32(ptr)
    ccall((:bt_get_be32, libbluetooth), UInt32, (Ptr{Cvoid},), ptr)
end

function bt_get_le16(ptr)
    ccall((:bt_get_le16, libbluetooth), UInt16, (Ptr{Cvoid},), ptr)
end

function bt_get_be16(ptr)
    ccall((:bt_get_be16, libbluetooth), UInt16, (Ptr{Cvoid},), ptr)
end

function bt_put_le64(val, ptr)
    ccall((:bt_put_le64, libbluetooth), Cvoid, (UInt64, Ptr{Cvoid}), val, ptr)
end

function bt_put_be64(val, ptr)
    ccall((:bt_put_be64, libbluetooth), Cvoid, (UInt64, Ptr{Cvoid}), val, ptr)
end

function bt_put_le32(val, ptr)
    ccall((:bt_put_le32, libbluetooth), Cvoid, (UInt32, Ptr{Cvoid}), val, ptr)
end

function bt_put_be32(val, ptr)
    ccall((:bt_put_be32, libbluetooth), Cvoid, (UInt32, Ptr{Cvoid}), val, ptr)
end

function bt_put_le16(val, ptr)
    ccall((:bt_put_le16, libbluetooth), Cvoid, (UInt16, Ptr{Cvoid}), val, ptr)
end

function bt_put_be16(val, ptr)
    ccall((:bt_put_be16, libbluetooth), Cvoid, (UInt16, Ptr{Cvoid}), val, ptr)
end

function bacmp(ba1, ba2)
    ccall((:bacmp, libbluetooth), Cint, (Ptr{bdaddr_t}, Ptr{bdaddr_t}), ba1, ba2)
end

function bacpy(dst, src)
    ccall((:bacpy, libbluetooth), Cvoid, (Ptr{bdaddr_t}, Ptr{bdaddr_t}), dst, src)
end

function baswap(dst, src)
    ccall((:baswap, libbluetooth), Cvoid, (Ptr{bdaddr_t}, Ptr{bdaddr_t}), dst, src)
end

function strtoba(str)
    ccall((:strtoba, libbluetooth), Ptr{bdaddr_t}, (Ptr{Cchar},), str)
end

function batostr(ba)
    ccall((:batostr, libbluetooth), Ptr{Cchar}, (Ptr{bdaddr_t},), ba)
end

function ba2str(ba, str)
    ccall((:ba2str, libbluetooth), Cint, (Ptr{bdaddr_t}, Ptr{Cchar}), ba, str)
end

function ba2strlc(ba, str)
    ccall((:ba2strlc, libbluetooth), Cint, (Ptr{bdaddr_t}, Ptr{Cchar}), ba, str)
end

function str2ba(str, ba)
    ccall((:str2ba, libbluetooth), Cint, (Ptr{Cchar}, Ptr{bdaddr_t}), str, ba)
end

function ba2oui(ba, oui)
    ccall((:ba2oui, libbluetooth), Cint, (Ptr{bdaddr_t}, Ptr{Cchar}), ba, oui)
end

function bachk(str)
    ccall((:bachk, libbluetooth), Cint, (Ptr{Cchar},), str)
end

function bt_malloc(size)
    ccall((:bt_malloc, libbluetooth), Ptr{Cvoid}, (Csize_t,), size)
end

function bt_free(ptr)
    ccall((:bt_free, libbluetooth), Cvoid, (Ptr{Cvoid},), ptr)
end

function bt_error(code)
    ccall((:bt_error, libbluetooth), Cint, (UInt16,), code)
end

function bt_compidtostr(id)
    ccall((:bt_compidtostr, libbluetooth), Ptr{Cchar}, (Cint,), id)
end

function bswap_128(src, dst)
    ccall((:bswap_128, libbluetooth), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), src, dst)
end

const BTPROTO_L2CAP = 0

const BTPROTO_HCI = 1

const BTPROTO_SCO = 2

const BTPROTO_RFCOMM = 3

const BTPROTO_BNEP = 4

const BTPROTO_CMTP = 5

const BTPROTO_HIDP = 6

const BTPROTO_AVDTP = 7

const SOL_HCI = 0

const SOL_L2CAP = 6

const SOL_SCO = 17

const SOL_RFCOMM = 18

const SOL_BLUETOOTH = 274

const BT_SECURITY = 4

const BT_SECURITY_SDP = 0

const BT_SECURITY_LOW = 1

const BT_SECURITY_MEDIUM = 2

const BT_SECURITY_HIGH = 3

const BT_SECURITY_FIPS = 4

const BT_DEFER_SETUP = 7

const BT_FLUSHABLE = 8

const BT_FLUSHABLE_OFF = 0

const BT_FLUSHABLE_ON = 1

const BT_POWER = 9

const BT_POWER_FORCE_ACTIVE_OFF = 0

const BT_POWER_FORCE_ACTIVE_ON = 1

const BT_CHANNEL_POLICY = 10

const BT_CHANNEL_POLICY_BREDR_ONLY = 0

const BT_CHANNEL_POLICY_BREDR_PREFERRED = 1

const BT_CHANNEL_POLICY_AMP_PREFERRED = 2

const BT_VOICE = 11

const BT_SNDMTU = 12

const BT_RCVMTU = 13

const BT_VOICE_TRANSPARENT = 0x0003

const BT_VOICE_CVSD_16BIT = 0x0060

const BT_PHY = 14

const BT_PHY_BR_1M_1SLOT = 0x00000001

const BT_PHY_BR_1M_3SLOT = 0x00000002

const BT_PHY_BR_1M_5SLOT = 0x00000004

const BT_PHY_EDR_2M_1SLOT = 0x00000008

const BT_PHY_EDR_2M_3SLOT = 0x00000010

const BT_PHY_EDR_2M_5SLOT = 0x00000020

const BT_PHY_EDR_3M_1SLOT = 0x00000040

const BT_PHY_EDR_3M_3SLOT = 0x00000080

const BT_PHY_EDR_3M_5SLOT = 0x00000100

const BT_PHY_LE_1M_TX = 0x00000200

const BT_PHY_LE_1M_RX = 0x00000400

const BT_PHY_LE_2M_TX = 0x00000800

const BT_PHY_LE_2M_RX = 0x00001000

const BT_PHY_LE_CODED_TX = 0x00002000

const BT_PHY_LE_CODED_RX = 0x00004000

const BDADDR_BREDR = 0x00

const BDADDR_LE_PUBLIC = 0x01

const BDADDR_LE_RANDOM = 0x02

end # module


text:
InputFolder = './Images/ObjColoc/';
OutputFolder = './Results/Images/ObjColoc_C01/';

@iA = '*_C01*.tif';

@fxm_mFilterObjSize [iA] > [M];
params.MinArea = 0;
params.MaxArea = Inf;
/endf

/show iA > M;
/keep M > tif;

ext: jl

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
function Dims_create!(nnodes::Integer, ndims::Integer, dims::MPIBuffertype{T}) where {T <: Integer}
    # int MPI_Dims_create(int nnodes, int ndims, int dims[])
    @mpichk ccall((:MPI_Dims_create, libmpi), Cint,
                  (Cint, Cint, Ptr{Cint}), nnodes, ndims, dims)
end

function Dims_create!(nnodes::Integer, dims::AbstractArray{T,N}) where {T<:Integer, N}
    cdims = Cint.(dims[:])
    ndims = length(cdims)
    Dims_create!(nnodes, ndims, cdims)
    dims[:] .= cdims
end

function Cart_create(comm_old::Comm, ndims::Integer, dims::MPIBuffertype{Cint}, periods::MPIBuffertype{Cint}, reorder)
    comm_cart = Comm()
    # int MPI_Cart_create(MPI_Comm comm_old, int ndims, const int dims[],
    #                     const int periods[], int reorder, MPI_Comm *comm_cart)
    @mpichk ccall((:MPI_Cart_create, libmpi), Cint,
                  (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Cint, Ptr{MPI_Comm}),
                  comm_old, ndims, dims, periods, reorder, comm_cart)
    if comm_cart.val != MPI_COMM_NULL
        refcount_inc()
        finalizer(free, comm_cart)
    end
    comm_cart
end

function Cart_create(comm_old::Comm, dims::AbstractArray{T,N}, periods::Array{T,N}, reorder) where T <: Integer where N
    cdims    = Cint.(dims[:])
    cperiods = Cint.(periods[:])
    ndims    = length(cdims)
    Cart_create(comm_old, ndims, cdims, cperiods, reorder)
end

"""
    Cart_get(comm::Comm, maxdims::Integer)

Obtain information on the Cartesian topology of dimension `maxdims` underlying the
communicator `comm`. This is specified by two `Cint` arrays of `maxdims` elements
for the number of processes and periodicity properties along each Cartesian dimension.
A third `Cint` array is returned, containing the Cartesian coordinates of the calling process.
"""
function Cart_get(comm::Comm, maxdims::Integer)
    # preallocate with nontrivial values
    dims    = Cint[-1 for i = 1:maxdims]
    periods = Cint[-1 for i = 1:maxdims]
    coords  = Cint[-1 for i = 1:maxdims]
    # int MPI_Cart_get(MPI_Comm comm, int maxdims, int dims[], int periods[], int coords[])
    @mpichk ccall((:MPI_Cart_get, libmpi), Cint,
                  (MPI_Comm, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}),
                  comm, maxdims, dims, periods, coords)
    return dims, periods, coords
end

"""
    Cartdim_get(comm::Comm)

Return number of dimensions of the Cartesian topology associated with the communicator `comm`.
"""
function Cartdim_get(comm::Comm)
    dims    = Cint[0]
    # int MPI_Cartdim_get(MPI_Comm comm, int *ndims)
    @mpichk ccall((:MPI_Cartdim_get, libmpi), Cint,
                  (MPI_Comm, Ptr{Cint}),
                  comm, dims)
    return Int(dims[1])
end

function Cart_coords!(comm::Comm, rank::Integer, maxdims::Integer, coords::MPIBuffertype{Cint})
    # int MPI_Cart_coords(MPI_Comm comm, int rank, int maxdims, int coords[])
    @mpichk ccall((:MPI_Cart_coords, libmpi), Cint,
                  (MPI_Comm, Cint, Cint, Ptr{Cint}),
                  comm, rank, maxdims, coords)
end

function Cart_coords(comm::Comm, maxdims::Integer)
    ccoords = Vector{Cint}(undef, maxdims)
    rank    = Comm_rank(comm)
    Cart_coords!(comm, rank, maxdims, ccoords)
    Int.(ccoords)
end

function Cart_shift(comm::Comm, direction::Integer, disp::Integer)
    rank_source = Ref{Cint}()
    rank_dest   = Ref{Cint}()
    # int MPI_Cart_shift(MPI_Comm comm, int direction, int disp,
    #                    int *rank_source, int *rank_dest)
    @mpichk ccall((:MPI_Cart_shift, libmpi), Cint,
                  (MPI_Comm, Cint, Cint, Ptr{Cint}, Ptr{Cint}),
                  comm, direction, disp, rank_source, rank_dest)
    Int(rank_source[]), Int(rank_dest[])
end

function Cart_sub(comm::Comm, remain_dims)
    comm_sub = Comm()
    remain_dims = [Cint(dim) for dim in remain_dims]
    # int MPI_Cart_sub(MPI_Comm comm, const int remain_dims[], MPI_Comm *comm_new)
    @mpichk ccall((:MPI_Cart_sub, libmpi), Cint,
                  (MPI_Comm, Ptr{Cint}, Ptr{MPI_Comm}),
                  comm, remain_dims, comm_sub)
    if comm_sub.val != MPI_COMM_NULL
        refcount_inc()
        finalizer(free, comm_sub)
    end
    comm_sub
end


text:
using Gaia, Tyndall, Jultran

Ts = 200:10:300
νs = 25:.1:40000

sunlight_νs, sunlight_Ws = sunlight()
println("sunlight done")
json = {
    :sunlight => {
        :values => sunlight_Ws,
        :νs => sunlight_νs
    },
    :νs => [νs[1],νs[end]],
    :ts => Ts,
    :data => Array(Dict, length(Ts))
}

earthlight_νs = 25:.1:10000
for (i,T) in enumerate(Ts)
    json[:data][i] = {
        :earthlight => {
            :values => [planck(ν, T) for ν in earthlight_νs],
            :νs => earthlight_νs
        }
    }

    for gas=[:co2,:ch4,:h2o] #,:n2o]
        κs = compute_κs(gas, νs, T, 1.e4)

        json[:data][i][gas] = {
            :values => nonzeros(κs),
            :νs => νs[findn(κs)[1]]
        }
    end
end

JSON.print(open("data/spectrum.json", "w"), json)


text:
using Distributions
using KernelEstimator

"""
Unpaired contingency table
"""
function contingencyTableUnpairedDiffPostMc(n1::Real, y1::Real, n2::Real, y2::Real)
    S = 10000
    rn1 = Beta(y1+1, n1-y1+1)
    θ1 = rand(rn1, S)
    rn2 = Beta(y2+1, n2-y2+1)
    θ2 = rand(rn2, S)
    diff = θ1-θ2
    den = kerneldensity(diff)
    (θ1, θ2, diff, den)
end


text:
@testset "Monomial compression test" begin

    R, (v1, v2, v3, v4) = Nemo.PolynomialRing(Nemo.QQ, ["v1", "v2", "v3", "v4"])
    tests = [
        v1 + v2 - 2,
        v4 + v4 * v1 - 3 * v2 * v4,
        (v1 + v2 + v2^2) * (v3 + v4),
        (v1 + v2^2) * (v3^3 + v4) - 7 * (v1 - 3 - v2) * (v3 - v4^2 - v3^2)
    ]

    for t in tests
        a, b = monomial_compress(t, [v1, v2])
        s = 0
        for (x, y) in zip(a, b)
            s += parent_ring_change(x, parent(t)) * parent_ring_change(y, parent(t))
        end
        @test s == t
    end

end


text:
"""
    default_name(x)

Return the lowercase of `nameof(x)` in `String`.
"""
default_name(x) = lowercase(string(nameof(x)))

"""
    default_name(x)

Return the lowercase of a given package name. It will
ignore the suffix if it ends with ".jl".
"""
function default_name(x::String)
    if endswith(x, ".jl")
        name = x[1:end-3]
    else
        name = x
    end
    return lowercase(name)
end

function docstring(x)
    return sprint(Base.Docs.doc(x); context = :color => true) do io, x
        show(io, MIME"text/plain"(), x)
    end
end

xcall(m::Module, name::Function, xs...; kwargs...) = xcall(m, nameof(name), xs...; kwargs...)
xcall(m::Module, name::Symbol, xs...; kwargs...) = xcall(GlobalRef(m, name), xs...; kwargs...)
xcall(name, xs...; kwargs...) = xcall(Parse, name, xs...; kwargs...)

function xcall(ref::GlobalRef, xs...; kwargs...)
    params = Expr(:parameters)
    for (key, value) in kwargs
        push!(params.args, Expr(:kw, key, value))
    end

    if isempty(kwargs)
        Expr(:call, ref, xs...)
    else
        return Expr(:call, ref, params, xs...)
    end
end

"""
    get_version(m::Module)

Get the version of a given module. It will try to find the version
of `Project.toml` if the given module is a project module. If fails,
it returns `v"0.0.0"`.
"""
function get_version(m::Module)
    # no version for scripts
    m === Main && return v"0.0.0"
    # project module
    path = pathof(m)
    if path !== nothing
        envpath = dirname(dirname(path))
        project = Pkg.Types.read_project(Pkg.Types.projectfile_path(envpath))
        if project.name == string(nameof(m)) && project.version !== nothing
            return project.version
        end
    end

    if hasproperty(m, :COMMAND_VERSION)
        return m.COMMAND_VERSION
    end
    return v"0.0.0"
end

"""
    cache_file([file=Base.PROGRAM_FILE])

Return the files that will be cached:

1. `cmd.jl`: a julia script that contains the generated CLI code.
2. `checksum`: a checksum file for checking if the generated CLI code matches the original file.
"""
function cachefile(file = Base.PROGRAM_FILE)
    dir = cachedir(file)
    return joinpath(dir, "cmd.jl"), joinpath(dir, "checksum")
end

"""
    cachedir([file=Base.PROGRAM_FILE])

Return the cache directory.
"""
function cachedir(file = Base.PROGRAM_FILE)
    name, _ = splitext(basename(file))
    dir = joinpath(dirname(file), "." * name * ".cmd")
    isabspath(file) || return joinpath(pwd(), dir)
    return dir
end

"""
    iscached([file=Base.PROGRAM_FILE])

Check if the given file is cached or not.
"""
function iscached(file = Base.PROGRAM_FILE)
    cache_file, crc = cachefile(file)
    isfile(crc) || return false
    isfile(cache_file) || return false
    if read(crc, String) == string(checksum(file), base = 16)
        return true
    end
    return false
end

# taken from Steven G Johnson
function checksum(filename, blocksize = 16384)
    crc = zero(UInt32)
    open(filename, "r") do f
        while !eof(f)
            crc = crc32c(read(f, blocksize), crc)
        end
    end
    return crc
end

"""
    create_cache(cmd[, file=Base.PROGRAM_FILE])

Create cache for given command `cmd` at file location `file`.
"""
function create_cache(cmd, file = Base.PROGRAM_FILE)
    isempty(file) && error("do not execute CLI script in -e, -E or REPL or disable_cache()")
    dir = cachedir(file)
    if !ispath(dir)
        mkpath(dir)
    end

    cache_file, crc = cachefile(file)
    write(cache_file, cmd)
    write(crc, string(checksum(file), base = 16))
    return
end


text:
using Pkg; Pkg.activate(first(splitdir(@__DIR__)))

using CUDA
using DLPack
using PyCall

hoomd = pyimport("hoomd")
md = pyimport("hoomd.md")
np = pyimport("numpy")

const kT = 0.596161
const dt = 0.02045

hoomd.context.initialize("--mode=gpu")
hoomd.option.set_num_threads(8)

# System Definition

## Set snapshot
snapshot = hoomd.data.make_snapshot(
    N = 14,
    box = hoomd.data.boxdim(Lx = 41, Ly = 41, Lz = 41),
    particle_types = ["C", "H"],
    bond_types = ["CC", "CH"],
    angle_types = ["CCC", "CCH", "HCH"],
    dihedral_types = ["CCCC", "HCCC", "HCCH"],
    pair_types = ["CCCC", "HCCC", "HCCH"],
    dtype = "double"
)

## Set particles types ids
let ids = PyArray(snapshot.particles."typeid")
    ids[[2:4; 6:7; 9:10; 12:14]] .= 1
end

## Compute initial positions
positions = [
    -2.990196   0.097881   0.000091;
    -2.634894  -0.911406   0.001002;
    -2.632173   0.601251  -0.873601;
    -4.060195   0.099327  -0.000736;
    -2.476854   0.823942   1.257436;
    -2.832157   1.833228   1.256526;
    -2.834877   0.320572   2.131128;
    -0.936856   0.821861   1.258628;
    -0.578833   1.325231   0.384935;
    -0.581553  -0.187426   1.259538;
    -0.423514   1.547922   2.515972;
    -0.781537   1.044552   3.389664;
    +0.646485   1.546476   2.516800;
    -0.778816   2.557208   2.515062;
]

reference_box_low_coords = [-22.206855  -19.677099  -19.241968]
box_low_coords = [-snapshot.box.Lx/2  -snapshot.box.Ly/2  -snapshot.box.Lz/2]
positions .+= box_low_coords .- reference_box_low_coords

PyArray(snapshot.particles."position") .= positions

## Set masses
let
    mC = 12.00
    mH = 1.008
    PyArray(snapshot.particles."mass") .= [
        mC, mH, mH, mH,
        mC, mH, mH,
        mC, mH, mH,
        mC, mH, mH, mH
    ]
end

## Set charges
let
    reference_charges = [
        -0.180000, 0.060000, 0.060000, 0.060000,
        -0.120000, 0.060000, 0.060000,
        -0.120000, 0.060000, 0.060000,
        -0.180000, 0.060000, 0.060000, 0.060000
    ]
    charge_conversion = 18.22262
    PyArray(snapshot.particles."charge") .= charge_conversion .* reference_charges
end

## Set bonds
snapshot.bonds.resize(13)

let bonds_ids = PyArray(snapshot.bonds."typeid")
    bonds_ids[[1:3; 5:6; 8:9; 11:13]] .= 1
end

PyArray(snapshot.bonds."group") .= [
    0   2;
    0   1;
    0   3;
    0   4;
    4   5;
    4   6;
    4   7;
    7   8;
    7   9;
    7   10;
    10  11;
    10  12;
    10  13;
]

## Set angles
snapshot.angles.resize(24)

let angles_ids = PyArray(snapshot.angles."typeid")
    angles_ids[[3; 5:8; 11:14; 17:21]] .= 1
    angles_ids[[1:2; 4; 10; 16; 22:24]] .= 2
end

PyArray(snapshot.angles."group") .= [
    1   0   2 ;
    2   0   3 ;
    2   0   4 ;
    1   0   3 ;
    1   0   4 ;
    3   0   4 ;
    0   4   5 ;
    0   4   6 ;
    0   4   7 ;
    5   4   6 ;
    5   4   7 ;
    6   4   7 ;
    4   7   8 ;
    4   7   9 ;
    4   7   10;
    8   7   9 ;
    8   7   10;
    9   7   10;
    7   10  11;
    7   10  12;
    7   10  13;
    11  10  12;
    11  10  13;
    12  10  13;
]

## Set dihedrals
snapshot.dihedrals.resize(27)

let dihedrals_ids = PyArray(snapshot.dihedrals."typeid")
    dihedrals_ids[[3; 6; 9:11; 15; 18:21]] .= 1
    dihedrals_ids[[1:2; 4:5; 7:8; 13:14; 16:17; 22:27]] .= 2
end

PyArray(snapshot.dihedrals."group") .= [
    2   0   4   5 ;
    2   0   4   6 ;
    2   0   4   7 ;
    1   0   4   5 ;
    1   0   4   6 ;
    1   0   4   7 ;
    3   0   4   5 ;
    3   0   4   6 ;
    3   0   4   7 ;
    0   4   7   8 ;
    0   4   7   9 ;
    0   4   7   10;
    5   4   7   8 ;
    5   4   7   9 ;
    5   4   7   10;
    6   4   7   8 ;
    6   4   7   9 ;
    6   4   7   10;
    4   7   10  11;
    4   7   10  12;
    4   7   10  13;
    8   7   10  11;
    8   7   10  12;
    8   7   10  13;
    9   7   10  11;
    9   7   10  12;
    9   7   10  13;
]

## Set pairs
snapshot.pairs.resize(27)

let pairs_ids = PyArray(snapshot.pairs."typeid")
    pairs_ids[2:11] .= 1
    pairs_ids[12:27] .= 2
end

PyArray(snapshot.pairs."group") .= [
    # CCCC
    0   10;
    # HCCC\n",
    0   8 ;
    0   9 ;
    5   10;
    6   10;
    1   7 ;
    2   7 ;
    3   7 ;
    11  4 ;
    12  4 ;
    13  4 ;
    # HCCH
    1   5 ;
    1   6 ;
    2   5 ;
    2   6 ;
    3   5 ;
    3   6 ;
    5   8 ;
    6   8 ;
    5   9 ;
    6   9 ;
    8   11;
    8   12;
    8   13;
    9   11;
    9   12;
    9   13;
]

## Materialize system from snapshot
system = hoomd.init.read_snapshot(snapshot)

# Interactions

## Neighbor list
nl_ex = md.nlist.cell()
nl_ex.reset_exclusions(exclusions = ["1-2", "1-3", "1-4"])

## Lennard-Jones interactions
lj = md.pair.lj(r_cut = 12.0, nlist = nl_ex)
lj.pair_coeff.set("C", "C", epsilon = 0.07, sigma = 3.55)
lj.pair_coeff.set("H", "H", epsilon = 0.03, sigma = 2.42)
lj.pair_coeff.set("C", "H", epsilon = sqrt(0.07 * 0.03), sigma = sqrt(3.55 * 2.42))

## Coulomb interactions
coulomb = md.charge.pppm(hoomd.group.charged(), nlist = nl_ex)
coulomb.set_params(Nx = 64, Ny = 64, Nz = 64, order = 6, rcut = 12.0)

## Harmonic interactions
harmonic = md.bond.harmonic()
harmonic.bond_coeff.set("CC", k = 2 * 268.0, r0 = 1.529)
harmonic.bond_coeff.set("CH", k = 2 * 340.0, r0 = 1.09)

## Angle potentials
angle = md.angle.harmonic()
angle.angle_coeff.set("CCC", k = 2 * 58.35, t0 = 112.7 * π / 180)
angle.angle_coeff.set("CCH", k = 2 * 37.5, t0 = 110.7 * π / 180)
angle.angle_coeff.set("HCH", k = 2 * 33.0, t0 = 107.8 * π / 180)

## Dihedral potentials
dihedral = md.dihedral.opls()
dihedral.dihedral_coeff.set("CCCC", k1 = 1.3, k2 = -0.05, k3 = 0.2, k4 = 0.0)
dihedral.dihedral_coeff.set("HCCC", k1 = 0.0, k2 = 0.0, k3 = 0.3, k4 = 0.0)
dihedral.dihedral_coeff.set("HCCH", k1 = 0.0, k2 = 0.0, k3 = 0.3, k4 = 0.0)

## Special Lennard-Jones interactions
lj_special_pairs = md.special_pair.lj()
lj_special_pairs.pair_coeff.set("CCCC", epsilon = 0.07, sigma = 3.55, r_cut = 12.0)
lj_special_pairs.pair_coeff.set("HCCH", epsilon = 0.03, sigma = 2.42, r_cut = 12.0)
lj_special_pairs.pair_coeff.set("HCCC",
    epsilon = sqrt(0.07 * 0.03), sigma = sqrt(3.55 * 2.42), r_cut = 12.0
)

## Special Coulomb interactions
coulomb_special_pairs = md.special_pair.coulomb()
coulomb_special_pairs.pair_coeff.set("CCCC", alpha = 0.5, r_cut = 12.0)
coulomb_special_pairs.pair_coeff.set("HCCC", alpha = 0.5, r_cut = 12.0)
coulomb_special_pairs.pair_coeff.set("HCCH", alpha = 0.5, r_cut = 12.0)

# Integration method

md.integrate.mode_standard(dt = dt)
integrator = md.integrate.nvt(group = hoomd.group.all(), kT = kT, tau = 100dt)
integrator.randomize_velocities(seed = 42)

# Advanced Sampling
include("SSAGES.jl")
#using .SSAGES

SSAGES.set_backend(:HOOMD)
simulation = hoomd.context.current
system_view = SSAGES.wrap(simulation)

ξ = SSAGES.DihedralAngle(1, 5, 8, 11)

grid = SSAGES.Grid(-π, π, 64, true)

sampler = SSAGES.ABF(grid, ξ, 100)(system_view)

hook = SSAGES.bind(simulation, sampler)

# Run simulation
hoomd.run(100_000)
hoomd.run(400_000)

# Data analysis
using Plots
using Statistics

θ = let #grid = sampler.grid
    Δ = only(grid.Δl)
    r = range(Δ / 2 + only(grid.r₀); step = Δ, length = only(grid.size))
    collect(r)
end

dA = let sm = sampler
    F = Array(vec(sm.ΣF) ./ max.(sm.hist, 1))
    dA = (F .- mean(F)) ./ std(F)
end

scatter(θ, (dA .- reverse(dA)) ./ 2; label = "")
scatter!(θ, dA; label = "")


text:
@testset "Orders in absolute number fields" begin
  println("NfOrd.jl")
  @time include("NfOrd/NfOrd.jl")

  println("Elem.jl")
  @time include("NfOrd/Elem.jl")
  println("Ideal.jl")
  @time include("NfOrd/Ideal.jl")
  println("FracIdl.jl")
  @time include("NfOrd/FracIdl.jl")
  println("ResidueRing.jl")
  @time include("NfOrd/ResidueRing.jl")
  println("Clgp.jl")
  @time include("NfOrd/Clgp.jl")
  println("RayClassGroup.jl")
  @time include("NfOrd/RayClassGroup.jl")
  println("ResidueRingMultGrp.jl")
  @time include("NfOrd/ResidueRingMultGrp.jl")
  println("Overorders.jl")
  @time include("NfOrd/Overorders.jl")
  println("LinearAlgebra.jl")
  @time include("NfOrd/LinearAlgebra.jl")
  #println("PicardGroup.jl")
  #@time include("NfOrd/PicardGroup.jl")
end


text:
global num_neur = 0
global num_prop = 0
global num_pool = 0

abstract Channel
abstract Prop
abstract Neuron
abstract Source
abstract NeuronPool
abstract Network
abstract NetworkS <: Network
abstract NetworkP <: Network
abstract Helper
abstract Puddle

if VERSION < v"0.5.0-"
    typealias RemoteChannel RemoteRef
end

immutable Pt3d
    x::Float64
    y::Float64
    z::Float64
    d::Float64
    arc::Float64 #normalized distance from 0 to end
end

Pt3d(pt::Pt3d,xyz::Array{Float64,1})=Pt3d(xyz[1],xyz[2],xyz[3],pt.d,pt.arc)

add(pt::Pt3d,x::Float64,y::Float64,z::Float64)=Pt3d(pt.x+x,pt.y+y,pt.z+z,pt.d,pt.arc)

immutable SegArea
    l::Float64
    r::Float64
    t::Float64
end

SegArea()=SegArea(50.0,50.0,100.0)

immutable SegRi
    l::Float64
    r::Float64
end

SegRi()=SegRi(0.0,0.0)

immutable Node
    area::SegArea
    ri::SegRi
    parent::Int64
    children::Array{Int64,1}
    internal::Bool
    pt3d::UnitRange{Int64}
end

Node(a::SegArea,r::SegRi,p::Int64,c::Array{Int64},m::UnitRange{Int64})=Node(a,r,p,c,true,m)

Node(p::Int64,c::Array{Int64,1},m::Int64)=Node(SegArea(),SegRi(),p,c,false,m:m)

Node(n::Node)=Node(n.area,n.ri,n.parent,n.children,n.internal,n.pt3d)

add_node(n::Neuron,nd::Node)=push!(n.nodes,nd)

type Prop0<:Prop
end

make_prop(a::Prop0,b::Int64)=Prop0()

function gen_prop_check(prop_array)

    global num_prop::Int64

    if method_exists(make_prop,(typeof(prop_array),Int))
    else
        num_prop+=1
        gen_prop(prop_array,num_prop)
    end

    nothing
end

function gen_prop(a,k::Int64)

    if issubtype(typeof(a),Tuple)
        num_fields=sum([length(fieldnames(a[i]))-1 for i=1:length(a)])
    else
        num_fields=sum(length(fieldnames(a)))
    end

    @eval begin
        type $(symbol("Prop_$k")) <: Prop
            p::Array{Float64,2}
        end

        make_prop(b::$(typeof(a)),n::Int64)=$(symbol("Prop_$k"))(zeros(Float64,$(num_fields),n))
    end

    c=make_prop(a,0)

    @eval begin
        make_prop(a::$(typeof(c)),n::Int64)=$(symbol("Prop_$k"))(zeros(Float64,$(num_fields),n))
    end

    nothing
end

immutable Section
    mtype::Int64
    pnode::UnitRange{Int64}
    child::Array{Int64,1}
    parent::Int64
    pt3d::Array{Pt3d,1}
    length::Float64
    prop::DataType
end

Section(s::Section,p::Int64)=Section(s.mtype,s.pnode,s.child,p,s.pt3d,s.length,s.prop)

Section(s::Section,p::DataType)=Section(s.mtype,s.pnode,s.child,s.parent,s.pt3d,s.length,p)

Section(s::Section,p::UnitRange{Int64})=Section(s.mtype,p,s.child,s.parent,s.pt3d,s.length,s.prop)

Section!(s::Array{Section,1},ind::Int64,p::Int64)=(s[ind]=Section(s[ind],p); nothing)

Section!(n::Neuron,ind::Int64,p::Int64)=Section!(n.secs,ind,p)

Section!(s::Array{Section,1},ind::Int64,p::UnitRange{Int64})=(s[ind]=Section(s[ind],p); nothing)

Section!(n::Neuron,ind::Int64,p::UnitRange{Int64})=Section!(n.secs,ind,p)

add_sec(neuron::Neuron, sec::Section)=push!(neuron.secs,sec)

type Extra_coeffs
    c::Array{Float64,1}
    inds::Array{Int64,1}
end

type Point <: Source
end

type Line <: Source
end

type Mixed <: Source
end

type Extracellular{S<:Source}
    xyz::Array{Float64,1}
    coeffs::Array{Extra_coeffs,1}
    v::Array{Float64,1}
end

Extracellular(xyz::Array{Float64,1})=Extracellular(Line(),xyz)

Extracellular(s::Source,xyz)=Extracellular{typeof(s)}(xyz,Array(Extra_coeffs,0),Array(Float64,0))

type Stim
    Is::Array{Float64,1}
    mtype::Int64
    neur::Int64
    node::Int64
    tstart::Float64
    tstop::Float64
end

Stim(Is::Float64,mtype::Int,neur::Int,node::Int,tstart::Float64,tstop::Float64)=Stim([Is],mtype,neur,node,tstart,tstop)

type Intracellular
    mtype::Int64
    neur::Int64
    node::Int64
    v::Array{Float64,1}
end

Intracellular(mtype::Int64,neur::Int64,node::Int64)=Intracellular(mtype,neur,node,Array(Float64,0))

type HelperS <: Helper
    flags::Array{Bool,1}
end

HelperS()=HelperS(falses(4))

type HelperP <: Helper
    flags::Array{Bool,1}
    dims::UnitRange{Int64}
    l::Int64
    i::Array{Intracellular,1}
    s::Array{Stim,1}
    e::Array{Extracellular,1}
end

HelperP()=HelperP(falses(4),0:0,0,Array(Intracellular,0),Array(Stim,0),Array(Extra_coeffs,0))

HelperP(dims::UnitRange{Int64},l::Int64)=HelperP(falses(4),dims,l,Array(Intracellular,0),Array(Stim,0),Array(Extra_coeffs,0))

type MyChannel <: AbstractChannel
    d::Array{Bool,1}
    cond_take::Condition    # waiting for data to become available
    MyChannel(l::Int64) = new(falses(l), Condition())
end

make_neuron()=nothing

function gen_neur_check(myprop)

    global num_neur::Int64

    if method_exists(make_neuron,(typeof(myprop),Array{Node,1},Array{Section,1},Array{Array{Int64,1}}))
    else
        num_neur+=1
        gen_neuron(myprop,num_neur)
    end

    nothing
end

function gen_neuron(prop,k::Int64)

    if issubtype(typeof(prop),Prop)
        myprop=(prop,prop,prop,prop)
    else
        myprop=prop
    end

    @eval begin
        type $(symbol("Neuron_$k")) <: Neuron
            soma::$(typeof(myprop[1]))
            axon::$(typeof(myprop[2]))
            dendrite::$(typeof(myprop[3]))
            apical::$(typeof(myprop[4]))
            secs::Array{Section,1}
            v::Array{Float64,1}
            a::Array{Float64,1}
            b::Array{Float64,1}
            d::Array{Float64,1}
            rhs::Array{Float64,1}
            Ra::Float64
            Cm::Float64
            dt::Float64
            nodes::Array{Node,1}
            i_vm::Array{Float64,1}
            divm::Array{Float64,1}
            diag_old::Array{Float64,1}
            internal_nodes::Array{Array{Int64,1},1}
            par::Array{Int64,1}
            v1::Array{Float64,1}
            i2::Array{Float64,1}
        end

        type $(symbol("Puddle_$k")) <: Puddle
            n::Array{$(symbol("Neuron_$k")),1}
            h::HelperP
            i::Int64
            imax::Int64
            r::RemoteChannel{MyChannel}
        end

        Puddle(n::Array{$(symbol("Neuron_$k")),1})=$(symbol("Puddle_$k"))(n,HelperP(),0,0,RemoteChannel(()->MyChannel(0),1))
        Puddle(n::Array{$(symbol("Neuron_$k")),1},h::HelperP)=$(symbol("Puddle_$k"))(n,h,0,0,RemoteChannel(()->MyChannel(0),1))
        Puddle(n::Array{$(symbol("Neuron_$k")),1},dims::UnitRange{Int64},l::Int64,t::Int64,mutex::RemoteChannel{MyChannel})=$(symbol("Puddle_$k"))(n,HelperP(dims,l),0,t,mutex)

        function make_neuron(prop::$(typeof(prop)),n::Array{Node,1},s::Array{Section,1},inter::Array{Array{Int64,1}})

            len=length(n)

            secs=deepcopy(s)
            nodes=deepcopy(n)
            inter_n=deepcopy(inter)

            (soma,axon,dendrite,apical)=make_props(prop,inter_n)

            n2 = $(symbol("Neuron_$k"))(soma,axon,dendrite,apical,secs,zeros(Float64,len),zeros(Float64,len),zeros(Float64,len),zeros(Float64,len),zeros(Float64,len),35.4,1.0,.025,nodes,zeros(Float64,len),zeros(Float64,len),zeros(Float64,len),inter_n,zeros(Float64,len),zeros(Float64,len),zeros(Float64,len))

            fillA!(n2)

            n2
        end
    end
end

add_inode(n::Neuron,i::Int64)=push!(n.internal_nodes[n.secs[i].mtype],length(n.nodes))

find_parents!(n::Neuron)=(for i=1:length(n.secs), j in n.secs[i].child; Section!(n,j,i); end)

function add(neuron::Neuron,prop_array...)

    num_arg=length(prop_array)

    if num_arg==1

        prop_array=prop_array[1]

        gen_prop_check(prop_array)

        myprop=make_prop(prop_array,0)
        gen_current(prop_array,myprop)

    elseif num_arg==4

        for i=1:4
            gen_prop_check(prop_array[i])
        end

        myprop=(make_prop(prop_array[1],0),make_prop(prop_array[2],0),make_prop(prop_array[3],0),make_prop(prop_array[4],0))

        for i=1:4
            gen_current(prop_array[i],myprop[i])
        end

    else
        #error checking
    end

    gen_neur_check(myprop)

    n1=make_neuron(myprop,neuron.nodes,neuron.secs,neuron.internal_nodes)

end

gen_neuron(Prop0(),0)

function make_props(prop::Prop,inter_n::Array{Array{Int64,1}})

    soma=make_prop(prop,length(inter_n[1]))
    axon=make_prop(prop,length(inter_n[2]))
    dendrite=make_prop(prop,length(inter_n[3]))
    apical=make_prop(prop,length(inter_n[4]))

    (soma,axon,dendrite,apical)
end

function make_props(prop::Tuple,inter_n::Array{Array{Int64,1}})

    soma=make_prop(prop[1],length(inter_n[1]))
    axon=make_prop(prop[2],length(inter_n[2]))
    dendrite=make_prop(prop[3],length(inter_n[3]))
    apical=make_prop(prop[4],length(inter_n[4]))

    (soma,axon,dendrite,apical)
end

function gen_pool_check(neur,par,ts)

    global num_pool::Int64

    if par==true
        if method_exists(make_ppool,(typeof(neur),))
        else
            gen_npool(neur,par)
            num_pool+=1
        end
        p=make_ppool(neur,ts)
    else
        if method_exists(make_spool,(typeof(neur),))
        else
            gen_npool(neur,par)
            num_pool+=1
        end
        p=make_spool(neur)
    end

    make_network(p,ts)
end

make_spool()=nothing
make_ppool()=nothing

findmyid(d)=parse(Int,split(string(d),"_")[end])

function findid(d)
    id=parse(Int,split(string(d),"_")[end])
    symbol("Puddle_$id")
end

function findin(d)
    id=parse(Int,split(string(d),"_")[end])
    "Neuron_$id"
end

function gen_npool(neur, par::Bool)

    global num_pool::Int64

    if issubtype(typeof(neur),Neuron)
        myfields=[:($(symbol("N_1"))::Array{($(typeof(neur))),1})]
    else
        typeinds=Array(DataType,0)
        mtypes=[typeof(neur[i]) for i=1:length(neur)]
        #Determine types of neurons
        for i=1:length(mtypes)
            if sum(mtypes[i].==typeinds)<1
                push!(typeinds,mtypes[i])
            end
        end

        if par==false
            myfields=[:($(symbol("N_$i"))::Array{($(typeinds[i])),1}) for i=1:length(typeinds)]
        else
            myfields=[:($(symbol("N_$i"))::DArray{($(findid(typeinds[i]))),1}) for i=1:length(typeinds)]
        end
    end

    if par==false

        @eval begin
            type $(symbol("Pool_$num_pool")) <: NeuronPool
                $(myfields...)
            end

            type $(symbol("Network_$num_pool")) <: NetworkS
                neur::$(symbol("Pool_$num_pool"))
                t::FloatRange{Float64}
                extra::Array{Extracellular,1}
                #extracellular Stimulation
                intra::Array{Intracellular,1}
                stim::Array{Stim,1}
                helper::HelperS
            end

            function make_spool(neur::$(typeof(neur)))

                if issubtype(typeof(neur),Neuron)
                    $(symbol("Pool_$num_pool"))([neur])
                else
                    inds=Array(Array{Int64,1},0)
                    mtypes=[typeof(neur[i]) for i=1:length(neur)]

                    for i=1:length(fieldnames($(symbol("Pool_$num_pool"))))
                        push!(inds,find(mtypes.==eltype(fieldtype($(symbol("Pool_$num_pool")),i))))
                    end

                    $(symbol("Pool_$num_pool"))([typeof(neur[inds[i]][1])[neur[inds[i]]...] for i=1:length(inds)]...)
                end
            end

            gen_net_func($(symbol("Pool_$num_pool")),"initialcon!")
            gen_net_func($(symbol("Pool_$num_pool")),"main")

            function make_network(p::$(symbol("Pool_$num_pool")),ts::Float64)
                $(symbol("Network_$num_pool"))(p,0.0:.025:ts,Array(Extracellular,0),Array(Intracellular,0),Array(Stim,0),HelperS())
            end

        end
    else
        @eval begin
            type $(symbol("Pool_$num_pool")) <: NeuronPool
                $(myfields...)
                c::Array{Array{UnitRange{Int64},1},1}
                t::Array{Int64,1}
                i::Array{Int64,1}
                s::Array{Int64,1}
            end

            type $(symbol("Network_$num_pool")) <: NetworkP
                neur::$(symbol("Pool_$num_pool"))
                t::FloatRange{Float64}
                extra::Array{Extracellular,1}
                #extracellular Stimulation
                intra::Array{Intracellular,1}
                stim::Array{Stim,1}
            end

            function make_ppool(neur::$(typeof(neur)),ts::Float64)

                inds=Array(Array{Int64,1},0)
                mtypes=[findmyid(typeof(neur[i])) for i=1:length(neur)]

                for i=1:length(fieldnames($(symbol("Pool_$num_pool"))))-4
                    nt=findmyid(eltype(fieldtype(($(symbol("Pool_$num_pool"))),i)))
                    push!(inds,find(mtypes .== nt))
                end

                c=[Array(UnitRange{Int64},0) for i=1:length(inds)]
                t=zeros(Int64,length(inds))
                a=[typeof(neur[inds[i][1]])[neur[inds[i]]...] for i=1:length(inds)]

                b=Array(Any,length(inds))

                w=length(workers())
                mutex=RemoteRef(()->MyChannel(w), 1)
                tl=length(0.0:.025:ts)

                for i=1:length(b)

                    if length(a[i])>w
                        dims=DistributedArrays.chunk_idxs(length(a[i]),w)[1]
                        l=length(a[i])
                        t[i]=l
                        c[i]=get_dims(dims)
                        mtype=typeof(Puddle([a[i][1]]))
                        b[i]=distribute(mtype[(k=j[1];Puddle(a[i][k],k,l,tl,mutex)) for j in dims])
                    else
                        l=length(a[i])
                        t[i]=l
                        c[i]=[1:l for i=1:1]
                        b[i]=distribute([Puddle(a[i],1:l,l,tl,mutex)])
                    end
                end

                $(symbol("Pool_$num_pool"))(b...,c,t,Array(Int64,0),Array(Int64,0))
            end

            gen_net_func_p($(symbol("Pool_$num_pool")),"initialcon!")
            gen_net_func_p($(symbol("Pool_$num_pool")),"main")
            gen_net_func_p($(symbol("Pool_$num_pool")),"reset_t!")

            function make_network(p::$(symbol("Pool_$num_pool")),ts::Float64)
                $(symbol("Network_$num_pool"))(p,0.0:.025:ts,Array(Extracellular,0),Array(Intracellular,0),Array(Stim,0))
            end
        end
    end
end

Network(n::Neuron,ts::Float64; par=false)=(gen_pool_check(n,par,ts))

Network(n::Tuple,ts::Float64; par=false)=(gen_pool_check(n,par,ts))

function gen_net_func(n::DataType,func::ASCIIString)

    a=length(fieldnames(n))

    @eval begin
        function $(symbol("$func"))(n::$(n))
            for j = 1 : $a
                for k=1:length(getfield(n,j))
                    $(symbol("$func"))(getfield(n,j)[k])
                end
            end
        end
    end
    nothing
end

function gen_net_func_p(n::DataType,func::ASCIIString)

    a=length(fieldnames(n))-4

    @eval begin
        function $(symbol("$func"))(n::$(n))
            for j = 1 : $a
                @sync for p in procs(getfield(n,j))
                    @async remotecall_wait((n)->($(symbol("$func"))(localpart(n))), p, getfield(n,j))
                end
            end
        end
    end
    nothing
end

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
"""
    treynor_ratio(asset_returns, benchmark_returns; multiplier=1.0, risk_free=0.0)

Calculates the Treynor ratio as the ratio of excess return divided by the CAPM beta. This ratio is similar to the Sharpe Ratio, but instead of dividing by the volatility, we devide by the CAPM beta as risk proxy.

# Formula

    TR = E[asset_returns - risk_free] / beta * multiplier

# Arguments
- `asset_returns`:      Vector of asset returns.
- `benchmark_returns`:  Vector of benchmark returns (e.g. market portfolio returns).
- `multiplier`:         Optional scalar multiplier, i.e. use `12` to annualize monthly returns, and use `252` to annualize daily returns.
- `risk_free`:          Optional vector or scalar value denoting the risk-free return (must have same frequency as the provided returns, e.g. daily).
"""
function treynor_ratio(asset_returns, benchmark_returns; multiplier=1.0, risk_free=0.0)
    α, β = capm(asset_returns, benchmark_returns; risk_free=risk_free)
    mean(asset_returns .- risk_free) / β * multiplier
end


text:
function reset!(maze::Maze)
	for i in 1 : maze.h
		for j in 1 : maze.w
			if get_mazecolor(maze, (j, i)) != 0
				set_mazecolor(maze, (j, i), 0xffffffff)
			end
		end
	end
end

function clutter!(maze::Maze, p::Float64)
	for i in 2 : maze.h - 1
		for j in 2 : maze.w - 1
			if rand(Float64) < p
				if get_mazecolor(maze, (j, i)) == 0
					set_mazecolor(maze, (j, i), 0xffffffff)
				else
					set_mazecolor(maze, (j, i), 0xff000000)
				end
			end
		end
	end
end

function make_random!(maze::Maze, p::Float64)
	for i in 2 : maze.h - 1
		for j in 2 : maze.w - 1
			if rand(Float64) < p
				set_mazecolor(maze, (j, i), 0xff000000)
			else
				set_mazecolor(maze, (j, i), 0xffffffff)
			end
		end
	end
end


text:
@testset "SparseMaskMatrix" begin
    using OptEnrichedSetCover: SparseMaskMatrix

    @testset "empty" begin
        @inferred SparseMaskMatrix()
        mask = SparseMaskMatrix()
        @test size(mask) == (0, 0)
        @test size(mask, 1) == 0
        @test size(mask, 2) == 0
        @test_throws ArgumentError size(mask, 0)
        @test_throws ArgumentError size(mask, 3)

        empty2 = SparseMaskMatrix(2, 3)
        @test size(empty2) == (2, 3)
        @test view(empty2, :, 1) == Int[]
        @test view(empty2, :, 2) == Int[]
        @test view(empty2, :, 3) == Int[]
        @test mask != empty2
        @test empty2 == SparseMaskMatrix(2, 3)
    end

    @testset "1x1" begin
        @inferred SparseMaskMatrix(1, 1, [1, 1], Int[])
        empty = SparseMaskMatrix(1, 1, [1, 1], Int[])
        @test size(empty) == (1, 1)
        @test size(empty, 1) == 1
        @test size(empty, 2) == 1
        @test empty[:, 1] == Int[]
        @test view(empty, :, 1) == Int[]

        nonempty = SparseMaskMatrix(1, 1, [1, 2], Int[1])
        @test size(nonempty) == (1, 1)
        @test size(nonempty, 1) == 1
        @test size(nonempty, 2) == 1
        @test nonempty[:, 1] == [1]
        @test view(nonempty, :, 1) == [1]
        @test nonempty != empty
    end

    @testset "SparseMaskMatrix(Collection{Set}, elm2ix)" begin
        sm = SparseMaskMatrix([Set([:a, :b]), Set{Symbol}(), Set([:c, :d]), Set([:a, :c])],
                              Dict(:a=>1, :b=>2, :c=>3, :d=>4, :e=>5))
        @test size(sm) == (5, 4)
        @test view(sm, :, 1) == [1, 2]
        @test view(sm, :, 2) == Int[]
        @test view(sm, :, 3) == [3, 4]
        @test view(sm, :, 4) == [1, 3]
    end

    @testset "SparseMaskMatrix(Vector{Vector{Int}})" begin
        sm = SparseMaskMatrix(4, [[1,2], Int[], [3,4]])
        @test size(sm) == (4, 3)
        @test view(sm, :, 1) == [1, 2]
        @test view(sm, :, 2) == Int[]
        @test view(sm, :, 3) == [3, 4]
    end

    @testset "SparseMaskMatrix(Matrix{Bool})" begin
        sm = SparseMaskMatrix([false false true; true false false])
        @test size(sm) == (2, 3)
        @test view(sm, :, 1) == [2]
        @test view(sm, :, 2) == Int[]
        @test view(sm, :, 3) == [1]
        @test sm == SparseMaskMatrix(2, [[2], Int[], [1]])
        @test sm != SparseMaskMatrix(3, [[2], Int[], [1]])
        @test sm != SparseMaskMatrix(2, [[1], Int[], [1]])
        @test sm != SparseMaskMatrix(2, [[2], [1], Int[]])
    end

    @testset "convert(Matrix{Bool}, SparseMaskMatrix)" begin
        sm1 = SparseMaskMatrix()
        @inferred convert(Matrix{Bool}, sm1)
        mtx1 = convert(Matrix{Bool}, sm1)
        @test mtx1 isa Matrix{Bool}
        @test mtx1 == Matrix{Bool}(undef, 0, 0)

        sm2 = SparseMaskMatrix([false false true; true false false])
        mtx2 = convert(Matrix{Bool}, sm2)
        @test mtx2 == [false false true; true false false]

        mtx3 = convert(Matrix, sm2)
        @test mtx3 isa Matrix{Bool}
        @test mtx3 == mtx2
    end

    @testset "getindex(SparseMaskMatrix)" begin
        sm = SparseMaskMatrix([false false true; true false false])

        @testset "SMS[:, ?]" begin
            @inferred sm[:, 1]
            @test sm[:, 1] == [2]
            @test sm[:, 2] == Int[]

            @inferred sm[:, [1, 2]]
            sm1 = sm[:, [1, 2]]
            @test sm1 isa SparseMaskMatrix
            @test convert(Matrix, sm1) == [false false; true false]
            @test size(sm1) == (2, 2)
            @inferred sm[:, [false, true, true]]
            @test sm[:, [2, 3]] == sm[:, [false, true, true]]
            @test sm[:, [2, 1, 3, 1]] == SparseMaskMatrix([false false true false; false true false true])
        end

        @testset "SMS[?, :]" begin
            @test_broken sm[1, :] == [false false true] # not implemented yet
            @test sm[[1], :] == SparseMaskMatrix([false false true])
            @inferred sm[[2,1,2], :]
            @test sm[[2,1,2], :] == SparseMaskMatrix([true false false; false false true; true false false])
        end

        @testset "SMS[?, ?]" begin
            @inferred sm[1, 2]
            @test sm[1, 1] == false
            @test sm[1, 3] == true
            @inferred sm[[2,1,2], [1]]
            @test sm[[2,1,2], [1]] == SparseMaskMatrix(reshape([true, false, true], 3, 1))
        end
    end

    @testset "permutedims(SMS, ?)" begin
        @test transpose(SparseMaskMatrix()) == SparseMaskMatrix()
        @test transpose(SparseMaskMatrix(2, 3)) == SparseMaskMatrix(3, 2)

        sm = SparseMaskMatrix([false false true; true false false])
        @inferred permutedims(sm, [1, 2])
        @test permutedims(sm, [1, 2]) == sm
        @inferred permutedims(sm, [2, 1])
        tsm = permutedims(sm, [2, 1])
        @test size(tsm) == (3, 2)
        @test tsm == SparseMaskMatrix([false true; false false; true false])
        @test transpose(sm) == tsm
        @test transpose(transpose(sm)) == sm
    end
end


text:
using Contour, Gadfly, DataFrames

using OneDimensional: goldenSectionSearch, fibSearch

f(x, y) = x^4 + y^4 - 5 * (x * y - x^2 * y^2)
fdx(x, y) = 4x^3 - 5y + 10 * x * y^2
fdy(x, y) = 4y^3 - 5x + 10 * y * x^2
▽(x, y) = [fdx(x, y), fdy(x, y)]
L(x, y, λ) = [x, y] .- (λ .* ▽(x, y))

function simpleGradientDescent(f, ▽, x0, y0, ɛ) # down gradient with exponentital slow step
  callCnt, gradCallCnt, xs, ys = 1, 0, [x0], [y0]

  λ = 0.05
  fprev = f(x0, y0)
  while true
    xk, yk = xs[end], ys[end]
    xnew, ynew = L(xk, yk, λ)
    gradCallCnt += 1
    if xnew < 0 || 1 < xnew || ynew < 0 || 1 < ynew
      λ *= 0.75
      continue
    end
    fnew = f(xnew, ynew)
    callCnt += 1
    if fnew > fprev
      λ *= 0.75
      continue
    end
    fprev = fnew
    push!(xs, xnew)
    push!(ys, ynew)
    if norm([xs[end] - xs[end - 1], ys[end] - ys[end - 1]]) < ɛ
      break
    end
  end
  xs, ys, callCnt, gradCallCnt
end

function fastGradienDescent(f, ▽, x0, y0, ɛ) # fastest descent gradient
  callCnt, gradCallCnt, xl, yl = 1, 0, [x0], [y0]

  left, right = 0, 1
  while true
    x_left, y_left = xl[end], yl[end]
    g = ▽(x_left, y_left)
    gradCallCnt += 1
    ax = [(left - x_left) / -g[1], (right - x_left) / -g[1]]
    ay = [(left - y_left) / -g[2], (right - y_left) / -g[2]]
    sort!(ax)
    sort!(ay)
    a = [max(ax[1], ay[1]), min(ax[2], ay[2])]
    λ, cnt, _, _, _ = fibSearch(arg -> f(x_left - arg * g[1], y_left - arg * g[2]), a[1], a[2], ɛ)
    callCnt += cnt
    xnew, ynew = x_left - λ * g[1], y_left - λ * g[2]
    push!(xl, xnew)
    push!(yl, ynew)
    if norm([xl[end] - xl[end - 1], yl[end] - yl[end - 1]]) < ɛ
      break
    end
  end
  xl, yl, callCnt, gradCallCnt
end

ɛ = 0.0001
x0, y0 = 0.2, 0.8
sgx, sgy, cnt1, gcnt1 = simpleGradientDescent(f, ▽, x0, y0, ɛ)
fgx, fgy, cnt2, gcnt2 = fastGradienDescent(f, ▽, x0, y0, ɛ)

@printf("Обычный градиентный спуск: %d вычислений функции, %d вычислений градиента\n", cnt1, gcnt1)
@printf("Наискорейший спуск: %d вычислений функции, %d вычислений градиента\n", cnt2, gcnt2)

Gadfly.set_default_plot_size(18cm, 18cm)
plot(layer(x=[0:0.01:1], y=[0:0.01:1], z=f, Geom.contour(levels=15)),
     layer(x=sgx, y=sgy, Theme(default_color=color("orange")), Geom.point, Geom.line),
     layer(x=fgx, y=fgy, Theme(default_color=color("red")), Geom.point, Geom.line(preserve_order=true)),
     Guide.xticks(ticks=[0:0.2:1.0]), Guide.yticks(ticks=[0:0.2:1.0]),
     Guide.title("Линии уровня функции x^4 + y^4 - 5(xy - x^2 * y^2)"))

max_stars_repo_name: ChShersh/university-courses

max_issues_repo_name: ChShersh/university-courses

max_forks_repo_name: ChShersh/university-courses


text:
# Helper function to modify the fairness tensor according to the values for sp2n, on2p, etc
# vals is a 2D array of the form [[sp2n, sn2p], [op2n, on2p]]
function _fairTensorLinProg(ft::FairTensor, vals)
	newftmat = Array{Union{VariableRef, Int, GenericAffExpr{Float64,VariableRef}}, 3}(undef, size(ft.mat)...)
	for i in 1:length(ft.labels)
		p2n, n2p = vals[i, :] #These vals are VariableRef from library JuMP
		p2p, n2n = 1-p2n, 1-n2p
		a = Array{Union{VariableRef, Int, GenericAffExpr{Float64,VariableRef}}, 2}(undef, 2, 2) # The numbers for modified fairness tensor values for a group
		a[1, 1] = ft.mat[i, 1, 1]*p2p + ft.mat[i, 2, 2]*n2p
		a[1, 2] = ft.mat[i, 1, 2]*p2p + ft.mat[i, 2, 1]*n2p
		a[2, 1] = ft.mat[i, 2, 1]*n2n + ft.mat[i, 1, 1]*p2n
		a[2, 2] = ft.mat[i, 2, 2]*n2n + ft.mat[i, 1, 2]*p2n
		newftmat[i, :, :] = a
	end
	return FairTensor(newftmat, ft.labels)
end

"""
    LinProgWrapper

It is a postprocessing algorithm that uses JuMP and Ipopt library to minimise error and satisfy the equality of specified specified measures for all groups at the same time.
Automatic differentiation and gradient based optimisation is used to find probabilities with which the predictions are changed for each group.
"""
struct LinProgWrapper{M<:MLJBase.Model} <: DeterministicComposite
	grp::Symbol
	classifier::M
	measures::Array{<:Measure}
end

"""
    LinProgWrapper(classifier=nothing, grp=:class, measure=nothing, measures=nothing)

Instantiates LinProgWrapper which wraps the classifier and containts the measure to optimised and the sensitive attribute(grp)
You can optimize the all fairness metrics in measures. You can optimize for only a single metric using keyword measure.
"""
function LinProgWrapper(; classifier::MLJBase.Model=nothing, grp::Symbol=:class, measure=nothing, measures=nothing)
	if measures==nothing measures=[measure] end
	model = LinProgWrapper(grp, classifier, measures)
	message = MLJBase.clean!(model)
	isempty(message) || @warn message
	return model
end

function MLJBase.clean!(model::LinProgWrapper)
    warning = ""
	model.measures[1]!=nothing || (warning *= "No Debiasing Measure specified\n")
	model.classifier!=nothing || (warning *= "No classifier specified in model\n")
    target_scitype(model) <: AbstractVector{<:Finite} || (warning *= "Only Binary Classifiers are supported\n")
    return warning
end

function MMI.fit(model::LinProgWrapper, verbosity::Int, X, y)
	grps = X[:, model.grp]
	n = length(levels(grps)) # Number of different values for sensitive attribute

	# As LinProgWrapper is a postprocessing algorithm, the model needs to be fitted first
	mch = machine(model.classifier, X, y)
	fit!(mch)
	ŷ = MMI.predict(mch, X)

	if typeof(ŷ[1]) <: MLJBase.UnivariateFinite
		ŷ = MLJBase.mode.(ŷ)
	end

	labels = levels(y)
	favLabel = labels[2]
	unfavLabel = labels[1]
	y = y.==favLabel
	ŷ = ŷ.==favLabel

	# Finding the probabilities of changing predictions is a Linear Programming Problem
	# JuMP and Ipopt Optimizer are used to for this Linear Programming Problem
	m = JuMP.Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level"=>0))

	@variable(m, 0<= p2p[1:n] <=1)
	@variable(m, 0<= p2n[1:n] <=1)
	@variable(m, 0<= n2p[1:n] <=1)
	@variable(m, 0<= n2n[1:n] <=1)

	@constraint(m, [i=1:n], p2p[i] == 1 - p2n[i])
	@constraint(m, [i=1:n], n2p[i] == 1 - n2n[i])

	ft = fair_tensor(categorical(ŷ), categorical(y), categorical(grps))

	vals = Array{Union{VariableRef, Int, GenericAffExpr{Float64,VariableRef}}, 2}(undef, n, 2)
	vals[: , 1] = p2n
	vals[: , 2] = n2p

	ft = _fairTensorLinProg(ft, vals)

	mat = reshape(ft.mat, (4n))
	@variable(m, aux[1:4n])
	@constraint(m,[i=1:4n], mat[i]==aux[i])

	register(m, :fpr, 4n, (x...)->fpr(Fairness.FairTensor{n}(reshape(collect(x), (n, 2, 2)), ft.labels)), autodiff=true)
	register(m, :fnr, 4n, (x...)->fnr(Fairness.FairTensor{n}(reshape(collect(x), (n, 2, 2)), ft.labels)), autodiff=true)
	@NLobjective(m, Min, fpr(aux...) + fnr(aux...))

	for j in 1:length(model.measures)
		measure = model.measures[j]
		register(m, Symbol("func$j"), 4n, (x...)->measure(Fairness.FairTensor{n}(reshape(collect(x), (n, 2, 2)), ft.labels), grp=levels(grps)[1]), autodiff=true)
		for i in 2:n
			fn_symbol = Symbol("func$(j)_$(i)")
			register(m, fn_symbol, 4n, (x...)->measure(Fairness.FairTensor{n}(reshape(collect(x), (n, 2, 2)), ft.labels), grp=levels(grps)[i]), autodiff=true)
			JuMP.add_NL_constraint(m, :($(Expr(:call, fn_symbol, aux...))==$(Expr(:call, Symbol("func$j"), aux...))))
			# TODO: Replace call to func1 with a pre-computed expression
		end
	end
	optimize!(m)

	# fitresult will provide the info we require in the predict function.
	fitresult = [[JuMP.value.(p2n), JuMP.value.(n2p)], mch.fitresult, labels]
	# Note: It was necessary to return the levels(y) value in fitreult because in predict there
	# is no way to infer the 2 possible values of labels/targets.
	# Main reason to return values is the edge case : Maybe all of the ŷ predictions are same and we don't get to know both labels

	return fitresult, nothing, nothing
end

# Corresponds to eq_odds function which uses mix_rates to modify results
function MMI.predict(model::LinProgWrapper, fitresult, Xnew)

	(p2n, n2p), classifier_fitresult, labels = fitresult

	ŷ = MMI.predict(model.classifier, classifier_fitresult, Xnew)

	if typeof(ŷ[1]) <: MLJBase.UnivariateFinite
		ŷ = MLJBase.mode.(ŷ)
	end

	favLabel = labels[2]
	unfavLabel = labels[1]

	grps = Xnew[:, model.grp]

	n = length(levels(grps)) # Number of different values for sensitive attribute

	for i in 1:n
		Class = levels(grps)[i]
		Grp = grps .== Class

		pp_indices = shuffle(findall((grps.==Grp) .& (ŷ.==favLabel))) # predicted positive for iᵗʰ class
		pn_indices = shuffle(findall((grps.==Class) .& (ŷ.==unfavLabel))) # predicted negative for iᵗʰ class

		# Note : arrays in julia start from 1
		p2n_indices = pp_indices[1:convert(Int, floor(length(pp_indices)*p2n[i]))]
		n2p_indices = pn_indices[1:convert(Int, floor(length(pn_indices)*n2p[i]))]

		ŷ[p2n_indices] .= unfavLabel
		ŷ[n2p_indices] .= favLabel
	end
	return ŷ
end

MMI.input_scitype(::Type{<:LinProgWrapper{M}}) where M = input_scitype(M)
MMI.target_scitype(::Type{<:LinProgWrapper{M}}) where M = AbstractVector{<:Finite{2}}


text:
abstract type AbstractSherlog <: AbstractFloat end

for SherlogN in (:Sherlog16, :Sherlog32, :Sherlog64)

@eval FloatN = $(Symbol("Float",string(SherlogN)[end-1:end]))

@eval begin

struct $SherlogN{T<:AbstractFloat,i} <: AbstractSherlog    # T is the bitpattern type for logging, i the logbook id
    val::$FloatN
end

# conversions back from Sherlog
Base.Float64(x::$SherlogN) = Float64(x.val)
Base.Float32(x::$SherlogN) = Float32(x.val)
Base.Float16(x::$SherlogN) = Float16(x.val)
Base.Int64(x::$SherlogN) = Int64(x.val)
Base.Int32(x::$SherlogN) = Int32(x.val)
Base.Int16(x::$SherlogN) = Int16(x.val)

# otherwise conversion to itself is ambigious
$SherlogN{T,i}(x::$SherlogN) where {T<:AbstractFloat,i} = $SherlogN{T,i}(x.val)
$SherlogN(x::$SherlogN{T,i}) where {T<:AbstractFloat,i} = $SherlogN{T,i}(x.val)

# generator functions from Float, Int and Bool
$SherlogN(x::AbstractFloat) = $SherlogN{Float16,1}(x)   # use Float16 and logbook #1 as standard
$SherlogN{T}(x::AbstractFloat) where T = $SherlogN{T,1}(x)      # use logbook #1 if no logbook provided
$SherlogN(x::Integer) = $SherlogN{Float16,1}(x)           # same for integers
$SherlogN{T}(x::Integer) where T = $SherlogN{T,1}(x)
$SherlogN{T}(x::Bool) where T = if x $SherlogN{T,1}(1) else $SherlogN{T,1}(0) end

Base.promote_rule(::Type{<:Integer},::Type{$SherlogN{T,i}}) where {T,i} = $SherlogN{T,i}
Base.promote_rule(::Type{Float64},::Type{$SherlogN{T,i}}) where {T,i} = $SherlogN{T,i}
Base.promote_rule(::Type{Float32},::Type{$SherlogN{T,i}}) where {T,i} = $SherlogN{T,i}
Base.promote_rule(::Type{Float16},::Type{$SherlogN{T,i}}) where {T,i} = $SherlogN{T,i}

Base.bitstring(x::$SherlogN) = bitstring(x.val)
Base.show(io::IO,x::$SherlogN) = print(io,$SherlogN,"(",string(x.val),")")

Base.eps(::Type{$SherlogN{T,i}}) where {T,i} = eps(Float64)
Base.eps(x::$SherlogN) = eps(x.val)

Base.typemin(::$SherlogN{T,i}) where {T,i} = $SherlogN{T,i}(typemin($FloatN))
Base.typemax(::$SherlogN{T,i}) where {T,i} = $SherlogN{T,i}(typemax($FloatN))
Base.floatmin(::$SherlogN{T,i}) where {T,i} = $SherlogN{T,i}(floatmin($FloatN))
Base.floatmax(::$SherlogN{T,i}) where {T,i} = $SherlogN{T,i}(floatmax($FloatN))
Base.precision(::$SherlogN{T,i}) where {T,i} = $SherlogN{T,i}(precision($FloatN))

end

for O in (:(+), :(-), :(*), :(/), :(^))
    @eval function Base.$O(x::$SherlogN{T,i},y::$SherlogN{T,i}) where {T,i}
        r = $O(x.val, y.val)
        log_it(T,r,i)
        return $SherlogN{T,i}(r)
    end
end

for O in (:(-), :(+),
          :sign,
          :prevfloat, :nextfloat,
          :round, :trunc, :ceil, :floor,
          :inv, :abs, :sqrt, :cbrt,
          :exp, :expm1, :exp2, :exp10,
          :log, :log1p, :log2, :log10,
          :rad2deg, :deg2rad, :mod2pi, :rem2pi,
          :sin, :cos, :tan, :csc, :sec, :cot,
          :asin, :acos, :atan, :acsc, :asec, :acot,
          :sinh, :cosh, :tanh, :csch, :sech, :coth,
          :asinh, :acosh, :atanh, :acsch, :asech, :acoth,
          :sinc, :sinpi, :cospi,
          :sind, :cosd, :tand, :cscd, :secd, :cotd,
          :asind, :acosd, :atand, :acscd, :asecd, :acotd
         )
    @eval function Base.$O(x::$SherlogN{T,i}) where {T,i}
        r = $O(x.val)
        log_it(T,r,i)
        return $SherlogN{T,i}(r)
    end
end

for O in (:(<), :(<=))
    @eval Base.$O(x::$SherlogN, y::$SherlogN) = $O(x.val, y.val)
end

@eval Random.rand(rng::Random.AbstractRNG, ::Random.Sampler{T}) where T<:$SherlogN= convert(T, rand(rng, $FloatN))
@eval Random.randn(rng::Random.AbstractRNG, ::Type{T}) where T<:$SherlogN = convert(T, randn(rng, $FloatN))
@eval Random.randexp(rng::Random.AbstractRNG, ::Type{T}) where T<:$SherlogN = convert(T, randexp(rng, $FloatN))

@eval round(::Type{T}, x::$SherlogN) where {T<:Integer} = round(T, x.val)
@eval rem(x::$SherlogN, y::$SherlogN) = Core.Intrinsics.rem_float(x.val, y.val)

end


text:
function filefromid(ID::Integer,DataDir)
    tarr = []
    filearr = readdir("$DataDir")
    for i in filearr
        fi = stat("$DataDir/$i")
        datei = Dates.unix2datetime(fi.mtime)
        hh = parse(Int,Dates.format(datei,"HH"))+1
        ti = parse(Int,Dates.format(datei,"yymmdd")*"$hh"*Dates.format(datei,"MM"))
        #println("$i : $ti")
        push!(tarr,ti)
    end

    i = findfirst(tarr,ID)

    if i==0
        error("No sample found")
    end

    FNAME = "$DataDir/$(filearr[i])"

    return FNAME
end

# ID = 1801161626
# FNAME = filefromid(ID)

max_stars_repo_name: akels/WeylDisc

max_issues_repo_name: akels/WeylDisc

max_forks_repo_name: akels/WeylDisc


text:
import Statistics
using LinearAlgebra
using OrderedCollections
using ChainRulesCore
using ChainRules
using NNlib
using Ghost
using Ghost: Tape, Variable, V, Call, mkcall, Constant, inputs
using Ghost: bound, compile, play!, isstruct
using Ghost: remove_first_parameter, kwfunc_signature, call_signature

include("helpers.jl")
include("drules.jl")
include("chainrules.jl")
include("grad.jl")
include("update.jl")
include("gradcheck.jl")

function __init__()
    update_chainrules_primitives!()
end


text:
# Quadratic discriminant
module QuadDiscClassifier
	export QuadDiscModel, QuadDisc_train, QuadDisc_exec
	using StatsFuns: softmax!

	# QuadDiscModel struct
	struct QuadDiscModel
		mean::Vector{Vector{Float64}}		# mean vector
		icm::Array{Float64}			# inverse covariance matrices for each class
		dcm::Vector{Float64}			# determinants of the covariance matrices
		priors::Dict{Int,Float64}		# class priors
		classes::Vector{Int}			# classes
		r1::Float64				# diagonal regularization weight
		r2::Float64				# trace regularization weight
	end

	# Printers
	Base.show(io::IO, m::QuadDiscModel) = print(io, "Quadratic Discriminant model, $(length(m.mean)) classes, r1=$(m.r1), r2=$(m.r2)")

	# Train methods (classification)
	QuadDisc_train(x::Matrix{Float64}, y::Vector{Int}, r1::Float64=0, r2::Float64=0) = begin

		m = size(x,1) 							# number of variables
		classes = sort(unique(y))					# classes
		C = length(classes)

		priors = Dict(yi=>sum(yi.==y)/length(y) for yi in classes)	# priors
		mv = [vec(mean(x[:,y.==yi],2)) for yi in classes]		# mean vectors

		cv = zeros(C,m,m)
		dcm = zeros(C)
		@inbounds for (i,yi) in enumerate(classes) # iterate through the classes

			# Calculate class covariance
			cv[i,:,:] = cov(x[:,y.==yi],2)

			# Calculate the log of the determinant of the cov matrix (for the fixed term)
			dcm[i] = -1/2*log(det(cv[i,:,:]))

			# Regularize and invert the covariance matrices
			cv[i,:,:] = inv( (1-r1-r2)*cv[i,:,:] + r1*Matrix(Diagonal(diag(cv[i,:,:]))) + r2*trace(cv[i,:,:])*eye(m) )

		end

		return QuadDiscModel(mv, cv, dcm, priors, classes, r1, r2)
	end

	# Execution methods (classification)
	QuadDisc_exec(m::QuadDiscModel, x::Matrix{Float64})=begin

		# Function for calculating the quadratic term: -1/2*xᵀ*inv(Σᵢ)*x (or diag(x'*(-1/2*m.icm[c,:,:])*x)', very slow)
		xTwx!(out,x,w) = begin
			for i in 1:size(x,2)
				v=view(x,:,i)
				@simd for c in 1:size(w,1)
					@inbounds out[c,i] *=v'*view(w,c,:,:)*v
				end
			end
			return out
		end

		# Initializations
		C = length(m.classes)
		out=fill(-0.5, size(m.icm,1), size(x,2))
		xTwx!(out, x, m.icm)

		@inbounds @simd for c in 1:C
					# Linear term: inv(Σᵢ)*μᵀ*x			# Constant term log(P(ωᵢ)) - ...
			out[c:c,:] += (m.icm[c,:,:]*m.mean[c])'*x .+ log(m.priors[m.classes[c]]).- 1/2*m.mean[c]'*m.icm[c,:,:]*m.mean[c] .+ m.dcm[c]
		end

		# Apply softmax
		@inbounds @simd for j in 1:size(out,2)
			softmax!(view(out,:,j))
		end

		return out
	end
end

##########################
# FunctionCell Interface #
##########################
"""
	quaddisc(r1=0.0, r2=0.0)

Constructs an untrained cell that when piped data inside, returns a quadratic liscriminant classifier trained
function cell based on the input data and labels.

# Arguments
  * `r1`::Float64 represents the diagonal covariance matrix regularization parameter
  * `r2`::Float64 represents the covariance matrix trace regularization parameter

For more information:
	[1] L. Kuncheva "Combining Pattern Classifiers 2'nd Ed." 2014, ISBN 978-1-118-31523-1
"""
quaddisc(r1::Float64=0.0, r2::Float64=0.0) = FunctionCell(quaddisc, (r1,r2), ModelProperties(), "Quadratic discriminant classifier: r1=$r1, r2=$r2")

############################
# DataCell/Array Interface #
############################
"""
	quaddisc(x, r1=0.0, r2=0.0)

Trains a quadratic discriminant classifcation model that using the data `x` and regularization parameters `r1`, `r2`.
"""
# Training
quaddisc(x::T where T<:CellDataL, r1::Float64=0.0, r2::Float64=0.0) = quaddisc((getx!(x), gety(x)), r1, r2)
quaddisc(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, r1::Float64=0.0, r2::Float64=0.0) = quaddisc((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), r1, r2)
quaddisc(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, r1::Float64=0.0, r2::Float64=0.0) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[quaddisc] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Transform labels first
	enc = labelencn(x[2])
	yenc = label2ind.(x[2],enc)

	# Train model
	qddata = QuadDiscClassifier.QuadDisc_train(getobs(x[1]), yenc, r1, r2)

	# Build model properties
	modelprops = ModelProperties(nvars(x[1]), length(enc.label), enc)

	FunctionCell(quaddisc, Model(qddata, modelprops), "Quadratic discriminant classifier: r1=$r1, r2=$r2")

end

# Execution
quaddisc(x::T where T<:CellData, model::Model{<:QuadDiscClassifier.QuadDiscModel}) =
	datacell(quaddisc(getx!(x), model), gety(x))
quaddisc(x::T where T<:AbstractVector, model::Model{<:QuadDiscClassifier.QuadDiscModel}) =
	quaddisc(mat(x, LearnBase.ObsDim.Constant{2}()), model)
quaddisc(x::T where T<:AbstractMatrix, model::Model{<:QuadDiscClassifier.QuadDiscModel}) =
	QuadDiscClassifier.QuadDisc_exec(model.data, getobs(x))


text:
# ============================================================
# Structs of Solution
# Solver stores data in structs of arrays (SoA)
# ============================================================

struct Solution{T1,T2,ND} <: AbstractSolution
    w::T1
    prim::T1
    ∇w::T2
end

struct Solution1F{T1,T2,T3,T4,ND} <: AbstractSolution
    w::T1
    prim::T1
    ∇w::T2
    f::T3
    ∇f::T4
end

struct Solution2F{T1,T2,T3,T4,ND} <: AbstractSolution
    w::T1
    prim::T1
    ∇w::T2
    h::T3
    b::T3
    ∇h::T4
    ∇b::T4
end

function Solution1D(w::AA, prim::AA)
    sw = deepcopy(w)
    if w[1] isa Number
        sw .= 0.0
    else
        for e in sw
            e .= 0.0
        end
    end

    return Solution{typeof(w),typeof(sw),1}(w, prim, sw)
end

function Solution1D(w::AA, prim::AA, f::AA)
    sw = deepcopy(w)
    sf = deepcopy(f)
    if w[1] isa Number
        sw .= 0.0
        sf .= 0.0
    else
        for e in sw
            e .= 0.0
        end
        for e in sf
            e .= 0.0
        end
    end

    return Solution1F{typeof(w),typeof(sw),typeof(f),typeof(sf),1}(w, prim, sw, f, sf)
end

function Solution1D(w::AA, prim::AA, h::AA, b::AA)
    sw = deepcopy(w)
    sh = deepcopy(h)
    sb = deepcopy(b)
    if w[1] isa Number
        sw .= 0.0
        sh .= 0.0
        sb .= 0.0
    else
        for e in sw
            e .= 0.0
        end
        for e in sh
            e .= 0.0
        end
        for e in sb
            e .= 0.0
        end
    end

    return Solution2F{typeof(w),typeof(sw),typeof(h),typeof(sh),1}(
        w,
        prim,
        sw,
        h,
        b,
        sh,
        sb,
    )
end

function Solution2D(w::AA, prim::AA)
    ET = ifelse(w[1] isa Number, eltype(w), eltype(w[1]))
    sw = zeros(ET, axes(w)[1], 2, axes(w)[2:end]...)
    return Solution{typeof(w),typeof(sw),2}(w, prim, sw)
end

function Solution2D(w::AA, prim::AA, f::AA)
    ET = ifelse(w[1] isa Number, eltype(w), eltype(w[1]))
    sw = zeros(ET, axes(w)[1], 2, axes(w)[2:end]...)
    sf = begin
        if ndims(f) - ndims(w) > 1
            zeros(ET, axes(f)[1:2]..., 2, axes(f)[3:end]...)
        else
            zeros(ET, axes(f)[1], 2, axes(f)[2:end]...)
        end
    end
    return Solution1F{typeof(w),typeof(sw),typeof(f),typeof(sf),2}(w, prim, sw, f, sf)
end

function Solution2D(w::AA, prim::AA, h::AA, b::AA)
    ET = ifelse(w[1] isa Number, eltype(w), eltype(w[1]))
    sw = zeros(ET, axes(w)[1], 2, axes(w)[2:end]...)
    sh = begin
        if ndims(h) - ndims(w) > 1
            zeros(ET, axes(h)[1:2]..., 2, axes(h)[3:end]...)
        else
            zeros(ET, axes(h)[1], 2, axes(h)[2:end]...)
        end
    end
    sb = zero(sh)
    return Solution2F{typeof(w),typeof(sw),typeof(h),typeof(sh),2}(
        w,
        prim,
        sw,
        h,
        b,
        sh,
        sb,
    )
end


text:
using Test

@testset "BoteSalvatICX" begin
    include("xione.jl")
end


text:
using Oceananigans

include("ConvergenceTests/ConvergenceTests.jl")

using .ConvergenceTests.ForcedFlowFreeSlip: setup_and_run_xy, setup_and_run_xz

# Run 4 simulations:
Nx = [32, 64, 128, 256]
stop_time = 0.01

h = π / maximum(Nx)
Δt = 0.01 * h^2
stop_iteration = round(Int, stop_time / Δt)
Δt = stop_time / stop_iteration

for N in Nx
    setup_and_run_xy(Nx=N, Δt=Δt, stop_iteration=stop_iteration)
    setup_and_run_xz(Nx=N, Δt=Δt, stop_iteration=stop_iteration)
end


text:
using DelimitedFiles

function read_matrix(dataset::String)::OBOPDataset
    C = readdlm(dataset,Int64)
    total_voters, total_itens = size(C)
    item_names = [string(i) for i in 1:total_itens]
    D = [C[i,j] - C[j,i] for i in 1:total_itens, j in 1:total_itens]
    interval_init = [i * (1/total_itens) for i in 0:total_itens-1]
    instance = OBOPDataset(dataset,
                           0,
                           total_itens,
                           item_names,
                           C,
                           D,
                           interval_init)
    return instance
end

function read_soc(dataset::String)::OBOPDataset
    instance = OBOPDataset(dataset,
                           0,
                           0,
                           Array{String}(undef,0),
                           Array{Int64}(undef,0,0),
                           Array{Int64}(undef,0,0),
                           Array{Float64}(undef,0))
    open(dataset) do file
        # Load total of itens
        instance.total_itens = parse(Int64,readline(file))
        # Load itens names
        for i in 1:instance.total_itens
            line = split(readline(file),",")
            push!(instance.item_names,strip(line[2]))
        end
        # Load total of ranks
        instance.total_voters = parse(Int64,split(readline(file),",")[3])
        # Create cumulative score matrix
        instance.C = zeros(Int64,instance.total_itens,instance.total_itens)
        found, bucketID = false, 1
        rank = Array{Int64}(undef, instance.total_itens)
        for r in 1:instance.total_voters
            line = parse.(Int64,split(readline(file),','))
            total = line[begin]
            rank = line[begin+1:end]
            # Update C matrix
            for i in 1:instance.total_itens-1
                for j in i+1:instance.total_itens
                    instance.C[rank[i],rank[j]] += total
                    instance.C[rank[j],rank[i]] -= total
                end
            end
        end
        # Create D matrix
        instance.D = [instance.C[i,j] - instance.C[j,i] for i in 1:instance.total_itens, j in 1:instance.total_itens]
        # Create interval vector
        instance.interval_init = [i * (1/instance.total_itens) for i in 0:instance.total_itens-1]
    end
    return instance
end

function read_soi(dataset::String)::OBOPDataset
    instance = OBOPDataset(dataset,
                           0,
                           0,
                           Array{String}(undef,0),
                           Array{Int64}(undef,0,0),
                           Array{Int64}(undef,0,0),
                           Array{Float64}(undef,0))
    open(dataset) do file
        # Load total of itens
        instance.total_itens = parse(Int64,readline(file))
        # Load itens names
        for i in 1:instance.total_itens
            line = split(readline(file),",")
            push!(instance.item_names,strip(line[2]))
        end
        # Load total of ranks
        instance.total_voters = parse(Int64,split(readline(file),",")[3])
        # Create cumulative score matrix
        instance.C = zeros(Int64,instance.total_itens,instance.total_itens)
        found, bucketID = false, 1
        rank = Array{Int64}(undef, instance.total_itens)
        for r in 1:instance.total_voters
            line = parse.(Int64,split(readline(file),','))
            total = line[begin]
            rank = line[begin+1:end]
            # Update C matrix
            for i in 1:size(rank,1)-1
                for j in i+1:size(rank,1)
                    instance.C[rank[i],rank[j]] += total
                    instance.C[rank[j],rank[i]] -= total
                end
            end
        end
        # Create D matrix
        instance.D = [instance.C[i,j] - instance.C[j,i] for i in 1:instance.total_itens, j in 1:instance.total_itens]
        # Create interval vector
        instance.interval_init = [i * (1/instance.total_itens) for i in 0:instance.total_itens-1]
    end
    return instance
end

function read_toc(dataset::String)::OBOPDataset
    instance = OBOPDataset(dataset,
                           0,
                           0,
                           Array{String}(undef,0),
                           Array{Int64}(undef,0,0),
                           Array{Int64}(undef,0,0),
                           Array{Float64}(undef,0))
    open(dataset) do file
        # Load total of itens
        instance.total_itens = parse(Int64,readline(file))
        # Load itens names
        for i in 1:instance.total_itens
            line = split(readline(file),",")
            push!(instance.item_names,strip(line[2]))
        end
        # Load total of ranks
        instance.total_voters = parse(Int64,split(readline(file),",")[3])
        # Create cumulative score matrix
        instance.C = zeros(Int64,instance.total_itens,instance.total_itens)
        rank = Array{Int64}(undef, instance.total_itens)
        found, bucketID = false, 1
        for r in 1:instance.total_voters
            fill!(rank,0)
            line = split(readline(file),',')
            total = parse(Int64,line[begin])
            line = line[begin+1:end]
            found, bucketID = false, 1
            for l in line
                if l[begin] != '{' && l[end] != '}'
                    value = parse(Int64,l)
                    rank[value] = bucketID
                    if found == false
                        bucketID += 1
                    end
                else
                    if l[begin] == '{'
                        found = true
                        value = parse(Int64,l[begin+1:end])
                        rank[value] = bucketID
                    elseif l[end] == '}'
                        found = false
                        value = parse(Int64,l[begin:end-1])
                        rank[value] = bucketID
                        bucketID += 1
                    end
                end
            end
            # Update C matrix
            for i in 1:instance.total_itens-1
                for j in i+1:instance.total_itens
                    if rank[i] == rank[j]
                        instance.C[i,j] += total
                        instance.C[j,i] += total
                    elseif rank[i] < rank[j]
                        instance.C[i,j] += total
                        instance.C[j,i] -= total
                    else
                        instance.C[i,j] -= total
                        instance.C[j,i] += total
                    end
                end
            end
        end
        # Create D matrix
        instance.D = [instance.C[i,j] - instance.C[j,i] for i in 1:instance.total_itens, j in 1:instance.total_itens]
        # Create interval vector
        instance.interval_init = [i * (1/instance.total_itens) for i in 0:instance.total_itens-1]
    end
    return instance
end

function read_toi(dataset::String)::OBOPDataset
    instance = OBOPDataset(dataset,
                           0,
                           0,
                           Array{String}(undef,0),
                           Array{Int64}(undef,0,0),
                           Array{Int64}(undef,0,0),
                           Array{Float64}(undef,0))
    open(dataset) do file
        # Load total of itens
        instance.total_itens = parse(Int64,readline(file))
        # Load itens names
        for i in 1:instance.total_itens
            line = split(readline(file),",")
            push!(instance.item_names,strip(line[2]))
        end
        # Load total of ranks
        instance.total_voters = parse(Int64,split(readline(file),",")[3])
        # Create cumulative score matrix
        instance.C = zeros(Int64,instance.total_itens,instance.total_itens)
        rank = Array{Int64}(undef, instance.total_itens)
        found, bucketID = false, 1
        for r in 1:instance.total_voters
            fill!(rank,0)
            line = split(readline(file),',')
            total = parse(Int64,line[begin])
            line = line[begin+1:end]
            found, bucketID = false, 1
            for l in line
                if l[begin] != '{' && l[end] != '}'
                    value = parse(Int64,l)
                    rank[value] = bucketID
                    if found == false
                        bucketID += 1
                    end
                else
                    if l[begin] == '{'
                        found = true
                        value = parse(Int64,l[begin+1:end])
                        rank[value] = bucketID
                    elseif l[end] == '}'
                        found = false
                        value = parse(Int64,l[begin:end-1])
                        rank[value] = bucketID
                        bucketID += 1
                    end
                end
            end
            # Update C matrix
            for i in 1:instance.total_itens-1
                for j in i+1:instance.total_itens
                    if rank[i] != 0 && rank[j] != 0
                        if rank[i] == rank[j]
                            instance.C[i,j] += total
                            instance.C[j,i] += total
                        elseif rank[i] < rank[j]
                            instance.C[i,j] += total
                            instance.C[j,i] -= total
                        else
                            instance.C[i,j] -= total
                            instance.C[j,i] += total
                        end
                    end
                end
            end
        end
        # Create D matrix
        instance.D = [instance.C[i,j] - instance.C[j,i] for i in 1:instance.total_itens, j in 1:instance.total_itens]
        # Create interval vector
        instance.interval_init = [i * (1/instance.total_itens) for i in 0:instance.total_itens-1]
    end
    return instance
end

"""
    read_dataset(dataset[, cost_matrix])

Read the `dataset` that can be either in Preflib.org formats or a pre-calculated Cumulative Cost Matrix.
If `cost_matrix` is defined as `true` it will use the cost matrix as input, else it will infer the file type considering the
following Prelib formats:

- `.soc` (Strict Order - Complete List)
- `.soi` (Strict Order - Incomplete List)
- `.toc` (Order with Ties - Complete List)
- `.toi` (Order with Ties - Incomplete List)
"""
function read_dataset(dataset::String; cost_matrix::Bool=false)::OBOPDataset
    instance = nothing
    if cost_matrix
        instance = read_matrix(dataset)
    else
        file_type = dataset[end-2:end]
        if file_type == "soc"
            instance = read_soc(dataset)
        elseif file_type == "soi"
            instance = read_soi(dataset)
        elseif file_type == "toc"
            instance = read_toc(dataset)
        elseif file_type == "toi"
            instance = read_toi(dataset)
        end
    end
    return instance
end

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop