text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

mutable struct CoreClient <: AbstractCoreClient
    channel::gRPCClient.gRPCChannel
    address::AbstractString
    port::Int
    core_stub::Core_TypeDBStub
    databaseMgr::CoreDatabaseManager
    sessions::Dict{Bytes, Optional{AbstractCoreSession}}
end

function Base.show(io::IO, core_client::CoreClient)
    print(io, "CoreClient(address: $(core_client.address):$(core_client.port))")
end

function CoreClient(address::AbstractString, port::Int = 1729)
    channel = gRPCChannel(address * ":" * string(port))
    stub = Core_TypeDBStub(channel)
    databaseMgr = CoreDatabaseManager()
    sessions = Dict{String, Union{<:AbstractCoreSession, Nothing}}()
    return CoreClient(channel,address, port, stub, databaseMgr,sessions)
end

is_cluster(client::AbstractCoreClient) = false

function remove_session(client::AbstractCoreClient, session::AbstractCoreSession)
    delete!(client.sessions, session.sessionID)
end

function close(client::AbstractCoreClient)
    safe_close.(collect(values(client.sessions)))
    return nothing
end


text:
using ProgressMeter

"""
    reduce_read_positions(sel::Array{<:Integer})

Reduces the individual read positions by finding sub-ranges.
Returns an array of read indices and an array with the number of entries per index.
"""
function reduce_read_positions(sel::Array{<:Integer})

    Npart     = size(sel,1)
    index     = zeros(Int64, Npart)
    n_to_read = zeros(Int64, Npart)

    index[1]     = sel[1] - 1
    n_to_read[1] = 1
    Nentries     = 1

    @inbounds for  i=2:Npart
        if (sel[i] == index[Nentries] + n_to_read[Nentries] + 1)
            n_to_read[Nentries] += 1
        else
            Nentries            += 1
            index[Nentries]     = sel[i] - 1
            n_to_read[Nentries] += 1
        end
    end

    resize!(index, Nentries)
    resize!(n_to_read, Nentries)

    return index, n_to_read
end

"""
    find_read_positions( snap_base::String, filter_function::Function)

Finds the number of particles and their storage location in a snapshot directory that pass the `filter_function`.
"""
function find_read_positions( snap_base::String, filter_function::Function;
                              verbose::Bool=true)

    if verbose
        @info "Getting number of particles to read..."
        t1 = time_ns()
    end

    # read the header of the zero'th file
    h = read_header(select_file(snap_base, 0))

    # number of particles that fulfill the filter criteria
    N_part = 0

    # storage dictionary
    d = Dict()

    for sub_snap = 0:(h.num_files-1)
        snap_file   = select_file(snap_base, sub_snap)
        sel         = filter_function(snap_file)
        N_this_file = size(sel,1)

        # store read positions of particles are in the file
        if N_this_file > 0
            index, n_to_read = reduce_read_positions(sel)

            d[sub_snap] = Dict( "index" => index, "n_to_read" => n_to_read)
        end

        if verbose
            @info "sub-snap $sub_snap: $N_this_file particles"
        end
        N_part   += N_this_file
    end

    if verbose
        t2 = time_ns()
        @info "Need to read $N_part particles"
        @info "  elapsed: $(output_time(t1,t2)) s"
    end

    d["N_part"] = N_part

    return d

end

"""
    find_read_positions( snap_base::String, geometry::AbstractGadgetGeometry;
                         parttype::Integer=0,
                         verbose::Bool=true )

Finds the number of particles and their storage location in a snapshot directory that are contained in the space defined by `geometry`.
"""
function find_read_positions( snap_base::String, geometry::AbstractGadgetGeometry;
                              parttype::Integer=0,
                              verbose::Bool=true)

    if verbose
        @info "Getting number of particles to read..."
        t1 = time_ns()
    end

    # read the header of the zero'th file
    h = read_header(select_file(snap_base, 0))

    # number of particles that fulfill the filter criteria
    N_part = 0

    # storage dictionary
    d = Dict()

    for sub_snap = 0:(h.num_files-1)
        snap_file = select_file(snap_base, sub_snap)

        # read header
        h_snap    = read_header(snap_file)

        # read the position block
        pos       = read_block(snap_file, "POS"; parttype)

        # find positions of particles that are contained in the geometry
        sel         = get_geometry_mask(geometry, pos)
        N_this_file = size(sel,1)

        # store read positions of particles are in the file
        if N_this_file > 0
            index, n_to_read = reduce_read_positions(sel)

            d[sub_snap] = Dict( "index" => index, "n_to_read" => n_to_read)
        end

        if verbose
            @info "sub-snap $sub_snap: $N_this_file / $(h_snap.npart[parttype+1]) particles"
        end
        N_part   += N_this_file
    end

    if verbose
        t2 = time_ns()
        N_all = Int64(h.nall[parttype+1])
        N_per = N_part / N_all * 100.0
        @info "Need to read $N_part / $N_all particles -> $(@sprintf("%0.2f", N_per)) %"
        @info "  elapsed: $(output_time(t1,t2)) s"
    end

    d["N_part"] = N_part

    return d

end

"""
    save_read_positions(read_positions_file::String, data)

Saves the relevant read-in positions to a binary file.
"""
function save_read_positions(read_positions_file::String, data)

    f = open(read_positions_file, "w")
    N_part = data["N_part"]

    delete!(data, "N_part")
    files = keys(data)
    N_files = length(keys(data))

    write(f, Int64(N_part))
    write(f, Int64(N_files))

    for file ∈ files
        N_entries = size(data[file]["index"],1)
        write(f, Int64(file))
        write(f, Int64(N_entries))
        write(f, Int64.(data[file]["index"]))
        write(f, Int64.(data[file]["n_to_read"]))
    end

    close(f)

end

"""
    load_read_positions(read_positions_file::String)

Loads the relevant read-in positions from a binary file.
"""
function load_read_positions(read_positions_file::String)

    read_positions = Dict()

    f = open(read_positions_file, "r")

    read_positions["N_part"] = read(f, Int64)

    N_files = read(f, Int64)

    for _ ∈ 1:N_files

        file      = read(f, Int64)
        N_entries = read(f, Int64)
        index     = read!(f, Array{Int64,1}(undef, N_entries))
        n_to_read = read!(f, Array{Int64,1}(undef, N_entries))

        read_positions[file] = Dict("index" => index, "n_to_read" => n_to_read)

    end

    close(f)

    return read_positions
end

ext: jl


text:
"""
#### Read the particle-data
- select variables
- limit to a maximum range
- print the name of each data-file before reading it
- toggle verbose mode

```julia
getparticles(       dataobject::InfoType;
                    lmax::Real=dataobject.levelmax,
                    vars::Array{Symbol,1}=[:all],
                    stars::Bool=true,
                    xrange::Array{<:Any,1}=[missing, missing],
                    yrange::Array{<:Any,1}=[missing, missing],
                    zrange::Array{<:Any,1}=[missing, missing],
                    center::Array{<:Any,1}=[0., 0., 0.],
                    range_unit::Symbol=:standard,
                    print_filenames::Bool=false,
                    verbose::Bool=verbose_mode)
```
#### Returns an object of type PartDataType, containing the particle-data table, the selected and the simulation ScaleType and summary of the InfoType
```julia
return PartDataType()

# get an overview of the returned fields:
# e.g.:
julia> info = getinfo(100)
julia> particles  = getparticles(info)
julia> viewfields(particles)
#or:
julia> fieldnames(particles)
```

#### Arguments
##### Required:
- **`dataobject`:** needs to be of type: "InfoType", created by the function *getinfo*
##### Predefined/Optional Keywords:
- **`lmax`:** not defined
- **`stars`:** not defined
- **`var(s)`:** the selected particle variables in arbitrary order: :all (default), :cpu, :mass, :vx, :vy, :vz, :birth :metals, ...
- **`xrange`:** the range between [xmin, xmax] in units given by argument `range_unit` and relative to the given `center`; zero length for xmin=xmax=0. is converted to maximum possible length
- **`yrange`:** the range between [ymin, ymax] in units given by argument `range_unit` and relative to the given `center`; zero length for ymin=ymax=0. is converted to maximum possible length
- **`zrange`:** the range between [zmin, zmax] in units given by argument `range_unit` and relative to the given `center`; zero length for zmin=zmax=0. is converted to maximum possible length
- **`range_unit`:** the units of the given ranges: :standard (code units), :Mpc, :kpc, :pc, :mpc, :ly, :au , :km, :cm (of typye Symbol) ..etc. ; see for defined length-scales viewfields(info.scale)
- **`center`:** in units given by argument `range_unit`; by default [0., 0., 0.]; the box-center can be selected by e.g. [:bc], [:boxcenter], [value, :bc, :bc], etc..
- **`print_filenames`:** print on screen the current processed particle file of each CPU
- **`verbose`:** print timestamp, selected vars and ranges on screen; default: set by the variable `verbose_mode`

### Defined Methods - function defined for different arguments
- getparticles( dataobject::InfoType; ...) # no given variables -> all variables loaded
- getparticles( dataobject::InfoType, var::Symbol; ...) # one given variable -> no array needed
- getparticles( dataobject::InfoType, vars::Array{Symbol,1}; ...)  # several given variables -> array needed

#### Examples
```julia
# read simulation information
julia> info = getinfo(420)

# Example 1:
# read particle data of all variables, full-box, all levels
julia> particles = getparticles(info)

# Example 2:
# read particle data of all variables
# data range 20x20x4 kpc; ranges are given in kpc relative to the box (here: 48 kpc) center at 24 kpc
julia> particles = getparticles( info,
                                  xrange=[-10., 10.],
                                  yrange=[-10., 10.],
                                  zrange=[-2., 2.],
                                  center=[24., 24., 24.],
                                  range_unit=:kpc )

# Example 3:
# give the center of the box by simply passing: center = [:bc] or center = [:boxcenter]
# this is equivalent to center=[24.,24.,24.] in Example 2
# the following combination is also possible: e.g. center=[:bc, 12., 34.], etc.
julia> particles = getparticles(    info,
                                    xrange=[-10.,10.],
                                    yrange=[-10.,10.],
                                    zrange=[-2.,2.],
                                    center=[33., bc:, 10.],
                                    range_unit=:kpc )

# Example 4:
# read particle data of the variables mass and the birth-time, full-box, all levels
julia> particles = getparticles( info, [:mass, :birth] ) # use array for the variables

# Example 5:
# read particle data of the single variable mass, full-box, all levels
julia> particles = getparticles( info, :mass ) # no array for a single variable needed
...
```

"""
function getparticles( dataobject::InfoType, var::Symbol;
                    lmax::Real=dataobject.levelmax,
                    stars::Bool=true,
                    xrange::Array{<:Any,1}=[missing, missing],
                    yrange::Array{<:Any,1}=[missing, missing],
                    zrange::Array{<:Any,1}=[missing, missing],
                    center::Array{<:Any,1}=[0., 0., 0.],
                    range_unit::Symbol=:standard,
                    print_filenames::Bool=false,
                    verbose::Bool=verbose_mode)

    return  getparticles( dataobject, vars=[var],
                        lmax=lmax,
                        stars=stars,
                        xrange=xrange,
                        yrange=yrange,
                        zrange=zrange,
                        center=center,
                        range_unit=range_unit,
                        print_filenames=print_filenames,
                        verbose=verbose)
