text:
import Statistics
using AlgebraicMultigrid
import IterativeSolvers
using SuiteSparse

function StokesSolvers(mesh, Kuu, Kup, fu, fp, solver)
    if solver==0
        # Coupled solve
        zero_p = spdiagm(mesh.nel, mesh.nel)
        K      = [Kuu Kup; -Kup' zero_p]
        f      = [fu; fp]
        xh     = lu(K)\f
        Vxh    = xh[1:mesh.nf]
        Vyh    = xh[mesh.nf+1:2*mesh.nf]
        Pe     = xh[2*mesh.nf+1:end]
        Pe     = Pe .- Statistics.mean(Pe)
    elseif solver==1
        # Decoupled solve
        coef  = 1e3.*mesh.ke./mesh.vole#*ones(mesh.nel)
        # file = matopen(string(@__DIR__,"/results/matrix_ppi.mat"), "w" )
        # write(file, "coef",    coef )
        # close(file)
        Kppi  = spdiagm(coef)
        Kpu   = -Kup'
        Kuusc = Kuu .- Kup*(Kppi*Kpu)
        PC    =  0.5*(Kuusc .+ Kuusc')
        t = @elapsed Kf    = cholesky(Hermitian(PC),check = false)
        @printf("Cholesky took = %02.2e s\n", t)
        u     = zeros(2*mesh.nf,1)
        ru    = zeros(2*mesh.nf, 1)
        fusc  = zeros(2*mesh.nf,1)
        p     = zeros(mesh.nel, 1)
        rp    = zeros(mesh.nel, 1)
        # Iterations
        for rit=1:20
            ru   .= fu .- Kuu*u .- Kup*p;
            rp   .= fp .- Kpu*u;
            nrmu = norm(ru)
            nrmp = norm(rp)
            @printf("  --> Powell-Hestenes Iteration %02d\n  Momentum res.   = %2.2e\n  Continuity res. = %2.2e\n", rit, nrmu/sqrt(length(ru)), nrmp/sqrt(length(rp)))
            if nrmu/sqrt(length(ru)) < 1e-10 && nrmp/sqrt(length(ru)) < 1e-10
                break
            end
            fusc .= fu  .- Kup*(Kppi*fp .+ p)
            u    .= Kf\fusc
            p   .+= Kppi*(fp .- Kpu*u)
        end
        # Post-process solve
        Vxh = u[1:mesh.nf]
        Vyh = u[mesh.nf+1:2*mesh.nf]
        Pe  = p[:]
    elseif solver==2
        # Decoupled solve
        coef  = 1e3.*mesh.ke./mesh.vole#*ones(mesh.nel)
        # file = matopen(string(@__DIR__,"/results/matrix_ppi.mat"), "w" )
        # write(file, "coef",    coef )
        # close(file)
        Kppi  = spdiagm(coef)
        Kpu   = -Kup'
        Kuusc = Kuu - Kup*(Kppi*Kpu)
        PC    =  0.5*(Kuusc + Kuusc')
        # t = @elapsed Kf    = cholesky(Hermitian(PC),check = false)
        # @time ml = ruge_stuben(PC) #pas mal
        @time ml = smoothed_aggregation(PC) #pas mal
        @time pc = aspreconditioner(ml)
        # @printf("Cholesky took = %02.2e s\n", t)
        u     = zeros(2*mesh.nf, 1)
        ru    = zeros(2*mesh.nf, 1)
        fusc  = zeros(2*mesh.nf, 1)
        p     = zeros( mesh.nel,  1)
        rp    = zeros( mesh.nel,  1)
        # Iterations
        for rit=1:20
            ru   .= fu .- Kuu*u .- Kup*p
            rp   .= fp .- Kpu*u
            @printf("  --> Powell-Hestenes Iteration %02d\n  Momentum res.   = %2.2e\n  Continuity res. = %2.2e\n", rit, norm(ru)/sqrt(length(ru)), norm(rp)/sqrt(length(rp)))
            if norm(ru)/sqrt(length(ru)) < 1e-12 && norm(rp)/sqrt(length(ru)) < 1e-12
                break
            end
            fusc .=  fu  - Kup*(Kppi*fp + p)
            # u    .= Kf\fusc
            # @time u = IterativeSolvers.gmres(Kuusc, fusc, Pl = pc, reltol=1e-6) #, maxiter=100
            # @time IterativeSolvers.cg!(u, Kuusc, fusc, Pl = pc, reltol=1e-3)
            @time IterativeSolvers.cg!(u, PC, fusc, Pl = pc, reltol=1e-6)
            p   .+= Kppi*(fp - Kpu*u)
        end
        # Post-process solve
        Vxh = u[1:mesh.nf]
        Vyh = u[mesh.nf+1:2*mesh.nf]
        Pe  = p[:]
    elseif solver==3
        # Decoupled solve
        coef  = 1e3.*mesh.ke./mesh.vole#*ones(mesh.nel)
        Kppi  = spdiagm(coef)
        Kpu   = .- Kup'
        Kuusc = Kuu .- Kup*(Kppi*Kpu)
        ndof  = size(Kuu,1)
        ndofx = Int64(ndof/2)
        Kxx   = Kuu[1:ndofx,1:ndofx]
        t = @elapsed Kxxf  = cholesky(Hermitian(Kxx),check = false)
        @printf("Cholesky took = %02.2e s\n", t)

        # @printf("Cholesky took = %02.2e s\n", t)
        u     = zeros(Float64, 2*mesh.nf)
        ru    = zeros(Float64, 2*mesh.nf)
        fusc  = zeros(Float64, 2*mesh.nf)
        p     = zeros(Float64, mesh.nel)
        rp    = zeros(Float64, mesh.nel)
        ######################################
        restart = 30
        f      = zeros(Float64, 2*mesh.nf)
        v      = zeros(Float64, 2*mesh.nf)
        s      = zeros(Float64, 2*mesh.nf)
        val    = zeros(Float64, restart)
        VV     = zeros(Float64, (2*mesh.nf, restart) )  # !!!!!!!!!! allocate in the right sense :D
        SS     = zeros(Float64, (2*mesh.nf, restart) )
        ######################################
        # Iterations
        for rit=1:20
            ru   .= fu[:] .- Kuu*u .- Kup*p
            rp   .= fp[:] .- Kpu*u
            nrmu = sqrt(mydotavx( ru,ru ) )#norm(v)
            nrmp = sqrt(mydotavx( rp,rp ) )#norm(v)
            @printf("  --> Powell-Hestenes Iteration %02d\n  Momentum res.   = %2.2e\n  Continuity res. = %2.2e\n", rit, nrmu/sqrt(length(ru)), nrmp/sqrt(length(rp)))
            if nrmu/sqrt(length(ru)) < 1e-10 && nrmp/sqrt(length(ru)) < 1e-10
                break
            end
            fusc .=  fu[:] .- Kup*(Kppi*fp .+ p)
            @time KSP_GCR_StokesFCFV!( u, Kuusc, fusc, 1e-10, 2, Kxxf, f, v, s, val, VV, SS, restart  )
            p   .+= Kppi*(fp .- Kpu*u)
        end
        # Post-process solve
        Vxh = u[1:mesh.nf]
        Vyh = u[mesh.nf+1:2*mesh.nf]
        Pe  = p[:]
    end
    return Vxh, Vyh, Pe
end

#--------------------------------------------------------------------#

function KSP_GCR_StokesFCFV!( x::Vector{Float64}, M::SparseMatrixCSC{Float64, Int64}, b::Vector{Float64}, eps::Float64, noisy::Int64, Kxxf::SuiteSparse.CHOLMOD.Factor{Float64}, f::Vector{Float64}, v::Vector{Float64}, s::Vector{Float64}, val::Vector{Float64}, VV::Matrix{Float64}, SS::Matrix{Float64}, restart::Int64 )
    # Initialise
    val .= 0.0
    s   .= 0.0
    v   .= 0.0
    VV  .= 0.0
    SS  .= 0.0
    # KSP GCR solver
    norm_r, norm0 = 0.0, 0.0
    N               = length(x)
    maxit           = 10*restart
    ncyc, its       = 0, 0
    i1, i2, success = 0, 0, 0
    # Initial residual
    @tturbo f     .= b .- M*x
    norm_r = sqrt(mydotavx( f, f ) )#norm(v)norm(f)
    norm0  = norm_r;
    ndof   = size(M,1)
    ndofu  = Int64(ndof/2)
    # ldiv!(P::SuiteSparse.CHOLMOD.Factor{Float64}, v) = (P \ v)
    # Solving procedure
     while ( success == 0 && its<maxit )
        for i1=1:restart
            # Apply preconditioner, s = PC^{-1} f
            s[1:ndofu]     .= Kxxf \ @view f[1:ndofu]
            s[ndofu+1:end] .= Kxxf \ @view f[ndofu+1:end]
            # Action of Jacobian on s: v = J*s
            @tturbo mul!(v, M, s)
            # Approximation of the Jv product
            for i2=1:i1
                val[i2] = mydotavx( v, view(VV, :, i2 ) )
            end
            # Scaling
            for i2=1:i1
                @tturbo v .-= val[i2] .* view(VV, :, i2 )
                @tturbo s .-= val[i2] .* view(SS, :, i2 )
            end
            # -----------------
            nrm_inv = 1.0 / sqrt(mydotavx( v, v ) )
            r_dot_v = mydotavx( f, v )  * nrm_inv
            # -----------------
            @tturbo v     .*= nrm_inv
            @tturbo s     .*= nrm_inv
            # -----------------
            @tturbo x     .+= r_dot_v.*s
            @tturbo f     .-= r_dot_v.*v
            # -----------------
            norm_r  = sqrt(mydotavx( f, f ) )
            if norm_r/sqrt(length(f)) < eps || its==23
                @printf("It. %04d: res. = %2.2e\n", its, norm_r/sqrt(length(f)))
                success = 1
                println("converged")
                break
            end
            # Store
            @tturbo VV[:,i1] .= v
            @tturbo SS[:,i1] .= s
            its      += 1
        end
        its  += 1
        ncyc += 1
    end
    if (noisy>1) @printf("[%1.4d] %1.4d KSP GCR Residual %1.12e %1.12e\n", ncyc, its, norm_r, norm_r/norm0); end
    return its
end
export KSP_GCR_Stokes!

#--------------------------------------------------------------------#

function mydotavx(A, B)
    s = zero(promote_type(eltype(A), eltype(B)))
    @tturbo for i in eachindex(A,B)
        s += A[i] * B[i]
    end
    s
end


text:
module BAI

include("general.jl")
include("reducers.jl")
include("response_oracles.jl")

end

ext: jl

max_stars_repo_name: tomaskocak/BAI

max_issues_repo_name: tomaskocak/BAI

max_forks_repo_name: tomaskocak/BAI


text:
@doc doc"""
Default length (bytes) upto which to pre-fetch block metadata.
(10 blocks of default size)
""" ->
const HDFS_READER_WINDOW_LENGTH = UInt64(64*1024*1024*10)

@doc doc"""
# HDFSFileReader
Provides Julia IO APIs for reading HDFS files.
Communicates with namenode for file metadata (through HDFSClient)
and to datanodes for file data (through HDFSBlockReader)
""" ->
type HDFSFileReader <: IO
    client::HDFSClient
    path::AbstractString
    size::UInt64
    block_sz::UInt64
    crc_check::Bool                         # whether to check CRC with file operations
    buffer::IOBuffer                        # to hold excess data read from HDFSBlockReader (should use a bucket brigade)
    fptr::UInt64                            # absolute offset in file we start (0:filesize. filesize => eof)
    blocks::Nullable{LocatedBlocksProto}    # block metadata around the current fptr (if connected)
    blk_reader::Nullable{HDFSBlockReader}   # current block reader (if connected)

    function HDFSFileReader(client::HDFSClient, path::AbstractString, offset::UInt64=zero(UInt64), check_crc::Bool=false)
        path = abspath(client, path)

        nfile_status = _get_file_info(client, path)
        isnull(nfile_status) && throw(HDFSException("File not found: $path"))
        file_status = get(nfile_status)
        (file_status.fileType == HdfsFileStatusProto_FileType.IS_FILE) || throw(HDFSException("Not a regular file: $path"))

        size = file_status.length
        block_sz = file_status.blocksize
        new(client, path, size, block_sz, check_crc, PipeBuffer(), offset, Nullable{LocatedBlocksProto}(), Nullable{HDFSBlockReader}())
    end
end

function URI(reader::HDFSFileReader, show_offset::Bool)
    ch = reader.client.nn_conn.channel
    user = username(ch.ugi)
    user_spec = isempty(user) ? user : "$(user)@"
    offset = show_offset ? "#$(reader.fptr)" : ""
    URI("hdfs://$(user_spec)$(ch.host):$(ch.port)$(reader.path)$(offset)")
end

show(io::IO, reader::HDFSFileReader) = println(io, "HDFSFileReader: $(URI(reader, true))")

#
# File pointer
eof(reader::HDFSFileReader) = (reader.fptr+1 > reader.size)
position(reader::HDFSFileReader) = reader.fptr
filesize(reader::HDFSFileReader) = reader.size
nb_available(reader::HDFSFileReader) = filesize(reader) - position(reader)

function seek(reader::HDFSFileReader, n::Integer)
    (n == reader.fptr) && return
    (0 <= n <= filesize(reader)) || throw(HDFSException("invalid seek position $n"))
    reader.fptr = n
    # TODO: do not disconnect if new fptr is still within the reach
    if isconnected(reader)
        reuse = false
        if !isnull(reader.blk_reader)
            blk_reader = get(reader.blk_reader)
            reuse = eof(blk_reader)
        end
        disconnect(reader, reuse)
    end
    nothing
end
seekend(reader::HDFSFileReader) = seek(reader, filesize(reader))
seekstart(reader::HDFSFileReader) = seek(reader, 0)
skip(reader::HDFSFileReader, n::Integer) = seek(reader, UInt64(n+position(reader)))

_block_has_offset(block::LocatedBlockProto, offset::UInt64, block_sz::UInt64) = (block.offset <= offset < (block.offset + block_sz))

function _find_block(blocks::LocatedBlocksProto, offset::UInt64, block_sz::UInt64)
    if isfilled(blocks, :blocks)
        blockarr = blocks.blocks
        for block in blockarr
            if _block_has_offset(block, offset, block_sz)
                block.corrupt && throw(HDFSException("Corrupt block found at offset $(block.offset)"))
                return Nullable(block)
            end
        end
    end
    Nullable{LocatedBlockProto}()
end

function _current_window(reader::HDFSFileReader)
    isnull(reader.blocks) && (return 0:-1)
    blocks = get(reader.blocks)

    block = blocks.blocks[end]
    window_max = block.offset + block.b.numBytes
    window_min = blocks.blocks[1].offset
    window_min:(window_max-1)
end

isconnected(reader::HDFSFileReader) = !isnull(reader.blk_reader)
function disconnect(reader::HDFSFileReader, reuse::Bool)
    @logmsg("disconnecting $(URI(reader,true)) reuse:$reuse")
    if !isnull(reader.blocks)
        reader.blocks = Nullable{LocatedBlocksProto}()
    end
    if !isnull(reader.blk_reader)
        blk_reader = get(reader.blk_reader)
        isconnected(blk_reader) && disconnect(blk_reader, reuse)
        reader.blk_reader = Nullable{HDFSBlockReader}()
    end
    reader.buffer = PipeBuffer()
    nothing
end

function connect(reader::HDFSFileReader)
    @logmsg("connect if required $(URI(reader,true))")
    eof(reader) && return
    isconnected(reader) && return

    @logmsg("connecting $(URI(reader,true))")
    if !(reader.fptr in _current_window(reader))
        @logmsg("moving block window for  $(URI(reader,true))")
        len = min(nb_available(reader), HDFS_READER_WINDOW_LENGTH)
        nblocks = _get_block_locations(reader.client, reader.path, reader.fptr, len)
        isnull(nblocks) && throw(HDFSException("Could not get block metadata for $(URI(reader,true))"))
        blocks = get(nblocks)
        reader.blocks = blocks
    else
        blocks = get(reader.blocks)
    end

    nblock = _find_block(blocks, reader.fptr, reader.block_sz)
    isnull(nblock) && throw(HDFSException("Could not get block metadata for $(URI(reader,true))"))
    block = get(nblock)
    block.corrupt && throw(HDFSException("Corrupt block found for $(URI(reader,true))"))

    nnode_info = _select_node_for_block(block)
    isnull(nnode_info) && throw(HDFSException("Could not get a valid datanode for $(URI(reader,true))"))
    node = get(nnode_info).id

    block_len = block.b.numBytes
    offset_into_block = (reader.fptr % reader.block_sz)
    len = block_len - offset_into_block
    reader.blk_reader = Nullable(HDFSBlockReader(node.ipAddr, node.xferPort, block, offset_into_block, len, reader.crc_check))
    nothing
end

function _read_and_buffer(reader::HDFSFileReader, out::Array{UInt8,1}, offset::UInt64, len::UInt64)
    requested::UInt64 = len
    try
        # while data remains to be copied
        while len > 0
            # set/reset block reader if required
            if !isnull(reader.blk_reader)
                blk_reader = get(reader.blk_reader)
                if eof(blk_reader)
                    disconnect(blk_reader, true)
                    reader.blk_reader = Nullable{HDFSBlockReader}()
                    connect(reader)
                end
            else
                connect(reader)
            end
            blk_reader = get(reader.blk_reader)

            # read packet from reader
            ret = read_packet!(blk_reader, out, offset)                     # try to read directly into output
            if ret < 0
                pkt_len = len + UInt64(abs(ret))                            # bytes in this packet
                buff = Array(UInt8, pkt_len)                                # allocate a temporary array
                @logmsg("allocated temporary array of size $pkt_len, len:$len, ret:$ret, offset:$offset, bufflen:$(length(buff)), outlen:$(length(out))")
                ret = read_packet!(blk_reader, buff, UInt64(1))             # read complete packet
                copy!(out, offset, buff, 1, len)                            # copy len bytes to output
                Base.write_sub(reader.buffer, buff, len+1, pkt_len-len)     # copy remaining data to buffer
                copylen = len
                len = 0                                                     # we definitely do not need to read any more
            else
                copylen = len - ret                                         # packet was read directly into output
                len = ret                                                   # we may still have to read more data
            end
            offset += copylen
            reader.fptr += copylen
        end
    catch ex
        channel = reader.channel
        @logmsg("exception receiving from $(channel.host):$(channel.port): $ex")
        disconnect(reader, false)
        rethrow(ex)
    end
    requested
end

function close(reader::HDFSFileReader)
    reuse = false
    if !isnull(reader.blk_reader)
        blk_reader = get(reader.blk_reader)
        reuse = eof(blk_reader)
    end
    disconnect(reader, reuse)
end

#
# File read
function read!(reader::HDFSFileReader, a::Vector{UInt8})
    remaining::UInt64 = length(a)*sizeof(UInt8)
    offset::UInt64 = 1
    # while not eof and remaining to be filled
    while !eof(reader) && (remaining > 0)
        if (reader.fptr == reader.size)
            reader.fptr += 1
            break
        end
        navlb = nb_available(reader.buffer)
        # if buffer empty, fill buffer
        if navlb == 0
            if (reader.fptr + remaining) > reader.size
                canread = reader.size - reader.fptr
                tb = Array(UInt8, Int(canread/sizeof(UInt8)))
                nbytes = _read_and_buffer(reader, tb, UInt64(1), canread)
                copy!(a, offset, tb, 1, length(tb))
            else
                nbytes = _read_and_buffer(reader, a, offset, remaining)
            end
        else
            nbytes = min(navlb, remaining)
            @logmsg("reading $nbytes from buffer. navlb:$navlb remaining:$remaining, offset:$offset")
            Base.read_sub(reader.buffer, a, offset, nbytes)
            reader.fptr += nbytes
        end
        # fill from buffer
        (nbytes == 0) && break
        remaining -= nbytes
        offset += nbytes
    end
    (remaining > 0) && throw(EOFError())
    a
end
const _a = Array(UInt8, 1)
read(reader::HDFSFileReader, ::Type{UInt8}) = (read!(reader, _a); _a[1])
readbytes(reader::HDFSFileReader, nb::Integer) = read!(reader, Array(UInt8, nb))
readall(reader::HDFSFileReader) = readbytes(reader, nb_available(reader))

@doc doc"""
# HDFSFileWriter
Provides Julia IO APIs for writing HDFS files.

Steps:
- get datanodes to write onto from namenode
- break bytes into packets
- write packets into first datanode's data channel (datanode will mirror packets to downstream datanodes)
- when complete, confirm namenode with a blockReceived

On close call NameNode.complete to:
- remove lease
    lease is a write lock for file modification. no leases are required for reading files.
    leases are managed by namenode
    lease is added on file create or append
    dfs client should start thread to renew leases periodically
- change file from under construction to complete
""" ->
type HDFSFileWriter <: IO
    client::HDFSClient
    path::AbstractString
    fptr::UInt64
    blk::Nullable{LocatedBlockProto}
    blk_writer::Nullable{HDFSBlockWriter}

    function HDFSFileWriter(client::HDFSClient, path::AbstractString)
        path = abspath(client, path)
        path_exists = exists(client, path)
        fs = _create_file(client, path, path_exists, zero(UInt32), zero(UInt64), DEFAULT_FILE_MODE, false)
        isnull(fs) && throw(HDFSException("Error openeing $path for write"))
        new(client, path, 0, Nullable{LocatedBlockProto}(), Nullable{HDFSBlockWriter}())
    end
end

function URI(writer::HDFSFileWriter, show_offset::Bool)
    ch = writer.client.nn_conn.channel
    user = username(ch.ugi)
    user_spec = isempty(user) ? user : "$(user)@"
    offset = show_offset ? "#$(writer.fptr)" : ""
    URI("hdfs://$(user_spec)$(ch.host):$(ch.port)$(writer.path)$(offset)")
end

show(io::IO, writer::HDFSFileWriter) = println(io, "HDFSFileWriter: $(URI(writer, true))")

function renewlease(writer::HDFSFileWriter)
    renewlease(writer.client)
end

function _write{T<:Union{UInt8,Vector{UInt8}}}(writer::HDFSFileWriter, data::T)
    rem_data = data
    L = rem_len = length(data)

    while rem_len > 0
        if isnull(writer.blk_writer)
            blk = _add_block(LocatedBlockProto, writer.client, writer.path, writer.blk)
            blk_writer = HDFSBlockWriter(blk, _get_server_defaults(writer.client))
            writer.blk_writer = Nullable(blk_writer)
            writer.blk = Nullable(blk)
        else
            blk_writer = get(writer.blk_writer)
        end

        rem_len -= write(blk_writer, rem_data)
        if rem_len > 0
            flush(blk_writer)
            disconnect(blk_writer, true)
            writer.blk_writer = Nullable{HDFSBlockWriter}()
            (T <: UInt8) || (rem_data = unsafe_wrap(Array, pointer(data, L-rem_len+1), rem_len))
        end
    end
    L
end
write(writer::HDFSFileWriter, data::UInt8) = _write(writer, data)
write(writer::HDFSFileWriter, data::Vector{UInt8}) = _write(writer, data)

function close(writer::HDFSFileWriter)
    if !isnull(writer.blk_writer)
        blk_writer = get(writer.blk_writer)
        flush(blk_writer)
        disconnect(blk_writer, true)
        _complete_file(writer.client, writer.path, Nullable(blk_writer.block.b))
        writer.blk_writer = Nullable{HDFSBlockWriter}()
    end
    nothing
end

#
# File open
function open(client::HDFSClient, path::AbstractString, rd::Bool, wr::Bool, ff::Bool;
                replication::UInt32=zero(UInt32), blocksz::UInt64=zero(UInt64), mode::UInt32=DEFAULT_FILE_MODE,
                offset::UInt64=zero(UInt64), crc::Bool=false)
    if rd
        return HDFSFileReader(client, path, offset, crc)
    elseif wr
        if ff
            error("Not implemented")
        else
            return HDFSFileWriter(client, path)
        end
    else
        throw(ArgumentError("invalid open mode. read:$rd, write:$wr, append:$ff"))
    end
end
open(file::HDFSFile, rd::Bool, wr::Bool, ff::Bool; opts...) = open(file.client, file.path, rd, wr, ff; opts...)
open(client::HDFSClient, path::AbstractString; opts...) = open(client, path, true, false, false; opts...)
open(file::HDFSFile; opts...) = open(file.client, file.path; opts...)

function open(client::HDFSClient, path::AbstractString, open_mode::AbstractString; opts...)
    open_mode == "r"  ? open(client, path, true , false, false; opts...) :
    open_mode == "w"  ? open(client, path, false, true , false; opts...) :
    open_mode == "a"  ? open(client, path, false, true , true;  opts...) :
    throw(ArgumentError("invalid open mode: $mode"))
end
open(file::HDFSFile, open_mode::AbstractString; opts...) = open(file.client, file.path, open_mode; opts...)

#
# File copy
function cp(frompath::Union{HDFSFile,AbstractString}, topath::Union{HDFSFile,AbstractString}; offset::UInt64=zero(UInt64), len::UInt64=zero(UInt64), crc::Bool=false)
    if isa(frompath, HDFSFile)
        fromfile = open(frompath.client, frompath.path, "r"; offset=offset, crc=crc)
    else
        fromfile = open(frompath, "r")
        (offset > 0) && skip(fromfile, offset)
    end

    if isa(topath, HDFSFile)
        tofile = open(topath.client, topath.path, "w")
    else
        tofile = open(topath, "w")
    end

    buff_sz = 64*1024*1024
    buff = Array(UInt8, buff_sz)
    brem = btot = (len == 0) ? (filesize(fromfile)-offset) : len
    while brem > 0
        bread = min(brem, buff_sz)
        (bread == length(buff)) || resize!(buff, bread)
        read!(fromfile, buff)
        write(tofile, buff)
        brem -= bread
        @logmsg("remaining $brem/$btot")
    end
    close(fromfile)
    close(tofile)
    nothing
end


text:
struct EulerCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
end

struct SplitEulerCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
end

function alg_cache(alg::SplitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  SplitEulerCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype))
end

u_cache(c::SplitEulerCache) = ()
du_cache(c::SplitEulerCache) = (c.k,c.fsalfirst)

struct SplitEulerConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::SplitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SplitEulerConstantCache()

function alg_cache(alg::Euler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  EulerCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype))
end

u_cache(c::EulerCache) = ()
du_cache(c::EulerCache) = (c.k,c.fsalfirst)

struct EulerConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Euler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = EulerConstantCache()

struct HeunCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  utilde::rateType
  fsalfirst::rateType
end

struct RalstonCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  utilde::rateType
  fsalfirst::rateType
end

u_cache(c::Union{HeunCache,RalstonCache}) = ()
du_cache(c::Union{HeunCache,RalstonCache}) = (c.k,c.fsalfirst,c.utilde)

function alg_cache(alg::Heun,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  HeunCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype),zero(rate_prototype))
end

function alg_cache(alg::Ralston,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  RalstonCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype),zero(rate_prototype))
end

struct HeunConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Heun,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = HeunConstantCache()

struct RalstonConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Ralston,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = RalstonConstantCache()

struct MidpointCache{uType,rateType,uEltypeNoUnits} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  atmp::uEltypeNoUnits
  fsalfirst::rateType
end

u_cache(c::MidpointCache) = (c.atmp,)
du_cache(c::MidpointCache) = (c.k,c.fsalfirst)

struct MidpointConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::Midpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u); atmp = similar(u, uEltypeNoUnits)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  MidpointCache(u,uprev,k,tmp,atmp,fsalfirst)
end

alg_cache(alg::Midpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = MidpointConstantCache()

struct RK4Cache{uType,rateType,uEltypeNoUnits} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k₂::rateType
  k₃::rateType
  k₄::rateType
  k::rateType
  tmp::uType
  atmp::uEltypeNoUnits
end

u_cache(c::RK4Cache) = (c.atmp,)
du_cache(c::RK4Cache) = (c.fsalfirst,c.k₂,c.k₃,c.k₄,c.k)

struct RK4ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::RK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k₁ = zero(rate_prototype)
  k₂ = zero(rate_prototype)
  k₃ = zero(rate_prototype)
  k₄ = zero(rate_prototype)
  k  = zero(rate_prototype)
  tmp = similar(u); atmp = similar(u, uEltypeNoUnits)
  RK4Cache(u,uprev,k₁,k₂,k₃,k₄,k,tmp,atmp)
end

alg_cache(alg::RK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = RK4ConstantCache()

struct CarpenterKennedy2N54Cache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  A2::Float64
  A3::Float64
  A4::Float64
  A5::Float64
  B1::Float64
  B2::Float64
  B3::Float64
  B4::Float64
  B5::Float64
  c2::Float64
  c3::Float64
  c4::Float64
  c5::Float64
end

u_cache(c::CarpenterKennedy2N54Cache) = ()
du_cache(c::CarpenterKennedy2N54Cache) = (c.k,c.fsalfirst)

struct CarpenterKennedy2N54ConstantCache <: OrdinaryDiffEqConstantCache
  A2::Float64
  A3::Float64
  A4::Float64
  A5::Float64
  B1::Float64
  B2::Float64
  B3::Float64
  B4::Float64
  B5::Float64
  c2::Float64
  c3::Float64
  c4::Float64
  c5::Float64
end

function alg_cache(alg::CarpenterKennedy2N54,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  A2 = -567301805773/1357537059087
  A3 = -2404267990393/2016746695238
  A4 = -3550918686646/2091501179385
  A5 = -1275806237668/842570457699
  B1 = 1432997174477/9575080441755
  B2 = 5161836677717/13612068292357
  B3 = 1720146321549/2090206949498
  B4 = 3134564353537/4481467310338
  B5 = 2277821191437/14882151754819
  c2 = 1432997174477/9575080441755
  c3 = 2526269341429/6820363962896
  c4 = 2006345519317/3224310063776
  c5 = 2802321613138/2924317926251
  CarpenterKennedy2N54Cache(u,uprev,k,tmp,fsalfirst,A2,A3,A4,A5,B1,B2,B3,B4,B5,c2,c3,c4,c5)
end

function alg_cache(alg::CarpenterKennedy2N54,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  A2 = -567301805773/1357537059087
  A3 = -2404267990393/2016746695238
  A4 = -3550918686646/2091501179385
  A5 = -1275806237668/842570457699
  B1 = 1432997174477/9575080441755
  B2 = 5161836677717/13612068292357
  B3 = 1720146321549/2090206949498
  B4 = 3134564353537/4481467310338
  B5 = 2277821191437/14882151754819
  c2 = 1432997174477/9575080441755
  c3 = 2526269341429/6820363962896
  c4 = 2006345519317/3224310063776
  c5 = 2802321613138/2924317926251
  CarpenterKennedy2N54ConstantCache(A2,A3,A4,A5,B1,B2,B3,B4,B5,c2,c3,c4,c5)
end

struct BS3Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::BS3Cache) = (c.atmp,c.utilde)
du_cache(c::BS3Cache) = (c.fsalfirst,c.k2,c.k3,c.k4)

function alg_cache(alg::BS3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  utilde = similar(u,axes(u))
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  BS3Cache(u,uprev,k1,k2,k3,k4,utilde,tmp,atmp,tab)
end

alg_cache(alg::BS3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct OwrenZen3Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::OwrenZen3Cache) = (c.atmp,c.utilde)
du_cache(c::OwrenZen3Cache) = (c.k1,c.k2,c.k3,c.k4)

function alg_cache(alg::OwrenZen3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  utilde = similar(u,axes(u))
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen3Cache(u,uprev,k1,k2,k3,k4,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct OwrenZen4Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::OwrenZen4Cache) = (c.atmp,c.utilde)
du_cache(c::OwrenZen4Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6)

function alg_cache(alg::OwrenZen4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  utilde = similar(u,axes(u))
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen4Cache(u,uprev,k1,k2,k3,k4,k5,k6,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct OwrenZen5Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::OwrenZen5Cache) = (c.atmp,c.utilde)
du_cache(c::OwrenZen5Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8)

function alg_cache(alg::OwrenZen5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  k8 = zero(rate_prototype)
  utilde = similar(u,axes(u))
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct BS5Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::BS5Cache) = (c.utilde,c.atmp)
du_cache(c::BS5Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.k8)

function alg_cache(alg::BS5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = BS5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  k8 = zero(rate_prototype)
  utilde = similar(u,axes(u))
  atmp = similar(u,uEltypeNoUnits,axes(u))
  tmp = similar(u)
  BS5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp,tab)
end

alg_cache(alg::BS5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = BS5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct Tsit5Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::Tsit5Cache) = (c.utilde,c.atmp)
du_cache(c::Tsit5Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7)

function alg_cache(alg::Tsit5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Tsit5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  utilde = similar(u,axes(u))
  atmp = similar(u,uEltypeNoUnits,axes(u))
  tmp = similar(u)
  Tsit5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
end

alg_cache(alg::Tsit5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Tsit5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct DP5Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  update::rateType
  bspl::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::DP5Cache) = (c.utilde,c.atmp)
du_cache(c::DP5Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.dense_tmp3,c.dense_tmp4,c.update,c.bspl)

function alg_cache(alg::DP5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = k2
  k7 = zero(rate_prototype) # This is FSAL'd to k1
  dense_tmp3 = k2
  dense_tmp4 = k5
  bspl = k3

  tmp = similar(u) # has to be separate for FSAL
  utilde = tmp

  if eltype(u) != uEltypeNoUnits || calck
    update = zero(rate_prototype)
    atmp = similar(u,uEltypeNoUnits,axes(u))
  else
    update = k7
    atmp = k3
  end

  tab = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  cache = DP5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp,tab)
  cache
end

alg_cache(alg::DP5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

struct DP5ThreadedCache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  update::rateType
  bspl::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::DP5ThreadedCache) = (c.utilde,c.atmp)
du_cache(c::DP5ThreadedCache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7,c.dense_tmp3,c.dense_tmp4,c.update,c.bspl)

function alg_cache(alg::DP5Threaded,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  dense_tmp3 = zero(rate_prototype)
  dense_tmp4 = zero(rate_prototype)
  update = zero(rate_prototype)
  bspl = zero(rate_prototype)
  utilde = similar(u,axes(u))
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits,axes(u))
  tab = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  DP5ThreadedCache(u,uprev,k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp,tab)
end

struct Anas5Cache{uType,uArrayType,rateType,uEltypeNoUnits,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  utilde::uArrayType
  tmp::uType
  atmp::uEltypeNoUnits
  tab::TabType
end

u_cache(c::Anas5Cache) = (c.atmp,c.utilde)
du_cache(c::Anas5Cache) = (c.k1,c.k2,c.k3,c.k4,c.k5,c.k6,c.k7)

function alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  utilde = similar(u,axes(u))
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  Anas5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
end

alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))


text:
# Find All Approximate Occurrences of a Pattern in a String

calc_hamming_distance(s1, s2) = count(((c1, c2),) -> c1 != c2, zip(s1, s2))

function find_approximate_pattern_occurrences(text, pattern, d)
    k = length(pattern)
    matching_indices = Int64[]
    for ix in 1:(length(text)-k+1)
        substring = text[ix:ix+k-1]
        if calc_hamming_distance(substring, pattern) <= d
            push!(matching_indices, ix - 1)
        end
    end
    matching_indices
end

function main()
    pattern = readline()
    text = readline()
    d = parse(Int64, readline())
    result = find_approximate_pattern_occurrences(text, pattern, d)
    println(join(result, " "))
end

main()


text:
#Runs a simple Metropolis-Hastings (ie MCMC) algoritm to simulate an
#exponential distribution, which has the probability density
#p(x)=exp(-x/m^2), where m>0.
#
#Author: H. Paul Keeler, 2019.
# Website: hpaulkeeler.com
# Repository: github.com/hpaulkeeler/posts

#Note: Need the .+ for adding a scalar to an array
#Also need . for sqrt, exp, cos, sin etc and assinging scalars to arrays

#clearconsole(); #for clearing Julia REPL console

using Distributions #for random simulations
using PyPlot #for plotting
using StatsBase #for histograms etc
using Random
using LinearAlgebra
PyPlot.close("all");  # close all PyPlot figures

#set random seed for reproducibility
#Random.seed!(1234)

numbSim=10^4; #number of random variables simulated
numbSteps=25; #number of steps for the Markov process
numbBins=50; #number of bins for histogram

sigma=1; #standard deviation for normal random steps
m=2; #parameter (ie mean) for distribution to be simulated

function fun_p(x)
    return ((exp.(-x./m)./m).*(x.>0));
end

xRand=rand(numbSim); #random intial values
probCurrent=fun_p(xRand); #current transition probabilities

for jj=1:numbSteps
    zRand= xRand.+sigma.*rand(Normal(),numbSim);#take a (normally distributed) random step
    #zRand= xRand +2*sigma*(rand(size(xRand))-0.5);#take a (uniformly distributed) random step
    probProposal=fun_p(zRand); #proposed probability

    #acceptance rejection step
    booleAccept=rand(numbSim) .< probProposal./probCurrent;
    #update state of random walk/Markov chain
    xRand[booleAccept]=zRand[booleAccept];
    #update transition probabilities
    probCurrent[booleAccept]=probProposal[booleAccept];

end

#histogram section: empirical probability density
histX=fit(Histogram, xRand,nbins=numbBins); #find histogram data
histX=normalize(histX,mode=:pdf); #normalize histogram
binEdges=histX.edges; #retrieve bin edges
xValues=(binEdges[1][2:end]+binEdges[1][1:end-1])./2; #mid-points of bins
pdfEmp=(histX.weights)
#analytic solution of probability density
pdfExact=fun_p(xValues);

# Plotting
PyPlot.plot(xValues, pdfExact)
PyPlot.scatter(xValues, pdfEmp, marker="x", c="r");
PyPlot.grid(true);
PyPlot.xlabel("x");
PyPlot.ylabel("p(x)");
PyPlot.show();

max_stars_repo_name: hpkeeler/posts

max_issues_repo_name: hpkeeler/posts

max_forks_repo_name: hpkeeler/posts


text:
"""
    orgGermanGT()
---
This is to manipulate the German genotype data, which were received on Apr. 04, 2020.
The data consist four platfomrs:
1. LD:  12 ID, ignored
2. MD:   5 ID, ignored, these are actually of platform v3, with a lot missings
3. V2:  70 ID, takes 73+ seconds to finish, will find a better way later.
4. V3: 733 ID, took 765s !!! too bad

# ToDo:
* A faster C++ version/improved Julia version, I/O adoped here is too slow
* Errror check need to be included
"""
function orgGermanGT()
    cd(work_dir)
    empty_dir("tmp")
    warning("Ignoring data in data/genotypes/german/ld")
    warning("Ignoring data in data/genotypes/german/md")

    dir = "data/genotypes/german"
    mdr = "data/maps/updated"
    den = "data/genotypes/step-1.plk"
    isdir(den) || mkdir(den)

    for ver in ["v2", "v3"]
        title("Dealing with the German data in $dir/$ver")
        list = readdir("$dir/$ver")
        fr2ped("$dir/$ver", list, "tmp/plink.ped", "Top")

        item("Creating a SNP dictionary with $mdr/$ver.map")
        dic = ref_map_dict("$mdr/$ver.map")
        done()

        item("Merge 'tmp/plink.ped' and 'tmp/plink.map' to '$den/german-$ver'")
        sample = list[1]
        create_plink_map("$dir/$ver/$sample", dic, "tmp/plink.map")
        ped_n_map_to_bed("tmp/plink.ped", "tmp/plink.map", "$den/german-$ver")
        done()
    end
end

"""
    orgDutchGT()
---
This is to manipulae the Dutch data, which were received much earlier.
I am doing them with `Julia` procedure.

# Simple description
1. platform 10690, 1994 ID
2. platform 10993,  102 ID
3. platform 11483,  109 ID
4. platform    v2,  343 ID
5. platform    v3,   44 ID
6. platform  777k,   10 ID
"""
function orgDutchGT()
    cd(work_dir)
    empty_dir("tmp")

    dir = "data/genotypes/dutch"
    mdr = "data/maps/updated"
    den = "data/genotypes/step-1.plk"
    isdir(den) || mkdir(den)

    for ver in ["d1", "d2", "d3", "v2", "v3", "v7"]
        title("Dealing with the Dutch data in $dir/$ver")
        list = readdir("$dir/$ver")
        fr2ped("$dir/$ver", list, "tmp/plink.ped", "Top")

        item("Creating a SNP dictionary with $mdr/$ver.map")
        dic = ref_map_dict("$mdr/$ver.map")
        done()

        item("Merge 'tmp/plink.ped' and 'tmp/plink.map' to '$den/dutch-$ver'")
        sample = list[1]
        create_plink_map("$dir/$ver/$sample", dic, "tmp/plink.map")
        ped_n_map_to_bed("tmp/plink.ped", "tmp/plink.map", "$den/dutch-$ver")
        done()
    end
end

"""
    autosome_subset(list)
---
Take subset of `bed` in `data/genotypes/step-1.plk/{list}`.
The subset derived consists the shared SNP between `data/maps/target.snp` and `bed`.
"""
function autosome_subset(list)
    cd(work_dir)
    empty_dir("tmp")
    title("Extract autosomal subset")
    item("Create the target set")
    target = begin
        snps = String[]
        p = r"^[1-9][0-9]*$"
        for line in eachline("data/maps/updated/v3.map")
            snp, chr, bp = split(line)
            match(p, chr) === nothing || push!(snps, snp)
        end
        Set(snps)
    end
    done()

    sdir = "data/genotypes/step-1.plk"
    tdir = "data/genotypes/step-2.plk"
    isdir(tdir) || mkdir(tdir)
    for plk in list
        item("Dealing with '$sdir/$plk'")
        open("tmp/sub.snp", "w") do io
            for line in eachline("$sdir/$plk.bim")
                snp = split(line)[2]
                ∈(snp, target) && write(io, snp, '\n')
            end
        end
        bed_snp_subset("$sdir/$plk", "tmp/sub.snp", "$tdir/$plk")
        done()
    end
end


text:
# Generate some test data
t = collect(0:0.1:5.0)
f_(t) = [sin(t); cos(t)]
df_(t) = [cos(t); -sin(t)]
X = hcat(map(i -> f_(i), t)...)
Y = randn(2,2)*X
DX = hcat(map(i -> df_(i), t)...)
p = [2.7]
u_(u, p, t) = [p[1] * u[1] - 3.0 * t]
U = hcat(map(i -> u_(X[:, i], p, t[i]), 1:length(t))...)

@testset "DiscreteProblem" begin
    p1 = DiscreteDataDrivenProblem(X)
    p2 = DiscreteDataDrivenProblem(X, t)
    p3 = DiscreteDataDrivenProblem(X, t, U)
    p4 = DiscreteDataDrivenProblem(X, t, u_, p = p)

    for p_ in [p1;p2;p3;p4]
        @test is_valid(p_)
        @test is_discrete(p_)
    end
    @test isequal(p3.U, p4.U)
    @test isequal(p4.p, p)
end

@testset "ContinuousProblem" begin
    p2 = ContinuousDataDrivenProblem(X, t, GaussianKernel())
    p3 = ContinuousDataDrivenProblem(X, DX)
    p4 = ContinuousDataDrivenProblem(X, t, U)
    p5 = ContinuousDataDrivenProblem(X, t, DX, U)
    p6 = ContinuousDataDrivenProblem(X, t, DX, u_, p = p)

    for p_ in [p2;p3;p4;p5;p6]
        @test is_valid(p_)
        @test is_continuous(p_)
    end

    @test isapprox(p3.X, p2.X, atol = 0.3)
    @test isapprox(p3.DX, p2.DX, atol = 0.3)
    @test isequal(p5.U, p6.U)
    @test is_parametrized(p6)
    @test isequal(p6.p, p)
    @test has_timepoints(p4)
end

@testset "DirectProblem" begin
    p1 = DirectDataDrivenProblem(X, Y)
    p2 = DirectDataDrivenProblem(X, t, Y)
    p3 = DirectDataDrivenProblem(X, t, Y, U)
    p4 = DirectDataDrivenProblem(X, t, Y, u_, p = p)

    for p_ in [p1;p2;p3;p4]
        @test is_valid(p_)
        @test is_direct(p_)
    end

    @test isequal(p3.U, p4.U)
    @test isequal(p4.p, p)
    @test !is_autonomous(p3)
end


text:
@testset "Properties" begin
    # the expected base methods for arrays also work on tensors
    t1 = Tensor(rand(3, 3))
    @test size(t1) == (3, 3)
    @test size(t1, 1) == 3
    @test ndims(t1) == 2
    @test length(t1) == 9
    @test eltype(t1) === Float64

    t2 = Tensor(rand(Int32, 3, 3, 3))
    @test size(t2) == (3, 3, 3)
    @test size(t2, 2) == 3
    @test ndims(t2) == 3
    @test length(t2) == 27
    @test eltype(t2) === Int32

    t3 = Tensor(3.0f0)
    @test size(t3) == ()
    @test size(t3, 3) == 1
    @test ndims(t3) == 0
    @test length(t3) == 1
    @test eltype(t3) === Float32
end

@testset "Equality" begin
    a1 = rand(2, 3, 3)
    t1 = Tensor(a1)

    # t1 and t2 share their underlying data, thus t1 === t2
    a2 = a1
    t2 = Tensor(a2)
    @test t1 === t2
    @test t1 == t2
    @test t1 ≈ t2

    # t1 and t2 don't share their data, but contain the same values
    a2 = copy(a1)
    t2 = Tensor(a2)
    @test t1 !== t2
    @test t1 == t2
    @test t1 ≈ t2

    # t1 and t2 differ very slightly in one value
    a2[1, 2, 3] += 1e-10
    t2 = Tensor(a2)
    @test t1 !== t2
    @test t1 != t2
    @test t1 ≈ t2

    # t1 and t2 differ in one value
    a2[1, 2, 3] = 4
    t2 = Tensor(a2)
    @test t1 !== t2
    @test t1 != t2
    @test !(t1 ≈ t2)
end

@testset "Types" begin
    # rand puts out Float64 per default
    t1 = Tensor(rand(3, 3))
    @test eltype(t1) === Float64

    t2 = Tensor{Float32}(t1)
    @test eltype(t2) === Float32

    t3 = convert(Tensor{Float16}, t1)
    @test eltype(t3) === Float16

    # like Base.convert throws an error when performing inexact conversions
    @test_throws InexactError Tensor{Int64}(t1)
    @test_throws InexactError convert(Tensor{Int64}, t1)

    # scalar tensors can be converted to the underlying scalar
    t4 = Tensor(4)
    @test convert(Int64, t4) == 4
    @test Int64(t4) == 4

    # not possible for non-scalar tensors
    @test_throws MethodError convert(Float64, t1)
    @test_throws MethodError Float64(t1)
end

@testset "Reshaping" begin
    t1 = Tensor(Vector(1:9))
    @test size(t1) == (9,)

    # reshaping does not alter the original tensor
    t2 = reshape(t1, 3, 3)
    @test size(t1) == (9,)
    @test size(t2) == (3, 3)
    @test t2 == Tensor([
        1 4 7
        2 5 8
        3 6 9
    ])

    # permutedims also does not alter the original tensor
    t3 = Tensor(rand(5, 4, 2))
    @test size(t3) == (5, 4, 2)
    t4 = permutedims(t3, (3, 1, 2))
    @test size(t3) == (5, 4, 2)
    @test size(t4) == (2, 5, 4)

    @test t3 != t4
    @test t3[:, :, 1] == t4[1, :, :]
    @test t3[:, :, 2] == t4[2, :, :]

    # on 1- and 2-dimensional tensors the permutation tuple can be omitted
    @test permutedims(Tensor([1, 2, 3])) == Tensor([1 2 3])
    @test permutedims(Tensor([1 2 3])) == Tensor(reshape([1, 2, 3], 3, 1))
    @test permutedims(Tensor([1 2 3; 4 5 6])) == Tensor([1 4; 2 5; 3 6])
end

@testset "Indexing" begin
    # indexing a tensor works like indexing the wrapped array
    t = Tensor(rand(1, 2, 3, 3))
    @test axes(t) == (Base.OneTo(1), Base.OneTo(2), Base.OneTo(3), Base.OneTo(3))
    @test eachindex(t) == Base.OneTo(18)
    @test CartesianIndices(t) == CartesianIndices(axes(t))

    t1 = t[1, :, :, :]
    @test t1 isa Tensor{Float64,3}
    @test size(t1) == (2, 3, 3)
    @test t1 == t[1, :, :, :]

    t2 = t[:, 2, :, 3]
    @test t2 isa Tensor{Float64,2}
    @test size(t2) == (1, 3)
    @test t2 == t[:, 2, :, 3]

    t3 = t[1, 2, 3, 2]
    @test t3 isa Tensor{Float64,0}
    @test size(t3) == ()
    @test t3 == t[1, 2, 3, 2]
    @test_throws BoundsError t3[1]

    t[1, 1, 1, :] = Tensor([1, 2, 3])
    @test t[1, 1, 1, :] == Tensor([1.0, 2.0, 3.0])

    # tensor elements can be assigned scalars and scalar tensors
    t[1, 1, 1, 2] = 4
    t[1, 1, 1, 3] = Tensor(9)
    @test t[1, 1, 1, :] == Tensor([1.0, 4.0, 9.0])

    # throws BoundsError if indices are out-of-range
    @test_throws BoundsError t[2, :, :, :]
    @test_throws BoundsError t[1, 3, 3, 3] = 2

    # checkbounds can also be used directly
    @test checkbounds(Bool, t, 1, 1, 1, 1)
    @test !checkbounds(Bool, t, 2, 1, 1, 1)
    @test_throws BoundsError checkbounds(t, 2, 1, 1, 1)
end

@testset "Printing" begin
    # test that tensors show methods work upon array show methods
    io = IOBuffer()
    # NOTE: for compatibility with Julia 1.0 which used slightly different; remove if
    # compatibility with this version is no longer needed/supported
    norm(s) = filter(x -> x != ' ', s)

    t1 = Tensor(2)
    show(io, t1)
    @test norm(String(take!(io))) == norm("Tensor{Int64, 0}(2)")
    show(io, "text/plain", t1)
    @test norm(String(take!(io))) == norm(join(
        [
            "scalar Tensor{Int64, 0}:"
            "2"
        ],
        "\n",
    ))

    t2 = Tensor(Vector{Float64}(1:4))
    show(io, t2)
    @test norm(String(take!(io))) == norm("Tensor{Float64, 1}($(t2.data))")
    show(io, "text/plain", t2)
    @test norm(String(take!(io))) == norm(join(
        [
            "4-element Tensor{Float64, 1}:"
            " 1.0"
            " 2.0"
            " 3.0"
            " 4.0"
        ],
        "\n",
    ))

    t3 = Tensor(reshape(Vector{Int32}(1:9), 3, 3))
    show(io, t3)
    @test norm(String(take!(io))) == norm("Tensor{Int32, 2}($(t3.data))")
    show(io, "text/plain", t3)
    @test norm(String(take!(io))) == norm(join(
        [
            "3×3 Tensor{Int32, 2}:"
            " 1  4  7"
            " 2  5  8"
            " 3  6  9"
        ],
        "\n",
    ))

    t4 = Tensor(reshape(Vector{Float32}(1:8), 2, 2, 2))
    show(io, t4)
    @test norm(String(take!(io))) == norm("Tensor{Float32, 3}($(t4.data))")
    show(io, "text/plain", t4)
    @test norm(String(take!(io))) == norm(
        join(
            [
                "2×2×2 Tensor{Float32, 3}:"
                "[:, :, 1] ="
                " 1.0  3.0"
                " 2.0  4.0"
                ""
                "[:, :, 2] ="
                " 5.0  7.0"
                " 6.0  8.0"
            ],
            "\n",
        ),
    )
end

@testset "Hashing" begin
    t = Tensor(rand(3, 3))

    # tensors are hashable and their hash is different than the hash of the underlying data
    @test hash(t) != hash(t.data)
    @test hash(t) == hash(hash(t.data) + hash("t"))

    # tensors can be used in a dict
    d = Dict(t => (3, 3))
    @test d[t] == (3, 3)
end

@testset "Delta" begin
    # the delta tensor is a square boolean tensor of dimension 2
    d3 = δ(3)
    @test size(d3) == (3, 3)
    @test eltype(d3) === Bool

    # the delta tensor acts like the neutral element of the inner product
    t = Tensor(rand(3, 4, 5))
    @test t ⋅ δ(5) == δ(3) ⋅ t == t
    @test δ(4) ⋅ δ(4) == δ(4)

    # delta tensor can also be created using a value type
    @test δ(Val(4)) == δ(4)
end

@testset "Epsilon" begin
    # the nth epsilon tensor is an n×...×n Int8 tensor of dimension n
    eps3 = ϵ(3)
    @test size(eps3) == (3, 3, 3)
    @test eltype(eps3) === Int8

    # the epsilon tensor can be used to generalize the cross product to tensors
    t1 = Tensor([1, 0, 0])
    t2 = Tensor([0, 1, 0])
    # TODO: adjust cross product definition if actual implementation changes
    @test -(t1 ⋅ ϵ(3) ⋅ t2) == Tensor([0, 0, 1]) # = t1 × t2

    # it also generalizes the tensors of higher dimensions than 1
    t3 = Tensor([
        1 0 0
        2 0 0
    ])
    @test -(t3 ⋅ ϵ(3) ⋅ t2) == Tensor([
        0 0 1
        0 0 2
    ])

    # it can also work on dimensions different than 3 with a little creativity
    t4 = Tensor([1, 0, 0, 0])
    t5 = Tensor([0, 1, 0, 0])
    t6 = Tensor([0, 0, 1, 0])
    @test ϵ(4) ⋅ t4 ⋅ t5 ⋅ t6 == Tensor([0, 0, 0, 1])
end


text:
# Getis-Ord statistic
struct GetisOrd <: AbstractLocalSpatialAutocorrelation
    n::Int
    G::Vector{Float64}
    p::Vector{Float64}
    Gperms::Matrix{Float64}
    Gpermsmean::Vector{Float64}
    Gpermsstd::Vector{Float64}
    z::Vector{Float64}
    q::Vector{Symbol}
    star::Bool
end

"""
    getisord(x, W)
Compute the Getis-Ord statistic.

# Optional Arguments
- `star=true`: compute the Gi* statistic, or the Gi if set to `false`.
- `permutations=9999`: number of permutations for the randomization test.
- `rng=default_rng()`: random number generator for the randomization test.
"""
function getisord(x::AbstractVector{T} where T<:Real, W::SpatialWeights; permutations::Int = 9999,
    star::Bool = true, rng::AbstractRNG = default_rng())::GetisOrd

    wt = wtransformation(W)
    wt == :row || wt == :binary || throw(ArgumentError("W must be row standardized or binary"))

    n = length(x)

    denon::Float64 = sum(x)

    # Auxiliar function to calculate Getis Ord
    function getisord_calc(xi::Number, wi::AbstractVector, xneighi::AbstractVector)::Float64
        return sum( wi .* xneighi) ./ (denon - xi)
    end

    function getisord_calc_star(xi::Number, wi::AbstractVector, xneighi::AbstractVector)::Float64
        if wt == :row
            wistar = 1 / (length(wi) + 1)
        elseif wt == :binary
            wistar = 1
        end
        return (wistar .* xi + sum(wistar .* xneighi)) ./ (denon)
    end

    if star
        getisord_calc_fun = getisord_calc_star
    else
        getisord_calc_fun = getisord_calc
    end

    # Getis-Ord
    G = zeros(n)
    for i in 1:n
        xi = x[i]
        wi = weights(W, i)
        xneighi = x[neighbors(W, i)]
        G[i] = getisord_calc_fun(xi, wi, xneighi)
    end

    # Conditional randomizatoin
    Gperms = crand_local(permutations, x, W, getisord_calc_fun, rng)

    larger = sum(Gperms .>= repeat(G, 1, permutations), dims = 2)
    low = (permutations .- larger) .< larger
    larger[low] .= permutations .- larger[low]
    p = (larger .+ 1) ./ (permutations + 1)
    p = vec(p)

    Gpermsstd = vec(std(Gperms, dims = 2, corrected = false))
    Gpermsmean = vec(mean(Gperms, dims = 2))
    zval = (G .- Gpermsmean)  ./ Gpermsstd
    zval = vec(zval)

    # Classification
    q = Array{Symbol}(undef, n)
    for i in 1:n
        if zval[i] > 0
            q[i] = :H
        else
            q[i] = :L
        end
    end

    return GetisOrd(n, G, p, Gperms, Gpermsmean, Gpermsstd, zval, q, star)

end

score(x::GetisOrd) = x.G;

scoreperms(x::GetisOrd) = x.Gperms;

mean(x::GetisOrd) = x.Gpermsmean;

std(x::GetisOrd) = x.Gpermsstd;

zscore(x::GetisOrd) = x.z;

pvalue(x::GetisOrd) = x.p;

assignments(x::GetisOrd) = x.q;

testname(x::GetisOrd) = x.star ? "Getis-Ord Gi*" : "Getis-Ord Gi";

function labelsorder(::GetisOrd)
    return [:ns, :H, :L]
end

function labelsnames(::GetisOrd)
    return ["Not Significant"; "High"; "Low"]
end

function labelcolor(::GetisOrd, x::String)

    if x == "H"
        catcolor = :red
        alpha = 1.0
    elseif x == "L"
        catcolor =:blue
        alpha = 1.0
    # not significant
    elseif x == "ns"
        catcolor = :lightgrey
        alpha = 0.4
    end

    return (catcolor, alpha)
end


text:
using JuLIP, JuLIP.Potentials, TightBinding
using JuLIP.Potentials: site_energy, site_energy_d
TB = TightBinding

@testset "Site Energy" begin

# test parameters
beta = 10.0        # temperature / smearing paramter: 10 to 50 for room temperature
n0 = 1            # site index where we compute the site energy
NQUAD = (4, 8, 16, 20)     # number of contour points
DIM = (1,2,3)

# define the model
tbm = TB.NRLTB.NRLTBModel(:Si, TB.FermiDiracSmearing(beta), bzquad = TB.GammaPoint() )

# check that the site energies add up to the total
at = bulk(:Si, pbc = (true, false, false), cubic=true) * DIM
TB.set_δNel!(tbm, at)
E = energy(tbm, at)
∑En = sum( site_energy(tbm, at, n) for n = 1:length(at) )
println("Testing that the spectral-decompositions site energy sums to total energy")
@show E - ∑En
@test abs.(E - ∑En) < 1e-10

# now the real system to test on
calc = TB.PEXSI(tbm, 5, [1])
at = DIM * bulk(:Si, pbc=(false,false,false), cubic=true)
JuLIP.rattle!(at, 0.02)

# calibrate the PEXSI calculator on a mini-system
print("calibrating . . . ")
TB.calibrate!(calc, at; at_train = bulk(:Si, pbc=true), npoles = 8)
println("done."); @test true

# output some useful info if we are watching the tests...
@show length(at)
@show TB.get_EminEmax(at)

# compute the site energy the old way and compare against the PEXSI calculation
Eold = TB.site_energy(tbm, at, n0)
println("Old Site Energy (via spectral decomposition): ", Eold)
println("Testing convergence of PEXSI site energy")
Enew = 0.0
for nquad in NQUAD
   TB.set_npoles!(calc, nquad)
   Enew = site_energy(calc, at, n0)
   println("nquad = ", nquad, "; error = ", abs(Enew - Eold))
end
@test abs(Enew - Eold) < 1e-5

println("Test consistency of site forces")
TB.set_npoles!(calc, 8)
X = positions(at) |> mat
Es = site_energy(calc, at, n0)
dEs = site_energy_d(calc, at, n0) |> mat
dEsh = []
errors = Float64[]
for p = 2:9
   h = 0.1^p
   dEsh = zeros(dEs)
   for n = 1:length(X)
      X[n] += h
      set_positions!(at, X)
      Esh = site_energy(calc, at, n0)
      dEsh[n] = (Esh - Es) / h
      X[n] -= h
   end
   println( " ", p, " | ", vecnorm(dEs-dEsh, Inf) )
   push!(errors, vecnorm(dEs-dEsh, Inf))
end
@test minimum(errors) < 1e-3 * maximum(errors)

println("Test consistency of ContourCalculator for multiple sites")
Is = unique(mod.(rand(Int, length(at) ÷ 3), length(at)) + 1)
Eold = sum( site_energy(tbm, at, n0) for n0 in Is )
Enew = 0.0
for nquad in NQUAD
   calc.nquad = nquad
   Enew = TB.partial_energy(calc, at, Is)
   println("nquad = ", nquad, "; rel-error = ", abs(Enew - Eold) / abs(Eold))
end
@test  abs(Enew - Eold) / abs(Eold) < 1e-5

println("Test consistency of multiple site forces")
calc.nquad = 8
X = copy( positions(at) |> mat )
Is = unique(mod.(rand(Int, length(at) ÷ 3), length(at)) + 1)
Es, dEs = TB.pexsi_partial_energy(calc, at, Is, true)
dEs = dEs |> mat
dEsh = []
errors = []
println(" p  |  error ")
for p = 2:9
   h = 0.1^p
   dEsh = zeros(dEs)
   for n = 1:length(X)
      X[n] += h
      set_positions!(at, X)
      Esh = TB.partial_energy(calc, at, Is)
      dEsh[n] = (Esh - Es) / h
      X[n] -= h
   end
   println( " ", p, " | ", vecnorm(dEs-dEsh, Inf) )
   push!(errors, vecnorm(dEs-dEsh, Inf))
end
@test minimum(errors) < 1e-3 * maximum(errors)

end


text:
#!/usr/bin/env julia

function fact(n::Int)::Int
    n >= 0 || error("fact argument must be positive")
    n <= 1 && return 1
    return n*fact(n - 1)
end

for n in 5:-1:-1
    print("fact($n) = ")
    println(fact(n))
end


text:
using DataTables
using CSV
using DecisionTree
# using MLBase

cd("/Users/pstey/projects_code/string_matching/")

include("tf_idf_utils.jl")

dt = CSV.read("AUCAdjudicatedCategories.csv", DataTable)

const STOPWORDS = [
    "the",
    "is",
    "at",
    "which",
    "on",
    "in",
    "for",
    "with"]

###
# Data Cleaning and Prepartion
###

# Recode a couple of minor errors
recode!(dt, :Category, Dict("Aortic Disease" => "aortic", "hypotension" => "acute"), String)

# Testing our TF/IDF functions
a = ["this is a string of words", "this is a sentence!", "this is it", "what about this?"]
tf_idf_datatable(a, STOPWORDS)

# Adding and ID column, and then make it the first
# column in the datatable.
dt[:id] = 1:(nrow(dt))
dt = dt[:, [3, 1, 2]]

# Get dataframe with TF/IDF values for each indication
dt_tf_idf = tf_idf_datatable(dt[:Indication], STOPWORDS)

###
# Fitting Random Forest
###

# Cast data to numeric arrays
labels = convert(Array{String,1}, dt[:Category])
features = convert(Array, dt_tf_idf)

## Split data into training and test set
function train_test_split(dat, pct_train; seed = rand(Int, 1))
    srand(seed)
    n = size(dat, 1)
    n_train = round(Int, pct_train * n)
    train_indcs = sample(1:n, n_train, replace = false)
    test_indcs = setdiff(1:n, train_indcs)
    return (train_indcs, test_indcs)
end

train, test = train_test_split(dt, 0.80, seed = 137)

labels_trn = labels[train]
labels_tst = labels[test]

features_trn = features[train, :]
features_tst = features[test, :]

# Run n-fold cross validation for forests using
# `mtry` random features, 100 trees, and 5 folds.
mtry = 250
ntrees = 1000
kfolds = 5

nfoldCV_forest(labels_trn, features_trn, mtry, ntrees, kfolds, 0.7)

# Build forest with meta-parameters we like from CV above
maxlabels = 5
fm1 = build_forest(labels_trn, features_trn, mtry, ntrees, maxlabels, 0.7)

# Get predicted category
yhat1 = apply_forest(fm1, features_tst)

# Write predictions to file
writecsv("/Users/pstey/Desktop/preds_labels.csv", hcat(yhat1, labels_tst))

# Get measures of model performance
fm1_acc = mean(yhat1 .== labels_tst)
cm1 = confusion_matrix(labels_tst, yhat1)
cm1.kappa


text:
using Test

include("dist.jl")
include("solving.jl")
include("scripts.jl")


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

const OUT_LEN = 32
const KEY_LEN = 32
const BLOCK_LEN = 64
const CHUNK_LEN = 1024

const CHUNK_START = UInt32(1) << 0
const CHUNK_END   = UInt32(1) << 1
const PARENT      = UInt32(1) << 2
const ROOT        = UInt32(1) << 3
const KEYED_HASH  = UInt32(1) << 4
const DERIVE_KEY_CONTEXT  = UInt32(1) << 5
const DERIVE_KEY_MATERIAL = UInt32(1) << 6

const IV = SVector{8, UInt32}(
    UInt32(0x6A09E667), UInt32(0xBB67AE85), UInt32(0x3C6EF372), UInt32(0xA54FF53A),
    UInt32(0x510E527F), UInt32(0x9B05688C), UInt32(0x1F83D9AB), UInt32(0x5BE0CD19)
)