text:
# Submodule to compute the coefficients of shearlet recovery and decoding in 2D

##############################################################################
# tools for doing efficient Fourier domain operations

"""
    unshearing!(X, neededShear, P, coeffs, padBy, j) where {N,T}
    unshearing!(X::AbstractArray{Complex{T},2}, neededShear, P,
                coeffs::AbstractArray{T,3}, padBy, j) where {N,T}
given a set of shearlet coefficients `coeffs`, add the inverse of the `j`th one
to X
"""
function unshearing!(X̂, neededShear, P, coeffs, padBy, j)
    for i in eachindex(view(X̂, 1, 1, axes(X̂)[3:end]...))
        cFreq = P * (pad(coeffs[:, :, j, i], padBy))
        X̂[:, :, i] = X̂[:, :, i] + cFreq .* neededShear
    end
end

function unshearing!(X̂::AbstractArray{S,2}, neededShear, P,
                     coeffs::AbstractArray{T,3}, padBy, j) where {T,
                                                                  S <: Complex}
    cFreq = P * (pad(coeffs[:, :, j], padBy))
    X̂[:, :] = X̂[:, :] + cFreq .* neededShear
end

"""
    baseDomain!(X, X̂, duals, P, used1, used2)
Multiply by the dual coordinates and convert to the space domain
"""
function baseDomain!(X::AbstractArray{S}, X̂::AbstractArray{T,N}, duals, P,
                     used1, used2) where {S<:Complex, T, N}
    for i in eachindex(view(X, 1, 1, axes(X)[3:end]...))
        X[:, :, i] = (P \ (X̂[:,:, i] ./ (duals)))[used1, used2]
    end
end

function baseDomain!(X::AbstractArray{S}, X̂::AbstractArray{T,2}, duals, P,
                     used1, used2) where {T, S<:Complex}
    X[:, :] = (P \ (X̂[:,:] ./ (duals)))[used1, used2]
end

function baseDomain!(X::AbstractArray{S}, X̂::AbstractArray{T,N}, duals, P,
                     used1, used2) where {S<:Real, T, N}
    for i in eachindex(view(X, 1, 1, axes(X)[3:end]...))
        X[:, :, i] = real.((P \ (X̂[:,:, i] ./ (duals)))[used1, used2])
    end
end

function baseDomain!(X::AbstractArray{S}, X̂::AbstractArray{T,2}, duals, P,
                     used1, used2) where {T, S<:Real}
    X[:, :] = real.((P \ (X̂[:,:] ./ (duals)))[used1, used2])
end

"""
    shearing!(X::AbstractArray{T,N}, neededShear, P, coeffs, padBy, used1,
                   used2, j) where {N,T}
    shearing!(X::AbstractArray{T,2}, neededShear, P, coeffs, padBy, used1,
                   used2, j) where T
fill in the `j`th coefficient and fill it into the matrix of coefficients
"""
function shearing!(X::AbstractArray{T,N}, neededShear, P, coeffs, padBy, used1,
                   used2, j) where {N,T}
    for i in eachindex(view(X, 1, 1, axes(X)[3:end]...))
        Xfreq = P * pad(X[:, :, i], padBy)
        tmpCoeff = (P \ (Xfreq .* neededShear))[used1, used2]
        if T <: Real
            coeffs[:, :, j, i] = real.(tmpCoeff)
        else
            coeffs[:, :, j, i] = tmpCoeff
        end
    end
end

function shearing!(X::AbstractArray{T,2}, neededShear, P, coeffs, padBy, used1,
                   used2, j) where T
    Xfreq =  P * pad(X[:, :], padBy)
    tmpCoeff = (P \ (Xfreq .* neededShear))[used1, used2]
    if T <: Real
        coeffs[:, :, j] = real.(tmpCoeff)
    else
        coeffs[:, :, j] = tmpCoeff
    end
end

function getAPlan(P, zeroMat, CT, padBy)
    if (typeof(P) <: Real) && (CT <: Real)
        P = plan_rfft(pad(zeroMat, padBy), flags = FFTW.PATIENT)
    elseif (typeof(P) <: Real) && (CT <: Complex)
        P = plan_fft!(pad(zeroMat, padBy), flags =
                      FFTW.PATIENT)
    end
    return P
end

##############################################################################
## Function that compute the coefficient matrix of the Shearlet Transform of
## some array
"""
    sheardec2D(X::AbstractArray{CT,N},
                    shearletSystem::Shearlab.Shearletsystem2D{T, CT};
                    P::{plan_(r)fft}=-1) where {T <: Number, CT, N}

Compute the coefficient matrix of the Shearlet transform across the first two
dimensions of X. If `shearletSystem.padded` is true, it pads the shearlet
coefficients by `padBy`, which defaults to the widest spatial component of any
of the shearlets (computed via `getPadBy`). If you're doing a lot of these, you
can precompute the Fourier transform plan and put it in `P` (note that it
should have the same element type as your data), and should have the size of X
*after* padding.

"""
function sheardec2D(X::AbstractArray{CT, N},
                    shearletSystem::Shearlab.Shearletsystem2D{T, CT}; P=-1) where {T <: Number, CT, N}
    @assert N >= 2

    padBy = shearletSystem.padBy

    if N>2
        coeffShape = (size(X)[1:2]..., size(shearletSystem.shearlets, 3),
                      size(X)[3:end]...)
    else
        coeffShape = (size(X)[1:2]..., size(shearletSystem.shearlets, 3))
    end

    coeffs = zeros(CT, coeffShape...)

    # if you didn't hand it a plan, adjust the kind made based on the inputs
    # data Type
    P = getAPlan(P, zeros(CT, size(X)[1:2]), CT, padBy)

    # compute shearlet coefficients at each scale
    # note that pointwise multiplication in the fourier domain equals
    # convolution in the time-domain
    used1 = padBy[1] .+ (1:size(X, 1))
    used2 = padBy[2] .+ (1:size(X, 2))
    for j = 1:shearletSystem.nShearlets
        # The fourier transform of X
        if j==1
            j = shearletSystem.nShearlets
        elseif j == shearletSystem.nShearlets
            j = 1
        end
        neededShear = shearletSystem.shearlets[:, :, j]
        shearing!(X, neededShear, P,  coeffs, padBy, used1, used2, j)
    end
    return coeffs
end # sheardec2D
⋆(shears::Shearletsystem2D, X) = sheardec(X, shears) #\star
⋆(X, shears::Shearletsystem2D) = sheardec(X, shears) #\star
import Base./
import Base.\
\(shears::Shearletsystem2D, X) = shearrec2D(X, shears) #\star
/(X, shears::Shearletsystem2D) = shearrecD2(X, shears) #\star

##############################################################################
## Function that recovers the image with the Shearlet transform with some shearletSystem
"""
    shearrec2D(coeffs::AbstractArray{complex{T},N},
               shearletSystem::Shearlab.Shearletsystem2D{T}; P=-1,
               padBy::Tuple{Int, Int} = getPadBy(shearletSystem)) where {T<:Real, N, B, C, D}

function that recovers the image with the Shearlet transform with some
shearletSystem. If you're doing a lot of these, you can precompute the Fourier
transform plan and put it in `P` (note that it is type sensitive), and should
have the size of X *after* padding.
"""

function shearrec2D(coeffs::AbstractArray{CT, N},
                    shearletSystem::Shearlab.Shearletsystem2D{T, CT}; P=-1) where {T<:Real, CT<:Union{T, Complex{T}}, N, B, C, D}

    padBy = shearletSystem.padBy
    # are there meta dimensions? The storage needs to adapt
    if N > 2
        X̂ = zeros(Complex{T}, size(shearletSystem.shearlets)[1:2]..., # this is X\^
                  size(coeffs)[4:end]...)
    else
        X̂ = zeros(Complex{T}, size(shearletSystem.shearlets)[1:2]...)
    end

    P = getAPlan(P, zeros(CT, shearletSystem.size), CT, padBy)

    # sum X in the Fourier domain over the coefficients at each scale/shearing
    used1 = padBy[1] .+ (1:shearletSystem.size[1])
    used2 = padBy[2] .+ (1:shearletSystem.size[2])
    for j = 1:shearletSystem.nShearlets
        if j==1
            j = shearletSystem.nShearlets
        elseif j== shearletSystem.nShearlets
            j = 1
        end

        neededShear = shearletSystem.shearlets[:, :, j]
        unshearing!(X̂, neededShear, P,  coeffs, padBy, j)
    end
    X = zeros(CT, size(coeffs)[1:end-1])
    baseDomain!(X, X̂, shearletSystem.dualFrameWeights, P, used1, used2)
    return X
end

##############################################################################
## Function that compute the array of the adjoint Shearlet Transform of
## some coefficients matrix
"""
...
    sheardecadjoint2D(coeffs,shearletSystem) compute the adjoint of the decomposition operator
...
"""
function sheardecadjoint2D(coeffs,shearletSystem)
		# Initialize reconstructed data
    gpu = shearletSystem.gpu;
    if gpu == 1
        X = AFArray(zeros(Complex{Float32},size(coeffs,1),size(coeffs,2)));
    else
        X = zeros(Complex{Float64},size(coeffs,1),size(coeffs,2));
    end
    for j = 1:shearletSystem.nShearlets
        X = X+fft(coeffs[:,:,j]).*shearletSystem.shearlets[:,:,j];
    end
    return real.(ifft((1 ./shearletSystem.dualFrameWeights).*X))

end # sheardecadjoint2D

##########################################################################################
## Function that compute the coefficient matrix of the adjoint inverse Shearlet Transform of
## some array
"""
...
    shearrecadjoint2D(X,shearletSystem)

compute the coefficient matrix of the adjoint inverse Shearlet transform of the array X
...
"""
function shearrecadjoint2D(X,shearletSystem)
    #Read the GPU info of the system
    gpu = shearletSystem.gpu;
    if gpu == 1
        coeffs = AFArray(zeros(Complex{Float32},size(shearletSystem.shearlets)));
    else
        coeffs = zeros(Complex{Float64},size(shearletSystem.shearlets));
    end
    # The fourier transform of X
    Xfreq = fft(X);
    #compute shearlet coefficients at each scale
    #not that pointwise multiplication in the fourier domain equals convolution
    #in the time-domain
    for j = 1:shearletSystem.nShearlets
        coeffs[:,:,j] = ifft(Xfreq.*shearletSystem.shearlets[:,:,j]);
    end
    return coeffs
end # shearrecadjoint2D


text:
struct PolygonDetection{T} <: AbstractObjectDetection{T}
    vertices::Vector{Point2{T}}
    core::ObjectDetectionCore

    function PolygonDetection(vertices::Vector{Point2{T}}, core_args...) where T
        if length(vertices) < 3
            throw(ArgumentError("Cannot create a polygon with less than 3 vertices."))
        end
        core = ObjectDetectionCore(core_args...)
        new{T}(vertices, core)
    end
end

vertices(detection::PolygonDetection) :: Vector{Point2} = detection.vertices

function centroid(detection::PolygonDetection{T}) :: Point2{T} where T
    return reduce((+), detection.vertices) / length(detection.vertices)
end

function bounding_box_detection(detection::PolygonDetection{T}) :: BoundingBoxDetection{T} where T
    return BoundingBoxDetection(detection.vertices, detection.core)
end

ext: jl


text:
# Type system for original files and header constants

# I'm using types as a enum here, consider changing this?
"Abstract class for representing matlab code fragments"
abstract type MatlabData end
"Type for representing Matlab strings"
struct MatStr <: MatlabData end
"Type for representing Matlab integers"
struct MatInt <: MatlabData end
"Type for representing Matlab floatingpoint numbers"
struct MatFloat <: MatlabData end
"Type for representing Matlab floatingpoint numbers"
struct OEDateTime <: MatlabData end
"Type for representing Matlab floatingpoint numbers"
struct MatVersion <: MatlabData end

"Exception type to indicate a malformed data file"
struct CorruptedException <: Exception
    message::String
end
CorruptedException() = CorruptedException("")
Base.showerror(io::IO, e::CorruptedException) = print(io, "Corrupted Exception: ", e.message)

### Constants for parsing header ###
const HEADER_N_BYTES = 1024
const HEADER_DATE_REGEX = r"^(\d{1,2}-\p{L}{3}-\d{4}) ([1-2]?\d)([0-5]\d)([1-5]?\d)"
const HEADER_DATEFORMAT = Dates.DateFormat("d-u-y")
const N_HEADER_LINE = 11

"""
    OriginalHeader(io::IOStream)
Data in the header of binary OpenEphys files.

Will throw [`CorruptedException`](@ref) if header is corrupt,
not an "OpenEphys" data format, or not version 0.4 of the data format.

# Fields

**`format`** is the name of the data format.

**`version`** is the version number of the data format.

**`headerbytes`** is the number of bytes in the header.

**`description`** is a description of the header.

**`created`** is the date and time the file was created.

**`channel`** is the name of the channel used to acquire this data.

**`channeltype`** is the type of channel used to acquire this data.

**`samplerate`** is the sampling rate in Hz.

**`blocklength`** is the length in bytes of each block of data within the file.

**`buffersize`** is the size of the buffer used during acquisition, in bytes.

**`bitvolts`** are the Volts per ADC bit.
"""
struct OriginalHeader
    "Data format"
    format::String
    "Version of data format"
    version::VersionNumber
    "Number of bytes in the header"
    headerbytes::Int
    "Description of the header"
    description::String
    "Time file created"
    created::DateTime
    "Channel name"
    channel::String
    "Channel type"
    channeltype::String
    "Sample rate for file"
    samplerate::Int
    "Length of data blocks in bytes"
    blocklength::Int
    "Size of buffer in bytes"
    buffersize::Int
    "Volts/bit of ADC values"
    bitvolts::Float64

    function OriginalHeader(
        format::String,
        version::VersionNumber,
        headerbytes::Int,
        description::String,
        created::DateTime,
        channel::String,
        channeltype::String,
        samplerate::Int,
        blocklength::Int,
        buffersize::Int,
        bitvolts::Float64
    )
        format == "Open Ephys Data Format" || throw(CorruptedException("Header is malformed"))
        version == v"0.4" || version == v"0.2" || throw(CorruptedException("Header is malformed"))
        return new(
            format,
            version,
            headerbytes,
            description,
            created,
            channel,
            channeltype,
            samplerate,
            blocklength,
            buffersize,
            bitvolts
        )
    end