end

function getparticles( dataobject::InfoType, vars::Array{Symbol,1};
                    lmax::Real=dataobject.levelmax,
                    stars::Bool=true,
                    xrange::Array{<:Any,1}=[missing, missing],
                    yrange::Array{<:Any,1}=[missing, missing],
                    zrange::Array{<:Any,1}=[missing, missing],
                    center::Array{<:Any,1}=[0., 0., 0.],
                    range_unit::Symbol=:standard,
                    print_filenames::Bool=false,
                    verbose::Bool=verbose_mode)

    return  getparticles( dataobject, vars=vars,
                                        lmax=lmax,
                                        stars=stars,
                                        xrange=xrange,
                                        yrange=yrange,
                                        zrange=zrange,
                                        center=center,
                                        range_unit=range_unit,
                                        print_filenames=print_filenames,
                                        verbose=verbose)
end

function getparticles( dataobject::InfoType;
                    lmax::Real=dataobject.levelmax,
                    vars::Array{Symbol,1}=[:all],
                    stars::Bool=true,
                    xrange::Array{<:Any,1}=[missing, missing],
                    yrange::Array{<:Any,1}=[missing, missing],
                    zrange::Array{<:Any,1}=[missing, missing],
                    center::Array{<:Any,1}=[0., 0., 0.],
                    range_unit::Symbol=:standard,
                    print_filenames::Bool=false,
                    verbose::Bool=verbose_mode)

    printtime("Get particle data: ", verbose)
    checkfortype(dataobject, :particles)

    #Todo: limit to a given lmax
    lmax=dataobject.levelmax # overwrite given lmax
    #checklevelmax(dataobject, lmax)
    isamr = checkuniformgrid(dataobject, lmax)
    #time = dataobject.time

    # create variabe-list and vector-mask (nvarh_corr) for getparticledata-function
    # print selected variables on screen
    nvarp_list, nvarp_i_list, nvarp_corr, read_cpu, used_descriptors = prepvariablelist(dataobject, :particles, vars, lmax, verbose)

    # convert given ranges and print overview on screen
    ranges = prepranges(dataobject, range_unit, verbose, xrange, yrange, zrange, center)

    # read particle-data of the selected variables
    if read_cpu
        if dataobject.descriptor.pversion == 0
            pos_1D, vars_1D, cpus_1D, identity_1D, levels_1D = getparticledata(  dataobject, length(nvarp_list), nvarp_corr, stars, lmax, ranges,
                                         print_filenames, verbose, read_cpu)
        elseif dataobject.descriptor.pversion > 0
            pos_1D, vars_1D, cpus_1D, identity_1D, family_1D, tag_1D, levels_1D = getparticledata(  dataobject, length(nvarp_list), nvarp_corr, stars, lmax, ranges,
                                         print_filenames, verbose, read_cpu)
        end
    else
        if dataobject.descriptor.pversion == 0
            pos_1D, vars_1D, identity_1D, levels_1D = getparticledata(  dataobject, length(nvarp_list), nvarp_corr, stars, lmax, ranges,
                                         print_filenames, verbose, read_cpu)
        elseif dataobject.descriptor.pversion > 0
            pos_1D, vars_1D, identity_1D, family_1D, tag_1D, levels_1D = getparticledata(  dataobject, length(nvarp_list), nvarp_corr, stars, lmax, ranges,
                                         print_filenames, verbose, read_cpu)
        end
    end

    if verbose
        @printf "Found %e particles\n" size(pos_1D)[2]
    end

    # prepare column names for the data table
    names_constr = preptablenames_particles(dataobject, dataobject.nvarp, nvarp_list, used_descriptors, read_cpu, lmax, dataobject.levelmin)

    if lmax != dataobject.levelmin # if AMR
        if dataobject.descriptor.pversion == 0
            Nkeys = [:level, :x, :y, :z, :id]
        elseif dataobject.descriptor.pversion > 0
            Nkeys = [:level, :x, :y, :z, :id, :family, :tag]
        end
    else # if uniform grid
        if dataobject.descriptor.pversion == 0
            Nkeys = [:x, :y, :z, :id]
        elseif dataobject.descriptor.pversion > 0
            Nkeys = [:x, :y, :z, :id, :family, :tag]
        end
    end

    # create data table
    # decouple pos_1D/vars_1D from ElasticArray with ElasticArray.data
    if read_cpu # read also cpu number related to particle
        if isamr
            if dataobject.descriptor.pversion == 0
                data = table( levels_1D[:],
                    pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data, identity_1D[:], cpus_1D[:],
                    [vars_1D[ nvarp_corr[i],: ].data for i in nvarp_i_list]...,
                    names=collect(names_constr), pkey=collect(Nkeys), presorted = false )
            elseif dataobject.descriptor.pversion > 0
                filter!(x->x≠6,nvarp_i_list)
                filter!(x->x≠5,nvarp_i_list)
                data = table( levels_1D[:],
                    pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data, identity_1D[:], family_1D[:], tag_1D[:], cpus_1D[:],
                    [vars_1D[ nvarp_corr[i],: ].data for i in nvarp_i_list]...,
                    names=collect(names_constr), pkey=collect(Nkeys), presorted = false )
            end

        else # if uniform grid
            if dataobject.descriptor.pversion == 0
                data = table(pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data, identity_1D[:], cpus_1D[:],
                    [vars_1D[ nvarp_corr[i],: ].data for i in nvarp_i_list]...,
                    names=collect(names_constr), pkey=collect(Nkeys), presorted = false )
            elseif dataobject.descriptor.pversion > 0
                filter!(x->x≠6,nvarp_i_list)
                filter!(x->x≠5,nvarp_i_list)
                data = table(pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data, identity_1D[:], family_1D[:], tag_1D[:], cpus_1D[:],
                    [vars_1D[ nvarp_corr[i],: ].data for i in nvarp_i_list]...,
                    names=collect(names_constr), pkey=collect(Nkeys), presorted = false )
            end
        end
    else
        if isamr
            if dataobject.descriptor.pversion == 0
                data = table( levels_1D[:],
                    pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data, identity_1D[:],
                    [vars_1D[ nvarp_corr[i],: ].data for i in nvarp_i_list]...,
                    names=collect(names_constr), pkey=collect(Nkeys), presorted = false )
            elseif dataobject.descriptor.pversion > 0
                filter!(x->x≠6,nvarp_i_list)
                filter!(x->x≠5,nvarp_i_list)
                data = table( levels_1D[:],
                    pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data, identity_1D[:], family_1D[:], tag_1D[:],
                    [vars_1D[ nvarp_corr[i],: ].data for i in nvarp_i_list]...,
                    names=collect(names_constr), pkey=collect(Nkeys), presorted = false )
            end
        else # if uniform grid
            if dataobject.descriptor.pversion == 0
                data = table(pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data, identity_1D[:],
                    [vars_1D[ nvarp_corr[i],: ].data for i in nvarp_i_list]...,
                    names=collect(names_constr), pkey=collect(Nkeys), presorted = false )
            elseif dataobject.descriptor.pversion > 0
                filter!(x->x≠6,nvarp_i_list)
                filter!(x->x≠5,nvarp_i_list)
                data = table(pos_1D[1,:].data, pos_1D[2,:].data, pos_1D[3,:].data, identity_1D[:], family_1D[:], tag_1D[:],
                    [vars_1D[ nvarp_corr[i],: ].data for i in nvarp_i_list]...,
                    names=collect(names_constr), pkey=collect(Nkeys), presorted = false )
            end

        end
    end

    printtablememory(data, verbose)

    partdata = PartDataType()
    partdata.data = data
    partdata.info = dataobject
    partdata.lmin = dataobject.levelmin
    partdata.lmax = lmax
    partdata.boxlen = dataobject.boxlen
    partdata.ranges = ranges
    partdata.selected_partvars = names_constr
    partdata.used_descriptors = used_descriptors
    partdata.scale = dataobject.scale
    return partdata

