text:
module MDPs

using GridInterpolations
using SparseArrays
using Printf
using DelimitedFiles
using PyCall
using TickTock
export MDP, Policy, policy!, value_iteration, get_belief, dimensions, RectangleGrid, __init__

function __init__()
    py"""
    import numpy as np
    from scipy.stats import bernoulli
    from fancyimpute import SoftImpute, BiScaler, NuclearNormMinimization

    def generate_mask(mask_prob=0.5):
        # generate mask for limited observation
        # change NUM_STATES for different planning tasks
        NUM_STATES = 10000
        NUM_ACTIONS = 1000
        RESIZE_SHAPE = (NUM_STATES, NUM_ACTIONS)
        mask = bernoulli.rvs(p=mask_prob, size=(RESIZE_SHAPE[0], RESIZE_SHAPE[1])).astype(float)
        mask[mask < 1] = np.nan

        return mask

    def bregman(image, mask, weight, eps=1e-3, max_iter=100):
        rows, cols, dims = image.shape
        rows2 = rows + 2
        cols2 = cols + 2
        total = rows * cols * dims
        shape_ext = (rows2, cols2, dims)

        u = np.zeros(shape_ext)
        dx = np.zeros(shape_ext)
        dy = np.zeros(shape_ext)
        bx = np.zeros(shape_ext)
        by = np.zeros(shape_ext)

        u[1:-1, 1:-1] = image
        u[0, 1:-1] = image[1, :]
        u[1:-1, 0] = image[:, 1]
        u[-1, 1:-1] = image[-2, :]
        u[1:-1, -1] = image[:, -2]

        i = 0
        rmse = np.inf
        lam = 2 * weight
        norm = (weight + 4 * lam)

        while i < max_iter and rmse > eps:
            rmse = 0
            for k in range(dims):
                for r in range(1, rows + 1):
                    for c in range(1, cols + 1):
                        uprev = u[r, c, k]
                        ux = u[r, c + 1, k] - uprev
                        uy = u[r + 1, c, k] - uprev
                        if mask[r - 1, c - 1]:
                            unew = (lam * (u[r + 1, c, k] + u[r - 1, c, k] +
                                           u[r, c + 1, k] + u[r, c - 1, k] +
                                           dx[r, c - 1, k] - dx[r, c, k] +
                                           dy[r - 1, c, k] - dy[r, c, k] -
                                           bx[r, c - 1, k] + bx[r, c, k] -
                                           by[r - 1, c, k] + by[r, c, k]
                                           ) + weight * image[r - 1, c - 1, k]
                                    ) / norm
                        else:
                            unew = (u[r + 1, c, k] + u[r - 1, c, k] +
                                    u[r, c + 1, k] + u[r, c - 1, k] +
                                    dx[r, c - 1, k] - dx[r, c, k] +
                                    dy[r - 1, c, k] - dy[r, c, k] -
                                    bx[r, c - 1, k] + bx[r, c, k] -
                                    by[r - 1, c, k] + by[r, c, k]
                                    ) / 4.0
                        u[r, c, k] = unew
                        rmse += (unew - uprev) ** 2
                        bxx = bx[r, c, k]
                        byy = by[r, c, k]

                        s = ux + bxx
                        if s > 1 / lam:
                            dxx = s - 1 / lam
                        elif s < -1 / lam:
                            dxx = s + 1 / lam
                        else:
                            dxx = 0
                        s = uy + byy
                        if s > 1 / lam:
                            dyy = s - 1 / lam
                        elif s < -1 / lam:
                            dyy = s + 1 / lam
                        else:
                            dyy = 0
                        dx[r, c, k] = dxx
                        dy[r, c, k] = dyy
                        bx[r, c, k] += ux - dxx
                        by[r, c, k] += uy - dyy

            rmse = np.sqrt(rmse / total)
            i += 1

        return np.squeeze(np.asarray(u[1:-1, 1:-1]))

    def tvm(input_array, keep_prob=0.8, lambda_tv=0.03):
        # expand input arr to 3-dims
        input_array = np.expand_dims(input_array, axis=2)
        mask = np.random.uniform(size=input_array.shape[:2])
        mask = mask < keep_prob
        return bregman(input_array, mask, weight=2.0/lambda_tv)

    def recover(mask, tmp_q_table, mask_prob=0.5):
        observed_q_table = mask * tmp_q_table
        recovered_q_table = SoftImpute(verbose=False).fit_transform(observed_q_table)

        return recovered_q_table
    """
end

const VTOL = 1e-6
const MAXITER = 10 #200
const GAMMA = 0.95
const SAVELOC = "../data/qvalue.csv"

mutable struct MDP
    S           :: RectangleGrid
    nstate      :: Int64
    A           :: Vector{Float64}
    naction     :: Int64
    transition  :: Function
    reward      :: Function

    function MDP(S, A, transition, reward)
        return new(S, length(S), A, length(A), transition, reward)
    end
end

mutable struct Policy
    Q       :: Matrix{Float64}
    A       :: Vector{Float64}
    naction :: Int64
    qvals   :: Vector{Float64}

    function Policy(Q, A)
        return new(Q, A, length(A), zeros(length(A)))
    end
end

function policy!(policy::Policy, belief::SparseMatrixCSC{Float64, Int64})
    fill!(policy.qvals, 0.0)
    for iaction in 1:policy.naction
        for ib in 1:length(belief.rowval)
            policy.qvals[iaction] +=
                belief.nzval[ib] * policy.Q[belief.rowval[ib], iaction]
        end
    end
    ibest = argmax(policy.qvals)
    return policy.A[ibest]
end

function saveq(Q::Matrix{Float64}, saveloc::String)
    writedlm(saveloc, Q, ',')
end

function value_iteration(mdp::MDP, save::Bool=false, saveloc::String=SAVE_LOC, verbose::Bool=true)
    nstate = length(mdp.S)
    naction = length(mdp.A)
    V = zeros(nstate)
    Q = zeros(nstate, naction)
    state = zeros(dimensions(mdp.S))
    cputime = 0.0

    println("Starting value iteration...")
    iter = 0
    for iter = 1:MAXITER
        tick()
        residual = 0.0
        prob = 0.2
        mask = py"generate_mask"(prob)

        # structured value-based planning
        for istate = 1:nstate
            for iaction = 1:naction
                if !isnan(mask[istate, iaction])
                    ind2x!(mdp.S, istate, state)
                    action = mdp.A[iaction]
                    snext = mdp.transition(state, action)
                    vnext = interpolate(mdp.S, V, snext)
                    Qprev = Q[istate, iaction]
                    Q[istate, iaction] = mdp.reward(state, action) + GAMMA * vnext
                    residual += (Q[istate, iaction] - Qprev)^2
                end
            end
        end

        # reconstruct via matrix estimation
        Q = py"recover"(mask, Q, prob)

        V = [maximum(Q[istate, :]) for istate = 1:nstate]

        iter_time = tok()
        cputime += iter_time
        residual /= (nstate * naction * prob)

        if verbose
            @printf("Iteration %d: average residual = %.2e, cputime = %.2e sec\n",
                    iter, residual, iter_time)
        end

        if residual < VTOL
            break
        end

        if iter == MAXITER
            println("Maximum number of iterations reached!")
        end
        print(iter)
    end
    @printf("Value iteration took %d iterations and %.2e sec\n", iter, cputime)

    if save
        saveq(Q, saveloc)
    end

    return Policy(Q, mdp.A)
end

function get_belief(mdp::MDP, state::Vector{Float64})
    belief = spzeros(mdp.nstate, 1)
    indices, weights = interpolants(mdp.S, state)
    belief[indices] = weights
    return belief
end

end

max_stars_repo_name: tsamtsam/SV-RL

max_issues_repo_name: tsamtsam/SV-RL

max_forks_repo_name: tsamtsam/SV-RL


text:
function slicerm(g, dim, x0, xargs; N=1000, w=0.5, m=10)
    #= multi-D slice sampler, g: log of distribution function,
    dim number of dimensions
    x0 start point,
    xargs extra arguments for distribution, N*dim samples to return,
        w "out step" distance, m max number of outward steps.

        Note that you only have to worry about log scaling in
	'vertical' comparisons
        in these Monte Carlo chain calculations.
    =#
    xs = zeros(N*dim,dim)  # array that will be returned
    xs[1,:] = x0  #x0 should be vector of length dim
    x1 = zeros(1,dim)
    L = zeros(1,dim)
    R = zeros(1,dim)
    way = zeros(1,dim)  # which axis to go along in space
    for i in 2:N*dim
        for d in 1:dim
            way = 0*way #clear it
            way[d] = 1.0 #set nonzero in direction we go for slicing
            y0 = g(x0,xargs)  #height of distribution at x0
            y = y0+log(rand()) # height for slice
            #start stepping out
            U = rand()   # between 0 and 1
            L = x0 - w*way*U
            R = L + w*way
            V = rand()
            J = floor(m*V)
            K = (m - 1) - J
            while J > 0 && y < g(L,xargs)
                L = L - w*way
                J = J - 1
            end
            while K > 0 && y < g(R,xargs)
                R = R + w*way
                K = K - 1
            end
            #now should be stepped out beyond distribution at slice level
            # work back in if no value found:
            Lbar, Rbar = L, R
            while true
                U = rand()
                x1 = Lbar + U * (Rbar - Lbar)  # vector subtraction should be
						# correct dir
                if y < g(x1,xargs)
                    break # exit while loop
                end
                if x1[d] < x0[d]
                    Lbar = x1
                    else
                    Rbar = x1
                end
            end
            xs[i,:] = x1 # found an acceptable point, record in chain (a row)
            x0 = x1 # set initial to new point for next round.
        end # for d
    end  #for N
    return xs
end

ext: jl

max_stars_repo_name: davidwcraig/slicer

max_issues_repo_name: davidwcraig/slicer

max_forks_repo_name: davidwcraig/slicer


text:
@testset "Dataset" begin
    res = fetch_dataset("WJP_STATE", vintage = "latest")
    @test res isa DataFrames.DataFrame
    @test "dates" in DataFrames.names(res)
    @test "vintage" in DataFrames.names(res)
    @test "series_id" in DataFrames.names(res)
    @test "values" in DataFrames.names(res)
    @test "State" in DataFrames.names(res)
    @test eltype(res.dates) <: Dates.Date
    @test eltype(res.vintage) <: Dates.Date
    @test eltype(res.series_id) <: String
    @test eltype(res.values) <: Union{Float64, Missing}

    res = fetch_dataset("WJP_STATE", restrictions = Dict("State" => ["TAS"]))
    @test res isa DataFrames.DataFrame
    @test "dates" in DataFrames.names(res)
    @test "vintage" in DataFrames.names(res)
    @test "series_id" in DataFrames.names(res)
    @test "values" in DataFrames.names(res)
    @test "State" in DataFrames.names(res)
    @test unique(res.State) == ["TAS"]

    res = fetch_dataset("WJP_STATE", restrictions = Dict("State" => ["TAS"]), dimensions = :name)
    @test res isa DataFrames.DataFrame
    @test "dates" in DataFrames.names(res)
    @test "vintage" in DataFrames.names(res)
    @test "series_id" in DataFrames.names(res)
    @test "values" in DataFrames.names(res)
    @test "State" in DataFrames.names(res)
    @test unique(res.State) == ["Tasmania"]

    @testset "Errors" begin
        @test_throws ErrorException("`*` is not a valid series code") fetch_dataset("*")
        @test_throws ErrorException("DataSet `valid_but_nonexistent` does not exist") fetch_dataset("valid_but_nonexistent")
        @test_throws ErrorException("Invalid vintage `abc`. Options are a date (in form YYYY-MM-DD) or one of `current` (alias `latest`) or `previous`") fetch_dataset("WJP_STATE", vintage = "abc")
        @test_throws ErrorException("`dimensions` can only be one of `:code` or `:name`") fetch_dataset("WJP_STATE", dimensions = :fizz)
    end