end

"""
    OriginalHeader(io::IOStream)
Reads the header of the open binary file `io`. Assumes that the stream
is at the beginning of the file.
"""
function OriginalHeader(io::IOStream)
    # Read the header from the IOStream and separate on semicolons
    head = read(io, HEADER_N_BYTES)
    length(head) == HEADER_N_BYTES || throw(CorruptedException("Header is malformed"))
    headstr = transcode(String, head)
    isvalid(headstr) || throw(CorruptedException("Header is malformed"))
    substrs =  split(headstr, ';', keepempty = false)
    resize!(substrs, N_HEADER_LINE)
    format = parseline(MatStr, substrs[1])
    version = parseline(MatVersion, substrs[2])
    headerbytes = parseline(MatInt, substrs[3])
    description = parseline(MatStr, substrs[4])
    created = parseline(OEDateTime, substrs[5])
    channel = parseline(MatStr, substrs[6])
    channeltype = parseline(MatStr, substrs[7])
    samplerate = parseline(MatInt, substrs[8])
    blocklength = parseline(MatInt, substrs[9])
    buffersize = parseline(MatInt, substrs[10])
    bitvolts = parseline(MatFloat, substrs[11])
    return OriginalHeader(
        format,
        version,
        headerbytes,
        description,
        created,
        channel,
        channeltype,
        samplerate,
        blocklength,
        buffersize,
        bitvolts
    )
end

"Parse a line of Matlab source code"
function parseline end
function parseline(::Type{M}, str::AbstractString) where {M<:MatlabData}
    parseto(parsetarget(M), matread(M, str))
end

parsetarget(::Type{MatStr}) = String
parsetarget(::Type{MatInt}) = Int
parsetarget(::Type{MatFloat}) = Float64
parsetarget(::Type{OEDateTime}) = DateTime
parsetarget(::Type{MatVersion}) = VersionNumber

"Convert a string to the desired type"
function parseto end
parseto(::Type{T}, str::AbstractString) where T = parse(T, str)
function parseto(::Type{DateTime}, str::AbstractString)
    m = match(HEADER_DATE_REGEX, str)
    if m === nothing
        throw(CorruptedException("Time created is improperly formatted"))
    end
    d = DateTime(m.captures[1], HEADER_DATEFORMAT)
    local hr, mn, sc
    try
        hr = parse(Int, m.captures[2])
        mn = parse(Int, m.captures[3])
        sc = parse(Int, m.captures[4])
    catch y
        if isa(y, ArgumentError)
            throw(CorruptedException("Time created is improperly formatted"))
        else
            rethrow(y)
        end
    end
    # Check for ambiguous time
    if mapreduce(length, +, m.captures[2:4]) == 5
        if 0 <= rem(hr, 10) <= 5 && 1 <= rem(mn, 10) <= 5 # a 'shifted' parsing would also be valid
            @warn("Header time ", str, " is ambiguous! Assigning the ambiguous digit to hours.")
        end
    end
    dt = d + Dates.Hour(hr) + Dates.Minute(mn) + Dates.Second(sc)
    return dt
end
parseto(::Type{VersionNumber}, str::AbstractString) = VersionNumber(str)
parseto(::Type{String}, str::AbstractString) = String(str)
parseto(::Type{T}, str::T) where T<:AbstractString = str

"read a Matlab source line"
function matread(::Type{T}, str::S) where {T<:MatlabData, S<:AbstractString}
    m = match(rx(T), str)
    m == nothing && throw(CorruptedException("Cannot parse header"))
    string(m[1])
end

### Matlab regular expressions ###
rx(::Type{MatStr}) = r" = '(.*)'$"
rx(::Type{MatInt}) = r" = (\d*)$"
rx(::Type{MatFloat}) = r" = ([-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)$"
rx(::Type{OEDateTime}) = rx(MatStr)
rx(::Type{MatVersion}) = rx(MatFloat)

function show(io::IO, a::O) where O<:OriginalHeader
    fields = fieldnames(O)
    for field in fields
        println(io, field, ": ", getfield(a, field))
    end
end
function showcompact(io::IO, header::OriginalHeader)
    show(IOContext(io, :compact => true), "channel: $(header.channel)")
end
function show(io::IO, headers::Vector{OriginalHeader})
    for header in headers
        println(io, showcompact(header))
    end
end


text:
#
# List of elements with properties
#
struct Element
 atomic_number::Int
 element::String
 pdb_name::String
 name::String
 mass::Float64
end
const elements = [
             Element( 0      ,  "X"  ,  "X"   ,   "NotFound"    ,  0.00000   )
             Element( 1      ,  "H"  ,  "H"   ,   "Hydrogen"    ,  1.00797   )
             Element( 2      ,  "He" ,  "HEL" ,   "Helium"      ,  4.00260   )
             Element( 3      ,  "Li" ,  "LI"  ,   "Lithium"     ,  6.941     )
             Element( 4      ,  "Be" ,  "BE"  ,   "Beryllium"   ,  9.01218   )
             Element( 5      ,  "B"  ,  "B"   ,   "Boron"       ,  10.81     )
             Element( 6      ,  "C"  ,  "C"   ,   "Carbon"      ,  12.011    )
             Element( 7      ,  "N"  ,  "N"   ,   "Nitrogen"    ,  14.0067   )
             Element( 8      ,  "O"  ,  "O"   ,   "Oxygen"      ,  15.9994   )
             Element( 9      ,  "F"  ,  "F"   ,   "Fluorine"    ,  18.998403 )
             Element( 10     ,  "Ne" ,  "NEO" ,   "Neon"        ,  20.179    )
             Element( 11     ,  "Na" ,  "NA"  ,   "Sodium"      ,  22.98977  )
             Element( 11     ,  "Na" ,  "SOD" ,   "Sodium"      ,  22.98977  )
             Element( 12     ,  "Mg" ,  "MG"  ,   "Magnesium"   ,  24.305    )
             Element( 13     ,  "Al" ,  "AL"  ,   "Aluminum"    ,  26.98154  )
             Element( 14     ,  "Si" ,  "SI"  ,   "Silicon"     ,  28.0855   )
             Element( 15     ,  "P"  ,  "P"   ,   "Phosphorus"  ,  30.97376  )
             Element( 16     ,  "S"  ,  "S"   ,   "Sulfur"      ,  32.06     )
             Element( 17     ,  "Cl" ,  "CL"  ,   "Chlorine"    ,  35.453    )
             Element( 17     ,  "Cl" ,  "CLA" ,   "Chlorine"    ,  35.453    )
             Element( 18     ,  "Ar" ,  "AR"  ,   "Argon"       ,  39.948    )
             Element( 19     ,  "K"  ,  "K"   ,   "Potassium"   ,  39.0983   )
             Element( 19     ,  "K"  ,  "POT" ,   "Potassium"   ,  39.0983   )
             Element( 20     ,  "Ca" ,  "CAL" ,   "Calcium"     ,  40.08     )
             Element( 21     ,  "Sc" ,  "SC"  ,   "Scandium"    ,  44.9559   )
             Element( 22     ,  "Ti" ,  "TI"  ,   "Titanium"    ,  47.90     )
             Element( 23     ,  "V"  ,  "V"   ,   "Vanadium"    ,  50.9415   )
             Element( 24     ,  "Cr" ,  "CR"  ,   "Chromium"    ,  51.996    )
             Element( 25     ,  "Mn" ,  "MN"  ,   "Manganese"   ,  54.9380   )
             Element( 26     ,  "Fe" ,  "FE"  ,   "Iron"        ,  55.847    )
             Element( 27     ,  "Co" ,  "CO"  ,   "Cobalt"      ,  58.9332   )
             Element( 28     ,  "Ni" ,  "NI"  ,   "Nickel"      ,  58.70     )
             Element( 29     ,  "Cu" ,  "CU"  ,   "Copper"      ,  63.546    )
             Element( 30     ,  "Zn" ,  "ZN"  ,   "Zinc"        ,  65.38     )
             Element( 31     ,  "Ga" ,  "GA"  ,   "Gallium"     ,  69.72     )
             Element( 32     ,  "Ge" ,  "GE"  ,   "Germanium"   ,  72.59     )
             Element( 33     ,  "As" ,  "AS"  ,   "Arsenic"     ,  74.9216   )
             Element( 34     ,  "Se" ,  "SE"  ,   "Selenium"    ,  78.96     )
             Element( 35     ,  "Br" ,  "BR"  ,   "Bromine"     ,  79.904    )
             Element( 36     ,  "Kr" ,  "KR"  ,   "Krypton"     ,  83.80     )
             Element( 37     ,  "Rb" ,  "RB"  ,   "Rubidium"    ,  85.4678   )
             Element( 38     ,  "Sr" ,  "SR"  ,   "Strontium"   ,  87.62     )
             Element( 39     ,  "Y"  ,  "Y"   ,   "Yttrium"     ,  88.9059   )
             Element( 40     ,  "Zr" ,  "ZR"  ,   "Zirconium"   ,  91.22     )
             Element( 41     ,  "Nb" ,  "NB"  ,   "Niobium"     ,  92.9064   )
             Element( 42     ,  "Mo" ,  "MO"  ,   "Molybdenum"  ,  95.94     )
             Element( 43     ,  "Tc" ,  "TC"  ,   "Technetium"  ,  98.       )
             Element( 44     ,  "Ru" ,  "RU"  ,   "Ruthenium"   ,  101.07    )
             Element( 45     ,  "Rh" ,  "RH"  ,   "Rhodium"     ,  102.9055  )
             Element( 46     ,  "Pd" ,  "PD"  ,   "Palladium"   ,  106.4     )
             Element( 47     ,  "Ag" ,  "AG"  ,   "Silver"      ,  107.868   )
             Element( 48     ,  "Cd" ,  "CAD" ,   "Cadmium"     ,  112.41    )
             Element( 49     ,  "In" ,  "IN"  ,   "Indium"      ,  114.82    )
             Element( 50     ,  "Sn" ,  "SN"  ,   "Tin"         ,  118.69    )
             Element( 51     ,  "Sb" ,  "SB"  ,   "Antimony"    ,  121.75    )
             Element( 52     ,  "Te" ,  "TE"  ,   "Tellurium"   ,  127.60    )
             Element( 53     ,  "I"  ,  "I"   ,   "Iodine"      ,  126.9045  )
             Element( 54     ,  "Xe" ,  "XE"  ,   "Xenon"       ,  131.30    )
             Element( 55     ,  "Cs" ,  "CES" ,   "Cesium"      ,  132.9054  )
             Element( 56     ,  "Ba" ,  "BA"  ,   "Barium"      ,  137.33    )
             Element( 57     ,  "La" ,  "LA"  ,   "Lanthanum"   ,  138.9055  )
             Element( 58     ,  "Ce" ,  "CER" ,   "Cerium"      ,  140.12    )
             Element( 59     ,  "Pr" ,  "PR"  ,   "Praseodymium",  140.9077  )
             Element( 60     ,  "Nd" ,  "ND"  ,   "Neodymium"   ,  144.24    )
             Element( 61     ,  "Pm" ,  "PM"  ,   "Promethium"  ,  145       )
             Element( 62     ,  "Sm" ,  "SM"  ,   "Samarium"    ,  150.4     )
             Element( 63     ,  "Eu" ,  "EU"  ,   "Europium"    ,  151.96    )
             Element( 64     ,  "Gd" ,  "GD"  ,   "Gadolinium"  ,  157.25    )
             Element( 65     ,  "Tb" ,  "TB"  ,   "Terbium"     ,  158.9254  )
             Element( 66     ,  "Dy" ,  "DY"  ,   "Dysprosium"  ,  162.50    )
             Element( 67     ,  "Ho" ,  "HLM" ,   "Holmium"     ,  164.9304  )
             Element( 68     ,  "Er" ,  "ER"  ,   "Erbium"      ,  167.26    )
             Element( 69     ,  "Tm" ,  "TM"  ,   "Thulium"     ,  168.9342  )
             Element( 70     ,  "Yb" ,  "YB"  ,   "Ytterbium"   ,  173.04    )
             Element( 71     ,  "Lu" ,  "LU"  ,   "Lutetium"    ,  174.967   )
             Element( 72     ,  "Hf" ,  "HAF" ,   "Hafnium"     ,  178.49    )
             Element( 73     ,  "Ta" ,  "TA"  ,   "Tantalum"    ,  180.9479  )
             Element( 74     ,  "W"  ,  "W"   ,   "Tungsten"    ,  183.85    )
             Element( 75     ,  "Re" ,  "RE"  ,   "Rhenium"     ,  186.207   )
             Element( 76     ,  "Os" ,  "OS"  ,   "Osmium"      ,  190.2     )
             Element( 77     ,  "Ir" ,  "IR"  ,   "Iridium"     ,  192.22    )
             Element( 78     ,  "Pt" ,  "PT"  ,   "Platinum"    ,  195.09    )
             Element( 79     ,  "Au" ,  "AU"  ,   "Gold"        ,  196.9665  )
             Element( 80     ,  "Hg" ,  "MER" ,   "Mercury"     ,  200.59    )
             Element( 81     ,  "Tl" ,  "TL"  ,   "Thallium"    ,  204.37    )
             Element( 82     ,  "Pb" ,  "PB"  ,   "Lead"        ,  207.2     )
             Element( 83     ,  "Bi" ,  "BI"  ,   "Bismuth"     ,  208.9804  )
             Element( 84     ,  "Po" ,  "PO"  ,   "Polonium"    ,  209       )
             Element( 85     ,  "At" ,  "AT"  ,   "Astatine"    ,  210       )
             Element( 86     ,  "Rn" ,  "RN"  ,   "Radon"       ,  222       )
             Element( 87     ,  "Fr" ,  "FR"  ,   "Francium"    ,  223       )
             Element( 88     ,  "Ra" ,  "RA"  ,   "Radium"      ,  226.0254  )
             Element( 89     ,  "Ac" ,  "AC"  ,   "Actinium"    ,  227.0278  )
             Element( 90     ,  "Th" ,  "TH"  ,   "Thorium"     ,  232.0381  )
             Element( 91     ,  "Pa" ,  "PA"  ,   "Protactinium",  231.0359  )
             Element( 92     ,  "U"  ,  "U"   ,   "Uranium"     ,  238.029   )
           ]