end

function preptablenames_particles(dataobject::InfoType, nvarp::Int, nvarp_list::Array{Int, 1}, used_descriptors::Dict{Any,Any}, read_cpu::Bool, lmax::Real, levelmin::Real)

    if read_cpu
        if lmax != levelmin # if AMR
            if dataobject.descriptor.pversion == 0
                names_constr = [:level, :x, :y, :z, :id, :cpu]
            elseif dataobject.descriptor.pversion > 0
                names_constr = [:level, :x, :y, :z, :id, :family, :tag, :cpu]
            end
        else # if uniform grid
            if dataobject.descriptor.pversion == 0
                names_constr = [:x, :y, :z, :id, :cpu]
            elseif dataobject.descriptor.pversion > 0
                names_constr = [:x, :y, :z, :id, :family, :tag, :cpu]
            end
        end
                    #, Symbol("x"), Symbol("y"), Symbol("z")
    else
        if lmax != levelmin # if AMR
            if dataobject.descriptor.pversion == 0
                names_constr = [:level, :x, :y, :z, :id]
            elseif dataobject.descriptor.pversion > 0
                names_constr = [:level, :x, :y, :z, :id, :family, :tag]
            end
        else # if uniform grid
            if dataobject.descriptor.pversion == 0
                names_constr = [:x, :y, :z, :id]
            elseif dataobject.descriptor.pversion > 0
                names_constr = [:x, :y, :z, :id, :family, :tag]
            end
        end
    end

    for i=1:nvarp
        if in(i, nvarp_list)
            #if  length(used_descriptors) == 0 || !haskey(used_descriptors, i)
            if dataobject.descriptor.pversion == 0
                if i == 1
                    append!(names_constr, [Symbol("vx")] )
                elseif i == 2
                    append!(names_constr, [Symbol("vy")] )
                elseif i == 3
                    append!(names_constr, [Symbol("vz")] )
                elseif i == 4
                    append!(names_constr, [Symbol("mass")] )
                elseif i == 5
                    append!(names_constr, [Symbol("birth")] )
                elseif i > 5
                    append!(names_constr, [Symbol("var$i")] )
                end
            elseif dataobject.descriptor.pversion > 0
                if i == 1
                    append!(names_constr, [Symbol("vx")] )
                elseif i == 2
                    append!(names_constr, [Symbol("vy")] )
                elseif i == 3
                    append!(names_constr, [Symbol("vz")] )
                elseif i == 4
                    append!(names_constr, [Symbol("mass")] )
                elseif i == 7
                    append!(names_constr, [Symbol("birth")] )
                elseif i == 8
                    append!(names_constr, [Symbol("metals")] )
                elseif i > 8
                    append!(names_constr, [Symbol("var$i")] )
                end
            end
            #else append!(names_constr, [used_descriptors[i]] )

            #end
        end
    end

    return names_constr
end


text:
"""
# Fit OLS Linear Regression with fitmodel

   + To fit linear regression with likelihood method, provide the following
   three information in `fitmodel`

   1.  `formula` Provide the equation. Eg. `@formula(y~x1+x2+...)`
   2.  `data` Provide training data as `DataFrame``
   3.  `modelClass`  : LinearRegression()

   ```Julia

   Julia> using RDatasets, StatsModels, StatsPlots

   Julia> df = dataset("datasets", "mtcars");

   Julia> model = @fitmodel(MPG ~ HP + WT+Gear,df,LinearRegression());

   Julia> model.fit

   ────────────────────────────────────────────────────────────────────────────
                     Coef.  Std. Error      t  Pr(>|t|)   Lower 95%   Upper 95%
   ────────────────────────────────────────────────────────────────────────────
   (Intercept)  32.0137     4.63226      6.91    <1e-06  22.5249     41.5024
   HP           -0.0367861  0.00989146  -3.72    0.0009  -0.0570478  -0.0165243
   WT           -3.19781    0.846546    -3.78    0.0008  -4.93188    -1.46374
   Gear          1.01998    0.851408     1.20    0.2410  -0.72405     2.76401
   ────────────────────────────────────────────────────────────────────────────

   ## Fitted model returns following informations

   Julia> model.sigma
   2.5741691724978977

   Julia> model.LogLike
   -73.52638935960971

   Julia> model.AIC
   157.05277871921942

   Julia> model.BIC
   164.38145823321804

   Julia> model.R_sqr
   0.8352309600685555

   Julia> model.Adjusted_R_sqr
   0.8175771343616149

   Julia> model.fittedResponse
   32-element Vector{Float64}:
   23.66884995233871
   22.85340824320635
   25.25355614074087
   20.746171762311327
   17.635570543830113
   ...

   Julia> model.residuals
   32-element Vector{Float64}:
   -2.668849952338711
   -1.8534082432063492
   -2.4535561407408686
   0.6538282376886713
   ...

   Julia> model.Cooks_distance
   32-element Vector{Float64}:
   0.013342034282302845
   0.00688728266731234
   0.015495847517058972
   ...

   Julia> plot(model.Cooks_distance)
   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LinearRegression)
   ans = linear_reg(formula,data);
   ans
end

"""
# Fit Bayesian Linear Regression with Ridge Prior with fitmodel

   # Prior
      v ~ InverseGamma(h,h)
      σ ~ InverseGamma(a0,b0)
      α|σ ~ Normal(0,v*σ)
      β|σ ~ Normal(0,v*σ)

      a0=b0=0.1
      h=0.01 (Default is 0.01, user can choose other values)

   # Data model /likelihood
      y ~ MvNormal(α .+ X * β, σ)

   Julia> model = @fitmodel(MPG ~ HP + WT+Gear,df,LinearRegression(),Prior_Ridge());

   ┌ Info: Found initial step size
   └   ϵ = 0.003125

   Julia> model.summaries
   Summary Statistics
   parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
         Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

            v    6.8184    3.8850     0.0388    0.0680   3316.5165    1.0003      266.6868
            σ    2.6816    0.3879     0.0039    0.0064   3526.9739    1.0003      283.6100
            α   28.4108    5.5547     0.0555    0.1075   2390.7583    1.0009      192.2450
         β[1]   -0.0402    0.0107     0.0001    0.0002   3503.1094    1.0008      281.6910
         β[2]   -2.6507    0.9780     0.0098    0.0189   2486.6245    1.0008      199.9537
         β[3]    1.6426    1.0086     0.0101    0.0189   2522.2853    1.0011      202.8213

   ## All rhat values are close to 1; indicates convergence of Markov Chain.
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LinearRegression,PriorMod::Prior_Ridge,h::Float64=0.01,sim_size::Int64=10000)
   ans = linear_reg(formula,data,Prior_Ridge(),h,sim_size);
   ans
end

"""
# Fit Bayesian Linear Regression with Laplace Prior with fitmodel

   # Prior
      v ~ InverseGamma(h,h)
      σ ~ InverseGamma(a0,b0)
      α|σ ~ Laplace(0,v*σ)
      β|σ ~ Laplace(0,v*σ)

      a0=b0=0.1
      h=0.01 (Default is 0.01, user can choose other values)

   # Data model /likelihood
      y ~ MvNormal(α .+ X * β, σ)

   # h is 0.01 as Default
   # sim_size is 10000

   Julia> model = @fitmodel(MPG ~ HP + WT+Gear,df,LinearRegression(),Prior_Laplace());

   Alternative:

   Julia> model = @fitmodel(MPG ~ HP + WT+Gear,df,LinearRegression(),Prior_Laplace(),0.01,10000);

   ┌ Warning: The current proposal will be rejected due to numerical error(s).
   │   isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)
   └ @ AdvancedHMC ~/.julia/packages/AdvancedHMC/w90s5/src/hamiltonian.jl:47
   ┌ Info: Found initial step size
   └   ϵ = 0.0015625

   Julia> model.summaries
   Summary Statistics
   parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
         Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

            v    4.3206    3.3165     0.0332    0.0567   3115.4478    1.0008      354.1489
            σ    2.6659    0.3791     0.0038    0.0063   4303.1480    0.9999      489.1609
            α   29.1723    5.0562     0.0506    0.0843   2889.0590    1.0000      328.4141
         β[1]   -0.0395    0.0104     0.0001    0.0002   4278.0975    0.9999      486.3132
         β[2]   -2.7379    0.9098     0.0091    0.0156   2979.1114    0.9999      338.6508
         β[3]    1.4906    0.9176     0.0092    0.0141   3140.5300    1.0000      357.0001

   ## All rhat values are close to 1; indicates convergence of Markov Chain.
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LinearRegression,PriorMod::Prior_Laplace,h::Float64=0.01,sim_size::Int64=10000)
   ans = linear_reg(formula,data,Prior_Laplace(),h,sim_size);
   ans
end

"""
# Fit Bayesian Linear Regression with Cauchy Prior with fitmodel

   #priors
      σ ~ Truncated(TDist(1),0,Inf)
      α|σ ~ σ*TDist(1)
      β|σ ~ σ*TDist(1)

   1. t-distribution with 1 df follows Cauchy distribution
   2. On σ we elicite Half-Cauchy distribution (see Gelman 2006)

   # Data model /likelihood
      y ~ MvNormal(α .+ X * β, σ)

   # sim_size is 10000 (default). In the following example we considered 20000.

   Julia> model = @fitmodel(MPG ~ HP + WT+Gear,df,LinearRegression(),Prior_Cauchy(),20000);

   ┌ Info: Found initial step size
   └   ϵ = 0.00078125

   Julia> model.summaries
   Summary Statistics
   parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
         Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

            σ    2.5919    0.3439     0.0024    0.0034   9909.9764    1.0000      771.7449
            α   30.2828    4.6736     0.0330    0.0677   4769.9124    1.0001      371.4596
         β[1]   -0.0395    0.0101     0.0001    0.0001   7610.6711    1.0000      592.6852
         β[2]   -2.8297    0.8653     0.0061    0.0119   5136.0610    1.0000      399.9736
         β[3]    1.2719    0.8446     0.0060    0.0117   5086.2652    1.0001      396.0957

   # All rhat values are close to 1; indicates convergence of Markov Chain.
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LinearRegression,PriorMod::Prior_Cauchy,sim_size::Int64=10000)
   ans = linear_reg(formula,data,Prior_Cauchy(),sim_size);
   ans