end


text:
@mpi_handle FileHandle MPI_File
const FILE_NULL = _FileHandle(MPI_FILE_NULL)
FileHandle() = FileHandle(FILE_NULL.val)

module File

import MPI: MPI, @mpichk, _doc_external, MPIPtr, libmpi, Finalized,
    Comm, MPI_Comm, FileHandle, MPI_File, Info, MPI_Info, FILE_NULL,
    Datatype, MPI_Datatype, MPI_Offset, Status, Buffer, Buffer_send
import Base: close

function open(comm::Comm, filename::AbstractString, amode::Cint, info::Info)
    file = FileHandle()
    # int MPI_File_open(MPI_Comm comm, const char *filename, int amode,
    #                   MPI_Info info, MPI_File *fh)
    @mpichk ccall((:MPI_File_open, libmpi), Cint,
                  (MPI_Comm, Cstring, Cint, MPI_Info, Ptr{MPI_File}),
                  comm, filename, amode, info, file)
    finalizer(close, file)
    file
end

"""
    MPI.File.open(comm::Comm, filename::AbstractString; keywords...)

Open the file identified by `filename`. This is a collective operation on `comm`.

Supported keywords are as follows:
 - `read`, `write`, `create`, `append` have the same behaviour and defaults as `Base.open`.
 - `sequential`: file will only be accessed sequentially (default: `false`)
 - `uniqueopen`: file will not be concurrently opened elsewhere (default: `false`)
 - `deleteonclose`: delete file on close (default: `false`)

Any additional keywords are passed via an [`Info`](@ref) object, and are implementation dependent.

# External links
$(_doc_external("MPI_File_open"))
"""
function open(comm::Comm, filename::AbstractString;
              read=nothing, write=nothing, create=nothing, append=nothing,
              sequential=false, uniqueopen=false, deleteonclose=false,
              infokws...)
    flags = Base.open_flags(read=read, write=write, create=create, truncate=nothing, append=append)
    amode = flags.read ? (flags.write ? MPI.MPI_MODE_RDWR : MPI.MPI_MODE_RDONLY) : (flags.write ? MPI.MPI_MODE_WRONLY : zero(Cint))
    if flags.write
        amode |= flags.create ? MPI.MPI_MODE_CREATE : MPI.MPI_MODE_EXCL
    end
    if flags.append
        amode |= MPI.MPI_MODE_APPEND
    end
    if sequential
        amode |= MPI.MPI_MODE_SEQUENTIAL
    end
    if uniqueopen
        amode |= MPI.MPI_MODE_UNIQUE_OPEN
    end
    if deleteonclose
        amode |= MPI.MPI_MODE_DELETE_ON_CLOSE
    end
    open(comm, filename, amode, Info(infokws...))
end

function close(file::FileHandle)
    if file.val != FILE_NULL.val && !Finalized()
        # int MPI_File_close(MPI_File *fh)
        @mpichk ccall((:MPI_File_close, libmpi), Cint,
                      (Ptr{MPI_File},), file)
    end
    return nothing
end

# View
"""
    MPI.File.set_view!(file::FileHandle, disp::Integer, etype::Datatype, filetype::Datatype, datarep::AbstractString; kwargs...)

Set the current process's view of `file`.

The start of the view is set to `disp`; the type of data is set to `etype`; the
distribution of data to processes is set to `filetype`; and the representation of data in
the file is set to `datarep`: one of `"native"` (default), `"internal"`, or `"external32"`.

# External links
$(_doc_external("MPI_File_set_view"))
"""
function set_view!(file::FileHandle, disp::Integer, etype::Datatype, filetype::Datatype, datarep::AbstractString, info::Info)
    # int MPI_File_set_view(MPI_File fh, MPI_Offset disp, MPI_Datatype etype,
    #                       MPI_Datatype filetype, const char *datarep, MPI_Info info)
    @mpichk ccall((:MPI_File_set_view, libmpi), Cint,
                  (MPI_File, MPI_Offset,  MPI_Datatype, MPI_Datatype, Cstring, MPI_Info),
                  file, disp, etype, filetype, datarep, info)
    return file
end

function set_view!(file::FileHandle, disp::Integer, etype::Datatype, filetype::Datatype, datarep="native"; infokwargs...)
    set_view!(file, disp, etype, filetype, datarep, Info(infokwargs...))
end

"""
    MPI.File.get_byte_offset(file::FileHandle, offset::Integer)

Converts a view-relative offset into an absolute byte position. Returns the absolute byte
position (from the beginning of the file) of `offset` relative to the current view of
`file`.

# External links
$(_doc_external("MPI_File_get_byte_offset"))
"""
function get_byte_offset(file::FileHandle, offset::Integer)
    r_displ = Ref{MPI_Offset}()
    # int MPI_File_get_byte_offset(MPI_File fh, MPI_Offset offset,
    #              MPI_Offset *disp)
    @mpichk ccall((:MPI_File_get_byte_offset, libmpi), Cint,
                  (MPI_File, MPI_Offset,  Ptr{MPI_Offset}),
                  file, offset, r_displ)
    r_displ[]
end

"""
    MPI.File.sync(fh::FileHandle)

A collective operation causing all previous writes to `fh` by the calling process to be
transferred to the storage device. If other processes have made updates to the storage
device, then all such updates become visible to subsequent reads of `fh` by the calling
process.

# External links
$(_doc_external("MPI_File_sync"))
"""
function sync(file::FileHandle)
    # int MPI_File_sync(MPI_File fh)
    @mpichk ccall((:MPI_File_sync, libmpi), Cint, (MPI_File,), file)
    return nothing
end

