text:
@doc "PathLoss Model"
abstract Model

@doc "Pathloss Calculation for model"
function pathloss(model::Model,dist::AbstractFloat)

  return pathloss_op(model,dist)

end

@doc "Pathloss Calculation for model"
function pathloss(model::Model,dists::Vector{AbstractFloat})
    nRows = size(dists,1)
    resp::Vector{AbstractFloat}=Vector{AbstractFloat}(nRows)

    @simd for i in 1:nRows
      @inbounds  resp[i] = pathloss_op(model,dists[i])
    end
    return resp
end


text:
module stochastic_master

export master, master_dynamic

using ...bases, ...states, ...operators
using ...operators_dense, ...operators_sparse
using ...timeevolution
import ...timeevolution: integrate_stoch, recast!
import ...timeevolution.timeevolution_master: dmaster_h, dmaster_nh, dmaster_h_dynamic, check_master

const DecayRates = Union{Vector{Float64}, Matrix{Float64}, Void}

"""
    stochastic.master(tspan, rho0, H, J, Js; <keyword arguments>)

Time-evolution according to a stochastic master equation.

For dense arguments the `master` function calculates the
non-hermitian Hamiltonian and then calls master_nh which is slightly faster.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `H`: Deterministic part of the Hamiltonian.
* `J`: Vector containing all deterministic
        jump operators which can be of any arbitrary operator type.
* `Js`: Vector containing the stochastic jump operators for a superoperator
        describing a measurement which has the form of the standard linear
        stochastic master equation, `Js[i]*rho + rho*Jsdagger[i]`.
* `Hs=nothing`: Vector containing additional stochastic terms of the Hamiltonian.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `rates_s=nothing`: Vector specifying the coefficients (decay rates)
        for the stochastic jump operators. If nothing is specified all rates
        are assumed to be 1.
* `Jdagger=dagger.(J)`: Vector containing the hermitian conjugates of the jump
        operators. If they are not given they are calculated automatically.
* `Jsdagger=dagger.(Js)`: Vector containing the hermitian conjugates of the
        stochastic jump operators.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master(tspan, rho0::DenseOperator, H::Operator,
                J::Vector, Js::Vector; Hs::Union{Void, Vector}=nothing,
                rates::DecayRates=nothing, rates_s::DecayRates=nothing,
                Jdagger::Vector=dagger.(J), Jsdagger::Vector=dagger.(Js),
                fout::Union{Function,Void}=nothing,
                kwargs...)

    tmp = copy(rho0)

    if isa(rates_s, Matrix{Float64})
        throw(ArgumentError("A matrix of stochastic rates is ambiguous! Please provide a vector of stochastic rates.
        You may want to use diagonaljumps."))
    end

    n = length(Js) + (isa(Hs, Void) ? 0 : length(Hs))
    dmaster_stoch(t::Float64, rho::DenseOperator, drho::DenseOperator, index::Int) = dmaster_stochastic(rho, Hs, rates_s, Js, Jsdagger, drho, tmp, index)

    isreducible = check_master(rho0, H, J, Jdagger, rates)
    check_master(rho0, H, Js, Jsdagger, rates_s)
    if !isreducible
        dmaster_h_determ(t::Float64, rho::DenseOperator, drho::DenseOperator) = dmaster_h(rho, H, rates, J, Jdagger, drho, tmp)
        integrate_master_stoch(tspan, dmaster_h_determ, dmaster_stoch, rho0, fout, n; kwargs...)
    else
        Hnh = copy(H)
        if typeof(rates) == Matrix{Float64}
            for i=1:length(J), j=1:length(J)
                Hnh -= 0.5im*rates[i,j]*Jdagger[i]*J[j]
            end
        elseif typeof(rates) == Vector{Float64}
            for i=1:length(J)
                Hnh -= 0.5im*rates[i]*Jdagger[i]*J[i]
            end
        else
            for i=1:length(J)
                Hnh -= 0.5im*Jdagger[i]*J[i]
            end
        end
        Hnhdagger = dagger(Hnh)

        dmaster_nh_determ(t::Float64, rho::DenseOperator, drho::DenseOperator) = dmaster_nh(rho, Hnh, Hnhdagger, rates, J, Jdagger, drho, tmp)
        integrate_master_stoch(tspan, dmaster_nh_determ, dmaster_stoch, rho0, fout, n; kwargs...)
    end
end
master(tspan, psi0::Ket, args...; kwargs...) = master(tspan, dm(psi0), args...; kwargs...)

"""
    stochastic.master_dynamic(tspan, rho0, f, fs; <keyword arguments>)

Time-evolution according to a stochastic master equation with a
dynamic Hamiltonian and J.

# Arguments
* `tspan`: Vector specifying the points of time for which output should be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `fdeterm`: Function `f(t, rho) -> (H, J, Jdagger)` or
        `f(t, rho) -> (H, J, Jdagger, rates)` giving the deterministic
        part of the master equation.
* `fstoch`: Function `f(t, rho) -> (Js, Jsdagger)` or
        `f(t, rho) -> (Js, Jsdagger, rates)` giving the stochastic superoperator
        of the form `Js[i]*rho + rho*Jsdagger[i]`.
* `rates=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the jump operators. If nothing is specified all rates are assumed
        to be 1.
* `rates_s=nothing`: Vector or matrix specifying the coefficients (decay rates)
        for the stochastic jump operators. If nothing is specified all rates are assumed
        to be 1.
* `fstoch_H=nothing`: Function `f(t, rho) -> Hs` providing a vector of operators
        that correspond to stochastic terms of the Hamiltonian.
* `fstoch_J=nothing`: Function `f(t, rho) -> (J, Jdagger)` or
        `f(t, rho) -> (J, Jdagger, rates)` giving a stochastic
        Lindblad term.
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `noise_processes=0`: Number of distinct white-noise processes in the equation.
        This number has to be equal to the total number of noise operators
        returned by `fstoch` and all optional functions. If unset, the number
        is calculated automatically from the function outputs. NOTE: Set this
        number if you want to avoid an initial calculation of function outputs!
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master_dynamic(tspan::Vector{Float64}, rho0::DenseOperator, fdeterm::Function, fstoch::Function;
                fstoch_H::Union{Function, Void}=nothing, fstoch_J::Union{Function, Void}=nothing,
                rates::DecayRates=nothing, rates_s::DecayRates=nothing,
                fout::Union{Function,Void}=nothing,
                noise_processes::Int=0,
                kwargs...)

    tmp = copy(rho0)

    if isa(rates_s, Matrix{Float64})
        throw(ArgumentError("A matrix of stochastic rates is ambiguous! Please provide a vector of stochastic rates.
        You may want to use diagonaljumps."))
    end

    if noise_processes == 0
        fs_out = fstoch(0, rho0)
        n = length(fs_out[1])
        if isa(fstoch_H, Function)
            n += length(fstoch_H(0, rho0))
        end
        if isa(fstoch_J, Function)
            n += length(fstoch_J(0, rho0)[1])
        end
    else
        n = noise_processes
    end

    dmaster_determ(t::Float64, rho::DenseOperator, drho::DenseOperator) = dmaster_h_dynamic(t, rho, fdeterm, rates, drho, tmp)
    if isa(fstoch_H, Void) && isa(fstoch_J, Void)
        dmaster_stoch_std(t::Float64, rho::DenseOperator, drho::DenseOperator, index::Int) =
            dmaster_stoch_dynamic(t, rho, fstoch, rates_s, drho, tmp, index)

        integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_std, rho0, fout, n; kwargs...)
    else
        dmaster_stoch_gen(t::Float64, rho::DenseOperator, drho::DenseOperator, index::Int) =
            dmaster_stoch_dynamic_general(t, rho, fstoch, fstoch_H, fstoch_J,
                    rates, rates_s, drho, tmp, index)

        integrate_master_stoch(tspan, dmaster_determ, dmaster_stoch_gen, rho0, fout, n; kwargs...)
    end
end
master_dynamic(tspan::Vector{Float64}, psi0::Ket, args...; kwargs...) = master_dynamic(tspan, dm(psi0), args...; kwargs...)

function dmaster_stochastic(rho::DenseOperator, H::Void, rates::Void,
            J::Vector, Jdagger::Vector, drho::DenseOperator, tmp::DenseOperator,
            index::Int)
    operators.gemm!(1, J[index], rho, 0, drho)
    operators.gemm!(1, rho, Jdagger[index], 1, drho)
    return drho
end
function dmaster_stochastic(rho::DenseOperator, H::Void, rates::Vector{Float64},
            J::Vector, Jdagger::Vector, drho::DenseOperator, tmp::DenseOperator,
            index::Int)
    operators.gemm!(rates[index], J[index], rho, 0, drho)
    operators.gemm!(rates[index], rho, Jdagger[index], 1, drho)
    return drho
end

function dmaster_stochastic(rho::DenseOperator, H::Vector, rates::Void,
            J::Vector, Jdagger::Vector, drho::DenseOperator, tmp::DenseOperator,
            index::Int)
    if index > length(J)
        operators.gemm!(-1.0im, H[index - length(J)], rho, 0.0, drho)
        operators.gemm!(1.0im, rho, H[index - length(J)], 1.0, drho)
    else
        operators.gemm!(1, J[index], rho, 0, drho)
        operators.gemm!(1, rho, Jdagger[index], 1, drho)
    end
    return drho
end
function dmaster_stochastic(rho::DenseOperator, H::Vector, rates::Vector{Float64},
            J::Vector, Jdagger::Vector, drho::DenseOperator, tmp::DenseOperator,
            index::Int)
    if index > length(J)
        operators.gemm!(-1.0im, H[index - length(J)], rho, 0.0, drho)
        operators.gemm!(1.0im, rho, H[index - length(J)], 1.0, drho)
    else
        operators.gemm!(rates[index], J[index], rho, 0, drho)
        operators.gemm!(rates[index], rho, Jdagger[index], 1, drho)
    end
    return drho
end

function dmaster_stoch_dynamic(t::Float64, rho::DenseOperator, f::Function, rates::DecayRates,
            drho::DenseOperator, tmp::DenseOperator, index::Int)
    result = f(t, rho)
    @assert 2 <= length(result) <= 3
    if length(result) == 2
        J, Jdagger = result
        rates_ = rates
    else
        J, Jdagger, rates_ = result
    end
    dmaster_stochastic(rho, nothing, rates_, J, Jdagger, drho, tmp, index)
end

function dmaster_stoch_dynamic_general(t::Float64, rho::DenseOperator, fstoch::Function,
            fstoch_H::Function, fstoch_J::Void, rates::DecayRates, rates_s::DecayRates,
            drho::DenseOperator, tmp::DenseOperator, index::Int)
    H = fstoch_H(t, rho)
    if index <= length(H)
        operators.gemm!(-1.0im, H[index], rho, 0.0, drho)
        operators.gemm!(1.0im, rho, H[index], 1.0, drho)
    else
        dmaster_stoch_dynamic(t, rho, fstoch, rates_s, drho, tmp, index-length(H))
    end
end
function dmaster_stoch_dynamic_general(t::Float64, rho::DenseOperator, fstoch::Function,
            fstoch_H::Void, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates,
            drho::DenseOperator, tmp::DenseOperator, index::Int)
    result_J = fstoch_J(t, rho)
    if index <= length(result_J[1])
        @assert 2 <= length(result_J) <= 3
        if length(result_J) == 2
            J_stoch, J_stoch_dagger = result_J
            rates_ = rates
        else
            J_stoch, J_stoch_dagger, rates_ = result_J
        end
        dlindblad(rho, rates_, J_stoch, J_stoch_dagger, drho, tmp, index)
    else
        dmaster_stoch_dynamic(t, rho, fstoch, rates_s, drho, tmp, index-length(result_J[1]))
    end
end
function dmaster_stoch_dynamic_general(t::Float64, rho::DenseOperator, fstoch::Function,
            fstoch_H::Function, fstoch_J::Function, rates::DecayRates, rates_s::DecayRates,
            drho::DenseOperator, tmp::DenseOperator, index::Int)
    H = fstoch_H(t, rho)
    result_J = fstoch_J(t, rho)
    if index <= length(H)
        operators.gemm!(-1.0im, H[index], rho, 0.0, drho)
        operators.gemm!(1.0im, rho, H[index], 1.0, drho)
    elseif length(H) < index <= length(H) + length(result_J[1])
        @assert 2 <= length(result_J) <= 3
        if length(result_J) == 2
            J_stoch, J_stoch_dagger = result_J
            rates_ = rates
        else
            J_stoch, J_stoch_dagger, rates_ = result_J
        end
        dlindblad(rho, rates_, J_stoch, J_stoch_dagger, drho, tmp, index-length(H))
    else
        dmaster_stoch_dynamic(t, rho, fstoch, rates_s, drho, tmp, index-length(H)-length(result_J[1]))
    end
end

function dlindblad(rho::DenseOperator, rates::Void, J::Vector, Jdagger::Vector,
    drho::DenseOperator, tmp::DenseOperator, i::Int)
    operators.gemm!(1, J[i], rho, 0, tmp)
    operators.gemm!(1, tmp, Jdagger[i], 1, drho)
    return drho
end
function dlindblad(rho::DenseOperator, rates::Vector{Float64}, J::Vector,
    Jdagger::Vector, drho::DenseOperator, tmp::DenseOperator, i::Int)
    operators.gemm!(rates[i], J[i], rho, 0, tmp)
    operators.gemm!(1, tmp, Jdagger[i], 1, drho)
    return drho
end

function integrate_master_stoch(tspan, df::Function, dg::Function,
                        rho0::DenseOperator, fout::Union{Void, Function},
                        n::Int;
                        kwargs...)
    tspan_ = convert(Vector{Float64}, tspan)
    x0 = reshape(rho0.data, length(rho0))
    state = copy(rho0)
    dstate = copy(rho0)
    integrate_stoch(tspan_, df, dg, x0, state, dstate, fout, n; kwargs...)
end

function recast!(x::SubArray{Complex128, 1}, rho::DenseOperator)
    rho.data = reshape(x, size(rho.data))
end
recast!(state::DenseOperator, x::SubArray{Complex128, 1}) = (x[:] = state.data)

end # module

ext: jl


text:
#=
-------------------------------------------------------------------------------------------------------------
This code computes the 1D Poisson Equation using the CG method.

This is the solution to Project 5a in Element-based Galerkin Methods
The algorithms here are described in Algorithms 8.1 and 8.2, and DSSed
using Alg. 5.4. in F.X. Giraldo's Introduction to Element-based Galerkin
Methods using Tensor-Product Bases: Analysis, Algorithms, and Applications.

Written by F.X. Giraldo on July 5, 2021.
Department of Applied Mathematics
Naval Postgraduate School
Monterey, CA 93943

The interpolation points used are the following:
ipoints=1: Lobatto
ipoints=2: Legendre
ipoints=3: Chebyshev
ipoints=4: Equi-spaced

The integration points used are the following:
qpoints=1: Lobatto
qpoints=2: Legendre
-------------------------------------------------------------------------------------------------------------
=#

using Plots

#Include all files used here
include("QuadraturePoints.jl")
include("create_grid.jl")
include("global_matrices.jl")
include("exact_solution.jl")
include("compute_norms.jl")

#Some Constants
DFloat = Float64
machine_zero=eps(DFloat)

function main()

    #-----------------------------Only Change these Input parameters---------------------------------#
    Ne=10 #Number of elements
    N=4 #Polynomial Order
    Q=N+1 #Quadrature Order
    ipoints=1 #interpolation point type
    qpoints=1 #integration point type
    case=1
    #-----------------------------Only Change these Input parameters---------------------------------#

    Np=N+1
    Nq=Q+1
    Npoin_dg=Ne*Np
    Npoin_cg=Ne*N + 1
    space_method="CG" #only capable of doing CG
    @show(N,Q,Ne,Npoin_cg,Npoin_dg,case)

    #Select Interpolation Points
    if ipoints == 1
        (ξ,ω) = QuadraturePoints.lobatto_gauss(Np,DFloat) #Lobatto
    elseif ipoints == 2
        (ξ,ω) = QuadraturePoints.legendre_gauss(Np,DFloat) #Legendre
    elseif ipoints == 3
        (ξ,ω) = QuadraturePoints.chebyshev_gauss(Np,DFloat) #Chebyshev
    elseif ipoints == 4
        (ξ,ω) = QuadraturePoints.equispaced_points(Np,DFloat) #Equi-spaced
    end

    #Select Quadrature Points
    if qpoints == 1
        (ξq,ωq) = QuadraturePoints.lobatto_gauss(Nq,DFloat) #Lobatto
    elseif qpoints == 2
        (ξq,ωq) = QuadraturePoints.legendre_gauss(Nq,DFloat) #Legendre
    end

    #Construct Lagrange Polynomials
    (ψ,dψ) = QuadraturePoints.lagrange_basis(Np,Nq,ξ,ξq,DFloat)

    #Construct Grid and (i,e)->I map
    (coord_cg,coord_dg,intma_cg,intma_dg,periodicity_cg,periodicity_dg) = create_grid(Np,Npoin_cg,Npoin_dg,Ne,ξ,DFloat)

    #CG or DG?
    if (space_method == "CG")
        Npoin=Npoin_cg
        coord=copy(coord_cg)
        intma=copy(intma_cg)
        periodicity=copy(periodicity_cg)
    elseif (space_method == "DG")
        Npoin=Npoin_dg
        coord=copy(coord_dg)
        intma=copy(intma_dg)
        periodicity=copy(periodicity_dg)
    end
    @show(space_method,Npoin,Ne,N,Q)

    #Construct Global Mass and Differentiation Matrices (only need Mass)
    #-----------------------------Students Add their Function here--------------------------#
    #-----------------------------Students Add their Function here--------------------------#
    (M,L) = global_matrices(ψ,dψ,ωq,intma,coord,Npoin,Ne,Np,Nq,DFloat)
    #-----------------------------Students Add their Function here--------------------------#
    #-----------------------------Students Add their Function here--------------------------#

    #Compute Exact Solution
    (qe,qeₓ,fe) = exact_solution(coord,Npoin,case,DFloat)

    #Build RHS vector and apply Dirichlet BC
    R=M*fe
    R[1]=qe[1]
    R[Npoin]=qe[Npoin]

    #Compute CG Solution
    q0=-L\R

    #Compute L2 Norm
    norms = compute_norms(q0,qe,Npoin,DFloat)
    println(" norm1 = ",norms[1]," norm2 = ",norms[2]," norm∞ = ",norms[3])

    #Plot Interpolation
    p1 = plot(coord,qe,linecolor = :blue,label = "Analytic Solution",legend = :bottomright,legendfont=font(5))
    plot!(p1,coord,q0,linestyle = :dashdot,linecolor = :red,label = "$space_method Solution")
    title!(p1,"1D Poisson Problem for Ne= $Ne, N=$N, Q=$Q")
    xlabel!(p1,"x")
    ylabel!(p1,"q(x)")
    plot!(p1,xlims = (first(coord),last(coord)+.1),ylims = (-1,1))

    figure_1 = plot(p1)
    display(figure_1)

#    savefig(plot_handle,"Project2.png")
    println("**Simulation Finished**") #output
end

#----------------------------------#
# Run the main function
#----------------------------------#
main()


text:
# FASTQ
# =====
#
# Reader and writer of the FASTQ file format.
#
# Cock, Peter JA, et al. "The Sanger FASTQ file format for sequences with
# quality scores, and the Solexa/Illumina FASTQ variants." Nucleic acids
# research 38.6 (2010): 1767-1771.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

"""
Metadata for FASTQ sequence records containing a `description` field,
and a `quality` string corresponding to the sequence.

Quality scores are stored as integer Phred scores.
"""
type FASTQMetadata
    description::StringField
    quality::Vector{Int8}
end

function FASTQMetadata()
    return FASTQMetadata(StringField(), Int8[])
end

function Base.:(==)(a::FASTQMetadata, b::FASTQMetadata)
    return a.description == b.description && a.quality == b.quality
end

function Base.copy(metadata::FASTQMetadata)
    return FASTQMetadata(copy(metadata.description), copy(metadata.quality))
end

"""
    FASTQSeqRecord(name, seq, quality[, description=""])

Create a sequence record for the FASTQ file format.
"""
typealias FASTQSeqRecord{S} SeqRecord{S,FASTQMetadata}

function (::Type{FASTQSeqRecord})(name::AbstractString,
                                  seq::Sequence,
                                  quality::Vector{Int8},
                                  description::AbstractString="")
    if length(seq) != length(quality)
        throw(ArgumentError("sequence and base quality must be the same length"))
    end
    return SeqRecord(name, seq, FASTQMetadata(description, quality))
end

include("quality.jl")
include("reader.jl")
include("parser.jl")
include("writer.jl")


text:
baremodule Restarts

export  restart, @define_restart

function with end
function restart end
macro define_restart end

module Internal

using ..Restarts: Restarts
import ..Restarts: @define_restart

using ContextVariablesX: @contextvar, with_context

include("internal.jl")

end  # module Internal

end  # baremodule Restarts


text:
"""
    _group_by_friends!(ata_model::AbstractModel)

# Description

Group the items by friend sets once the friend sets have been added to the `ATA.AbstractModel` by [`add_friends!`](#ATA.add_friends!-Tuple{ATA.AbstractModel}) function and as **last operation** of the assembly.

# Arguments

- **`ata_model::AbstractModel`** : Required. The model built with `start_ata()`, settings loaded by [`start_ata`](#ATA.start_ata) function and friend sets added by [`add_friends!`](#ATA.add_friends!-Tuple{ATA.AbstractModel}) function.
"""
function _group_by_friends!(ata_model::AbstractModel) #last
    try
        n_items = ata_model.settings.n_items
        #only works for categorical variables and item use, all the other contraitns need expansion by fs_items
        if ata_model.settings.n_fs == 0
            error!(ata_model, "No friend sets found, run start_ata() before.")
            return nothing
        end
        if any([
            size(ata_model.constraints[t].constr_b, 1) > 0 for t = 1:ata_model.settings.T
        ])
            A = Vector{Matrix{Float64}}(undef, ata_model.settings.T)
            b = Vector{Vector{Float64}}(undef, ata_model.settings.T)
            A_new = Vector{Matrix{Float64}}(undef, ata_model.settings.T)
            for t = 1:ata_model.settings.T
                A[t] = ata_model.constraints[t].constr_A
                b[t] = ata_model.constraints[t].constr_b
            end
            n_fs = ata_model.settings.n_fs
            for t = 1:ata_model.settings.T
                A_new[t] = Matrix{Float64}(undef, size(A[t], 1), n_fs)
            end
            for fs = 1:n_fs
                for t = 1:ata_model.settings.T
                    A_new[t][:, fs] =
                        sum(A[t][:, ata_model.settings.fs.items[fs]], dims = 2)
                end
            end
            for t = 1:ata_model.settings.T
                DelimitedFiles.writedlm("opt/A_$t.csv", A_new[t])
                DelimitedFiles.writedlm("opt/b_$t.csv", b[t])
                ata_model.constraints[t].constr_A = A_new[t]
            end
            open("opt/fs_items.jl", "w") do f
                write(f, "fs_items = Vector{Vector{Int64}}(undef, $n_fs)\n")
                for fs = 1:n_fs
                    write(
                        f,
                        string("fs_items[$fs] =", ata_model.settings.fs.items[fs], "\n"),
                    )
                end
            end
            # open("opt/Settings.jl", "a") do f
            #     write(f, "group_by_fs = true\n\n")
            # end
            #transform forced0
            x_forced0 = ata_model.settings.forced0
            x_forced0_new = Vector{Vector{Bool}}(undef, ata_model.settings.T)
            for v = 1:ata_model.settings.T
                x_forced0_new[v] = fill(true, n_fs)
                for i = 1:ata_model.settings.n_items
                    if x_forced0[v][i] == false
                        for fs = 1:n_fs
                            if any(i .== ata_model.settings.fs.items[fs])
                                x_forced0_new[v][fs] = false
                            end
                        end
                    end
                end
            end
            #update ata_model.forced0
            ata_model.settings.forced0 = x_forced0_new
            DelimitedFiles.writedlm("opt/x_forced0.txt", x_forced0_new)
            success!(
                ata_model,
                "Categorical and linear quantitative constraints grouped by friend sets.",
            )
        end
        #item use
        if ata_model.settings.to_apply[1]
            item_use_max = ata_model.settings.iu.max
            item_use_max_new = zeros(Int, n_fs)
            for fs = 1:n_fs
                item_use_max_new[fs] =
                    Int(minimum(ata_model.settings.iu.max[ata_model.settings.fs.items[fs]]))
            end
            ata_model.settings.iu.max = item_use_max_new
            open("opt/Settings.jl", "a") do f
                write(f, "item_use_max = $item_use_max_new\n\n")
            end
            success!(ata_model, "Maximum item use grouped.")
        end
        if ata_model.settings.to_apply[2]
            item_use_min = ata_model.settings.iu.min
            item_use_min_new = zeros(Int, n_fs)
            for fs = 1:n_fs
                item_use_min_new[fs] =
                    Int(maximum(ata_model.settings.iu.min[ata_model.settings.fs.items[fs]]))
            end
            ata_model.settings.iu.min = item_use_min_new
            open("opt/Settings.jl", "a") do f
                write(f, "item_use_min = $item_use_min_new\n\n")
            end
            success!(ata_model, "Minimum item use grouped.")
        end
    catch e
        error!(ata_model, string("", sprint(showerror, e)))
    end
    return nothing
end


text:
# Plot the visibilities

module plot_vis

export plot_RawModelVis

using visibilities
using constants

import PyPlot
import PyPlot.plt
using LaTeXStrings

FSF(str) = PyPlot.matplotlib[:ticker][:FormatStrFormatter](str)

function plot_RawModelVis(vis::RawModelVis)
    fig, ax = plt.subplots(nrows=2, figsize=(5, 8))

    FullVis = fillModelVis(vis)
    VV = FullVis.VV
    uu = FullVis.uu
    vv = FullVis.vv

    rvis = real(VV)
    cvis = imag(VV)

    ext = (minimum(uu), maximum(uu), minimum(vv), maximum(vv))

    rimg = ax[1][:imshow](rvis, interpolation="none", origin="upper", cmap=plt.get_cmap("Greys"), extent=ext)
    ax[1][:set_title]("real")
    ax[1][:set_xlabel](L"uu [k$\lambda$]")
    ax[1][:set_ylabel](L"vv [k$\lambda$]")
    ax[1][:xaxis][:set_major_formatter](FSF("%.0f"))
    ax[1][:yaxis][:set_major_formatter](FSF("%.0f"))

    #[left, bottom, width, height]
    cax = fig[:add_axes]([0.84, 0.65, 0.03, 0.25])
    cb = fig[:colorbar](rimg, cax=cax)

        #ticks = np.linspace(0, np.max(cov), num=6)
        #cb.set_ticks(ticks)
        #cb.set_ticks(MaxNLocator(nbins=5))

    cimg = ax[2][:imshow](cvis, interpolation="none", origin="upper", cmap=plt.get_cmap("Greys"), extent=ext)
    ax[2][:set_title]("imag")
    ax[2][:set_xlabel](L"uu [k$\lambda$]")
    ax[2][:set_ylabel](L"vv [k$\lambda$]")
    ax[2][:xaxis][:set_major_formatter](FSF("%.0f"))
    ax[2][:yaxis][:set_major_formatter](FSF("%.0f"))

    cax = fig[:add_axes]([0.84, 0.15, 0.03, 0.25])
    cb = fig[:colorbar](cimg, cax=cax)

    fig[:subplots_adjust](left=0.15, right=0.85, hspace=0.2)
    plt.savefig("plots/visibilities.png")

end

end

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IPConfigurationPropertiesFormat <: SwaggerModel
    privateIPAddress::Any # spec type: Union{ Nothing, String } # spec name: privateIPAddress
    privateIPAllocationMethod::Any # spec type: Union{ Nothing, String } # spec name: privateIPAllocationMethod
    subnet::Any # spec type: Union{ Nothing, Subnet } # spec name: subnet
    publicIPAddress::Any # spec type: Union{ Nothing, PublicIPAddress } # spec name: publicIPAddress
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function IPConfigurationPropertiesFormat(;privateIPAddress=nothing, privateIPAllocationMethod=nothing, subnet=nothing, publicIPAddress=nothing, provisioningState=nothing)
        o = new()
        validate_property(IPConfigurationPropertiesFormat, Symbol("privateIPAddress"), privateIPAddress)
        setfield!(o, Symbol("privateIPAddress"), privateIPAddress)
        validate_property(IPConfigurationPropertiesFormat, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        setfield!(o, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        validate_property(IPConfigurationPropertiesFormat, Symbol("subnet"), subnet)
        setfield!(o, Symbol("subnet"), subnet)
        validate_property(IPConfigurationPropertiesFormat, Symbol("publicIPAddress"), publicIPAddress)
        setfield!(o, Symbol("publicIPAddress"), publicIPAddress)
        validate_property(IPConfigurationPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type IPConfigurationPropertiesFormat

const _property_map_IPConfigurationPropertiesFormat = Dict{Symbol,Symbol}(Symbol("privateIPAddress")=>Symbol("privateIPAddress"), Symbol("privateIPAllocationMethod")=>Symbol("privateIPAllocationMethod"), Symbol("subnet")=>Symbol("subnet"), Symbol("publicIPAddress")=>Symbol("publicIPAddress"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_IPConfigurationPropertiesFormat = Dict{Symbol,String}(Symbol("privateIPAddress")=>"String", Symbol("privateIPAllocationMethod")=>"String", Symbol("subnet")=>"Subnet", Symbol("publicIPAddress")=>"PublicIPAddress", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ IPConfigurationPropertiesFormat }) = collect(keys(_property_map_IPConfigurationPropertiesFormat))
Swagger.property_type(::Type{ IPConfigurationPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IPConfigurationPropertiesFormat[name]))}
Swagger.field_name(::Type{ IPConfigurationPropertiesFormat }, property_name::Symbol) =  _property_map_IPConfigurationPropertiesFormat[property_name]

const _allowed_IPConfigurationPropertiesFormat_privateIPAllocationMethod = ["Static", "Dynamic"]

function check_required(o::IPConfigurationPropertiesFormat)
    true
end

function validate_property(::Type{ IPConfigurationPropertiesFormat }, name::Symbol, val)
    if name === Symbol("privateIPAllocationMethod")
        Swagger.validate_param(name, "IPConfigurationPropertiesFormat", :enum, val, _allowed_IPConfigurationPropertiesFormat_privateIPAllocationMethod)
    end
end


text:
#################### Chains ####################

## Constructors ##

# Default name map.
const DEFAULT_MAP = Dict{Symbol, Vector{Any}}(:parameters => [])

# Constructor to handle a vector of vectors.
Chains(val::AbstractVector{<:AbstractVector{<:Union{Missing, Real}}}, args...; kwargs...) =
	Chains(copy(reduce(hcat, val)'), args...; kwargs...)

# Constructor to handle a 1D array.
Chains(val::AbstractVector{<:Union{Missing, Real}}, args...; kwargs...) =
	Chains(reshape(val, :, 1, 1), args...; kwargs...)

# Constructor to handle a 2D array
Chains(val::AbstractMatrix{<:Union{Missing, Real}}, args...; kwargs...) =
	Chains(reshape(val, size(val, 1), size(val, 2), 1), args...; kwargs...)

# Generic chain constructor.
function Chains(
        val::AbstractArray{<:Union{Missing, Real},3},
        parameter_names::Vector{String} = map(i->"Param$i", 1:size(val, 2)),
        name_map_original = copy(DEFAULT_MAP);
        start::Int=1,
        thin::Int=1,
        evidence = missing,
        info::NamedTuple=NamedTuple(),
        sorted::Bool=false)
    # If we received an array of pairs, convert it to a dictionary.
    name_map = if typeof(name_map_original) <: Dict
        # Copying can avoid state mutation.
        deepcopy(name_map_original)
    elseif typeof(name_map_original) <: Array
        Dict(deepcopy(name_map_original))
    elseif typeof(name_map_original) <: NamedTuple
        _namedtuple2dict(name_map_original)
    end

    # Make sure that we have a :parameters index.
    if !in(:parameters, keys(name_map))
        name_map[:parameters] = []
    end

    # Preclean the name_map of names that aren't in the
    # parameter_names vector.
    for section in keys(name_map)
        filter!(x -> x ∈ parameter_names, name_map[section])
    end

    if length(keys(name_map)) == 1
        name_map[first(keys(name_map))] = parameter_names
    else
        # Store unassigned variables.
        unassigned = Set([])

        # Check that all parameters are assigned.
        for param in parameter_names
            found = false
            for (_, assigned) in name_map
                if param in assigned
                    found = true
                    break
                end
            end

            # Assign to :parameters by default, or :internals if it starts
            # or ends with an underscore.
            if !found
                push!(unassigned, param)
            end
        end

        # Assign all unassigned parameter names.
        for param in unassigned
            push!(name_map[:parameters], param)
        end
    end

    # Make the name_map NamedTuple.
    name_map_tupl = _dict2namedtuple(name_map)

    # Construct the AxisArray.
    arr = AxisArray(val;
                    iter = range(start, step=thin, length=size(val, 1)),
                    var = parameter_names,
                    chain = 1:size(val, 3))

    if sorted
        return sort(
            Chains{eltype(val), typeof(evidence), typeof(name_map_tupl), typeof(info)}(
                arr, evidence, name_map_tupl, info)
        )
    else
        return Chains{eltype(val), typeof(evidence), typeof(name_map_tupl), typeof(info)}(
                arr, evidence, name_map_tupl, info)
    end
end

# Retrieve a new chain with only a specific section pulled out.
Chains(::Chains, ::Tuple{}; kwargs...) = throw(ArgumentError("no section specified"))
Chains(c::Chains, section; kwargs...) = Chains(c, (section,); kwargs...)
function Chains(c::Chains, section::AbstractArray; kwargs...)
	Chains(c, ntuple(i -> section[i], length(section)); kwargs...)
end
function Chains(c::Chains, section::NTuple{<:Any,Symbol}; sorted::Bool=false)
    # Make sure the section exists first.
	section ⊆ keys(c.name_map) ||
		throw(ArgumentError("$section is not a subset of the chain's name map"))

	# Create the new section map.
    name_map = NamedTuple{section}(getindex.(Ref(c.name_map), section))

    # Extract wanted values.
    value = c.value[:, mapreduce(collect, vcat, name_map), :]

    # Create the new chain.
    chain = Chains{eltype(c.value),typeof(c.logevidence),typeof(name_map),typeof(c.info)}(
		value,
        c.logevidence,
        name_map,
        c.info)

    if sorted
        return sort(chain)
    else
        return chain
    end
end

#################### Indexing ####################
function _sym2index(c::Chains, v::Union{Vector{Symbol}, Vector{String}}; sorted::Bool = false)
    v_str = string.(v)
    idx = indexin(v_str, names(c))
    syms = []
    for i in eachindex(idx)
        value = v_str[i]
        if idx[i] == nothing
            append!(syms,
                collect(Iterators.filter(k -> startswith(string(k), value*"["), names(c))))
        else
            push!(syms, value)
        end
    end
    return sorted ? sort!(syms, lt=natural) : syms
end

Base.getindex(c::Chains, i1::T) where T<:Union{AbstractUnitRange, StepRange} = c[i1, :, :]
Base.getindex(c::Chains, i1::Integer) = c[i1:i1, :, :]
Base.getindex(c::Chains, v::Symbol) = c[[v]]
Base.getindex(c::Chains, v::String) = c[:, [v], :]
Base.getindex(c::Chains, v::Vector{String}) = c[:, v, :]

function Base.getindex(c::Chains, v::Vector{Symbol})
    syms = _sym2index(c, v)
    return c[:, syms, :]
end

function Base.getindex(c::Chains, i...)
    # Make sure things are in array form to preserve the axes.
    ind = [typeof(i[1]) <: Integer ? (i[1]:i[1]) : i[1],
           typeof(i[2]) <: Union{AbstractArray, Colon} ?  i[2] : [i[2]],
           typeof(i[3]) <: Union{AbstractArray, Colon} ?  i[3] : [i[3]]
    ]

    # Check to see if we received a symbol in i[2].
    if ind[2] != Colon() && typeof(first(ind[2])) <: Symbol
        ind[2] = _sym2index(c, ind[2])
    end

    newval = getindex(c.value, ind...)
    names = newval.axes[2].val
    new_name_map = _trim_name_map(names, c.name_map)
    return Chains{eltype(c.value),typeof(c.logevidence),typeof(new_name_map),
        typeof(c.info)}(newval, c.logevidence, new_name_map, c.info)
end

Base.setindex!(c::Chains, v, i...) = setindex!(c.value, v, i...)
Base.lastindex(c::Chains) = lastindex(c.value, 1)
Base.lastindex(c::Chains, d::Integer) = lastindex(c.value, d)

"""
    Base.get(c::Chains, v::Symbol; flatten=false)
    Base.get(c::Chains, vs::Vector{Symbol}; flatten=false)

Returns a `NamedTuple` with `v` as the key, and matching paramter
names as the values.

Passing `flatten=true` will return a `NamedTuple` with keys ungrouped.

Example:

```julia
x = get(c, :param1)
x = get(c, [:param1, :param2])
```
"""
Base.get(c::Chains, v::Symbol; flatten = false) = get(c, [v], flatten=flatten)
function Base.get(c::Chains, vs::Vector{Symbol}; flatten = false)
    pairs = Dict()
    for v in vs
        syms = _sym2index(c, [v])
        len = length(syms)
        val = ()
        if len > 1
            val = ntuple(i -> c.value[:,syms[i],:], length(syms))
        elseif len == 1
            val = c.value[:,syms[1],:]
        else
            continue
        end

        if flatten
            for i in eachindex(syms)
                pairs[syms[i]] = val[i]
            end
        else
            pairs[v] = val
        end
    end
    return _dict2namedtuple(pairs)
end

"""
    get(c::Chains; section::Union{Vector{Symbol}, Symbol; flatten=false}

Returns all parameters in a given section(s) as a `NamedTuple`.

Passing `flatten=true` will return a `NamedTuple` with keys ungrouped.

Example:

```julia
x = get(chn, section = :parameters)
x = get(chn, section = [:internals, :parameters])
```
"""
function Base.get(c::Chains;
        section::Union{Vector{Symbol}, Symbol},
        flatten = false)
    section = section isa Symbol ? [section] : section
    not_found = Symbol[]
    names = Set(String[])
    for v in section
        if v in keys(c.name_map)
            # If the name contains a bracket,
            # split it so get can group them correctly.
            nms = flatten ?
                c.name_map[v] :
                map(n -> String(split(n, "[")[1]), c.name_map[v])
            push!(names, nms...)
        else
            push!(not_found, v)
        end
    end

    if length(not_found) > 0
        throw(ArgumentError("$not_found not found in chains name map."))
    end

    return get(c, Symbol.(names), flatten = flatten)
end

"""
    get_params(c::Chains; flatten=false)

Returns all parameters packaged as a `NamedTuple`. Variables with a bracket
in their name (as in "P[1]") will be grouped into the returned value as P.

Passing `flatten=true` will return a `NamedTuple` with keys ungrouped.

Example:

```julia
x = get_params(chn)
x.P
```
"""
get_params(c::Chains; flatten = false) = get(c, section = sections(c), flatten=flatten)

#################### Base Methods ####################

function Base.show(io::IO, c::Chains)
    print(io, "Object of type Chains, with data of type $(summary(c.value.data))\n\n")
    println(io, header(c))

    # Show summary stats.
    show(io, describe(c))
end

Base.keys(c::Chains) = names(c)
Base.size(c::Chains) = size(c.value)
Base.size(c::Chains, ind) = size(c)[ind]
Base.length(c::Chains) = length(range(c))
Base.first(c::Chains) = first(c.value[Axis{:iter}].val)
Base.step(c::Chains) = step(c.value[Axis{:iter}].val)
Base.last(c::Chains) = last(c.value[Axis{:iter}].val)

Base.convert(::Type{Array}, chn::Chains) = convert(Array, chn.value)

#################### Auxilliary Functions ####################

"""
    range(c::Chains)

Returns the range used in a `Chains` object.
"""
function Base.range(c::Chains)
    return c.value[Axis{:iter}].val
end

"""
    chains(c::Chains)

Returns the names or symbols of each chain in a `Chains` object.
"""
function chains(c::Chains)
    return c.value[Axis{:chain}].val
end

"""
    names(c::Chains, sections)

Return the parameter names in a `Chains` object.
"""
function Base.names(c::Chains)
    return c.value[Axis{:var}].val
end

"""
    names(c::Chains, sections::Union{Symbol, Vector{Symbol}})

Return the parameter names in a `Chains` object, given an array of sections.
"""
function Base.names(c::Chains, sections::Union{Symbol, Vector{Symbol}})
    # Check that sections is an array.
    sections = typeof(sections) <: AbstractArray ?
        sections :
        [sections]

    nms = []

    for i in sections
        push!(nms, c.name_map[i]...)
    end
    return nms
end

"""
    get_sections(c::Chains, sections::Vector = [])

Returns multiple `Chains` objects, each containing only a single section.
"""
function get_sections(c::Chains, sections::Vector = [])
    sections = length(sections) == 0 ? collect(keys(c.name_map)) : sections
    return [Chains(c, section) for section in sections]
end

# Return a new chain for each section.
function get_sections(c::Chains, section::Union{Symbol, String})
    return get_sections(c, [section])
end

"""
    sections(c::Chains)

Retrieve a list of the sections in a chain.
"""
sections(c::Chains) = collect(keys(c.name_map))

"""
    header(c::Chains; section=missing)

Returns a string containing summary information for a `Chains` object.
If the `section` keyword is used, this function prints only the relevant section
header.

Example:

```julia
# Printing the whole header.
header(chn)