end

"""
# Fit Bayesian Linear Regression with t(ν) distributed Prior with fitmodel

   #priors
      ν ~ InverseGamma(h,h)
      σ ~ InverseGamma(a0,b0)
      α|σ ~ σ*TDist(ν)
      β|σ ~ σ*TDist(ν)

   1. t-distribution with ν df follows Cauchy distribution
   2. On σ we elicite InverseGamma(a0,b0) distribution (see Gelman 2006)

   # Data model /likelihood
      y ~ MvNormal(α .+ X * β, σ)

   # sim_size is 10000 (default).

   Julia> model = @fitmodel(MPG ~ HP + WT+Gear,df,LinearRegression(),Prior_TDist());

   ┌ Info: Found initial step size
   └   ϵ = 0.0001953125

   Julia> model.summaries
   Summary Statistics
   parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
         Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

            ν    1.0513    0.5599     0.0056    0.0082   5053.9417    1.0000      366.6794
            σ    2.6200    0.3526     0.0035    0.0050   5835.7064    0.9999      423.3989
            α   30.2533    4.7311     0.0473    0.0918   2619.3110    1.0006      190.0393
         β[1]   -0.0394    0.0100     0.0001    0.0001   4267.6480    1.0002      309.6313
         β[2]   -2.8329    0.8630     0.0086    0.0155   2854.0611    0.9999      207.0711
         β[3]    1.2791    0.8614     0.0086    0.0165   2659.5009    1.0005      192.9552

   # All rhat values are close to 1; indicates convergence of Markov Chain.
   # We estimate the posterior estimates of degrees of freedom is 1.0513.
   So posterior distribution of β is close to Cauchy.
   Hence it is better to estimate the posterior median. Instead of posterior mean.

   Julia> model.quantiles
   Quantiles
   parameters      2.5%     25.0%     50.0%     75.0%     97.5%
         Symbol   Float64   Float64   Float64   Float64   Float64

            ν    0.3849    0.6656    0.9170    1.2905    2.4945
            σ    2.0373    2.3724    2.5847    2.8248    3.4096
            α   20.5143   27.2162   30.2498   33.3969   39.3762
         β[1]   -0.0591   -0.0460   -0.0394   -0.0328   -0.0199
         β[2]   -4.5112   -3.4040   -2.8377   -2.2691   -1.0928
         β[3]   -0.3783    0.7007    1.2674    1.8364    3.0627

   ## Check if trace-plots are stationary !
   plot(model.chain)

"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LinearRegression,PriorMod::Prior_TDist,h::Float64=2.0,sim_size::Int64=10000)
   ans = linear_reg(formula,data,Prior_TDist(),h,sim_size);
   ans
end

"""
# Fit Bayesian Linear Regression with Uniform Prior
# distributed Prior with fitmodel

   #priors
   v=1/h;
   σ ~ Uniform(0,v)
   α ~ Uniform(-v*σ,v*σ)
   β ~ filldist(Uniform(-v,v), predictors)

   # Data model /likelihood
      y ~ MvNormal(α .+ X * β, σ)

   # sim_size is 10000 (default).

   Julia> model = @fitmodel(MPG ~ HP + WT+Gear,df
                           ,LinearRegression()
                           ,Prior_Uniform());

   ## Check if trace-plots are stationary !
   Julia> plot(model.chain)

"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LinearRegression,PriorMod::Prior_Uniform,h::Float64=0.01,sim_size::Int64=10000)
   ans = linear_reg(formula,data,Prior_TDist(),h,sim_size);
   ans
end

"""
# Logistic Regression with Logit link : GLM with Scoring Method

    ## Link Function
    z = X*β

    prob = Logit.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end
   ```
   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                  ,turnout,LogisticRegression(),Logit());

   Julia> model.fit

   ────────────────────────────────────────────────────────────────────────────
                     Coef.  Std. Error      z  Pr(>|z|)   Lower 95%   Upper 95%
   ────────────────────────────────────────────────────────────────────────────
   (Intercept)  -3.03426    0.325927    -9.31    <1e-19  -3.67307    -2.39546
   Age           0.0283543  0.00346034   8.19    <1e-15   0.0215722   0.0351365
   Race: white   0.250798   0.146457     1.71    0.0868  -0.0362521   0.537847
   Income        0.177112   0.0271516    6.52    <1e-10   0.123896    0.230328
   Educate       0.175634   0.0203308    8.64    <1e-17   0.135786    0.215481
   ────────────────────────────────────────────────────────────────────────────

   Julia> model.LogLike
   -1011.9906318515576

   Julia> model.AIC
   2033.9812637031152

   Julia> model.BIC
   2061.985776000826
   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LogisticRegression,Link::Logit)
   ans = logistic_reg(formula,data,"LogitLink");
   ans
end

"""
# Logistic Regression with Probit link : GLM with Scoring Method

    ## Link Function
    z = X*β

    prob = Probit.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                  ,turnout,LogisticRegression(),Probit());
   Julia> model.fit

   ────────────────────────────────────────────────────────────────────────────
                     Coef.  Std. Error      z  Pr(>|z|)   Lower 95%   Upper 95%
   ────────────────────────────────────────────────────────────────────────────
   (Intercept)  -1.76141    0.188556    -9.34    <1e-20  -2.13097    -1.39185
   Age           0.0164973  0.00199897   8.25    <1e-15   0.0125794   0.0204152
   Race: white   0.162856   0.0876885    1.86    0.0633  -0.0090108   0.334722
   Income        0.0963117  0.0149675    6.43    <1e-09   0.066976    0.125647
   Educate       0.10417    0.0116713    8.93    <1e-18   0.0812949   0.127046
   ────────────────────────────────────────────────────────────────────────────
   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LogisticRegression,Link::Probit)
   ans = logistic_reg(formula,data,"ProbitLink");
   ans
end

"""
# Logistic Regression with CloglogLink link : GLM with Scoring Method

    ## Link Function
    z = X*β

    prob = CloglogLink.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia
   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                  ,turnout,LogisticRegression(),Cloglog());
   Julia> model.fit

   ─────────────────────────────────────────────────────────────────────────────
                     Coef.  Std. Error       z  Pr(>|z|)   Lower 95%   Upper 95%
   ─────────────────────────────────────────────────────────────────────────────
   (Intercept)  -1.94617    0.184123    -10.57    <1e-25  -2.30704    -1.58529
   Age           0.0147857  0.00184088    8.03    <1e-15   0.0111776   0.0183937
   Race: white   0.185139   0.087101      2.13    0.0335   0.014424    0.355854
   Income        0.0768268  0.0126411     6.08    <1e-08   0.0520506   0.101603
   Educate       0.0983976  0.0108857     9.04    <1e-18   0.077062    0.119733
   ─────────────────────────────────────────────────────────────────────────────
   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LogisticRegression,Link::Cloglog)
   ans = logistic_reg(formula,data,"CloglogLink");
   ans
end

"""
# Logistic Regression with Cauchit link : GLM with Scoring Method

    ## Link Function
    z = X*β

    prob = Cauchit.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

    ```Julia
      Julia> turnout = dataset("Zelig", "turnout")

      Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                     ,turnout,LogisticRegression(),Cauchit());
      Julia> model.fit

      ────────────────────────────────────────────────────────────────────────────
                        Coef.  Std. Error      z  Pr(>|z|)   Lower 95%   Upper 95%
      ────────────────────────────────────────────────────────────────────────────
      (Intercept)  -3.16889    0.384429    -8.24    <1e-15  -3.92235    -2.41542
      Age           0.0304105  0.00413473   7.35    <1e-12   0.0223066   0.0385144
      Race: white   0.181839   0.144766     1.26    0.2091  -0.101898    0.465576
      Income        0.235267   0.038152     6.17    <1e-09   0.16049     0.310043
      Educate       0.169276   0.0240098    7.05    <1e-11   0.122217    0.216334
      ────────────────────────────────────────────────────────────────────────────
   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LogisticRegression,Link::Cauchit)
   ans = logistic_reg(formula,data,"CauchitLink");
   ans
end

"""
# Logistic Regression with Ridge prior using HMC method in Turing

    ## priors
    λ ~ InverseGamma(h,h)
    β ~ Normal(0,λ)

    ## Link Function
    z = X*β

    prob = Logit.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                  ,turnout,LogisticRegression(),Cauchit());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LogisticRegression,Link::Logit,PriorMod::Prior_Ridge, h::Real=0.1, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Logit(),Prior_Ridge(),h,sim_size)
   ans
end

"""
# Logistic Regression with Ridge prior using HMC method in Turing

    ## priors
    λ ~ InverseGamma(h,h)
    β ~ Normal(0,λ)

    ## Link Function
    z = X*β

    prob = Probit.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                  ,turnout,LogisticRegression()
                  ,Probit(),Prior_Ridge());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LogisticRegression,Link::Probit,PriorMod::Prior_Ridge, h::Real=0.1, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Probit(),Prior_Ridge(),h,sim_size)
   ans
end

"""
# Logistic Regression with Ridge prior using HMC method in Turing

    ## priors
    λ ~ InverseGamma(h,h)
    β ~ Normal(0,λ)

    ## Link Function
    z = X*β

    prob = Cloglog.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Cloglog()
                           ,Prior_Ridge());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LogisticRegression,Link::Cloglog,PriorMod::Prior_Ridge, h::Real=0.1, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Cloglog(),Prior_Ridge(),h,sim_size)
   ans
end

"""
# Logistic Regression with Ridge prior using HMC method in Turing

    ## priors
    λ ~ InverseGamma(h,h)
    β ~ Normal(0,λ)

    ## Link Function
    z = X*β

    prob = Cauchit.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Cauchit()
                           ,Prior_Ridge());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame,modelClass::LogisticRegression,Link::Cauchit,PriorMod::Prior_Ridge, h::Real=0.1, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Cauchit(),Prior_Ridge(),h,sim_size)
   ans
end

"""
# Logistic Regression with Laplace prior using HMC method in Turing

    ## priors
    λ ~ InverseGamma(h,h)
    β ~ Laplace(0,λ)

    ## Link Function
    z = X*β

    prob = Logit.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Logit()
                           ,Prior_Laplace());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Logit, PriorMod::Prior_Laplace, h::Real=0.1, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Logit(),Prior_Laplace(),h,sim_size)
   ans
end

"""
# Logistic Regression with Laplace prior using HMC method in Turing

    ## priors
    λ ~ InverseGamma(h,h)
    β ~ Laplace(0,λ)

    ## Link Function
    z = X*β

    prob = Probit.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Probit()
                           ,Prior_Laplace());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Probit, PriorMod::Prior_Laplace, h::Real=0.1, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Probit(),Prior_Laplace(),h,sim_size)
   ans
end

"""
# Logistic Regression with Laplace prior using HMC method in Turing

    ## priors
    λ ~ InverseGamma(h,h)
    β ~ Laplace(0,λ)

    ## Link Function
    z = X*β

    prob = Cloglog.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Cloglog()
                           ,Prior_Laplace());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Cloglog, PriorMod::Prior_Laplace, h::Real=0.1, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Cloglog(),Prior_Laplace(),h,sim_size)
   ans
end

"""
# Logistic Regression with Laplace prior using HMC method in Turing

    ## priors
    λ ~ InverseGamma(h,h)
    β ~ Laplace(0,λ)

    ## Link Function
    z = X*β

    prob = Cauchit.(z)

    #likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Cauchit()
                           ,Prior_Laplace());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Cauchit, PriorMod::Prior_Laplace, h::Real=0.1, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Cauchit(),Prior_Laplace(),h,sim_size)
   ans
end

"""
# Logistic Regression with Cauchy prior using HMC method in Turing

   ## priors

    λ   ~ Truncated(TDist(1),0,Inf) # Half-Cauchy prior
    β|λ ~ λ*TDist(1)  # Cauchy prior

    ## Link Function
    z = X*β

    prob = Logit.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Logit()
                           ,Prior_Cauchy());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Logit, PriorMod::Prior_Cauchy, h::Real=1.0, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Logit(),Prior_Cauchy(),h,sim_size)
   ans
end

"""
# Logistic Regression with Cauchy prior using HMC method in Turing

   ## priors

    λ   ~ Truncated(TDist(1),0,Inf) # Half-Cauchy prior
    β|λ ~ λ*TDist(1)  # Cauchy prior

    ## Link Function
    z = X*β

    prob = Probit.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Probit()
                           ,Prior_Cauchy());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Probit, PriorMod::Prior_Cauchy, h::Real=1.0, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Probit(),Prior_Cauchy(),h,sim_size)
   ans
end

"""
# Logistic Regression with Cauchy prior using HMC method in Turing

   ## priors

    λ   ~ Truncated(TDist(1),0,Inf) # Half-Cauchy prior
    β|λ ~ λ*TDist(1)  # Cauchy prior

    ## Link Function
    z = X*β

    prob = Cloglog.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Cloglog()
                           ,Prior_Cauchy());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Cloglog, PriorMod::Prior_Cauchy, h::Real=1.0, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Cloglog(),Prior_Cauchy(),h,sim_size)
   ans
end

"""
# Logistic Regression with Cauchy prior using HMC method in Turing

   ## priors

    λ   ~ Truncated(TDist(1),0,Inf) # Half-Cauchy prior
    β|λ ~ λ*TDist(1)  # Cauchy prior

    ## Link Function
    z = X*β

    prob = Cauchit.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Cauchit()
                           ,Prior_Cauchy());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Cauchit, PriorMod::Prior_Cauchy, h::Real=1.0, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Cauchit(),Prior_Cauchy(),h,sim_size)
   ans
end

"""
# Logistic Regression with T-Dist prior using HMC method in Turing

   ## priors

    λ ~ InverseGamma(h,h)
    ν ~ InverseGamma(h,h)
    β|λ,ν ~ λ*TDist(ν)

    ## Link Function
    z = X*β

    prob = Logit.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Logit()
                           ,Prior_TDist());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Logit, PriorMod::Prior_TDist, h::Real=1.0, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Logit(),Prior_TDist(),h,sim_size)
   ans
end

"""
# Logistic Regression with T-Dist prior using HMC method in Turing

   ## priors

    λ ~ InverseGamma(h,h)
    ν ~ InverseGamma(h,h)
    β|λ,ν ~ λ*TDist(ν)

    ## Link Function
    z = X*β

    prob = Probit.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Probit()
                           ,Prior_TDist());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Probit, PriorMod::Prior_TDist, h::Real=1.0, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Probit(),Prior_TDist(),h,sim_size)
   ans
end

"""
# Logistic Regression with T-Dist prior using HMC method in Turing

   ## priors

    λ ~ InverseGamma(h,h)
    ν ~ InverseGamma(h,h)
    β|λ,ν ~ λ*TDist(ν)

    ## Link Function
    z = X*β

    prob = Cloglog.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Cloglog()
                           ,Prior_TDist());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Cloglog, PriorMod::Prior_TDist, h::Real=1.0, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Cloglog(),Prior_TDist(),h,sim_size)
   ans
end

"""
# Logistic Regression with T-Dist prior using HMC method in Turing

   ## priors

    λ ~ InverseGamma(h,h)
    ν ~ InverseGamma(h,h)
    β|λ,ν ~ λ*TDist(ν)

    ## Link Function
    z = X*β

    prob = Cauchit.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Cauchit()
                           ,Prior_TDist());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Cauchit, PriorMod::Prior_TDist, h::Real=1.0, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Cauchit(),Prior_TDist(),h,sim_size)
   ans
end

"""
# Logistic Regression with Uniform prior using HMC method in Turing

   ## priors
    v ~ InverseGamma(h,h)
    β ~ filldist(Uniform(-v,v), p)

    ## Link Function
    z = X*β

    prob = Cauchit.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Logit()
                           ,Prior_Uniform());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Logit, PriorMod::Prior_Uniform, h::Real=0.01, level::Real=0.95, sim_size::Int64=10000)
  ans = logistic_reg(formula,data,Logit(),Prior_Uniform(),h,sim_size)
  ans
end

"""
# Logistic Regression with Uniform prior using HMC method in Turing

   ## priors
    v ~ InverseGamma(h,h)
    β ~ filldist(Uniform(-v,v), p)

    ## Link Function
    z = X*β

    prob = Probit.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Probit()
                           ,Prior_Uniform());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Probit, PriorMod::Prior_Uniform, h::Real=0.01, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Probit(),Prior_Uniform(),h,sim_size)
   ans
end

"""
# Logistic Regression with Uniform prior using HMC method in Turing

   ## priors
    v ~ InverseGamma(h,h)
    β ~ filldist(Uniform(-v,v), p)

    ## Link Function
    z = X*β

    prob = Cloglog.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Cloglog()
                           ,Prior_Uniform());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Cloglog, PriorMod::Prior_Uniform, h::Real=0.01, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Cloglog(),Prior_Uniform(),h,sim_size)
   ans
end

"""
# Logistic Regression with Uniform prior using HMC method in Turing

   ## priors
    v ~ InverseGamma(h,h)
    β ~ filldist(Uniform(-v,v), p)

    ## Link Function
    z = X*β

    prob = Cauchit.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Bernoulli(prob[i])
    end

   ```Julia

   Julia> turnout = dataset("Zelig", "turnout")

   Julia> model = @fitmodel(Vote ~ Age + Race +Income + Educate
                           ,turnout
                           ,LogisticRegression()
                           ,Cauchit()
                           ,Prior_Uniform());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::LogisticRegression, Link::Cauchit, PriorMod::Prior_Uniform, h::Real=0.01, level::Real=0.95, sim_size::Int64=10000)
   ans = logistic_reg(formula,data,Cauchit(),Prior_Uniform(),h,sim_size)
   ans