# I/O using individual file pointers
"""
    MPI.File.read!(file::FileHandle, data)

Reads current view of `file` into `data`. `data` can be a [`Buffer`](@ref), or
any object for which `Buffer(data)` is defined.

# See also
- [`MPI.File.set_view!`](@ref) to set the current view of the file
- [`MPI.File.read_all!`](@ref) for the collective operation

# External links
$(_doc_external("MPI_File_read"))
"""
function read!(file::FileHandle, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_File_read(MPI_File fh, void *buf,
    #                   int count, MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_read, libmpi), Cint,
                  (MPI_File, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
read!(file::FileHandle, data) = read!(file, Buffer(data))

"""
    MPI.File.read_all!(file::FileHandle, data)

Reads current view of `file` into `data`. `data` can be a [`Buffer`](@ref), or
any object for which `Buffer(data)` is defined. This is a collective operation, so must be
called on all ranks in the communicator on which `file` was opened.

# See also
- [`MPI.File.set_view!`](@ref) to set the current view of the file
- [`MPI.File.read!`](@ref) for the noncollective operation

# External links
$(_doc_external("MPI_File_read_all"))
"""
function read_all!(file::FileHandle, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_File_read_all(MPI_File fh, void *buf,
    #                       int count, MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_read_all, libmpi), Cint,
                  (MPI_File, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
read_all!(file::FileHandle, data) = read_all!(file, Buffer(data))

"""
    MPI.File.write(file::FileHandle, data)

Writes `data` to the current view of `file`. `data` can be a [`Buffer`](@ref),
or any object for which `Buffer_send(data)` is defined.

# See also
- [`MPI.File.set_view!`](@ref) to set the current view of the file
- [`MPI.File.write_all`](@ref) for the collective operation

# External links
$(_doc_external("MPI_File_write"))
"""
function write(file::FileHandle, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_File_write(MPI_File fh, const void *buf,
    #                    int count, MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_write, libmpi), Cint,
                  (MPI_File, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
write(file::FileHandle, data) = write(file, Buffer_send(data))

"""
    MPI.File.write_all(file::FileHandle, data)

Writes `data` to the current view of `file`. `data` can be a [`Buffer`](@ref),
or any object for which `Buffer_send(data)` is defined. This is a collective
operation, so must be called on all ranks in the communicator on which `file` was opened.

# See also
- [`MPI.File.set_view!`](@ref) to set the current view of the file
- [`MPI.File.write`](@ref) for the noncollective operation

# External links
$(_doc_external("MPI_File_write_all"))
"""
function write_all(file::FileHandle, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_File_write_all(MPI_File fh, const void *buf,
    #                    int count, MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_write_all, libmpi), Cint,
                  (MPI_File, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
write_all(file::FileHandle, data) = write_all(file, Buffer_send(data))

# Explicit offsets
"""
    MPI.File.read_at!(file::FileHandle, offset::Integer, data)

Reads from `file` at position `offset` into `data`. `data` can be a [`Buffer`](@ref), or
any object for which `Buffer(data)` is defined.

# See also
- [`MPI.File.read_at_all!`](@ref) for the collective operation

# External links
$(_doc_external("MPI_File_read_at"))
"""
function read_at!(file::FileHandle, offset::Integer, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_File_read_at(MPI_File fh, MPI_Offset offset, void *buf, int count,
    #                      MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_read_at, libmpi), Cint,
                  (MPI_File, MPI_Offset, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, offset, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
read_at!(file::FileHandle, offset::Integer, data) = read_at!(file, offset, Buffer(data))

"""
    MPI.File.read_at_all!(file::FileHandle, offset::Integer, data)

Reads from `file` at position `offset` into `data`. `data` can be a [`Buffer`](@ref), or
any object for which `Buffer(data)` is defined. This is a collective operation, so must be
called on all ranks in the communicator on which `file` was opened.

# See also
- [`MPI.File.read_at!`](@ref) for the noncollective operation

# External links
$(_doc_external("MPI_File_read_at_all"))
"""
function read_at_all!(file::FileHandle, offset::Integer, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)

    # int MPI_File_read_at_all(MPI_File fh, MPI_Offset offset, void *buf,
    #                          int count, MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_read_at_all, libmpi), Cint,
                  (MPI_File, MPI_Offset, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, offset, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
read_at_all!(file::FileHandle, offset::Integer, data) = read_at_all!(file, offset, Buffer(data))

"""
    MPI.File.write_at(file::FileHandle, offset::Integer, data)

Writes `data` to `file` at position `offset`. `data` can be a [`Buffer`](@ref), or any
object for which [`Buffer_send(data)`](@ref) is defined.

# See also
- [`MPI.File.write_at_all`](@ref) for the collective operation

# External links
$(_doc_external("MPI_File_write_at"))
"""
function write_at(file::FileHandle, offset::Integer, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_File_write_at(MPI_File fh, MPI_Offset offset, const void *buf,
    #                       int count, MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_write_at, libmpi), Cint,
                  (MPI_File, MPI_Offset, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, offset, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
write_at(file::FileHandle, offset::Integer, data) = write_at(file, offset, Buffer_send(data))

"""
    MPI.File.write_at_all(file::FileHandle, offset::Integer, data)

Writes from `data` to `file` at position `offset`. `data` can be a [`Buffer`](@ref), or
any object for which [`Buffer_send(data)`](@ref) is defined. This is a collective
operation, so must be called on all ranks in the communicator on which `file` was opened.

# See also
- [`MPI.File.write_at`](@ref) for the noncollective operation

# External links
$(_doc_external("MPI_File_write_at_all"))
"""
function write_at_all(file::FileHandle, offset::Integer, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_File_write_at_all(MPI_File fh, MPI_Offset offset, const void *buf,
    #                           int count, MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_write_at_all, libmpi), Cint,
                  (MPI_File, MPI_Offset, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, offset, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
write_at_all(file::FileHandle, offset::Integer, data) = write_at_all(file, offset, Buffer_send(data))

# Shared file pointer
"""
    MPI.File.read_shared!(file::FileHandle, data)

Reads from `file` using the shared file pointer into `data`.  `data` can be a
[`Buffer`](@ref), or any object for which `Buffer(data)` is defined.

# See also
- [`MPI.File.read_ordered!`](@ref) for the collective operation

# External links
$(_doc_external("MPI_File_read_shared"))
"""
function read_shared!(file::FileHandle, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_File_read_shared(MPI_File fh, void *buf, int count,
    #              MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_read_shared, libmpi), Cint,
                  (MPI_File, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
read_shared!(file::FileHandle, data) = read_shared!(file, Buffer(data))

"""
    MPI.File.write_shared(file::FileHandle, data)

Writes to `file` using the shared file pointer from `data`. `data` can be a
[`Buffer`](@ref), or any object for which `Buffer(data)` is defined.

# See also
- [`MPI.File.write_ordered`](@ref) for the collective operation

# External links
$(_doc_external("MPI_File_write_shared"))
"""
function write_shared(file::FileHandle, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_File_write_shared(MPI_File fh, const void *buf, int count,
    #          MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_write_shared, libmpi), Cint,
                  (MPI_File, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
write_shared(file::FileHandle, buf) = write_shared(file, Buffer_send(buf))

# Shared collective operations
"""
    MPI.File.read_ordered!(file::FileHandle, data)

Collectively reads in rank order from `file` using the shared file pointer into `data`.
`data` can be a [`Buffer`](@ref), or any object for which `Buffer(data)` is defined. This
is a collective operation, so must be called on all ranks in the communicator on which
`file` was opened.

# See also
- [`MPI.File.read_shared!`](@ref) for the noncollective operation

# External links
$(_doc_external("MPI_File_read_ordered"))
"""
function read_ordered!(file::FileHandle, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_File_read_ordered(MPI_File fh, void *buf, int count,
    #              MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_read_ordered, libmpi), Cint,
                  (MPI_File, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
read_ordered!(file::FileHandle, data) = read_ordered!(file, Buffer(data))

"""
    MPI.File.write_ordered(file::FileHandle, data)

Collectively writes in rank order to `file` using the shared file pointer from `data`.
`data` can be a [`Buffer`](@ref), or any object for which `Buffer(data)` is defined. This
is a collective operation, so must be called on all ranks in the communicator on which
`file` was opened.

# See also
- [`MPI.File.write_shared`](@ref) for the noncollective operation

# External links
$(_doc_external("MPI_File_write_ordered"))
"""
function write_ordered(file::FileHandle, buf::Buffer)
    stat_ref = Ref{Status}(MPI.STATUS_EMPTY)
    # int MPI_File_write_ordered(MPI_File fh, const void *buf, int count,
    #              MPI_Datatype datatype, MPI_Status *status)
    @mpichk ccall((:MPI_File_write_ordered, libmpi), Cint,
                  (MPI_File, MPIPtr, Cint, MPI_Datatype, Ptr{Status}),
                  file, buf.data, buf.count, buf.datatype, stat_ref)
    return stat_ref[]
end
write_ordered(file::FileHandle, buf) = write_ordered(file, Buffer_send(buf))

# seek
@enum Seek begin
    SEEK_SET = MPI.MPI_SEEK_SET
    SEEK_CUR = MPI.MPI_SEEK_CUR
    SEEK_END = MPI.MPI_SEEK_END
end

"""
    MPI.File.seek_shared(file::FileHandle, offset::Integer, whence::Seek=SEEK_SET)

Updates the shared file pointer according to `whence`, which has the following possible
values:

- `MPI.File.SEEK_SET` (default): the pointer is set to `offset`
- `MPI.File.SEEK_CUR`: the pointer is set to the current pointer position plus `offset`
- `MPI.File.SEEK_END`: the pointer is set to the end of file plus `offset`

This is a collective operation, and must be called with the same value on all processes in
the communicator.

# External links
$(_doc_external("MPI_File_seek_shared"))
"""
function seek_shared(file::FileHandle, offset::Integer, whence::Seek=SEEK_SET)
    # int MPI_File_seek_shared(MPI_File fh, MPI_Offset offset, int whence)
    @mpichk ccall((:MPI_File_seek_shared, libmpi), Cint,
                  (MPI_File, MPI_Offset, Cint),
                  file, offset, whence)
end

"""
    MPI.File.get_position_shared(file::FileHandle)

The current position of the shared file pointer (in `etype` units) relative to the current view.

# External links
$(_doc_external("MPI_File_get_position_shared"))
"""
function get_position_shared(file::FileHandle)
    r = Ref{MPI_Offset}()
    # int MPI_File_get_position_shared(MPI_File fh, MPI_Offset *offset)
    @mpichk ccall((:MPI_File_get_position_shared, libmpi), Cint,
                  (MPI_File, Ptr{MPI_Offset}), file, r)
    return r[]
end

end # module


text:
# This implements the variant when A is dense and m > n

using LinearAlgebra

struct IndGraphSkinny{T <: RealOrComplex} <: IndGraph
  m::Int
  n::Int
  A::Array{T, 2}
  AA::Array{T, 2}
  F::Cholesky{T, Array{T, 2}} #LL factorization
  tmp::Array{T, 1}
end

function IndGraphSkinny(A::Array{T, 2}) where {T <: RealOrComplex}
  m, n = size(A)
  AA = A' * A

  F = cholesky(AA + I)
  #normrows = vec(sqrt.(sum(abs2.(A), 2)))

  #The tmp vector assumes that difference between m and n is not drastic.
  # If someone will decide to solve it using m >> n, then tmp vector might be
  # considered to have only n position required to prox esitmation and indicator
  # calculation might be converted to less efficient.
  tmp = Array{T, 1}(undef, m)
  IndGraphSkinny(m, n, A, AA, F, tmp)
end

function (f::IndGraphSkinny)(x::AbstractArray{T}, y::AbstractArray{T}) where
    {T <: RealOrComplex}
  # the tolerance in the following line should be customizable
  mul!(f.tmp, f.A, x)
  f.tmp .-= y
  if norm(f.tmp, Inf) <= 1e-10
    return 0.0
  end
  return +Inf
end

function prox!(
    x::AbstractArray{T},
    y::AbstractArray{T},
    f::IndGraphSkinny,
    c::AbstractArray{T},
    d::AbstractArray{T},
    gamma=1.0
    ) where {T <: RealOrComplex}

  # x[:] = f.F \ (c + f.A' * d)
  mul!(x, adjoint(f.A), d)
  x .+= c
  ldiv!(f.F, x)

  mul!(y, f.A, x)
  return 0.0
end

fun_name(f::IndGraphSkinny) = "Indicator of an operator graph defined by dense full column rank matrix"
# fun_dom(f::IndGraph) = "Array{Real,1}, Array{Complex,1}"
# fun_expr(f::IndGraph) = "x,y ↦ 0 if Ax = y, +∞ otherwise"
# fun_params(f::IndGraph) =
#   string( "A = ", typeof(f.A), " of size ", size(f.A))

function prox_naive(
    f::IndGraphSkinny,
    c::AbstractArray{T},
    d::AbstractArray{T},
    gamma=1.0
    ) where {T <: RealOrComplex}

  x = f.F \ (c + f.A' * d)
  return x, f.A * x, 0.0
end

function (f::IndGraphSkinny)(xy::AbstractVector{T}) where
  {T <: RealOrComplex}
  x, y = splitinput(f, xy)
  return f(x, y)
end

(f::IndGraphSkinny)(xy::Tuple{AbstractVector{T}, AbstractVector{T}}) where
  {T <: RealOrComplex} = f(xy[1], xy[2])


text:
@testset "Testing read_dataset.jl" begin

    @testset "Reading matrix file" begin
        instance = read_dataset("datasets/movie.dat",cost_matrix = true)
        @test instance.file_name == "datasets/movie.dat"
        @test instance.total_itens == 10
        @test instance.total_voters == 0
        @test instance.interval_init == [i * (1/10) for i in 0:9]
        @test size(instance.item_names) == (10,)
        @test instance.item_names == string.(collect(1:10))
        @test size(instance.C) == (10,10)
        @test instance.C == [0   3   3   5   5   5   5   5   5  5;
                            3   0   5   5   5   5   5   5   5  5;
                            -1   1   0   1   3   5   3   5   5  5;
                            1  -1   1   0   5   5   5   5   5  5;
                            -5  -5  -1  -1   0   5   5   5   5  5;
                            -5  -5  -3  -5  -5   0  -5   5   5  5;
                            -5  -5  -1  -3   3   5   0   5   5  5;
                            -5  -5  -5  -5  -5   3  -5   0   5  5;
                            -5  -5  -5  -5  -5  -1  -5  -1   0  5;
                            -5  -5  -5  -5  -5  -5  -5  -5  -1  0]
        @test size(instance.D) == (10,10)
        @test instance.D == [0    0     4    4   10   10   10   10  10  10;
                            0    0     4    6   10   10   10   10  10  10;
                            -4   -4     0    0    4    8    4   10  10  10;
                            -4   -6     0    0    6   10    8   10  10  10;
                            -10  -10   -4   -6    0   10    2   10  10  10;
                            -10  -10   -8  -10  -10    0  -10    2   6  10;
                            -10  -10   -4   -8   -2   10    0   10  10  10;
                            -10  -10  -10  -10  -10   -2  -10    0   6  10;
                            -10  -10  -10  -10  -10   -6  -10   -6   0   6;
                            -10  -10  -10  -10  -10  -10  -10  -10  -6   0]
    end;

    @testset "Reading soc file" begin
        instance = read_dataset("datasets/movie.soc")
        @test instance.file_name == "datasets/movie.soc"
        @test instance.total_itens == 10
        @test instance.total_voters == 5
        @test instance.interval_init == [i * (1/10) for i in 0:9]
        @test size(instance.item_names) == (10,)
        @test instance.item_names == ["The Shawshank Redemption (1994)",
                                      "The Godfather (1972)",
                                      "The Godfather: Part II (1974)",
                                      "Schindler's List (1993)",
                                      "Pulp Fiction (1994)",
                                      "Rocky (1976)",
                                      "Her (2013)",
                                      "Mission: Impossible II (2000)",
                                      "Cats (2019)",
                                      "Disaster Movie (2008)"]
        @test size(instance.C) == (10,10)
        @test instance.C == [0   3   3   5   5   5   5   5   5  5;
                            -3   0   5   5   5   5   5   5   5  5;
                            -3  -5   0   1   3   5   3   5   5  5;
                            -5  -5  -1   0   5   5   5   5   5  5;
                            -5  -5  -3  -5   0   5   5   5   5  5;
                            -5  -5  -5  -5  -5   0  -5   5   5  5;
                            -5  -5  -3  -5  -5   5   0   5   5  5;
                            -5  -5  -5  -5  -5  -5  -5   0   5  5;
                            -5  -5  -5  -5  -5  -5  -5  -5   0  5;
                            -5  -5  -5  -5  -5  -5  -5  -5  -5  0]
        @test size(instance.D) == (10,10)
        @test instance.D == [0    6    6   10   10   10   10   10   10  10;
                            -6    0   10   10   10   10   10   10   10  10;
                            -6  -10    0    2    6   10    6   10   10  10;
                           -10  -10   -2    0   10   10   10   10   10  10;
                           -10  -10   -6  -10    0   10   10   10   10  10;
                           -10  -10  -10  -10  -10    0  -10   10   10  10;
                           -10  -10   -6  -10  -10   10    0   10   10  10;
                           -10  -10  -10  -10  -10  -10  -10    0   10  10;
                           -10  -10  -10  -10  -10  -10  -10  -10    0  10;
                           -10  -10  -10  -10  -10  -10  -10  -10  -10   0]
    end;

    @testset "Reading soi file" begin
        instance = read_dataset("datasets/movie.soi")
        @test instance.file_name == "datasets/movie.soi"
        @test instance.total_itens == 10
        @test instance.total_voters == 5
        @test instance.interval_init == [i * (1/10) for i in 0:9]
        @test size(instance.item_names) == (10,)
        @test instance.item_names == ["The Shawshank Redemption (1994)",
                                      "The Godfather (1972)",
                                      "The Godfather: Part II (1974)",
                                      "Schindler's List (1993)",
                                      "Pulp Fiction (1994)",
                                      "Rocky (1976)",
                                      "Her (2013)",
                                      "Mission: Impossible II (2000)",
                                      "Cats (2019)",
                                      "Disaster Movie (2008)"]
        @test size(instance.C) == (10,10)
        @test instance.C == [0   3   4   3   2   3   4   3   3  3;
                            -3   0   4   2   2   4   3   4   3  3;
                            -4  -4   0  -1   2   4   2   4   4  4;
                            -3  -2   1   0   2   2   3   2   2  2;
                            -2  -2  -2  -2   0   2   2   2   1  1;
                            -3  -4  -4  -2  -2   0  -3   4   3  3;
                            -4  -3  -2  -3  -2   3   0   3   3  3;
                            -3  -4  -4  -2  -2  -4  -3   0   3  3;
                            -3  -3  -4  -2  -1  -3  -3  -3   0  4;
                            -3  -3  -4  -2  -1  -3  -3  -3  -4  0]
        @test size(instance.D) == (10,10)
        @test instance.D == [0   6   8   6   4   6   8   6   6  6
                            -6   0   8   4   4   8   6   8   6  6
                            -8  -8   0  -2   4   8   4   8   8  8
                            -6  -4   2   0   4   4   6   4   4  4
                            -4  -4  -4  -4   0   4   4   4   2  2
                            -6  -8  -8  -4  -4   0  -6   8   6  6
                            -8  -6  -4  -6  -4   6   0   6   6  6
                            -6  -8  -8  -4  -4  -8  -6   0   6  6
                            -6  -6  -8  -4  -2  -6  -6  -6   0  8
                            -6  -6  -8  -4  -2  -6  -6  -6  -8  0]
    end;

    @testset "Reading toc file" begin
        instance = read_dataset("datasets/movie.toc")
        @test instance.file_name == "datasets/movie.toc"
        @test instance.total_itens == 10
        @test instance.total_voters == 5
        @test instance.interval_init == [i * (1/10) for i in 0:9]
        @test size(instance.item_names) == (10,)
        @test instance.item_names == ["The Shawshank Redemption (1994)",
                                    "The Godfather (1972)",
                                    "The Godfather: Part II (1974)",
                                    "Schindler's List (1993)",
                                    "Pulp Fiction (1994)",
                                    "Rocky (1976)",
                                    "Her (2013)",
                                    "Mission: Impossible II (2000)",
                                    "Cats (2019)",
                                    "Disaster Movie (2008)"]
        @test size(instance.C) == (10,10)
        @test instance.C == [0   3   3   5   5   5   5   5   5  5;
                            3   0   5   5   5   5   5   5   5  5;
                            -1   1   0   1   3   5   3   5   5  5;
                            1  -1   1   0   5   5   5   5   5  5;
                            -5  -5  -1  -1   0   5   5   5   5  5;
                            -5  -5  -3  -5  -5   0  -5   5   5  5;
                            -5  -5  -1  -3   3   5   0   5   5  5;
                            -5  -5  -5  -5  -5   3  -5   0   5  5;
                            -5  -5  -5  -5  -5  -1  -5  -1   0  5;
                            -5  -5  -5  -5  -5  -5  -5  -5  -1  0]
        @test size(instance.D) == (10,10)
        @test instance.D == [0    0     4    4   10   10   10   10  10  10;
                            0    0     4    6   10   10   10   10  10  10;
                            -4   -4     0    0    4    8    4   10  10  10;
                            -4   -6     0    0    6   10    8   10  10  10;
                            -10  -10   -4   -6    0   10    2   10  10  10;
                            -10  -10   -8  -10  -10    0  -10    2   6  10;
                            -10  -10   -4   -8   -2   10    0   10  10  10;
                            -10  -10  -10  -10  -10   -2  -10    0   6  10;
                            -10  -10  -10  -10  -10   -6  -10   -6   0   6;
                            -10  -10  -10  -10  -10  -10  -10  -10  -6   0]
    end;

    @testset "Reading toi file" begin
        instance = read_dataset("datasets/movie.toi")
        @test instance.file_name == "datasets/movie.toi"
        @test instance.total_itens == 6
        @test instance.total_voters == 5
        @test instance.interval_init == [i * (1/6) for i in 0:5]
        @test size(instance.item_names) == (6,)
        @test instance.item_names == ["The Shawshank Redemption (1994)",
                                    "The Godfather (1972)",
                                    "Schindler's List (1993)",
                                    "Pulp Fiction (1994)",
                                    "Rocky (1976)",
                                    "Cats (2019)"]
        @test size(instance.C) == (6,6)
        @test instance.C == [0   2   2   4  2  3;
                            2   0   2   5  3  3;
                            0   2   0   4  2  3;
                            -2  -5  -2   0  1  3;
                            -2  -3  -2   1  0  1;
                            -3  -3  -3  -3  1  0]
        @test size(instance.D) == (6,6)
        @test instance.D == [ 0    0   2   6  4  6;
                            0    0   0  10  6  6;
                            -2    0   0   6  4  6;
                            -6  -10  -6   0  0  6;
                            -4   -6  -4   0  0  0;
                            -6   -6  -6  -6  0  0]
    end;

end;

max_stars_repo_name: LuizHNLorena/AbrkgaObop

max_issues_repo_name: LuizHNLorena/AbrkgaObop

max_forks_repo_name: LuizHNLorena/AbrkgaObop


text:
using OMOPCommonDataModel

using DataFrames
using Documenter
using PrettyPrint
using StructArrays
using Tables
using Test

@testset "OMOPCommonDataModel.jl" begin
    @testset "Unit tests" begin
        @testset "version() and cdm_version()" begin
            @test OMOPCommonDataModel.version() isa VersionNumber
            @test OMOPCommonDataModel.version() > v"0"
            @test OMOPCommonDataModel.cdm_version() isa VersionNumber
            @test OMOPCommonDataModel.cdm_version() > v"0"
            @test OMOPCommonDataModel.cdm_version() == OMOPCommonDataModel.CDM_VERSION
        end
        @testset "String casing" begin
            @testset "pascalcase" begin
                @test OMOPCommonDataModel.Generate.pascalcase("HelloWorld") == "HelloWorld"
                @test OMOPCommonDataModel.Generate.pascalcase("helloWorld") == "HelloWorld"
                @test OMOPCommonDataModel.Generate.pascalcase("Hello World") == "HelloWorld"
                @test OMOPCommonDataModel.Generate.pascalcase("Hello_World") == "HelloWorld"
                @test OMOPCommonDataModel.Generate.pascalcase("HelloWorld") == "HelloWorld"
                @test OMOPCommonDataModel.Generate.pascalcase("hello world") == "HelloWorld"
                @test OMOPCommonDataModel.Generate.pascalcase("hello_world") == "HelloWorld"
                @test OMOPCommonDataModel.Generate.pascalcase("hello_World") == "HelloWorld"
                @test OMOPCommonDataModel.Generate.pascalcase("Hello_World") == "HelloWorld"
                @test OMOPCommonDataModel.Generate.pascalcase("Hello_world") == "HelloWorld"
                @test OMOPCommonDataModel.Generate.pascalcase("hello-world") == "HelloWorld"
            end
            @testset "snakecase" begin
                @test OMOPCommonDataModel.Generate.snakecase("HelloWorld") == "hello_world"
                @test OMOPCommonDataModel.Generate.snakecase("helloWorld") == "hello_world"
                @test OMOPCommonDataModel.Generate.snakecase("Hello World") == "hello_world"
                @test OMOPCommonDataModel.Generate.snakecase("Hello_World") == "hello_world"
                @test OMOPCommonDataModel.Generate.snakecase("HelloWorld") == "hello_world"
                @test OMOPCommonDataModel.Generate.snakecase("hello world") == "hello_world"
                @test OMOPCommonDataModel.Generate.snakecase("hello_world") == "hello_world"
                @test OMOPCommonDataModel.Generate.snakecase("hello_World") == "hello_world"
                @test OMOPCommonDataModel.Generate.snakecase("Hello_World") == "hello_world"
                @test OMOPCommonDataModel.Generate.snakecase("Hello_world") == "hello_world"
                @test OMOPCommonDataModel.Generate.snakecase("hello-world") == "hello_world"
            end
            @testset "spacecase" begin
                @test OMOPCommonDataModel.Generate.spacecase("HelloWorld") == "hello world"
                @test OMOPCommonDataModel.Generate.spacecase("helloWorld") == "hello world"
                @test OMOPCommonDataModel.Generate.spacecase("Hello World") == "hello world"
                @test OMOPCommonDataModel.Generate.spacecase("Hello_World") == "hello world"
                @test OMOPCommonDataModel.Generate.spacecase("HelloWorld") == "hello world"
                @test OMOPCommonDataModel.Generate.spacecase("hello world") == "hello world"
                @test OMOPCommonDataModel.Generate.spacecase("hello_world") == "hello world"
                @test OMOPCommonDataModel.Generate.spacecase("hello_World") == "hello world"
                @test OMOPCommonDataModel.Generate.spacecase("Hello_World") == "hello world"
                @test OMOPCommonDataModel.Generate.spacecase("Hello_world") == "hello world"
                @test OMOPCommonDataModel.Generate.spacecase("hello-world") == "hello world"
            end
        end
        @testset "String casing helper function" begin
            @testset "_remove_duplicate_spaces" begin
                @test OMOPCommonDataModel.Generate._remove_duplicate_spaces("foo    \n    \t                      bar") == "foo bar"
            end
            @testset "" begin
            end
        end
        @testset "Generate helper functions" begin
            @testset "_generate_export_statement" begin
                @test OMOPCommonDataModel.Generate._generate_export_statement("Foo", true) == "export Foo"
                @test OMOPCommonDataModel.Generate._generate_export_statement("Foo", false) == ""
            end
            @testset "_cdm_type_to_julia_partialtype" begin
                @test_throws ArgumentError OMOPCommonDataModel.Generate._cdm_type_to_julia_partialtype("foo_bar")
            end
            @testset "_generate_full_fieldtype" begin
                @test OMOPCommonDataModel.Generate._generate_full_fieldtype("SomeType", "NULL"; make_all_fields_optional = true) == "Union{SomeType, Missing} = missing # optional"
                @test OMOPCommonDataModel.Generate._generate_full_fieldtype("SomeType", "NOT NULL"; make_all_fields_optional = true) == "Union{SomeType, Missing} = missing # required"
                @test OMOPCommonDataModel.Generate._generate_full_fieldtype("SomeType", "NULL"; make_all_fields_optional = false) == "Union{SomeType, Missing} = missing # optional"
                @test OMOPCommonDataModel.Generate._generate_full_fieldtype("SomeType", "NOT NULL"; make_all_fields_optional = false) == "SomeType # required"
                @test_throws ArgumentError OMOPCommonDataModel.Generate._generate_full_fieldtype("SomeType", "FOO BAR"; make_all_fields_optional = false)
            end
        end
        @testset "Generate.generate()" begin
            generate_temp_dir = mktempdir()
            atexit(() -> rm(generate_temp_dir; force = true, recursive = true))
            generate_output_file = joinpath(abspath(generate_temp_dir), "autogenerated.jl")
            string_cases = [
                OMOPCommonDataModel.Generate.snakecase,
                OMOPCommonDataModel.Generate.pascalcase,
                identity,
                lowercase,
                uppercase,
                titlecase,
            ]
            for case in string_cases
                rm(generate_output_file; force = true, recursive = true)
                @test !ispath(generate_output_file)
                @test !isfile(generate_output_file)
                OMOPCommonDataModel.Generate.generate(; string_casing = case,
                                                        output_file = generate_output_file)
                @test ispath(generate_output_file)
                @test isfile(generate_output_file)
                generate_output_file_contents = read(generate_output_file, String)
                @test length(strip(generate_output_file_contents)) > 1_000
                # make sure that the generated file is parseable as Julia code
                @test_nowarn Base.parse_input_line(String(generate_output_file_contents))
                @test Base.parse_input_line(String(generate_output_file_contents)) isa Expr
            end
            rm(generate_temp_dir; force = true, recursive = true)
        end
    end
    @testset "Doctests" begin
        doctest(OMOPCommonDataModel)
    end
end


text:
using SpectralKit, PGFPlotsX

basis = smolyak_basis(Chebyshev, EndpointGrid(), Val(3),
                      let b = BoundedLinear(-1,1); (b, b); end)
g = grid(basis)
p = @pgf Axis({scale_only_axis = true,
               width = "6cm",
               height = "6cm" },
              Plot({ only_marks, red }, Coordinates(Tuple.(g))))
pgfsave("smolyak_grid.png", p; dpi = 600)


text:
using BenchmarkTools

# Get the number of iterations from the command line
n, = size(ARGS)
if n < 1
    println("Usage: TestGaussLegendreQuadrature.jl N")
    println("       ---> Please specify the number of grid points.")
    exit()
end

N = parse(Int, ARGS[1])

"""
    gauss(a, b, N)

Finds the approximate integral over the region [a, b]
with N iterations.
"""
function gauss(a, b, N)
    lambda, Q = eig(SymTridiagonal(zeros(N), [n/sqrt(4n^2 - 1) for n = 1:N-1]))
    return (lambda+1)*(b-a)/2 + a, [2*Q[1, i]^2 for i = 1:N]*(b-a)/2
end

println("--------------------------")
println(@sprintf "Gauss-Legendre Quadrature %d" N)
println("--------------------------")

@btime begin
   x, w = gauss(-3, 3, $N)
   quad = sum(exp.(x) .* w)
end

#exact = exp(3) - exp(-3)
#println(quad)
#println(exact)


text:
include("ispc_toolchain.jl")

# A basic ISPC kernel:
code = """
export void simple(uniform float vin[], uniform float vout[],
                   uniform int count) {
    foreach (index = 0 ... count) {
        float v = vin[index];
        if (v < 0.5)
            v = v * v;
        else
            v = sqrt(v);
        vout[index] = v;
    }
}
"""

println(code)

# Compile the code and get a function pointer to our kernel:
lib = load_ispc(code, `--target=avx1-i32x8`)
fptr = Libc.Libdl.dlsym(lib, "simple")

# Call the kernel:
vin = rand(Float32, 1000);
vout = zeros(Float32, 1000);
ccall(fptr, Cvoid, (Ref{Float32}, Ref{Float32}, UInt64), vin, vout, length(vout))

println(vout[1:10])


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()

# ## Setup
using Dojo
using JLD2
include("methods/simulation.jl")
include("methods/process.jl")
include("methods/benchmark.jl")

# ## Visualizer
vis = Visualizer()
open(vis)

################################################################################
## Astronaut
################################################################################
## humanoid
## multiple bodies
## no initial linear velocities
## no gravity
## no spring and damper
## random control
################################################################################

function ctrl!(mechanism, t)
	nu = input_dimension(mechanism)
	u = 0.3 * [zeros(6); mechanism.timestep; zeros(nu-7)]
	set_input!(mechanism, u)
	return
end

# ## Initialize
Random.seed!(0)
ϵ = 1.0e-14

# ## Mechanism
mech = get_mechanism(:humanoid,
	timestep=0.01,
	gravity=0.0,
	spring=0.0,
	damper=0.0,
	contact=false)

# ## Simulate
initialize!(mech, :humanoid)
storage = simulate!(mech, 1.0, ctrl!,
	record=true,
	opts=SolverOptions(rtol=ϵ, btol=ϵ))

# ## Visualize
visualize(mech, storage,
	vis=vis)

################################################################################
## Momentum
################################################################################

# ## Simulate
timestep= [0.10,]
mech, storage, tcompute = astronaut_simulation(;
	Nsim=1,
	timestep=1e-1,
	tsim=10.0,
	tctrl=1.0,
	seed=0,
	control_amplitude=0.05)

# ## Visualize
visualize(mech, storage[1], vis=vis)

# ## Benchmark
timestep= [0.10, 0.03, 0.01, 0.003, 0.001]
speed_dj, ener_dj = benchmark_energy(timestep; Nsim=1, tsim=2.0, tctrl=1.0, seed=0, ϵ=1e-14)
speed_dj, plin_dj, pang_dj = benchmark_momentum(timestep; Nsim=1, tsim=1.0, seed=0, ϵ=1e-14)

# ## Saving results
path = joinpath(@__DIR__, "results/astronaut_momentum.jld2")
jldsave(path, speed=speed_dj, plin=plin_dj, pang=pang_dj)

# ## Recover results
file = jldopen(joinpath(@__DIR__, "../MuJoCoBenchmark.jl/results/astronaut_momentum.jld2"))
speed_mj, plin_mj, pang_mj, ener_mj = file["speed"], file["plin"], file["pang"], file["ener"]
file = jldopen(joinpath(@__DIR__, "results/astronaut_momentum.jld2"))
speed_dj, plin_dj, pang_dj = file["speed"], file["plin"], file["pang"]

plt = plot(layout=(1,2), legend=false, xlims=(1e-2,2e4), ylims=(1e-15,1e2))
plot!(plt[1,1], speed_dj, plin_dj, axis=:log, yaxis=:flip, linewidth=3, markersize=6)
plot!(plt[1,2], speed_dj, pang_dj, axis=:log, yaxis=:flip, linewidth=3, markersize=6)
scatter!(plt[1,1], speed_dj, plin_dj, axis=:log, yaxis=:flip, linewidth=3, markersize=6)
scatter!(plt[1,2], speed_dj, pang_dj, axis=:log, yaxis=:flip, linewidth=3, markersize=6)

plot!(plt[1,1], speed_mj, plin_mj, axis=:log, yaxis=:flip, linewidth=3, markersize=6)
plot!(plt[1,2], speed_mj, pang_mj, axis=:log, yaxis=:flip, linewidth=3, markersize=6)
scatter!(plt[1,1], speed_mj, plin_mj, axis=:log, yaxis=:flip, linewidth=3, markersize=6)
scatter!(plt[1,2], speed_mj, pang_mj, axis=:log, yaxis=:flip, linewidth=3, markersize=6)

################################################################################
## Energy
################################################################################

# ## Simulate
timestep= [0.10, 0.02, 0.005]
ener_traj = []
for i = 1:3
	mech, storage, tcompute = astronaut_simulation(;
		Nsim=1,
		timestep=timestep[i],
		tsim=100.0,
		tctrl=1.0,
		seed=0,
		control_amplitude=0.02)
	stride = Int(floor(1 / timestep[i]))
	ener_t = [kinetic_energy(mech, storage[1], i) for i in stride+2:stride:length(storage[1])]
	push!(ener_traj, ener_t .- ener_t[1])
	@show i
end

# ## Saving results
path = joinpath(@__DIR__, "results/astronaut_energy.jld2")
jldsave(path, ener_traj=ener_traj)

# ## Recover results
file = jldopen(joinpath(@__DIR__, "../MuJoCoBenchmark.jl/results/astronaut_energy.jld2"))
ener_traj_mj = file["ener_traj"]
file = jldopen(joinpath(@__DIR__, "results/astronaut_energy.jld2"))
ener_traj_dj = file["ener_traj"]

# ## Plot
plt = plot(layout=(3,1), legend=false, xlims=(0,100))
for i = 1:3
	plot!(plt[i,1],
		range(0,stop=100,length=length(ener_traj_mj[i])),
		ener_traj_mj[i], linewidth=3, markersize=6)
	plot!(plt[i,1],
		range(0,stop=100,length=length(ener_traj_dj[i])),
		ener_traj_dj[i], linewidth=3, markersize=6)
end
display(plt)


text:
using Pkg.TOML, Libdl

config_toml = joinpath(first(DEPOT_PATH), "prefs", "MPI.toml")
mkpath(dirname(config_toml))

if !isfile(config_toml)
    touch(config_toml)
end

config = TOML.parsefile(config_toml)
update_config = false

# MPI.toml has 4 keys
#  binary   = "" (default) | "system" | "MPICH_jll" | "OpenMPI_jll" | "MicrosoftMPI_jll"
#  path     = "" (default) | top-level directory location
#  library  = "" (default) | library name/path | list of library names/paths
#  abi      = "" (default) | "MPICH" | "OpenMPI" | "MicrosoftMPI" | "unknown"
#  mpiexec  = "" (default) | executable name/path | [executable name/path, extra args...]

# 1. check if environment variables have changed
if haskey(ENV, "JULIA_MPI_BINARY")
    config["binary"] = ENV["JULIA_MPI_BINARY"]
    update_config = true
elseif haskey(ENV, "JULIA_MPI_LIBRARY") || haskey(ENV, "JULIA_MPI_PATH")
    config["binary"] = "system"
end

if haskey(ENV, "JULIA_MPI_PATH")
    config["path"] = ENV["JULIA_MPI_PATH"]
    update_config = true
end

if haskey(ENV, "JULIA_MPI_LIBRARY")
    config["library"] = ENV["JULIA_MPI_LIBRARY"]
    update_config = true
end

if haskey(ENV, "JULIA_MPI_ABI")
    config["abi"] = ENV["JULIA_MPI_ABI"]
    update_config = true
end

if haskey(ENV, "JULIA_MPIEXEC")
    if ENV["JULIA_MPIEXEC"] == ""
        delete!(config, "mpiexec")
    else
        config["mpiexec"] = ENV["JULIA_MPIEXEC"]
    end
    update_config = true
end

if update_config
    open(config_toml, write=true) do f
        TOML.print(f, config)
    end
end

binary = get(config, "binary", "")

# 2. generate deps.jl
if binary == "system"
    @info "using system MPI"
    library = get(config, "library", "")
    if library == ""
        library = ["libmpi", "libmpi_ibm", "msmpi", "libmpich"]
    end
    path    = get(config, "path", "")
    mpiexec = get(config, "mpiexec", path == "" ? "mpiexec" : joinpath(path, "bin", "mpiexec"))
    abi     = get(config, "abi", "")

    const libmpi = find_library(library, path == "" ? [] : [joinpath(path, "lib"), joinpath(path, "lib64")])
    if libmpi == ""
        error("libmpi could not be found")
    end

    const mpiexec_cmd = Cmd(mpiexec isa String ? [mpiexec] : mpiexec)

    _doc_external(fname) = ""
    include(joinpath("..","src","implementations.jl"))

    @info "Using implementation" libmpi mpiexec_cmd MPI_LIBRARY_VERSION_STRING

    if abi === ""
        # 3. check ABI
        impl, version = identify_implementation()
        if (impl == MPICH && version >= v"3.1" ||
            impl == IntelMPI && version > v"2014" ||
            impl == MVAPICH && version >= v"2" ||
            impl == CrayMPICH && version >= v"7")
            # https://www.mpich.org/abi/
            abi = "MPICH"
        elseif impl == OpenMPI || impl == IBMSpectrumMPI
            abi = "OpenMPI"
        elseif impl == MicrosoftMPI
            abi = "MicrosoftMPI"
        else
            abi = "unknown"
        end
        @info "MPI implementation detected" impl version abi
    else
        @info "MPI implementation config" abi
    end
    if abi == "MPICH"
        abi_incl = :(include("consts_mpich.jl"))
    elseif abi == "OpenMPI"
        abi_incl = :(include("consts_openmpi.jl"))
    elseif abi == "MicrosoftMPI"
        abi_incl = :(include("consts_microsoftmpi.jl"))
    else
        include("gen_consts.jl")
        abi_incl = :(include("consts.jl"))
    end

    deps = quote
        const libmpi = $libmpi
        const mpiexec_cmd = $mpiexec_cmd
        const mpiexec_path = mpiexec_cmd[1]
        _mpiexec(fn) = fn(mpiexec_cmd)
        $abi_incl

        using Requires

        function __init__deps()
            if Get_library_version() != $MPI_LIBRARY_VERSION_STRING
                error("MPI library has changed, please re-run `Pkg.build(\"MPI\")`.")
            end

            # error if a jll is loaded.
            @require(MPICH_jll        = "7cb0a576-ebde-5e09-9194-50597f1243b4",
                     error("MPICH_jll cannot be loaded: MPI.jl is configured to use the system MPI library"))
            @require(MicrosoftMPI_jll = "9237b28f-5490-5468-be7b-bb81f5f5e6cf",
                     error("MicrosoftMPI_jll cannot be loaded: MPI.jl is configured to use the system MPI library"))
            @require(OpenMPI_jll      = "fe0851c0-eecd-5654-98d4-656369965a5c",
                     error("OpenMPI_jll cannot be loaded: MPI.jl is configured to use the system MPI library"))
        end
    end
elseif binary == ""
    @info "using default MPI jll"
    deps = quote
        if Sys.iswindows()
            using MicrosoftMPI_jll
            const _mpiexec = MicrosoftMPI_jll.mpiexec
            const mpiexec_path = MicrosoftMPI_jll.mpiexec_path
            include("consts_microsoftmpi.jl")
        else
            using MPICH_jll
            const _mpiexec = MPICH_jll.mpiexec
            const mpiexec_path = MPICH_jll.mpiexec_path
            include("consts_mpich.jl")
        end

        function __init__deps()
            if (haskey(ENV, "SLURM_JOBID") || haskey(ENV, "PBS_JOBID") || haskey(ENV, "LSB_JOBID")) &&
                !(get(ENV, "JULIA_MPI_CLUSTER_WARN", "") == "n")
                @warn """
                    You appear to be using MPI.jl with the default MPI binary on a cluster.
                    We recommend using the system-provided MPI, see the Configuration section of the MPI.jl docs.
                    """
            end
        end
    end
elseif binary ==  "MPICH_jll"
    @info "using MPICH_jll"
    deps = quote
        using MPICH_jll
        include("consts_mpich.jl")
        const _mpiexec = MPICH_jll.mpiexec
        const mpiexec_path = MPICH_jll.mpiexec_path
        __init__deps() = nothing
    end
elseif binary ==  "OpenMPI_jll"
    @info "using OpenMPI_jll"
    deps = quote
        using OpenMPI_jll
        include("consts_openmpi.jl")
        const _mpiexec = OpenMPI_jll.mpiexec
        const mpiexec_path = OpenMPI_jll.mpiexec_path

        function __init__deps()
            # Required for OpenMPI relocateable binaries
            # TODO: this should be done in OpenMPI_jll package
            # https://github.com/JuliaPackaging/Yggdrasil/issues/390
            ENV["OPAL_PREFIX"] = OpenMPI_jll.artifact_dir
        end
    end
elseif binary ==  "MicrosoftMPI_jll"
    @info "using MicrosoftMPI_jll"
    deps = quote
        using MicrosoftMPI_jll
        include("consts_microsoftmpi.jl")
        const _mpiexec = MicrosoftMPI_jll.mpiexec
        const mpiexec_path = MicrosoftMPI_jll.mpiexec_path

        __init__deps() = nothing
    end
else
    error("Unknown binary $binary")
end

remove_line_numbers(x) = x
function remove_line_numbers(ex::Expr)
    if ex.head == :macrocall
        ex.args[2] = nothing
    else
        ex.args = [remove_line_numbers(arg) for arg in ex.args if !(arg isa LineNumberNode)]
    end
    return ex
end

# only update deps.jl if it has changed.
# allows users to call Pkg.build("MPI") without triggering another round of precompilation
deps_str = string(remove_line_numbers(deps))

if !isfile("deps.jl") || deps_str != read("deps.jl", String)
    write("deps.jl", deps_str)
end


text:
module XSV
#    using
#        DataFrames,

    export
        # types

        # methods
        iterxsv,
        iterxsvh,
        iterxsvht,
        fiterxsv,
        fiterxsvh,
        fiterxsvht,
        freadxsv,
        genxsv,
        readxsv

    include("genxsv.jl")
    include("readxsv.jl")
end


text:
"""
    _reflect(a,b,s,e)

reflection of the Brownian motion `B` where `B_{t-1} = a` and  `B_{t} = b`
on the hypercube `[s,e]^d` where `d = size(a,1)`
"""
function _reflect(a::T, b::T, s::T, e::T) where T <: Vector
    r = 2; n = zeros(size(a))
    # first checking if b is in the hypercube
    all((a .>= s) .& (a .<= e)) ? nothing : error("a = $a not in hypercube")
    size(a) == size(b) ? nothing : error("a not same dim as b")
    for i in 1:length(a)
        if b[i] < s[i]
            rtemp = (a[i] - s[i]) / (a[i] - b[i])
            if rtemp < r
                r = rtemp
                n .= 0
                n[i] = -1
            end
        elseif  b[i] > e[i]
            rtemp =  (e[i] - a[i]) / (b[i]- a[i])
            if rtemp < r
                r = rtemp
                n .= 0
                n[i] = 1
            end
        end
    end
    while r < 1
        c = a + r * ( b - a )
        # dat = hcat(a,c)
        # Plots.plot3d!(dat[1,:],dat[2,:],dat[3,:],label = "",color="blue")
        a = c
        b = b - 2 * n * ( dot(b-c,n))
        r = 2;
        for i in 1:length(a)
            if b[i] < s[i]
                rtemp = (a[i] - s[i]) / (a[i] - b[i])
                if rtemp < r
                    r = rtemp
                    n .= 0
                    n[i] = -1
                end
            elseif  b[i] > e[i]
                rtemp =  (e[i] - a[i]) / (b[i]- a[i])
                if rtemp < r
                    r = rtemp
                    n .= 0
                    n[i] = 1
                end
            end
        end
    end
    # dat = hcat(a,b)
    # Plots.plot3d!(dat[1,:],dat[2,:],dat[3,:],label = "",color="blue")
    return b
end

function _reflect(a::T, b::T, s, e) where T <: CuArray
    @assert all((a .>= s) .& (a .<= e)) "a = $a not in hypercube"
    @assert size(a) == size(b) "a not same dim as b"
    out1 = b .< s
    out2 = b .> e
    out = out1 .| out2
    n = similar(a)
    n .= 0
    # Allocating
    while any(out)
        rtemp1 = @. (s - a) #left
        rtemp2 = @. (e - a) #right
        div = @. (out * (b-a) + !out)
        rtemp = (rtemp1 .* out1 .+ rtemp2 .* out2) ./ div .+ (.!(out1 .| out2))
        rmin = minimum(rtemp,dims=1)
        n .= rtemp .== minimum(rtemp;dims=1)
        c = @. (a + (b-a) * rmin)
        b = @. ( b - 2 * n * (b-c) )
        a = c
        @. out1 = b < s
        @. out2 = b > e
        @. out = out1 | out2
    end
    return b
end


text:
# Code used to generate Lovasz Theta function problems. Save in Testproblems/Lovasz with an accurante
# solution from MOSEK

# Original problem has the following format:
# min_X  Tr(XJ)
# s.t.   Tr(X) = 1
#        X_ij = 1, if (i,j) ∈ E
#        X ⪴ 0

# workspace()
# using JLD, JuMP, Mosek

function createNonzeroExtractor(E)
  n = size(E,1)
  nz = nnz(sparse(E))
  S = spzeros(nz,n^2)
  e = vec(E)
  counter = 1
  for iii=1:length(e)
    if e[iii] == 1.
      S[counter,iii] = 1
      counter+=1
    end
  end
  return sparse(S)
end

# creates a random symmetric edge matrix with zeros on the diagonal and 1 as nonzero entries
function createRandomEdgeMatrix(rng,n,density)
  E  = full(Symmetric(sprand(rng,n,n,density)))
  map!(x->(x!=0 ? x=1. : x=0.),E,E)
  for iii=1:n
    E[iii,iii] = 0
  end
  return sparse(E)
end

rng = MersenneTwister(12345)
dirPath = "../DataFiles/Julia/Lovasz/"
!ispath(dirPath) && mkdir(dirPath)

nn = 25

 for iii =1:1:nn
  # choose size of problem
  n = rand(rng,10:30)
  dmin = 0.1
  dmax = 0.9
  density = rand(rng)*(dmax-dmin)+dmin

  # create random edge matrix E for graph
  E = createRandomEdgeMatrix(rng,n,density)

  # Test: create E for pentagon graph
  # n = 5
  # E = zeros(5,5)
  # E[1,2] = E[2,3] = E[3,4] = E[4,5] = E[5,1] = 1
  # E[1,5] = E[2,1] = E[3,2] = E[4,3] = E[5,4] = 1
  # correctVal = sqrt(5)

  nz = nnz(sparse(E))
  J = sparse(ones(n,n))

  # put problem into solver format
  q = -vec(J)
  P = spzeros(n^2,n^2)
  r = 0.
  A1 = vec(speye(n))'
  A2 = createNonzeroExtractor(E)
  A3 = -speye(n^2)

  Aa = [A1;A2;A3]
  ba = [1;spzeros(nz+n^2)]

  # specify cone
  Kf = 1+nz
  Kl = 0
  Kq = []
  Ks = [n^2]

  # # solve accurately once with mosek
  model = Model(solver=MosekSolver())
  @variable(model, X[1:n,1:n], SDP)
  @objective(model, Max ,trace(J*X))
  @constraint(model, trace(X) == 1)
  x = vec(X)
  @constraint(model, A2*x .== 0)
  status = JuMP.solve(model)

  # objective value wit minus since it's a max-objective
  objTrue = -getobjectivevalue(model)
  solTrue = getvalue(x)
  # println("Objective value: ",objTrue)
  # println("Corrected objective value: ", objCorrected)
  # println("x = ", getvalue(x))

  nr = "$(iii)"
  if iii < 10
    nr = "0$(iii)"
  end
  problemType = "Lovasz Theta Function"
  problemName = "Lovasz"*nr

  fn = "Lovasz"*nr*".jld"
  JLD.save(dirPath*fn,"n",size(Aa,2),"m",size(Aa,1),"A",Aa,"b",ba,"P",P,"q",q,"r",r,"E",E,"objTrue",objTrue,"solTrue",solTrue,"problemType",problemType,"problemName",problemName,"Kf",Kf,"Kl",Kl,"Kq",Kq,"Ks",Ks)
  println("$(iii)/$(nn) completed!")
end


text:
# Demignot, S., D., D.,
# Effect of prosthetic sugar groups on the pharmacokinetics of glucose-oxidase
# https://pubmed.ncbi.nlm.nih.gov/2855567/
using Logging

using StructuralIdentifiability

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

ode = @ODEmodel(
    x0'(t) = a1 * (x1(t) - x0(t)) - (ka * n * x0(t)) / (kc * ka + kc * x2(t) + ka * x0(t)),
    x1'(t) = a2 * (x0(t) - x1(t)),
    x2'(t) = b1 * (x3(t) - x2(t)) - (kc * n * x2(t)) / (kc * ka + kc * x2(t) + ka * x0(t)),
    x3'(t) = b2 * (x2(t) - x3(t)),
    y1(t) = x0(t)
)

@time println(assess_global_identifiability(ode))


text:
function _convert(ctx::ConverterContext, decl::CLMacroDefinition)
	isbuiltin(decl) && return
	# HACK:  builtin check doesnt seem to work on all builtins, but code locations with unknown file seems to do the trick
	CodeLocation(decl).file == "<unknown>" && return

	loc = string(CodeLocation(decl))
	name = _convertName(ctx, decl)
	func = _gensym(ctx, name)*"_"
	tokens = map(t -> t.text, tokenize(decl))

	args = []
	if isfunctionlike(decl) != 0
		tokens[2] == "(" || error("Macro was expected to be function-like.")

		i = 3
		while i <= length(tokens) && tokens[i] != ")"
			tokens[i] != "," && push!(args, _convertName(ctx, tokens[i]))
			i += 1
		end
		i += 1  # skip the closing paren ')'

		if !isempty(args) && args[end] == :(...)
			args[end] = "$(_gensym(ctx, "args"))..."
		end

		tokens = tokens[i:end]
	else
		tokens = tokens[2:end]
	end

	literal = nothing
	if !isempty(tokens) && tokens[1] == "(" && tokens[end] == ")"
		tokens = tokens[2:end-1]
	end

	# parse C string literal
	if isnothing(literal) && length(tokens) == 1
		literal = _parseString(tokens[1])
	end

	# parse C character literal
	if isnothing(literal) && length(tokens) == 1
		literal = _parseChar(tokens[1])
	end

	# parse C integer literal
	if isnothing(literal) && length(tokens) in (1, 2)
		literal = _parseInteger(join(tokens))
	end

	# parse C floating point literal
	if isnothing(literal) && length(tokens) == 1
		literal = _parseFloat(tokens[1])
	end

	# parse C symbol/token literal
	if isnothing(literal) && length(tokens) == 1
		literal = _parseSymbol(tokens[1])
		literal = isnothing(literal) ? literal : Symbol(_convertName(ctx, literal))
	end

	# parse C expressions
	if isnothing(literal)
		literal = _parseExpr(join(tokens, ' '))
	end

	# parse C expressions
	if isnothing(literal)
		literal = _parseExpr(join(tokens))
	end

	if isempty(tokens)
		body = "nothing"
	elseif !isnothing(literal)
		body = literal isa Symbol || literal isa Expr ? replace(repr(literal), r"#=.*=# " => "") : ":($(repr(literal)))"
	else
		body = """:(error("Macro not converted:  "*$(repr(join(tokens, ' ')))))"""
	end
	body = """function($(join(args, ", ")))
		$(body)
	end"""

	if !("@$(name)" in ctx.oneofs)
		push!(ctx.oneofs, "@$(name)")

		_export(ctx, "@$(name)")
		binding = "const $(func) = Ref{Function}($(body))"
		expr = "macro $(name)($(join(args, ", "))) ; $(func)[]($(join(args, ", "))) ; end"
		push!(ctx.converted, JuliaizedC(
			decl,
			binding,
			:atcompile,
		))
		push!(ctx.converted, JuliaizedC(
			decl,
			expr,
			:atcompile,
		))
	else
		expr = "$(func)[] = $(body)"
		push!(ctx.converted, JuliaizedC(
			decl,
			expr,
			:atcompile,
		))
	end
end


text:
# ---
# title: 866. Prime Palindrome
# id: problem866
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/prime-palindrome/description/>
# hidden: true
# ---
#
# Find the smallest prime palindrome greater than or equal to `N`.
#
# Recall that a number is _prime_ if it's only divisors are 1 and itself, and it
# is greater than 1.
#
# For example, 2,3,5,7,11 and 13 are primes.
#
# Recall that a number is a _palindrome_ if it reads the same from left to right
# as it does from right to left.
#
# For example, 12321 is a palindrome.
#
#
#
# **Example 1:**
#
#
#
#     Input: 6
#     Output: 7
#
#
# **Example 2:**
#
#
#
#     Input: 8
#     Output: 11
#
#
# **Example 3:**
#
#
#
#     Input: 13
#     Output: 101
#
#
#
# **Note:**
#
#   * `1 <= N <= 10^8`
#   * The answer is guaranteed to exist and be less than `2 * 10^8`.
#
#
## @lc code=start
using LeetCode

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


text:
@test has_cuda(true)
@test has_cuda_gpu(true)

# the API shouldn't have been initialized
@test CuCurrentContext() == nothing
@not_if_memcheck @test CuCurrentDevice() == nothing

task_cb = Any[nothing for tid in 1:Threads.nthreads()]
CUDA.attaskswitch() do
    task_cb[Threads.threadid()] = current_task()
end

device_switch_cb = Any[nothing for tid in 1:Threads.nthreads()]
CUDA.atdeviceswitch() do
    device_switch_cb[Threads.threadid()] = (dev=device(), ctx=context())
end

device_reset_cb = Any[nothing for tid in 1:Threads.nthreads()]
CUDA.atdevicereset() do dev
    device_reset_cb[Threads.threadid()] = dev
end

function reset_cb()
    fill!(task_cb, nothing)
    fill!(device_switch_cb, nothing)
    fill!(device_reset_cb, nothing)
end

# now cause initialization
ctx = context()
dev = device()
@test CuCurrentContext() == ctx
@test CuCurrentDevice() == dev
@test task_cb[1] == current_task()
@test device_switch_cb[1].ctx == ctx
@test device_switch_cb[1].dev == dev

reset_cb()

# ... on a different task
task = @async begin
    context()
end
@test ctx == fetch(task)
@test task_cb[1] == task
@test device_switch_cb[1].ctx == ctx
@test device_switch_cb[1].dev == dev

reset_cb()

# ... back to the main task
ctx = context()
dev = device()
@test task_cb[1] == current_task()
@test device_switch_cb[1] == nothing

device!(CuDevice(0))
device!(CuDevice(0)) do
    nothing
end

context!(ctx)
context!(ctx) do
    nothing
end

@test_throws ErrorException device!(0, CUDA.CU_CTX_SCHED_YIELD)

reset_cb()

device_reset!()

@test device_reset_cb[1] == CuDevice(0)

reset_cb()

device!(0, CUDA.CU_CTX_SCHED_YIELD)
@test task_cb[1] == nothing
@test device_switch_cb[1].dev == CuDevice(0)

# reset on a different task
let ctx = context()
    @test CUDA.isvalid(ctx)
    @test ctx == fetch(@async context())

    @sync @async device_reset!()

    @test CUDA.isvalid(context())
    @test ctx != context()
end

# test the device selection functionality
if length(devices()) > 1
    device!(0)
    device!(1) do
        @test device() == CuDevice(1)
    end
    @test device() == CuDevice(0)

    device!(1)
    @test device() == CuDevice(1)
end

# test that each task can work with devices independently from other tasks
if length(devices()) > 1
    device!(0)
    @test device() == CuDevice(0)

    task = @async begin
        device!(1)
        @test device() == CuDevice(1)
    end
    fetch(task)

    @test device() == CuDevice(0)

    # reset on a task
    task = @async begin
        device!(1)
        device_reset!()
    end
    fetch(task)
    @test device() == CuDevice(0)

    # tasks on multiple threads
    Threads.@threads for d in 0:1
        for x in 1:100  # give threads a chance to trample over each other
            device!(d)
            yield()
            @test device() == CuDevice(d)
            yield()

            sleep(rand(0.001:0.001:0.01))

            device!(1-d)
            yield()
            @test device() == CuDevice(1-d)
            yield()
        end
    end
    @test device() == CuDevice(0)
end

@test deviceid() >= 0
@test deviceid(CuDevice(0)) == 0
if length(devices()) > 1
    @test deviceid(CuDevice(1)) == 1
end


text:
blockchars = [
:lower1_8 => '\U2581',
:lower1_4 => '\U2582',
:lower3_8 => '\U2583',
:lower1_2 => '\U2584',
:lower5_8 => '\U2585',
:lower3_4 => '\U2586',
:lower7_8 => '\U2587',
:lowerfull => '\U2588',
]


text:
"""
```
parse_model(network::Network)
```

Construct a `state <: Vector{Int}` and `model <: ReactionSystem` from a `Network` object.

Use this in `simulate` to avoid constructing internal data structures over multiple samples.
"""
function parse_model(network::Network)
    state = [ x.population for (key, x) in species_list(network) ]
    model = ReactionSystem(network)

    return state, model
end

##### user-facing interface #####

"""
```
simulate(network::Network, algname::SimulationAlgorithm;
            [tfinal = 0.0],
            [rates_cache = HasRates],
            [save_points = nothing],
            [save_function = save_state],
            [ntrials = nothing])
```

Simulate a `Network` of interacting populations.
Note that simulations may terminate early if the cumulative intensity reaches `0` (that is, reaching an absorbing state).

#### Keyword Arguments

- `tfinal`: The final simulation time.
- `rates_cache`: Indicates the type of information stored in a `rates` cache. If `HasRates` is chosen, store the actual rates. If `HasSums` is chosen, store partial sums of the `rates`. This effectively toggles between linear and binary searches, provided the algorithm supports the option.
- `save_points`: An indexable collection indicating time points to sample and record system state. The default `nothing` forces saving after every reaction event.
- `save_function`: A three argument function that maps the three arguments `(simulator, state, model)` to data recorded in a `SamplePath`. Default behavior is to record each species.
- `ntrials`: Specifies the number of times to simulate the given model. The default `nothing` generates a single sample path.
"""
function simulate(network::Network, algname::SimulationAlgorithm;
    tfinal=0.0,
    rates_cache = HasRates,
    save_points = nothing,
    save_function::funcT = save_state,
    ntrials = nothing
    ) where funcT <: Function
    # build the internal representation of our stochastic process
    initial_state, model = parse_model(network)

    # feedforward down the chain...
    return simulate(initial_state, model, algname, tfinal, rates_cache, save_points, save_function, ntrials)
end

"""
```
simulate(state, model, algname::SimulationAlgorithm;
            [tfinal = 0.0],
            [rates_cache = HasRates],
            [save_points = nothing],
            [save_function = save_state],
            [ntrials = nothing])
```

Simulate a `model` with the given initial `state`.
Note that simulations may terminate early if the cumulative intensity reaches `0` (that is, reaching an absorbing state).

#### Arguments

For well-mixed systems:

- `state`: A `Vector{Int}` constructed by `parse_model`.
- `model`: A `ReactionSystem` type constructed by `parse_model`.

For lattice-based systems:

- `state`: A `Lattice` type whose topology comptabile is compatible with `model`.
- `model`: An `InteractingParticleSystem` type built from `@enumerate_with_sclass`.

#### Keyword Arguments

- `tfinal`: The final simulation time.
- `rates_cache`: Indicates the type of information stored in a `rates` cache. If `HasRates` is chosen, store the actual rates. If `HasSums` is chosen, store partial sums of the `rates`. This effectively toggles between linear and binary searches, provided the algorithm supports the option.
- `save_points`: An indexable collection indicating time points to sample and record system state. The default `nothing` forces saving after every reaction event.
- `save_function`: A three argument function that maps the three arguments `(simulator, state, model)` to data recorded in a `SamplePath`. Default behavior is to store population counts for well-mixed systems or `Configuration` objects for lattice-based systems.
- `ntrials`: Specifies the number of times to simulate the given model. The default `nothing` generates a single sample path.
"""
function simulate(initial_state, model, algname::SimulationAlgorithm;
    tfinal = 0.0,
    rates_cache = HasRates,
    save_points = nothing,
    save_function::funcT = save_state,
    ntrials = nothing
    ) where funcT <: Function
    # feedforward down the chain...
    return simulate(initial_state, model, algname, tfinal, rates_cache, save_points, save_function, ntrials)
end

##### internals #####

function simulate(initial_state, model, algname, tfinal, rates_cache, save_points, save_function, ntrials)
    # copy state
    state = copy(initial_state)

    # build the simulator
    simulator = build_simulator(algname, state, model, rates_cache)
    # TODO: figure out a way of initializing so that output is built correctly
    initialize!(simulator, state, model, tfinal)

    # build the output data
    output = build_output(save_function, simulator, state, model, ntrials)

    # feedforward down the chain...
    simulate!(simulator, state, initial_state, model, tfinal, output, save_points, save_function, ntrials)
end

# case: single trajectory
function simulate!(simulator, state, initial_state, model, tfinal, output, save_points, save_function, ::Nothing)
    initialize_datastructs!(state, initial_state, model)
    simulate!(simulator, state, model, tfinal, output, save_points, save_function)

    return output
end

# case: ensemble simulation
function simulate!(simulator, state, initial_state, model, tfinal, output, save_points, save_function, ntrials::Integer)
    for k in 1:ntrials
        initialize_datastructs!(state, initial_state, model)
        sample_path = output[k]
        simulate!(simulator, state, model, tfinal, sample_path, save_points, save_function)
    end

    return output
end

function simulate!(simulator, state, model, tfinal, output, save_points, save_function)
    initialize!(simulator, state, model, tfinal)

    in_timespan = simulator.t < tfinal
    transient = cumulative_intensity(simulator) > 0
    running = in_timespan && transient

    while running
        tnew = get_new_time(simulator)

        if tnew <= tfinal
            step!(simulator, state, model)
        else
            simulator.t = tfinal
        end
        # update output
        t = simulator.t
        update_samplepath!(save_function, output, simulator, t, state, model, save_points)

        # update flags
        in_timespan = simulator.t < tfinal
        transient = cumulative_intensity(simulator) > 0
        running = in_timespan && transient
    end

    # add end point in case of extinction
    if !transient
        update_samplepath!(save_function, output, simulator, tfinal, state, model, save_points)
    end

    return output
end

##### temporary hacks

initialize_datastructs!(state, initial_state, model) = copyto!(state, initial_state)

function initialize_datastructs!(lattice::Lattice, initial_lattice::Lattice, model::InteractingParticleSystem)
    # copy fields for lattice...
    empty!(lattice.site)
    for s in initial_lattice.site
        new_site = Site(s.label, State(s.state.ptype), tuple(s.coord...))
        push!(lattice.site, new_site)
    end

    empty!(lattice.coord_order)
    foreach(s -> push!(lattice.coord_order, s), lattice.site)
    sort!(lattice.coord_order,
    alg = Base.Sort.DEFAULT_STABLE,
    by = coordinates,
    lt = <)

    foreach(nb -> empty!(nb), lattice.neighbors)

    # unpack information
    number_init = number_sites(lattice)
    site = lattice.site
    nbtype = topology(lattice)
    enum = model.enumeration
    class = enum.class
    isactive = enum.isactive
    pair_to_classes = enum.pair_to_classes
    dummy_composition = enum.dummy_composition

    # make sure sample classes are empty
    foreach(c -> empty!(c), class)

    # add missing open sites
    for i in 1:number_init
        x = site[i]
        for coord in eachdir(nbtype, x)
            if !istracked(lattice, coord)
                y = spawn_new_site(lattice, coord)
                add_site!(lattice, y)
            end
        end
    end

    # build neighborhoods
    build_neighborhoods!(nbtype, lattice)

    # assign sites to classes in the enumeration
    for x in site
        # build the neighborhood composition around x
        build_composition!(dummy_composition, lattice, x)

        # what sample classes does x belong to?
        l = get_ptype(x)
        k = get_nbclass_index(enum, dummy_composition)

        if isactive[l]
            classes = pair_to_classes[(l,k)]

            # add x to each assigned class
            for s in classes
                add_class_member!(class[s], x)
            end
        end

        # change the neighborhood class of x
        change_neighbor_class!(x, k)
    end

    # clean up
    fill!(dummy_composition, 0)
    sort!(lattice.coord_order,
    alg = Base.Sort.DEFAULT_STABLE,
    by = coordinates,
    lt = <)

    return nothing
end


text:
function R_tree(tree::Tree ; file = "tree.Rdata")
    edges, branch_lengths, temp_rates, tip_labels = make_ape(tree)
    ntip = (tree.n_nodes + 1)/2

    @rput edges
    @rput branch_lengths
    @rput temp_rates
    @rput ntip
    @rput file
    @rput tip_labels

    reval("""
        temp_tree = list(edge = edges, Nnode = ntip - 1, edge.length = branch_lengths, tip.label = tip_labels)
        class(temp_tree) = "phylo"
    """)
end

"""
    save_ClaDS_in_R(co::CladsOutput, path::String ; maxit = Inf, ...)
sss
Save the output of a ClaDS run as a Rdata file.

# Arguments
- `co::CladsOutput`: the result of a run of `infer_ClaDS`.
- `path::String`: the path the file should be saved to.

# .Rdata file
The .Rdata file contains an object called `CladsOutput` that is a list with the following fields :
- `tree`: the phylogeny, saved as a `phylo` object
- `chains`: the MCMC chains
- `rtt_chains`: the MCMCs with the rate through time information
- `sig_map`: the MAP estimate for the σ parameter
- `al_map`: the MAP estimate for the α parameter
- `eps_map`: the MAP estimate for the ε parameter
- `lambda0_map`: the MAP estimate for the λ0 parameter
- `lambdai_map`: the MAP estimate for the branch-specific speciation rates
- `lambdatip_map`: the MAP estimate for the tip speciation rates
- `DTT_mean`: the point estimate for the diversity through time
- `RTT_map` : the rate through time estimates
- `time_points` : the times at which `DTT_mean` and `RTT_map` are computed
- `enhanced_trees` : a sample of the complete tree distribution, as a list of `phylo` objects
- `gelm` : the gelman parameter
"""
function save_ClaDS_in_R(co::CladsOutput, path::String)

    chains = co.chains
    rtt_chains = co.rtt_chains
    sig_map = co.σ_map
    al_map = co.α_map
    eps_map = co.ε_map
    lambda0_map = co.λ0_map
    lambdai_map = co.λi_map
    lambdatip_map = co.λtip_map
    DTT_mean = co.DTT_mean
    RTT_map = co.RTT_map
    time_points = co.time_points
    enhanced_trees = co.enhanced_trees
    gelm_id, gelm_stat  = co.gelm

    R_tree(co.tree)
    reval("""tree = temp_tree""")

    reval("""enhanced_trees = list()""")
    for i in 1:length(co.enhanced_trees)
        @rput i
        R_tree(co.enhanced_trees[i])
        reval("""enhanced_trees[[i]] = list(tree = temp_tree, rates = temp_rates)""")
    end

    @rput chains
    @rput rtt_chains
    @rput sig_map
    @rput al_map
    @rput eps_map
    @rput lambda0_map
    @rput lambdai_map
    @rput lambdatip_map
    @rput DTT_mean
    @rput RTT_map
    @rput time_points
    @rput gelm_id
    @rput gelm_stat

    @rput path
    reval("""
        CladsOutput = list(
            tree = tree,
            rtt_chains = rtt_chains,
            chains = chains,
            sig_map = sig_map,
            al_map = al_map,
            eps_map = eps_map,
            lambda0_map = lambda0_map,
            lambdai_map = lambdai_map,
            lambdatip_map = lambdatip_map,
            DTT_mean = DTT_mean,
            RTT_map = RTT_map,
            time_points = time_points,
            enhanced_trees = enhanced_trees,
            gelm = c(gelm_id, gelm_stat)
        )

        save(CladsOutput, file = path)
    """)
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA