text:
function hirasaki_model(;log_k=[0.30, 1.74, 1.62, 0.14, 0.5, 2.23, 1.0, 0.09, -0.64],
                        dh=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    return """
    SURFACE_MASTER_SPECIES
        Chalk_a    Chalk_aOH-0.75
        Chalk_c    Chalk_cOH+0.75
    SURFACE_SPECIES
        Chalk_aOH-0.75 = Chalk_aOH-0.75
        log_k   0.0
        Chalk_cOH+0.75 = Chalk_cOH+0.75
        log_k   0.0
        # calcium site
        Chalk_aOH-0.75 + H+ = Chalk_aOH2+0.25
        log_k   $(log_k[1]) # 0.30
        delta_h      $(dh[1])
        Chalk_aOH-0.75 + Ca+2 = Chalk_aOHCa+1.25
        log_k   $(log_k[2]) # 1.74
        delta_h      $(dh[2])
        Chalk_aOH-0.75 + Mg+2 = Chalk_aOHMg+1.25
        log_k   $(log_k[3]) # 1.62
        delta_h      $(dh[3])
        Chalk_aOH-0.75 + Na+ = Chalk_aOHNa+0.25
        log_k   $(log_k[4]) # 0.14
        delta_h      $(dh[4])
        # carbonate site
        Chalk_cOH+0.75 + OH- = Chalk_cO-0.25 + H2O
        log_k   $(log_k[5]) # 0.5
        delta_h      $(dh[5])
        Chalk_cOH+0.75 + CO3-2 = Chalk_cOHCO3-1.25
        log_k   $(log_k[6]) # 2.23
        delta_h      $(dh[6])
        Chalk_cOH+0.75 + SO4-2 = Chalk_cOHSO4-1.25
        log_k   $(log_k[7]) # 1.0
        delta_h      $(dh[7])
        Chalk_cOH+0.75 + HCO3- = Chalk_cOH2CO3-0.25
        log_k   $(log_k[8]) # 0.09
        delta_h      $(dh[8])
        Chalk_cOH+0.75 + Cl- = Chalk_cOHCl-0.25
        log_k   $(log_k[9]) # -0.64
        delta_h      $(dh[9])
    END
    """
end

function modified_heberling(;log_k=[0.501, 0.56, 1.68, 1.48, -0.05, 0.04, -0.07, 0.4, 0.56, 1.68, 1.48],
                        dh=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    return """
    SURFACE_MASTER_SPECIES
    Chalk_a Chalk_aOH1.5
    Chalk_c Chalk_cH0.5

    SURFACE_SPECIES
    Chalk_aOH1.5 = Chalk_aOH1.5
    -cd_music 0 0 0
    log_k 0.0
    Chalk_cH0.5 = Chalk_cH0.5
    -cd_music 0 0 0
    log_k 0.0

    #dummy reactions producing surface species with fractional charges
    # are these necessary? Is there a way around it?
    Chalk_aOH1.5 = Chalk_aOH-0.5 + 0.5 H+
    -cd_music -0.5 0 0
    log_k 20
    Chalk_cH0.5 = Chalk_c-0.5 + 0.5 H+
    -cd_music -0.5 0 0
    log_k 20

    #inner sphere reactions
    Chalk_aOH-0.5 + H+ = Chalk_aOH2+0.5
    -cd_music 1 0 0
    log_K $(log_k[1]) # 0.501
    delta_h      $(dh[1])

    #outer sphere reactions
    Chalk_aOH-0.5 + Na+ = Chalk_aOHNa+0.5
    -cd_music 0 1 0
    log_K $(log_k[2]) # 0.56
    delta_h      $(dh[2])
    Chalk_aOH-0.5 + Ca+2 = Chalk_aOHCa+1.5
    -cd_music 0 2 0
    log_K $(log_k[3]) # 1.68
    delta_h      $(dh[3])
    Chalk_aOH-0.5 + Mg+2 = Chalk_aOHMg+1.5
    -cd_music 0 2 0
    log_K $(log_k[4]) # 1.48
    delta_h      $(dh[4])
    Chalk_aOH2+0.5 + Cl- = Chalk_aOH2Cl-0.5
    -cd_music 0 -1 0
    log_K $(log_k[5]) # -0.05
    delta_h      $(dh[5])
    Chalk_aOH2+0.5 + HCO3- = Chalk_aOH3CO3-0.5
    -cd_music 0 -1 0
    log_K $(log_k[6]) # 0.04
    delta_h      $(dh[6])
    Chalk_aOH2+0.5 + HCO3- = Chalk_aOH2CO3-1.5+ H+
    -cd_music 0 -2 0
    log_K $(log_k[7]) # -7.07
    delta_h      $(dh[7])
    Chalk_aOH2+0.5 + SO4-2 = Chalk_aOH2SO4-1.5
    -cd_music 0 -2 0
    log_K $(log_k[8]) # 0.4
    delta_h      $(dh[8])
    Chalk_c-0.5 + Na+ = Chalk_cNa+0.5
    -cd_music 0 1 0
    log_K $(log_k[9]) #  0.56
    delta_h      $(dh[9])
    Chalk_c-0.5 + Ca+2 = Chalk_cCa+1.5
    -cd_music 0 2 0
    log_K $(log_k[10]) # 1.68
    delta_h      $(dh[10])
    Chalk_c-0.5 + Mg+2 = Chalk_cMg+1.5
    -cd_music 0 2 0
    log_K $(log_k[11]) # 1.48
    delta_h      $(dh[11])
    """
end

function modified_so(;log_k=[-3.52108, -2.61239, -2.72276, 14.5456, -24.73, 9.26094, 1.37606, 0.205842, -9.8144, -8.8383],
    dh=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    return """
    SURFACE_MASTER_SPECIES
        Chalk_a Chalk_aOH-0.667
        Chalk_c Chalk_cH+0.667
    SURFACE_SPECIES
        Chalk_cH+0.667 = Chalk_cH+0.667
        log_k 0.0
        -cd_music 0 0 0 0 0
        Chalk_aOH-0.667 = Chalk_aOH-0.667
        log_k 0.0
        -cd_music 0 0 0 0 0
        Chalk_cH+0.667 = Chalk_c-0.333 + H+
        log_k $(log_k[1]) # -3.58
        delta_h      $(dh[1])
        -cd_music -1 0 0 0 0
        Chalk_cH+0.667 + Ca+2 = Chalk_cCa+1.667 + H+
        log_k $(log_k[2]) # -2.8
        delta_h      $(dh[2])
        -cd_music -1 2 0 0 0
        Chalk_cH+0.667 + Mg+2 = Chalk_cMg+1.667 + H+
        log_k $(log_k[3]) # -2.6
        delta_h      $(dh[3])
        -cd_music -1 2 0 0 0
        Chalk_aOH-0.667 + H+ = Chalk_aOH2+0.333
        log_k $(log_k[4]) # 12.85
        delta_h      $(dh[4])
        -cd_music 1 0 0 0 0
        Chalk_aOH-0.667 = Chalk_aO-1.667 + H+
        log_k $(log_k[5]) # -24.73
        delta_h      $(dh[5])
        -cd_music -1 0 0 0 0
        Chalk_aOH-0.667 + CO3-2 + H+ = Chalk_aHCO3-0.667 + OH-
        log_k $(log_k[6]) # 10.15
        delta_h      $(dh[6])
        -cd_music 0.6 -0.6 0 0 0
        Chalk_aOH-0.667 + CO3-2 = Chalk_aCO3-1.667 + OH-
        log_k $(log_k[7]) # 1.55
        delta_h      $(dh[7])
        -cd_music 0.6 -1.6 0 0 0
        Chalk_aOH-0.667 + SO4-2 = Chalk_aSO4-1.667 + OH-
        log_k $(log_k[8]) # 1.55
        delta_h      $(dh[8])
        -cd_music 0.6 -1.6 0 0 0
        Chalk_cH+0.667 + Na+ = Chalk_cNa+0.667+H+
        log_k $(log_k[9])  #-9.11
        delta_h      $(dh[9])
        -cd_music -1 1 0 0 0
        Chalk_aOH-0.667 + Cl- = Chalk_aCl-0.667 + OH-
        log_k $(log_k[10])  #-5.0
        delta_h      $(dh[10])
        -cd_music 1 -1 0 0 0

    """
end

function modified_wolthers(;log_k=[-3.58, -2.8, -2.2, 12.85, -24.73, 10.15, 1.55, 0.35],
    dh=[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
    return """
    SURFACE_MASTER_SPECIES
        Chalk_a Chalk_aOH-0.667
        Chalk_c Chalk_cH+0.667
    SURFACE_SPECIES
        Chalk_cH+0.667 = Chalk_cH+0.667
        log_k 0.0
        -cd_music 0 0 0 0 0
        Chalk_aOH-0.667 = Chalk_aOH-0.667
        log_k 0.0
        -cd_music 0 0 0 0 0
        Chalk_cH+0.667 = Chalk_c-0.333 + H+
        log_k $(log_k[1]) # -3.58
        delta_h      $(dh[1])
        -cd_music -1 0 0 0 0
        Chalk_cH+0.667 + Ca+2 = Chalk_cCa+1.667 + H+
        log_k $(log_k[2]) # -2.8
        delta_h      $(dh[2])
        -cd_music -1 2 0 0 0
        Chalk_cH+0.667 + Mg+2 = Chalk_cMg+1.667 + H+
        log_k $(log_k[3]) # -2.6
        delta_h      $(dh[3])
        -cd_music -1 2 0 0 0
        Chalk_aOH-0.667 + H+ = Chalk_aOH2+0.333
        log_k $(log_k[4]) # 12.85
        delta_h      $(dh[4])
        -cd_music 1 0 0 0 0
        Chalk_aOH-0.667 = Chalk_aO-1.667 + H+
        log_k $(log_k[5]) # -24.73
        delta_h      $(dh[5])
        -cd_music -1 0 0 0 0
        Chalk_aOH-0.667 + CO3-2 + H+ = Chalk_aHCO3-0.667 + OH-
        log_k $(log_k[6]) # 10.15
        delta_h      $(dh[6])
        -cd_music 0.6 -0.6 0 0 0
        Chalk_aOH-0.667 + CO3-2 = Chalk_aCO3-1.667 + OH-
        log_k $(log_k[7]) # 1.55
        delta_h      $(dh[7])
        -cd_music 0.6 -1.6 0 0 0
        Chalk_aOH-0.667 + SO4-2 = Chalk_aSO4-1.667 + OH-
        log_k $(log_k[8]) # 1.55
        delta_h      $(dh[8])
        -cd_music 0.6 -1.6 0 0 0

    """
end

max_stars_repo_name: simulkade/calciteSCM

max_issues_repo_name: simulkade/calciteSCM

max_forks_repo_name: simulkade/calciteSCM


text:
export conv, conv!, ∇conv_data, ∇conv_data!, ∇conv_filter, ∇conv_filter!, depthwiseconv,
        depthwiseconv!, ∇depthwiseconv_data, ∇depthwiseconv_data!, ∇depthwiseconv_filter,
        ∇depthwiseconv_filter!

## Convolution API
#
#  We provide the following generic methods, for 3d, 4d, and 5d tensors, calculating 1d,
#  2d and 3d convolutions, based on the rank of the input tensors, in both mutating and
#  non-mutating auto-allocating variants:
#   - Convolution:
#     - conv(x, w, cdims)
#     - conv!(y, x, w, cdims)
#   - Convolution data backpropagation
#     - ∇conv_data(dy, w, cdims)
#     - ∇conv_data!(dx, dy, w, cdims)
#   - Convolution filter backpropagation
#     - ∇conv_filter(x, dy, cdims)
#     - ∇conv_filter!(dw, x, dy, cdims)
#
#   All methods require a `ConvDims` object to define the dimensions and optional
#   elements of the convolution (padding, stride, dilation, kernel-flipping, etc...),
#   which is easily constructable through something like `DenseConvDims(x, w)`.  All
#   methods take in the `ConvDims` of the associated normal, forward-pass convolution,
#   that is, the following is legal:
#
#       cdims = ConvDims(x, w; stride=2, dilation=(3,2))
#       dx = ∇conv_data(conv(x, w, cdims), w, cdims)

#   The computational flow, starting from the user facing functions,
#   goes through the following steps:
#
#   STEP 1:
#       use ConvDims objects (only for `conv` and `depthwiseconv`)
#   STEP 2:
#        define autoallocating version (frontend and implementations)
#   STEP 3:
#        reshape to 3d convolutions (frontend and implementions)
#   STEP 4:
#        choose implementation

# TODO: should we also add
#   STEP X:
#        use homogeneus datatypes
# to handle etherogeneus inputs now handled by conv_direct?

########## STEP 1 ############
"""
    conv(x, w; stride=1, pad=0, dilation=1, flipped=false)

Apply convolution filter `w` to input `x`. `x` and `w` are 3d/4d/5d tensors
in 1d/2d/3d convolutions respectively.
"""
function conv(x, w::AbstractArray{T, N}; stride=1, pad=0, dilation=1, flipped=false) where {T, N}
    stride = expand(Val(N-2), stride)
    pad = expand(Val(N-2), pad)
    dilation = expand(Val(N-2), dilation)
    cdims = DenseConvDims(x, w; stride=stride, padding=pad, dilation=dilation, flipkernel=flipped)
    return conv(x, w, cdims)
end

"""
    depthwiseconv(x, w; stride=1, pad=0, dilation=1, flipped=false)

Depthwise convolution operation with filter `w` on input `x`. `x` and `w`
are 3d/4d/5d tensors in 1d/2d/3d convolutions respectively.
"""
function depthwiseconv(x, w::AbstractArray{T, N}; stride=1, pad=0, dilation=1, flipped=false) where {T, N}
    stride = expand(Val(N-2), stride)
    pad = expand(Val(N-2), pad)
    dilation = expand(Val(N-2), dilation)
    cdims = DepthwiseConvDims(x, w; stride=stride, padding=pad, dilation=dilation, flipkernel=flipped)
    return depthwiseconv(x, w, cdims)
end
##############################

########### STEP 2 ###################
# Let's generate auto-allocating versions of all our functions, for all backends.
# We `@timeit` these methods separately, as we want to know how much time is spent in
# allocation.  :P
for backend in (Symbol(), :_direct, :_im2col, :_nnpack)
    # First make auto-allocating versions of the conv()-like calls:
    for name in (:conv, :depthwiseconv)
        @eval begin
            function $(Symbol("$(name)$(backend)"))(
                            x::AbstractArray{xT,N}, w::AbstractArray{wT,N},
                            cdims::ConvDims; kwargs...) where {xT, wT, N}
                y = similar(x, promote_type(xT, wT), output_size(cdims)...,
                               channels_out(cdims), size(x,N))
                return $(Symbol("$(name)$(backend)!"))(y, x, w, cdims; kwargs...)
            end
        end
    end

    for name in (:∇conv_data, :∇depthwiseconv_data)
        @eval begin
            function $(Symbol("$(name)$(backend)"))(
                            dy::AbstractArray{yT,N}, w::AbstractArray{wT,N},
                            cdims::ConvDims; kwargs...) where {yT, wT, N}
                dx = similar(dy, input_size(cdims)..., channels_in(cdims),
                                                        size(dy, N))
                return $(Symbol("$(name)$(backend)!"))(dx, dy, w, cdims; kwargs...)
            end
        end
    end

    # We do the conv/depthwiseconv filter backprops separately, as the shape calculation
    # for `w` is slightly different for depthwise than for normal dense convolution.
    @eval begin
        function $(Symbol("∇conv_filter$(backend)"))(
                        x::AbstractArray{xT,N}, dy::AbstractArray{yT,N},
                        cdims::ConvDims; kwargs...) where {xT, yT, N}
            dw = similar(dy, kernel_size(cdims)..., channels_in(cdims),
                                                    channels_out(cdims))
            return $(Symbol("∇conv_filter$(backend)!"))(dw, x, dy, cdims; kwargs...)
        end
    end

    @eval begin
        function $(Symbol("∇depthwiseconv_filter$(backend)"))(
                        x::AbstractArray{xT,N}, dy::AbstractArray{yT,N},
                        cdims::ConvDims; kwargs...) where {xT, yT, N}
            dw = similar(dy, kernel_size(cdims)..., channel_multiplier(cdims),
                                                    channels_in(cdims))
            return $(Symbol("∇depthwiseconv_filter$(backend)!"))(dw, x, dy, cdims;
                                                                 kwargs...)
        end
    end
end
##########################################

########## STEP 3 ############

# Our strategy for 1d and 2d convolution is to reshape to 3d convolutions, which
# makes things MUCH EASIER for us on the backend side, and is in general pretty fast,
# since we can specialize on sizes.
for front_name in (:conv, :∇conv_data, :∇conv_filter,
                   :depthwiseconv, :∇depthwiseconv_data, :∇depthwiseconv_filter)
    for backend in (Symbol(), :_direct, :_im2col) ## NNPACK  is only for 2d conv
        for N in (3, 4)
            @eval begin
                function $(Symbol("$(front_name)$(backend)!"))(
                                y::AbstractArray{yT,$N}, x::AbstractArray{xT,$N},
                                w::AbstractArray{wT,$N}, cdims::ConvDims;
                                kwargs...) where {yT, xT, wT}
                    $(Symbol("$(front_name)$(backend)!"))(
                        insert_singleton_spatial_dimension(y, $(5 - N)),
                        insert_singleton_spatial_dimension(x, $(5 - N)),
                        insert_singleton_spatial_dimension(w, $(5 - N)),
                        insert_singleton_spatial_dimension(cdims, $(5 - N));
                        kwargs...
                    )

                    # We explicitly return `y` here, because the backend call
                    # itself may return a reshaped view, which we don't want.
                    return y
                end
            end
        end
    end
end
#######################################

########### STEP 4 ############

# First, we will define mappings from the generic API names to our accelerated backend
# implementations. For homogeneous-datatype 1, 2 and 3d convolutions, we default to using
# im2col + GEMM.  Do so in a loop, here:
for (front_name, backend) in (
        # This maps from public, front-facing name, to internal backend name
        :conv                   => :im2col,
        :∇conv_data             => :im2col,
        :∇conv_filter           => :im2col,
        :depthwiseconv          => :im2col,
        :∇depthwiseconv_data    => :im2col,
        :∇depthwiseconv_filter  => :im2col,
    )

    # These are the GEMM types we will accelerate with `im2col`
    G = Union{[x[2] for x in gemm_datatype_mappings]...}

    # We only define 3d conv primitives, we reshape lower down to get 1d and 2d convolution
    @eval begin
        # im2col-accelerated function forwarding definition
        function $(Symbol("$(front_name)!"))(
                        out::AbstractArray{T,5}, in1::AbstractArray{T,5},
                        in2::AbstractArray{T,5}, cdims::ConvDims; kwargs...) where {T <: $G}
            $(Symbol("$(front_name)_$(backend)!"))(out, in1, in2, cdims; kwargs...)
        end
    end
end

# We always support a fallback, non-accelerated path, where we use the direct, but
# slow, implementations.  These should not typically be used, hence the `@warn`,
# but let's go ahead and define them first:
for front_name in (:conv, :∇conv_data, :∇conv_filter,
                   :depthwiseconv, :∇depthwiseconv_data, :∇depthwiseconv_filter)
    @eval begin
        function $(Symbol("$(front_name)!"))(
                        y::AbstractArray{yT,N}, in1::AbstractArray{T1,N},
                        in2::AbstractArray{T2,N}, cdims::ConvDims;
                        kwargs...) where {yT, T1, T2, N}
            @warn string("Slow fallback implementation invoked for ", $(string(front_name)), "!  ",
                          "You probably don't want this; check your datatypes.") yT T1 T2 maxlog=1
            $(Symbol("$(front_name)_direct!"))(y, in1, in2, cdims; kwargs...)
        end
    end
end

# Use NNPACK if it is available and the operation is supported
# commented out 'till proper benchmarking and more correctness test are performed
# if is_nnpack_available()
#     function conv(x::Array{Float32, 4}, w::Array{Float32, 4},
#                   cdims::DenseConvDims{2, K, C_in, C_out, (1, 1), P, (1, 1), F};
#                   kwargs...) where {K, C_in, C_out, P, F}
#         return conv_nnpack(x, w, cdims; kwargs...)
#     end

#     function ∇conv_data(dy::Array{Float32, 4}, w::Array{Float32, 4},
#                 cdims::DenseConvDims{2, K, C_in, C_out, (1, 1), P, (1, 1), F};
#                 kwargs...) where {K, C_in, C_out, P, F}
#         return ∇conv_data_nnpack(dy, w, cdims; kwargs...)
#     end

#     function ∇conv_filter(x::Array{Float32, 4}, dy::Array{Float32, 4},
#                 cdims::DenseConvDims{2, K, C_in, C_out, (1, 1), P, (1, 1), F};
#                 kwargs...) where {K, C_in, C_out, P, F}
#         return ∇conv_filter_nnpack(x, dy, cdims; kwargs...)
#     end
# end
########################################################

ext: jl


text:
using SFML

# This program draws the SFML.jl logo and applies a blur shader to it

window = RenderWindow("Shader Test", 800, 600)

event = Event()

texture = Texture(joinpath(dirname(@__FILE__),"..","..","assets","sfmljl_logo.png"))
sprite = Sprite()
set_texture(sprite, texture)
set_position(sprite, Vector2f(400, 300))
scale(sprite, Vector2f(0.25, 0.25))
size = get_size(texture)
set_origin(sprite, Vector2f(size.x/2, size.y/2))

vert =
"""
uniform float wave_phase;
uniform vec2 wave_amplitude;

void main()
{
    vec4 vertex = gl_Vertex;
    vertex.x += cos(gl_Vertex.y * 0.02 + wave_phase * 3.8) * wave_amplitude.x
              + sin(gl_Vertex.y * 0.02 + wave_phase * 6.3) * wave_amplitude.x * 0.3;
    vertex.y += sin(gl_Vertex.x * 0.02 + wave_phase * 2.4) * wave_amplitude.y
              + cos(gl_Vertex.x * 0.02 + wave_phase * 5.2) * wave_amplitude.y * 0.3;

    gl_Position = gl_ModelViewProjectionMatrix * vertex;
    gl_TexCoord[0] = gl_TextureMatrix[0] * gl_MultiTexCoord0;
    gl_FrontColor = gl_Color;
}
"""

frag =
"""
uniform sampler2D texture;
uniform float blur_radius;

void main()
{
    vec2 offx = vec2(blur_radius, 0.0);
    vec2 offy = vec2(0.0, blur_radius);

    vec4 pixel = texture2D(texture, gl_TexCoord[0].xy)               * 4.0 +
                 texture2D(texture, gl_TexCoord[0].xy - offx)        * 2.0 +
                 texture2D(texture, gl_TexCoord[0].xy + offx)        * 2.0 +
                 texture2D(texture, gl_TexCoord[0].xy - offy)        * 2.0 +
                 texture2D(texture, gl_TexCoord[0].xy + offy)        * 2.0 +
                 texture2D(texture, gl_TexCoord[0].xy - offx - offy) * 1.0 +
                 texture2D(texture, gl_TexCoord[0].xy - offx + offy) * 1.0 +
                 texture2D(texture, gl_TexCoord[0].xy + offx - offy) * 1.0 +
                 texture2D(texture, gl_TexCoord[0].xy + offx + offy) * 1.0;

    gl_FragColor =  gl_Color * (pixel / 16.0);
}
"""

shader = ShaderFromMemory(vert, frag)
set_parameter(shader, "texture")

state = RenderStates(shader)

clock = Clock()

while isopen(window)
    while pollevent(window, event)
        if get_type(event) == EventType.CLOSED
            close(window)
        end
    end

    x = get_mousepos(window).x / get_size(window).x
    y = get_mousepos(window).y / get_size(window).y

    set_parameter(shader, "blur_radius", (x + y) * 0.008)
    set_parameter(shader, "wave_phase", get_elapsed_time(clock) |> as_seconds)
    set_parameter(shader, "wave_amplitude", Vector2f(x * 40, y * 40))

    clear(window, SFML.black)
    draw(window, sprite, state)
    display(window)
end


text:
# From https://github.com/JuliaInterop/RCall.jl/blob/master/src/types.jl

abstract type Sxp end # SEXPREC
const SxpPtrInfo = UInt32 # sxpinfo_struct

"R Sxp header: a pointer to this is used for unknown types."
struct SxpHead <: Sxp # SEXPREC_HEADER
    info::SxpPtrInfo
    attrib::Ptr{Cvoid}
    gc_next::Ptr{Cvoid}
    gc_prev::Ptr{Cvoid}
end
const UnknownSxp = SxpHead

abstract type VectorSxp <: Sxp end
abstract type VectorAtomicSxp <: VectorSxp end
abstract type VectorNumericSxp <: VectorAtomicSxp end
abstract type VectorListSxp <: VectorSxp end
abstract type PairListSxp <: Sxp end
abstract type FunctionSxp <: Sxp end

"R NULL value"
struct NilSxp <: PairListSxp   # type tag 0
    head::SxpHead
end

"R pairs (cons) list cell"
struct ListSxp <: PairListSxp  # type tag 2
    head::SxpHead
    car
    cdr
    tag
end

"R function closure"
struct ClosSxp <: FunctionSxp  # type tag 3
    head::SxpHead
    formals::Ptr{ListSxp}
    body
    env
end

"R environment"
struct EnvSxp <: Sxp  # type tag 4
    head::SxpHead
    frame
    enclos
    hashtab
end

"R promise"
struct PromSxp <: Sxp  # type tag 5
    head::SxpHead
    value
    expr
    env
end

"R function call"
struct LangSxp <: PairListSxp  # type tag 6
    head::SxpHead
    car
    cdr
    tag
end

"R special function"
struct SpecialSxp <: FunctionSxp  # type tag 7
    head::SxpHead
end

"R built-in function"
struct BuiltinSxp <: FunctionSxp  # type tag 8
    head::SxpHead
end

"R character string"
struct CharSxp <: VectorAtomicSxp     # type tag 9
    head::SxpHead
    length::Int
    truelength::Int
end

"R symbol"
struct SymSxp <: Sxp   # type tag 1
    head::SxpHead
    name::Ptr{CharSxp}
    value
    internal
end

"R logical vector"
struct LglSxp <: VectorNumericSxp     # type tag 10
    head::SxpHead
    length::Int
    truelength::Int
end

"R integer vector"
struct IntSxp <: VectorNumericSxp     # type tag 13
    head::SxpHead
    length::Int
    truelength::Int
end

"R real vector"
struct RealSxp <: VectorNumericSxp    # type tag 14
    head::SxpHead
    length::Int
    truelength::Int
end

"R complex vector"
struct CplxSxp <: VectorNumericSxp    # type tag 15
    head::SxpHead
    length::Int
    truelength::Int
end

"R vector of character strings"
struct StrSxp <: VectorListSxp     # type tag 16
    head::SxpHead
    length::Int
    truelength::Int
end

"R dot-dot-dot object"
struct DotSxp <: Sxp     # type tag 17
    head::SxpHead
end

"R \"any\" object"
struct AnySxp <: Sxp     # type tag 18
    head::SxpHead
end

"R list (i.e. Array{Any,1})"
struct VecSxp <: VectorListSxp     # type tag 19
    head::SxpHead
    length::Int
    truelength::Int
end

"R expression vector"
struct ExprSxp <: VectorListSxp    # type tag 20
    head::SxpHead
    length::Int
    truelength::Int
end

"R byte code"
struct BcodeSxp <: Sxp   # type tag 21
    head::SxpHead
end

"R external pointer"
struct ExtPtrSxp <: Sxp  # type tag 22
    head::SxpHead
    ptr::Ptr{Cvoid}
    prot::Ptr{Cvoid}
    tag
end

"R weak reference"
struct WeakRefSxp <: Sxp  # type tag 23
    head::SxpHead
end

"R byte vector"
struct RawSxp <: VectorAtomicSxp      # type tag 24
    head::SxpHead
    length::Int
    truelength::Int
end

"R S4 object"
struct S4Sxp <: Sxp      # type tag 25
    head::SxpHead
end


text:
"""
    TransformedKernel(k::Kernel, t::Transform)

Kernel derived from `k` for which inputs are transformed via a [`Transform`](@ref) `t`.

It is preferred to create kernels with input transformations with [`transform`](@ref)
instead of  `TransformedKernel` directly since [`transform`](@ref) allows optimized
implementations for specific kernels and transformations.

# Definition

For inputs ``x, x'``, the transformed kernel ``\\widetilde{k}`` derived from kernel ``k`` by
input transformation ``t`` is defined as
```math
\\widetilde{k}(x, x'; k, t) = k\\big(t(x), t(x')\\big).
```
"""
struct TransformedKernel{Tk<:Kernel,Tr<:Transform} <: Kernel
    kernel::Tk
    transform::Tr
end

@functor TransformedKernel

(k::TransformedKernel)(x, y) = k.kernel(k.transform(x), k.transform(y))

# Optimizations for scale transforms of simple kernels to save allocations:
# Instead of a multiplying every element of the inputs before evaluating the metric,
# we perform a scalar multiplcation of the distance of the original inputs, if possible.
function (k::TransformedKernel{<:SimpleKernel,<:ScaleTransform})(
    x::AbstractVector{<:Real}, y::AbstractVector{<:Real}
)
    return kappa(k.kernel, _scale(k.transform, metric(k.kernel), x, y))
end

function _scale(t::ScaleTransform, metric::Euclidean, x, y)
    return first(t.s) * evaluate(metric, x, y)
end
function _scale(t::ScaleTransform, metric::Union{SqEuclidean,DotProduct}, x, y)
    return first(t.s)^2 * evaluate(metric, x, y)
end
_scale(t::ScaleTransform, metric, x, y) = evaluate(metric, t(x), t(y))

"""
    transform(k::Kernel, t::Transform)

Create a [`TransformedKernel`](@ref) for kernel `k` and transform `t`.
"""
transform(k::Kernel, t::Transform) = TransformedKernel(k, t)
function transform(k::TransformedKernel, t::Transform)
    return TransformedKernel(k.kernel, t ∘ k.transform)
end

"""
    transform(k::Kernel, ρ::Real)

Create a [`TransformedKernel`](@ref) for kernel `k` and inverse lengthscale `ρ`.
"""
transform(k::Kernel, ρ::Real) = transform(k, ScaleTransform(ρ))

"""
    transform(k::Kernel, ρ::AbstractVector)

Create a [`TransformedKernel`](@ref) for kernel `k` and inverse lengthscales `ρ`.
"""
transform(k::Kernel, ρ::AbstractVector) = transform(k, ARDTransform(ρ))

kernel(κ) = κ.kernel

Base.show(io::IO, κ::TransformedKernel) = printshifted(io, κ, 0)

function printshifted(io::IO, κ::TransformedKernel, shift::Int)
    printshifted(io, κ.kernel, shift)
    return print(io, "\n" * ("\t"^(shift + 1)) * "- $(κ.transform)")
end

# Kernel matrix operations

function kernelmatrix_diag!(K::AbstractVector, κ::TransformedKernel, x::AbstractVector)
    return kernelmatrix_diag!(K, κ.kernel, _map(κ.transform, x))
end

function kernelmatrix!(K::AbstractMatrix, κ::TransformedKernel, x::AbstractVector)
    return kernelmatrix!(K, kernel(κ), _map(κ.transform, x))
end

function kernelmatrix!(
    K::AbstractMatrix, κ::TransformedKernel, x::AbstractVector, y::AbstractVector
)
    return kernelmatrix!(K, kernel(κ), _map(κ.transform, x), _map(κ.transform, y))
end

function kernelmatrix_diag(κ::TransformedKernel, x::AbstractVector)
    return kernelmatrix_diag(κ.kernel, _map(κ.transform, x))
end

function kernelmatrix(κ::TransformedKernel, x::AbstractVector)
    return kernelmatrix(kernel(κ), _map(κ.transform, x))
end

function kernelmatrix(κ::TransformedKernel, x::AbstractVector, y::AbstractVector)
    return kernelmatrix(kernel(κ), _map(κ.transform, x), _map(κ.transform, y))
end


text:
@testset "Spin--orbit interaction" begin
    @testset "Quantum numbers" begin
        orbitals = sos"1[s]"
        @test Atoms.mⱼ(orbitals[1]) == half(1)
        @test Atoms.mⱼ(orbitals[2]) == half(-1)

        orbitals = sos"2[p]"
        @test Atoms.mⱼ(orbitals[1]) == half(-1)
        @test Atoms.mⱼ(orbitals[2]) == half(-3)
        @test Atoms.mⱼ(orbitals[3]) == half(1)
        @test Atoms.mⱼ(orbitals[4]) == half(-1)
        @test Atoms.mⱼ(orbitals[5]) == half(3)
        @test Atoms.mⱼ(orbitals[6]) == half(1)
    end
end


text:
# interfacing with other packages

using LinearAlgebra
using LinearAlgebra: BlasComplex, BlasFloat, BlasReal

function mv_wrapper(transa::SparseChar, alpha::Number, A::CuSparseMatrix{T}, X::DenseCuVector{T},
                    beta::Number, Y::CuVector{T}) where {T}
    mv!(transa, alpha, A, X, beta, Y, 'O')
end

function mm_wrapper(transa::SparseChar, transb::SparseChar, alpha::Number,
                    A::CuSparseMatrix{T}, B::CuMatrix{T}, beta::Number, C::CuMatrix{T}) where {T}
    if version() <= v"10.3.1"
        # Generic mm! doesn't support transposed B on CUDA10
        return mm2!(transa, transb, alpha, A, B, beta, C, 'O')
    end
    mm!(transa, transb, alpha, A, B, beta, C, 'O')
end

tag_wrappers = ((identity, identity),
                (T -> :(HermOrSym{T, <:$T}), A -> :(parent($A))))
op_wrappers = (
    (identity, 'N', identity),
    (T -> :(Transpose{<:T, <:$T}), 'T', A -> :(parent($A))),
    (T -> :(Adjoint{<:T, <:$T}), 'C', A -> :(parent($A)))
)
for (taga, untaga) in tag_wrappers, (wrapa, transa, unwrapa) in op_wrappers
    TypeA = wrapa(taga(:(CuSparseMatrix{T})))

    @eval begin
        function LinearAlgebra.mul!(C::CuVector{T}, A::$TypeA, B::DenseCuVector{T},
                                    alpha::Number, beta::Number) where {T <: BlasFloat}
            mv_wrapper($transa, alpha, $(untaga(unwrapa(:A))), B, beta, C)
        end
    end

    for (tagb, untagb) in tag_wrappers, (wrapb, transb, unwrapb) in op_wrappers
        TypeB = wrapb(tagb(:(DenseCuMatrix{T})))

        isadjoint(expr) = expr.head == :curly && expr.args[1] == :Adjoint
        if isadjoint(TypeA) || isadjoint(TypeB)
            # CUSPARSE defines adjoints only for complex inputs. For real inputs we run
            # adjoints as tranposes so that we can still support the whole API surface of
            # LinearAlgebra.
            @eval begin
                function LinearAlgebra.mul!(C::CuMatrix{T}, A::$TypeA, B::$TypeB,
                                            alpha::Number, beta::Number) where {T <: BlasComplex}
                    mm_wrapper($transa, $transb, alpha, $(untaga(unwrapa(:A))),
                               $(untagb(unwrapb(:B))), beta, C)
                end
            end

            transa_real = transa == 'C' ? 'T' : transa
            transb_real = transb == 'C' ? 'T' : transb
            @eval begin
                function LinearAlgebra.mul!(C::CuMatrix{T}, A::$TypeA, B::$TypeB,
                                            alpha::Number, beta::Number) where {T <: BlasReal}
                    mm_wrapper($transa_real, $transb_real, alpha, $(untaga(unwrapa(:A))),
                               $(untagb(unwrapb(:B))), beta, C)
                end
            end
        else
            @eval begin
                function LinearAlgebra.mul!(C::CuMatrix{T}, A::$TypeA, B::$TypeB,
                                            alpha::Number, beta::Number) where {T <: BlasFloat}
                    mm_wrapper($transa, $transb, alpha, $(untaga(unwrapa(:A))),
                               $(untagb(unwrapb(:B))), beta, C)
                end
            end
        end
    end
end

Base.:(+)(A::Union{CuSparseMatrixCSR,CuSparseMatrixCSC},
          B::Union{CuSparseMatrixCSR,CuSparseMatrixCSC}) = geam(A,B,'O','O','O')
Base.:(-)(A::Union{CuSparseMatrixCSR,CuSparseMatrixCSC},
          B::Union{CuSparseMatrixCSR,CuSparseMatrixCSC}) = geam(A,-one(eltype(A)),B,'O','O','O')

# triangular

## direct
for (t, uploc, isunitc) in ((:LowerTriangular, 'L', 'N'),
                            (:UnitLowerTriangular, 'L', 'U'),
                            (:UpperTriangular, 'U', 'N'),
                            (:UnitUpperTriangular, 'U', 'U'))
    @eval begin
        # Left division
        LinearAlgebra.ldiv!(A::$t{T,<:AbstractCuSparseMatrix},
                            B::DenseCuVector{T}) where {T<:BlasFloat} =
            sv2!('N', $uploc, $isunitc, one(T), parent(A), B, 'O')

        LinearAlgebra.ldiv!(A::$t{T,<:AbstractCuSparseMatrix},
                            B::DenseCuMatrix{T}) where {T<:BlasFloat} =
            sm2!('N', 'N', $uploc, $isunitc, one(T), parent(A), B, 'O')
    end
end

## adjoint/transpose ('uploc' reversed)
for (t, uploc, isunitc) in ((:LowerTriangular, 'U', 'N'),
                            (:UnitLowerTriangular, 'U', 'U'),
                            (:UpperTriangular, 'L', 'N'),
                            (:UnitUpperTriangular, 'L', 'U'))
    @eval begin
        # Left division with vectors
        LinearAlgebra.ldiv!(A::$t{<:Any,<:Transpose{T,<:AbstractCuSparseMatrix}},
                            B::DenseCuVector{T}) where {T<:BlasFloat} =
            sv2!('T', $uploc, $isunitc, one(T), parent(parent(A)), B, 'O')

        LinearAlgebra.ldiv!(A::$t{<:Any,<:Adjoint{T,<:AbstractCuSparseMatrix}},
                            B::DenseCuVector{T}) where {T<:BlasReal} =
            sv2!('T', $uploc, $isunitc, one(T), parent(parent(A)), B, 'O')

        LinearAlgebra.ldiv!(A::$t{<:Any,<:Adjoint{T,<:AbstractCuSparseMatrix}},
                            B::DenseCuVector{T}) where {T<:BlasComplex} =
            sv2!('C', $uploc, $isunitc, one(T), parent(parent(A)), B, 'O')

        # Left division with matrices
        LinearAlgebra.ldiv!(A::$t{<:Any,<:Transpose{T,<:AbstractCuSparseMatrix}},
                            B::DenseCuMatrix{T}) where {T<:BlasFloat} =
            sm2!('T', 'N', $uploc, $isunitc, one(T), parent(parent(A)), B, 'O')

        LinearAlgebra.ldiv!(A::$t{<:Any,<:Adjoint{T,<:AbstractCuSparseMatrix}},
                            B::DenseCuMatrix{T}) where {T<:BlasReal} =
            sm2!('T', 'N', $uploc, $isunitc, one(T), parent(parent(A)), B, 'O')

        LinearAlgebra.ldiv!(A::$t{<:Any,<:Adjoint{T,<:AbstractCuSparseMatrix}},
                            B::DenseCuMatrix{T}) where {T<:BlasComplex} =
            sm2!('C', 'N', $uploc, $isunitc, one(T), parent(parent(A)), B, 'O')
    end
end


text:
using EarthDataLab
using Test
import DataFrames: DataFrame,aggregate
import Dates: year
import Statistics: cor, mean

@testset "Dataframe representation" begin
c=Cube()

d = getCubeData(c,variable=["air_temperature_2m","gross_primary_productivity"],longitude=(10,11),latitude=(50,51),
                time=(Date("2002-01-01"),Date("2008-12-31")))

dmem=readcubedata(d)

function docor(xout,xin)
    #Inside this function, xin is now a data frame
    @test isa(xin,DataFrame)
    xout[1]=cor(xin[:air_temperature_2m],xin[:gross_primary_productivity])
end
indims = InDims(TimeAxis,VariableAxis,artype = AsDataFrame())
outdims = OutDims()
o = mapCube(docor,dmem,indims=indims,outdims=outdims)

@test all(isapprox.(o.data,[cor(dmem.data[i,j,:,1],dmem.data[i,j,:,2]) for i=1:4, j=1:4]))

function annMean(xout,xin)
    #xin is now a DataFrame where time is added as the third column
    #We derive the year and add it to the dataframe
    xin[!,:year] = year.(xin[!,:Time])
    #Now we do the annual aggregation, note that we have to exclude the time column, because we can't aggregate here
    x2 = aggregate(xin[!,[1,2,4]],:year,a->mean(skipmissing(a)))
    #We copy the results to our output
    xout[:,1] = x2[!,2]
    xout[:,2] = x2[!,3]
end
indims = InDims("time","var",artype=AsDataFrame(true))
outdims = OutDims(RangeAxis("Year",2002:2008),"var")

o = mapCube(annMean,dmem,indims=indims,outdims=outdims)

@test all(isapprox.(o.data[1,1,:,:],mean(dmem.data[:,:,1:46,1],dims=3)))
end


text:
const TimeNS = typeof(time_ns())
const TIME_NOT_AVAILABLE = typemax(TimeNS)
const TIME_MISSING = TIME_NOT_AVAILABLE - 1
# TODO: use a wrapper type instead of the manual sentinel values

@inline function interpret_time(x::TimeNS)
    if x === TIME_NOT_AVAILABLE
        return nothing
    elseif x === TIME_MISSING
        return missing
    else
        return x
    end
end

const NOTAG = Symbol("")

"""
A mutable object whose identity (pointer) identifies a location.
"""
mutable struct LocationProxy
    line::LineNumberNode
    _module::Module
    tag::Symbol
end

macro location()
    QuoteNode(LocationProxy(__source__, __module__, NOTAG))
end

# Used in test:
_check_location1() = @location
_check_location2() = @location
_check_location3() = @location

event_id_counters(n = Threads.nthreads()) = RecordArrays.unsafe_zeros(UInt, n)

const EVENT_IDS_REF = Ref{typeof(event_id_counters())}()

function init_event_ids!()
    EVENT_IDS_REF[] = event_id_counters()
end

init_event_ids!()
atexit() do
    EVENT_IDS_REF[] = event_id_counters(0)
end

@inline event_id_ref(tid) = view(EVENT_IDS_REF[], tid)

struct RawEventID
    threadid::typeof(Threads.threadid())
    localid::UInt
end

function event_id()
    tid = Threads.threadid()
    ref = event_id_ref(tid)
    lid = ref[] += 1
    return RawEventID(tid, lid)
end

struct IntervalRecord
    location::UInt
    event::RawEventID
    approx_taskid::UInt
    start::TimeNS
    stop::TimeNS
end
# TODO: move `tag` here (requires GC object support in RecordArrays)

function interval_record_stores(blocksize = 2^8)
    BlockLinkedList(RecordVector{IntervalRecord}, blocksize)
end

const INTERVAL_RECORD_STORES = typeof(interval_record_stores())[]

function init_interval_record_stores!()
    empty!(INTERVAL_RECORD_STORES)
    append!(
        INTERVAL_RECORD_STORES,
        (interval_record_stores() for _ in 1:Threads.nthreads()),
    )
end

init_interval_record_stores!()
atexit() do
    empty!(INTERVAL_RECORD_STORES)
end

function interval_begin(location, event)
    handle = alloclast!(INTERVAL_RECORD_STORES[Threads.threadid()])
    handle.location[] = UInt(pointer_from_objref(location))
    handle.event[] = event
    handle.approx_taskid[] = UInt(pointer_from_objref(current_task()))
    handle.start[] = time_ns()
    handle.stop[] = TIME_MISSING
    return handle
end
# `approx_taskid` is only an approximation of task identity since the same
# memory location can be reused. Note that `objectid` has the same problem (and
# also a few ns slower).

function interval_end!(handle)
    handle.stop[] = time_ns()
end

function _record_point(location, event)
    handle = interval_begin(location, event)
    handle.stop[] = TIME_NOT_AVAILABLE
    return handle
end

const RawHandle = let
    location = @location
    event = event_id()
    handle = interval_begin(location, event)
    typeof(handle)
end

"""
    @recordinterval tag code
    @recordinterval code

Record the start and stop times for executing `code`.

`tag` must be a literal symbol; i.e., `@recordinterval :my_tag ...` instead of
`@recordinterval my_tag ...`.

# Examples

```jldoctest
julia> using EventTrackerBase

julia> @recordinterval begin
           sleep(0.01)
           1
       end
1
```

The timing can now be printed with `using EventTracker` then
`EventTracker.summary()`.
"""
macro recordinterval(tag::QuoteNode, code)
    tag.value isa Symbol || error("`tag` must be a symbol; got: $(tag.value)")
    return recordinterval_impl(__source__, __module__, tag.value, code)
end

macro recordinterval(code)
    return recordinterval_impl(__source__, __module__, NOTAG, code)
end

macro recordinterval(tag::Symbol, code)
    error("use `@recordinterval :$tag ...`")
end

"""
    @recordpoint tag
    @recordpoint

Record the times this part of code is executed.
"""
macro recordpoint(tag::QuoteNode)
    tag.value isa Symbol || error("`tag` must be a symbol; got: $(tag.value)")
    return recordpoint_impl(__source__, __module__, tag.value)
end

macro recordpoint()
    return recordpoint_impl(__source__, __module__, NOTAG)
end

macro recordpoint(tag::Symbol)
    error("use `@recordpoint :$tag`")
end

# Rooting is probably unnecessary and is done by the compiler automatically
# already. But maybe it's also nice to iterate over instrumented locations
# anyway?
# (TODO: Iteration over instrumented location won't work for precompiled
# locations. It's possible to support it by injecting `__init__` function that
# populates `LOCATIONS_ROOT` upon import.)
const LOCATIONS_ROOT = LocationProxy[]

function recordinterval_impl(__source__, __module__, tag::Symbol, code)
    # Injecting `location` objection in the AST so that the compiler will root
    # it (I guess?).
    location = LocationProxy(__source__, __module__, tag)
    # So, this is probably not needed (but see above):
    push!(LOCATIONS_ROOT, location)
    return quote
        local location = $(QuoteNode(location))
        local event = $event_id()
        local handle = $interval_begin(location, event)
        local ans = $(esc(Expr(:block, __source__, code)))
        $interval_end!(handle)
        ans
    end
end

function recordpoint_impl(__source__, __module__, tag::Symbol)
    location = LocationProxy(__source__, __module__, tag)
    push!(LOCATIONS_ROOT, location)
    return quote
        local location = $(QuoteNode(location))
        local event = $event_id()
        local handle = $_record_point(location, event)
        nothing
    end
end

"""
    EventTracker.clear()

Clear all records. The caller must ensure that no recordings are happening
concurrently.
"""
function clear()
    init_interval_record_stores!()
    return
end


text:
### Process examples
using Pkg
Pkg.add(Pkg.PackageSpec(; url="https://github.com/JuliaGaussianProcesses/JuliaGPsDocs.jl")) # While the package is unregistered, it's a workaround

using JuliaGPsDocs

using KernelFunctions
using PDMats, Kronecker  # we have to load all optional packages to generate the full API documentation

JuliaGPsDocs.generate_examples(KernelFunctions)

### Build documentation
using Documenter

# Doctest setup
DocMeta.setdocmeta!(
    KernelFunctions,
    :DocTestSetup,
    quote
        using KernelFunctions
    end;  # we have to load all packages used (implicitly) within jldoctest blocks in the API docstrings
    recursive=true,
)

makedocs(;
    sitename="KernelFunctions.jl",
    format=Documenter.HTML(),
    modules=[KernelFunctions],
    pages=[
        "Home" => "index.md",
        "userguide.md",
        "kernels.md",
        "transform.md",
        "metrics.md",
        "create_kernel.md",
        "API" => "api.md",
        "Design" => "design.md",
        "Examples" => JuliaGPsDocs.find_generated_examples(KernelFunctions),
    ],
    strict=true,
    checkdocs=:exports,
    doctestfilters=JuliaGPsDocs.DOCTEST_FILTERS,
)

deploydocs(;
    repo="github.com/JuliaGaussianProcesses/KernelFunctions.jl.git", push_preview=true
)


text:
function IllinoisSolver()
    return IllinoisSolver(())
end

function IllinoisSolver(arg0::jdouble)
    return IllinoisSolver((jdouble,), arg0)
end

function IllinoisSolver(arg0::jdouble, arg1::jdouble)
    return IllinoisSolver((jdouble, jdouble), arg0, arg1)
end

function IllinoisSolver(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return IllinoisSolver((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end


text:
# Activating Julia environment
using Pkg
Pkg.activate(".")

# Packages
using JuMP, MathOptInterface, Gurobi, Random, Distributions, LinearAlgebra, CSV, DataFrames
MOI = MathOptInterface

# CRUNCHING THE DATA (just like in lecture)
# NODES (I = International supplier; R = Regional supplier; L = Local market (both supply and deliver); D = delivery point)
N = []
N_I = []     # set of international suppliers
N_R = []     # set of regional suppliers
N_L = []     # set of local markets
N_D = []     # delivery points
dem = Dict() # set of demands
file = CSV.File("syria_nodes.csv")
for row in file
    push!(N, row.Name)
    if !ismissing(row.Demand)
        dem[row.Name] = row.Demand
    end
    if row.Type == "I"
        push!(N_I, row.Name)
    elseif row.Type == "R"
        push!(N_R, row.Name)
    elseif row.Type == "L"
        push!(N_L, row.Name) # Note: local markets supply and deliver goods at the given cost.
    elseif row.Type == "D"
        push!(N_D, row.Name)
    else
        throw(ErrorException("rowType $(row.Type) not supported."))
    end
end

# EDGES
hc = DataFrame(CSV.File("syria_edges.csv"))

# FOOD NUTRITION AND INTERNATIONAL COSTS
fooddata = DataFrame(CSV.File("syria_foodnutrition.csv"))
intfoodcosts = select(fooddata, [:Food, :InternationalPrice])
commodities = sort(Array(intfoodcosts.Food)) # Commodities
select!(fooddata, Not([14,15]))
fooddata = Dict(fooddata.Food .=> eachrow(fooddata))

# FOOD COST ($/metric ton for regional suppliers)
pc = DataFrame(CSV.File("syria_foodcost.csv"))
pc.Food = convert.(String63, pc.Food)
for int_supply_node in N_I # adding international prices to pc for easier processing
    for row in eachrow(intfoodcosts)
        append!(pc, DataFrame(:A => N_I, :Food => row.Food, :Price => row.InternationalPrice))
    end
end
pc = unique(pc)
international_items = DataFrame([r for r in eachrow(pc) if r.A in N_I])
regional_items =  DataFrame([r for r in eachrow(pc) if r.A in N_R])

# FOOD REQUIREMENTS (avg. per person per day)
foodreqs = DataFrame(CSV.File("syria_foodreq.csv"))
select!(foodreqs, Not(:Type))
nutrients = String.(propertynames(foodreqs))
foodreqs = Dict(string(pptname) => foodreqs[1, pptname] for pptname in propertynames(foodreqs))

function WFP_model(slack::Union{Nothing, Real} = 1, pof::Real = 0.5)
    # FINALLY CREATING THE MODEL
    m = Model(Gurobi.Optimizer)
    m[:pof] = pof
    normalquant = quantile(Normal(0,1), pof)

    # Procurement and delivery
    procurement_links = unique([row.A => row.Food for row in eachrow(pc)])       # all places where we can procure food
    @variable(m, procurement[A = N, Food = commodities; (A => Food) in procurement_links] >= 0) # procurement in tons
    @variable(m, delivery[N_D, commodities] >= 0)                                               # delivery in tons
    # Total procurement cost
    @variable(m, procurement_cost >= 0)
    if normalquant == 0
        @constraint(m, procurement_cost >= sum(r[:Price] * procurement[r.A, r.Food] for r in eachrow(pc)))
    else
        @variable(m, proccert[1:size(pc, 1)])
        for i=1:size(pc,1)
            r = eachrow(pc)[i]
            if r.A in N_I
                @constraint(m, proccert[i] == 0.05*r.Price*procurement[r.A, r.Food])
            else
                @constraint(m, proccert[i] == 0.30*r.Price*procurement[r.A, r.Food] +
                                    0.05*sum(cr.Price*procurement[cr.A, cr.Food] for cr in eachrow(pc) if cr.A == r.A))
            end
        end
        @constraint(m, [1/normalquant*(procurement_cost - sum(r[:Price] * procurement[r.A, r.Food] for r in eachrow(pc))),
                        proccert...] in MOI.SecondOrderCone(length(proccert) + 1))
        # Equivalent constraint in understandable format
        # @constraint(m, procurement_cost >= normalquant * norm(proccert, 2) +
        #         sum(r[:Price] * procurement[r.A, r.Food] for r in eachrow(pc)))
    end

    # Transportation
    transportation_links = unique([row.A => row.B for row in eachrow(hc)])
    @variable(m, transportation[A = N, B = N; (A => B) in transportation_links] >= 0)
    @variable(m, F[A = N, B = N, W = commodities; (A => B) in transportation_links] >= 0) # connections F
    for r in eachrow(hc) # Linking transportation cost to total food hauled
        @constraint(m, transportation[r.A, r.B] == sum(F[r.A, r.B, commodity] for commodity in commodities))
    end
    # Total transportation cost
    @variable(m, transportation_cost >= 0)
    @constraint(m, transportation_cost >= sum(r.tCost * transportation[r.A, r.B] for r in eachrow(hc)))

    # Flow constraints
    for node in N
        valid_sources = [link.first for link in transportation_links if link.second == node]
        valid_sinks = [link.second for link in transportation_links if link.first == node]
        for commodity in commodities
            if (node =>  commodity) in procurement_links
                @constraint(m, procurement[node, commodity] + sum(F[source, node, commodity] for source in valid_sources) ==
                                    sum(F[node, sink, commodity] for sink in valid_sinks))
            elseif node in N_D
                @assert length(valid_sinks) == 0
                @constraint(m, delivery[node, commodity] == sum(F[source, node, commodity] for source in valid_sources))
            else
                @constraint(m, sum(F[source, node, commodity] for source in valid_sources) ==
                                    sum(F[node, sink, commodity] for sink in valid_sinks))
            end
        end
    end

    # Other helpers
    @variable(m, international_p_costs >= 0)
    @constraint(m, international_p_costs == sum(procurement[r.A, r.Food] * r.Price for r in eachrow(international_items)))

    # Serving demand
    @variable(m, ration_pp[commodities] >= 0) # Rations (kg/person) of commodities
    @variable(m, nutrients_pp[nutrients] >= 0) # Total nutrients per person

    # Making sure the rations are good nutritionally
    @variable(m, nutrient_slack >= 0)
    if slack isa Real
        @constraint(m, nutrient_slack == slack)
    end
    for nutrient in nutrients # Note the factor of 10 for conversion of 100g to kg (since rations are in kg/pp)
        @constraint(m, nutrients_pp[nutrient] == 10 * sum(ration_pp[commodity] * fooddata[commodity][nutrient] for commodity in commodities))
        @constraint(m, nutrients_pp[nutrient] >= nutrient_slack * foodreqs[nutrient])
    end
    for node in N_D
        for commodity in commodities
            @constraint(m, 1000*delivery[node, commodity] >= dem[node] * ration_pp[commodity])
        end
    end

    # Diet constraints
    # Achieving a greater than 4:1 ratio by mass of carbohydrates to protein, and a greater than 4:1 ratio by mass of carbohydrates to fats
    # for the same total energy intake. The energy stored in a gram of carbohydrate, protein and fat are 4kcal/g, 4kcal/g and 9kcal/g respectively.
    @variable(m, carbs_pp >= 0)
    @constraint(m, 4*carbs_pp == nutrients_pp["Energy(kcal)"] - 4*nutrients_pp["Protein(g)"] - 9*nutrients_pp["Fat(g)"])
    @constraint(m, carbs_pp >= 4*nutrients_pp["Protein(g)"])
    @constraint(m, carbs_pp >= 4*nutrients_pp["Fat(g)"])

    # Setting objectives
    @objective(m, Min, procurement_cost + transportation_cost)
    return m
end

function report_results(m)
    return DataFrame("PP" => m[:pof],
                     "TotalCost" => round(value(m[:procurement_cost] + m[:transportation_cost]), sigdigits = 5),
                     "Proc/Total" => round(value(m[:procurement_cost]) / value(m[:procurement_cost] + m[:transportation_cost]), sigdigits=3),
                     "Trans/Total" => round(value(m[:transportation_cost]) / value(m[:procurement_cost] + m[:transportation_cost]), sigdigits=3),
                     "Intl/TotalProc" => round(value(m[:international_p_costs]) / value(m[:procurement_cost]), sigdigits=3),
                     "Slack" =>  round(value(m[:nutrient_slack]), sigdigits = 3),
                     "Cost/Person" => round(value(m[:procurement_cost] + m[:transportation_cost]) /
                                        (value(m[:nutrient_slack]) * sum(values(dem))), sigdigits=5),
                     "NumActiveProc" => sum(values(value.(m[:procurement])) .> 1e-3),
                     "NumActiveTrans" => sum(values(value.(m[:transportation])) .> 1e-3))
end

# Solving nominal problem
m = WFP_model(1)
optimize!(m);
report_results(m)

# Maximizing nutrition (optimal solution just scales the nutrition)
nm = WFP_model(nothing, 0.5)
# set_optimizer_attribute_(nm, "OptimalityTol", 1e-8)
@constraint(nm, nm[:procurement_cost] + nm[:transportation_cost] <= 6000)
@objective(nm, Max, nm[:nutrient_slack])
optimize!(nm);
report_results(nm)

# Robust solutions with a budget limit
pofs = [75, 85, 90, 92, 95, 96, 97, 98, 99, 99.5]./100
robustdata_lim = report_results(nm)
for pof in pofs
    m_lim = WFP_model(nothing, pof)
    @constraint(m_lim, m_lim[:procurement_cost] + m_lim[:transportation_cost] <= 6000)
    @objective(m_lim, Max, m_lim[:nutrient_slack])
    optimize!(m_lim);
    df = report_results(m_lim)
    append!(robustdata_lim, df)
end

# Robust solutions without a budget limit
m_unlim = WFP_model(1)
optimize!(m_unlim);
robustdata_unlim = report_results(m_unlim)
for pof in pofs
    m_unlim = WFP_model(1, pof)
    optimize!(m_unlim);
    df = report_results(m_unlim)
    append!(robustdata_unlim, df)
end


text:
#=
Copyright 2020 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

__precompile__(true)
module openBF

    export Vessel, Heart, Blood, runSimulation
    using YAML
    using ArgParse
    using StaticArrays
    using DelimitedFiles
    using LinearAlgebra
    using Printf

    """
    Heart type
    """
    struct Heart
        inlet_type :: String
        cardiac_T :: Float64
        input_data :: Array{Float64,2}
        inlet_number :: Int
    end

    """
    Blood type
    """
    struct Blood
        mu  :: Float64
        rho :: Float64
        rho_inv :: Float64
    end

    """
    Vessel type
    """
    mutable struct Vessel
        label :: String

        #Topology
        ID :: Int
        sn :: Int
        tn :: Int

        #Inlet
        inlet :: Bool
        heart :: Heart

        #Numerical constants
        M       :: Int
        dx      :: Float64
        invDx   :: Float64
        halfDx  :: Float64

        #Physical constants
        beta            :: Array{Float64,1}
        gamma           :: Array{Float64,1}
        s_15_gamma      :: Array{Float64,1}
        gamma_ghost     :: Array{Float64,1}
        A0              :: Array{Float64,1}
        s_A0            :: Array{Float64,1}
        inv_A0          :: Array{Float64,1}
        s_inv_A0        :: Array{Float64,1}
        Pext            :: Float64
        viscT           :: Float64
        wallE           :: Array{Float64,1}
        wallVa          :: Array{Float64,1}
        wallVb          :: Array{Float64,1}

        #Iterative solution
        A :: Array{Float64,1}
        Q :: Array{Float64,1}
        u :: Array{Float64,1}
        c :: Array{Float64,1}
        P :: Array{Float64,1}

        A_t :: Array{Float64,2}
        Q_t :: Array{Float64,2}
        u_t :: Array{Float64,2}
        c_t :: Array{Float64,2}
        P_t :: Array{Float64,2}

        A_l :: Array{Float64,2}
        Q_l :: Array{Float64,2}
        u_l :: Array{Float64,2}
        c_l :: Array{Float64,2}
        P_l :: Array{Float64,2}

        #Riemann invariants
        W1M0 :: Float64
        W2M0 :: Float64

        #Ghost cells
        U00A :: Float64
        U00Q :: Float64
        U01A :: Float64
        U01Q :: Float64

        UM1A :: Float64
        UM1Q :: Float64
        UM2A :: Float64
        UM2Q :: Float64

        #Result file names
        last_P_name :: String
        last_Q_name :: String
        last_A_name :: String
        last_c_name :: String
        last_u_name :: String

        out_P_name :: String
        out_Q_name :: String
        out_A_name :: String
        out_c_name :: String
        out_u_name :: String

        #Saving nodes
        node2 :: Int
        node3 :: Int
        node4 :: Int

        #Peripheral boundary condition
        Rt :: Float64
        R1 :: Float64
        R2 :: Float64
        Cc :: Float64
        Pc :: Float64

        #Slope
        slope :: Array{Float64,1}

        #MUSCLArrays
        flux :: Array{Float64,2}
        uStar :: Array{Float64,2}

        vA :: Array{Float64,1}
        vQ :: Array{Float64,1}

        dU :: Array{Float64,2}

        slopesA :: Array{Float64,1}
        slopesQ :: Array{Float64,1}

        Al :: Array{Float64,1}
        Ar :: Array{Float64,1}

        Ql :: Array{Float64,1}
        Qr :: Array{Float64,1}

        Fl :: Array{Float64,2}
        Fr :: Array{Float64,2}

        #Outlet type
        outlet :: String
    end

    include("initialise.jl")

    include("boundary_conditions.jl")
    include("solver.jl")

    include("junctions.jl")
    include("conjunctions.jl")
    include("bifurcations.jl")
    include("anastomosis.jl")

    include("IOutils.jl")
    include("check_convergence.jl")

    include("program.jl")
end


text:
"""
barychebinit(n)

Get Chebyshev nodes of first kind and corresponding barycentric Lagrange interpolation weights.
Reference: Berrut, J.P. and Trefethen, L.N., 2004. Barycentric lagrange interpolation. SIAM review, 46(3), pp.501-517.

# Arguments
- `n`: order of the Chebyshev interpolation

# Returns
- Chebyshev nodes
- Barycentric Lagrange interpolation weights
"""
function barychebinit(n)
    x = zeros(Float64, n)
    w = similar(x)
    for i in 1:n
        c = (2i - 1)π / (2n)
        x[n - i + 1] = cos(c)
        w[n - i + 1] = (-1)^(i - 1) * sin(c)
    end
    return x, w
end

"""
function barycheb(n, x, f, wc, xc)

Barycentric Lagrange interpolation at Chebyshev nodes
Reference: Berrut, J.P. and Trefethen, L.N., 2004. Barycentric lagrange interpolation. SIAM review, 46(3), pp.501-517.

# Arguments
- `n`: order of the Chebyshev interpolation
- `x`: coordinate to interpolate
- `f`: array of size n, function at the Chebyshev nodes
- `wc`: array of size n, Barycentric Lagrange interpolation weights
- `xc`: array of size n, coordinates of Chebyshev nodes

# Returns
- Interpolation result
"""
function barycheb(n, x, f, wc, xc)
    for j in 1:n
        if x == xc[j]
            return f[j]
        end
    end

    num, den = 0.0, 0.0
    for j in 1:n
        q = wc[j] / (x - xc[j])
        num += q * f[j]
        den += q
    end
    return num / den
end

function barycheb2(n, x, f, wc, xc)
    for j in 1:n
        if x == xc[j]
            return f[j, :]
        end
    end

    den = 0.0
    num = zeros(eltype(f), size(f)[2])
    for j in 1:n
        q = wc[j] / (x - xc[j])
        num += q .* f[j, :]
        den += q
    end
    return num ./ den
end