end

"""
# Poisson Regression - Likelihood Method

    ## Link Function
    z = X*β

    μ = exp.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Poisson(μ[i])
    end

   ```Julia

   Julia> sanction = dataset("Zelig", "sanction");

   Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                           , sanction
                           , PoissonRegression());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::PoissonRegression)
   ans = Poisson_Reg(formula,data)
   ans
end

"""
# Poisson Regression with Ridge Prior using HMC with Turing

    ## priors
    λ ~ InverseGamma(h,h)
    α|λ ~ λ*Normal(0,λ)
    β|λ ~ λ*Normal(0,λ)

    ## Link Function
    z = X*β

    μ = exp.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Poisson(μ[i])
    end

   ```Julia

   Julia> sanction = dataset("Zelig", "sanction");

   Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                           , sanction
                           , PoissonRegression()
                           , Prior_Ridge());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::PoissonRegression,PriorMod::Prior_Ridge,h::Float64=0.1,sim_size::Int64=10000)
   ans = Poisson_Reg(formula,data,Prior_Ridge(),h,sim_size)
   ans
end

"""
# Poisson Regression with Laplace Prior using HMC with Turing

    ## priors
    λ ~ InverseGamma(h,h)
    α|λ ~ λ * Laplace(0,λ)
    β|λ ~ λ * Laplace(0,λ)

    ## Link Function
    z = X*β

    μ = exp.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Poisson(μ[i])
    end

   ```Julia

   Julia> sanction = dataset("Zelig", "sanction");

   Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                           , sanction
                           , PoissonRegression()
                           , Prior_Laplace());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::PoissonRegression,PriorMod::Prior_Laplace,h::Float64=0.1,sim_size::Int64=10000)
   ans = Poisson_Reg(formula,data,Prior_Laplace(),h,sim_size)
   ans
end

"""
# Poisson Regression with Cauchy Prior using HMC with Turing

    ## priors
    λ~InverseGamma(h,h)
    α ~ TDist(1)*λ
    β ~ filldist(TDist(1)*λ, p)

    ## Link Function
    z = X*β

    μ = exp.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Poisson(μ[i])
    end

   ```Julia

   Julia> sanction = dataset("Zelig", "sanction");

   Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                           , sanction
                           , PoissonRegression()
                           , Prior_Cauchy());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::PoissonRegression,PriorMod::Prior_Cauchy,h::Float64=1.0,sim_size::Int64=10000)
   ans = Poisson_Reg(formula,data,Prior_Cauchy(),h,sim_size)
   ans
end

"""
# Poisson Regression with TDist Prior using HMC with Turing

    ## priors
    λ~InverseGamma(h,h)
    α|ν ~ TDist(ν)*λ
    β|ν ~ TDist(ν)*λ

    ## Link Function
    z = X*β

    μ = exp.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Poisson(μ[i])
    end

   ```Julia

   Julia> sanction = dataset("Zelig", "sanction");

   Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                           , sanction
                           , PoissonRegression()
                           , Prior_TDist());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::PoissonRegression,PriorMod::Prior_TDist,h::Float64=2.0,sim_size::Int64=10000)
   ans = Poisson_Reg(formula,data,Prior_TDist(),h,sim_size)
   ans
end

"""
# Poisson Regression with Uniform Prior using HMC with Turing

    ## priors
    λ~InverseGamma(h,h)
    α ~ Uniform(-λ,λ)
    β ~ Uniform(-λ,λ)

    ## Link Function
    z = X*β

    μ = exp.(z)

    ## likelihood
    for i = 1:n
        y[i] ~ Poisson(μ[i])
    end

   ```Julia

   Julia> sanction = dataset("Zelig", "sanction");

   Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                           , sanction
                           , PoissonRegression()
                           , Prior_Uniform());

   ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::PoissonRegression,PriorMod::Prior_Uniform,h::Float64=1.0,sim_size::Int64=10000)
   ans = Poisson_Reg(formula,data,Prior_Uniform(),h,sim_size)
   ans
end

"""
## NegBinomial Regression with Ridge Prior

    #priors
    λ ~ InverseGamma(h,h)
    α|λ ~ Normal(0,λ)
    β|λ ~ Normal(0,λ)

    ## link
    z = α .+ X * β
    mu = exp.(z)

    #likelihood
    for i = 1:n
      y[i] ~ NegativeBinomial2(mu[i],λ)
    end

    h::Float64 = 0.1 (default)
    sim_size::Int64 = 10000 (default)
    ```Julia
    Julia> using RDatasets
    Julia> sanction = dataset("Zelig", "sanction");
    Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                            , sanction
                            , NegBinomRegression()
                            , Prior_Ridge());

    ┌ Info: Found initial step size
    └   ϵ = 0.025

    Julia> model.summaries
    Summary Statistics
    parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

            λ    2.0350    0.4392     0.0044    0.0052   7933.2864    1.0000      639.1111
            α   -1.0802    0.5194     0.0052    0.0082   3565.0359    1.0007      287.2018
            β[1]   -0.0053    0.1637     0.0016    0.0022   5160.0061    1.0003      415.6937
            β[2]    1.0617    0.1320     0.0013    0.0017   5259.8195    1.0001      423.7348
            β[3]   -0.1775    0.5511     0.0055    0.0057   7967.5892    0.9999      641.8746
            β[4]    1.2770    0.3160     0.0032    0.0038   7812.4049    0.9999      629.3728
            β[5]    0.1518    0.2839     0.0028    0.0036   6062.4239    1.0006      488.3931

    ```
"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::NegBinomRegression,PriorMod::Prior_Ridge,h::Float64=0.1,sim_size::Int64=10000)
    ans = NegBinom_Reg(formula,data,Prior_Ridge(),h,sim_size)
    ans
 end

 """
 ## NegBinomial Regression with Laplace Prior

     #priors
     λ ~ InverseGamma(h,h)
     α|λ ~ Laplace(0,λ)
     β|λ ~ Laplace(0,λ)

     ## link
     z = α .+ X * β
     mu = exp.(z)

     #likelihood
     for i = 1:n
       y[i] ~ NegativeBinomial2(mu[i],λ)
     end

     h::Float64 = 0.1 (default)
     sim_size::Int64 = 10000 (default)

    ```Julia

    Julia> using RDatasets
    Julia> sanction = dataset("Zelig", "sanction");
    Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                            , sanction
                            , NegBinomRegression()
                            , Prior_Laplace());

    ```

 """
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::NegBinomRegression,PriorMod::Prior_Laplace,h::Float64=0.1,sim_size::Int64=10000)
    ans = NegBinom_Reg(formula,data,Prior_Laplace(),h,sim_size)
    ans
end

"""
 ## NegBinomial Regression with Cauchy Prior

     #priors
     λ~InverseGamma(h,h)
     α|λ ~ λ*TDist(1)
     β|λ ~ λ*TDist(1)

     ## link
     z = α .+ X * β
     mu = exp.(z)

     #likelihood
     for i = 1:n
       y[i] ~ NegativeBinomial2(mu[i],λ)
     end

     h::Float64 = 1.0 (default)
     sim_size::Int64 = 10000 (default)

    ```Julia

    Julia> using RDatasets
    Julia> sanction = dataset("Zelig", "sanction");
    Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                            , sanction
                            , NegBinomRegression()
                            , Prior_Cauchy());

    ```

"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::NegBinomRegression,PriorMod::Prior_Cauchy,h::Float64=1.0,sim_size::Int64=10000)
    ans = NegBinom_Reg(formula,data,Prior_Cauchy(),h,sim_size)
    ans
end

"""
 ## NegBinomial Regression with TDist Prior

     #priors
     λ ~ InverseGamma(h,h)
     ν ~ InverseGamma(h,h)
     α|λ,μ ~ λ*TDist(ν)
     β|λ,μ ~ λ*TDist(ν)

     ## link
     z = α .+ X * β
     mu = exp.(z)

     #likelihood
     for i = 1:n
       y[i] ~ NegativeBinomial2(mu[i],λ)
     end

     h::Float64 = 1.0 (default)
     sim_size::Int64 = 10000 (default)

    ```Julia

    Julia> using RDatasets
    Julia> sanction = dataset("Zelig", "sanction");
    Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                            , sanction
                            , NegBinomRegression()
                            , Prior_TDist());

    ```

"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::NegBinomRegression,PriorMod::Prior_TDist,h::Float64=1.0,sim_size::Int64=10000)
    ans = NegBinom_Reg(formula,data,Prior_TDist(),h,sim_size)
    ans
end

"""
## NegBinomial Regression with Uniform Prior

    Ibrahim and Laud (JASA, 1990) showed that the uniform flat priors
    for GLM's can lead to improper posterior distributions thus making
    them undesirable. In such cases, the Markov Chain struggles to converge.
    Even if it converges, results are unreliable.

     #priors
     λ~InverseGamma(h,h)
     α|λ ~ Uniform(-λ,λ)
     β|λ ~ Uniform(-λ,λ)

     ## link
     z = α .+ X * β
     mu = exp.(z)

     #likelihood
     for i = 1:n
       y[i] ~ NegativeBinomial2(mu[i],λ)
     end

     h::Float64 = 1.0 (default)
     sim_size::Int64 = 10000 (default)

    ```Julia

    Julia> using RDatasets
    Julia> sanction = dataset("Zelig", "sanction");
    Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                            , sanction
                            , NegBinomRegression()
                            , Prior_Uniform());

    ```

"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::NegBinomRegression,PriorMod::Prior_Uniform,h::Float64=0.1,sim_size::Int64=10000)
    ans = NegBinom_Reg(formula,data,Prior_Uniform(),h,sim_size)
    ans
end

"""
## NegBinomial Regression - Likelihood method

     ## link
     z = α .+ X * β
     mu = exp.(z)

     #likelihood
     for i = 1:n
       y[i] ~ NegativeBinomial2(mu[i],λ)
     end

    ```Julia

    Julia> using RDatasets
    Julia> sanction = dataset("Zelig", "sanction");
    Julia> model = @fitmodel(Num ~ Target + Coop + NCost
                            , sanction
                            , NegBinomRegression());

    ```

"""
function fitmodel(formula::FormulaTerm, data::DataFrame, modelClass::NegBinomRegression)
    ans = NegBinom_Reg(formula,data,"LogLink")
    ans
end

"""
```julia
   @fitmodel(formula, args...)
```

Macro for calling `fitmodel` without using `@formula` to create the formula.

# Example
```julia
using CRRao, RDatasets
sanction = dataset("Zelig", "sanction")
model = @fitmodel(Num ~ Target + Coop + NCost, sanction, NegBinomRegression())
```
"""
macro fitmodel(formula, args...)
   quote
      fitmodel(@formula($formula), $(map(esc, args)...))
   end
end


text:
import BenchmarkConfigSweeps
import BenchmarkTools
using DataFrames
using FileIO
using JSON
using Statistics
using VegaLite

results = only(BenchmarkTools.load(joinpath(@__DIR__, "build/results.json")))
df_raw = DataFrame(BenchmarkConfigSweeps.flattable(results))

access_param = "Access density"

begin
    df = select(df_raw, Not(:trial))
    df = select(df, Not(r"JULIA_.*"))
    df[:, :ms] = map(t -> mean(t).time, df_raw.trial) ./ 1e6
    df[:, :Implementation] = map(df.alg) do alg
        if alg === :base_seq || alg == :base_par
            Symbol("Base.Dict + Divide-and-Conquer")
        elseif alg === :cdict_seq || alg === :cdict_par
            :ConcurrentDict
        else
            error("unknown alg = ", alg)
        end
    end
    datasize = 2^19
    df[:, access_param] = datasize ./ df.nkeys
    df
end
#-

df_speedup = combine(groupby(df, Not([:ms, :ntasks, :alg, :Implementation]))) do g
    baseline = only(g.ms[g.alg.===:base_seq])
    hcat(g, DataFrame((; speedup = baseline ./ g.ms)))
end
#-

function parallel_algorithms(df)
    idx = df.alg .∈ Ref((:base_par, :cdict_par))
    return df[idx, :]
end

plt = @vlplot(
    facet = {column = {field = :Implementation}},
    spec = {
        layer = [
            {
                # :line,
                mark = {:line, point = true},
                encoding = {
                    x = {:ntasks, type = :quantitative, title = "Number of Tasks"},
                    y = {
                        :speedup,
                        type = :quantitative,
                        title = "Speedup wrt sequential program",
                    },
                    color = {field = access_param, type = :ordinal},
                },
            },
            {mark = :rule, encoding = {y = {datum = 1}}},
        ],
    },
    data = parallel_algorithms(df_speedup),
)

save(joinpath(@__DIR__, "build/results.png"), plt)
save(joinpath(@__DIR__, "build/results.svg"), plt)

plt


text:
module MbedTLSStream

using MbedTLS

const roots = ObjectIdDict()
push_root(obj) = (roots[obj] = nothing)
pop_root(obj) = delete!(roots, obj)

function set_error(cls::LibGit2.Error.Class, msg::String)
    ccall((:giterr_set_str, :libgit2), Void, (Cint, Cstring), Cint(cls), msg)
    return
end

function ssl_set_error(error)
    @assert error != MbedTLS.MBEDTLS_ERR_SSL_WANT_READ
    @assert error != MbedTLS.MBEDTLS_ERR_SSL_WANT_WRITE

    errstr = error != 0 ? MbedTLS.strerror(error) : "unknown SSL error"

    set_error(LibGit2.Error.SSL, "SSL error: $(num2hex(Int16(error))) - $errstr")
    return Cint(error)
end

function memcpy(src::Ptr{UInt8}, sz)
    dest = convert(Ptr{UInt8}, Libc.malloc(sizeof(UInt8) * sz))
    unsafe_copy!(dest, src, sz)
    return dest
end
memcpy(src::Cstring, sz) = memcpy(convert(Ptr{UInt8}, src), sz)

function init_ssl()
    info("Initializing MbedTLS...")

    global tls_entropy = MbedTLS.Entropy()
    global tls_rng = MbedTLS.CtrDrbg()
    MbedTLS.seed!(tls_rng, tls_entropy)

    global tls_conf = MbedTLS.SSLConfig()
    MbedTLS.rng!(tls_conf, tls_rng)
    MbedTLS.config_defaults!(tls_conf)
    MbedTLS.authmode!(tls_conf, MbedTLS.MBEDTLS_SSL_VERIFY_REQUIRED)
    MbedTLS.ca_chain!(tls_conf, MbedTLS.crt_parse_file(ENV["SSL_CERT_FILE"]))
end

@enum GIT_CERT CERT_NONE CERT_X509 GIT_CERT_HOSTKEY_LIBSSH2 GIT_CERT_STRARRAY

immutable GitCert
    cert_type::Cint
end
GitCert(ctype::GIT_CERT) = GitCert(Cint(ctype))

immutable GitCertX509
    parent::GitCert
    data::Ptr{Void} # Pointer to the X.509 certificate data
    len::Csize_t    # Length of the memory block pointed to by `data`.
end
GitCertX509(data::Ptr{Void}, len::Int) = GitCertX509(GitCert(CERT_X509), data, Csize_t(len))
GitCertX509() = GitCertX509(C_NULL, 0)

immutable GitStream
    version::Cint

    encrypted::Cint
    proxy_support::Cint

    connect::Ptr{Void}     # int (*connect)(struct git_stream *);
    certificate::Ptr{Void} # int (*certificate)(git_cert **, struct git_stream *);
    set_proxy::Ptr{Void}   # int (*set_proxy)(struct git_stream *, const char *proxy_url);
    read::Ptr{Void}        # ssize_t (*read)(struct git_stream *, void *, size_t);
    write::Ptr{Void}       # ssize_t (*write)(struct git_stream *, const char *, size_t, int);
    close::Ptr{Void}       # int (*close)(struct git_stream *);
    free::Ptr{Void}        # void (*free)(struct git_stream *);
end
GitStream(; encrypted::Cint = zero(Cint),
            proxy_support::Cint = zero(Cint),
            connect::Ptr{Void} = C_NULL,
            certificate::Ptr{Void} = C_NULL,
            set_proxy::Ptr{Void} = C_NULL,
            read::Ptr{Void} = C_NULL,
            write::Ptr{Void} = C_NULL,
            close::Ptr{Void} = C_NULL,
            free::Ptr{Void} = C_NULL) =
    GitStream(one(Cint),
              encrypted,
              proxy_support,
              connect,
              certificate,
              set_proxy,
              read,
              write,
              close,
              free)

type TLSStream
    sock::TCPSocket
    host::String
    port::String
    ctx::MbedTLS.SSLContext # mbedtls_ssl_context *ssl;
    cert::GitCertX509 # git_cert_x509 cert_info;
end

type TLSWrapper
    parent::GitStream
    stream::Ptr{TLSStream}
end

function getstream(tls_ptr::Ptr{TLSWrapper})
    tls = unsafe_load(tls_ptr,1)
    stream = unsafe_pointer_to_objref(tls.stream)
end

function bio_write(ctx_stm::Ptr{Void}, c_msg::Ptr{UInt8}, sz::Csize_t)
    print("bio_write: $sz => ")
    jl_ctx = unsafe_pointer_to_objref(ctx_stm)
    print("B: $jl_ctx ")
    n = unsafe_write(jl_ctx, c_msg, sz)
    println("A: $jl_ctx => $n")
    return Cint(n)
end

function bio_read(ctx_stm::Ptr{Void}, c_msg::Ptr{UInt8}, sz::Csize_t)
    print("bio_read: $sz <= ")
    jl_ctx = unsafe_pointer_to_objref(ctx_stm)
    print("B: $jl_ctx")
    n = unsafe_read(jl_ctx, c_msg, sz)
    n = n != nothing ? n : sz
    println("A: $jl_ctx => $n")
    return Cint(n)
end

function tls_stream_connect(tls_ptr::Ptr{TLSWrapper})
    stream = getstream(tls_ptr)
    try
        # println(stream.host, " ", stream.port)

        Base.connect!(stream.sock, stream.host, parse(UInt16, stream.port))

        MbedTLS.hostname!(stream.ctx, stream.host)
        MbedTLS.set_bio!(stream.ctx, stream.sock)
        # stream.ctx.bio = stream.sock
        # MbedTLS.set_bio!(stream.ctx,
        #                  pointer_from_objref(stream.sock),
        #                  cfunction(bio_write, Cint, (Ptr{Void}, Ptr{UInt8}, Csize_t)),
        #                  cfunction(bio_read, Cint, (Ptr{Void}, Ptr{UInt8}, Csize_t)))
        MbedTLS.handshake(stream.ctx)
    catch ex
        ret = if isa(ex, MbedTLS.MbedException)
            ssl_set_error(ex.ret)
        else
            set_error(LibGit2.Error.SSL, "$ex")
            Cint(-1)
        end
        return ret
    end

    return zero(Cint)
end

function tls_stream_certificate(out::Ptr{Ptr{Void}}, stream::Ptr{Void})
    # out[] = Ref(GitCert(CERT_X509))
    println("tls_stream_certificate: ")
    return zero(Cint)
end

function tls_stream_setproxy(stream::Ptr{Void}, proxyurl::Cstring)
    println("tls_stream_setproxy: ")
    return zero(Cint)
end

function tls_stream_read(tls_ptr::Ptr{TLSWrapper}, data::Ptr{UInt8}, len::Csize_t)
    stream = getstream(tls_ptr)
    ret = ccall((:mbedtls_ssl_read, MbedTLS.MBED_TLS), Cint,
                (Ptr{Void}, Ptr{Void}, Csize_t),
                stream.ctx.data, data, len)
    ret <= 0 && ssl_set_error(ret)
    return Cint(ret)
end

function tls_stream_write(tls_ptr::Ptr{TLSWrapper}, data::Ptr{UInt8}, len::Csize_t, flags::Cint)
    stream = getstream(tls_ptr)
    ret = ccall((:mbedtls_ssl_write, MbedTLS.MBED_TLS), Cint,
              (Ptr{Void}, Ptr{Void}, Csize_t),
              stream.ctx.data, data, len)
    ret <= 0 && ssl_set_error(ret)
    return Cint(ret)
end

function tls_stream_close(tls_ptr::Ptr{TLSWrapper})
    stream = getstream(tls_ptr)

    # close SSL context
    try
        if isopen(stream.sock)
            close(stream.ctx)
        end
    catch
        return Cint(-1)
    end

    # close socket
    Base.uvfinalize(stream.sock)

    # println("$(stream.sock)")

    return zero(Cint)
end

function tls_stream_free(tls_ptr::Ptr{TLSWrapper})
    stream = getstream(tls_ptr)

    # remove anchor
    pop_root(stream)

    # clear wrapper memory
    Libc.free(tls_ptr)

    return
end

function tls_stream(out::Ptr{Ptr{Void}}, host::Cstring, port::Cstring)

    # Create socket: socket.jl:243
    sock = TCPSocket(Libc.malloc(Base._sizeof_uv_tcp))
    Base.associate_julia_struct(sock.handle,sock)
    err = ccall(:uv_tcp_init,Cint,(Ptr{Void},Ptr{Void}),
                Base.eventloop(),sock.handle)
    if err != 0
        Libc.free(sock.handle)
        sock.handle = C_NULL
        return Cint(-1)
    end
    sock.status = Base.StatusInit

    # init ssl context
    ctx = MbedTLS.SSLContext()
    MbedTLS.setup!(ctx, tls_conf)

    io = TLSStream(
        sock,
        unsafe_string(host),
        unsafe_string(port),
        ctx,
        GitCertX509()
    )
    push_root(io)

    parent = GitStream(
        encrypted   = one(Cint),
        connect     = cfunction(tls_stream_connect, Cint, (Ptr{TLSWrapper},)),
        certificate = cfunction(tls_stream_certificate, Cint, (Ptr{Ptr{Void}}, Ptr{Void})),
        set_proxy   = cfunction(tls_stream_setproxy, Cint, (Ptr{Void}, Cstring)),
        read        = cfunction(tls_stream_read, Cint, (Ptr{TLSWrapper}, Ptr{UInt8}, Csize_t)),
        write       = cfunction(tls_stream_write, Cint, (Ptr{TLSWrapper}, Ptr{UInt8}, Csize_t, Cint)),
        close       = cfunction(tls_stream_close, Cint, (Ptr{TLSWrapper},)),
        free        = cfunction(tls_stream_free, Void, (Ptr{TLSWrapper},))
    )

    tls = TLSWrapper(
        parent,
        pointer_from_objref(io)
    )

    # duplicate tls stream structure (to avoid GC)
    data = memcpy(convert(Ptr{UInt8}, pointer_from_objref(tls)), sizeof(tls))
    unsafe_store!(out, data, 1)

    return zero(Cint)
end

function __init__()
    init_ssl()
    const tls_stream_cb = cfunction(tls_stream, Cint, (Ptr{Ptr{Void}}, Cstring, Cstring))
    ccall((:git_stream_register_tls, :libgit2), Cint, (Ptr{Void},), tls_stream_cb)
end

end # module


text:
module MachineGa

struct Chromosome
    genes::Array{Float64,1}
    cost::Float64
end

struct MCGA
    popsize::Int
    chsize::Int
    costfunction::Function
    mutationprobability::Float64
    crossoverprobability::Float64
    lowerbounds::Array{Float64,1}
    upperbounds::Array{Float64,1}
end

include("./byteworks.jl")
include("./rawoperators.jl")
include("./floatoperators.jl")
include("./gaoperators.jl")

# exports
export MCGA
export Chromosome

export ByteWorks
export RawOperators
export FloatOperators
export GaOperators

end # module


text:
tdf1 = TimeDataFrame(DataFrame(randn(4, 3), :auto), SemesterDate(2000,2))
tdf2 = TimeDataFrame(DataFrame(randn(4, 3), :auto), SemesterDate(2001,1))
tdf3 = TimeDataFrame(DataFrame(randn(4, 3), :auto), SemesterDate(1999,2))
tdf4 = TimeDataFrame(DataFrame(randn(5, 3), :auto), SemesterDate(2002,2))
tdf5 = TimeDataFrame(DataFrame(randn(5, 3), :auto), SemesterDate(2004,2))
tdf6 = TimeDataFrame(DataFrame(randn(5, 3), :auto), SemesterDate(1996,2))

tdf1_orig = copy(tdf1)
tdf2_orig = copy(tdf2)
tdf3_orig = copy(tdf3)
tdf4_orig = copy(tdf4)
tdf5_orig = copy(tdf5)
tdf6_orig = copy(tdf6)
TimeDataFrames.align!(tdf1, tdf1)
@test tdf1 == tdf1_orig

TimeDataFrames.align!(tdf1, tdf2)
@test size(tdf1) == (5, 3)
@test size(tdf2) == (5, 3)
@test periods(tdf1)[1] == SemesterDate(2000,2)
@test periods(tdf1)[end] == SemesterDate(2002,2)
@test periods(tdf2)[1] == SemesterDate(2000,2)
@test periods(tdf2)[end] == SemesterDate(2002,2)
@test tdf1[1, 1] == tdf1_orig[1, 1]
@test tdf1[4, 1] == tdf1_orig[4, 1]
@test ismissing.(tdf1[5, 1])
@test tdf2[2, 1] == tdf2_orig[1, 1]
@test tdf2[5, 1] == tdf2_orig[4, 1]
@test ismissing.(tdf2[1, 1])

tdf1 = copy(tdf1_orig)
tdf2 = copy(tdf2_orig)
tdf = tdf1.x1 + tdf2.x2
@test ismissing(tdf[1,1])
@test ismissing(tdf[5,1])
@test tdf[2,1] ≈ tdf1[2,1] + tdf2[1, 2]
@test tdf[4,1] ≈ tdf1[4,1] + tdf2[3, 2]

tdf1 = copy(tdf1_orig)
tdf2 = copy(tdf2_orig)
tdf = tdf1.x1 - tdf2.x2
@test ismissing(tdf[1,1])
@test ismissing(tdf[5,1])
@test tdf[2,1] ≈ tdf1[2,1] - tdf2[1, 2]
@test tdf[4,1] ≈ tdf1[4,1] - tdf2[3, 2]

x1 = tdf1[!, :x1]
lx1 = lag(x1)
@test ismissing(lx1[1])
@test lx1[2] == x1[1]
@test collect(skipmissing(lag(x1, 2))) == collect(skipmissing(lag(lag(x1))))

x1 = tdf1[!, :x1]
lx1 = lead(x1)
@test ismissing(lx1[end])
@test lx1[1] == x1[2]
@test collect(skipmissing(lead(x1, 2))) == collect(skipmissing(lead(lead(x1))))

tdf2 = index(tdf1, SemesterDate(2001,1))
df1 = dataframe(tdf1)
df2 = dataframe(tdf2)
for (i, c) in enumerate(eachcol(df2))
    @test c .* df1[2, i]/100 ≈ df1[!, i]
end