text:
module JGUI

import Base: show
import Base: getindex, setindex!, length,
             push!, append!, prepend!, insert!, splice!, delete!, shift!, unshift!, pop!,
             findin
import Base: size, endof, ndims
import Base: connect, notify

export properties

export getValue, setValue, setIcon
export replace!

export disconnect

export window,
       destroy, raise, lower

export labelframe,
       box, hbox, vbox, addstretch, addstrut, addspacing,
       formlayout,
       notebook,
       children,
       grid,
       row_minimum_height, column_minimum_width, row_stretch, column_stretch

export label, separator, button, lineedit, textedit,
       checkbox, radiogroup, buttongroup, combobox,
       slider, slider2d, spinbox,
       listview, storeview, treeview,
       imageview,
       icon

export Store, TreeStore

export treestore, expand_node, collapse_node, node_to_path, path_to_node, update_node

export filedialog, messagebox, confirmbox, dialog

export action, menubar, menu,
       addMenu, addAction

## adjust this to pick a gtoolkit
function pick_toolkit()
    toolkits = ["Gtk", "Tk", "Qt"]
    toolkit_map = {"Gtk"=>"Gtk", "Tk"=>"Tk", "Qt"=>"PySide"}

    ## filter these

    println("Which toolkit to use:")
    for (i, kit) in enumerate(toolkits)
        println("\t[$i] $kit")
    end
    i = parseint(chomp(readline()))

    toolkit = toolkits[i]
    ENV["toolkit"] = toolkit
end

## which toolkit?
#__init__() = haskey(ENV, "toolkit") || pick_toolkit()
if !haskey(ENV, "toolkit")
    ENV["toolkit"] = "Gtk"      # default to Gtk for now
end

export manipulate

isqt() = lowercase(ENV["toolkit"]) == lowercase("Qt")
istk() = lowercase(ENV["toolkit"]) == lowercase("Tk")
isgtk() = lowercase(ENV["toolkit"]) == lowercase("Gtk")

if istk()
    using Tk
elseif isqt()
    using PyCall
    using PySide
    import PySide: raise, setFocus, setIcon
elseif isgtk()
    using Gtk, Cairo
end

using Docile
@docstrings

include("types.jl")

include("methods.jl")
include("icons.jl")
include("models.jl")
include("containers.jl")
include("widgets.jl")
include("dialogs.jl")
include("menu.jl")
include("manipulate.jl")        # code depends on Tk or Qt

if istk()
    default_toolkit = MIME("application/x-tcltk")
    include("tk.jl")
    ENV["WINSTON_OUTPUT"] = :tk
    export cairographic
    include("winston.jl")
elseif isqt()
    default_toolkit = MIME("application/x-qt")
    include("qt.jl")
    export pyplotgraphic
elseif isgtk()
    default_toolkit = MIME("application/x-gtk")
    include("gtk.jl")
    ENV["WINSTON_OUTPUT"] = :gtk
    export cairographic
    include("winston.jl")
end

end


text:
##############################################################################
#
# InvertedGamma distribution, cf Bayesian Theory (Barnardo & Smith) p 119
#
#  Note that B&S parametrize in terms of shape/rate, but this is uses
#  shape/scale so that it is consistent with the implementation of Gamma
#
##############################################################################

immutable InvertedGamma <: ContinuousUnivariateDistribution
    shape::Float64 # location
    scale::Float64 # scale
    function InvertedGamma(sh::Real, sc::Real)
    	sh > zero(sh) && sc > zero(sc) || error("Both shape and scale must be positive")
    	new(float64(sh), float64(sc))
    end
end

scale(d::InvertedGamma) = d.scale
rate(d::InvertedGamma) = 1.0 / d.scale

insupport(::InvertedGamma, x::Real) = zero(x) <= x < Inf
insupport(::Type{InvertedGamma}, x::Real) = zero(x) <= x < Inf

mean(d::InvertedGamma) = d.shape > 1.0 ? (1.0 / d.scale) / (d.shape - 1.0) : Inf
var(d::InvertedGamma) = d.shape > 2.0 ? (1.0 / d.scale)^2 / ((d.shape - 1.0)^2 * (d.shape - 2.0)) : Inf
skewness(d::InvertedGamma) = d.shape > 3.0 ? (4.0 * sqrt(d.shape - 2.0)) / (d.shape - 3.0) : NaN
kurtosis(d::InvertedGamma) = d.shape > 4.0 ? (30.0 * d.shape - 66.0) / ((d.shape - 3.0) * (d.shape - 4.0)) : NaN

mode(d::InvertedGamma) = (1.0 / d.scale) / (d.shape + 1.0)
modes(d::InvertedGamma) = [mode(d)]

cdf(d::InvertedGamma, x::Real) = ccdf(Gamma(d.shape, d.scale), 1.0 / x)
ccdf(d::InvertedGamma, x::Real) = cdf(Gamma(d.shape, d.scale), 1.0 / x)
logcdf(d::InvertedGamma, x::Real) = logccdf(Gamma(d.shape, d.scale), 1.0 / x)
logccdf(d::InvertedGamma, x::Real) = logcdf(Gamma(d.shape, d.scale), 1.0 / x)

quantile(d::InvertedGamma, p::Real) = 1.0 / cquantile(Gamma(d.shape, d.scale), p)
cquantile(d::InvertedGamma, p::Real) = 1.0 / quantile(Gamma(d.shape, d.scale), p)
invlogcdf(d::InvertedGamma, p::Real) = 1.0 / invlogccdf(Gamma(d.shape, d.scale), p)
invlogccdf(d::InvertedGamma, p::Real) = 1.0 / invlogcdf(Gamma(d.shape, d.scale), p)

function entropy(d::InvertedGamma)
    a, b = d.shape, d.scale
    a - log(b) + lgamma(a) - (1.0 + a) * digamma(a)
end

function mgf(d::InvertedGamma, t::Real)
    a, b = d.shape, d.scale
    (2 * (-b * t)^(a / 2)) / gamma(a) * besselk(a, sqrt(-4.0 * b * t))
end

function cf(d::InvertedGamma, t::Real)
    a, b = d.shape, d.scale
    (2 * (-im * b * t)^(a / 2)) / gamma(a) * besselk(a, sqrt(-4.0 * im * b * t))
end

pdf(d::InvertedGamma, x::Real) = exp(logpdf(d, x))

function logpdf(d::InvertedGamma, x::Real)
    a, b = d.shape, d.scale
    -(a * log(b)) - lgamma(a) - ((a + 1.0) * log(x)) - 1.0/ (b * x)
end

rand(d::InvertedGamma) = 1.0 / rand(Gamma(d.shape, d.scale))

function rand!(d::InvertedGamma, A::Array{Float64})
    rand!(Gamma(d.shape, d.scale), A)
    for i in 1:length(A)
    	A[i] = 1.0 / A[i]
    end
    A
end

ext: jl


text:
using Test, LinearAlgebra
using CVChannel

@testset "./src/operations.jl" begin

@testset "isPPT" begin
    @test isPPT( Matrix(I,4,4), 2, [2,2] )
    @test !isPPT( [1 0 0 1;0 0 0 0;0 0 0 0;1 0 0 1], 2, [2,2] )
end

@testset "swapOperator" begin
    @test isapprox(swapOperator(2),[1 0 0 0 ; 0 0 1 0 ; 0 1 0 0 ; 0 0 0 1], atol = 1e-6)
    swap3 = [
        1 0 0 0 0 0 0 0 0 ;
        0 0 0 1 0 0 0 0 0 ;
        0 0 0 0 0 0 1 0 0 ;
        0 1 0 0 0 0 0 0 0 ;
        0 0 0 0 1 0 0 0 0 ;
        0 0 0 0 0 0 0 1 0 ;
        0 0 1 0 0 0 0 0 0 ;
        0 0 0 0 0 1 0 0 0 ;
        0 0 0 0 0 0 0 0 1
        ]
    @test isapprox(swapOperator(3),swap3, atol =1e-6)
end

@testset "permuteSubsystems" begin
    ρ = [1:1:4;];
    @test permuteSubsystems(ρ,[2,1],[2,2]) == [1;3;2;4]
    ρ = [1:1:8;]; target = [1;3;5;7;2;4;6;8];
    @test permuteSubsystems(ρ,[3,1,2],[2,2,2])==target
    ρ = copy(reshape([1:1:16;],(4,4))');
    target = [ 1 3 2 4 ; 9 11 10 12 ; 5 7 6 8; 13 15 14 16];
    @test permuteSubsystems(ρ,[2,1],[2,2])==target
    ρ = copy(reshape([1:1:64;],(8,8))')
    target = [
        1 5 2 6 3 7 4 8;
        33 37 34 38 35 39 36 40;
        9 13 10 14 11 15 12 16;
        41 45 42 46 43 47 44 48;
        17 21 18 22 19 23 20 24;
        49 53 50 54 51 55 52 56;
        25 29 26 30 27 31 28 32;
        57 61 58 62 59 63 60 64;
    ]
    @test permuteSubsystems(ρ,[2,3,1],[2,2,2])==target
    ρ = copy(reshape([1:1:144;],(12,12))')
    target = [
    1    3    5    7    9   11    2    4    6    8   10   12
    25   27   29   31   33   35   26   28   30   32   34   36
    49   51   53   55   57   59   50   52   54   56   58   60
    73   75   77   79   81   83   74   76   78   80   82   84
    97   99  101  103  105  107   98  100  102  104  106  108
    121  123  125  127  129  131  122  124  126  128  130  132
    13   15   17   19   21   23   14   16   18   20   22   24
    37   39   41   43   45   47   38   40   42   44   46   48
    61   63   65   67   69   71   62   64   66   68   70   72
    85   87   89   91   93   95   86   88   90   92   94   96
    109  111  113  115  117  119  110  112  114  116  118  120
    133  135  137  139  141  143  134  136  138  140  142  144
    ]
    @test permuteSubsystems(ρ,[3,1,2],[2,3,2]) == target
    @test_throws DomainError permuteSubsystems([1 2 3; 4 5 6], [2,3,1],[1,2,3])
end

@testset "shiftOperator" begin
    @test shiftOperator(2) == [0 1 ; 1 0]
    @test shiftOperator(3)*[0;1;0] == [0;0;1]
    @test shiftOperator(4)*[0;0;0;1] == [1;0;0;0]
end

@testset "discreteWeylOperator" begin
    #This guarantees it generates the Weyl basis in 2d
    @test discreteWeylOperator(0,0,2) == [1 0 ; 0 1]
    @test discreteWeylOperator(0,1,2) == [0 1 ; 1 0]
    @test isapprox(discreteWeylOperator(1,0,2),[1 0 ; 0 -1],atol=1e-6)
    @test isapprox(discreteWeylOperator(1,1,2), [0 -1 ; 1 0],atol=1e-6)
    #This checks that it is doing the right thing as the dimension grows
    @test discreteWeylOperator(0,0,3) == [1 0 0 ; 0 1 0 ; 0 0 1]
    @test discreteWeylOperator(0,2,3) == [0 1 0 ; 0 0 1 ; 1 0 0]
    @test isapprox(discreteWeylOperator(2,0,3),[1 0 0 ; 0 exp(pi*4im/3) 0 ; 0 0 exp(pi*8im/3)], atol=1e-6)
end
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
__precompile__(true)

module CLBlast

using Compat
using OpenCL.cl

depsfile = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if isfile(depsfile)
    include(depsfile)
else
    error("CLBlast not properly installed. Please run Pkg.build(\"CLBlast\") then restart Julia.")
end

include("constants.jl")
include("error.jl")
include("auxiliary_functions.jl")

include("L1/L1.jl")
include("L2/L2.jl")
include("L3/L3.jl")

end # module


text:
"""
    MSSSIM([kernel], [W]; num_scales = length(W)) <: FullReferenceIQI
    assess(iqi::MSSSIM, img, ref)
    assess_msssim(img, ref; num_scales = 5)

Computes the multi-scale structural similarity (MS-SSIM) between two images.

!!! tip

    The default parameters comes from [1]. For benchmark usage, it is recommended to not
    change the parameters, because most other SSIM implementations follows the same settings.

# Examples

For benchmark usage, it is recommended to use `assess_msssim(img, ref)`. One could also create
a custom `MSSSIM` instance and then pass it to `assess` or use it as a function.

```julia
iqi = MSSSIM(KernelFactors.gaussian(2.5, 17), (1.0, 1.0, 2.0)./4)
assess(iqi, img, ref)
iqi(img, ref) # both usages are equivalent
```

# Reference

[1] Wang, Z., Simoncelli, E. P., & Bovik, A. C. (2003, November). Multiscale structural similarity for image quality assessment.In *The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers*, 2003 (Vol. 2, pp. 1398-1402). Ieee.

[2] Wang, Z. IW-SSIM: Information Content Weighted Structural Similarity Index for Image Quality Assessment. Retrived July 9, 2020, from https://ece.uwaterloo.ca/~z70wang/research/iwssim/

[3] Jorge Pessoa. pytorch-msssim. Retrived July 9, 2020, from https://github.com/jorge-pessoa/pytorch-msssim

"""
struct MSSSIM{A, N} <: FullReferenceIQI
    kernel::A
    W::NTuple{N, NTuple{3, Float64}}

    function MSSSIM(kernel=SSIM_KERNEL, W=MSSSIM_W; num_scales::Integer=length(W))
        ndims(kernel) == 1 || throw(ArgumentError("only 1-d kernel is valid"))
        issymetric(kernel) || @warn "MSSSIM kernel is assumed to be symmetric"
        all(length.(W) .== 3) || throw(ArgumentError("(α, β, γ) required for all scales, instead it's $(W)"))
        (num_scales ∈ 1:length(W)) || throw(ArgumentError("`num_scales` $(num_scales) should be positive integer between [1, $(length(W))]"))
        all(x-> x>=0, flatten(W)) || throw(ArgumentError("α, β, γ should be non-negative for all scales, instead it's $(W)"))
        sum(flatten(W)) == 0 && throw(ArgumentError("MS-SSIM must have at least one weight > 0"))

        kernel = centered(kernel)
        if num_scales ≠ length(W)
            W ≠ MSSSIM_W && @warn "truncate MS-SSIM weights to scale $num_scales"
        end
        W = W[1:num_scales]
        sw = reduce((w1, w2) -> w1.+w2, W[1:num_scales])
        if !all(x -> x≈1.0, sw)
            W ≠ MSSSIM_W[1:num_scales] && @warn "normalize MS-SSIM weights so that (∑α, ∑β, ∑γ) == (1.0, 1.0, 1.0)"
            W = map(w->w./sw, W)
        end

        new{typeof(kernel), Int(num_scales)}(kernel, W)
    end
end
# Weights for α, β, γ in [1]
const MSSSIM_W = (
    (0.0448, 0.0448, 0.0448), # α₁, β₁, γ₁
    (0.2856, 0.2856, 0.2856), # α₂, β₂, γ₂
    (0.3001, 0.3001, 0.3001), # α₃, β₃, γ₃
    (0.2363, 0.2363, 0.2363), # α₄, β₄, γ₄
    (0.1333, 0.1333, 0.1333), # α₅, β₅, γ₅
)
# shorthand for αᵢ=βᵢ=γᵢ for all scales i
function MSSSIM(kernel, W::NTuple{N, <:Real}; num_scales=length(W)) where N
    MSSSIM(kernel, map(x->(x, x, x), W); num_scales=num_scales)
end

Base.:(==)(ia::MSSSIM, ib::MSSSIM) = ia.kernel == ib.kernel && ia.W == ib.W

# SSIM does not allow for user specifying peakval and K, so we don't allow it here either
(iqi::MSSSIM)(x, ref) = _msssim(iqi, x, ref)
assess_msssim(x, ref; num_scales=5) = MSSSIM(num_scales=num_scales)(x, ref)

# Implementation details
function _msssim(iqi::MSSSIM,
                 x::Union{GenericGrayImage, AbstractArray{<:AbstractRGB}},
                 ref::Union{GenericGrayImage, AbstractArray{<:AbstractRGB}};
                 peakval=1.0,
                 K=SSIM_K)
    if size(x) ≠ size(ref)
        err = ArgumentError("images should be the same size, instead they're $(size(x))-$(size(ref))")
        throw(err)
    end
    if axes(x) ≠ axes(ref)
        x = OffsetArrays.no_offset_view(x)
        ref = OffsetArrays.no_offset_view(ref)
    end
    C₁, C₂ = @. (peakval * K)^2
    C₃ = C₂/2

    num_scales = length(iqi.W)
    T = promote_type(float(eltype(ref)), float(eltype(x)))
    x = of_eltype(T, x)
    ref = of_eltype(T, ref)

    mean_lcs = NTuple{3, Float64}[]
    for i in 1:num_scales
        # instead of the original implementation that calculates mean(c .* s)
        # here we use a more general version mean(c) * mean(s)
        lcs = mean.(__ssim_map_general(x, ref, iqi.kernel, C₁, C₂, C₃; crop=true))
        push!(mean_lcs, lcs)

        # downsampling
        i == num_scales && break
        x = _average_pooling(x)
        ref = _average_pooling(ref)
    end

    # weighted them up, but only l of the last scale is used according to equation (7) in [1]
    mean_lcs = map(mean_lcs, iqi.W) do lcs, w
        # similar to SSIM, here we ensure that negative numbers in s are not being raised to powers
        # less than 1
        s = w[3] < 1 ? max(0, lcs[3]) : lcs[3]
        (lcs[1], lcs[2], s) .^ w
    end
    if length(mean_lcs) == 1
        return prod(first(mean_lcs))
    else
        return mapreduce(x->x[2]*x[3], *, mean_lcs[1:end-1]) * prod(mean_lcs[end])
    end
end

_msssim(iqi::MSSSIM, x::AbstractArray{<:Color3}, ref::AbstractArray{<:Color3}) =
    _msssim(iqi, of_eltype(RGB, x), of_eltype(RGB, ref))

function _average_pooling(x::GenericImage; kernel=ones(2)./2)
    # note that this is slightly slower than two-fold method `ImageTransformations.restrict`
    window = kernelfactors(Tuple(repeated(kernel, ndims(x))))
    R = ntuple(i->first(axes(x, i)) : 2 : last(axes(x, i)), ndims(x))
    return imfilter(x, window, "symmetric")[R...]
end


text:
import ..Abstract.AbstractMolecularContainer

"""
A ``Atom`` (immutable) object stores all the information
about an atom.
The coordinates are stored as a 2D matrix (N by 3),
which each row represents a trajectory frame and
each column represents one Cartesian dimension.
+--------------------+-------------------------+-----------------------------------+
| field              | type                    | note                              |
+====================+=========================+===================================+
| :atom_id           | Integer                 | atom ID                           |
+--------------------+-------------------------+-----------------------------------+
| :residue_name      | AbstractString          | residue name                      |
+--------------------+-------------------------+-----------------------------------+
| :residue_id        | Integer                 | residue ID                        |
+--------------------+-------------------------+-----------------------------------+
| :segment_name      | AbstractString          | segment name                      |
+--------------------+-------------------------+-----------------------------------+
| :chain_name        | AbstractString          | chain name                        |
+--------------------+-------------------------+-----------------------------------+
| :system_name       | AbstractString          | molecular system name             |
+--------------------+-------------------------+-----------------------------------+
| :atom_name         | AbstractString          | atom name (no extra spaces)       |
+--------------------+-------------------------+-----------------------------------+
| :atom_type         | AbstractString          | atom type                         |
+--------------------+-------------------------+-----------------------------------+
| :raw_atom_name     | AbstractString          | atom name with spaces preserved   |
+--------------------+-------------------------+-----------------------------------+
| :alt_location      | Char                    | alternative location              |
+--------------------+-------------------------+-----------------------------------+
| :insertion         | Char                    | insertion code                    |
+--------------------+-------------------------+-----------------------------------+
| :occupancy         | AbstractFloat           | atomic occupancy                  |
+--------------------+-------------------------+-----------------------------------+
| :beta              | AbstractFloat           | temperature factor                |
+--------------------+-------------------------+-----------------------------------+
| :element           | AbstractString          | chemical element                  |
+--------------------+-------------------------+-----------------------------------+
| :charge            | Integer                 | atomic charge                     |
+--------------------+-------------------------+-----------------------------------+
| :coordinate        | Array{AbstractFloat,1}  | atomic coordinate                 |
+--------------------+-------------------------+-----------------------------------+
| :topology          | TopologicAtom           | topological information           |
+--------------------+-------------------------+-----------------------------------+

"""

immutable Atom <: AbstractMolecularContainer
    index::Integer # sequential index
    atom_name::AbstractString # name without spaces
    atom_type::AbstractString # atom type name
    alt_location::Char # 17
    insertion::Char # 27
    occupancy::AbstractFloat # 55-60
    beta::AbstractFloat # 61-66
    element::AbstractString # 77-78
    charge::Integer # 79-80

    residue_name::AbstractString
    residue_id::Integer
    segment_name::AbstractString
    chain_name::AbstractString
    system_name::AbstractString

    coordinate::Array{AbstractFloat,1}
end

# constructor with keyword arguments
include("constructor.jl")


text:
using UncertainData
import Interpolations: Linear, Flat, OnGrid
using Interpolations

# Define the InterpolateAndBin instance
intp_bin = InterpolateAndBin(mean, 0:10:100, Linear(), 0:0.01:100, Flat(OnGrid()))

ib_test = InterpolateBinTest(CrossMappingTest(), intp_bin)
@test ib_test isa InterpolateBinTest
@test ib_test.n == 1

ib_test = InterpolateBinTest(CrossMappingTest(), intp_bin, 10)
@test ib_test isa InterpolateBinTest
@test ib_test.n == 10

# Some test data
N = 100
sys = ar1_unidir(c_xy = 0.2)
X, Y = example_uncertain_indexvalue_datasets(sys, N, (1, 2),
    d_xval = Uniform(0.001, 0.2), d_yval = Uniform(0.001, 0.2));

# Define a causality test
k, l, m = 1, 1, 1 # embedding parameters
n_subdivisions = floor(Int, N^(1/(k + l + m + 1)))
state_space_binning = RectangularBinning(n_subdivisions)

ηmax = 3
ηs = -ηmax:ηmax
te_test = VisitationFrequencyTest(k = k, l = l, m = m,
            binning = state_space_binning,
            ηs = ηs, estimator = VisitationFrequency(b = 2)) # use base-2 logarithms
pa_test = PredictiveAsymmetryTest(predictive_test = te_test)

# Define interpolation grid over the range of available index values
tmin = max(minimum(mean.(X.indices)), minimum(mean.(Y.indices)))
tmax = max(maximum(mean.(X.indices)), maximum(mean.(Y.indices)))
intp_grid = tmin:0.01:tmax

# Define binning grid
left_bin_edges = tmin:1.5:tmax

# Define the Interpolate AndBin instance
intp_bin = InterpolateAndBin(mean, left_bin_edges, Linear(), intp_grid, Flat(OnGrid()))

# Define interpolate-and-bin test
ib_test = InterpolateBinTest(pa_test, intp_bin, 2)

res_xy = causality(X, Y, ib_test)

@test res_xy isa Vector{Vector{T}} where T
@test length(res_xy) == ib_test.n


text:
using Parameters

abstract type AbstractHenryLawConstant end
export AbstractHenryLawConstant

struct EmptyHenryLawConstant <: AbstractHenryLawConstant end
(eh::EmptyHenryLawConstant)(;T::N) where {N<:Number} = Inf
export EmptyHenryLawConstant

@with_kw struct ConstantHenryLawConstant{N<:Number} <: AbstractHenryLawConstant
    kH::N
end
(ch::ConstantHenryLawConstant)(;T::N) where {N<:Number} = ch.kH
export ConstantHenryLawConstant

struct TemperatureDependentHenryLawConstant{N<:Function} <: AbstractHenryLawConstant
    kH::N
end

function TemperatureDependentHenryLawConstant(;Ts::N1,kHs::N2) where {N1<:AbstractArray,N2<:AbstractArray}
    return TemperatureDependentHenryLawConstant(getspline(Ts,kHs))
end

(th::TemperatureDependentHenryLawConstant)(;T::N) where {N<:Number} = th.kH(T)
export TemperatureDependentHenryLawConstant

abstract type AbstractLiquidVolumetricMassTransferCoefficient end
export AbstractLiquidVolumetricMassTransferCoefficient

struct EmptyLiquidVolumetricMassTransferCoefficient <: AbstractLiquidVolumetricMassTransferCoefficient end
(eh::EmptyLiquidVolumetricMassTransferCoefficient)(;T::N) where {N<:Number} = 0
export EmptyLiquidVolumetricMassTransferCoefficient

@with_kw struct ConstantLiquidVolumetricMassTransferCoefficient{N<:Number} <: AbstractLiquidVolumetricMassTransferCoefficient
    kLA::N
end
(cv::ConstantLiquidVolumetricMassTransferCoefficient)(;T::N) where {N<:Number} = cv.kLA
export ConstantLiquidVolumetricMassTransferCoefficient

struct TemperatureDependentLiquidVolumetricMassTransferCoefficient{N<:Function} <: AbstractLiquidVolumetricMassTransferCoefficient
    kLA::N
end

function TemperatureDependentLiquidVolumetricMassTransferCoefficient(;Ts::N1,kLAs::N2) where {N1<:AbstractArray,N2<:AbstractArray}
    return TemperatureDependentLiquidVolumetricMassTransferCoefficient(getspline(Ts,kLAs))
end

(tv::TemperatureDependentLiquidVolumetricMassTransferCoefficient)(;T::N) where {N<:Number} = tv.kLA(T)
export TemperatureDependentLiquidVolumetricMassTransferCoefficient


text:
const Nspecies = 47 # number of complexes in surface + endosome + free ligand
const halfL = 19 # number of complexes on surface alone
const internalFrac = 0.5 # Same as that used in TAM model
const recIDX = [1, 2, 3, 10, 17]
const recIDXint = [ii + halfL for ii in recIDX]
const ligIDX = [39, 40, 41]
const activeSpec = [8, 9, 15, 16, 19]
const pSTATidx = [43, 44, 45]

const Nparams = 36 # number of unknowns for the full model
const kfbnd = 0.60 # Assuming on rate of 10^7 M-1 sec-1
const internalV = 623.0 # Same as that used in TAM model

# p[1:4] is kfwd, k1rev, k2rev, k4rev
# p[5:8] is k5rev, k10rev, k11rev, k13rev
# p[9:12] is k14rev, k16rev, k17rev, k22rev
# p[13:14] is k23rev, k24rev

function dYdT(du, u, p, ILs)
    # IL2Ra, IL2Rb, gc, IL2_IL2Ra, IL2_IL2Rb, IL2_IL2Ra_IL2Rb, IL2_IL2Ra_gc, IL2_IL2Rb_gc, IL2_IL2Ra_IL2Rb_gc
    # IL15Ra, IL15_IL15Ra, IL15_IL2Rb, IL15_IL15Ra_IL2Rb, IL15_IL15Ra_gc, IL15_IL2Rb_gc, IL15_IL15Ra_IL2Rb_gc

    # IL2/15
    for i = 0:1
        pp = view(p, (1 + i * 6):(8 + i * 6))
        k12rev = pp[2] * pp[7] / pp[3] # Detailed balance
        k8rev = pp[6] * k12rev / pp[5] # Detailed balance
        k9rev = pp[6] * pp[7] / pp[4] # Detailed balance

        du[2] +=
            -kfbnd * u[2] * ILs[i + 1] + pp[3] * u[5 + 7 * i] - p[1] * u[2] * u[7 + 7 * i] + k9rev * u[9 + 7 * i] - p[1] * u[2] * u[4 + 7 * i] +
            pp[7] * u[6 + 7 * i]
        du[3] +=
            -p[1] * u[5 + 7 * i] * u[3] + p[5 + 7 * i] * u[8 + 7 * i] - p[1] * u[4 + 7 * i] * u[3] + p[4 + 7 * i] * u[7 + 7 * i] -
            p[1] * u[6 + 7 * i] * u[3] + p[6 + 7 * i] * u[9 + 7 * i]
        du[1 + 9 * i] =
            -kfbnd * u[1 + 9 * i] * ILs[i + 1] + pp[2] * u[4 + 7 * i] - p[1] * u[1 + 9 * i] * u[8 + 7 * i] + k8rev * u[9 + 7 * i] -
            p[1] * u[1 + 9 * i] * u[5 + 7 * i] + k12rev * u[6 + 7 * i]
        du[4 + 7 * i] =
            -p[1] * u[4 + 7 * i] * u[2] + pp[7] * u[6 + 7 * i] - p[1] * u[4 + 7 * i] * u[3] +
            pp[4] * u[7 + 7 * i] +
            kfbnd * ILs[i + 1] * u[1 + 9 * i] - pp[2] * u[4 + 7 * i]
        du[5 + 7 * i] =
            -p[1] * u[5 + 7 * i] * u[1 + 9 * i] + k12rev * u[6 + 7 * i] - p[1] * u[5 + 7 * i] * u[3] +
            p[5 + 7 * i] * u[8 + 7 * i] +
            kfbnd * ILs[i + 1] * u[2] - pp[3] * u[5 + 7 * i]
        du[6 + 7 * i] =
            -p[1] * u[6 + 7 * i] * u[3] + pp[6] * u[9 + 7 * i] + p[1] * u[4 + 7 * i] * u[2] - pp[7] * u[6 + 7 * i] +
            p[1] * u[5 + 7 * i] * u[1 + 9 * i] - k12rev * u[6 + 7 * i]
        du[9 + 7 * i] =
            p[1] * u[8 + 7 * i] * u[1 + 9 * i] - k8rev * u[9 + 7 * i] + p[1] * u[7 + 7 * i] * u[2] - k9rev * u[9 + 7 * i] +
            p[1] * u[6 + 7 * i] * u[3] - pp[6] * u[9 + 7 * i]
        du[7 + 7 * i] = -p[1] * u[7 + 7 * i] * u[2] + k9rev * u[9 + 7 * i] + p[1] * u[4 + 7 * i] * u[3] - pp[4] * u[7 + 7 * i]
        du[8 + 7 * i] = -p[1] * u[8 + 7 * i] * u[1 + 9 * i] + k8rev * u[9 + 7 * i] + p[1] * u[3] * u[5 + 7 * i] - pp[5] * u[8 + 7 * i]
    end

    du[3] += -p[1] * u[3] * u[18] + p[15] * u[19]
    du[17] = -kfbnd * u[17] * ILs[3] + p[14] * u[18]
    du[18] = kfbnd * u[17] * ILs[3] - p[14] * u[18] - p[1] * u[3] * u[18] + p[15] * u[19]
    du[19] = p[1] * u[3] * u[18] - p[15] * u[19]

    return nothing
end

function fullDeriv(du, u, p, t)
    fill!(du, 0.0)

    ILs = view(p, 1:3)
    pSurf = view(p, 4:18)
    pEndo = p[19]
    trafP = view(p, 20:30)

    # Calculate cell surface and endosomal reactions
    dYdT(du, u, pSurf, ILs)

    # Add STAT5 reactions
    activR = sum(u[activeSpec]) + internalFrac * sum(u[activeSpec .+ halfL])
    du[42] = p[36] * u[47] - p[31] * u[42] * activR # STAT5
    du[43] = p[31] * u[42] * activR - p[32] * u[43]^2 # pSTAT5
    du[44] = 0.5 * p[32] * u[43]^2 - p[33] * u[44] # pSTAT5d
    du[45] = p[33] * u[44] - p[34] * u[45] # pSTAT5nd
    du[46] = p[34] * u[45] - 0.5 * p[35] * u[46] # STAT5nd
    du[47] = p[35] * u[46] - p[36] * u[47] # STAT5n

    dYdT(view(du, (halfL + 1):(2 * halfL)), view(u, (halfL + 1):(2 * halfL)), pEndo * pSurf, view(u, ligIDX))

    # Handle endosomal ligand balance.
    # Must come before trafficking as we only calculate this based on reactions balance
    du[39] = -sum(view(du, (halfL + 4):(halfL + 9))) / internalV
    du[40] = -sum(view(du, (halfL + 11):(halfL + 16))) / internalV
    du[41] = -sum(view(du, (halfL + 18):(halfL + 19))) / internalV

    # Actually calculate the trafficking
    for ii = 1:halfL
        if ii ∈ activeSpec
            du[ii] -= u[ii] * (trafP[1] + trafP[2]) # Endo
            du[ii + halfL] += u[ii] * (trafP[1] + trafP[2]) / internalFrac - trafP[5] * u[ii + halfL] # Endo, deg
        else
            du[ii] += -u[ii] * trafP[1] + trafP[4] * (1.0 - trafP[3]) * u[ii + halfL] * internalFrac # Endo, recycling
            du[ii + halfL] += u[ii] * trafP[1] / internalFrac - trafP[4] * (1.0 - trafP[3]) * u[ii + halfL] - (trafP[5] * trafP[3]) * u[ii + halfL] # Endo, rec, deg
        end
    end

    # Expression: IL2Ra, IL2Rb, gc, IL15Ra, IL7Ra
    du[recIDX] .+= view(trafP, 6:10)

    # Degradation does lead to some clearance of ligand in the endosome
    du[ligIDX] .-= view(u, ligIDX) * trafP[5]

    return nothing
end

" Takes in parameters and solves for steady state expression and initial species states "
function solveAutocrine(rIn::Vector)
    @assert all(rIn .>= 0.0)
    r = view(rIn, 20:30)
    @assert r[3] < 1.0

    # r is endo, activeEndo, sortF, kRec, kDeg, Rexpr*8
    y0 = zeros(eltype(r), Nspecies)

    # Add base amount of STAT5
    y0[42] = r[11]

    # Expand out trafficking terms
    kRec = r[4] * (1.0 - r[3])
    kDeg = r[5] * r[3]

    # Assuming no autocrine ligand, so can solve steady state
    # Add the species
    y0[recIDXint] = r[6:10] / kDeg / internalFrac
    y0[recIDX] = (r[6:10] + kRec * y0[recIDXint] * internalFrac) / r[1]

    return y0
end


text:
const bondoperators = ['@', '-', '\\', '/', '=', '#', '\$']
const specialoperators = [ '.', ']', '[', '(', ')']
const operators = [ specialoperators;  bondoperators ]

const bonds = Dict( '-' => 1, '\\' => 1, '/' => 1, '=' => 2, '#' => 3, '\$' => 4 )
const valence = Dict( "B" => 3, "C" => 4, "N" => 3, "O" => 2, "P" => 3,
                      "S" => 2, "F" => 1, "Cl" => 1, "Br" => 1, "I" => 1,
                      "H" => 1, "Po" => 4, "As" => 5, "Co" => 6)

const aliphatics  = ["B", "C", "N", "O", "S", "P", "F", "Cl", "Br", "I"]
const aromatics   = ["b", "c", "n", "o", "s", "p"]

const bracket = [ "H", "He" ,"Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na",
                  "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr",
                  "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br",
                  "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd",
                  "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "Hf",
                  "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi",
                  "Po", "At", "Rn", "Fr", "Ra", "Rf", "Db", "Sg", "Bh", "Hs", "Mt",
                  "Ds", "Rg", "Cn", "Fl", "Lv", "La", "Ce", "Pr", "Nd", "Pm", "Sm",
                  "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Ac", "Th",
                  "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md",
                  "No", "Lr", "Mg"]

const bracket_aromatic = ["b", "c", "n", "o", "p", "s", "se", "as"]

isoperator(x::AbstractChar) = x ∈ operators
isbondoperator(x::AbstractChar) = x ∈ bondoperators
isspecialoperator(x::AbstractChar) = x ∈ specialoperators
ispm(x::AbstractChar) = x ∈ ('+', '-')
FindNumerics( S::Vector{Char} ) = findall( isnumeric.( S ) )
FindPMs( S::Vector{Char} ) = findall( ispm.( S ) )

mutable struct Element
    symbol::String
    isotope::Union{Int16, Nothing}
    aromatic::Bool
    ringID::Array{Int16,1}
    explicithydrogens::Int8
    implicithydrogens::Int8
    charge::Int8
end

function Element(symbol::String)
    if islowercase(symbol[1])
        return Element(uppercasefirst(symbol), nothing, true, [], 0, 0, 0 )
    else
        return Element(symbol, nothing, false, [], 0, 0, 0 )
    end
end

Element(e::Element) = Element(e.symbol, e.isotope, e.aromatic, copy(e.ringID), e.explicithydrogens, e.implicithydrogens, e.charge)

function Base.:(==)(a::Element, b::Element)
    return all( [ a.symbol == b.symbol, a.isotope == b.isotope,
                a.aromatic == b.aromatic, a.ringID == b.ringID,
                a.explicithydrogens == b.explicithydrogens,
                a.charge == b.charge])
end

Base.copy(e::Element) = Element(e)

#immutable struct to hold elemental information - less memory then "Element"
struct GraphElement
    symbol::String
    isotope::Union{Int16, Nothing}
    aromatic::Bool
    hydrogens::Int8
    charge::Int8
end

GraphElement(E::Element) = GraphElement(E.symbol, E.isotope, E.aromatic,
                                        E.hydrogens, E.charge)

abbreviation(E::Union{GraphElement,Element}) = E.symbol
implicitH(E::Element) = E.implicithydrogens
explicitH(E::Element) = E.explicithydrogens
H(E::Element) = E.implicithydrogens + E.explicithydrogens
Hydrogens(E::Element) = E.implicithydrogens + E.explicithydrogens
charge(E::Union{GraphElement,Element}) = E.charge
isotope(E::Union{GraphElement,Element}) = E.isotope


text:
# # Visual

# The function `visual` can be used to give data-independent visual information about the plot
# (plotting function or attributes).

# The available plotting functions are documented [here](http://makie.juliaplots.org/dev/examples/plotting_functions/).
# Refer to plotting functions using upper CamelCase for `visual`'s first argument (e.g. `visual(Scatter), visual(BarPlot)`).
# See the documentation of each plotting function to discover the available attributes.

## Examples

using AlgebraOfGraphics, CairoMakie
set_aog_theme!()

df = (x=randn(1000), y=randn(1000))
plt = data(df) * mapping(:x, :y) * AlgebraOfGraphics.density(npoints=50)
draw(plt * visual(Heatmap)) # plot as heatmap (the default)

#

draw(plt * visual(colormap=:viridis)) # set a different colormap

#

draw(plt * visual(Contour)) # plot as contour

#

draw(plt * visual(Contour, linewidth=2)) # plot as contour with thicker lines


text:
using Revise
using AdFem
using ADCME
using Random
using Mamba
using ProgressMeter
using PyPlot
using DelimitedFiles

n = 10
m = 15
h = 0.1

σ0 = 0.1
N = 500

if length(ARGS)==2
    σ0 = parse(Float64, ARGS[1])
    N = parse(Int64, ARGS[2])
end
DIR = "sigma$(σ0)-$N-$(randstring(10))"
if !isdir(DIR)
    mkdir(DIR)
end

Random.seed!(233)
obs_id = rand(interior_node("all", m, n, h), 10)
obs_id = [obs_id; obs_id .+ (m+1)*(n+1)]
function solve_elasticity(c)
    c = constant(c)
    E, ν = c[1], c[2]
    H = 1e6*E/(1+ν)/(1-2ν)*tensor([
             1-ν ν 0.0
             ν 1-ν 0.0
             0.0 0.0 (1-2ν)/2
         ])
    K = compute_fem_stiffness_matrix(H, m, n, h)
    bdedge = bcedge("right", m, n, h)
    bdnode = bcnode("left", m, n, h)
    K, _ = fem_impose_Dirichlet_boundary_condition(K, bdnode, m, n, h)
    F = compute_fem_traction_term([zeros(4*m*n) 1e6*ones(4*m*n)], bdedge, m, n, h)
    sol = K\F
    sol = sol[obs_id]
    return sol
end

sol = solve_elasticity([2.0;0.35])
sess = Session(); init(sess)
sol_ = run(sess, sol)
# # visualize_displacement(sol_*1e3, m, n, h)

# error()
function invsigma(x)
    return log(x/(1-x))
end

c = placeholder(rand(2))
sol = solve_elasticity(c)

function visualize(sim)
    close("all")
    figure(figsize=(10, 5))
    subplot(121)
    title("E")
    p = hist(sim[:,1], bins=50, density=true)
    PyPlot.plot(ones(100)*2.0, LinRange(0, maximum(p[1]), 100), "--")
    subplot(122)
    title("\$\\mu\$")
    p = hist(sim[:,2], bins=50, density=true)
    PyPlot.plot(ones(100)*0.35, LinRange(0, maximum(p[1]), 100), "--")
    savefig("$DIR/hist.png")
    writedlm("$DIR/data.txt", sim)
end

#------------------------------------------------------------------------------------------
# RWMC
function logf(x)
    b0 = sigmoid(x[1])*10
    b1 = sigmoid(x[2])*0.499
    r = sol_ - run(sess, sol, c=>[b0;b1])
    -sum(r.^2)/(2 * σ0^2)
end
sess = Session(); init(sess)

burnin = div(N, 5)
sim = Chains(N, 2, names=["E", "mu"])
θ = RWMVariate([invsigma(2.0/10),invsigma(0.35/0.499)], [0.1,0.1], logf, proposal = SymUniform)

@showprogress for i = 1:N
    sample!(θ)
    sim[i, :, 1] = [ sigmoid(θ[1])*10; sigmoid(θ[2])*0.499]
end
visualize(sim.value[burnin+1:end,:,1])

# #------------------------------------------------------------------------------------------
# # HMC

# r = sol_ - sol
# lf = -sum(r^2)/(2σ0^2)
# lfg = gradients(lf, c)

# function logfgrad(x)
#     b0 = sigmoid(x[1])*10
#     b1 = sigmoid(x[2])*0.499
#     r, dr = run(sess, [lf,lfg], c=>[b0;b1])
#     r, dr
# end
# sess = Session(); init(sess)
# N = 500
# burnin = 100
# sim = Chains(N, 2, names=["E", "mu"])

# ε = 0.1
# L = 50
# θ = HMCVariate([invsigma(2.0/10),invsigma(0.35/0.499)],
#     ε, L, logfgrad)

# @showprogress for i = 1:N
#     sample!(θ)
#     sim[i, :, 1] = [ sigmoid(θ[1])*10; sigmoid(θ[2])*0.499]
# end

# describe(sim)

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
#Copyright (c) 2015: Julia Computing Inc. All rights reserved.
module JuliaInXL

export xldate, parse_and_eval, jlsetvar
#import Base.include_from_node1; export include_from_node1

using Reexport
@reexport using JuliaWebAPI
using Logging
using Dates
using ZMQ

global const SECONDS_PER_MINUTE = 60
global const MINUTES_PER_HOUR = 60
global const HOURS_PER_DAY = 24
global const SECONDS_PER_DAY = (HOURS_PER_DAY * MINUTES_PER_HOUR * SECONDS_PER_MINUTE)
global const DAY_MILLISECONDS = SECONDS_PER_DAY * 1000
#global const threadid = Vector{Int}(128)
global const threadid = Vector(undef,128)

#The xldata function is adapted from the code in the Apache POI project
function xldate(date::Number; use1904windowing=false, roundtoSeconds=false)
  wholeDays = floor(Int, date)
  millisInDay = round(Int, (date-wholeDays)*DAY_MILLISECONDS)
  startYear = 1900
  dayAdjust = -1 #Excel thinks 2/29/1900 is a valid date, which it isn't
  if (use1904windowing)
    startYear = 1904
    dayAdjust = 1 #// 1904 date windowing uses 1/2/1904 as the first day
  elseif (wholeDays < 61)
    dayAdjust = 0
  end

  d = DateTime(startYear, 1, 1)
  d = d+Day(wholeDays + dayAdjust - 1)
  if roundtoSeconds
    millisInDay = round(Int, millisInDay/1000)*1000
  end
  d = d + Millisecond(millisInDay)
end

parse_and_eval(arg::String) = Main.eval(Meta.parse(arg))

function jlsetvar(name::String, arg)
    y = Any[arg]
    JuliaWebAPI.narrow_args!(y)
    Main.eval(Main.eval(Meta.parse(":(global $name = $(y[1]))")))
    "Set global variable $name"
end

# Overriding base include method to avoid serializing issue
# Issue is `include` returns the last thing that it encounters in the file. Which may be something that is not serializable. To avoid the error, we add a `nothing` at the end

function include(x)
    Base.MainInclude.include(x)
    nothing
end

# entry point for new thread
function heartbeat_thread(sock::Ptr{Nothing})
    ccall((:zmq_device,ZMQ.zmq), Cint, (Cint, Ptr{Nothing}, Ptr{Nothing}),
          ZMQ.QUEUE, sock, sock)
    nothing
end

function start_heartbeat(port=9998)
	ctx=Context()
	sock=Socket(ctx, REP)
	ZMQ.bind(sock, "tcp://*:$port")
    heartbeat_c = cfunction(heartbeat_thread, Nothing, (Ptr{Nothing},))
    ccall(:uv_thread_create, Cint, (Ptr{Int}, Ptr{Nothing}, Ptr{Nothing}),
          threadid, heartbeat_c, sock.data)
end

function start_async_server(port=9999)
  global_logger(SimpleLogger(stdout, Logging.Warn))
  transport = ZMQTransport("tcp://127.0.0.1:$port", REP, true)
  msgformat = JSONMsgFormat()
  api = APIResponder(transport, msgformat, nothing, true)
  map(fn->register(api, fn), [+, -, *, /, include, parse_and_eval, jlsetvar])
  conn = process(api; async=true)

  Base.eval(Main, :(conn=$conn))
  #start_heartbeat(port-1)
  println("JuliaInXL server bound to variable 'conn' at tcp://127.0.0.1:$port")
  println("JuliaInXL server Connection accessible from tcp://localhost:$port inside of Excel")
end

end # module


text:
module HMCExamples
using Base: String
using StatsPlots
using Turing, TuringCallbacks, TensorBoardLogger, ArgParse, CSV, DataFrames, JSON, LinearAlgebra, DifferentiableStateSpaceModels, Random, DelimitedFiles, Logging
import ArgParse.parse_item

# auto-generated code from /deps/generate_models.jl
include("generated_models/rbc_1.jl")
include("generated_models/rbc_2.jl")
include("generated_models/FVGQ20_1.jl")
include("generated_models/FVGQ20_2.jl")

# Utilities
include("utilities.jl")
include("experiment_results.jl")

# Specific Models and CLI drivers
include("turing_models.jl")
include("estimate_rbc_1_kalman.jl")
# include("estimate_rbc_1_joint.jl")
include("estimate_rbc_2_joint.jl")
include("estimate_FVGQ_1_kalman.jl")
include("estimate_FVGQ_1_joint.jl")
include("estimate_FVGQ_2_joint.jl")

end #module


text:
# Autogenerated wrapper script for FastHadamardStructuredTransforms_jll for x86_64-apple-darwin
export libfasttransforms

JLLWrappers.@generate_wrapper_header("FastHadamardStructuredTransforms")
JLLWrappers.@declare_library_product(libfasttransforms, "@rpath/libfasttransforms.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libfasttransforms,
        "lib/libfasttransforms.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#Returns the path as determined by llvm-config, or throws an error if LLVM
#could not be found
function find_llvm()
    if haskey(ENV, "TRAVIS")
        return readchomp(`which llvm-config-3.4`)
    end
    tmp = joinpath(JULIA_HOME, "llvm-config")
    isfile(tmp) && return tmp
    #Get path from llvm-config
    try
        tmp = readchomp(`which llvm-config`)
    end
    isfile(tmp) && return tmp
    error("No useable llvm-config found")
end

#On OSX, install llvm automatically through Homebrew.jl if no llvm installation
#is found
@static if is_apple()
    using Homebrew

    try
        find_llvm()
    catch #Not installed on system
        info("Installing llvm through Homebrew.jl")
        if !Homebrew.installed("llvm")
            Homebrew.add("llvm")
        end
    end
end

if !haskey(ENV, "LLVM_CONFIG")
    ENV["LLVM_CONFIG"] = find_llvm()
end

# on travis and some other systems, the llvm lib may not be in the
# default path.
push!(Libdl.DL_LOAD_PATH, readchomp(`$(ENV["LLVM_CONFIG"]) --libdir`))

cd(joinpath(dirname(@__FILE__), "src"))
try
    run(`make`)
catch
    has_build_clang_flag=false
    julia_make_user = joinpath(JULIA_HOME, "..", "..", "Make.user")

    if isfile(julia_make_user)
        open(julia_make_user) do f
            for l in eachline(f)
                if startswith(l, "BUILD_LLVM_CLANG=1")
                    has_build_clang_flag = true
                    break
                end
            end
        end
    end

    if !has_build_clang_flag
        error("""clang headers could not be found and BUILD_LLVM_CLANG=1 not specified in Make.user
Try setting BUILD_LLVM_CLANG=1 in Make.user in $JULIA_HOME and try again.""")
    else
        exit()
    end
end

if (!ispath("../usr"))
    run(`mkdir ../usr`)
end
if (!ispath("../usr/lib"))
    run(`mkdir ../usr/lib`)
end

run(`mv libwrapclang.$(Libdl.dlext) ../usr/lib`)


text:
function show{FSA <: FixedVector}(io::IO, a::Vector{FSA})
	println(io, length(a), "-element", typeof(a), ":")
	for elem in a
		print(io, " ")
		showcompact(io, elem)
		println(io)
	end
end
immutable MaxFun <: Func{2} end
immutable MinFun <: Func{2} end
call(::MaxFun, a, b) = max(a, b)
call(::MinFun, a, b) = min(a, b)
minimum{T <: FixedArray}(a::Vector{T}) = reduce(MinFun(), a)
maximum{T <: FixedArray}(a::Vector{T}) = reduce(MaxFun(), a)
function isapprox{FSA <: FixedArray}(a::FSA, b::Array)
    for i=1:length(a)
        !isapprox(a[i], b[i]) && return false
    end
    true
end

(.+){T<:FixedArray, ND}(a::Array{T, ND}, x::T) = T[elem .+ x for elem in a]
(+){T<:FixedArray, ND}(a::Array{T, ND}, x::T) = a .+ x


text:
@testset "FairEvaluate Benchmarking" begin
	model = ConstantClassifier()
	X, y = @load_toydata
	results = fairevaluate([model, model], X, y, grp=:Sex, priv_grps=["M"])
	@test length(results)==5
	@test size(results["measures"])==(1,)
	@test size(results["pvalues"])==(1, 2, 2)
	@test size(results["classifier_names"])==(2,)
	@test size(results["results"])==(1, 2, 6)
	@test size(results["tstats"])==(1, 2, 2)
end


text:
@test_throws MethodError GtkBuilderAid.arguments(5)
@test_throws ErrorException GtkBuilderAid.arguments(:(hello::There))

nullio = IOBuffer(0)

function test_macro_throws(error_type, macroexpr)
  expansion = macroexpand(macroexpr)
  if expansion.head != :error
    warn("Expansion did not throw error")
    println(expansion)
    @test false
  end
  if !isa(expansion.args[1], error_type)
    warn("Expansion threw wrong error")
    warn("Expected: $error_type")
    warn("Got: $(expansion.args[1])")
    dump(expansion)
    @test false
  end
end

test_macro_throws(TypeError, quote
@GtkBuilderAid helloSymbol
end)

test_macro_throws(ArgumentError, quote
@GtkBuilderAid helloExpr()
end)

test_macro_throws(ArgumentError, quote
@GtkBuilderAid
end)

test_macro_throws(ErrorException, quote
@GtkBuilderAid symbolDirective begin
end
end)

test_macro_throws(ErrorException, quote
@GtkBuilderAid badexpr::Directive begin
end
end)

test_app = @GtkApplication("com.github.test_gtkbuilderaid", 0)

# Check that an empty builder doesn't crash
@GtkBuilderAid function_name(empty_builder) begin

end
empty_builder("resources/nothing.ui", test_app; wpipe=nullio)

# Check that a poorly made builder doesn't crash
@GtkBuilderAid function_name(poor_builder) begin

function click_ok()
end

function close_window()
end

end
poor_builder("resources/nothing.ui", test_app; wpipe=nullio)

type InternalData
  counter::Int
end

@GtkBuilderAid function_name(long_builder) begin

@guarded function click_ok(
    widget,
    data)
  data.counter += 1
  return nothing
end

@guarded function close_window(
    widget,
    window)
  destroy(window)
  return nothing
end

end

@test_throws MethodError long_builder()

type ActivateData
  app
  data
end

@guarded function activate_long_builder_app(app, data)
  app = GObject(app)
  built = long_builder("resources/nothing.ui", data)

  w = Gtk.GAccessor.object(built, "main_window")
  push!(app, w)
  showall(w)

  # Try triggering the clicked signal
  ccall(
    (:g_signal_emit_by_name, Gtk.libgtk),
    Void,
    (Ptr{GObject}, Ptr{Int8}),
    Gtk.GAccessor.object(built, "ok_button"), "clicked")

  ccall(
    (:g_signal_emit_by_name, Gtk.libgtk),
    Void,
    (Ptr{GObject}, Ptr{Int8}),
    Gtk.GAccessor.object(built, "cancel_button"), "clicked")

  quit(app)
  return nothing
end

data = InternalData(0)
@test data.counter == 0
signal_connect(activate_long_builder_app, test_app, "activate", Void, (), false, data)
run(test_app)
@test data.counter == 1

# Show the expanded macro
# Mostly check that this succeeds
@GtkBuilderAid function_name(builder) userdata(test_app) "resources/nothing.ui" begin

function click_ok(
    widget,
    user_info)
  println("OK clicked!")
  return nothing
end

function close_window(
    widget,
    window)
  destroy(window)
  return nothing
end

end

builder()
# Also check that the unbound form works
builder("$(Pkg.dir("GtkBuilderAid"))/test/resources/nothing.ui")

@GtkBuilderAid function_name(tuple_builder) userdata_tuple(test_app::GtkApplication) "resources/nothing.ui" begin

function click_ok(
    widget,
    user_info)
  println("OK clicked!")
  return nothing
end

function close_window(
    widget,
    window)
  destroy(window)
  return nothing
end

end

# check again but with an explicit name for the builder function
@GtkBuilderAid function_name(build_nothing) userdata(test_app::GtkApplication) "resources/nothing.ui" begin

function close_window(
    widget,
    window)
  destroy(window)
  return nothing
end

function click_ok(
    widget,
    user_info)
  println("OK clicked!")
  return nothing
end

end
build_nothing()

# Test non-string file arguments
@GtkBuilderAid function_name(base_method_builder) begin

function close_window(
    widget,
    window)
  window = GObject(window)
  destroy(window)
  return nothing
end

function click_ok(
    widget,
    user_info)
  println("OK clicked!")
  return nothing
end

end

@test_throws MethodError base_method_builder()
# Should succeed
base_method_builder("resources/nothing.ui")

# Test non-existent files
test_macro_throws(ErrorException, quote
@GtkBuilderAid "resources/nonexistentfile.ui" begin
function close_window(
    widget,
    window)
  destroy(window)
  return nothing
end

end
end)

expanding_builder = @GtkBuilderAid begin

@guarded function close_window(
    widget,
    window)
  destroy(window)
  return nothing::Void
end

@guarded function click_ok(
    widget,
    user_info)
  println("OK clicked!")
  return nothing::Void
end

end
expanding_builder("resources/nothing.ui")

shorthand_builder = @GtkBuilderAid begin

close_window(
    widget,
    window) = begin
  destroy(window)
  return nothing::Void
end

click_ok(
    widget,
    user_info) = begin
  println("OK clicked!")
  return nothing::Void
end

end
shorthand_builder("resources/nothing.ui")

# Test using an external builder
external_builder = @GtkBuilderAid begin

close_window(
    widget,
    window) = begin
  destroy(window)
  return nothing::Void
end

click_ok(
    widget,
    user_info) = begin
  println("OK clicked!")
  return nothing::Void
end

end
builder_obj = @GtkBuilder(filename = "resources/nothing.ui")
external_builder(builder_obj)

# create a drag-drop interface
# TODO complete this and use it as ex3
immutable GtkTargetEntry
  target::Ptr{UInt8}
  flags::Cuint
  info::Cuint
end

@GtkBuilderAid function_name(drag_drop_builder) begin

@guarded function close_window(widget, window)
  destroy(window)
  nothing
end

@guarded function text_view_realize_cb(
    widget,
    userinfo)
  # enable drag-drop stuff
  targets = GtkTargetEntry[]
  ccall(
      (:gtk_drag_dest_set, Gtk.libgtk),
      Void,
      (Ptr{GObject}, Cint, Ptr{GtkTargetEntry}, Cint, Cint),
      widget,
      Gtk.GConstants.GtkDestDefaults.DROP,
      targets,
      length(targets),
      Gtk.GConstants.GdkDragAction.COPY)
  nothing
end

@guarded function text_view_drag_data_received_cb(
    widget,
    context,
    x,
    y,
    data,
    info,
    time,
    userinfo)
  println("Drag Received")
  nothing
end

@guarded Cint(0) function text_view_drag_drop_cb(
    widget,
    context,
    x,
    y,
    time,
    userinfo)
  # drag-drop check area
  println("Drag Drop")
  Cint(1)
end

end

built = drag_drop_builder("resources/multiparameter.ui")

# Based off of GTK's Custom Drawing source code

# Test using canvas features of DrawingArea
type CanvasData
  surface::CairoSurface
end

# Clear the surface to white
function clear_surface(surf::CairoSurface)
  ctx = CairoContext(surf)
  set_source_rgb(ctx, 1, 1, 1)
  paint(ctx)
  destroy(ctx)
end

function draw_brush(
    canvas,
    x,
    y,
    surf)

  px = Int(round(x - 3))
  py = Int(round(y - 3))
  ctx = CairoContext(surf)
  rectangle(ctx, px, py, 6, 6)
  fill(ctx)
  destroy(ctx)

  # Force the widget to update that part of the drawing
  ccall(
      (:gtk_widget_queue_draw_area, Gtk.libgtk),
      Void,
      (Ptr{GObject}, Cint, Cint, Cint, Cint),
      canvas,
      px,
      py,
      6,
      6)
end

@GtkBuilderAid function_name(canvas_builder) begin

@guarded Cint(0) function configure_event_cb(
    canvas,
    configure_event,
    userdata)
  destroy(userdata.surface)
  w = width(canvas)
  h = height(canvas)
  userdata.surface = CairoSurface(
    ccall((:gdk_window_create_similar_surface, Gtk.libgtk),
      Ptr{Void},
      (Ptr{Void}, Gtk.GEnum, Cint, Cint),
      Gtk.gdk_window(canvas),
      Cairo.CONTENT_COLOR_ALPHA,
      w,
      h),
    w,
    h)
  clear_surface(userdata.surface)
  return Cint(1)
end

@guarded Cint(1) function draw_cb(
    canvas,
    ctx_ptr,
    userdata)

  ctx = CairoContext(ctx_ptr)

  set_source_surface(ctx, userdata.surface, 0, 0)
  paint(ctx)
  return Cint(0)
end

@guarded Cint(0) function button_press_event_cb(
    canvas,
    event_ptr,
    userdata)

  event = Gtk.GdkEvent(event_ptr)

  # Make sure the surface is valid
  if userdata.surface.ptr == C_NULL
    return Cint(0)
  end

  if event.button == 1
    draw_brush(canvas, event.x, event.y, userdata.surface)
  elseif event.button == 3
    clear_surface(userdata.surface)
    ccall(
        (:gtk_widget_queue_draw, Gtk.libgtk),
        Void,
        (Ptr{GObject}, ),
        canvas)
  end

  return Cint(1)
end

@guarded Cint(0) function motion_notify_event_cb(
    canvas,
    event_ptr,
    userdata)

  event = Gtk.GdkEvent(event_ptr)

  # Make sure the surface is valid
  if userdata.surface.ptr == C_NULL
    return Cint(0)
  end

  if (event.state & Gtk.GdkModifierType.BUTTON1) != 0
    draw_brush(canvas, event.x, event.y, userdata.surface)
  end

  return Cint(1)

end

# Destroy the cairo surface when the window closes
@guarded function close_window(
    window,
    userdata)
  destroy(userdata.surface)
  return nothing
end

end

# Start with an under-defined cairo surface
canvas_builder("resources/draw.ui", CanvasData(CairoSurface(C_NULL, -1, -1)))

@test_throws ErrorException builder("resources/nonexistant.ui")

# run(test_app)


text:
using NewtonsMethod
using Test

@testset "NewtonsMethod.jl" begin
    # test 1: make sure it is working
        g(x)=(x-1)^3
        g_prime(x)=3*(x-1)^2
        iv = 1.1
        sol = newtonroot(g, g_prime, iv)
        @test sol≈ 1.0000001545213348
    # test 2: make sure that the second dispach works well
        solu = newtonroot(g, iv)
    @test solu ≈ 1.0000001545213348

    # test 3: make sure tolerance works
    solu = newtonroot(g, g_prime, iv, tolerance = 1)
    solut= newtonroot(g, iv, tolerance =  1)
    @test !( solu  ≈ 1.0000001545213348)
    @test !( solut ≈ 1.00000013565)

    # test 4: make sure maxiter works
    soluti  = newtonroot(g, g_prime, iv,maxiter= 5 )
    solutio = newtonroot(g, iv, maxiter= 5 )
    @test soluti  == nothing
    @test solutio == nothing
    #test 5: Non-Convergence
    h(x)=x^2+2
    h_prime(x)=2x
    @test newtonroot(h, h_prime, iv)== nothing
    @test newtonroot(h,  iv)== nothing
    #test BiFloat
    iv=BigFloat(3)
    @test newtonroot(g,g_prime, iv)≈ 1.0000001808754536763325638480
    @test newtonroot(g, iv)≈ 1.0000001808754536763325638480
    # test Int
    iv=3
    @test newtonroot(g,g_prime, iv)≈ 1.0000001808754536763325638480
    @test newtonroot(g, iv)≈ 1.0000001808754536763325638480
    # test complicated function
    k(x)=exp(x)-3+sin(x)
    k_prime(x)=exp(x)+cos(x)
    @test newtonroot(k,k_prime, iv)≈ 0.8194425797965913
    @test newtonroot(k, iv)≈ 0.8194425797965913
end


text:
# script to install all deps
import Pkg; Pkg.add("JSON")
import Pkg; Pkg.add("SHA")
import Pkg; Pkg.add("Dates")
import Bisquit; Pkg.add("Bisquit")
import Genie, Genie.Router, Genie.Renderer.Json, Genie.Requests; Pkg.add("Genie")
import HTTP; Pkg.add("HTTP")
import UUIDs; Pkg.add("UUIDs")

max_stars_repo_name: timo-cmd/Bisquit


text:
@testset "logistic" begin
    N, d = 20, 2
    k = 2.0 * SqExponentialKernel() ∘ ScaleTransform(10.0)
    X, f = generate_f(N, d, k)
    y = f .> 0
    floattypes = [Float64]
    tests_likelihood(
        LogisticLikelihood(),
        LogisticLikelihood,
        Dict(
            "VGP" => Dict("AVI" => true, "QVI" => true, "MCVI" => false),
            "SVGP" => Dict("AVI" => true, "QVI" => true, "MCVI" => false),
            "OSVGP" => Dict("AVI" => true, "QVI" => false, "MCVI" => false),
            "MCGP" => Dict("Gibbs" => true, "HMC" => false),
        ),
        floattypes,
        "Classification",
        1,
        X,
        f,
        y,
        k,
    )
end


text:
##################################################################################
# This file is part of ModelBaseEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

using Lazy: @forward

# export Parameters, ParamAlias, ParamLink, peval, @alias, @link
export Parameters, ModelParam, peval, @alias, @link

abstract type ParamAbstract end

"""
    struct Parameters <: AbstractDict{Symbol, Any}

Container for model parameters. It functions as a `Dict` where the keys are the
parameter names. Simple parameter values are stored directly. Special parameters
depend on other parameters are are wrapped in the appropriate data structures to
keep track of such dependencies. There are two types of special parameters -
aliases and links.

Individual parameters can be accessed in two different ways - dot and bracket
notation.

Read access by dot notation calls [`peval`](@ref) while bracket notation
doesn't. This makes no difference for simple parameters. For special parameters,
access by bracket notation returns its internal structure, while access by dot
notation returns its current value depending on other parameters.

Write access is the same in both dot and bracket notation. The new parameter
value is assigned directly in the case of simple parameter. To create an alias
parameter, use the [`@alias`](@ref) macro. To create a link parameter use the
[`@link`](@ref) macro.

See also: [`ModelParam`](@ref), [`peval`](@ref), [`@alias`](@ref),
[`@link`](@ref), [`update_links!`](@ref).
"""
struct Parameters <: AbstractDict{Symbol,ParamAbstract}
    mod::Ref{Module}
    contents::Dict{Symbol,ParamAbstract}
end

"""
    Parameters([mod::Module])

When creating an instance of `Parameters`, optionally one can specify the module
in which parameter expressions will be evaluated. This only matters if there are
any link parameters that depend on custom functions or global
variables/constants. In this case, the `mod` argument should be the module in
which these definitions exist.
"""
Parameters(mod::Module=@__MODULE__) = Parameters(Ref(mod), Dict{Symbol,ParamAbstract}())

"""
    params = @parameters

When called without any arguments, return an empty [`Parameters`](@ref)
container, with its evaluation module set to the module in which the macro is
being called.
"""
macro parameters()
    return :( Parameters($__module__) )
end

# To deepcopy() Parameters, we make a new Ref to the same module and a deepcopy of contents.
Base.deepcopy_internal(p::Parameters, stackdict::IdDict) = Parameters(Ref(p.mod[]), Base.deepcopy_internal(p.contents, stackdict))

# The following functionality is forwarded to the contents
# iteration
@forward Parameters.contents Base.keys, Base.values, Base.pairs
@forward Parameters.contents Base.iterate, Base.length
# bracket notation read access
@forward Parameters.contents Base.getindex
# dict access
@forward Parameters.contents Base.get, Base.get!

"""
    mutable struct ModelParam

Contains a model parameter. For a simple parameter it simply stores its value.
For a link or an alias, it stores the link information and also caches the
current value for speed.
"""
mutable struct ModelParam <: ParamAbstract
    depends::Set{Symbol}  # stores the names of parameter that depend on this one
    link::Union{Nothing,Symbol,Expr}
    value
end
ModelParam() = ModelParam(Set{Symbol}(), nothing, nothing)
ModelParam(value) = ModelParam(Set{Symbol}(), nothing, value)
ModelParam(value::Union{Symbol,Expr}) = ModelParam(Set{Symbol}(), value, nothing)

"""
    @alias name

    Create a parameter alias. Use `@alias` in the [`@parameters`](@ref) section of your
model definition.
```
@parameters model begin
    a = 5
    b = @alias a
end
```
"""
macro alias(arg)
    return arg isa Symbol ? ModelParam(arg) : :( throw(ArgumentError("`@alias` requires a symbol. Use `@link` with expressions.")) )
end

"""
    @link expr

Create a parameter link. Use `@link` in the [`@parameters`](@ref) section of
your model definition.

If your parameter depends on other parameters, then you use `@link` to declare
that. The expression can be any valid Julia code.
```
@parameters model begin
    a = 5
    b = @link a + 1
end
```

When a parameter the link depends on is assigned a new value, the link that
depends on it gets updated automatically.

!!! note "Important note"
    There are two cases in which the value of a link does not get updated automatically.
    If the parameter it depends on is mutable, e.g. a `Vector`, it is possible for it to get
    updated in place. The other case is when the link contains global variable or custom function.

    In such case, it is necessary to call [`update_links!`](@ref).
"""
macro link(arg)
    return arg isa Union{Symbol,Expr} ? ModelParam(arg) : :( throw(ArgumentError("`@link` requires an expression.")) )
end

Base.show(io::IO, p::ModelParam) = begin
    p.link === nothing ? print(io, p.value) :
    p.link isa Symbol ? print(io, "@alias ", p.link) :
                        print(io, "@link ", p.link)
end

Base.:(==)(l::ModelParam, r::ModelParam) = l.link === nothing && l.value == r.value || l.link == r.link

###############################
# setindex!

_value(val) = val isa ModelParam ? val.value : val
_link(val) = val isa ModelParam ? val.link : nothing

function _rmlink(params, p, key)
    # we have to remove `key` from the `.depends` of all parameters we depend on
    MacroTools.postwalk(p.link) do e
        if e isa Symbol
            ep = get(params.contents, e, nothing)
            ep !== nothing && delete!(ep.depends, key)
        end
        return e
    end
    return
end

function _check_circular(params, val, key)
    # we have to check for circular dependencies
    deps = Set{Symbol}()
    while true
        MacroTools.postwalk(_link(val)) do e
            if e isa Symbol
                if e == key
                    throw(ArgumentError("Circular dependency of $(key) and $(e) in redefinition of $(key)."))
                end
                if e in keys(params.contents)
                    push!(deps, e)
                end
            end
            return e
        end
        if isempty(deps)
            return
        end
        k = pop!(deps)
        val = params[k]
    end
end

function _addlink(params, val, key)
    if !isa(val, ModelParam) || val.link === nothing
        return
    end
    # we have add `key` to the `.depends` of all parameters we depend on
    MacroTools.postwalk(val.link) do e
        if e isa Symbol
            ep = get(params.contents, e, nothing)
            ep !== nothing && push!(ep.depends, key)
        end
        return e
    end
    return
end

"""
    peval(params, what)

Evaluate the given expression in the context of the given parameters.

If `what` is a `ModelParam`, its current value is returned. If there's a chance
it might be out of date, call [`update_links!`](@ref).

If `what` is a Symbol or an Expr, all mentions of parameter names are
substituted by their values and the the expression is evaluated.

If `what is any other value, it is returned unchanged.`

See also: [`Parameters`](@ref), [`@alias`](@ref), [`@link`](@ref),
[`ModelParam`](@ref), [`update_links!`](@ref).

"""
peval(params, val) = val
peval(params, par::ModelParam) = par.value
peval(params, sym::Symbol) = sym ∈ keys(params) ? peval(params, params[sym]) : sym
function peval(params, expr::Expr)
    ret = Expr(expr.head)
    ret.args = [peval(params, a) for a in expr.args]
    params.mod[].eval(ret)
end

function _update_values(params, p, key)
    # update my own value
    if p.link !== nothing
        p.value = peval(params, p.link)
    end
    # update values that depend on me
    deps = copy(p.depends)
    while !isempty(deps)
        pk = params.contents[pop!(deps)]
        pk.value = peval(params, pk.link)
        if !isempty(pk.depends)
            push!(deps, pk.depends...)
        end
    end
end

"""
    iterate(params::Parameters)

Iterates the given Parameters collection in the order of dependency.
Specifically, each parameter comes up only after all parameters it depends on
have already been visited. The order within that is alphabetical.
"""
function Base.iterate(params::Parameters, done=Set{Symbol}())
    if length(done) == length(params.contents)
        return nothing
    end
    for k in sort(collect(keys(params.contents)))
        if k in done
            continue
        end
        v = params.contents[k]
        if v.link === nothing
            return k => v, push!(done, k)
        end
        ready = true
        MacroTools.postwalk(v.link) do e
            if e ∈ keys(params) && e ∉ done
                ready = false
            end
            e
        end
        if ready
            return k => v, push!(done, k)
        end
    end
end

export update_links!

"""
    update_links!(model)
    update_links!(params)

Recompute the current values of all parameters.

Typically when a new value of a parameter is assigned, all parameter links and
aliases that depend on it are updated recursively. If a parameter mutable, e.g.
a Vector or another collection, its value can be updated in place without
re-assigning it, thus the automatic updated does not happen. In this case, it is
necessary to call `update_links!`.
"""
update_links!(m::AbstractModel) = update_links!(parameters(m))
function update_links!(params::Parameters)
    for (k, v) in params
        if v.link !== nothing
            v.value = peval(params, v.link)
        end
    end
end

function Base.setindex!(params::Parameters, val, key)
    if key in fieldnames(typeof(params))
        throw(ArgumentError("Invalid parameter name: $key."))
    end
    # If param[key] is not a link, then key doesn't appear in anyone's depends
    # Invariant: my.depends always contains the names of parameters that depend on me.
    _check_circular(params, val, key)
    p = get!(params.contents, key, ModelParam())
    _rmlink(params, p, key)
    _addlink(params, val, key)
    p.link = _link(val)
    p.value = _value(val)
    _update_values(params, p, key)
    params
end

Base.propertynames(params::Parameters) = tuple(keys(params)...)

function Base.setproperty!(params::Parameters, key::Symbol, val)
    if key ∈ fieldnames(typeof(params))
        return setfield!(params, key, val)
    else
        return setindex!(params, val, key)
    end
end

function Base.getproperty(params::Parameters, key::Symbol)
    if key ∈ fieldnames(typeof(params))
        return getfield(params, key)
    end
    par = get(params, key, nothing)
    if par === nothing
        throw(ArgumentError("Unknown parameter $key."))
    else
        return peval(params, par)
    end
end


text:
# Don't change the value below, it's used in an `include` test
55


text:
function Base.copy(x::Camera)
    Camera(ntuple(7) do i
        getfield(x, i)
    end...)
end

function Base.:(==)(a::Camera, b::Camera)
    to_value(a.view) == to_value(b.view) &&
    to_value(a.projection) == to_value(b.projection) &&
    to_value(a.resolution) == to_value(b.resolution)
end

function Base.show(io::IO, camera::Camera)
    println(io, "Camera:")
    println(io, "  $(length(camera.steering_nodes)) steering observables connected")
    println(io, "  pixel_space: ", camera.pixel_space[])
    println(io, "  view: ", camera.view[])
    println(io, "  projection: ", camera.projection[])
    println(io, "  projectionview: ", camera.projectionview[])
    println(io, "  resolution: ", camera.resolution[])
    println(io, "  eyeposition: ", camera.eyeposition[])
end

function disconnect!(c::Camera)
    for obsfunc in c.steering_nodes
        off(obsfunc)
    end
    empty!(c.steering_nodes)
    return
end

function disconnect!(c::EmptyCamera)
    return
end

function disconnect!(observables::Vector)
    for obs in observables
        disconnect!(obs)
    end
    empty!(observables)
    return
end

struct CameraLift{F, Args}
    f::F
    args::Args
end

function (cl::CameraLift{F, Args})(val) where {F, Args}
    cl.f(map(to_value, cl.args)...)
end

"""
    on(f, c::Camera, observables::Observable...)

When mapping over observables for the camera, we store them in the `steering_node` vector,
to make it easier to disconnect the camera steering signals later!
"""
function Observables.on(f::Function, camera::Camera, observables::AbstractObservable...; priority=Int8(0))
    # PriorityObservables don't implement on_any, because that would replace
    # the method in Observables. CameraLift acts as a workaround for now.
    cl = CameraLift(f, observables)
    for n in observables
        obs = if n isa PriorityObservable
            on(cl, n, priority=priority)
        else
            on(cl, n)
        end
        push!(camera.steering_nodes, obs)
    end
    return f
end

function Camera(px_area)
    pixel_space = lift(px_area) do window_size
        nearclip = -10_000f0
        farclip = 10_000f0
        w, h = Float32.(widths(window_size))
        return orthographicprojection(0f0, w, 0f0, h, nearclip, farclip)
    end
    view = Observable(Mat4f(I))
    proj = Observable(Mat4f(I))
    proj_view = map(*, proj, view)
    Camera(
        pixel_space,
        view,
        proj,
        proj_view,
        lift(a-> Vec2f(widths(a)), px_area),
        Observable(Vec3f(1)),
        ObserverFunction[]
    )
end

function set_proj_view!(camera::Camera, projection, view)
    # hack, to not double update projectionview
    # TODO, this makes code doing on(view), not work correctly...
    # But nobody should do that, right?
    camera.view.val = view
    camera.projection[] = projection
end

function is_mouseinside(scene, target)
    scene === target && return false
    Vec(scene.events.mouseposition[]) in pixelarea(scene)[] || return false
    for child in r.children
        is_mouseinside(child, target) && return true
    end
    return false
end

"""
    is_mouseinside(scene)

Returns true if the current mouseposition is inside the given scene.
"""
function is_mouseinside(scene)
    return Vec(scene.events.mouseposition[]) in pixelarea(scene)[]
    # Check that mouse is not inside any other screen
    # for child in scene.children
    #     is_mouseinside(child) && return false
    # end
end


text:
#=
Model: MNA_1 (578 variables, 9 inputs)
=#
using Reachability, MAT, SparseArrays

mna1(o::Pair{Symbol, <:Any}...) = mna1(Options(Dict{Symbol, Any}(o)))

function mna1(input_options::Options)
    # =====================
    # Problem specification
    # =====================
    # dynamics
    file = matopen(@relpath "mna1.mat")
    A = sparse(read(file, "A"))
    n = size(A, 1)
    b = sparsevec(570:578, [fill(-0.1, 5); fill(-0.2, 4)], n)

    # initial set
    X0 = Hyperrectangle([fill(0.00125, 2); zeros(576)],
                        [fill(0.00025, 2); zeros(576)])

    # instantiate continuous LTI system
    S = IVP(CACS(A, b, nothing), X0)

    # property: x1 < 0.5
    property = SafeStatesProperty(HalfSpace(sparsevec([1], [1.0], n), 0.5))

    # =======================
    # Problem default options
    # =======================
    partition = [(2*i-1:2*i) for i in 1:289] # 2D blocks

    if input_options[:mode] == "reach"
        problem_options = Options(:vars => [1],
                                  :partition => partition,
                                  :plot_vars => [0, 1])
    elseif input_options[:mode] == "check"
        problem_options = Options(:vars => [1],
                                  :partition => partition,
                                  :property => property)
    end

    return (S, merge(problem_options, input_options))
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
using DataFrames, CSV, Missings

dirpath = "d:/data/fannie_mae/"
dirpaths = readdir(dirpath)
filepath = joinpath(dirpath, "Performance_2000Q1.txt")

const types = [
    String,                     Union{String, Missing},     Union{String, Missing},     Union{Float64, Missing},    Union{Float64, Missing},
    Union{Float64, Missing},    Union{Float64, Missing},    Union{Float64, Missing},    Union{String, Missing},     Union{String, Missing},
    Union{String, Missing},     Union{String, Missing},     Union{String, Missing},     Union{String, Missing},     Union{String, Missing},
    Union{String, Missing},     Union{String, Missing},     Union{Float64, Missing},    Union{Float64, Missing},    Union{Float64, Missing},
    Union{Float64, Missing},    Union{Float64, Missing},    Union{Float64, Missing},    Union{Float64, Missing},    Union{Float64, Missing},
    Union{Float64, Missing},    Union{Float64, Missing},    Union{Float64, Missing},    Union{String, Missing},     Union{Float64, Missing},
    Union{String, Missing}]

@time perf = CSV.read(filepath, delim='|', header = false, types = types, weakrefstrings=false) #45;

@time Feather.write("julia_feather.feather", perf)
using JLD2, FileIO, DataFrames
df = DataFrame(a=[1])

# @time @save "df.jld2" perf # too slow

const ufm = Union{Float64, Missing}
# parse Fannie Mae performance data
function parsefm(T, s)
    if T == ufm
        if s  == ""
            return missing
        end
        return parse(Float64, s)
    end
    String(s)
end

strings_to_parse = ["100007365142|11/01/2001||8|73753.13|22|338|337|01/2030|00000|0|N|||||||||||||||||||",
"100007365142|12/01/2001|JPMORGAN CHASE BANK, NA|8|73694.5|23|337|336|01/2030|00000|0|N|||||||||||||||||||"]

df = DataFrame(a = strings_to_parse)

CSV.write("mwe.txt", df, header = false)

df1 = CSV.read("mwe.txt", header=false)
# 2×1 DataFrames.DataFrame
# │ Row │ Column1                                                                             │
# ├─────┼─────────────────────────────────────────────────────────────────────────────────────┤
# │ 1   │ 100007365142|11/01/2001||8|73753.13|22|338|337|01/2030|00000|0|N||||||||||||||||||| │
# │ 2   │ 100007365142|12/01/2001|JPMORGAN CHASE BANK

df2 = CSV.read("mwe.txt", header=false, delim = '|')

split.(strings_to_parse, '|')

function convert2tuple(r)
    ([parsefm(T, s) for (T,s) in zip(types, split(r, '|'))]...)
end

@time hehe = convert2tuple.(perfraw[:Column1]);

for i =1:size(perfraw,1)
    println(i)
    convert2tuple(perfraw[:Column1][i])
end

@time abc(perfraw)

    pv = ["LOAN_ID", "Monthly.Rpt.Prd", "Servicer.Name", "LAST_RT", "LAST_UPB", "Loan.Age", "Months.To.Legal.Mat" , "Adj.Month.To.Mat",
"Maturity.Date", "MSA", "Delq.Status", "MOD_FLAG", "Zero.Bal.Code", "ZB_DTE", "LPI_DTE", "FCC_DTE","DISP_DT", "FCC_COST", "PP_COST", "AR_COST",
"IE_COST", "TAX_COST", "NS_PROCS", "CE_PROCS", "RMW_PROCS", "O_PROCS", "NON_INT_UPB", "PRIN_FORG_UPB_FHFA", "REPCH_FLAG", "PRIN_FORG_UPB_OTH",
"TRANSFER_FLAG"]

@time perf = CSV.read(filepath, delim='|', header = false, types = types) #45;

using RCall
using Feather
function rfread(path)
    R"""
    feather::write_feather(data.table::fread($path), 'tmp.feather')
    gc()
    """
    Feather.read("tmp.feather")
end

@time perf = rfread(filepath) # 45

@time Feather.read("tmp.feather")

abctp = ([abc(T,s) for (T,s) in zip(types, split(perf[1,1],'|'))]...)

DataFrame(abctp)

function()
    df = abc(perf[1,1])
for r in eachrow(perf)
    print(([abc(T,s) for (T,s) in zip(types, split(r[1],'|'))]...))
    break
end

using DataFrames
const N = 100_000_000; const K = 100
srand(1)
@time df = DataFrame(idstr = rand("id".*dec.(1:N÷K,10), N)
        , id = rand(1:K, N)
        , val = rand(1:5,N))
using CSV
@time CSV.write("df.csv", df);

for i in 1:ncol(perf)
    println(pv[i],":", sum(ismissing.(perf[i])))
end

@time perf = CSV.read(filepath, header=false)

@time perf1 = readtable(filepath, separator = '|', header=false) # 221 # 261

using Feather
@time af  = Feather.read("d:/p.feather")

using TextParse
@time iperf = TextParse.csvread(filepath, '|', type_detect_rows = 11857)

using uCSV

uCSV.read(filepath, delim = '|', typedetectrows = 11857)

using FileIO, CSVFiles, DataFrames, FeatherFiles
df = DataFrame(load(File(format"CSV", filepath), '|'))

using Feather

using RCall
R"""
library(feather)
system.time(read_feather("d:/p.feather"))
"""

df = Feather.read("d:/p.feather")

R"""
memory.limit(2^31-1)
library(data.table)
library(feather)
list(
system.time(a <- fread($filepath)),
system.time(write_feather(a, "d:/p.feather")))
"""


text:
"""
    get_times(x::AbstractVector{<:Real})

Get the times associated with the inputs.
"""
get_times(x::AbstractVector{<:Real}) = x

"""
    inputs_to_time_form(x::AbstractVector)

Get the time-form representation of a vector of inputs.
Outputs an `AbstractVector` whose length is `length(get_times(x))`,
and whose elements comprise all of those associated with each time point.

For single-output time series problems, `x` is usually returned as-is, however, for
multi-output and spatio-temporal problems there are usually multiple inputs associated with
each point in time.
"""
inputs_to_time_form(x::AbstractVector{<:Real}) = Fill(nothing, length(x))

"""
    merge_inputs(x1::AbstractVector, x2::AbstractVector)

Merge two collections of inputs `x1` and `x2` if possible.
A notable example of inputs which cannot be merged are two rectilinear grids whose
spatial components differ.
"""
merge_inputs(x1::AbstractVector{<:Real}, x2::AbstractVector{<:Real}) = vcat(x1, x2)

"""
    sort_in_time(x::AbstractVector)

Sort the input `x` temporally.
This is different from sorting the input in general, but agrees for real-valued inputs.
Returns both the sorting indices and the sorted inputs.
"""
function sort_in_time(x::AbstractVector{<:Real})
    idx = sortperm(x)
    return idx, x[idx]
end

"""
    observations_to_time_form(x::AbstractVector, y::AbstractVector{<:Real})

Get the time-form representation of a vector of observations `y` of length
`length(time_points(x))`, which is useful for state-space models.
The precise arrangements of the elements of `y` in the output depends on `x`.

For example, single-output time-series problems typically just return `y`, but multi-output
and spatio-temporal problems typically have multiple elements of `y` associated with a
single element of `x`.
"""
function observations_to_time_form(
    x::AbstractVector{<:Real},
    y::AbstractVector{<:Union{Missing, Real}},
)
    return y
end

"""
    noise_var_to_time_form(x::AbstractVector, S::AbstractMatrix{<:Real})

Get the time-form representation of an `AbstractMatrix` of observation noise variances `S`,
based on the `AbstractVector` of inputs `x`.

For example, single-output time problems will typically require `S` to be `Diagonal`, and
will just return its diagonal.
Multi-output and spatio-temporal problems can potentially have a block-diagonal `S`, as
it's typically only necessary to require that the observation variance at each point in time
is independent of other times.
"""
noise_var_to_time_form(x::AbstractVector{<:Real}, S::Diagonal{<:Real}) = diag(S)

"""
    destructure(x::AbstractVector, ys::AbstractVector)

Construct a flattened representation of `ys` in accordance with the structure in `x`.
"""
destructure(x::AbstractVector{<:Real}, ys::AbstractVector) = ys


text:
"""
$(TYPEDSIGNATURES)

Return [`AnyElement`](@ref) wraping a `tag` symbol and `PnmlDict` holding
a well-formed XML node.

See [`ToolInfo`](@ref) for one intended use-case and [`unclaimed_label`](@ref) for when
a pnml label is expected but does not have a parser or the tag appears in an
unexpected place.
"""
function anyelement end
anyelement(node; kw...) =  anyelement(node, PnmlCore(); kw...)
function anyelement(node, pntd; kw...)::AnyElement
    d = _harvest_any!(node, pntd, _harvest_any!; kw...)
    return AnyElement(Symbol(nodename(node)) => d, node)
end

"""
$(TYPEDSIGNATURES)

Return `tag` => `PnmlDict` holding a pnml label and its children.

The main use-case is to wrap in a [`PnmlLabel`](@ref), [`Structure`](@ref),
[`Term`](@ref) or other specialized label. These wrappers add type to the
nested dictonary holding the contents of the label.

Differs from `AnyElement` in that any id attribute of `node`
will be registered in [`IDRegistry`](@ref) as a unique identifier.
"""
function unclaimed_label end
unclaimed_label(node; kw...) = unclaimed_label(node, PnmlCore(); kw...)
function unclaimed_label(node, pntd; kw...)::Pair{Symbol,PnmlDict}
    @assert haskey(kw, :reg)
    # ID attributes can appear in various places. Each is unique and added to the registry.
    #EzXML.haskey(node, "id") && register_id!(kw[:reg], node["id"])
    # Children may be claimed.
    return Symbol(nodename(node)) => _harvest_any!(node, pntd, _harvest_any!; kw...)
end

"""
$(TYPEDSIGNATURES)

Return `PnmlDict` holding a well-formed XML `node`.

If element `node` has any children, each is placed in the dictonary with the
child's tag name symbol as the key, repeated tags produce a vector as the value.
Any XML attributes found are added as as key,value pairs.

Descend the well-formed XML using `parser` on child nodes.

Note the assumption that "children" and "content" are mutually exclusive.
Content is always a leaf element. However XML attributes can be anywhere in
the hiearchy. And neither children nor content nor attribute may be present.
"""
function _harvest_any!(node::XMLNode, pntd::PNTD, parser; kw...) where {PNTD<:PnmlType}
    @assert haskey(kw, :reg)
    # Extract XML attributes. Register IDs as symbols.
    dict = PnmlDict()
    for a in eachattribute(node)
        dict[Symbol(a.name)] = a.name == "id" ? register_id!(kw[:reg], a.content) : a.content
    end

    # Extract children or content
    children = elements(node)
    if !isempty(children)
        merge!(dict, _anyelement_content(children, pntd, parser; kw...))
    elseif !isempty(nodecontent(node))
        # <tag> </tag> will have nodecontent, though the whitespace is discarded.
        dict[:content] = strip(nodecontent(node))
    else
        # <tag/> and <tag></tag> will not have any nodecontent.
        #TODO Force dict[:content] = "" ?
        #TODO <tag/> serves as a flag. (is key present?)
    end
    #@show dict
    return dict
end

"""
$(TYPEDSIGNATURES)

Apply `parser` to each node in `nodes`.
Return PnmlDict with values that are vectors when there
are multiple instances of a tag in `nodes` and scalar otherwise.
"""
function _anyelement_content(nodes::Vector{XMLNode}, pntd::PNTD, parser; kw...) where {PNTD<:PnmlType}
    namevec = [nodename(node) => node for node in nodes] # Not yet turned into Symbols.
    tagnames = unique(map(first, namevec))
    dict = PnmlDict()
    foreach(tagnames) do tagname
        tags = filter(x->x.first===tagname, namevec)
        dict[Symbol(tagname)] = if length(tags) > 1 # Now its a symbol.
            parser.(map(x->x.second, tags), Ref(pntd), Ref(parser); kw...) # vector
        else
            parser(tags[1].second, pntd, parser; kw...) # scalar
        end
    end
    return dict
end


text:
clibrary(:colorcet)
plot(x, [sin.(x), cos.(x)], m = ([5 5], [:d :rect], stroke(1)),
    c = [:fire :plasma], markerz =  x, linez = x,
    colorbar = false, ls = [:dash :dashdot], leg =:false,
    ylim = (-1.5, 1.5), size = (600, 300))
plot!(x, fill(1.2,20), c = :viridis, markerz =  x, linez = x,
    m = ([10, 5, 4], [:d, :rect, :utriangle], stroke(0)))
plot!(x, fill(-1.2,20), c = :kdc, markerz = x, linez = x, ls = :dashdotdot,
    m = ([10,8,6,4],[:d,:rect,:utriangle,:circle],stroke(1)))
savefig("figmixmarkers.pdf")
# second plot
tronLegacy = ["#FF410DFF", "#6EE2FFFF", "#F7C530FF", "#95CC5EFF",
    "#D0DFE6FF", "#F79D1EFF", "#748AA6FF"]
plot(x, [sin.(x), cos.(x)], m = ([5 5], [:d :rect], stroke(1)),
    c = [tronLegacy[1] tronLegacy[2]], ls = [:dash :dashdot],
    colorbar = false,  ylim = (-1.5, 1.5), leg =:false, bg=:black,
    size = (600, 400))
plot!(x, fill(1.2,20), m = ([10, 5, 4], [:d, :rect, :utriangle], stroke(0)),
    c = tronLegacy[3])
plot!(x, fill(-1.2,20), c = tronLegacy[5],ls = :dashdotdot,
    m = ([10, 8, 6, 4], [:d, :rect, :utriangle, :circle], stroke(1)))
savefig("figmixmarkerstron.pdf")

max_stars_repo_name: lazarusA/helpfulJulia

max_issues_repo_name: lazarusA/helpfulJulia

max_forks_repo_name: lazarusA/helpfulJulia


text:
@doc raw"""
    CenteredMatrices{m,n,𝔽} <: AbstractDecoratorManifold{𝔽}

The manifold of $m × n$ real-valued or complex-valued matrices whose columns sum to zero, i.e.
````math
\bigl\{ p ∈ 𝔽^{m × n}\ \big|\ [1 … 1] * p = [0 … 0] \bigr\},
````
where $𝔽 ∈ \{ℝ,ℂ\}$.

# Constructor
    CenteredMatrices(m, n[, field=ℝ])

Generate the manifold of `m`-by-`n` (`field`-valued) matrices whose columns sum to zero.
"""
struct CenteredMatrices{M,N,𝔽} <: AbstractDecoratorManifold{𝔽} end

function CenteredMatrices(m::Int, n::Int, field::AbstractNumbers=ℝ)
    return CenteredMatrices{m,n,field}()
end

active_traits(f, ::CenteredMatrices, args...) = merge_traits(IsEmbeddedSubmanifold())

@doc raw"""
    check_point(M::CenteredMatrices{m,n,𝔽}, p; kwargs...)

Check whether the matrix is a valid point on the
[`CenteredMatrices`](@ref) `M`, i.e. is an `m`-by-`n` matrix whose columns sum to
zero.

The tolerance for the column sums of `p` can be set using `kwargs...`.
"""
function check_point(M::CenteredMatrices{m,n,𝔽}, p; kwargs...) where {m,n,𝔽}
    if !isapprox(sum(p, dims=1), zeros(1, n); kwargs...)
        return DomainError(
            p,
            string(
                "The point $(p) does not lie on $(M), since its columns do not sum to zero.",
            ),
        )
    end
    return nothing
end

"""
    check_vector(M::CenteredMatrices{m,n,𝔽}, p, X; kwargs... )

Check whether `X` is a tangent vector to manifold point `p` on the
[`CenteredMatrices`](@ref) `M`, i.e. that `X` is a matrix of size `(m, n)` whose columns
sum to zero and its values are from the correct [`AbstractNumbers`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#number-system).
The tolerance for the column sums of `p` and `X` can be set using `kwargs...`.
"""
function check_vector(M::CenteredMatrices{m,n,𝔽}, p, X; kwargs...) where {m,n,𝔽}
    if !isapprox(sum(X, dims=1), zeros(1, n); kwargs...)
        return DomainError(
            X,
            "The vector $(X) is not a tangent vector to $(p) on $(M), since its columns do not sum to zero.",
        )
    end
    return nothing
end

embed(::CenteredMatrices, p) = p
embed(::CenteredMatrices, p, X) = X

get_embedding(::CenteredMatrices{m,n,𝔽}) where {m,n,𝔽} = Euclidean(m, n; field=𝔽)

@doc raw"""
    manifold_dimension(M::CenteredMatrices{m,n,𝔽})

Return the manifold dimension of the [`CenteredMatrices`](@ref) `m`-by-`n` matrix `M` over the number system
`𝔽`, i.e.

````math
\dim(\mathcal M) = (m*n - n) \dim_ℝ 𝔽,
````
where $\dim_ℝ 𝔽$ is the [`real_dimension`](https://juliamanifolds.github.io/ManifoldsBase.jl/stable/types.html#ManifoldsBase.real_dimension-Tuple{ManifoldsBase.AbstractNumbers}) of `𝔽`.
"""
function manifold_dimension(::CenteredMatrices{m,n,𝔽}) where {m,n,𝔽}
    return (m * n - n) * real_dimension(𝔽)
end

@doc raw"""
    project(M::CenteredMatrices, p)

Projects `p` from the embedding onto the [`CenteredMatrices`](@ref) `M`, i.e.

````math
\operatorname{proj}_{\mathcal M}(p) = p - \begin{bmatrix}
1\\
⋮\\
1
\end{bmatrix} * [c_1 \dots c_n],
````
where $c_i = \frac{1}{m}\sum_{j=1}^m p_{j,i}$ for $i = 1, \dots, n$.
"""
project(::CenteredMatrices, ::Any)

project!(::CenteredMatrices, q, p) = copyto!(q, p .- mean(p, dims=1))

@doc raw"""
    project(M::CenteredMatrices, p, X)

Project the matrix `X` onto the tangent space at `p` on the [`CenteredMatrices`](@ref) `M`, i.e.

````math
\operatorname{proj}_p(X) = X - \begin{bmatrix}
1\\
⋮\\
1
\end{bmatrix} * [c_1 \dots c_n],
````
where $c_i = \frac{1}{m}\sum_{j=1}^m x_{j,i}$  for $i = 1, \dots, n$.
"""
project(::CenteredMatrices, ::Any, ::Any)

project!(::CenteredMatrices, Y, p, X) = (Y .= X .- mean(X, dims=1))

@generated representation_size(::CenteredMatrices{m,n,𝔽}) where {m,n,𝔽} = (m, n)

function Base.show(io::IO, ::CenteredMatrices{m,n,𝔽}) where {m,n,𝔽}
    return print(io, "CenteredMatrices($(m), $(n), $(𝔽))")
end


text:
# Sorting.

mutable struct OrderNode <: TabularNode
    over::Union{SQLNode, Nothing}
    by::Vector{SQLNode}

    OrderNode(; over = nothing, by) =
        new(over, by)
end

OrderNode(by...; over = nothing) =
    OrderNode(over = over, by = SQLNode[by...])

"""
    Order(; over = nothing, by)
    Order(by...; over = nothing)

`Order` sorts the input rows `by` the given key.

```sql
SELECT ...
FROM \$over
ORDER BY \$by...
```

# Examples

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :year_of_birth]);

julia> q = From(person) |>
           Order(Get.year_of_birth);

julia> print(render(q))
SELECT "person_1"."person_id", "person_1"."year_of_birth"
FROM "person" AS "person_1"
ORDER BY "person_1"."year_of_birth"
```
"""
Order(args...; kws...) =
    OrderNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Order), pats::Vector{Any}) =
    dissect(scr, OrderNode, pats)

function PrettyPrinting.quoteof(n::OrderNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(Order))
    if isempty(n.by)
        push!(ex.args, Expr(:kw, :by, Expr(:vect)))
    else
        append!(ex.args, quoteof(n.by, ctx))
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::OrderNode) =
    label(n.over)

rebase(n::OrderNode, n′) =
    OrderNode(over = rebase(n.over, n′), by = n.by)


text:
# Based on the Julia implementation
# InlineStringN are all `isbitstype` so they work with JDF automatically
# However for other languages, we may still need to explicitly support them and conver them
# to string is not available in those languages
using JDF
using InlineStrings, DataFrames, CSV

using Random:randstring

@testset "Test InlineStrings get loaded and saved properly" begin
    a = DataFrames.DataFrame(a = [randstring(254) |> String255 for i in 1:100])

    path = tempdir()

    CSV.write(joinpath(path, "tmp.csv"), a)

    a1 = CSV.read(joinpath(path, "tmp.csv"), DataFrame)

    a1.a = a1.a .|> String255

    JDF.save(joinpath(path, "tmp.jdf"), a1)

    a2 = JDF.load(joinpath(path, "tmp.jdf")) |> DataFrame

    @test eltype(a2.a) == String255

    # clean up
    rm("tmp.csv"; force=true)
    rm("tmp.jdf"; force=true, recursive=true)
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
function cr3bpEomsCondition(out, x, t, integrator)

    # Get Requirements
    TU  = integrator.p.crp.TU
    LU  = integrator.p.crp.LU
    MU  = integrator.p.crp.MU
    R1  = integrator.p.crp.R1
    R2  = integrator.p.crp.R2
    μ   = integrator.p.crp.μ
    isp = integrator.p.sp.isp
    m0  = integrator.p.sp.initMass
    mp  = integrator.p.sp.initProp

    ϵ   = integrator.p.ϵ
    utype = integrator.p.utype

    # Compute Scaling
    cSc = isp*9.81*TU / (LU * 1000.0)

    # Compute required
    λv = sqrt(x[11]^2 + x[12]^2 + x[13]^2)
    S  = computeS(x, λv, cSc)
    xpμ = x[1] + μ
    r1  = sqrt(xpμ^2 + x[2]^2 + x[3]^2)
    r2  = sqrt((xpμ - 1)^2 + x[2]^2 + x[3]^2)

    # Switching condition
    if ϵ != 0.0
        if  utype == 0
            out[1] = S - ϵ
        elseif utype == 2
            out[1] = S + ϵ
        else
            out[1] = abs(S) - ϵ
        end
    else
        out[1] = S
    end

    # Termination conditions
    out[2] = x[7]*MU - (m0 - mp)
    out[3] = r1*LU - R1
    out[4] = r2*LU - R2
end

function cr3bpEomsAffect!(integrator, idx)

    # Get Requirements
    TU  = integrator.p.crp.TU
    LU  = integrator.p.crp.LU
    isp = integrator.p.sp.isp

    ϵ   = integrator.p.ϵ
    utype = integrator.p.utype

    # Switching affect
    if idx == 1
        if ϵ != 0.0
            if utype != 1
                integrator.p.utype = 1
            else
                cSc = isp*9.81*TU / (LU*1000.0)
                λv = sqrt(integrator.u[11]*integrator.u[11] +
                          integrator.u[12]*integrator.u[12] +
                          integrator.u[13]*integrator.u[13])
                S = computeS(integrator.u, λv, cSc)
                if S < 0.0
                    integrator.p.utype = 2
                else
                    integrator.p.utype = 0
                end
            end

        else # Need to prop STM across switching
            if utype == 0
                integrator.p.utype = 2
            else
                integrator.p.utype = 0
            end
        end

        # Reset step
        auto_dt_reset!(integrator)

    # Termination affect
    else
        #println("Terminating!")
        terminate!(integrator)
    end
end

function cr3bpEomsConditionNoTerm(x, t, integrator)

    # Get Requirements
    TU  = integrator.p.crp.TU
    LU  = integrator.p.crp.LU
    isp = integrator.p.sp.isp

    ϵ   = integrator.p.ϵ
    utype = integrator.p.utype

    # Compute Scaling
    cSc = isp*9.81*TU / (LU * 1000.0)

    # Compute required
    λv = sqrt(x[11]^2 + x[12]^2 + x[13]^2)
    S  = computeS(x, λv, cSc)

    # Switching condition
    out = 0.0
    if ϵ != 0.0
        if  utype == 0
            out = S - ϵ
        elseif utype == 2
            out = S + ϵ
        else
            out = abs(S) - ϵ
        end
    else
        out = S
    end

    return out
end

function cr3bpEomsAffectNoTerm!(integrator)

    # Get Requirements
    TU  = integrator.p.crp.TU
    LU  = integrator.p.crp.LU
    isp = integrator.p.sp.isp

    ϵ   = integrator.p.ϵ
    utype = integrator.p.utype

    # Switching affect
    if ϵ != 0.0
        if utype != 1
            integrator.p.utype = 1
        else
            cSc = isp*9.81*TU / (LU*1000.0)
            λv = sqrt(integrator.u[11]*integrator.u[11] +
                      integrator.u[12]*integrator.u[12] +
                      integrator.u[13]*integrator.u[13])
            S = computeS(integrator.u, λv, cSc)
            if S < 0.0
                integrator.p.utype = 2
            else
                integrator.p.utype = 0
            end
        end

    else # Need to prop STM across switching
        if utype == 0
            integrator.p.utype = 2
        else
            integrator.p.utype = 0
        end
    end

    # Reset step
    auto_dt_reset!(integrator)
end

function cr3bpEomsAffectNoTermWithSTM!(integrator)

    # Get Requirements
    TU  = integrator.p.crp.TU
    LU  = integrator.p.crp.LU
    isp = integrator.p.sp.isp

    ϵ   = integrator.p.ϵ
    utype = integrator.p.utype

    # Switching affect
    if ϵ != 0.0
        if utype != 1
            integrator.p.utype = 1
        else
            cSc = isp*9.81*TU / (LU*1000.0)
            λv = sqrt(integrator.u[11]*integrator.u[11] +
                      integrator.u[12]*integrator.u[12] +
                      integrator.u[13]*integrator.u[13])
            S = computeS(integrator.u, λv, cSc)
            if S < 0.0
                integrator.p.utype = 2
            else
                integrator.p.utype = 0
            end
        end
    else
        propSTM!(integrator.u, integrator.p)
    end

    # Reset step
    auto_dt_reset!(integrator)
end

function propSTM!(u::AbstractVector, ps::CR3BPIndirectWithSTMParams)
    @inbounds begin
        # Get requirements
        TU  = ps.crp.TU
        LU  = ps.crp.LU
        MU  = ps.crp.MU
        μ   = ps.crp.μ
        c   = ps.sp.c

        # Scale Requirements
        tMaxSc = ps.sp.tMax * TU * TU / (MU*LU*1000.0)
        cSc = c*TU / (LU*1000.0)

        # Compute Requirements
        λv = sqrt(u[11]^2 + u[12]^2 + u[13]^2)

        # Get state/co-state derivatives before updating utype
        dy⁻ = cr3bpEomIndirect(view(u, 1:14), ps, 0.0, MEMF())

        # Update utype
        ps.utype == 0 ? ps.utype = 2 : ps.utype = 0

        # Get state/co-state derivatives after updating utype
        dy⁺ = cr3bpEomIndirect(view(u, 1:14), ps, 0.0, MEMF())

        # Compute dy difference
        dyDiff = @SVector [dy⁺[1]  - dy⁻[1],
                           dy⁺[2]  - dy⁻[2],
                           dy⁺[3]  - dy⁻[3],
                           dy⁺[4]  - dy⁻[4],
                           dy⁺[5]  - dy⁻[5],
                           dy⁺[6]  - dy⁻[6],
                           dy⁺[7]  - dy⁻[7],
                           dy⁺[8]  - dy⁻[8],
                           dy⁺[9]  - dy⁻[9],
                           dy⁺[10] - dy⁻[10],
                           dy⁺[11] - dy⁻[11],
                           dy⁺[12] - dy⁻[12],
                           dy⁺[13] - dy⁻[13],
                           dy⁺[14] - dy⁻[14]]

        # Compute 1 / Sdot and ∂S/∂y / Sdot
        dSInv   = (u[8] - 2.0*u[12])*u[11] + (u[9] + 2.0*u[11])*u[12] + u[10]*u[13]
        dSInv   = λv*u[7]/(cSc*dSInv)
        temp    = -cSc*dSInv / (u[7]*λv)
        ∂SdSInv = @SVector [0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                            λv*cSc*dSInv / (u[7]^2),
                            0.0, 0.0, 0.0,
                            u[11]*temp,
                            u[12]*temp,
                            u[13]*temp,
                            -1.0*dSInv]

        # Compute Ψ using allocated jac matrix in parameters
        Ψ = ps.m1
        vecOuterProdPlusI!(Ψ, dyDiff, ∂SdSInv)

        # Propagate STM
        Φ⁺ = reshape(view(u, 15:210), (14,14))
        Φ⁻ = ps.m2; Φ⁻ .= Φ⁺
        Octavian.matmul!(Φ⁺, Ψ, Φ⁻)
    end
    return nothing
end