#
# Retrive index of element in elements list from name. Returns 1 (element "X" of list) if not found
#
function element_index(name::String)
  # Try to find if there is any exact match
  i = findfirst(el -> (el.element == name ||
                       el.name == name ||
                       el.pdb_name == name ), elements)
  (i != nothing) && return i
  #
  # Try to find by PDB name, note that NT2 and 2NT2 must match N, for example
  #
  # First, check if the first char is a number
  i0 = try
    parse(Int,name[1:1])
    2
  catch
    1
  end
  # Now check if the start of this name matches some PDB code
  for (iel,el) in pairs(elements)
    lpdb = length(el.pdb_name)
    # If the name is shorter than the PDB name, no match
    (length(name) - i0 + 1) < lpdb && continue
    if name[i0:i0+lpdb-1] == el.pdb_name
      return iel
    end
  end
  return 1
end


text:
using Spid
using CSV
using Test
using BioSequences
using FASTX

open("test.pileup") do pileup_f
    open("test.pileup.ref.fa") do ref_f
        open("test.pileup.out.fa") do out_f
            ref_fasta = FASTA.Reader(ref_f)
            @test pileup_consensus_fasta(
                pileup_f, ref_fasta, 10, .9, typemax(Int)) == [
                    rec for rec in FASTA.Reader(out_f)]
        end
    end
end

rm("scratch", force=true, recursive=true)
mkdir("scratch")

merge_alignments(
    "scratch/test", ["Sample1.fa", "Sample2.fa", "Sample3.fa"],
    Regex("(?:.*/)?([^/]*)\\.fa")
)

align_assembly("MT-human.fa.gz", "MT-orang.fa",
               "scratch/aln.bam", "scratch/aln.fa"; preset="asm20")

align_assembly("MT-human.fa", "MT-orang.fa",
               "scratch/aln2.bam", "scratch/aln2.fa"; preset="asm20")

for fname in ["test.fa", "test.core.fa",
              "test.variantsOnly.fa", "test.core.variantsOnly.fa",
              "test.fa.pairwise_diffs.csv", "test.core.fa.pairwise_diffs.csv",
              "aln.fa", "aln2.fa"]
    open(fname) do f1
        open("scratch/$fname") do f2
            @test read(f1, String) == read(f2, String)
        end
    end
end

# test empty core genome doesn't crash
merge_alignments(
    "scratch/test", ["Sample1.fa", "Sample2.fa", "Sample3.fa", "SampleN.fa"],
    Regex("(?:.*/)?([^/]*)\\.fa")
)


text:
module BinaryBuilderBase

using Pkg, Pkg.Artifacts, Random, Libdl
using Base.BinaryPlatforms
using Downloads
using JSON, OutputCollectors

# Re-export useful stuff from Base.BinaryPlatforms:
export HostPlatform, platform_dlext, valid_dl_path, arch, libc,
       libgfortran_version, libstdcxx_version, cxxstring_abi,
       detect_libgfortran_version, detect_libstdcxx_version, detect_cxxstring_abi,
       call_abi, wordsize, triplet, select_platform, platforms_match,
       Platform, AnyPlatform

export AbstractSource, AbstractDependency, SetupSource, PatchSource,
    resolve_jlls, coerce_dependency, coerce_source, Runner,
    generate_compiler_wrappers!, preferred_runner, CompilerShard, UserNSRunner,
    DockerRunner, choose_shards, exeext, preferred_libgfortran_version,
    preferred_cxxstring_abi, gcc_version, available_gcc_builds, getversion,
    getpkg, aatriplet,
    nbits, proc_family, storage_dir, extract_kwargs, extract_fields,
    download_source, setup_workspace, setup_dependencies, update_registry,
    getname, cleanup_dependencies, compress_dir, prepare_for_deletion,
    run_interactive, sourcify, dependencify, with_logfile, get_concrete_platform

include("compat.jl")

include("ArchiveUtils.jl")
include("Sources.jl")
include("Dependencies.jl")
include("Prefix.jl")
include("Products.jl")
include("Platforms.jl")

include("Runner.jl")
include("Rootfs.jl")
include("squashfs_utils.jl")
include("UserNSRunner.jl")
include("DockerRunner.jl")

include("utils.jl")

# This is the location that all binary builder-related files are stored under.
# downloads, unpacked .tar.gz shards, mounted shards, ccache cache, etc....
function storage_dir(args::AbstractString...)
    global storage_cache
    dir = joinpath(storage_cache, args...)
    mkpath(dirname(dir))
    return dir
end
ccache_dir() = storage_dir("ccache")

"""
    enable_apple_file()

Return the path to file that, if exists, indicates that the user accepts to
download macOS SDK.  The file is automatically created when the package is
loaded if the environment variable `BINARYBUILDER_AUTOMATIC_APPLE` is set to
`"true"`.
"""
enable_apple_file() = storage_dir("enable_apple")

# These globals store important information such as where we're downloading
# the rootfs to, and where we're unpacking it.  These constants are initialized
# by `__init__()` to allow for environment variable overrides from the user.
storage_cache = ""
use_squashfs = false
allow_ecryptfs = false
use_ccache = false
bootstrap_list = Symbol[]

function __init__()
    global runner_override, use_squashfs, allow_ecryptfs
    global use_ccache, storage_cache

    # Allow the user to override the default value for `storage_dir`
    storage_cache = get(ENV, "BINARYBUILDER_STORAGE_DIR",
                        abspath(joinpath(@__DIR__, "..", "deps")))

    # If the user has signalled that they really want us to automatically
    # accept apple EULAs, do that.
    if get(ENV, "BINARYBUILDER_AUTOMATIC_APPLE", "") == "true"
        touch(enable_apple_file())
    end

    # If the user has overridden our runner selection algorithms, honor that
    runner_override = lowercase(get(ENV, "BINARYBUILDER_RUNNER", ""))
    if runner_override == "unprivileged"
        runner_override = "userns"
    end
    if !(runner_override in ["", "userns", "privileged", "docker"])
        @warn("Invalid runner value $runner_override, ignoring...")
        runner_override = ""
    end

    # If the user has asked for squashfs mounting instead of tarball mounting,
    # use that here.  Note that on Travis, we default to using squashfs, unless
    # BINARYBUILDER_USE_SQUASHFS is set to "false", which overrides this
    # default. If we are not on Travis, we default to using tarballs and not
    # squashfs images as using them requires `sudo` access.
    if get(ENV, "BINARYBUILDER_USE_SQUASHFS", "") == "false"
        use_squashfs = false
    elseif get(ENV, "BINARYBUILDER_USE_SQUASHFS", "") == "true"
        use_squashfs = true
    else
        # If it hasn't been specified, but we're on Travis, default to "on"
        if get(ENV, "TRAVIS", "") == "true"
            use_squashfs = true
        end

        # If it hasn't been specified but we're going to use the docker runner,
        # then set `use_squashfs` to `true` here.
        if preferred_runner() == DockerRunner
            # Conversely, if we're dock'ing it up, don't use it.
            use_squashfs = false
        elseif runner_override == "privileged"
            # If we're forcing a privileged runner, go ahead and default to squashfs
            use_squashfs = true
        end
    end

    # If the user has signified that they want to allow mounting of ecryptfs
    # paths, then let them do so at their own peril.
    if get(ENV, "BINARYBUILDER_ALLOW_ECRYPTFS", "") == "true"
        allow_ecryptfs = true
    end

    # If the user has enabled `ccache` support, use it!
    if get(ENV, "BINARYBUILDER_USE_CCACHE", "false") == "true"
        use_ccache = true
    end
end

end # module BinaryBuilderBase


text:
#==============================================================================
    Generate DataFrame of Column Names and Type from a DataFrame

    author:     mahiki@users.noreply.github.com
=============================================================================#

using DataFrames, Dates

"""
    printschema(df::DataFrame)

Return a DataFrame of ordered columns and column type for input dataFrame df

# Examples
```jldoctest
julia> df = DataFrame(a = 1:2, b = [1.0, π], c = [Date("20191002", "yyyymmdd"), Date("21120101", "yyyymmdd")], d = [1//2, missing]);
julia> printschema(df)
4×2 DataFrame
 Row │ variable  value
     │ String    Type
─────┼───────────────────────────────────────────
   1 │ a         Int64
   2 │ b         Float64
   3 │ c         Date
   4 │ d         Union{Missing, Rational{Int64}}
```
"""
function printschema(df::DataFrame)
    dfcols = mapcols(eltype, df)
    dfcolnames = names(dfcols)
    stack(dfcols, dfcolnames)
end

max_stars_repo_name: mahiki/DesertIslandDisk

max_issues_repo_name: mahiki/DesertIslandDisk

max_forks_repo_name: mahiki/DesertIslandDisk


text:
using HDF5
using DataFrames
# Arguments:
# ARGS[1] = option
# ARGS[2] = table store_sales path e.g."/home/whassan/tmp/csv/q25/store_sales_sanitized.csv"
# ARGS[5] = table web_sales path e.g. "/home/whassan/tmp/csv/q25/web_sales_sanitized.csv"

function example_dataset(ss_customer_sk, ss_ticket_number, ss_sold_date_sk, ss_net_paid,
       ws_bill_customer_sk, ws_order_number, ws_sold_date_sk, ws_net_paid)
    # ss_customer_sk   ss_ticket_number  ss_sold_date_sk  ss_net_paid
    # 1                0                 37580            101.0
    # 1                0                 37600            3.2
    # 1                3                 34000            24.0
    # 2                5                 35000            3.5
    # 2                6                 36000            50.0
    append!(ss_customer_sk, [1,1,1,2,2])
    append!(ss_ticket_number, [0,0,3,5,6])
    append!(ss_sold_date_sk, [37580,37600,34000,35000,36000])
    append!(ss_net_paid, [101.0,3.2,24.0,3.5,50.0])

    # ws_bill_customer_sk ws_order_number ws_sold_date_sk ss_net_paid
    # 1                1                 37581            100.0
    # 1                1                 32610            30.2
    # 1                3                 34000            240.0
    # 2                5                 35000            35.0
    # 2                5                 36000            50.3
    # 3                4                 37620            30.1

    append!(ws_bill_customer_sk, [1,1,1,2,2,3])
    append!(ws_order_number, [0,0,3,5,5,4])
    append!(ws_sold_date_sk, [37581,37610,34000,35000,36000,37620])
    append!(ws_net_paid, [100.0,30.2,240.0,35.0,50.3,30.1])

    # cid   recency   frequency  totalspend
    # 1     1.0       4          468.2
    # 2     0.0       3          138.8
    # 3     1.0       1          30.1
end

function generate_dataset(ss_customer_sk, ss_ticket_number, ss_sold_date_sk, ss_net_paid,
        ws_bill_customer_sk, ws_order_number, ws_sold_date_sk, ws_net_paid,
        table_ss_path, table_ws_path)

    ss_df = readtable(open(table_ss_path))
    ws_df = readtable(open(table_ws_path))
    # order of columns in generate-dataset.sh
    #replace NA values in Dataframes with 2147483648
    append!(ss_customer_sk, convert(Array, ss_df[2], typemax(Int32)))
    append!(ss_ticket_number, convert(Array, ss_df[3], typemax(Int32)))
    append!(ss_sold_date_sk, convert(Array, ss_df[1], typemax(Int32)))
    append!(ss_net_paid, convert(Array, ss_df[4], typemax(Float32)))

    append!(ws_bill_customer_sk, ws_df[2])
    append!(ws_order_number, ws_df[3])
    append!(ws_sold_date_sk, ws_df[1])
    append!(ws_net_paid, ws_df[4])
end