# Print only one section's header.
header(chn, section = :parameter)
```
"""
function header(c::Chains; section=missing)
    rng = range(c)

    # Function to make section strings.
    section_str(sec, arr) = string(
        "$sec",
        repeat(" ", 18 - length(string(sec))),
        "= $(join(map(string, arr), ", "))\n"
    )

    # Set up string array.
    section_strings = String[]

    # Get section lines.
    if section isa Missing
        for (sec, nms) in pairs(c.name_map)
            section_string = section_str(sec, nms)
            push!(section_strings, section_string)
        end
    else
        section in keys(c.name_map) ||
            throw(ArgumentError("$section not found in name map."))
        section_string = section_str(section, c.name_map[section])
        push!(section_strings, section_string)
    end

    #

    # Return header.
    return string(
        ismissing(c.logevidence) ? "" : "Log evidence      = $(c.logevidence)\n",
        "Iterations        = $(first(c)):$(last(c))\n",
        "Thinning interval = $(step(c))\n",
        "Chains            = $(join(map(string, chains(c)), ", "))\n",
        "Samples per chain = $(length(range(c)))\n",
        section_strings...
    )
end

function indiscretesupport(c::Chains,
                           bounds::Tuple{Real, Real}=(0, Inf))
  nrows, nvars, nchains = size(c.value)
  result = Array{Bool}(undef, nvars * (nrows > 0))
  for i in 1:nvars
    result[i] = true
    for j in 1:nrows, k in 1:nchains
      x = c.value[j, i, k]
      if !isinteger(x) || x < bounds[1] || x > bounds[2]
        result[i] = false
        break
      end
    end
  end
  result
end

function link(c::Chains)
  cc = copy(c.value.data)
  for j in 1:length(c.names)
    x = cc[:, j, :]
    if minimum(x) > 0.0
      cc[:, j, :] = maximum(x) < 1.0 ? logit.(x) : log.(x)
    end
  end
  cc
end

"""
    _trim_name_map(names::Vector, name_map::NamedTuple)

This is an internal function used to remove values from a name map
and return a new name_map.
"""
function _trim_name_map(names::Vector, name_map::NamedTuple)
    n = _namedtuple2dict(name_map)
    for (key, values) in n
        intersection = values ∩ names
        if length(intersection) > 0
            n[key] = intersection
        else
            delete!(n, key)
        end
    end
    return _dict2namedtuple(n)
end

### Chains specific functions ###
"""
    sort(c::Chains)

Returns a new column-sorted version of `c`, using natural sort order.
"""
function Base.sort(c::Chains)
    v = c.value
    x, y, z = size(v)
    unsorted = collect(zip(1:y, v.axes[2].val))
    sorted = sort(unsorted, by = x -> string(x[2]), lt=natural)
    new_axes = (v.axes[1], Axis{:var}([n for (_, n) in sorted]), v.axes[3])
    new_v = copy(v.data)
    for i in eachindex(sorted)
        new_v[:, i, :] = v[:, sorted[i][1], :]
    end

    # Sort the name map too:
    name_dict = Dict()
    for (name, values) in pairs(c.name_map)
        name_dict[name] = sort(values, by=x -> string(x), lt=natural)
    end
    new_name_map = _dict2namedtuple(name_dict)

    aa = AxisArray(new_v, new_axes...)
    return Chains{eltype(c.value),typeof(c.logevidence),typeof(new_name_map),
        typeof(c.info)}(aa, c.logevidence, new_name_map, c.info)
end

"""
    setinfo(c::Chains, n::NamedTuple)

Returns a new `Chains` object with a `NamedTuple` type `n` placed in the `info` field.

Example:

```julia
new_chn = setinfo(chn, NamedTuple{(:a, :b)}((1, 2)))
```
"""
function setinfo(c::Chains, n::NamedTuple)
    return Chains{eltype(c.value),typeof(c.logevidence),typeof(c.name_map),typeof(n)}(
        c.value, c.logevidence, c.name_map, n)
end

set_section(c::Chains, nt::NamedTuple) = set_section(c, _namedtuple2dict(nt))

"""
    set_section(c::Chains, nt::Dict)

Changes a chains name mapping to a provided dictionary. This also supports a NamedTuple.
Any parameters in the chain that are unassigned will be placed into
the :parameters section.
"""
function set_section(c::Chains, d::Dict)
    # Add :parameters if it's not there.
    if !(:parameters in keys(d))
        d[:parameters] = []
    end

    # Make sure all the names are in the new name map.
    nms = Set([])
    for values in values(d)
        for val in values
            push!(nms, val)
        end
    end
    missing_names = setdiff(names(c), nms)

    # Assign everything to :parameters if anything's missing.
    if length(missing_names) > 0
        @warn "Section mapping does not contain all parameter names, " *
            "$missing_names assigned to :parameters."
        push!(d[:parameters], missing_names...)
    end

    nt = _dict2namedtuple(d)
    return Chains{eltype(c.value),typeof(c.logevidence),typeof(nt),typeof(c.info)}(
        c.value, c.logevidence, nt, c.info)
end

function _use_showall(c::Chains, section::Symbol)
    if section == :parameters && !in(:parameters, keys(c.name_map))
        return true
    end
    return false
end

function _clean_sections(c::Chains, sections::Union{Vector{Symbol}, Symbol})
    sections = sections isa AbstractArray ? sections : [sections]
    ks = collect(keys(c.name_map))
    return ks ∩ sections
end

#################### Concatenation ####################

Base.cat(c::Chains, cs::Chains...; dims = Val(1)) = _cat(dims, c, cs...)
Base.cat(c::T, cs::T...; dims = Val(1)) where T<:Chains = _cat(dims, c, cs...)

Base.vcat(c::Chains, cs::Chains...) = _cat(Val(1), c, cs...)
Base.vcat(c::T, cs::T...) where T<:Chains = _cat(Val(1), c, cs...)

Base.hcat(c::Chains, cs::Chains...) = _cat(Val(2), c, cs...)
Base.hcat(c::T, cs::T...) where T<:Chains = _cat(Val(2), c, cs...)

AbstractMCMC.chainscat(c::Chains, cs::Chains...) = _cat(Val(3), c, cs...)

_cat(dim::Int, cs::Chains...) = _cat(Val(dim), cs...)

function _cat(::Val{1}, c1::Chains, args::Chains...)
    # check inputs
    thin = step(c1)
    all(c -> step(c) == thin, args) || throw(ArgumentError("chain thinning differs"))
    nms = names(c1)
    all(c -> names(c) == nms, args) || throw(ArgumentError("chain names differ"))
    chns = chains(c1)
    all(c -> chains(c) == chns, args) || throw(ArgumentError("sets of chains differ"))

    # concatenate all chains
    data = mapreduce(c -> c.value.data, vcat, args; init = c1.value.data)
    value = AxisArray(data;
                      iter = range(first(c1); length = size(data, 1), step = thin),
                      var = nms,
                      chain = chns)

    return Chains(value, missing, c1.name_map, c1.info)
end

function _cat(::Val{2}, c1::Chains, args::Chains...)
    # check inputs
    rng = range(c1)
    all(c -> range(c) == rng, args) || throw(ArgumentError("chain ranges differ"))
    chns = chains(c1)
    all(c -> chains(c) == chns, args) || throw(ArgumentError("sets of chains differ"))

    # combine names and sections of parameters
    nms = names(c1)
    n = length(nms)
    for c in args
        nms = union(nms, names(c))
        n += length(names(c))
        n == length(nms) || throw(ArgumentError("non-unique parameter names"))
    end

    name_map = mapreduce(c -> c.name_map, merge_union, args; init = c1.name_map)

    # concatenate all chains
    data = mapreduce(c -> c.value.data, hcat, args; init = c1.value.data)
    value = AxisArray(data; iter = rng, var = nms, chain = chns)

    return Chains(value, missing, name_map, c1.info)
end

function _cat(::Val{3}, c1::Chains, args::Chains...)
    # check inputs
    rng = range(c1)
    all(c -> range(c) == rng, args) || throw(ArgumentError("chain ranges differ"))
    nms = names(c1)
    all(c -> names(c) == nms, args) || throw(ArgumentError("chain names differ"))

    # concatenate all chains
    data = mapreduce(c -> c.value.data, (x, y) -> cat(x, y; dims = 3), args;
                     init = c1.value.data)
    value = AxisArray(data; iter = rng, var = nms, chain = 1:size(data, 3))

    return Chains(value, missing, c1.name_map, c1.info)
end

function pool_chain(c::Chains)
    data = c.value.data
    pool_data = reshape(permutedims(data, [1, 3, 2]), :, size(data, 2), 1)
    return Chains(pool_data, names(c), c.name_map; info=c.info)
end

function set_names(c::Chains, d::Dict; sorted::Bool=true)
    # Set new parameter names.
    params = names(c)
    new_params = replace(params, d...)

    # Iterate through each pair in the name map to
    # create a new one.
    new_map = Dict()
    for (section, parameters) in pairs(c.name_map)
        new_map[section] = replace(parameters, d...)
    end

    # Return a new chains object.
    return Chains(
        c.value.data,
        new_params,
        new_map;
        info=c.info,
        evidence=c.logevidence,
        sorted=sorted,
    )
end


text:
"""
    struct JohnsonState{T, U}

An [`AbstractPathState`](@ref) designed for Johnson shortest-paths calculations.
"""
struct JohnsonState{T <: Real,U <: Integer} <: AbstractPathState
    dists::Matrix{T}
    parents::Matrix{U}
end

"""
    johnson_shortest_paths(g, distmx=weights(g))

Use the [Johnson algorithm](https://en.wikipedia.org/wiki/Johnson%27s_algorithm)
to compute the shortest paths between all pairs of vertices in graph `g` using an
optional distance matrix `distmx`.

### Implementation Notes
If the parameter parallel is set true, dijkstra_shortest_paths will run in parallel.
Parallel bellman_ford_shortest_paths is currently unavailable
Return a [`LightGraphs.JohnsonState`](@ref) with relevant
traversal information.
Behaviour in case of negative cycle depends on bellman_ford_shortest_paths.
Throws NegativeCycleError() if a negative cycle is present.

### Performance
Complexity: O(|V|*|E|)
If distmx is not mutable or of type, DefaultDistance than a sparse matrix will be produced using distmx.
In the case that distmx is immutable, to reduce memory overhead,
if edge (a, b) does not exist in g then distmx[a, b] should be set to 0.
"""
function johnson_shortest_paths(g::AbstractGraph{U},
    distmx::AbstractMatrix{T}=weights(g)) where T <: Real where U <: Integer

    nvg = nv(g)
    type_distmx = typeof(distmx)
    #Change when parallel implementation of Bellman Ford available
    wt_transform = bellman_ford_shortest_paths(g, vertices(g), distmx).dists

    if !type_distmx.mutable && type_distmx !=  LightGraphs.DefaultDistance
        distmx = sparse(distmx) #Change reference, not value
    end

    #Weight transform not needed if all weights are positive.
    if type_distmx !=  LightGraphs.DefaultDistance
        for e in edges(g)
            distmx[src(e), dst(e)] += wt_transform[src(e)] - wt_transform[dst(e)]
        end
    end

    dists = Matrix{T}(undef, nvg, nvg)
    parents = Matrix{U}(undef, nvg, nvg)
    for v in vertices(g)
        dijk_state = dijkstra_shortest_paths(g, v, distmx)
        dists[v, :] = dijk_state.dists
        parents[v, :] = dijk_state.parents
    end

    broadcast!(-, dists, dists, wt_transform)
    for v in vertices(g)
        dists[:, v] .+= wt_transform[v] #Vertical traversal prefered
    end

    if type_distmx.mutable
        for e in edges(g)
            distmx[src(e), dst(e)] += wt_transform[dst(e)] - wt_transform[src(e)]
        end
    end

    return JohnsonState(dists, parents)
end

function enumerate_paths(s::JohnsonState{T,U}, v::Integer) where T <: Real where U <: Integer
    pathinfo = s.parents[v, :]
    paths = Vector{Vector{U}}()
    for i in 1:length(pathinfo)
        if (i == v) || (s.dists[v, i] == typemax(T))
            push!(paths, Vector{U}())
        else
            path = Vector{U}()
            currpathindex = i
            while currpathindex != 0
                push!(path, currpathindex)
                currpathindex = pathinfo[currpathindex]
            end
            push!(paths, reverse(path))
        end
    end
    return paths
end

enumerate_paths(s::JohnsonState) = [enumerate_paths(s, v) for v in 1:size(s.parents, 1)]
enumerate_paths(st::JohnsonState, s::Integer, d::Integer) = enumerate_paths(st, s)[d]


text:
using MinimaxEstimation
using Documenter

makedocs(;
    modules=[MinimaxEstimation],
    authors="Olle Kjellqvist",
    repo="https://github.com/kjellqvist/MinimaxEstimation.jl/blob/{commit}{path}#L{line}",
    sitename="MinimaxEstimation.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://kjellqvist.github.io/MinimaxEstimation.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => Any[
            "Library" => "lib.md",
            "Example" => "example.md"
        ]
    ],
)

deploydocs(;
    devbranch = "main",
    repo="github.com/kjellqvist/MinimaxEstimation.jl",
)


text:
module HaskellGC

export retain, release

refs = Any[]
stack = Int64[]

# retain a reference
function retain(r)
  if length(stack) == 0
    push!(refs, r)
    length(refs) :: Int64
  else
    ind = pop!(stack)
    refs[ind] = r
    ind :: Int64
  end
end

# release a held reference
function release(i::Int64)
  refs[i] = Union{}
  push!(stack, i)
end

# utility to show a julia value in haskell
function show(x)
  Base.unsafe_convert(Cstring, repr(x))
end

function finalize_hs(p::Ptr{Void})
  ccall(:hs_free_stable_ptr, Void, (Ptr{Void},), p)
end

end

max_stars_repo_name: tkonolige/inline-julia

max_issues_repo_name: tkonolige/inline-julia

max_forks_repo_name: tkonolige/inline-julia


text:
module BasicTypes
  include(joinpath(@__DIR__, "testcommon.jl"))

  struct ImmutableBits
    a::Float64
    b::Float64
  end

  mutable struct MutableBits
    a::Float64
    b::Float64
  end

  struct A
    x :: Float32
    y :: Float32
  end

  @wrapmodule CxxWrap.CxxWrapCore.libbasic_types

  function __init__()
    @initcxx
  end
end

using CxxWrap
using Test

result = Ref(0.0)
function test_boxed_struct(x)
    result[] = x.a + x.b
    return
end

@testset "$(basename(@__FILE__)[1:end-3])" begin

let a = BasicTypes.A(2,3)
  @test BasicTypes.f(a) == 5.0
  @test BasicTypes.g(a) == 5.0
  @test BasicTypes.g(Ref(a)) == 5.0
  @test BasicTypes.h(Ref(a)) == 5.0
  @test BasicTypes.h(C_NULL) == 0.0
end

let imm = BasicTypes.ImmutableBits(1.0, 2.0)
  @test BasicTypes.increment_immutable(Ref(imm)) == BasicTypes.ImmutableBits(2.0, 3.0)
end

let f = Float32(5.0), a = [f]
  @test BasicTypes.twice_val(f) == 10.0
  @test BasicTypes.twice_cref(f) == 10.0
  @test BasicTypes.twice_ref(Ref(f)) == 10.0
  @test BasicTypes.twice_cptr(Ref(f)) == 10.0
  @test BasicTypes.twice_ptr(Ref(f)) == 10.0
  @test BasicTypes.twice_ptr(pointer(a)) == 10.0
  @test BasicTypes.twice_ptr(a) == 10.0
  BasicTypes.twice_ptr_mut(a)
  @test a[1] == 10.0
end

let f = Ref(Float32(2.0))
  BasicTypes.twice_ref_mut(f)
  @test f[] == 4.0
end

let a = Ref(BasicTypes.A(1,2))
  b = BasicTypes.return_a_ptr(a)
  @test a[].x == 2
  @test a[].y == 3
  @test unsafe_load(b) == a[]

  b = BasicTypes.return_a_cptr(a)
  @test a[].x == 3
  @test a[].y == 4
  @test unsafe_load(b) == a[]

  b = BasicTypes.return_a_ref(a)
  @test a[].x == 4
  @test a[].y == 5
  @test b[] == a[]

  b = BasicTypes.return_a_cref(a)
  @test a[].x == 5
  @test a[].y == 6
  @test b[] == a[]
end

let s = "abc"
  @test BasicTypes.strlen_cchar(s) == length(s)
  @test BasicTypes.strlen_char(s) == length(s)
end

let s = StdString("hello")
  @test BasicTypes.strlen_str(s) == 5
  @test BasicTypes.strlen_strcref(s) == 5
  @test BasicTypes.strlen_strref(s) == 5
  @test BasicTypes.strlen_strptr(CxxPtr(s)) == 5
  @test BasicTypes.strlen_strcptr(ConstCxxPtr(s)) == 5
end

let s = BasicTypes.StringHolder("hello")
  # Check reference return values
  @test BasicTypes.str_return_val(s) == "hello"
  @test BasicTypes.str_return_cref(s)[] == "hello"
  strref = BasicTypes.str_return_ref(s)
  @test strref[] == "hello"
  strptr = BasicTypes.str_return_ptr(s)
  @test strptr[] == "hello"
  @test BasicTypes.str_return_cptr(s)[] == "hello"

  # Modification through reference
  BasicTypes.replace_str_val!(strref, "world")
  @test strref[] == "world"
  @test strptr[] == "world"
  @test BasicTypes.str_return_val(s) == "world"

  # Modification through pointer
  BasicTypes.replace_str_val!(CxxRef(strptr), "bye!")
  @test strref[] == "bye!"
  @test strptr[] == "bye!"
  @test BasicTypes.str_return_val(s) == "bye!"

  # Modification through value
  strval = BasicTypes.str_return_val(s)
  BasicTypes.replace_str_val!(strval, "no really, bye!")
  @test strval == "no really, bye!"
  @test strptr[] == "bye!"

  # Check that const is respected
  strcref = BasicTypes.str_return_cref(s)
  strcptr = BasicTypes.str_return_cptr(s)
  @test_throws MethodError BasicTypes.replace_str_val!(strcref, "can't work")
  @test_throws MethodError BasicTypes.replace_str_val!(strcptr, "can't work")
end

let cfunc = @safe_cfunction(test_boxed_struct, Cvoid, (Any,))
  BasicTypes.boxed_mirrored_type(cfunc)
  @test result[] == 3.0

  BasicTypes.boxed_mutable_mirrored_type(cfunc)
  @test result[] == 5.0
end

let a = CxxChar(3), b = CxxWchar(2)
  @test typeof(a*b) == CxxWrap.CxxWrapCore.julia_int_type(CxxWchar)
  @test typeof(a+b) == CxxWrap.CxxWrapCore.julia_int_type(CxxWchar)
  @test typeof(a/b) == Float64
  @test a*b == 6
  @test a+b == 5
  @test typeof(a/a) == Float64
  @test typeof(a+a) == CxxWrap.CxxWrapCore.julia_int_type(CxxChar)
  @test typeof(a*a) == CxxWrap.CxxWrapCore.julia_int_type(CxxChar)
end

let buf = IOBuffer()
  show(buf,CxxBool(true))
  @test String(take!(buf)) == "true"
  show(buf,CxxLong(42))
  @test String(take!(buf)) == "42"
end

@test BasicTypes.test_for_each_type() == (sizeof(Float32) + sizeof(Float64))

@test BasicTypes.strict_method(CxxChar(1)) == "char"
@test BasicTypes.strict_method(CxxLong(1)) == "long"
@test_throws MethodError BasicTypes.strict_method(Int16(1))
@test BasicTypes.strict_method(CxxBool(true)) == "bool"
@test BasicTypes.strict_method(true) == "bool"
@test BasicTypes.loose_method(Int32(3)) == "int"
@test BasicTypes.loose_method(false) == "bool"

let (intnames, inttypes) = BasicTypes.julia_integer_mapping()
  lmax = maximum(length.(intnames))
  for (name, type) in zip(intnames, inttypes)
    padding = repeat(' ', lmax-length(name))
    println("$name$padding -> $type")
  end
end

end


text:
@defstruct ReplicationLayer Layer (
    name :: AbstractString = "power",
    (num_copies :: Int = 1, num_copies > 0),
    (tops :: Vector{Symbol} = Symbol[], length(tops) == 1),
    (bottoms :: Vector{Symbol} = Symbol[], length(bottoms) == length(tops)),
)
@characterize_layer(ReplicationLayer,
    can_do_bp => true
)

type ReplicationLayerState <: LayerState
    layer            :: ReplicationLayer
    blobs            :: Vector{Blob}
    blobs_diff :: Vector{Blob}
end

function setup(backend::Backend, layer::ReplicationLayer, inputs::Vector{Blob},
                diffs::Vector{Blob})
    data_type = eltype(inputs[1])
    output_dim = (layer.num_copies, size(inputs[1])...)
    blobs = [make_blob(backend, data_type, output_dim)]

    if all(map(b -> isa(b, NullBlob), diffs))
        blobs_diff = Blob[NullBlob()]
    else
        blobs_diff = [make_blob(backend, data_type, output_dim)]
    end

    state = ReplicationLayerState(layer, blobs, blobs_diff)
end

function shutdown(backend::Backend, state::ReplicationLayerState)
    map(destroy, state.blobs)
    map(destroy, state.blobs_diff)
end

function forward(backend::CPUBackend, state::ReplicationLayerState,
            inputs::Vector{Blob})
    input = inputs[1]
    output = state.blobs[1]

    p, n = size(input)
    for i in 1:n, j in 1:p
        output.data[:, j, i] = input.data[j, i]
    end
end

function backward(backend::CPUBackend, state::ReplicationLayerState,
        inputs::Vector{Blob}, diffs::Vector{Blob})

    if !isa(diffs[1], NullBlob)
        data_type = eltype(inputs[1])
        p, n = size(inputs[1])

        for i in 1:n, j in 1:p
            diffs[1].data[j, i] = sum(state.blobs_diff[1].data[:, j, i])
        end
    end
end


text:
module Tyndall
using Jultran
using Sundials
using JSON
export planck, dry_adiabat, isothermal, OLR, compute_κ_grid, dry_adiabat_with_strat

# constants
p_s = 1.013e5 # p_s == surface pressure; ps == pressures
R = 287
c_p = 1006
σ = 5.67e-8
h = 6.626075540e-34    #Planck's constant
c = 2.99792458e8       #Speed of light
k =1.38065812e-23      #Boltzman thermodynamic constant
g = 9.8
cos_θ = .5

# settings

ν_min = 10
ν_max = 2500
no_of_ν_bands = 80
ν_step = (ν_max - ν_min) / no_of_ν_bands
νs = ν_min:ν_step:ν_max
ν_bands = [ν:(ν+ν_step) for ν in νs[1:end-1]]
no_of_γs = 15
γs = 1:no_of_γs

no_of_ps = 4
p_step = 1.013e5/(no_of_ps - 1)
ps = 0:p_step:1.013e5

T_min = 150
T_max = 340
no_of_Ts = 4
T_step = (T_max-T_min)/(no_of_Ts-1)
Ts = 150:T_step:340

κ_grid = reshape(JSON.parse(open("$(Pkg.dir())/Tyndall/data/κ_grid.json")), (no_of_ps, no_of_Ts, no_of_ν_bands, no_of_γs))

# ν in cm-1
planck(ν,T) = 100 * 2h * c^2 * (100ν)^3 / expm1(c*h*100ν/(k*T))

I_plus(ν_i, γ, T, q) = function(p, I_plus, I_plus_dot)
    ν = mean(ν_bands[ν_i])
    I_plus_dot[1] = q * κ(ν_i, γ, p, T(p)) * (I_plus[1] - π * planck(ν,T(p))) / (g * cos_θ)
end

function κ(ν_i, γ, p, T)
    if p >= 0
        p_weight = (p % p_step) / p_step
        p_i = fld(p,p_step) + 1
        if p_weight < 1.e-8
            p_weight = 1
            p_i -= 1
        end

        T_weight = (T % T_step) / T_step
        T_i = fld(T-T_min,T_step) + 1
        sum(T_weight * (κ_grid[p_i:p_i+1,T_i,ν_i,γ] .* [p_weight,1-p_weight]) + (1 - T_weight) * (κ_grid[p_i:p_i+1,T_i+1,ν_i,γ] .* [p_weight,1-p_weight]))
    else
        0
    end
end

function OLR(T,q; by_ν=false)
    if by_ν
        map(1:no_of_ν_bands) do ν_i
            sum(map(1:no_of_γs) do γ
                Sundials.ode(I_plus(ν_i,γ,T,q), [π * planck(mean(ν_bands[ν_i]),T(p_s))], [p_s, 0])[2]
            end)/no_of_γs
        end
    else
        sum(OLR(T,q,by_ν=true)) * ν_step
    end
end

# temperature profiles
isothermal(T_s) = p -> T_s
# dry adiabat
dry_adiabat(T_s) = p -> T_s * (p/p_s)^(R/c_p)
# dry_adiabat_with_strat(T_s) = p -> p > 1.e4 ? T_s * (p/p_s)^(R/c_p) : T_s * (1.e4/p_s)^(R/c_p)
dry_adiabat_with_strat(T_s) = p -> p > p_s * (160.0/T_s) ^(c_p/R) ? T_s * (p/p_s)^(R/c_p) : 160.0

function compute_κ_grid()
    κ_grid = Array(Float64, (no_of_ps, no_of_Ts, no_of_ν_bands, no_of_γs))

    for (p_i,p)=enumerate(ps), (T_i,T)=enumerate(Ts)
        println("p is $p")
        println("T is $T")
        ν_i = 0
        for ν_band=ν_bands
            println("ν_i is $(ν_i)")
            ν_i += 1
            νs_in_band = ν_band[1]:.01:ν_band[end]
            sorted_κs = sort(compute_κs(:co2, νs_in_band, T, p)) # compute_κs from Jultran
            γ_size = length(νs_in_band) / no_of_γs
            for γ=γs
                κ_grid[p_i,T_i,ν_i,γ] = mean(filter(x->!isnan(x),sorted_κs[(1:γ_size) + floor(γ_size * (γ-1))]))
            end
        end
    end

    open("$(Pkg.dir())/Tyndall/data/κ_grid.json", "w+") do f
        JSON.print(f, reshape(κ_grid,(no_of_ps * no_of_Ts * no_of_ν_bands * no_of_γs)))
    end

    κ_grid
end

end # module


text:
#=
Basic Types and Aliases
=#

const Mb = 1_000_000
const GENOMESIZE = 3235Mb

# ALIAS NEW TYPES FOR INCREASED CODE READABILITY
if GENOMESIZE < typemax(UInt32)
   const CoordInt = UInt32
else
   const CoordInt = UInt64
end

const CoordTuple = Tuple{Vararg{CoordInt}}
const CoordArray = Vector{CoordInt}
const CoordTree  = IntervalTree{CoordInt,IntervalTrees.Interval{CoordInt}}
const CoordSet   = Set{Tuple{CoordInt,CoordInt}}

const ExonInt    = UInt16
const ExonTree   = IntervalTree{ExonInt,IntervalTrees.Interval{ExonInt}}

const SGSequence = BioSequences.BioSequence{DNAAlphabet{4}}

const GeneName    = String
const SeqName     = String
const RefSeqId    = String
const ASCIIString = String

const GeneMeta = Tuple{GeneName, SeqName, Char}

const BufOut = BufferedStreams.BufferedOutputStream

struct GeneInfo
   gene::GeneName
   name::SeqName
   strand::Bool # pos is true, neg is false
end

GeneInfo( gene::S, name::S, strand::Char ) where {S <: AbstractString} =
                               GeneInfo( convert(GeneName, gene),
                                         convert(SeqName, name),
                                         strand == '+' ? true :
                                                         false )

struct TxInfo
   name::RefSeqId
   txstart::CoordInt
   txend::CoordInt
   exnum::CoordInt
end

#=
Basic IO Functions
=#

fixpath( str::String ) = abspath( expanduser( str ) )

isgzipped( filename::String ) = splitext(filename)[2] == ".gz"

function hasextension( filename::String, ext::String )
   restr = "\\.$ext\$"
   return occursin(restr, filename)
end

function increment!( dict::Dict{K,V}, key::K, val::V=one(V) ) where {K,V}
   if haskey( dict, key )
      dict[key] += val
   else
      dict[key] = val
   end
   dict
end


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

struct InvalidCharError <: Exception
    char::Char
end
struct CodePointError <: Exception
    code::Integer
end
@noinline invalid_char(c::Char) = throw(InvalidCharError(c))
@noinline code_point_err(u::UInt32) = throw(CodePointError(u))

function ismalformed(c::Char)
    u = reinterpret(UInt32, c)
    l1 = leading_ones(u) << 3
    t0 = trailing_zeros(u) & 56
    (l1 == 8) | (l1 + t0 > 32) |
    (((u & 0x00c0c0c0) ⊻ 0x00808080) >> t0 != 0)
end

@inline is_overlong_enc(u::UInt32) = (u >> 24 == 0xc0) | (u >> 24 == 0xc1) | (u >> 21 == 0x0704) | (u >> 20 == 0x0f08)

function isoverlong(c::Char)
    u = reinterpret(UInt32, c)
    is_overlong_enc(u)
end

function UInt32(c::Char)
    # TODO: use optimized inline LLVM
    u = reinterpret(UInt32, c)
    u < 0x80000000 && return u >> 24
    l1 = leading_ones(u)
    t0 = trailing_zeros(u) & 56
    (l1 == 1) | (8l1 + t0 > 32) |
    ((((u & 0x00c0c0c0) ⊻ 0x00808080) >> t0 != 0) | is_overlong_enc(u)) &&
        invalid_char(c)::Union{}
    u &= 0xffffffff >> l1
    u >>= t0
    (u & 0x0000007f >> 0) | (u & 0x00007f00 >> 2) |
    (u & 0x007f0000 >> 4) | (u & 0x7f000000 >> 6)
end

function decode_overlong(c::Char)
    u = reinterpret(UInt32, c)
    l1 = leading_ones(u)
    t0 = trailing_zeros(u) & 56
    u &= 0xffffffff >> l1
    u >>= t0
    (u & 0x0000007f >> 0) | (u & 0x00007f00 >> 2) |
    (u & 0x007f0000 >> 4) | (u & 0x7f000000 >> 6)
end

function Char(u::UInt32)
    u < 0x80 && return reinterpret(Char, u << 24)
    u < 0x00200000 || code_point_err(u)::Union{}
    c = ((u << 0) & 0x0000003f) | ((u << 2) & 0x00003f00) |
        ((u << 4) & 0x003f0000) | ((u << 6) & 0x3f000000)
    c = u < 0x00000800 ? (c << 16) | 0xc0800000 :
        u < 0x00010000 ? (c << 08) | 0xe0808000 :
                         (c << 00) | 0xf0808080
    reinterpret(Char, c)
end

function (T::Union{Type{Int8},Type{UInt8}})(c::Char)
    i = reinterpret(Int32, c)
    i ≥ 0 ? ((i >>> 24) % T) : T(UInt32(c))
end

function Char(b::Union{Int8,UInt8})
    0 ≤ b ≤ 0x7f ? reinterpret(Char, (b % UInt32) << 24) : Char(UInt32(b))
end

convert(::Type{Char}, x::Number) = Char(x)
convert(::Type{T}, x::Char) where {T<:Number} = T(x)

rem(x::Char, ::Type{T}) where {T<:Number} = rem(UInt32(x), T)

typemax(::Type{Char}) = reinterpret(Char, typemax(UInt32))
typemin(::Type{Char}) = reinterpret(Char, typemin(UInt32))

size(c::Char) = ()
size(c::Char,d) = convert(Int, d) < 1 ? throw(BoundsError()) : 1
ndims(c::Char) = 0
ndims(::Type{Char}) = 0
length(c::Char) = 1
firstindex(c::Char) = 1
lastindex(c::Char) = 1
getindex(c::Char) = c
getindex(c::Char, i::Integer) = i == 1 ? c : throw(BoundsError())
getindex(c::Char, I::Integer...) = all(x -> x == 1, I) ? c : throw(BoundsError())
first(c::Char) = c
last(c::Char) = c
eltype(::Type{Char}) = Char

start(c::Char) = false
next(c::Char, state) = (c, true)
done(c::Char, state) = state
isempty(c::Char) = false
in(x::Char, y::Char) = x == y

==(x::Char, y::Char) = reinterpret(UInt32, x) == reinterpret(UInt32, y)
isless(x::Char, y::Char) = reinterpret(UInt32, x) < reinterpret(UInt32, y)
hash(x::Char, h::UInt) =
    hash_uint64(((reinterpret(UInt32, x) + UInt64(0xd4d64234)) << 32) ⊻ UInt64(h))
widen(::Type{Char}) = Char

-(x::Char, y::Char) = Int(x) - Int(y)
-(x::Char, y::Integer) = Char(Int32(x) - Int32(y))
+(x::Char, y::Integer) = Char(Int32(x) + Int32(y))
+(x::Integer, y::Char) = y + x

print(io::IO, c::Char) = (write(io, c); nothing)

const hex_chars = UInt8['0':'9';'a':'z']

function show(io::IO, c::Char)
    if c <= '\\'
        b = c == '\0' ? 0x30 :
            c == '\a' ? 0x61 :
            c == '\b' ? 0x62 :
            c == '\t' ? 0x74 :
            c == '\n' ? 0x6e :
            c == '\v' ? 0x76 :
            c == '\f' ? 0x66 :
            c == '\r' ? 0x72 :
            c == '\e' ? 0x65 :
            c == '\'' ? 0x27 :
            c == '\\' ? 0x5c : 0xff
        if b != 0xff
            write(io, 0x27, 0x5c, b, 0x27)
            return
        end
    end
    if isoverlong(c) || ismalformed(c)
        write(io, 0x27)
        u = reinterpret(UInt32, c)
        while true
            a = hex_chars[((u >> 28) & 0xf) + 1]
            b = hex_chars[((u >> 24) & 0xf) + 1]
            write(io, 0x5c, 'x', a, b)
            (u <<= 8) == 0 && break
        end
        write(io, 0x27)
    elseif isprint(c)
        write(io, 0x27, c, 0x27)
    else # unprintable, well-formed, non-overlong Unicode
        u = UInt32(c)
        write(io, 0x27, 0x5c, c <= '\x7f' ? 0x78 : c <= '\uffff' ? 0x75 : 0x55)
        d = max(2, 8 - (leading_zeros(u) >> 2))
        while 0 < d
            write(io, hex_chars[((u >> ((d -= 1) << 2)) & 0xf) + 1])
        end
        write(io, 0x27)
    end
    return
end

function show(io::IO, ::MIME"text/plain", c::Char)
    show(io, c)
    if !ismalformed(c)
        print(io, ": ")
        if isoverlong(c)
            print(io, "[overlong] ")
            u = decode_overlong(c)
            c = Char(u)
        else
            u = UInt32(c)
        end
        h = hex(u, u ≤ 0xffff ? 4 : 6)
        print(io, (isascii(c) ? "ASCII/" : ""), "Unicode U+", h)
    else
        print(io, ": Malformed UTF-8")
    end
    abr = Unicode.category_abbrev(c)
    str = Unicode.category_string(c)
    print(io, " (category ", abr, ": ", str, ")")
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
using Test
using MLStyle, MacroTools

@testset "search_and_replace" begin

ex = @search_and_replace sin(x) sin(x) cos(x)
@test ex == :( cos(x) )

ex = @search_and_replace sin(y) sin($x) cos($x)
@test ex == :( cos(y) )

ex = @search_and_replace sin(sin(x)) sin($x) cos($x +1)
@test ex == :( cos(cos(x+1) +1) )

end


text:
using ImageShow, ImageCore
using TestImages, FileIO
using Test

function check_summary(n, msg)
    function generate_summary_regex(n)
        summary_regex = raw"""Frame: \d+/\d+ FPS: \d+\.\d+\s+\nexit: "q" play\/pause: "space-bar" seek: "arrow keys"\n"""
        Regex(mapreduce(i->summary_regex, (x,y)->x*raw".*"*y, 1:n))
    end
    function _check_summary(n, msg)
        return match(generate_summary_regex(n), msg) !== nothing
    end

    return _check_summary(n, msg) && !_check_summary(n+1, msg)
end

@testset "multipage" begin
    @testset "play" begin
        img = RGB.(testimage("mri-stack"))
        framestack = [img[:, :, i] for i in 1:size(img, 3)]

        mktempdir() do workdir
            tmpfile = joinpath(workdir, "multipage.png")
            save(tmpfile, framestack[1])
            frame_size = stat(tmpfile).size

            # Case: quit immediately -- render once
            summary_output_io = IOBuffer()
            key_input_io = IOBuffer(UInt8['q'])
            open(tmpfile, "w") do fn
                ImageShow._play(framestack; fps=1, paused=true, quit_after_play=true, display_io=fn, summary_io=summary_output_io, keyboard_io=key_input_io)
            end
            summary_msg = String(take!(summary_output_io))
            @test check_summary(1, summary_msg)
            @test RGB.(load(tmpfile)) == framestack[1]
            @test 1.0 == stat(tmpfile).size/frame_size

            # Case: quit after one forward -- render twice
            fn = open(tmpfile, "w")
            summary_output_io = IOBuffer()
            key_input_io = IOBuffer(UInt8['f', 'q'])
            open(tmpfile, "w") do fn
                ImageShow._play(framestack; fps=5, paused=true, quit_after_play=true, display_io=fn, summary_io=summary_output_io, keyboard_io=key_input_io)
            end
            summary_msg = String(take!(summary_output_io))
            @test check_summary(2, summary_msg)
            # NOTE: show method will append to the given file handler, however, PNG reader will only
            #        read the first valid png data; all extra data block are discarded.
            # This somehow still proves that we're writing more than one image to the display_io
            @test RGB.(load(tmpfile)) == framestack[1]
            @test 1.8 < stat(tmpfile).size/frame_size <= 2.2

            # Case: quit after resume play -- render twice
            summary_output_io = IOBuffer()
            key_input_io = IOBuffer(UInt8[' ', 'q'])
            open(tmpfile, "w") do fn
                ImageShow._play(framestack; fps=5, paused=true, quit_after_play=true, display_io=fn, summary_io=summary_output_io, keyboard_io=key_input_io)
            end
            summary_msg = String(take!(summary_output_io))
            @test check_summary(2, summary_msg)
            # NOTE: show method will append to the given file handler, however, PNG reader will only
            #        read the first valid png data; all extra data block are discarded.
            # This somehow still proves that we're writing more than one image to the display_io
            @test RGB.(load(tmpfile)) == framestack[1]
            @test 1.8 < stat(tmpfile).size/frame_size <= 2.2

            # Case: quit after a useless control -- only render once
            summary_output_io = IOBuffer()
            key_input_io = IOBuffer(UInt8['b', '?', 'q'])
            open(tmpfile, "w") do fn
                ImageShow._play(framestack; fps=5, paused=true, quit_after_play=true, display_io=fn, summary_io=summary_output_io, keyboard_io=key_input_io)
            end
            summary_msg = String(take!(summary_output_io))
            @test check_summary(1, summary_msg)
            @test RGB.(load(tmpfile)) == framestack[1]
            @test 1.0 == stat(tmpfile).size/frame_size

            # Case: quit after all play
            summary_output_io = IOBuffer()
            key_input_io = IOBuffer()
            # use a small fps to make sure each frame are actually written
            open(tmpfile, "w") do fn
                ImageShow._play(framestack; fps=15, paused=false, quit_after_play=true, display_io=fn, summary_io=summary_output_io, keyboard_io=key_input_io)
            end
            summary_msg = String(take!(summary_output_io))
            @test check_summary(length(framestack), summary_msg)
            # NOTE: show method will append to the given file handler, however, PNG reader will only
            #        read the first valid png data; all extra data block are discarded.
            # This somehow still proves that we're writing more than one image to the display_io
            @test RGB.(load(tmpfile)) == framestack[1]
            @test 0.8*length(framestack) < stat(tmpfile).size/frame_size < length(framestack)
        end
    end

    @testset "utils" begin
        # Although it's a no-op, it gets blocked at fps=2, which is about 0.5 second
        f() = nothing
        f() # precompile it
        t = @elapsed ImageShow.fixed_fps(f, 0.5)
        @test 1.5 < t < 2.5
    end
end


text:
filename = Pkg.dir("DataFrames", "test", "data", "big_data.csv")

minibatch_sizes = [1, 5, 25, 100, 1_000, 10_000]

for f in (colmeans, colvars, cor)
	for minibatch_size in minibatch_sizes
		ds = DataStream(filename, minibatch_size)
		N = 3
		df = benchmark(() -> apply(f, (ds,)),
			           "DataStream Functions",
			           join({
			           	      string(f),
			           	      "w/ minibatches of",
			           	      minibatch_size,
			           	      "rows"
			           	    }, " "),
			           N)
		# TODO: Keep permanent record
		printtable(df, header=false)
	end
end


text:
using Flux
using Flux: logitcrossentropy, onehotbatch, onecold, @epochs
using Flux.Data: DataLoader
using MLDatasets: MNIST
using Random
using Zygote

function get_dataloaders(batch_size::Int, shuffle::Bool)
    train_x, train_y = MNIST.traindata(Float32)
    test_x, test_y = MNIST.testdata(Float32)

    train_y, test_y = onehotbatch(train_y, 0:9), onehotbatch(test_y, 0:9)

    train_loader = DataLoader(train_x, train_y, batchsize=batch_size, shuffle=shuffle)
    test_loader = DataLoader(test_x, test_y, batchsize=batch_size, shuffle=shuffle)

    return train_loader, test_loader
end

function accuracy(data_loader, model)
    acc_correct = 0
    for (x_batch, y_batch) in data_loader
        batch_size = size(x_batch)[end]
        acc_correct += sum(onecold(model(x_batch)) .== onecold(y_batch)) / batch_size
    end
    return acc_correct / length(data_loader)
end

struct FFNetwork
    fc_1
    dropout
    fc_2
    FFNetwork(
        input_dim::Int, hidden_dim::Int, dropout_ratio::Float32, num_classes::Int
    ) = new(
        Dense(input_dim, hidden_dim, relu),
        Dropout(dropout_ratio),
        Dense(hidden_dim, num_classes),
    )
end

function (net::FFNetwork)(x)
    x = Flux.flatten(x)
    return net.fc_2(net.dropout(net.fc_1(x)))
end

function cross_entropy_loss(model, x, y)
    ŷ = model(x)
    return logitcrossentropy(model(x), y)
end

function main(num_epochs, batch_size, shuffle, η)
    train_loader, test_loader = get_dataloaders(batch_size, shuffle)

    model = FFNetwork(28*28, 128, 0.2f0, 10)
    trainable_params = Flux.params(model.fc_1, model.fc_2)
    optimiser = ADAM(η)

    for epoch = 1:num_epochs
        acc_loss = 0.0
        for (x_batch, y_batch) in train_loader
            loss, back = pullback(trainable_params) do
                cross_entropy_loss(model, x_batch, y_batch)
            end
            # Feed the pullback 1 to obtain the gradients and update then model parameters
            gradients = back(1f0)
            Flux.Optimise.update!(optimiser, trainable_params, gradients)
            acc_loss += loss
        end
        avg_loss = acc_loss / length(train_loader)
        println("Epoch: $epoch - loss: $avg_loss")
    end

    testmode!(model)
    @show accuracy(train_loader, model)
    @show accuracy(test_loader, model)
    println("Complete!")
end

if abspath(PROGRAM_FILE) == @__FILE__
    batch_size = 64
    shuffle_data = true
    η = 0.0001
    num_epochs = 5
    main(num_epochs, batch_size, shuffle_data, η)
end

max_stars_repo_name: alecokas/ml-framework-comparison

max_issues_repo_name: alecokas/ml-framework-comparison

max_forks_repo_name: alecokas/ml-framework-comparison


text:
using PackageCompiler

if length(ARGS) == 1
    sysimage_path = ARGS[1]
else
    sysimage_path = "sys_optmap.so"
end

println("Saving image as $sysimage_path")

#Add to this list other solver packages as needed
packages_to_compile = [:MultivariateStats, :Optim, :JuMP, :Clp, :Cbc, :GLPK, :Ipopt, :PyCall, :ColorSchemes, :ColorTypes, :Images, :Luxor, :CSV, :DataFrames, :Glob]

create_sysimage(packages_to_compile,
                sysimage_path=sysimage_path,
                precompile_execution_file="test_optmap.jl")

max_stars_repo_name: nmaterise/optmap

max_issues_repo_name: nmaterise/optmap

max_forks_repo_name: nmaterise/optmap