function main()

    option = parse(Int,ARGS[1])
    if  option == 1 && length(ARGS) >= 3
        table_ss_path=ARGS[2]
        table_ws_path=ARGS[3]
    end

    ss_customer_sk = Int64[]
    ss_ticket_number = Int64[]
    ss_sold_date_sk = Int64[]
    ss_net_paid = Float64[]

    ws_bill_customer_sk = Int64[]
    ws_order_number = Int64[]
    ws_sold_date_sk = Int64[]
    ws_net_paid = Float64[]

    if option == 0
        example_dataset(ss_customer_sk, ss_ticket_number, ss_sold_date_sk, ss_net_paid, ws_bill_customer_sk, ws_order_number, ws_sold_date_sk, ws_net_paid)
    end

    if option == 1
        generate_dataset(ss_customer_sk, ss_ticket_number, ss_sold_date_sk, ss_net_paid, ws_bill_customer_sk, ws_order_number, ws_sold_date_sk, ws_net_paid, table_ss_path, table_ws_path)
    end
    println(":D Done reading into arrays whole dataset")
    file_name = "test_q25.hdf5"

    if isfile(file_name)
        rm(file_name)
    end
    h5write(file_name ,"/ss_customer_sk", ss_customer_sk)
    h5write(file_name ,"/ss_ticket_number", ss_ticket_number)
    h5write(file_name ,"/ss_sold_date_sk", ss_sold_date_sk)
    h5write(file_name ,"/ss_net_paid", ss_net_paid)
    h5write(file_name ,"/ws_bill_customer_sk", ws_bill_customer_sk)
    h5write(file_name ,"/ws_order_number", ws_order_number)
    h5write(file_name ,"/ws_sold_date_sk", ws_sold_date_sk)
    h5write(file_name ,"/ws_net_paid", ws_net_paid)
end
main()


text:
struct Differential <: AbstractOperator
    x::Union{Variable,Operation}
    order::Int
end
Differential(x) = Differential(x,1)

Base.show(io::IO, D::Differential) = print(io,"($(D.x),$(D.order))")
Base.Expr(D::Differential) = :($(Symbol("D_$(D.x.name)_$(D.order)")))

function Derivative end
Base.:*(D::Differential,x::Operation) = Operation(Derivative,Expression[x,D])
function Base.:*(D::Differential,x::Variable)
    if D.x === x
        return Constant(1)
    elseif D.x ∉ x.dependents
        return Constant(0)
    else
        return Variable(x,D)
    end
end
Base.:(==)(D1::Differential, D2::Differential) = D1.order == D2.order && D1.x == D2.x

Variable(x::Variable,D::Differential) = Variable(x.name,x.value,x.value_type,
                x.subtype,D,x.dependents,x.description,x.flow,x.domain,
                x.size,x.context)

function expand_derivatives(O::Operation)
    if O.op == Derivative
        D = O.args[2]
        o = O.args[1]
        simplify_constants(sum(i->Derivative(o,i)*expand_derivatives(D*o.args[i]),1:length(o.args)))
    else
        for i in 1:length(O.args)
            O.args[i] = expand_derivatives(O.args[i])
        end
    end
end
expand_derivatives(x::Variable) = x

# Don't specialize on the function here
function Derivative(O::Operation,idx)
    # This calls the Derivative dispatch from the user or pre-defined code
    Derivative(O.op,O.args,Val{idx})
end

# Pre-defined derivatives
import DiffRules, SpecialFunctions, NaNMath
for (modu, fun, arity) in DiffRules.diffrules()
    if arity ==  1 && !(fun in (:-, :+)) # :+ and :- are both unary and binary operators
        @eval begin
            function Derivative(::typeof($modu.$fun), arg, ::Type{Val{1}})
                M, f = $(modu, fun)
                @assert length(arg) == 1 "$M.$f is a unary function!"
                dx = DiffRules.diffrule(M, f, arg[1])
                parse(Operation,dx)
            end
        end
    elseif arity ==  2
        for i in 1:2
            @eval begin
                function Derivative(::typeof($modu.$fun), args, ::Type{Val{$i}})
                    M, f =  $(modu, fun)
                    if f in (:-, :+)
                        @assert length(args) in (1, 2) "$M.$f is a unary or a binary function!"
                    else
                        @assert length(args) == 2 "$M.$f is a binary function!"
                    end
                    dx = DiffRules.diffrule(M, f, args[1], args[2])[$i]
                    parse(Operation,dx)
                end
            end
        end
    end
end

function count_order(x)
    @assert !(x isa Symbol) "The variable $x must have an order of differentiation that is greater or equal to 1!"
    n = 1
    while !(x.args[1] isa Symbol)
        n = n+1
        x = x.args[1]
    end
    n, x.args[1]
end

function _differetial_macro(x)
    ex = Expr(:block)
    lhss = Symbol[]
    x = flatten_expr!(x)
    for di in x
        @assert di isa Expr && di.args[1] == :~ "@Deriv expects a form that looks like `@Deriv D''~t E'~t`"
        lhs = di.args[2]
        rhs = di.args[3]
        order, lhs = count_order(lhs)
        push!(lhss, lhs)
        expr = :($lhs = Differential($rhs, $order))
        push!(ex.args,  expr)
    end
    push!(ex.args, Expr(:tuple, lhss...))
    ex
end

macro Deriv(x...)
    esc(_differetial_macro(x))
end

function calculate_jacobian(eqs,vars)
    Expression[Differential(vars[j])*eqs[i] for i in 1:length(eqs), j in 1:length(vars)]
end

export Differential, Derivative, expand_derivatives, @Deriv, calculate_jacobian


text:
using Documenter, Graft

makedocs(
   # Options
   modules=[Graft]
)

custom_deps() = run(`pip install --user pygments mkdocs mkdocs-material`)

deploydocs(
   #options
   deps = custom_deps,
   repo = "github.com/pranavtbhat/Graft.jl.git",
   julia = "0.5"
)


text:
# # Josh: I changed update_get! -> fit_get! and update! -> fit!.  Does everything
# # still work?  I can't follow everything here.
#
# fit_get!(o::OnlineStat, args...) = (fit!(o, args...); o)
# fit_get!(o::OnlineStat, t::Tuple) = fit_get!(o, t...)
#
# function handlePipeExpr(lhs, rhs)
#
#   # println("!!PIPE!! ", lhs, " ::: ", rhs)
#   if isa(rhs, Expr) && rhs.head == :tuple
#
#     # special handling... apply this pipe to each element in the tuple, and return a tuple of the gensyms
#     blk = Expr(:block)
#     gensyms = Symbol[]
#     for rhsitem in rhs.args
#       # create a pipe expression "gensym() = lhs |> rhsitem" and add it to the block
#       gs = gensym()
#       push!(gensyms, gs)
#       pipeexpr = handlePipeExpr(lhs, rhsitem)
#       push!(blk.args, :($gs = $pipeexpr))
#     end
#
#     # now the final item in the block returns the tuple of gensyms
#     gsexpr = Expr(:tuple)
#     gsexpr.args = gensyms
#     push!(blk.args, gsexpr)
#
#     # finally done... return the block
#     return blk
#   end
#
#   # if we got here, it's a normal pipe operation
#   gs = gensym()
#   lhs = buildStreamExpr(lhs)[1]
#
#   # if lhs is a tuple, splat the arguments directly into the call to fit_get!
#   if isa(lhs, Expr) && lhs.head == :tuple
#     return :(fit_get!($rhs, $(lhs.args...)))
#   end
#
#   # must be a single argument
#   return :(fit_get!($rhs, $lhs))
#
#   # quote
#   #   $gs = $lhs
#   #   fit_get!($rhs, $gs)
#   # end
# end
#
# getMaxArg(s) = 0
# function getMaxArg(expr::Expr)
#   if expr.head == :$
#     inputnum = expr.args[1]
#     isa(inputnum, Int) || error("Cannot use dollar sign in stream macro unless referring to input (i.e. \$2 refers to the 2nd input): $expr")
#     return inputnum
#   else
#     return maximum(map(getMaxArg, expr.args))
#   end
# end
#
#
# replaceUnderscore(sym, gs::Symbol) = (sym == :_ ? gs : sym)
# function replaceUnderscore(expr::Expr, gs::Symbol)
#   # println("!!UNDER!! ", expr, " ::: ", gs)
#   expr.args = map(ex->replaceUnderscore(ex,gs), expr.args)
#   buildStreamExpr(expr)[1]
# end
#
# # create a new expression which replaces any symbol "_" in the rhs with a gensym of the lhs
# function handleCurryingExpr(lhs, rhs)
#
#   # first we process the lhs, and figure out whether there is a "mapping pipeline"
#   lhs, ismapped = buildStreamExpr(lhs)
#
#   if ismapped
#
#     # special handling if the lhs is a "mapping pipeline"
#     # for every gensym in the return tuple, inject the rhs but replacing the underscore with that element's gensym
#     # dump(lhs)
#     lhs_return_tuple = lhs.args[end]
#     for (i,lhsgs) in enumerate(lhs_return_tuple.args)
#       lhs_return_tuple.args[i] = replaceUnderscore(copy(rhs), lhsgs)
#     end
#
#     # we've now injected this expression into the return tuple of the lhs...
#     # return that, plus tell the next curry to continue mapping
#     return lhs, ismapped
#
#   else
#
#     # not mapped, so return a new expr
#     gs = gensym()
#     rhs =  replaceUnderscore(rhs, gs)
#     blk = quote
#       $gs = $lhs
#       $rhs
#     end
#     return blk, false
#   end
# end
#
# type StreamParamInfo
#   numInputs::Int
# end
# const STREAMPARAMS = StreamParamInfo(1)
#
#
# buildStreamExpr(x) = x, false
#
# function buildStreamExpr(expr::Expr)
#   head = expr.head
#   if head == :block || head == :tuple
#
#     # map build to each arg in the block
#     expr.args = map(ex -> buildStreamExpr(ex)[1], expr.args)
#     return expr, false
#
#   elseif head == :line
#
#     # keep this as-is
#     return expr, false
#
#   elseif head == :$
#
#     # if it's an integer i, replace with INPUT[i], (or INPUT for the 1-arg case)
#     inputnum = expr.args[1]
#     isa(inputnum, Int) || error("Cannot use dollar sign in stream macro unless referring to input (i.e. \$2 refers to the 2nd input): $expr")
#     @assert inputnum <= STREAMPARAMS.numInputs
#     return symbol(string("INPUT", inputnum)), false
#     # if STREAMPARAMS.numInputs > 1
#     #   return :(INPUT[$inputnum]), false
#     # else
#     #   return :INPUT, false
#     # end
#
#   elseif head == :(=)
#
#     # local variable... make sure we have a symbol on the lhs, then
#     # recursively call this on the rhs
#     # println("!!! = !!! ", expr)
#     @assert isa(expr.args[1], Symbol)
#     expr.args[2] = buildStreamExpr(expr.args[2])[1]
#     return expr, false
#
#   elseif head == :call
#
#     # operators and functions
#     fname = expr.args[1]
#     @assert isa(fname, Symbol)
#     if fname == :|>
#
#       # pipe symbol could be streaming or currying
#       lhs, rhs = expr.args[2:3]
#       if isa(rhs, Symbol)
#         return handlePipeExpr(lhs, rhs), false
#       elseif isa(rhs, Expr) && rhs.head == :tuple
#         return handlePipeExpr(lhs, rhs), true
#       elseif isa(rhs, Expr)
#         return handleCurryingExpr(lhs, rhs)
#       else
#         error("Unexpected rhs in buildStreamExpr pipe: $expr")
#       end
#
#     else
#
#       # normal expression... recursively call this on the arguments
#       # expr.args[2:end] = map(buildStreamExpr, expr.args[2:end])
#       expr.args = map(ex -> buildStreamExpr(ex)[1], expr.args)
#       return expr, false
#
#     end
#   else
#     # dump(expr,20)
#     # warn("Unexpected expr in buildStreamExpr: $expr")
#
#     # map build to each arg in the block
#     expr.args = map(ex -> buildStreamExpr(ex)[1], expr.args)
#     return expr, false
#   end
# end
#
# """
# Define pipelines of streaming data that include updating OnlineStats.
# The @stream macro creates an anonymous function to calculate a full
# data pipeline on a new data point.  Wrap in begin/end blocks for complex
# processing.  Use "\$i" to refer to the ith argument, and "_" (underscore)
# to curry results through a pipeline.
# Some features:
#   - update univariate and multivariate OnlineStats, with chaining:
#         myMean = Mean()
#         f = @stream mean(\$1 |> myMean)
#         runningMean = map(f, 0.:10)
#         # should create the series: runningMean == 0.0 : 0.5 : 5.0
#         myRegression = Adagrad(10)
#         f = @stream (\$1,\$2) |> myRegression
#         # now call f(x,y) to update the regression
#   - apply arbitrary functions and control flow:
#         myMean1 = Mean()
#         myMean2 = Mean()
#         f = @stream begin
#           if \$1 > 0.0
#             return log(mean(\$2 |> myMean1)) - 1.0
#           end
#           mean(\$2 |> myMean2)
#         end
#   - mapping pipelines allow you to update many OnlineStats with the same expression,
#     plus continue the mapping into curried results:
#       This example:
#         reg1 = Adagrad(p)
#         reg2 = SGD(p)
#         f = @stream (\$1,\$2) |> (reg1, reg2) |> \$2 - predict(_, \$1)
#       is equivalent to creating:
#         function f(x,y)
#           tmp1 = fit_get!(reg1, x, y)
#           tmp2 = fit_get!(reg2, x, y)
#           (y - predict(tmp1, x)), y - predict(tmp2, x)))
#         end
# """
# macro stream(expr::Expr)
#
#   # figure out the biggest $i arg
#   STREAMPARAMS.numInputs = getMaxArg(expr)
#   # println("GOT: ", STREAMPARAMS.numInputs)
#   if STREAMPARAMS.numInputs > 10
#     error("too many inputs in stream: maxInputs=$(STREAMPARAMS.numInputs)  expr: $expr")
#   end
#
#   # generate a unique function name, and spell out the args: ##streamed##1232(INPUT1, INPUT2)
#   fname = gensym("streamed")
#   # println(fname)
#   fargs = [symbol(string("INPUT",i)) for i in 1:STREAMPARAMS.numInputs]
#   # println(fargs)
#
#   # now create the function body
#   # print("Before: "); dump(expr, 20)
#   fbody,_ = buildStreamExpr(expr)
#   # print("After : "); dump(fbody, 20)
#   # println(fbody)
#
#
#   blk = esc(quote
#     function $fname($(fargs...))
#       $fbody
#     end
#     $fname
#   end)
#   println(blk)
#   blk
# end
#
#
# # -------------------------------------------------------------------------------
# # NOTE: everything below is an old implementation leveraging Reactive.jl... it was too slow, so I abandoned it
#
# # # some helper methods to create Input types
# # RealInput{T<:Real}(::Type{T}) = Input(zero(T))
# # FloatInput() = RealInput(Float64)
# # IntInput() = RealInput(Int)
# # VecInput{T<:Real}(::Type{T}) = Input{AVec{T}}(zeros(T,0))
# # VecInput() = VecInput(Float64)
# # RegressionInput{T<:Real}(::Type{T}) = Input{ Tuple{AVec{T},T}}((zeros(T,0),zero(T)))
# # RegressionInput() = RegressionInput(Float64)
#
#
# # function liftexpr(lhs, rhs, f = :nop)
# #   gs = gensym()
# #   fgs = f == :nop ? gs : :($f($gs))
# #   quote
# #     lift($gs -> (fit!($rhs, $fgs...); $rhs), $lhs; init = $rhs)
# #   end
# # end
#
# # applyPipe(sym::Symbol) = sym
#
# # function applyPipe(expr::Expr)
# #   @assert expr.head == :call
# #   fname = expr.args[1]
# #   @assert fname == :|>
#
# #   lhs, rhs = expr.args[2:3]
# #   @assert isa(rhs, Symbol)
# #   # TODO: we assume rhs is an OnlineStat... can we assert this and error now?
#
# #   if isa(lhs, Symbol)
# #     # both sides are symbols... return the lift expression
# #     return liftexpr(lhs, rhs)
# #   end
#
# #   # if we're here, then we have a Symbol on the RHS (presumably an OnlineStat) but something else on the LHS
# #   # NOTE: to keep it simple, assume there is a single function with a single input, which is a Signal
# #   @assert isa(lhs, Expr)
# #   @assert lhs.head == :call
# #   @assert length(lhs.args) == 2  # should be a function call with 1 param... function symbols is args[1], param is args[2]
# #   f = lhs.args[1]
# #   lhs = applyPipe(lhs.args[2])
# #   return liftexpr(lhs, rhs, f)
# # end
#
# # # pass one to many symbols/expressions that either refer to an OnlineStat object,
# # # or have it as the first argument to a function call.
# # # returns a Reactive.Input{inputType} object which you should push! the inputs to
# # macro stream(expr::Expr)
# #   expr = applyPipe(expr)
# #   # println(expr)
# #   esc(expr)
# # end


text:
"""
    load_trajs(fpath; trajrange=missing, timerange=missing)

Loads trajectories of a MCWF simulation from a .jld file at the provided `fpath`.

# Arguments
* `fpath`: savefile path (e.g. `some/valid/path/filename.jld2`). Directory and
file must exist.
* `trajrange`: MCWF trajectories to be loaded. Defaults to `1:Ntrajectories`.
* `timerange`: Timerange to be loaded. Defaults to the full timerange.

See also: [`save_trajs`](@ref)
"""
function load_trajs(fpath::String; trajrange::Union{UnitRange{T},Missing}=missing,
        timerange::Union{UnitRange{T},Missing}=missing) where T<:Integer
    @assert isfile(fpath) "ERROR: accessing "*fpath*": No such file or directory"
    times = nothing
    sols = nothing;
    jldopen(fpath, "r") do file
        fulltrajrange = 1:length(keys(file["trajs"]));
        if ismissing(trajrange)
            trajrange = fulltrajrange;
        else
            @assert trajrange[end] <= fulltrajrange[end] "The 'trajrange' argument limits must not overpass those of the saved data"
        end
        fulltimerange = 1:length(file["t"]);
        if ismissing(timerange)
            timerange = fulltimerange;
        else
            @assert timerange[end] <= fulltimerange[end] "The 'timerange' argument limits must not overpass those of the saved data"
        end
        times = file["t"][timerange]
        sols = [file["trajs/" * string(i)][t] for t in timerange, i in trajrange];
    end
    return times, sols
end;

"""
    ParallelMCWF.split_last_integer(s)

For any input string `"string_i"`, where `i` is an integer, return `("string_", i)`.
Used internally by [`ParallelMCWF.safe_fpaths`](@ref)
"""
function split_last_integer(s::String)
    num::String = ""
    for c::Char in s
        try
            n = parse(Int, c);
            num *= c;
        catch
            num="";
        end
    end
    return SubString(s,1,length(s)-length(num)), num;
end;
"""
    ParallelMCWF.safe_fpath(fpath)

For any input file path `"some/path/file.extension"`, check whether the file exists.
If so return a free file path by appending or incrementing a number at the end
of the original filename. Otherwise return `fpath`.
Used internally by [`save_trajs`](@ref)
"""
function safe_fpath(fpath::String)
    path = dirname(fpath);
    @assert isdir(path) "ERROR: accessing "*path*": No such directory"
	fpath, ext = splitext(fpath);
    while isfile(fpath*ext)
        nfpath::String, num::String = split_last_integer(fpath);
        nfpath *= num == "" ? "_1" : string(parse(Int, num)+1);
        fpath = nfpath;
    end
    return fpath*ext;
end;

"""
    save_trajs(fpath, res=missing; additional_data=missing)

Saves trajectories of a MCWF simulation to a .jld file at the provided `fpath`.

# Arguments
* `fpath`: savefile path (e.g. `some/valid/path/filename.jld2`).
Directory must pre-exist. If the file already exists, either an integer is
appended to the name or is incremented to prevent overwritting any file.
* `res=missing`: MCWF solutions on under the form `Tuple(times, sols)`.
* `additional_data=missing`: If given a `Dict`, entries are added to the
savefile.

See also: [`load_trajs`](@ref)
"""
function save_trajs(fpath::String, res::Union{Tuple{Vector{T1},Vector{T2}},Missing}=missing;
        additional_data::Union{Dict{String,T3},Missing}=missing) where {T1<:Real,T2,T3}
    if isdir(dirname(fpath))
		if isfile(fpath)
			nfpath = safe_fpath(fpath);
			@warn "$(basename(fpath)) already exists at $(dirname(fpath)).\nSaving to $nfpath..."
		end
        jldopen(nfpath, "a+") do file
            if !ismissing(res)
                times, trajs = res;
                file["t"] = times;
                for i in 1:length(trajs)
                    file["trajs/" * string(i)] = [trajs[i][t] for t in 1:length(times)];
                end
            end
            if !ismissing(additional_data)
                for (key, val) in additional_data
                    file[key] = val;
                end
            end
        end
    else
        path = dirname(fpath);
        fname = basename(fpath);
        printstyled("ERROR: accessing "*path*": No such file or directory\n",color=:light_red)
        println("Please type a valid path (like \"my/valid/path/\") or enter \"exit\"");
        npath = readline(stdin);
        if npath == "exit"
        elseif isdirpath(npath)
            save_trajs(npath*fname,res;additional_data=additional_data);
        else
            println("Entered new path not found. Create it? [y/n]");
            ans = readline(stdin);
            if ans == "y"
                mkdir(npath)
                save_trajs(npath*fname,res;additional_data=additional_data);
            end
        end
    end
end;


text:
using Test
using GridWorld

# TODO: Add Unit Tests

max_stars_repo_name: Joshua-Robison/GridWorld

max_issues_repo_name: Joshua-Robison/GridWorld

max_forks_repo_name: Joshua-Robison/GridWorld


text:
#
using Test
using Arena
using LinearAlgebra

@testset "[TEST] basics" begin

  as = ArcPointsRangeSolve([0.0;1.0;0.0],[1.0/sqrt(2.0);1.0/sqrt(2.0);0.0],[1.0;0.0;0.0],1.0)
  findaxiscenter!(as)
  @show as.center, as.axis, as.angle
  @test norm(as.center-[0.0;0.0;0.0]) < 1e-4
  @test norm(as.axis-[0.0;0.0;-1.0]) < 1e-4
  @test -1e-4 < as.angle-pi/2 < 1e-4

  as = ArcPointsRangeSolve([2.0;0;0.0],[0.5;1.5;0.0],[-1.0;0;0.0],1.5)
  findaxiscenter!(as)
  @show as.center, as.axis, as.angle
  @test norm(as.center-[0.5;0.0;0.0]) < 1e-4
  @test norm(as.axis-[0.0;0.0;1.0]) < 1e-4
  @test -1e-4 < as.angle-pi < 1e-4

  as = ArcPointsRangeSolve([0;-1.0;0],[0;0.5;1.5],[0;2.0;0.0],1.5)
  findaxiscenter!(as)
  @show as.center, as.axis, as.angle
  @test norm(as.center-[0;0.5;0]) < 1e-4
  @test norm(as.axis-[-1.0;0.0;0.0]) < 1e-4
  @test -1e-4 < as.angle-pi < 1e-4

end


text:
@testset "draw.base" begin
    @testset "atomnotation" begin
        ASSET_DIR = joinpath(dirname(@__FILE__), "..", "..", "assets", "test")
        mol = sdftomol(open(joinpath(ASSET_DIR, "demo.mol")))

        @test atomcolor(mol)[6] == Color(0, 192, 0)
        @test atomcolor(mol)[23] == Color(0, 192, 192)  # Default color
        @test isatomvisible(mol)[1] === false
        @test isatomvisible(mol)[8] === true
    end

    @testset "bondstyle" begin
        ASSET_DIR = joinpath(dirname(@__FILE__), "..", "..", "assets", "test")
        mol = sdftomol(open(joinpath(ASSET_DIR, "demo.mol")))
        coords_, styles_ = coords2d(mol)
        @test styles_[1] == 0  # 2 ニ 1
        @test styles_[21] == 1 # 17 ニ 18
        @test styles_[8] == 2  # 10 = 13
        @test styles_[25] == 0 # 27 ≡ 28
        @test styles_[15] == 1 # 3 ◀ 7
        @test styles_[5] == 6  # 4 ◁ 8
        @test styles_[28] == 4 # 5 ~ 11
    end

    @testset "atomhtml" begin
        @test atomhtml(:C, 0, 3, :right) == "CH3"
        @test atomhtml(:N, 1, 4, :left) == "+H4N"
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GTK3_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GTK3")
JLLWrappers.@generate_main_file("GTK3", UUID("77ec8976-b24b-556a-a1bf-49a033a670a6"))
end  # module GTK3_jll


text:
function vector_to_symmetric_matrix_converter()
    E1 = [
        1.0 0.0
        0.0 0.0
        0.0 1.0
    ]
    E2 = [
        0.0 0.0
        0.0 1.0
        1.0 0.0
    ]
    return [E1, E2]
end

function make_row_matrix(matrix, vals)
    return hcat([v * matrix for v in vals]...)
end

function interpolation_matrix(vals, ndofs)
    return make_row_matrix(diagm(ones(ndofs)), vals)
end

function transform_gradient(gradf, jacobian)
    return gradf / Diagonal(jacobian)
end

function bilinear_form(basis, quad, stiffness, jacobian)
    dim = dimension(basis)
    nf = number_of_basis_functions(basis)
    ndofs = dim * nf
    matrix = zeros(ndofs, ndofs)
    vectosymmconverter = vector_to_symmetric_matrix_converter()
    detjac = prod(jacobian)
    for (p, w) in quad
        grad = transform_gradient(gradient(basis, p), jacobian)
        NK = zeros(3, 2nf)
        for k = 1:dim
            NK .+= make_row_matrix(vectosymmconverter[k], grad[:, k])
        end
        matrix .+= NK' * stiffness * NK * detjac * w
    end
    return matrix
end

function bilinear_form(basis, quad, stiffness, cellmap::CellMap)
    return bilinear_form(basis, quad, stiffness, jacobian(cellmap))
end

function mass_matrix(basis, quad, detjac::R, ndofs) where {R<:Real}
    nf = number_of_basis_functions(basis)
    totaldofs = ndofs * nf
    matrix = zeros(totaldofs, totaldofs)
    for (p, w) in quad
        vals = basis(p)
        N = interpolation_matrix(vals, ndofs)
        matrix .+= N' * N * detjac * w
    end
    return matrix
end

function mass_matrix(basis, quad, cellmap::CellMap, ndofs)
    detjac = determinant_jacobian(cellmap)
    return mass_matrix(basis, quad, detjac, ndofs)
end

function mass_matrix(basis, quad, scale::V, ndofs) where {V<:AbstractVector}
    nf = number_of_basis_functions(basis)
    totaldofs = ndofs * nf
    matrix = zeros(totaldofs, totaldofs)
    @assert length(scale) == length(quad)
    for (idx, (p, w)) in enumerate(quad)
        vals = basis(p)
        N = interpolation_matrix(vals, ndofs)
        matrix .+= N' * N * scale[idx] * w
    end
    return matrix
end

function component_mass_matrix(basis, quad, component, scale)
    nf = number_of_basis_functions(basis)
    dim = dimension(basis)
    @assert dim == length(component)
    totaldofs = dim * nf
    matrix = zeros(totaldofs, totaldofs)
    projector = component * component'
    for (p, w) in quad
        vals = basis(p)
        NI = interpolation_matrix(vals, dim)
        NP = make_row_matrix(projector, vals)
        matrix .+= NI' * NP * scale * w
    end
    return matrix
end

function face_traction_operator(basis, quad, normal, stiffness, facescale, cellmap)

    dim = dimension(basis)
    @assert length(normal) == dim
    nf = number_of_basis_functions(basis)
    ndofs = dim * nf
    matrix = zeros(ndofs, ndofs)

    vectosymmconverter = vector_to_symmetric_matrix_converter()
    N = sum([normal[k] * vectosymmconverter[k]' for k = 1:dim])
    jac = jacobian(cellmap)

    for (p, w) in quad
        grad = transform_gradient(gradient(basis, p), jac)
        vals = basis(p)

        NK = zeros(3, 2nf)
        for k = 1:dim
            NK .+= make_row_matrix(vectosymmconverter[k], grad[:, k])
        end
        NI = interpolation_matrix(vals, dim)
        matrix .+= NI' * N * stiffness * NK * facescale * w
    end
    return matrix
end

function face_traction_component_operator(
    basis,
    quad,
    component,
    normal,
    stiffness,
    facescale,
    cellmap,
)
    dim = dimension(basis)
    @assert length(normal) == dim
    @assert length(component) == dim

    nf = number_of_basis_functions(basis)
    ndofs = dim * nf
    matrix = zeros(ndofs, ndofs)

    vectosymmconverter = vector_to_symmetric_matrix_converter()
    N = sum([normal[k] * vectosymmconverter[k]' for k = 1:dim])
    projector = component * component'
    jac = jacobian(cellmap)

    for (p, w) in quad
        grad = transform_gradient(gradient(basis, p), jac)
        vals = basis(p)

        NK = sum([make_row_matrix(vectosymmconverter[k], grad[:, k]) for k = 1:dim])
        NI = interpolation_matrix(vals, dim)

        matrix .+= NI' * projector * N * stiffness * NK * facescale * w
    end
    return matrix
end

function component_linear_form(rhsfunc, basis, quad, component, cellmap, detjac)
    dim = dimension(basis)
    @assert length(component) == dim
    nf = number_of_basis_functions(basis)
    rhs = zeros(dim * nf)
    for (p, w) in quad
        rhsval = rhsfunc(cellmap(p))
        NI = interpolation_matrix(basis(p), dim)
        rhs .+= rhsval * NI' * component * detjac * w
    end
    return rhs
end

function component_linear_form(rhsval, basis, quad, component, detjac)
    dim = dimension(basis)
    @assert length(component) == dim
    nf = number_of_basis_functions(basis)
    rhs = zeros(dim * nf)
    for (p, w) in quad
        NI = interpolation_matrix(basis(p), dim)
        rhs .+= rhsval * NI' * component * detjac * w
    end
    return rhs
end

function linear_form(rhsfunc, basis, quad, cellmap)
    dim = dimension(basis)
    nf = number_of_basis_functions(basis)
    rhs = zeros(dim * nf)
    detjac = determinant_jacobian(cellmap)
    for (p, w) in quad
        vals = rhsfunc(cellmap(p))
        N = interpolation_matrix(basis(p), dim)
        rhs .+= N' * vals * detjac * w
    end
    return rhs
end

function linear_form(rhsfunc, basis, quad, cellmap, detjac)
    dim = dimension(basis)
    nf = number_of_basis_functions(basis)
    rhs = zeros(dim * nf)
    for (p, w) in quad
        vals = rhsfunc(cellmap(p))
        N = interpolation_matrix(basis(p), dim)
        rhs .+= N' * vals * detjac * w
    end
    return rhs
end

function constant_linear_form(rhsval, basis, quad, detjac)
    ndofs = length(rhsval)
    nf = number_of_basis_functions(basis)
    rhs = zeros(ndofs * nf)
    for (p, w) in quad
        vals = basis(p)
        NI = interpolation_matrix(basis(p), ndofs)
        rhs .+= NI' * rhsval * detjac * w
    end
    return rhs
end


text:
# write.jl: Silo variable writing
export silowrite

function silowrite{T<:Real}(file::silofile, name::String, x_values::Vector{T}, y_values::Vector{T}; optlist::Ref{DBoptlist}=Ref{DBoptlist}())
  N = length(x_values)
  if length(y_values) != N
    error("silowrite requires x and y-values to be of the same length.")
  end
  y_values_eltype = eltype(y_values)
  if haskey(DBdatatype.JuliaSilotypemap, y_values_eltype)
    datatype = DBdatatype.JuliaSilotypemap[y_values_eltype]
  else
    error("Datatype $y_values_eltype currently unsupported.")
  end
  DBPutCurve(file.dbfile, name, x_values, y_values, datatype, N, optlist)
end
silowrite{T<:Real}(file::silofile, name::String, y_values::Vector{T}) = silowrite(file, name, convert(Vector{eltype(y_values)},1:length(y_values)), y_values)


text:
try
    Pkg.installed("SymEngine")
catch
    Pkg.add("SymEngine")
end

using Documenter, Dolang

makedocs(
    modules=[Dolang],
    authors="Spencer Lyon, Pablo Winant, and contributors",
    clean=false,
    format=:html,
    sitename="Dolang.jl",
    doctest=true,
    pages=[
        "Guide" => [
            "index.md",
            "symbolic.md",
            "compiler.md",
            "examples.md",
        ],
        "Developer Documentation" => [
            "dev/compiler.md"
        ]
    ],
)

deploydocs(
    repo   = "github.com/EconForge/Dolang.jl.git",
    julia  = "0.6",
    target = "build",
    deps   = nothing,
    make   = nothing,
)

max_issues_repo_name: EconForge/Dolang


text:
function selecteph2jld(sseph::TaylorInterpolant, bodyind::AbstractVector{Int}, tspan::Number, N::Int)
    nast = N - 11 # number of asteroids in sseph
    indvec = nbodyind(N, bodyind)
    nastout = length(bodyind) - 11 # number of asteroids to be saved in file
    @assert nastout <= nast
    sgn_yrs = sign(tspan) == 1.0 ? "p" : "m"
    nyrs_int = Int(abs(tspan))
    # write output to jld file
    ss16ast_fname = "sseph$(lpad(nast,3,'0'))ast$(lpad(nastout,3,'0'))_"*sgn_yrs*"$(nyrs_int)y_et.jld"
    ss16ast_eph = TaylorInterpolant(sseph.t0, sseph.t, sseph.x[:, union(indvec, 6N+1:6N+13)])
    jldopen(ss16ast_fname, "w") do file
        addrequire(file, TaylorSeries)
        addrequire(file, PlanetaryEphemeris)
        write(file, "ss16ast_eph", ss16ast_eph)
    end
    #check that written output is equal to original variable `ss16ast_eph`
    recovered_sol_i = load(ss16ast_fname, "ss16ast_eph")
    @show recovered_sol_i == ss16ast_eph
    return nothing
end

function propagate(maxsteps::Int, jd0::T, tspan::T;
        output::Bool=true, dense::Bool=false, ephfile::String="sseph.jld",
        dynamics::Function=NBP_pN_A_J23E_J23M_J2S!, nast::Int=343,
        quadmath::Bool=false, ss16ast::Bool=true,
        bodyind::AbstractVector{Int}=1:(11+nast), order::Int=order,
        abstol::T=abstol, parse_eqs::Bool=true) where {T<:Real}

    # total number of bodies
    N = 11+nast
    # get initial conditions (6N translational + 6 lunar mantle physical librations + 6 lunar core + TT-TDB)
    _q0 = initialcond(N, jd0) ### <--- length(_q0) == 6N+13
    # set initial time equal to zero (improves accuracy in data reductions)
    _t0 = zero(jd0)
    # final time (years)
    @show _tmax = zero(_t0)+tspan*yr

    if quadmath
        # use quadruple precision
        q0 = Float128.( _q0 )
        t0 = Float128(_t0)
        tmax = Float128(_tmax)
        _abstol = Float128(abstol)
        _jd0 = Float128(jd0)
    else
        q0 = _q0
        t0 = _t0
        tmax = _tmax
        _abstol = abstol
        _jd0 = jd0
    end

    params = (N, _jd0)

    # do integration
    if dense
        # @time sol_ = taylorinteg(dynamics, q0, t0, tmax, order, _abstol, params, maxsteps=maxsteps, dense=dense)
        @time sol_ = taylorinteg_threads(dynamics, q0, t0, tmax, order, _abstol, params, maxsteps=maxsteps, dense=dense, parse_eqs=parse_eqs)
        if quadmath
            et0 = (jd0-J2000)*daysec
            etv = Float64.( sol_.t[:]*daysec )
            sseph_x_et = map( x->x(Taylor1(order)/daysec), map(x->Taylor1(Float64.(x.coeffs)), sol_.x[:,:]) )
        else
            et0 = (jd0-J2000)*daysec
            etv = sol_.t[:]*daysec
            sseph_x_et = map(x->x(Taylor1(order)/daysec), sol_.x[:,:])
        end
        sseph = TaylorInterpolant(et0, etv, sseph_x_et)
        sol = (sseph=sseph,)
    else
        # @time sol_ = taylorinteg(dynamics, q0, t0, tmax, order, _abstol, params, maxsteps=maxsteps, dense=dense)
        @time sol_ = taylorinteg_threads(dynamics, q0, t0, tmax, order, _abstol, params, maxsteps=maxsteps, dense=dense, parse_eqs=parse_eqs)
        sol = (t=sol_[1][:], x=sol_[2][:,:])
    end

    #write solution to .jld files
    if output
        if dense && ss16ast
            selecteph2jld(sseph, bodyind, tspan, N)
        else
            println("Saving solution to file: $ephfile")
            jldopen(ephfile, "w") do file
                addrequire(file, TaylorSeries)
                addrequire(file, PlanetaryEphemeris)
                # write variables to jld file
                for ind in eachindex(sol)
                    varname = string(ind)
                    println("Saving variable: ", varname)
                    write(file, varname, sol[ind])
                end
            end
            #check that recovered variables are equal to original variables
            for ind in eachindex(sol)
                varname = string(ind)
                #read varname from jld file and assign recovered variable to recovered_sol_i
                recovered_sol_i = load(ephfile, varname)
                #check that recovered variable is equal to original variable
                @show recovered_sol_i == sol[ind]
            end
        end
        println("Saved solution")
        return nothing
    else
        return sol
    end
end


text:
using Reexport
@reexport using Entropies
using Entropies: NearestNeighborEntropyEstimator

import DelayEmbeddings: AbstractDataset, Dataset
export Dataset

using Distances, SpecialFunctions

const Vector_or_Dataset{D, T} = Union{AbstractVector{T}, AbstractDataset{D, T}} where {D, T}
const Est = Union{ProbabilitiesEstimator, EntropyEstimator}


text:
include("./representation.jl")
include("./parent_selection.jl")
include("./recombination.jl")
include("./mutation.jl")
include("./metrics.jl")
include("./functions_collection.jl")
include("./utils.jl")
include("./survivor_selection.jl")
include("./graphing_utilities.jl")

using .Representation, .ParentSelection, .Recombination, .Mutation, .Metrics, .FunctionsCollection
using .SurvivorSelection, .Utils, Statistics, .GraphingUtilities

parameters =  Utils.get_algo_params()
x_values = Utils.generate_linespace(parameters["x_range"][1], parameters["x_range"][2], parameters["population_size"])
true_values = parameters["evaluate_function"](x_values)
population = Representation.generate_population(parameters["population_size"], parameters["n_kernels"], parameters["σ_initial"])
extra_params_error = Dict("threshold"=>parameters["threshold"],
                        "lower_weight"=>parameters["lower_weight"],
                        "upper_weight"=>parameters["upper_weight"])

Metrics.calc_population_fitness!(population,
                                 true_values,
                                 parameters["error_function"],
                                 extra_params_error)

all_fitness_values = Array{Array{Float64, 1}, 1}[]
all_hit_ratios = Array{Array{Float64, 1}, 1}[]
for run in 1:parameters["num_runs"]
    fitness_values = Array{Float64, 1}[]
    hit_ratios = Array{Float64, 1}[]
    for i in 1:parameters["num_generations"]
        println("Generation: ", i)
        println("Selecting parents...")
        parents_groups = ParentSelection.random_parent_selection(population, parameters["λ_children"],
                                                                parameters["μ_parents"])
        println("Applying recombination...")
        offspring = Recombination.apply_recombination(parents_groups, parameters["μ_parents"],
                                                    parameters["recombination_type"])
        println("Applying mutation")
        offspring_population = Representation.Population(offspring)
        mutated_population = parameters["mutation_type"](offspring_population)
        println("Calculating fitness...")
        Metrics.calc_population_fitness!(mutated_population,
                                         true_values,
                                         parameters["error_function"],
                                         extra_params_error)
        println("Selecting survivors...")
        next_generation = nothing
        if parameters["survivor_selection"] == "children"
            next_generation = SurvivorSelection.rank_selection(mutated_population, parameters["μ_parents"])
        elseif parameters["survivor_selection"] == "all"
            next_generation = SurvivorSelection.rank_selection(mutated_population, population, parameters["μ_parents"])
        else
            throw(ArgumentError("Survivor selection requested not available. Options are: all, children"))
        end
        push!(fitness_values, [member.fitness for member in next_generation.members])
        push!(hit_ratios, [member.hit_ratio for member in next_generation.members])
        global population = next_generation
    end
    push!(all_fitness_values, fitness_values)
    push!(all_hit_ratios, hit_ratios)
end
println("Finished training...")
mean_runs_fitness = mean(all_fitness_values)
mean_hit_ratio = mean(all_hit_ratios)

output_name = nothing
if parameters["output_file"] != nothing
    output_name = "outputs/" * parameters["output_file"]
end

GraphingUtilities.plot_all(mean_runs_fitness,
                            parameters["num_generations"],
                            population,
                            true_values,
                            x_values,
                            mean_hit_ratio,
                            parameters["population_size"],
                            output_name)

open(output_name * ".txt", "w") do file
    write(file, string([mean(values) for values in mean_runs_fitness]))
    write(file, "\n")
    write(file, string([mean(values) for values in mean_hit_ratio]))
    write(file, "\n")
end


text:
module Blackscholes

using GPUBenchmarks, BenchmarkTools
import CUDAdrv

using CUDAnative
const cu = CUDAnative

const description = """
Blackschole is a nice benchmark for broadcasting performance.
It's a medium heavy calculation per array element, where the calculation is completely
independant from each other.
The CuArray package is a bit slower here compared to GPUArrays, which should be straightforward to fix.
I suspect that it's due to more promotions between integer types in the indexing code.
"""

function blackscholes(
        sptprice,
        strike,
        rate,
        volatility,
        time
    )
    logterm = log( sptprice / strike)
    powterm = .5f0 * volatility * volatility
    den = volatility * sqrt(time)
    d1 = (((rate + powterm) * time) + logterm) / den
    d2 = d1 - den
    NofXd1 = cndf2(d1)
    NofXd2 = cndf2(d2)
    futureValue = strike * exp(-rate * time)
    c1 = futureValue * NofXd2
    call_ = sptprice * NofXd1 - c1
    put  = call_ - futureValue + sptprice
    return put
end
function blackscholes_af(
        sptprice,
        strike,
        rate,
        volatility,
        time
    )
    logterm = log( sptprice / strike)
    powterm = .5f0 * volatility * volatility
    den = volatility * sqrt(time)
    d1 = (((rate + powterm) * time) + logterm) / den
    d2 = d1 - den
    NofXd1 = cndf2(d1)
    NofXd2 = cndf2(d2)
    futureValue = strike * exp(-rate * time)
    c1 = futureValue * NofXd2
    call_ = sptprice * NofXd1 - c1
    put  = call_ - futureValue + sptprice
    return put
end

@inline function cndf2(x)
    0.5f0 + 0.5f0 * erfc(0.707106781f0 * x)
end

function cu_blackscholes(sptprice, strike, rate, volatility, time)
    logterm = cu.log( sptprice / strike)
    powterm = .5f0 * volatility * volatility
    den = volatility * cu.sqrt(time)
    d1 = (((rate + powterm) * time) + logterm) / den
    d2 = d1 - den
    NofXd1 = cu_cndf2(d1)
    NofXd2 = cu_cndf2(d2)
    futureValue = strike * cu.exp(- rate * time)
    c1 = futureValue * NofXd2
    call_ = sptprice * NofXd1 - c1
    put  = call_ - futureValue + sptprice
    return put
end

function cu_cndf2(x)
    0.5f0 + 0.5f0 * cu.erfc(0.707106781f0 * x)
end

function execute_broadcast(f, res, a, b, c, d, e)
    res .= f.(a, b, c, d, e)
    synchronize(res)
end

function execute_broadcast_af(f, res, a, b, c, d, e)
    res .= f.(a, b, c, d, e)
    synchronize(res)
end

function makeresult(bench, N, device, hardware, mdiff)
    BenchResult(
        "Blackscholes",
        bench,
        N,
        Complex64,
        string(device),
        hardware,
        @__FILE__,
        mdiff
    )
end

function execute(device)
    hardware, array_type = init(device)
    T = Float32
    results = BenchResult[]
    func = is_cudanative(device) ? cu_blackscholes : blackscholes
    func = is_arrayfire(device) ? blackscholes_af : func
    for x in 1:7
        N = 10^x
        jl_arrays = (
            zeros(T, N),
            T[42.0 for i = 1:N],
            T[40.0 + (i / N) for i = 1:N],
            T[0.5 for i = 1:N],
            T[0.2 for i = 1:N],
            T[0.5 for i = 1:N]
        )
        gpu_arrays = array_type.(jl_arrays)
        res, a, b, c, d, e = gpu_arrays
        br_func = is_arrayfire(device) ? execute_broadcast_af : execute_broadcast
        bench = @benchmark $(br_func)($func, $res, $a, $b, $c, $d, $e)
        jl_arrays[1] .= blackscholes.(jl_arrays[2:end]...)
        mdiff = meandifference(jl_arrays[1], res)
        push!(results, makeresult(bench, N, device, hardware, mdiff))
        for elem in gpu_arrays
            free(elem)
            println("freeing element: ", CUDAdrv.Mem.used() / 10^7)
        end
        gc()
    end
    return results
end

end


text:
@testset "StandardModel utilities" begin
    model = load_model(StandardModel, model_paths["e_coli_core.json"])

    # FBA
    fluxes = flux_balance_analysis_dict(
        model,
        Tulip.Optimizer;
        modifications = [change_objective("BIOMASS_Ecoli_core_w_GAM")],
    )

    # bounds setting
    cbm = make_optimization_model(model, Tulip.Optimizer)
    ubs = cbm[:ubs]
    lbs = cbm[:lbs]
    glucose_index = first(indexin(["EX_glc__D_e"], reactions(model)))
    o2_index = first(indexin(["EX_o2_e"], reactions(model)))
    atpm_index = first(indexin(["ATPM"], reactions(model)))
    set_optmodel_bound!(glucose_index, cbm; ub = -1.0, lb = -1.0)
    @test normalized_rhs(ubs[glucose_index]) == -1.0
    @test normalized_rhs(lbs[glucose_index]) == 1.0
    set_optmodel_bound!(o2_index, cbm; ub = 1.0, lb = 1.0)
    @test normalized_rhs(ubs[o2_index]) == 1.0
    @test normalized_rhs(lbs[o2_index]) == -1.0
end


text:
using Test
using MPI
using CLIMA
using CLIMA.ConfigTypes
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.VTK
using Logging
using Printf
using LinearAlgebra
using CLIMA.DGmethods: DGModel, init_ode_state, LocalGeometry, courant
using CLIMA.DGmethods.NumericalFluxes:
    Rusanov, CentralNumericalFluxGradient, CentralNumericalFluxDiffusive
using CLIMA.Courant
using CLIMA.Atmos:
    AtmosModel,
    AtmosAcousticLinearModel,
    RemainderModel,
    FlatOrientation,
    NoReferenceState,
    ReferenceState,
    DryModel,
    NoRadiation,
    NoPrecipitation,
    Gravity,
    HydrostaticState,
    IsothermalProfile,
    ConstantViscosityWithDivergence,
    vars_state,
    soundspeed
using CLIMA.Atmos
using CLIMA.ODESolvers

using CLIMAParameters
using CLIMAParameters.Planet: kappa_d
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

using CLIMA.MoistThermodynamics:
    air_density, total_energy, internal_energy, soundspeed_air

using CLIMA.VariableTemplates: Vars
using StaticArrays

const p∞ = 10^5
const T∞ = 300.0

function initialcondition!(bl, state, aux, coords, t)
    FT = eltype(state)

    translation_speed::FT = 150
    translation_angle::FT = pi / 4
    α = translation_angle
    u∞ = SVector(
        FT(translation_speed * coords[1]),
        FT(translation_speed * coords[1]),
        FT(0),
    )
    _kappa_d::FT = kappa_d(param_set)

    u = u∞
    T = FT(T∞)
    # adiabatic/isentropic relation
    p = FT(p∞) * (T / FT(T∞))^(FT(1) / _kappa_d)
    ρ = air_density(T, p, bl.param_set)

    state.ρ = ρ
    state.ρu = ρ * u
    e_kin = u' * u / 2
    state.ρe = ρ * total_energy(e_kin, FT(0), T, bl.param_set)

    nothing
end

let
    # boiler plate MPI stuff
    CLIMA.init()
    ArrayType = CLIMA.array_type()
    mpicomm = MPI.COMM_WORLD

    # Mesh generation parameters
    N = 4
    Nq = N + 1
    Neh = 10
    Nev = 4

    @testset "$(@__FILE__) DGModel matrix" begin
        for FT in (Float64, Float32)
            for dim in (2, 3)
                if dim == 2
                    brickrange = (
                        range(FT(0); length = Neh + 1, stop = 1),
                        range(FT(1); length = Nev + 1, stop = 2),
                    )
                elseif dim == 3
                    brickrange = (
                        range(FT(0); length = Neh + 1, stop = 1),
                        range(FT(0); length = Neh + 1, stop = 1),
                        range(FT(1); length = Nev + 1, stop = 2),
                    )
                end
                μ = FT(2)
                topl = StackedBrickTopology(mpicomm, brickrange)

                grid = DiscontinuousSpectralElementGrid(
                    topl,
                    FloatType = FT,
                    DeviceArray = ArrayType,
                    polynomialorder = N,
                )

                model = AtmosModel{FT}(
                    AtmosLESConfigType;
                    ref_state = NoReferenceState(),
                    turbulence = ConstantViscosityWithDivergence(μ),
                    moisture = DryModel(),
                    source = Gravity(),
                    boundarycondition = (),
                    init_state = initialcondition!,
                    param_set = param_set,
                )

                dg = DGModel(
                    model,
                    grid,
                    Rusanov(),
                    CentralNumericalFluxDiffusive(),
                    CentralNumericalFluxGradient(),
                )

                Δt = FT(1 // 200)

                Q = init_ode_state(dg, FT(0))

                Δx = min_node_distance(grid, EveryDirection())
                Δx_v = min_node_distance(grid, VerticalDirection())
                Δx_h = min_node_distance(grid, HorizontalDirection())

                translation_speed = FT(norm([150.0, 150.0, 0.0]))
                diff_speed_h =
                    FT(μ / air_density(FT(T∞), FT(p∞), model.param_set))
                diff_speed_v =
                    FT(μ / air_density(FT(T∞), FT(p∞), model.param_set))
                c_h =
                    Δt * (
                        translation_speed +
                        soundspeed_air(FT(T∞), model.param_set)
                    ) / Δx_h
                c_v = Δt * (soundspeed_air(FT(T∞), model.param_set)) / Δx_v
                d_h = Δt * diff_speed_h / Δx_h^2
                d_v = Δt * diff_speed_v / Δx_v^2
                simtime = FT(0)

                # tests for non diffusive courant number
                @test courant(
                    nondiffusive_courant,
                    dg,
                    model,
                    Q,
                    Δt,
                    simtime,
                    HorizontalDirection(),
                ) ≈ c_h rtol = 1e-4
                @test courant(
                    nondiffusive_courant,
                    dg,
                    model,
                    Q,
                    Δt,
                    simtime,
                    VerticalDirection(),
                ) ≈ c_v rtol = 1e-4

                # tests for diffusive courant number
                @test courant(
                    diffusive_courant,
                    dg,
                    model,
                    Q,
                    Δt,
                    simtime,
                    HorizontalDirection(),
                ) ≈ d_h
                @test courant(
                    diffusive_courant,
                    dg,
                    model,
                    Q,
                    Δt,
                    simtime,
                    VerticalDirection(),
                ) ≈ d_v
            end
        end
    end
end

nothing

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
"""
    solve_btpde_midpoint(model, matrices, experiment)

Solve the Bloch-Torrey partial differential equation using P1 finite elements.
This function uses a manual time stepping scheme (theta-rule), that requires a degree of
implicitness `θ` and a time step `Δt`.
    `θ = 0.5`: Crank-Nicolson (second order)
    `θ = 1.0`: Implicit Euler (first order)
"""
function solve_btpde_midpoint(model::Model, matrices, experiment::Experiment)

    # Measure function evalutation time
    starttime = Base.time()

    # Extract input parameters
    @unpack mesh, D, T₂, ρ = model
    @unpack M, S, R, Mx, Q, M_cmpts = matrices
    @unpack directions, sequences, values, values_type = experiment.gradient
    @unpack θ, timestep = experiment.btpde_midpoint

    # Deduce sizes
    ncompartment = length(ρ)
    npoint_cmpts = size.(mesh.points, 2)
    inds_cmpts = cumsum([0; npoint_cmpts])
    ndirection = size(directions, 2)
    nsequence = length(sequences)
    namplitude = length(values)

    qvalues, bvalues = get_values(experiment.gradient)

    # Create initial conditions (enforce complex values)
    ρ = vcat(fill.(complex.(ρ), npoint_cmpts)...)

    # Allocate output arrays
    signal = zeros(ComplexF64, ncompartment, namplitude, nsequence, ndirection)
    signal_allcmpts = zeros(ComplexF64, namplitude, nsequence, ndirection)
    magnetization =
        Array{Matrix{ComplexF64},4}(undef, ncompartment, namplitude, nsequence, ndirection)
    time = Array{Vector{Float64},3}(undef, namplitude, nsequence, ndirection)
    itertimes = zeros(namplitude, nsequence, ndirection)

    Jac(q, fₜ, g) = -(S + Q + R + im * fₜ * q * (g' * Mx))

    # Iterate over gradient amplitudes, time profiles and directions
    for idir = 1:ndirection, iseq = 1:nsequence, iamp = 1:namplitude

        # Measure solving time
        itertime = Base.time()

        # Gradient direction
        dir = directions[:, idir]

        # Time profile
        f = sequences[iseq]
        ivals = intervals(f)

        # Gradient amplitude
        q = qvalues[iamp, iseq]
        b = bvalues[iamp, iseq]

        # Display state of iterations
        @printf "Solving BTPDE with size %d\n" (sum(npoint_cmpts))
        @printf "  Direction %d of %d: g = [%.2f, %.2f, %.2f]\n" idir ndirection dir...
        @printf "  Sequence  %d of %d: f = %s\n" iseq nsequence f
        @printf "  Amplitude %d of %d: q = %g, b = %g\n" iamp namplitude q b

        # Crank-Nicolson time stepping
        t = 0.0
        y = copy(ρ)
        Ey = copy(y)
        for i = 1:length(ivals)-1
            @printf "    Solving for interval [%g, %g]\n" ivals[i] ivals[i+1]

            # Adjust time step to obtain divide interval uniformly
            ival_length = ivals[i+1] - ivals[i]
            nt = round(Int, ival_length / timestep)
            Δt = ival_length / nt

            # Build matrices for interval
            J = Jac(q, f(ivals[i]), dir)
            F = lu(M .- Δt .* θ .* J)
            # F = factorize(M .- Δt .* θ .* J)
            E = @. complex(M) + Δt * (1 - θ) * J

            # Advance one interval
            # for it = 1:nt
            while t + Δt ≤ ivals[i+1]
                # @show t
                mul!(Ey, E, y)
                ldiv!(y, F, Ey)
                t += Δt
            end
        end

        # Extract solution
        time[iamp, iseq, idir] = [t]
        mag = [y;]

        # Split solution into compartments
        for icmpt = 1:ncompartment
            inds = inds_cmpts[icmpt]+1:inds_cmpts[icmpt+1]

            # Store magnetization in compartment
            magnetization[icmpt, iamp, iseq, idir] = mag[inds, :]

            # Integrate final magnetization over compartment
            signal[icmpt, iamp, iseq, idir] =
                sum(M_cmpts[icmpt] * mag[inds, end], dims = 1)[1]
        end

        itertimes[iamp, iseq, idir] = Base.time() - itertime
    end

    signal_allcmpts = sum(signal, dims = 1)[1, :, :, :]

    totaltime = Base.time() - starttime

    (; magnetization, signal, signal_allcmpts, time, itertimes, totaltime)
end


text:
module Writers

using Avro.DataFile
using Avro.DataFile.Codecs
using Avro.Generic
using Avro.Io
using Avro.Schemas

export wrap

mutable struct Writer{CodecName}
    schema::Schemas.Schema
    output_encoder::BinaryEncoder
    buffer_encoder::BinaryEncoder

    # Users may modify these settings
    codec::Codec{CodecName}
    sync_marker::Vector{UInt8}
    sync_interval::Int

    # Internal, mutable state
    block_count::Int
end

function Writer(
        schema::Schemas.Schema,
        output_encoder::BinaryEncoder,
        buffer_encoder::BinaryEncoder,
        codec::Codec{CodecName},
        sync_marker::Vector{UInt8},
        sync_interval::Int) where CodecName
    Writer(
        schema,
        output_encoder,
        buffer_encoder,
        codec,
        sync_marker,
        sync_interval,
        0)
end

function wrap(
        schema::Schemas.Schema,
        output::IO;
        codec_name::String = "null",
        sync_marker::Vector{UInt8} = _generate_sync_marker(),
        sync_interval::Integer = 256)

    # Initialize the output and buffer encoders
    output_encoder = BinaryEncoder(output)
    buffer_encoder = BinaryEncoder(IOBuffer())
    codec = Codecs.create(codec_name)

    # Initialize the file writer
    file_writer = Writer(
        schema,
        output_encoder,
        buffer_encoder,
        codec,
        sync_marker,
        sync_interval)

    # Write the header
    _write_header(file_writer)

    # Return the file writer
    file_writer
end

function Base.write(file_writer::Writer, datum)
    buffer_encoder = file_writer.buffer_encoder
    schema = file_writer.schema

    # Write the datum to the buffer encoder using given schema
    write(buffer_encoder, schema, datum)

    # Increment the number of records in the block
    file_writer.block_count += 1

    # Flush the buffer to the output if we have reached the sync interval
    if file_writer.block_count >= file_writer.sync_interval
        write_block(file_writer)
    end
end

function Base.close(file_writer::Writer)
    write_block(file_writer)

    close(file_writer.output_encoder.stream)
end

function write_block(file_writer::Writer)
    if file_writer.block_count > 0
        buffer_encoder = file_writer.buffer_encoder
        output_encoder = file_writer.output_encoder
        codec = file_writer.codec

        # Extract the buffer data as bytes, compress the buffer data, and
        # reset the buffer state
        buffer_data = take!(buffer_encoder.stream)
        buffer_data = Codecs.compress(codec, buffer_data)

        # Write number of records, the blocks size (bytes), the data, and then
        # the sync marker
        encode_long(output_encoder, file_writer.block_count)
        encode_long(output_encoder, length(buffer_data))
        encode_fixed(output_encoder, buffer_data)
        encode_fixed(output_encoder, file_writer.sync_marker)

        # Reset the block counter
        file_writer.block_count = 0
    end
end

"""
Generate a 16 byte synchronization marker.
"""
_generate_sync_marker() = rand(UInt8, 16)

"""
Generates a header for an object container file.
"""
function _generate_header(schema, sync_marker, codec::Codec)
    GenericRecord(
        METADATA_SCHEMA,
        [
            GenericFixed(MAGIC_SCHEMA, OBJECT_CONTAINER_FILE_MAGIC),
            Dict(
                META_CODEC_KEY => Codecs.name(codec),
                META_SCHEMA_KEY => string(schema),
            ),
            GenericFixed(SYNC_SCHEMA, sync_marker)
        ]
    )
end

function _write_header(file_writer::Writer)
    header = _generate_header(
        file_writer.schema,
        file_writer.sync_marker,
        file_writer.codec)

    # Write the header directly the output encoder
    write(file_writer.output_encoder, METADATA_SCHEMA, header)
end

end


text:
using LazyUtils
using Test
using Pkg

@testset "LazyUtils" begin
    @testset "PkgUtils" begin
        environment = joinpath(Pkg.envdir(), "v$(VERSION.major).$(VERSION.minor)")
        projectfile = joinpath(environment, "Project.toml")
        environments = get_environments()
        projectfiles = get_projectfiles()
        @test any(occursin.(environment, environments)) || environments === nothing
        @test any(occursin.(projectfile, projectfiles)) || projectfiles === nothing
        @test get_pkg_version("Test") === nothing
        @test get_pkg_version("TOML").major == 1
    end
end


text:
using Test
import IPFitting
import IPFitting.FIO

##
println("===================================")
println("Testing IPFitting.FIO")
println("===================================")

for filetype in [".h5", ".json"]

println("-----------------------------")
println("Testing FIO with $filetype.")
println("-----------------------------")
fname1 = tempname() * filetype
fname2 = tempname() * filetype

D = Dict( "a" => rand(),
          "b" => rand(3,4,10),
          "c" => Dict( "d" => rand(10), "e" => "mystring"),
          "f" => nothing )

if filetype == ".json"
   D["b"] = D["b"][:]
end

println("Save and Load Entire Dict")
FIO.save(fname1, D)
D1 = FIO.load(fname1)
println(@test D1 == D)

println("Load an element of D")
a = FIO.load(fname1, "a")
println(@test a == D["a"])

println("Save elements by name")
FIO.save(fname2, collect(D)...)
D2 = FIO.load(fname2)
println(@test D2 == D)

if filetype == ".h5"
   println("Read from multi-dim array")
   bsub = FIO.load(fname1, ("b", (:,:,2:5)) )
   println(@test bsub == D["b"][:,:,2:5])
end

println("Delete temporary files")
run(`rm $fname1 $fname2`)

end

println("===================================")


text:
using LoopVectorization: @tvectorize

function eval_exp(N)
    a = range(0, stop=2*pi, length=N)
    A = Matrix{ComplexF64}(undef, N, N)

    @inbounds Threads.@threads for j in 1:N
    for i in 1:N
        A[i,j] = exp((100+im)*im*sqrt(a[i]^2 + a[j]^2))
    end
    end

    return A

end

function eval_exp_tvectorize(N)
    a = range(0, stop=2*pi, length=N)
    A = Matrix{ComplexF64}(undef, N, N)

    _A = reinterpret(reshape, Float64, A)
    @tvectorize for j in 1:N, i in 1:N
        x = sqrt(a[i]^2 + a[j]^2)
        prefac = exp(-x)
        s, c = sincos(100*x)

        _A[1,i,j] = prefac * c
        _A[2,i,j] = prefac * s
    end
    return A
end

eval_exp(5)
print(string("running loop on ", Threads.nthreads(), " threads \n"))
for N in 1000:1000:10000
    @time begin
    A = eval_exp(N)
    end

    @time begin
    A = eval_exp_tvectorize(N)
    end
end

max_stars_repo_name: mdmaas/julia-numpy-fortran-test

max_issues_repo_name: mdmaas/julia-numpy-fortran-test

max_forks_repo_name: mdmaas/julia-numpy-fortran-test


text:
conn = DBInterface.connect(MySQL.Connection, "127.0.0.1", "root", "root")

DBInterface.execute(conn, "CREATE DATABASE db_cabida")
DBInterface.execute(conn, "use db_cabida")

DBInterface.execute(conn, """CREATE TABLE datosCabidaNormativa
                 (
                     ID INT NOT NULL AUTO_INCREMENT,
                     SEPMIN DOUBLE,
                     ANTEJARDIN DOUBLE,
                     RASANTE DOUBLE,
                     ALTURAMAX DOUBLE,
                     MAXPISOS INT,
                     COEFOCUPACION DOUBLE,
                     SUBPREDIALMIN DOUBLE,
                     DENSIDADMAX DOUBLE,
                     FLAGDENSIDADBRUTA DOUBLE,
                     COEFCONSTRUCTIBILIDAD DOUBLE,
                     ESTACIONAMIENTOSPORVIV DOUBLE,
                     PORCADICESTACVISITAS DOUBLE,
                     SUPPORESTACIONAMIENTO DOUBLE,
                     ESTBICICLETAPOREST DOUBLE,
                     BICICLETASPOREST DOUBLE,
                     FLAGCAMBIOESTPORBICICLETA BIT(1),
                     MAXSUBTE INT,
                     COEFOCUPACIONEST DOUBLE,
                     SEPESTMIN DOUBLE,
                     REDUCCIONESTPORDISTMETRO DOUBLE,
                     PRIMARY KEY (ID)
                 );""")

DBInterface.execute(conn, """INSERT INTO datosCabidaNormativa (SEPMIN, ANTEJARDIN, RASANTE, ALTURAMAX, MAXPISOS, COEFOCUPACION, SUBPREDIALMIN, DENSIDADMAX, FLAGDENSIDADBRUTA, COEFCONSTRUCTIBILIDAD, ESTACIONAMIENTOSPORVIV, PORCADICESTACVISITAS, SUPPORESTACIONAMIENTO, ESTBICICLETAPOREST, BICICLETASPOREST, FLAGCAMBIOESTPORBICICLETA, MAXSUBTE, COEFOCUPACIONEST, SEPESTMIN, REDUCCIONESTPORDISTMETRO)
                 VALUES
                 (5, 5, 2.75, 60, 30, .7, 1000, 2500, 1, 4, 1, .15, 34, .25, 3, 1, 10, .8, 7, 0);
              """)

queryData = DBInterface.execute(conn, "select * from datosCabidaNormativa")

row = DBInterface.execute(conn, "select * from datosCabidaNormativa") |> columntable
A = [getproperty(row, prop) for (i, prop) in enumerate(propertynames(row))]


text:
module image_transform
export image_to_vector

using Images

function image_to_vector(filepath::String)
    image = load(filepath)
    inverted_bw = broadcast(x -> 1.0 - x, Gray.(image)) # black on white to white on black
    not_black = broadcast(!=, inverted_bw, Gray{N0f8}(0.0)) # pixels which are not black
    bounding_box = (
        up = findall(any(not_black, dims=2)[:, 1])[begin],
        down = findall(any(not_black, dims=2)[:, 1])[end],
        left = findall(any(not_black, dims=1)[1, :])[begin],
        right = findall(any(not_black, dims=1)[1, :])[end]
    )
    canvas_dims = (28, 28)
    canvas = zeros((canvas_dims))

    bounding_box_dims = (bounding_box.down - bounding_box.up, bounding_box.right - bounding_box.left)
    rescaled_height, rescaled_width = round.(Int, (bounding_box_dims) .* (20 / maximum(bounding_box_dims)))

    # image within bounding box is scaled (maintaining aspect ratio) to fit in 20x20 and pasted in the center of the canvas
    canvas[
        canvas_dims[1]÷2 - rescaled_height÷2 : canvas_dims[1]÷2 + (rescaled_height - rescaled_height÷2) - 1,
        canvas_dims[2]÷2 - rescaled_width÷2 : canvas_dims[2]÷2 + (rescaled_width - rescaled_width÷2) - 1,
        ] = imresize(
        inverted_bw[bounding_box.up:bounding_box.down, bounding_box.left:bounding_box.right]
            |> dilate
            |> dilate
            |> dilate,
        (rescaled_height, rescaled_width)
    )

    return reduce(vcat, Float64.(canvas)')
end
end

max_stars_repo_name: ShakyaMajumdar/neural-network

max_issues_repo_name: ShakyaMajumdar/neural-network

max_forks_repo_name: ShakyaMajumdar/neural-network