text:
using LinearAlgebra
using DataFrames
using CSV

#This file contains various renditions of  which perform numerical
#quadrature given integration endpoints and a function handle
#More importantly, it provides a composite type which stores nodes and weight.

"""
Data-type which stores quadrature nodes and weights
"""

mutable struct nodesWeights{T1<:Array{Float64}, T2<:Array{Float64}}
    nodes::T1
    weights::T2
end

"""
Apply affine transformation to nodes to integration over range r (linear change of variables)
"""
function affineTransformNodes(nw::nodesWeights{O}, r::Array{<:Real}) where O<:Union{Array{Float64, 2}, Array{Float64, 1}}
    b = r[2]; a = r[1]
    nw.nodes = (b-a)/2 .* nw.nodes .+ (b + a)/2
end

"""
Uses the Golub-Welsch eigenvalue method to compute
Gauss-Legendre quadrature nodes and weights for the domain [-1, 1]

INPUTS: n, the number of desired nodes
OUTPUTS: x and w, the nodes and weights
"""
function gausslegpts(n)
    x = zeros(n-1)
    for i = 1:n-1
        x[i] = i/sqrt(4*i^2-1)
    end
    x = SymTridiagonal(zeros(n), x) #Jacobi matrix
    d, V = eigen(x)
    w = 2*V[1,:].^2 #factor of 2 comes from the fact that integral of measure is 2
    return d, w
end

#precomute various integration nodes/weights
dm0, wm0 = gausslegpts(24)
dm, wm = gausslegpts(12)
dm2, wm2 = gausslegpts(50)
dm3, wm3 = gausslegpts(100)

function getGaussQuadraturedata()
    nodesWeights(dm, wm)
end

"""
INPUTS:
f is a function handle, a and b are integration endpoints (a<b), and n
is the number of quadrature nodes (defaults to 10). Note that f can be
multi-output, in which case int1D integrates each entry individually. This way,
the function lends itself to being composed with Gauss-Turan integration.
"""
function int1D(f, arr, num="1")
    if num=="1"
        x, w= dm, wm
    elseif num=="2"
        x, w= dm2, wm2
    elseif num=="0"
        x, w = dm0, wm0
    else
        x, w = dm3, wm3
    end
    sample = f(arr[1])
    res = zeros(size(sample))
    for i = 1:length(x)
        expr = (arr[2]-arr[1])/2 * x[i] .+ (arr[1]+arr[2])/2
        #println(res)
        #println(f(expr) .* w[i])
        res = res .+ f(expr) .* w[i]
    end
    res = (arr[2]-arr[1])/2 .* res
end

function int1D_print(f, arr, num="0")
    if num=="1"
        x, w= dm, wm
    elseif num=="2"
        x, w= dm2, wm2
    elseif num=="0"
        x, w = dm0, wm0
    else
        x, w = dm3, wm3
    end
    sample = f(arr[1])
    res = zeros(size(sample))
    for i = 1:length(x)
        expr = (arr[2]-arr[1])/2 * x[i] .+ (arr[1]+arr[2])/2
        println(res)
        println(f(expr) .* w[i])
        res = res .+ f(expr) .* w[i]
    end
    res = (arr[2]-arr[1])/2 .* res
end

"""
Integrates function over 2D domain using Gauss-Legendre Quadrature
arr[:, 1] are lower bounds, while arr[:, 2] are upper bounds
"""
function int2D(f, arr)
    d = size(arr, 1) #number of dimensions
    int = 0.0
    x, w = dm, wm
    avg1 = (arr[1, 2] - arr[1, 1])/2
    avg2 = (arr[2, 2] - arr[2, 1])/2
    for i = 1:length(x)
        for j = 1:length(x)
           expr1 = avg1*x[i] + (arr[1, 2]+arr[1, 1])/2
           expr2 = avg2*x[j] + (arr[2, 2]+arr[2, 1])/2
           int = int + w[i]*w[j]*f(expr1, expr2)
        end
    end
    int = int*avg1*avg2
end

data50 = DataFrame(CSV.File("nodes_weights_50.csv", header=0))
nodes50 = convert(Array, data50[:,1]) #integration nodes for Gauss-Turan Quadrature
weights50 = convert(Matrix, data50[:, 2:end]) #integration weights

data20 = DataFrame(CSV.File("nodes_weights_20.csv", header=0))
nodes20  = convert(Array, data20[:,1])
weights20 = convert(Array, data20[:,2:end])

data12 = DataFrame(CSV.File("nodes_weights_12.csv", header=0))
nodes12  = convert(Array, data12[:,1])
weights12 = convert(Array, data12[:,2:end])

function getTuranQuadratureData()
    return nodesWeights(nodes12, weights12)
end

"""
Gauss-Turan integration with 2 derivatives
INPUTS:
arr contains enpoints of integration interval
"""
function Gauss_Turan(F, arr, nodes = nodes12, weights=weights12)
    a = arr[1]
    b = arr[2]
    nodes = (b-a)/2 .* nodes .+ (b+a)/2
    T = zeros(length(nodes), 3)
    time = @elapsed begin
    for i = 1:length(nodes)
        (T[i, 1], T[i, 2], T[i, 3]) = F(nodes[i])
    end
    end
    println(time)
    fn = T[:, 1]
    dfn = T[:, 2]
    d2fn = T[:, 3]
    return (fn'*weights[:, 1] + dfn'*weights[:, 2] + d2fn'*weights[:, 3])*(b-a)/2
end

"""
n-dimensional Gauss-Legendre integration
"""
function intnD(f, arr, n)

end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
export Airspace, box_spherical_polygon, isinside

struct Airspace
    name::String
    polygon::Vector{Point{Float64}}
    bounding_box::Vector{Point{Float64}}
    function Airspace(name::String, polygon::Vector{Point{Float64}})
        box = box_spherical_polygon(polygon)
        new(name, polygon, box)
    end
end

function box_spherical_polygon(polygon::Vector{Point{Float64}})
    # box = Array{Float64, 2}(undef, 2, 2)
    lats_max = Array{Float64,1}()
    lats_min = Array{Float64,1}()
    previous_point = Point(NaN, NaN)
    for point in polygon
        # println(point)
        if isnan(previous_point.λ)
            previous_point = point
        else
            append!(lats_max, maximum(latitude_options(previous_point, point)))
            append!(lats_min, minimum(latitude_options(previous_point, point)))
            previous_point = point
        end
    end
    lat_max = maximum(lats_max)
    lat_min = minimum(lats_min)

    # Limiting longitudes
    lon_generator = (x.λ for x in polygon)
    lon_max = maximum(lon_generator)
    lon_min = minimum(lon_generator)

    box = Vector{Point}()
    box = vcat(box, Point(lat_min, lon_min))
    box = vcat(box, Point(lat_min, lon_max))
    box = vcat(box, Point(lat_max, lon_max))
    box = vcat(box, Point(lat_max, lon_min))
    box = vcat(box, Point(lat_min, lon_min))
    return box
end

function latitude_options(pos₁::Point, pos₂::Point)
    brg = bearing(pos₁, pos₂)
    closest_point = closest_point_to_pole(pos₁, brg)
    if distance(pos₁, pos₂) > distance(pos₁, closest_point)
        return pos₁.ϕ, pos₂.ϕ, closest_point.ϕ
    else
        return pos₁.ϕ, pos₂.ϕ
    end
end

function isinside(pos₁::Point, polygon::Vector{Point{Float64}})
    south_pole = Point(-90.0, 0.0)
    poly₁ = polygon[1]
    count_intersections = 0
    for poly₂ in polygon[2:end]
        isnan(intersection_point(pos₁, south_pole, poly₁, poly₂).ϕ) ? false : count_intersections += 1
        poly₁ = poly₂
    end
    return count_intersections%2 != 0
end

ext: jl


text:
module SteadyStateDiffEq

using Reexport
@reexport using DiffEqBase

using NLsolve, DiffEqCallbacks
using LinearAlgebra

include("algorithms.jl")
include("solve.jl")

export SSRootfind, DynamicSS

end # module


text:
module AMQPClient

import Base: write, read, read!, close, convert, show, isopen

using Sockets

const DEBUG = false

# 0.7: use builtin logging by enabling following statement
# using Logging; Logging.global_logger(Logging.ConsoleLogger(stderr, Logging.Debug))
macro debug(s)
    esc(:(DEBUG && Base.@debug($s)))
end

# Client property info that gets sent to the server on connection startup
const CLIENT_IDENTIFICATION = Dict{String,Any}(
    "product" => "Julia AMQPClient",
    "product_version" => string(VERSION),
    "capabilities" => Dict{String,Any}()
)

include("types.jl")
include("spec.jl")
include("message.jl")
include("auth.jl")
include("protocol.jl")
include("convert.jl")
include("show.jl")

export connection, channel, CloseReason
export exchange_declare, exchange_delete, exchange_bind, exchange_unbind, default_exchange_name
export queue_declare, queue_bind, queue_unbind, queue_purge, queue_delete
export tx_select, tx_commit, tx_rollback
export basic_qos, basic_consume, basic_cancel, basic_publish, basic_get, basic_ack, basic_reject, basic_recover
export confirm_select
export EXCHANGE_TYPE_DIRECT, EXCHANGE_TYPE_FANOUT, EXCHANGE_TYPE_TOPIC, EXCHANGE_TYPE_HEADERS
export read, read!, close, convert, show
export Message, set_properties, PERSISTENT, NON_PERSISTENT

end # module


text:
# A Modeling Language Package for Solving Data Envelopment Analysis Problems
module FrontierEfficiencyAnalysis

# An algebraic modeling language for Julia
import JuMP
# Solver-independent functions and low-level interface for Mathematical Programming
import MathProgBase
# An example solver
import Gurobi

# function
export solveDEA

# solving dea problems by a smart algorithm
# it is computationally efficient and can solve large scale DEA problems via any solvers (even with LP size limitation)
include("solveDEA.jl")

end


text:
function initial_position(i,j,ls,λ,shift,posx,posy)
    xs=steady_position(i,ls)
    ys=steady_position(j,ls)
    dx=xs-posx
    dy=ys-posy
    dr=sqrt(dx^2+dy^2)
    ds=(-1/λ^2)*exp(-(dr/λ)^2/2)*shift
    xs+dx*ds,ys+dy*ds
end

function update_force!(fx::Array{Float64,2},xc::Array{Float64,2},ks)
    nsx,nsy=size(xc)
    @inbounds for j ∈ 2:nsy-1
        for i ∈ 2:nsx-1
            fx[i,j] = -ks*(4xc[i,j]
            -xc[i,j-1]
            -xc[i-1,j]
            -xc[i+1,j]
            -xc[i,j+1])
        end
    end
end

function update_force!(fx::AbstractArray{Float64,2},xc::AbstractArray{Float64,2},ks)
    ns=size(xc,1)
    r0=1:ns-2
    r1=2:ns-1
    r2=3:ns

    @views @. fx[r1,r1] = -ks*(4xc[r1,r1]
                            -xc[r1,r0]    # Down
                            -xc[r0,r1]    # Left
                            -xc[r2,r1]    # Right
                            -xc[r1,r2])   # Up
end

#V is the type of the returned arrays (type are first class citizens !)
function initialize_arrays(sp,ip,ap,V=Array{Float64,2})
    ls,ms,ks,ns=getvalues(sp)
    λ,shift,pos=getvalues(ip)
    δt,nδt,nδtperframe=getvalues(ap)

    #2D Array comprehensions
    #  -> cxs,xys,cxc,cyc are 2D arrays of size (ns,ns)
    cxs=[steady_position(i,ls) for i ∈ 1:ns, j ∈ 1:ns]
    cys=[steady_position(j,ls) for i ∈ 1:ns, j ∈ 1:ns]
    cxc=[initial_position(i,j,ls,λ,shift,pos,pos)[1] for i ∈ 1:ns, j ∈ 1:ns]
    cyc=[initial_position(i,j,ls,λ,shift,pos,pos)[2] for i ∈ 1:ns, j ∈ 1:ns]

    #Broadcast of V Ctor equivalent to xs=V(cxs);ys=V(cys);xc=V(cxc);yc=V(cyc)
    xs,ys,xc,yc=V.((cxs,cys,cxc,cyc))
    #Assume zero initial velocities : xp(revious)=xc(urrent)
    fx,xt,xp=zero(xc),zero(xc),copy(xc)
    fy,yt,yp=zero(yc),zero(yc),copy(yc)

    xs,xc,xp,xt,fx,ys,yc,yp,yt,fy
end

function display_perf(td,xc,nδt)
    nflops=20 # Manual count floating point ops
    nRW=12    # Manual count Array Reads and Writes
    T=eltype(xc) # i.e. Float64
    floatsize_inbytes=sizeof(T) #Float64 -> 8 Bytes
    nbytes=nRW*floatsize_inbytes # Float64 -> 96 Bytes

    mupns=length(xc)*nδt/(td*1.e9) # mass pos updates per ns
    println("$(round((mupns*nflops),sigdigits=3))\t GFLOPS")
    println("$(round((mupns*nbytes),sigdigits=3))\t GB/s")
end

do_synchronize(V) = nothing #no need for sync for CPU Arrays

function animate_spring2D(sp,ip,ap,V=Array{Float64,2})
    δt,nδt,nδtperframe=getvalues(ap)
    xs,xc,xp,xt,fx,ys,yc,yp,yt,fy=initialize_arrays(sp,ip,ap,V)

    dc=zero(xc)    # Displacement (type V)
    cdc=Array(xc)  # Displacement (type Array{Float64,2})
    @. dc = sqrt((xc - xs)^2+(yc - ys)^2) #displacement norm
    mdc=maximum(dc)

    xaxis=Array(xs[:,1]) # 1D Array
    yaxis=Array(ys[1,:]) # 1D Array
    nf=nδt÷nδtperframe # We assume a null remainder

    t=0.0 # Simulation time
    tdynamic=0.0 # Computing time for the dynamic
    anim = @animate for i ∈ 1:nf

        tdynamic+= @elapsed begin
            xc,xp,xt=advance_nδtpf(xc,xp,xt,fx,sp,ap)
            yc,yp,yt=advance_nδtpf(yc,yp,yt,fy,sp,ap)
            do_synchronize(V) #sync if V is a GPU Array
        end
        @. dc = sqrt((xc - xs)^2+(yc - ys)^2) #displacement norm
        @. cdc = dc  # Copy dc (type V) to cdc (type Array{Float64,2})

        t+=nδtperframe*δt
        contour(xaxis,yaxis,cdc,clims=(0,mdc/4), #level lines graph
            thickness_scaling = 1.4,size=(600,600),
            xlabel=L"x_s",ylabel=L"y_s",
            right_margin = 10Plots.PlotMeasures.mm,
            title=L"\sqrt{(x_c-x_s)^2+(y_c-y_s)^2}",aspect_ratio=:equal)
    end

    display_perf(tdynamic,xc,nδt) #Display GFLOPS and GB/s
    gif(anim,"toto.gif",fps=15) # save the animation
end


text:
using JuliaDB
using Base.Test

@testset "rechunk" begin
    t = NDSparse([1:10;], ones(10), ones(10))
    d = distribute(t, 2)
    d2 = compute(distribute(d, [2,3,4,1]))
    @test get.(map(JuliaDB.nrows, d2.domains)) == [2,3,4,1]
end

@testset "rechunk by" begin
    t = distribute(NDSparse(Columns([1,1,1,2,2,2], [1,1,2,2,3,3]), [1:6;]), [2,2,2])
    t1=rechunk(t, (1, 2=>x->x.>=2), closed=true, chunks=3)
    @test collect(t1) == selectkeys(collect(t), (1,2=>x->x>=2))
    @test collect(t1.chunks[1]) == NDSparse(Columns([1,1],[false,false]), [1,2])
  # @test collect(t1.chunks[2]) == NDSparse(Columns([1],[2]), [3])
  # @test collect(t1.chunks[3]) == NDSparse(Columns([2,2,2],[2,3,3]), [4,5,6])
end


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

JLLWrappers.@generate_main_file_header("mimalloc")
JLLWrappers.@generate_main_file("mimalloc", UUID("38cef826-3163-5c31-a333-3957420b3fb7"))
end  # module mimalloc_jll


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -36.7700820821577 <= q <= 37.6713142442229)

                     add_NL_constraint(m, :(softplus(0.9266287961121966 + 0.9630341445616355*softplus(0.9080159995418975 + 0.10769169268185319*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.8197021439185934*softplus(-0.9264374759333807 + -0.7028105280156569*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.8564198456770362*softplus(-0.7118497063916944 + 0.36586202274128254*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.7868467121279954*softplus(0.8429853377330745 + -0.6377098680254001*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) + softplus(-0.9209290943615045 + -0.04560727661330999*softplus(0.9080159995418975 + 0.10769169268185319*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.04187997574232316*softplus(-0.9264374759333807 + -0.7028105280156569*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.009897273687461627*softplus(-0.7118497063916944 + 0.36586202274128254*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.4138607422831795*softplus(0.8429853377330745 + -0.6377098680254001*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) + softplus(0.41055131347137674 + 0.5074173875637427*softplus(0.9080159995418975 + 0.10769169268185319*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.7865169798440914*softplus(-0.9264374759333807 + -0.7028105280156569*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.6322955384333278*softplus(-0.7118497063916944 + 0.36586202274128254*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.03598168048021799*softplus(0.8429853377330745 + -0.6377098680254001*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) + softplus(-0.7045693885767101 + 0.765322888364711*softplus(0.9080159995418975 + 0.10769169268185319*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.2597058016512084*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.8597063403720888*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.17950756837409587*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + 0.9087169409161477*softplus(-0.9264374759333807 + -0.7028105280156569*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.6656745963467641*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.9099166873549671*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.9082646170966706*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.29202574365958966*softplus(-0.7118497063916944 + 0.36586202274128254*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + -0.9129427463398758*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + -0.5191517354247042*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.42477634545139953*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5]))) + -0.14636407330293322*softplus(0.8429853377330745 + -0.6377098680254001*softplus(-0.02399021036854254 + -0.4943608904262171*$(x[1]) + -0.3028702711871105*$(x[2]) + 0.7730701583831108*$(x[3]) + -0.5656579689281509*$(x[4]) + 0.08835909373115314*$(x[5])) + 0.04468698484435718*softplus(-0.6908874409241781 + -0.9493883109743697*$(x[1]) + 0.1034335229123764*$(x[2]) + 0.3901179347519834*$(x[3]) + -0.16073609758033847*$(x[4]) + 0.052390589356133166*$(x[5])) + 0.3607403638969684*softplus(0.7510327333371474 + 0.48933101371416265*$(x[1]) + 0.5299590321169751*$(x[2]) + 0.05546767568522437*$(x[3]) + -0.4663274814822995*$(x[4]) + -0.6567285643500202*$(x[5])) + 0.4004266402265344*softplus(-0.27157848416188113 + 0.560003468778171*$(x[1]) + 0.2933473130688298*$(x[2]) + 0.8808820816540863*$(x[3]) + 0.6726368959629845*$(x[4]) + -0.19953521813356856*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#
# Representing data shapes and pipeline signatures.
#

import Base:
    convert,
    getindex,
    eltype,
    show

#
# Shapes of data.
#

"""
    AbstractShape

Describes the structure of column-oriented data.
"""
abstract type AbstractShape end

quoteof(shp::AbstractShape) =
    Expr(:call, nameof(typeof(shp)))

quoteof(p::Pair{Symbol,<:AbstractShape}) =
    Expr(:call, :(=>), labelquote(p.first), quoteof(p.second))

quoteof_inner(p::Pair{Symbol,<:AbstractShape}) =
    Expr(:call, :(=>), labelquote(p.first), quoteof_inner(p.second))

labelquote(lbl::Symbol) =
    quoteof(labelsyntax(lbl))

syntaxof(shp::AbstractShape) =
    nameof(typeof(shp))

syntaxof(p::Pair{Symbol,<:AbstractShape}) =
    Expr(:(=), labelsyntax(p.first), syntaxof(p.second))

labelsyntax(lbl::Symbol) =
    Base.isidentifier(lbl) ? lbl : string(lbl)

show(io::IO, shp::AbstractShape) =
    print_expr(io, quoteof(shp))

#
# Concrete shapes.
#

abstract type DataShape <: AbstractShape end

"""
    AnyShape()

Nothing is known about the data.
"""
struct AnyShape <: DataShape
end

quoteof_inner(::AnyShape) =
    :Any

syntaxof(::AnyShape) =
    :Any

eltype(::AnyShape) = Any

"""
    NoShape()

Inconsistent constraints on the data.
"""
struct NoShape <: DataShape
end

syntaxof(::NoShape) =
    :Bottom

eltype(::NoShape) = Union{}

"""
    ValueOf(::Type)

Regular Julia vector with elements of the given type.
"""

struct ValueOf <: DataShape
    ty::Type
end

ValueOf(::Type{Any}) =
    AnyShape()

ValueOf(::Type{Union{}}) =
    NoShape()

convert(::Type{AbstractShape}, ty::Type) =
    ValueOf(ty)

quoteof(shp::ValueOf) =
    Expr(:call, nameof(ValueOf), shp.ty)

quoteof_inner(shp::ValueOf) =
    shp.ty

syntaxof(shp::ValueOf) = shp.ty

eltype(shp::ValueOf) = shp.ty

"""
    TupleOf([lbls::Vector{Symbol},] cols::Vector{AbstractShape})
    TupleOf(cols::AbstractShape...)
    TupleOf(lcols::Pair{<:Union{Symbol,AbstractString},<:AbstractShape}...)

Shape of a `TupleVector`.
"""

struct TupleOf <: DataShape
    lbls::Vector{Symbol}
    cols::Vector{AbstractShape}
end

TupleOf() =
    TupleOf(Symbol[], AbstractShape[])

TupleOf(cols::Vector{<:AbstractShape}) =
    TupleOf(Symbol[], cols)

TupleOf(cols::Union{AbstractShape,Type}...) =
    TupleOf(Symbol[], AbstractShape[col for col in cols])

TupleOf(lcols::Pair{<:Union{Symbol,AbstractString},<:Union{AbstractShape,Type}}...) =
    TupleOf(Symbol[first(lcol) for lcol in lcols], AbstractShape[last(lcol) for lcol in lcols])

quoteof(shp::TupleOf) =
    if isempty(shp.lbls)
        Expr(:call, nameof(TupleOf), Any[quoteof_inner(col) for col in shp.cols]...)
    else
        Expr(:call, nameof(TupleOf), Any[quoteof_inner(lbl => col) for (lbl, col) in zip(shp.lbls, shp.cols)]...)
    end

syntaxof(shp::TupleOf) =
    if isempty(shp.lbls)
        Expr(:tuple, Any[syntaxof(col) for col in shp.cols]...)
    else
        Expr(:tuple, Any[syntaxof(lbl => col) for (lbl, col) in zip(shp.lbls, shp.cols)]...)
    end

labels(shp::TupleOf) = shp.lbls

function ordinal_label(k::Int)
    lbl = ""
    if k == 1
        lbl = 'A' * lbl
    else
        while k > 1
            lbl = ('A' + (k - 1) % 26) * lbl
            k = 1 + (k - 1) ÷ 26
        end
    end
    return Symbol('#' * lbl)
end

function label(shp::TupleOf, k::Int)
    !isempty(shp.lbls) ? shp.lbls[k] : ordinal_label(k)
end

width(shp::TupleOf) = length(shp.cols)

columns(shp::TupleOf) = shp.cols

column(shp::TupleOf, j::Int) = shp.cols[j]

column(shp::TupleOf, lbl::Symbol) =
    column(shp, findfirst(isequal(lbl), shp.lbls))

function replace_column(shp::TupleOf, j::Int, f)
    col = shp.cols[j]
    col′ = f isa AbstractShape ? f : f(col)
    cols′ = copy(shp.cols)
    cols′[j] = col′
    TupleOf(shp.lbls, cols′)
end

function eltype(shp::TupleOf)
    t = Tuple{Any[eltype(col) for col in shp.cols]...}
    if isempty(shp.lbls)
        t
    else
        NamedTuple{(shp.lbls...,),t}
    end
end

"""
    BlockOf(elts::AbstractShape, card::Cardinality=x0toN)

Shape of a `BlockVector`.
"""
struct BlockOf <: DataShape
    elts::AbstractShape
    card::Cardinality
end

BlockOf(elts) =
    BlockOf(elts, x0toN)

quoteof(shp::BlockOf) =
    if shp.card == x0toN
        Expr(:call, nameof(BlockOf), quoteof_inner(shp.elts))
    else
        Expr(:call, nameof(BlockOf), quoteof_inner(shp.elts), quoteof(shp.card))
    end

syntaxof(shp::BlockOf) =
    Expr(:call, :×, syntaxof(shp.card), syntaxof(shp.elts))

elements(shp::BlockOf) = shp.elts

function replace_elements(shp::BlockOf, f)
    elts′ = f isa AbstractShape ? f : f(shp.elts)
    BlockOf(elts′, shp.card)
end

cardinality(shp::BlockOf) = shp.card

ismandatory(shp::BlockOf) = ismandatory(shp.card)

issingular(shp::BlockOf) = issingular(shp.card)

function eltype(shp::BlockOf)
    t = eltype(shp.elts)
    if shp.card == x1to1
        t
    elseif shp.card == x0to1
        Union{t,Missing}
    else
        Vector{t}
    end
end

#
# Annotations.
#

abstract type Annotation <: AbstractShape end

syntaxof(shp::Annotation) =
    syntaxof(shp.sub)

eltype(shp::Annotation) =
    eltype(subject(shp))

deannotate(shp::AbstractShape) =
    shp

deannotate(shp::Annotation) =
    deannotate(subject(shp))

"""
    sub |> IsLabeled(::Symbol)

The shape has an attached label.
"""
struct IsLabeled <: Annotation
    sub::AbstractShape
    lbl::Symbol
end

IsLabeled(sub::Union{AbstractShape,Type}, lbl::Union{Symbol,AbstractString}) =
    IsLabeled(convert(AbstractShape, sub), Symbol(lbl))

IsLabeled(lbl::Union{Symbol,AbstractString}) =
    sub -> IsLabeled(sub, lbl)

quoteof(shp::IsLabeled) =
    Expr(:call, nameof(|>), quoteof_inner(shp.sub), Expr(:call, nameof(IsLabeled), labelquote(shp.lbl)))

subject(shp::IsLabeled) = shp.sub

replace_subject(shp::IsLabeled, f) =
    IsLabeled(f isa AbstractShape ? f : f(shp.sub), shp.lbl)

label(shp::IsLabeled) =
    shp.lbl

"""
    sub |> IsFlow

The annotated `BlockVector` holds the output flow.
"""

struct IsFlow <: Annotation
    sub::BlockOf
end

quoteof(shp::IsFlow) =
    Expr(:call, nameof(|>), quoteof_inner(shp.sub), nameof(IsFlow))

subject(shp::IsFlow) = shp.sub

replace_subject(shp::IsFlow, f) =
    IsFlow(f isa AbstractShape ? f : f(shp.sub))

elements(shp::IsFlow) =
    elements(shp.sub)

replace_elements(shp::IsFlow, f) =
    replace_subject(shp, sub -> replace_elements(sub, f))

cardinality(shp::IsFlow) =
    cardinality(shp.sub)

"""
    sub |> IsScope

The annotated `TupleVector` holds the scoping context.
"""

struct IsScope <: Annotation
    sub::TupleOf
end

quoteof(shp::IsScope) =
    Expr(:call, nameof(|>), quoteof_inner(shp.sub), nameof(IsScope))

subject(shp::IsScope) = shp.sub

replace_subject(shp::IsScope, f) =
    IsScope(f isa AbstractShape ? f : f(shp.sub))

column(shp::IsScope) =
    column(shp.sub, 1)

replace_column(shp::IsScope, f) =
    replace_subject(shp, sub -> replace_column(sub, 1, f))

context(shp::IsScope) =
    column(shp.sub, 2)

#
# Partial order on shapes.
#

"""
    fits(x::T, y::T) :: Bool

Checks if constraint `x` implies constraint `y`.
"""
function fits end

fits(c1::Cardinality, c2::Cardinality) =
    (c1 | c2) == c2

fits(::AbstractShape, ::AbstractShape) = false

fits(::DataShape, ::AnyShape) = true

fits(::NoShape, ::DataShape) = true

fits(::NoShape, ::AnyShape) = true

fits(shp1::ValueOf, shp2::ValueOf) =
    shp1.ty <: shp2.ty

fits(shp1::TupleOf, shp2::TupleOf) =
    length(shp1.cols) == length(shp2.cols) &&
    all(Bool[fits(col1, col2) for (col1, col2) in zip(shp1.cols, shp2.cols)]) &&
    (isempty(shp2.lbls) || shp1.lbls == shp2.lbls)

fits(shp1::BlockOf, shp2::BlockOf) =
    fits(shp1.elts, shp2.elts) &&
    fits(shp1.card, shp2.card)

fits(shp1::IsLabeled, shp2::IsLabeled) =
    fits(shp1.sub, shp2.sub) && shp1.lbl == shp2.lbl

fits(shp1::IsLabeled, shp2::AbstractShape) =
    fits(shp1.sub, shp2)

fits(shp1::IsFlow, shp2::IsFlow) =
    fits(shp1.sub, shp2.sub)

fits(shp1::IsFlow, shp2::AbstractShape) =
    fits(shp1.sub, shp2)

fits(shp1::IsScope, shp2::IsScope) =
    fits(shp1.sub, shp2.sub)

fits(shp1::IsScope, shp2::AbstractShape) =
    fits(shp1.sub, shp2)

#
# Guessing the shape of a vector.
#

shapeof(v::AbstractVector) =
    ValueOf(eltype(v))

shapeof(tv::TupleVector) =
    TupleOf(labels(tv), AbstractShape[shapeof(col) for col in columns(tv)])

shapeof(bv::BlockVector) =
    BlockOf(shapeof(elements(bv)), cardinality(bv))

fits(v::AbstractVector, shp::AbstractShape) =
    fits(shapeof(v), shp)

#
# Signature of a pipeline.
#

"""
    Signature(::AbstractShape, ::AbstractShape)

Shapes of a pipeline source and tagret.
"""
struct Signature
    src::AbstractShape
    tgt::AbstractShape
end

Signature() = Signature(NoShape(), AnyShape())

quoteof(sig::Signature) =
    Expr(:call, nameof(Signature), quoteof(sig.src), quoteof(sig.tgt))

syntaxof(sig::Signature) =
    Expr(:(->), syntaxof(sig.src), syntaxof(sig.tgt))

show(io::IO, sig::Signature) =
    print_expr(io, quoteof(sig))

source(sig::Signature) = sig.src

target(sig::Signature) = sig.tgt

#
# Rendering as a tree.
#

print_graph(shp::AbstractShape) =
    print_graph(stdout, shp)

function print_graph(io::IO, shp::AbstractShape; indent=0)
    gr = graphof(shp)
    w = 0
    for (ind, name, descr) in gr
        w = max(w, ind*2 + textwidth(name))
    end
    for (k, (ind, name, descr)) in enumerate(gr)
        print(io, " " ^ indent)
        for j = 1:ind
            more = false
            k′ = k + 1
            while k′ <= length(gr)
                ind′ = gr[k′][1]
                if ind′ == j
                    more = true
                end
                ind′ > j || break
                k′ += 1
            end
            print(io, j < ind ? (more ? "│ " : "  ") : (more ? "├╴" : "└╴"))
        end
        print(io, name)
        if !isempty(descr)
            print(io, " " ^ (2 + w - ind*2 - textwidth(name)))
            print(io, descr)
        end
        println(io)
    end
end

function graphof(shp::AbstractShape)
    gr = Tuple{Int,String,String}[]
    graphof!(gr, shp, nothing, 0, "#")
    gr
end

function graphof!(gr, shp::AbstractShape, card, ind, name)
    descr = "$(syntaxof(shp))"
    if card !== nothing
        descr = "$card × $descr"
    end
    push!(gr, (ind, name, descr))
end

function graphof!(gr, shp::BlockOf, card, ind, name)
    card′ = "$(syntaxof(cardinality(shp)))"
    if card !== nothing
        card′ = "$card × $card′"
    end
    graphof!(gr, elements(shp), card′, ind, name)
end

function graphof!(gr, shp::TupleOf, card, ind, name)
    descr = card !== nothing ? card : ""
    push!(gr, (ind, name, descr))
    for j = 1:width(shp)
        graphof!(gr, column(shp, j), nothing, ind+1, string(label(shp, j)))
    end
    gr
end

graphof!(gr, shp::IsLabeled, card, ind, name) =
    graphof!(gr, subject(shp), card, ind, string(label(shp)))

graphof!(gr, shp::Annotation, card, ind, name) =
    graphof!(gr, subject(shp), card, ind, name)


text:
_pr(t::Dict, x::Int, y::Int, z::Int) = join((rstrip(join(t[(n, m)] for n in range(0, 3+x+z))) for m in reverse(range(0, 3+y+z))), "\n")

function cuboid(x::Int, y::Int, z::Int)
    t = Dict((n, m) => " " for n in range(0, 3 + x + z), m in range(0, 3 + y + z))
    xrow = vcat("+", collect("$(i % 10)" for i in range(0, x)), "+")
    for (i, ch) in enumerate(xrow) t[(i, 0)] = t[(i, 1+y)] = t[(1+z+i, 2+y+z)] = ch end
    yrow = vcat("+", collect("$(j % 10)" for j in range(0, y)), "+")
    for (j, ch) in enumerate(yrow) t[(0, j)] = t[(x+1, j)] = t[(2+x+z, 1+z+j)] = ch end
    zdep = vcat("+", collect("$(k % 10)" for k in range(0, y)), "+")
    for (k, ch) in enumerate(xrow) t[(k, 1+y+k)] = t[(1+x+k, 1+y+k)] = t[(1+x+k, k)] = ch end

    return _pr(t, x, y, z)
end

for (x, y, z) in [(2, 3, 4), (3, 4, 2), (4, 2, 3), (5, 5, 6)]
    println("\nCUBOID($x, $y, $z)\n")
    println(cuboid(x, y, z))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export solarposition, zenith, zenith_azimuth, sines_and_cosines

# Calculate the position of the sun as (zenith,azimuth) in degrees for a location
# on earth given as (latitude,longitude) in degrees. The DateTime should be given
# as Universal Time, e.g. UTC.
function solarposition_wiki(datetime::DateTime, latitude, longitude)
	# Position of the Sun as seen from Earth.
	# Ported from Matlab, original helpstring below.

	#   This is the most basic algorithm. It is documented in Seinfeld &
	#   Pandis, Duffie & Beckman and Wikipedia.
	#
	# [ANGLES,PROJECTION] = SOLARPOSITION(DATE,TIME,LATITUDE,LONGITUDE,TIME_ZONE)
	# returns ZENITH & AZIMUTH for all DATE & TIME pairs at LATITUDE, LONGITUDE.
	# ANGLES = [ZENITH,AZIMUTH] and PROJECTION = [PHI_X, PHI_Y]
	# PHI_X is projection on x-z plane & PHI_Y is projection on y-z plane.
	# DATETIME can be string, vector [YEAR, MONTH, DAY, HOURS, MINUTES, SECONDS],
	#   cellstring or matrix N x [YEAR, MONTH, DAY, HOURS, MINUTES, SECONDS] for N
	#   times.
	# LATITUDE [degrees] and LONGITUDE [degrees] are the coordinates of the site.
	# TIME_ZONE [hours] of the site.
	# ROTATION [degrees] clockwise rotation of system relative to north.
	# DST [logical] flag for daylight savings time, typ. from March to November
	#   in the northern hemisphere.
	#
	# References:
	# http://en.wikipedia.org/wiki/Solar_azimuth_angle
	# http://en.wikipedia.org/wiki/Solar_elevation_angle
	#
	# Mark A. Mikofski
	# Copyright (c) 2013

	# Equation of Time
	# "Alternative calculation" on wikipedia:
	# https://en.wikipedia.org/wiki/Equation_of_time#Calculating_the_equation_of_time
	# Supposedly accurate to 6 seconds (0.1 minutes)

	year = Dates.year(datetime)
	# date in days starting at zero on 1 January (the subtraction produces a difference in milliseconds)
	d_ms = datetime - DateTime(year,1,1)
	d = d_ms.value/1000/3600/24

	# Earth's mean angular orbital velocity [degrees/day]
	w = 360/365.24
	# angle the earth moves on its orbit at its average speed from the December solstice
	# 10 is the approximate number of days from the December solstice to January 1
	a = w * (d+10)
	# The angle the Earth moves from the solstice to date D, including a first-order correction
	# for the Earth's orbital eccentricity, 0.0167. The number 2 is the number of days from 1 January
	# to the date of the Earth's perihelion.
	b = a + (360/pi*0.0167)*sind(w*(d-2))
	# C is the difference between the angles moved at mean speed, and at the corrected speed projected
	# onto the equatorial plane, and divided by 180 to get the difference in "half turns". The value 23.44°
	# is the obliquity (tilt) of the Earth's axis. The subtraction gives the conventional sign to the equation of time.
	c = (a - atand(tand(b)/cosd(23.44))) / 180
	# For any given value of x, arctan x (sometimes written as tan−1 x) has multiple values, differing from each other
	# by integer numbers of half turns. This may cause C to be wrong by an integer number of half turns. The excess
	# half turns are removed in the next step of the calculation to give the equation of time:
	ET = 720*(c - round(c))

	# approximate solar time [hours]
	solarTime = 24*(d - floor(d)) + longitude*24/360 + ET/60
	t_h = 15*(solarTime - 12) # [degrees] hour angle

	# declination [degrees]
	# accurate to 0.2 degrees, see https://en.wikipedia.org/wiki/Position_of_the_Sun#Calculations
	delta = -asind(sind(23.44)*cosd(b))
	zenith = acosd(sind(latitude)*sind(delta) + cosd(latitude)*cosd(delta)*cosd(t_h)) # [degrees] zenith

	# azimuth [0, 180], absolute value measured from due south, so east = west = 90,
	# south = 0, north = 180
	cos_phi = (cosd(zenith)*sind(latitude) - sind(delta)) ./ (sind(zenith)*cosd(latitude)) # cosine(azimuth)
	phi_south = acosd(clamp(cos_phi,-1,1))
	#phi_south = acosd(cos_phi);

	# azimuth [0, 360], measured clockwise from due north,
	# so east = 90, south = 180, and west = 270 degrees
	azimuth = 180 + sign(t_h)*phi_south # Shift domain to 0-360 deg

	#angles = [theta, phi]; # [degrees] zenith, azimuth
	return zenith, azimuth, solarTime, delta
end

function solarposition_wiki_faster(datetime::DateTime, latitude, longitude)
	year = Dates.year(datetime)
	d_ms = datetime - DateTime(year,1,1)
	d = d_ms.value/1000/3600/24
	w = 360/365.24
	a = w * (d+10)
	b = a + (360/pi*0.0167)*sind(w*(d-2))
	c = (a - atand(tand(b)/cosd(23.44))) / 180
	ET = 720*(c - round(c))
	solarTime = 24*(d - floor(d)) + longitude*24/360 + ET/60
	t_h = 15*(solarTime - 12) # [degrees] hour angle
	delta = -asind(sind(23.44)*cosd(b))
	slat, clat = sind(latitude), cosd(latitude)
	sdel, cdel = sind(delta), cosd(delta)
	czen = slat*sdel + clat*cdel*cosd(t_h)
	zenith = acosd(czen) # [degrees] zenith

	# azimuth [0, 180], absolute value measured from due south, so east = west = 90,
	# south = 0, north = 180
	cos_phi = (czen*slat - sdel) / (sind(zenith)*clat) # cosine(azimuth)
	phi_south = acosd(clamp(cos_phi,-1,1))
	#phi_south = acosd(cos_phi);

	# azimuth [0, 360], measured clockwise from due north,
	# so east = 90, south = 180, and west = 270 degrees
	azimuth = 180 + sign(t_h)*phi_south # Shift domain to 0-360 deg

	#angles = [theta, phi]; # [degrees] zenith, azimuth
	return zenith, azimuth, solarTime, delta
end

# Splits the DateTime into (year, month, day) (as integers) and the hour (as a floating point)
function splitdatetime(datetime::DateTime)
	# date = floor(datetime, Dates.Day)
	y, m, d = Dates.yearmonthday(datetime)
	return y, m, d, (datetime - DateTime(y,m,d)).value/1000/3600
end

# time correction term, difference between Terrestrial Time and Universal Time, see Grena (2012) section 2.1.
# Conversion to TT is necessary since UT follows the rotation of the earth, which is not uniform.
# t [days] days from January 1, 2060 (in UT)
# returns Δτ [seconds]
differenceTTUT(t) = 96.4 + 0.00158*t

# Calculate time variable needed for Grena (2012) algorithms, see section 3.1.
# The time t is the input DateTime expressed as days from January 1, 2060 (in UT).
# returns t [days]
calctime(datetime::DateTime) = (datetime-DateTime(2060,1,1)).value/1000/3600/24

# Same result as above, but this is the function in Grena's paper. The Julia version is faster.
function calctime_grena(datetime::DateTime)
	y, m, d, h = splitdatetime(datetime)
	if m <= 2
		m += 12
		y -= 1
	end
	return trunc(365.25*(y-2000)) + trunc(30.6001*(m+1)) - trunc(0.01*y) + d + 0.0416667*h - 21958
end

# Full version of algorithm 1 of Grena (2012), see section 3.2.
# Calculate the position of the sun as (zenith,azimuth) in radians for a location
# on earth given as (latitude,longitude) in degrees. The DateTime should be given
# as Universal Time, e.g. UTC.
function solarposition_grena1(datetime::DateTime, latitude, longitude)
	ω = 0.017202786 # [1/day]
	t = calctime(datetime)
	t_e = t + 1.1574e-5 * differenceTTUT(t)
	s1, c1 = sin(ω*t_e), cos(ω*t_e)
	s2, c2 = 2*s1*c1, (c1+s1)*(c1-s1)
	α = mod(-1.38880 + 1.72027920e-2*t_e + 3.199e-2*s1 - 2.65e-3*c1 + 4.050e-2*s2 + 1.525e-2*c2, 2π)
	δ = 6.57e-3 + 7.347e-2*s1 - 3.9919e-1*c1 + 7.3e-4*s2 - 6.60e-3*c2
	H = mod(1.75283 + 6.3003881*t + deg2rad(longitude) - α + π, 2π) - π
	sϕ, cϕ = sind(latitude), cosd(latitude)
	sδ, cδ = sin(δ), cos(δ)
	sH, cH = sin(H), cos(H)
	se0 = sϕ*sδ + cϕ*cδ*cH
	ep = asin(se0) - 4.26e-5*sqrt(1-se0^2)  #cos(e0)
	azimuth = atan(sH, cH*sϕ - sδ*cϕ/cδ)	#azimuth = 0 towards south and positive direction towards west
	zenith = π/2 - ep
	return zenith, azimuth				# radians
	# return zenith, azimuth, δ, H		# radians
	# return rad2deg(zenith), rad2deg(mod(azimuth+π,2π)), rad2deg(H)/15+12, rad2deg(δ) # degrees, same azimuth origin as solarposition0
end

# First half of algorithm 1 of Grena (2012), see section 3.2.
# Calculate the position of the sun as (declination, hour angle) in radians for a location
# on earth given by its longitude in degrees. The DateTime should be given
# as Universal Time, e.g. UTC.
function solarposition(datetime::DateTime, longitude)
	ω = 0.017202786 # [1/day]
	t = calctime(datetime)
	t_e = t + 1.1574e-5 * differenceTTUT(t)
	s1, c1 = sin(ω*t_e), cos(ω*t_e)
	s2, c2 = 2*s1*c1, (c1+s1)*(c1-s1)
	α = mod(-1.38880 + 1.72027920e-2*t_e + 3.199e-2*s1 - 2.65e-3*c1 + 4.050e-2*s2 + 1.525e-2*c2, 2π)
	δ = 6.57e-3 + 7.347e-2*s1 - 3.9919e-1*c1 + 7.3e-4*s2 - 6.60e-3*c2
	H = mod(1.75283 + 6.3003881*t + deg2rad(longitude) - α + π, 2π) - π
	return δ, H
end

# Full version of algorithm 1 of Grena (2012), see section 3.2.
# Calculate the position of the sun as (zenith, azimuth) in degrees for a location
# on earth given as (latitude,longitude) in degrees. The DateTime should be given
# as Universal Time, e.g. UTC.
function solarposition(datetime::DateTime, latitude, longitude)
	δ, H = solarposition(datetime, longitude)
	zenith, azimuth = zenith_azimuth(latitude, sines_and_cosines(δ, H)...)
	return rad2deg(zenith), rad2deg(azimuth)
end

sines_and_cosines(δ, H) = sin(δ), cos(δ), sin(H), cos(H)

function zenith_azimuth(latitude, sδ, cδ, sH, cH)
	sϕ, cϕ = sind(latitude), cosd(latitude)
	se0 = sϕ*sδ + cϕ*cδ*cH
	azimuth = atan(sH, cH*sϕ - sδ*cϕ/cδ)	# azimuth = 0 towards south and positive direction towards west
	zenith = π/2 - asin(se0) + 4.26e-5*sqrt(1-se0^2)  # sqrt(1-se0^2) = cos(e0), but slightly faster
	return zenith, azimuth		# radians
end

function zenith(latitude, sδ, cδ, sH, cH)
	sϕ, cϕ = sind(latitude), cosd(latitude)
	se0 = sϕ*sδ + cϕ*cδ*cH
	return π/2 - asin(se0) + 4.26e-5*sqrt(1-se0^2)  # sqrt(1-se0^2) = cos(e0), but slightly faster
end

const solarconstant = 1361  # W/m2

# Annual variation of solar insolation outside the atmosphere (W/m2):
# https://www.itacanet.org/the-sun-as-a-source-of-energy/part-2-solar-energy-reaching-the-earths-surface/#2.1.-The-Solar-Constant
solarinsolation(dt::DateTime) = solarconstant * (1 + 0.034*cos(2*π*(dt - DateTime(year(dt))).value/1000/3600/24/365.25))

# dt = DateTime(2019,4,15,13,46)
# δ, H = solarposition(dt, 12)	# return rad2deg(zenith), rad2deg(mod(azimuth+π,2π)), rad2deg(H)/15+12, rad2deg(δ)

# pos = sines_and_cosines(δ, H)
# zen = zenith(58, pos...)

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
#!/usr/bin/env julia

module RdRand

export rdrand16, rdrand32, rdrand64, rdrand

const librdrand = joinpath(dirname(@__FILE__), "..", "deps", "librdrand")

rdrand_cpucheck() = ccall( (:cpucheck,librdrand), Int, ())

if rdrand_cpucheck() != 0
    error("RdRand: unsupported cpu")
end

rdrand16() = ccall( (:rdrand16,librdrand), UInt16, ())
rdrand32() = ccall( (:rdrand32,librdrand), UInt32, ())
rdrand64() = ccall( (:rdrand64,librdrand), UInt64, ())

#rdrand64 should be uniform from 0 to 2^64
#but rdrand() isn't quite uniform on [0,1].
#What does it mean for a floating point number to be
#"uniform" from [0,1]? Here the smallest non-zero
#number is 1/2^64, which is about 5x10^-20.
#Since floating point numbers can represent much
#smaller quantities, should some greater care be taken?
#This function does pass diehard suite of tests.
rdrand() = rdrand64()/typemax(UInt64)

end


text:
using Aligator
using Base.Test

# write your own tests here
include("testpetkovsek.jl")


text:
using Documenter, StatsModels

DocMeta.setdocmeta!(StatsModels, :DocTestSetup, :(using StatsModels, StatsBase); recursive=true)

using Pkg
Pkg.precompile()

makedocs(
    sitename = "StatsModels.jl",
    pages = [
        "Introduction" => "index.md",
        "Modeling tabular data" => "formula.md",
        "Internals and extending the `@formula`" => "internals.md",
        "Contrast coding categorical variables" => "contrasts.md",
        "Temporal variables and Time Series Terms" => "temporal_terms.md",
        "API documentation" => "api.md"
    ],
    modules = [StatsModels]
)

deploydocs(
    repo = "github.com/JuliaStats/StatsModels.jl.git",
    push_preview = true,
)


text:
#=
   Ruin problem.

   This version use a fixed max time (length) of the sequences, here
   maxT.

   Summary Statistics
   parameters      mean       std   naive_se      mcse          ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64      Float64   Float64       Float64

      ...
         len    4.6035    4.5269     0.0453    0.0431    9406.7545    1.0000     3560.4673
    amount_t    0.9586    2.0299     0.0203    0.0203   10001.5751    1.0003     3785.6075
      ruin_t    0.7785    0.4153     0.0042    0.0040    9748.3505    1.0000     3689.7618

   Distributions of variable len (num:0)
   1.00000 =>    5004  (0.500400)
   12.00000 =>    2215  (0.221500)
   3.00000 =>    1255  (0.125500)
   5.00000 =>     660  (0.066000)
   7.00000 =>     386  (0.038600)
   9.00000 =>     298  (0.029800)
   11.00000 =>     182  (0.018200)

   Distributions of variable amount_t (num:0)
   0.00000 =>    7785  (0.778500)
   4.00000 =>     801  (0.080100)
   2.00000 =>     658  (0.065800)
   6.00000 =>     527  (0.052700)
   8.00000 =>     195  (0.019500)
   10.00000 =>      32  (0.003200)
   12.00000 =>       2  (0.000200)

   Distributions of variable ruin_t (num:0)
   1.00000 =>    7785  (0.778500)
   0.00000 =>    2215  (0.221500)

.  Cf ~/blog/ruin_problem.blog
      ~/webppl/ruin_problem.wppl

=#
using Turing, StatsPlots
include("jl_utils.jl")

@model function ruin_problem(maxT=12,start=1,win=1,loose=1)
    head = 1
    tail = 2
    coins = [head, tail]

    # Draw a coin at timestep t
    draw = tzeros(maxT)
    for t in 1:maxT
        if t == 1
            draw[t] = head
        else
            draw[t] ~ Categorical([0.5,0.5])
        end
    end

    # What is the score at time t?
    # After ruin there's no way back at the game again.
    amount = tzeros(maxT)
    for t in 1:maxT
        if t == 1
            amount[t] = start
        elseif amount[t-1] <= 0
            amount[t] = 0
        elseif draw[t] == head
            amount[t] = amount[t-1] + win
        else
            if amount[t-1] - loose < 0
                amount[t] =  0
            else
                amount[t] = amount[t-1] - loose
            end
        end
    end

    # Is the player busted at time t?
    # This is the first time the player bused
    ruin = tzeros(maxT)
    for t in 1:maxT
        if amount[t] <= 0
            ruin[t] = 1
        else
            ruin[t] = 0
        end
    end

    len ~ Dirac(length(filter(i -> i > 0, amount)))
    amount_t ~ Dirac(amount[maxT])
    ruin_t ~ Dirac(ruin[maxT])

    # true ~ Dirac(ruin[5]==true)
    # true ~ Dirac(amount[maxT]==8)

end

maxT=12
start=1
win=1
loose=1
model = ruin_problem(maxT,start,win,loose)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

show_var_dist_pct(chns, :len)
show_var_dist_pct(chns, :amount_t)
show_var_dist_pct(chns, :ruin_t)

# chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
# genq = generated_quantities(model, chains_params)
# show_var_dist_pct(genq)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Distributed

workers = addprocs(20)

@everywhere begin
    using HMMs
    using Brobdingnag
    using SpecialFunctions
    using Random
    using Distributed
    using Models
end

@everywhere function fithmm( n, nhmm, seed )
    y = rand( nhmm, 1:n, seed=seed  );
    hmm = deepcopy( nhmm )
    HMMs.setobservations( hmm, y )
    HMMs.em( hmm, debug=2 )
    p = Models.getcompressedparameters( hmm )
    C = Models.sandwich( hmm )
    return (p, C)
end

hmmtype = HMMs.HMM{1,HMMs.GenTDist,Brob,Float64,Int}
nhmm = rand( hmmtype )
n = 1_000_000

N = 100
data1 = pmap( fithmm, fill(n,N), fill(nhmm,N), 1:N )

using StatsBase

p = Models.getcompressedparameters( nhmm )
vs = Vector{Float64}[]
for i = 1:N
    (p1, C1) = data1[i]
    sC1 = real.(sqrt(C1))
    @assert( maximum(abs.(C1 - sC1 * sC1')) < 1e-10 )
    push!( vs, inv(sC1) * (p - p1) )
end

hmmtype = HMMs.HMM{2,HMMs.GenTDist,Brob,Float64,Int}
nhmm = rand( hmmtype )
n = 1_000_000

N = 100
data2 = pmap( fithmm, fill(n,N), fill(nhmm,N), 1:N )

p = Models.getcompressedparameters( nhmm )
vs = Vector{Float64}[]
for i = 1:N
    (p1, C1) = data1[i]
    sC1 = real.(sqrt(C1))
    @assert( maximum(abs.(C1 - sC1 * sC1')) < 1e-10 )
    push!( vs, inv(sC1) * (p - p1) )
end


text:
#! /usr/bin/julia
#
#	csv_delete.jl
#
#					Jul/05/2018
# --------------------------------------------------------------------
include("text_manipulate.jl")
#
# --------------------------------------------------------------------
println("*** 開始 ***")
fname_csv = ARGS[1]
key_in = ARGS[2]
println(fname_csv)
println(key_in)
#
dict_aa = csv_to_dict_proc(fname_csv)
#
if ( haskey(dict_aa,key_in))
	delete!(dict_aa, key_in)
	dict_to_csv_proc(dict_aa,fname_csv)
end
#
println("*** 終了 ***")
# --------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
# export Subotnik_A

# @doc raw"""
#     Model A according to Coffmann & Subotnik, Phys.Chem.Chem.Phys., 20, 9847 (2018).
#     V_0 =  0.5*m*ω^2 r^2
#     V_1 =  0.5*m*ω^2 (r-g)^2 + ΔG
#     Molecule-metal coupling:
#     Γ(x) = Γ_1(Γ_0 (K(r-d)^2)/(1+K(r-d)^2))
#     d is the minimum of the PES, d = g/2 +ΔG/(mω^2g)
#     In the original model, K = 0.1 and Γ_0=0.01; Γ_1 is used as an energy unit
#     The numbers given as default correspond to the one listed in the paper.
# """
# @with_kw struct Subotnik_A <: DiabaticModel
#     n_states::UInt8 = 2
#     m::Float64 = 2000
#     om::Float64 = 2e-4
#     g::Float64 = 20.6097
#     DG::Float64 = -0.0038
#     G1::Float64 = 1e-4
#     G0::Float64 = 1e-2
#     K:: Float64 = 0.1
#     d::Float64 = 8
# end

# """
# Diabatic potenial matrix
# """
# function potential!(model::Subotnik_A, V::Hermitian, R::AbstractMatrix)
#     # Update d
#     d = model.g/2 + model.DG/(model.m*model.om^2*model.g)

#     # Higher potential well
#     V_0(r) = (model.m*model.om^2*r^2)/2.0
#     # Loweer energy potential well
#     V_1(r) = (model.m*model.om^2*(r-model.g)^2)/2.0 + model.DG
#     #println(R[1]," " , V_0(R[1]), " ", V_1(R[1]))
#     # weighted metal-molecule coupling between the diabates
#     G(r)=model.G1*(model.G0 + (model.K*(r - model.d)^2)/(1.0 + model.K*(r -model.d)^2))

#     #V = Hermitian(zeros(2,2))
#     # Set up the Hamiltonian
#     V[1,1] = V_0(R[1])
#     V[2,2] = V_1(R[1])
#     V.data[1,2] = G(R[1])
#     V.data[2,1] = G(R[1])
# end

# """
# Elementwise position derivative of the above potential
# """
# function derivative!(model::Subotnik_A, derivative::AbstractMatrix{<:Hermitian}, R::AbstractMatrix)

#     q = R[1]
#     derivative[1][1,1] = model.m*model.om^2*q
#     derivative[1][2,2] =model.m*model.om^2*(q-model.g)
#     derivative[1].data[1,2] = model.G1*((2*model.K*(q-model.d)*(1+model.K*(q-model.d)^2)-2*model.K^2*(q-model.d)^2)/
#     (1+model.K*(q-model.d)^2)^2)

#     derivative[1].data[2,1] = model.G1*((2*model.K*(q-model.d)*(1+model.K*(q-model.d)^2)-2*model.K^2*(q-model.d)^2)/
#     (1+model.K*(q-model.d)^2)^2)

# end

@doc raw"""
    Model MiaoSubotnik according to Miao & Subotnik, J. Chem Phys., 150, 041711 (2019).
    V_0 =  0.5*m*ω^2 x^2
    V_1 =  0.5*m*ω^2 (x-g)^2 + ΔG
    Molecule-metal coupling:
    V = sqrt(\Gamma*2*W/(2*\pi*M))
    d is the minimum of the PES, d = g/2 +ΔG/(mω^2g)
    2W is the bandwidth
    In the original model, m = 2000, ω = 0.0002, g = 20.6097, ΔG = -0.0038 and
    kBT=0.00095. In the paper, Γ, 2W and M were varied and tested.
    The numbers given as default correspond to the one listed in the paper.
    The model is defined such that h=(V_0-V_1)d^t d + V_1, however, the diabatic surfaces
    that they show in the paper in Fig. 1 correspond to the picture implemented below
"""
Parameters.@with_kw struct MiaoSubotnik <: LargeDiabaticModel
    m::Float64 = 2000
    ω::Float64 = 2e-4
    g::Float64 = 20.6097
    ΔG::Float64 = -3.8e-3
    Γ::Float64 = 6.4e-3
    W::Float64 = 5Γ
    M::Int = 40
    ρ::Float64 = M/2W
    increment::Float64 = 1/ρ
    n_states::Int = M+1
end

NQCModels.ndofs(::MiaoSubotnik) = 1
NQCModels.nstates(model::MiaoSubotnik) = model.n_states
NQCModels.nelectrons(model::MiaoSubotnik) = div(model.M, 2)

function NQCModels.potential!(model::MiaoSubotnik, V::Hermitian, R::Real)
    Parameters.@unpack m, ω, g, ΔG, Γ, ρ, increment, n_states = model

    U0(x) = (m*ω^2*x^2)/2
    U1(x) = (m*ω^2*(x-g)^2)/2 + ΔG

    V[1,1] = U1(R) - U0(R)

    Vₙ = sqrt(Γ/(2π*ρ))
    V.data[1,2:end] .= Vₙ

    V[2,2] = -model.W
    for i=3:n_states
        V[i,i] = V[i-1,i-1] + increment
    end

    return V
end

function NQCModels.derivative!(model::MiaoSubotnik, derivative::Hermitian, R::Real)
    Parameters.@unpack m, ω, g, n_states = model

    dU0(x) = m*ω^2*x
    dU1(x) = m*ω^2*(x-g)

    derivative[1,1] = dU1(R) - dU0(R)

    return derivative
end

function NQCModels.state_independent_potential(model::MiaoSubotnik, r::AbstractMatrix)
    (;m, ω) = model
    return 1/2 * m * ω^2 * r[1]^2
end

function NQCModels.state_independent_derivative!(model::MiaoSubotnik, derivative::AbstractMatrix, r::AbstractMatrix)
    (;m, ω) = model
    derivative[1] = m * ω^2 * r[1]
    return derivative
end


text:
using DifferentialEquations
using LinearAlgebra
using DiffEqSensitivity
using Zygote
using Statistics
using Distributions

# include("./AssimilationProblem.jl")

"""
    timeindex(t, Δt)

For a given time value and time step, return the corresponding array index.
"""
function timeindex(t, Δt)
    return Int(round(t/Δt)) + 1
end

"""
    update!(sol, t_now, t_next, Δt, uₐ)

Update solution array uₐ with the values from the solution object sol at specified time step Δt
"""
function update!(sol, t_now, t_next, Δt, uₐ)
    idx_now = timeindex(t_now, Δt)
    idx_next = timeindex(t_next, Δt)

    uₐ[:, idx_now:idx_next] = hcat(sol(t_now:Δt:t_next)...)
end
Zygote.@nograd update!  # tell zygote to ignore update! when computing derivatives

"""
    model_forward!(uₖ, uₐ, t_now, t_next, prob)

Given value uₖ and time range t_now and t_next, integrate the ODE problem forward and update solution object.
"""
function model_forward!(uₖ, uₐ, t_now, t_next, prob)
    ode_prob = ODEProblem(prob.f, uₖ, (t_now, t_next), prob.p)
    sol = solve(ode_prob, Tsit5(), reltol=1e-6,abstol=1e-6)

    update!(sol, t_now, t_next, prob.Δt, uₐ)

    return sol[:, end]
end

"""
    EKFsolve(prob::AssimilationProblem)

Given an AssimilationProblem `prob`, obtain the analysis uₐ and error covariance matrix B using the Extended Kalman Filter (EKF) method.
"""
function EKFsolve(prob::AssimilationProblem)
    # 1. allocate solution arrays for u(t) and B(t)
    uₐ = zeros(size(prob.u₀, 1), size(prob.tspan[1]:prob.Δt:prob.tspan[2], 1))
    B = zeros(size(prob.u₀,1), size(prob.u₀,1), size(prob.tspan[1]:prob.Δt:prob.tspan[2], 1))

    # 2. Set initial values
    uₐ[:,1] = prob.u₀
    B[:,:,1] = prob.B₀

    # 3. define time ranges
    times = prob.tspan[1]:prob.Δt:prob.tspan[2]

    n_obs = size(prob.w, 2)
    times_obs = prob.Δt_obs:prob.Δt_obs:(prob.Δt_obs * (n_obs-1))

    for k ∈ 1:(length(times)-1)

        # 4. Perform the forecast to obtain the background

        t_now = times[k]
        t_next = times[k+1]

        idx_now = timeindex(t_now, prob.Δt)
        idx_next = timeindex(t_next, prob.Δt)

        u_now = uₐ[:, idx_now]

        u_next, DM = Zygote.withjacobian(model_forward!, u_now, uₐ, t_now, t_next, prob)
        DM = DM[1]

        # update error covariance matrix
        B[:,:,idx_next] = DM * B[:,:,idx_now] * DM' + prob.Q

        # 5. Perform the assimilation step
        if t_next ∈ times_obs
            n = size(uₐ, 1)
            idx_h = timeindex(t_next, prob.Δt_obs) - 1 # -1 since this doesn't include 0

            # compute the Kalman gain matrix
            DH = prob.Dh(uₐ[:,idx_next])

            denom = DH*B[:,:, idx_next]*DH' + prob.R

            K = B[:,:, idx_next]*DH'*inv(denom)

            uₐ[:,idx_next] .= uₐ[:, idx_next] + K*(prob.w[:, idx_h] - prob.h(uₐ[:, idx_next]))

            B[:,:,idx_next] = (I(n) - K*DH)*B[:,:,idx_next]
        end

    end

    # 6. Return the resulting analysis uₐ, and error covariance matrix B
    return uₐ, B
end


text:
using Test
using LinearAlgebra
using Distributions
using DataDeps
using ParameterEstimocean

using Oceananigans
using Oceananigans.Units
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity, MixingLength

using ParameterEstimocean.Parameters: unconstrained_prior, transform_to_constrained

@testset "Forward map evaluation with LESbrary observations" begin
    @testset "Convective-adjustment-based ensemble simulations" begin

        data_path = datadep"two_day_suite_4m/free_convection_instantaneous_statistics.jld2"
        times = [1hours, 2hours]
        observations = SyntheticObservations(data_path; field_names=:b, times, transformation=ZScore())
        Nensemble = 10

        simulation = ensemble_column_model_simulation(observations; Nensemble, tracers=:b,
                                                      closure = ConvectiveAdjustmentVerticalDiffusivity())

        Qᵘ = simulation.model.velocities.u.boundary_conditions.top.condition
        Qᵇ = simulation.model.tracers.b.boundary_conditions.top.condition
        N² = simulation.model.tracers.b.boundary_conditions.bottom.condition

        simulation.Δt = 10.0

        Qᵘ .= observations.metadata.parameters.momentum_flux
        Qᵇ .= observations.metadata.parameters.buoyancy_flux
        N² .= observations.metadata.parameters.N²_deep

        priors = (convective_κz = ScaledLogitNormal(bounds=(0.1, 1.0)),
                  background_κz = ScaledLogitNormal(bounds=(0.0, 1e-5)))

        free_parameters = FreeParameters(priors)
        calibration = InverseProblem(observations, simulation, free_parameters)

        θ1 = (convective_κz=0.25, background_κz=5e-5)
        θ2 = (convective_κz=1.00, background_κz=1e-4)

        y = observation_map(calibration)
        G = forward_map(calibration, [θ1, θ1, θ2])

        @test G[:, 1] == G[:, 2]
        @test G[:, 2] != G[:, 3]
        @test all(G[:, 3] == G[:, i] for i = 4:Nensemble)
    end

    @testset "CATKE-based ensemble simulations" begin
        data_path = datadep"two_day_suite_4m/weak_wind_strong_cooling_instantaneous_statistics.jld2"
        times = [1hours, 2hours]
        field_names = (:b, :u, :v, :e)
        transformation = (b=ZScore(), u=ZScore(), v=ZScore(), e=RescaledZScore(0.1))
        observations = SyntheticObservations(data_path; field_names, times, transformation)

        catke = CATKEVerticalDiffusivity()

        simulation = ensemble_column_model_simulation(observations;
                                                      Nensemble = 10,
                                                      architecture = CPU(),
                                                      tracers = (:b, :e),
                                                      closure = catke)

        Qᵘ = simulation.model.velocities.u.boundary_conditions.top.condition
        Qᵇ = simulation.model.tracers.b.boundary_conditions.top.condition
        N² = simulation.model.tracers.b.boundary_conditions.bottom.condition

        simulation.Δt = 10.0

        Qᵘ .= observations.metadata.parameters.momentum_flux
        Qᵇ .= observations.metadata.parameters.buoyancy_flux
        N² .= observations.metadata.parameters.N²_deep

        Cᴰ_default = catke.Cᴰ

        priors = (; Cᴰ = lognormal(mean=Cᴰ_default, std=Cᴰ_default/100))

        free_parameters = FreeParameters(priors)
        calibration = InverseProblem(observations, simulation, free_parameters)

        Nensemble = simulation.model.grid.Nx
        names = free_parameters.names
        Nθ = length(names)
        unconstrained_priors = NamedTuple(name => unconstrained_prior(priors[name]) for name in names)

        # Initial sample
        X = [rand(unconstrained_priors[i]) for i=1:Nθ, k=1:Nensemble]

        @test size(X) == (Nθ, Nensemble)

        θ = [transform_to_constrained(priors, X[:, k]) for k=1:Nensemble]

        @test length(θ) == Nensemble
        @test all(length(θ[k]) == Nθ for k=1:Nensemble)

        G = forward_map(calibration, θ)
        y = observation_map(calibration)

        @test size(G, 1) == length(y)
        @test size(G, 2) == Nensemble

        θ1 = θ[1]
        G = forward_map(calibration, [θ1])
        @test all(G[:, 1] == G[:, i] for i = 2:Nensemble)

        θ2 = θ[2]
        @test values(θ1) != values(θ2)

        G = forward_map(calibration, [θ1, θ2])
        @test G[:, 1] != G[:, 2]
        @test all(G[:, 2] == G[:, i] for i = 3:Nensemble)

        G = forward_map(calibration, [θ1, θ1, θ2])
        @test G[:, 1] == G[:, 2]
        @test all(G[:, 3] == G[:, i] for i = 4:Nensemble)

        G = forward_map(calibration, [θ1, θ1, θ2, θ2])
        @test G[:, 1] == G[:, 2]
        @test G[:, 2] != G[:, 3]
        @test G[:, 3] == G[:, 4]
        @test all(G[:, 4] == G[:, i] for i = 5:Nensemble)
    end
end


text:
module LapackUtil

const liblapack = Base.liblapack_name

import LinearAlgebra.BLAS.@blasfunc

import LinearAlgebra: BlasFloat, BlasInt, LAPACKException,
    DimensionMismatch, SingularException, PosDefException, chkstride1, checksquare

using Base: iszero, has_offset_axes

export tgsyl!, lanv2, ladiv, lag2, lacn2!

function chklapackerror(ret::BlasInt)
    if ret == 0
        return
    elseif ret < 0
        throw(ArgumentError("invalid argument #$(-ret) to LAPACK call"))
    else # ret > 0
        throw(LAPACKException(ret))
    end
end

for (fn, elty, relty) in ((:dtgsyl_, :Float64, :Float64),
                   (:stgsyl_, :Float32, :Float32),
                   (:ztgsyl_, :ComplexF64, :Float64),
                   (:ctgsyl_, :ComplexF32, :Float32))
    @eval begin
        function tgsyl!(trans::AbstractChar, A::AbstractMatrix{$elty}, B::AbstractMatrix{$elty}, C::AbstractMatrix{$elty},
                        D::AbstractMatrix{$elty}, E::AbstractMatrix{$elty}, F::AbstractMatrix{$elty})
            @assert !has_offset_axes(A, B, C, D, E, F)
            chkstride1(A, B, C, D, E, F)
            m, n = checksquare(A, B)
            lda = max(1, stride(A, 2))
            ldb = max(1, stride(B, 2))
            ldc = max(1, stride(C, 2))
            m1, n1 = size(C)
            (m == m1 && n == n1) ||
                throw(DimensionMismatch("dimensions of A($m,$m),  B($n,$n), and C($m1,$n1) must match"))
            m2, n2 = checksquare(D, E)
            m == m2 || throw(DimensionMismatch("dimensions of A($m,$m) and D($m2,$m2) must match"))
            n == n2 || throw(DimensionMismatch("dimensions of B($n,$n) and E($n2,$n2), must match"))
            ldd = max(1, stride(D, 2))
            lde = max(1, stride(E, 2))
            ldf = max(1, stride(F, 2))
            m3, n3 = size(F)
            (m2 == m3 && n2 == n3) ||
                throw(DimensionMismatch("dimensions of D($m,$m),  E($n,$n), and F($m3,$n3) must match"))
            dif = Vector{$relty}(undef, 1)
            scale = Vector{$relty}(undef, 1)
            info  = Ref{BlasInt}()
            ijob = 0
            work = Vector{$elty}(undef, 1)
            lwork = 1
            iwork = Vector{BlasInt}(undef,m+n+6)
            #SUBROUTINE DTGSYL( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
            #       LDD, E, LDE, F, LDF, SCALE, DIF, WORK, LWORK,
            #       IWORK, INFO )
            ccall((@blasfunc($fn), liblapack), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                 Ptr{$relty}, Ptr{$relty}, Ptr{$relty}, Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
                 trans, ijob, m, n,
                 A, lda, B, ldb, C, ldc, D, ldd, E, lde, F, ldf,
                 scale, dif, work, lwork, iwork, info)
            chklapackerror(info[])
            C, F, scale[1]
        end
    end
end

"""
    tgsyl!(A, B, C, D, E, F) -> (C, F, scale)

Solve the Sylvester system of matrix equations

      AX - YB = scale*C
      DX - YE = scale*F ,

where `X` and `Y` are unknown matrices, the pairs `(A, D)`, `(B, E)` and  `(C, F)`
have the same sizes, and the pairs `(A, D)` and `(B, E)` are in
generalized (real) Schur canonical form, i.e. `A`, `B` are upper quasi
triangular and `D`, `E` are upper triangular.
Returns `X` (overwriting `C`), `Y` (overwriting `F`) and `scale`.

    tgsyl!(trans, A, B, C, D, E, F) -> (C, F, scale)

Solve for `trans = 'T'` and real matrices or for `trans = 'C'` and complex
matrices,  the (adjoint) Sylvester system of matrix equations

      A'X + D'Y = scale*C
      XB' + YE' = scale*(-F) .

`tgsyl!('N', A, B, C, D, E, F)` corresponds to the call `tgsyl!(A, B, C, D, E, F)`.

Interface to the LAPACK subroutines DTGSYL/STGSYL/ZTGSYL/CTGSYL.
"""
tgsyl!(trans::AbstractChar,A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix, E::AbstractMatrix, F::AbstractMatrix)

tgsyl!(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix, E::AbstractMatrix, F::AbstractMatrix) =
tgsyl!('N',A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, D::AbstractMatrix, E::AbstractMatrix, F::AbstractMatrix)

for (fn, elty) in ((:dlanv2_, :Float64),
                   (:slanv2_, :Float32))
    @eval begin
        function lanv2(A::$elty, B::$elty, C::$elty, D::$elty)
           """
           SUBROUTINE DLANV2( A, B, C, D, RT1R, RT1I, RT2R, RT2I, CS, SN )

           DOUBLE PRECISION A, B, C, CS, D, RT1I, RT1R, RT2I, RT2R, SN
           """
           RT1R = Ref{$elty}(1.0)
           RT1I = Ref{$elty}(1.0)
           RT2R = Ref{$elty}(1.0)
           RT2I = Ref{$elty}(1.0)
           CS = Ref{$elty}(1.0)
           SN = Ref{$elty}(1.0)
           ccall((@blasfunc($fn), liblapack), Cvoid,
                 (Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty},
                 Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty}),
                 A, B, C, D,
                 RT1R, RT1I, RT2R, RT2I, CS, SN)
           return RT1R[], RT1I[], RT2R[], RT2I[], CS[], SN[]
        end
    end
end
"""
    lanv2(A, B, C, D) -> (RT1R, RT1I, RT2R, RT2I, CS, SN)

Compute the Schur factorization of a real 2-by-2 nonsymmetric matrix `[A,B;C,D]` in
standard form. `A`, `B`, `C`, `D` are overwritten on output by the corresponding elements of the
standardised Schur form. `RT1R+im*RT1I` and `RT2R+im*RT2I` are the resulting eigenvalues.
`CS` and `SN` are the parameters of the rotation matrix.
Interface to the LAPACK subroutines DLANV2/SLANV2.
"""
lanv2(A::Union{Float32,Float64}, B::Union{Float32,Float64}, C::Union{Float32,Float64}, D::Union{Float32,Float64})

for (fn, elty) in ((:dlag2_, :Float64),
                   (:slag2_, :Float32))
    @eval begin
        function lag2(A::StridedMatrix{$elty}, B::StridedMatrix{$elty}, SAFMIN::$elty)
           """
           SUBROUTINE DLAG2( A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI )

           INTEGER            LDA, LDB
           DOUBLE PRECISION   SAFMIN, SCALE1, SCALE2, WI, WR1, WR2
           DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
           """
           LDA = stride(A,2)
           LDB = stride(B,2)
           SCALE1 = Ref{$elty}(1.0)
           SCALE2 = Ref{$elty}(1.0)
           WR1 = Ref{$elty}(1.0)
           WR2 = Ref{$elty}(1.0)
           WI = Ref{$elty}(1.0)
           ccall((@blasfunc($fn), liblapack), Cvoid,
                (Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt}, Ref{$elty},
                Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty}),
                A, LDA, B, LDB, SAFMIN,
                SCALE1, SCALE2, WR1, WR2, WI)
           return SCALE1[], SCALE2[], WR1[], WR2[], WI[]
        end
    end
end
"""
    lag2(A, B, SAFMIN) -> (SCALE1, SCALE2, WR1, WR2, WI)

Compute the eigenvalues of a 2-by-2 generalized real eigenvalue problem for
the matrix pair `(A,B)`, with scaling as necessary to avoid over-/underflow.
`SAFMIN` is the smallest positive number s.t. `1/SAFMIN` does not overflow.
If `WI = 0`, `WR1/SCALE1` and `WR2/SCALE2` are the resulting real eigenvalues, while
if `WI <> 0`, then `(WR1+/-im*WI)/SCALE1` are the resulting complex eigenvalues.
Interface to the LAPACK subroutines DLAG2/SLAG2.
"""
lag2(A::StridedMatrix{Union{Float32,Float64}}, B::StridedMatrix{Union{Float32,Float64}}, SAFMIN::Union{Float32,Float64})

for (fn, elty) in ((:dladiv_, :Float64),
                   (:sladiv_, :Float32))
    @eval begin
        function ladiv(A::$elty, B::$elty, C::$elty, D::$elty)
           """
           SUBROUTINE DLADIV( A, B, C, D, P, Q )

           DOUBLE PRECISION   A, B, C, D, P, Q
           """
           P = Ref{$elty}(1.0)
           Q = Ref{$elty}(1.0)
           ccall((@blasfunc($fn), liblapack), Cvoid,
                (Ref{$elty},Ref{$elty},Ref{$elty},Ref{$elty},
                 Ref{$elty},Ref{$elty}),
                 A, B, C, D, P, Q)
           return P[], Q[]
        end
    end
end
"""
    ladiv(A, B, C, D) -> (P, Q)

Perform the complex division in real arithmetic

  ``P + iQ = \\displaystyle\\frac{A+iB}{C+iD}``

by avoiding unnecessary overflow.
Interface to the LAPACK subroutines DLADIV/SLADIV.
"""
ladiv(A::Union{Float32,Float64}, B::Union{Float32,Float64}, C::Union{Float32,Float64}, D::Union{Float32,Float64})

for (fn, elty) in ((:dlacn2_, :Float64),
                   (:slacn2_, :Float32))
    @eval begin
        function lacn2!(V::AbstractVector{$elty}, X::AbstractVector{$elty}, ISGN::AbstractVector{BlasInt},
                        EST::$elty, KASE::BlasInt, ISAVE::AbstractVector{BlasInt})
            @assert !has_offset_axes(V, X)
            chkstride1(V,X)
            n = length(V)
            if n != length(X) || n != length(ISGN)
                throw(DimensionMismatch("dimensions of V,  X, and ISIGN must be equal"))
            end
            """
            *       SUBROUTINE DLACN2( N, V, X, ISGN, EST, KASE, ISAVE )
            *
            *       .. Scalar Arguments ..
            *       INTEGER            KASE, N
            *       DOUBLE PRECISION   EST
            *       ..
            *       .. Array Arguments ..
            *       INTEGER            ISGN( * ), ISAVE( 3 )
            *       DOUBLE PRECISION   V( * ), X( * )
            """
            KASE1 = Array{BlasInt,1}(undef,1)
            KASE1[1] = KASE
            EST1 = Vector{$elty}(undef, 1)
            EST1[1] = EST
            ccall((@blasfunc($fn), liblapack), Cvoid,
                (Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
                Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
                n, V, X, ISGN, EST1, KASE1, ISAVE)
            return EST1[1], KASE1[1]
        end
    end
end
"""
    lacn2!(V, X, ISGN, EST, KASE, ISAVE ) -> (EST, KASE )
Estimate the 1-norm of a `real` linear operator `A`, using reverse communication
by applying the operator or its transpose/adjoint to a vector.
`KASE` is a parameter to control the norm evaluation process as follows.
On the initial call, `KASE` should be `0`. On an intermediate return,
`KASE` will be `1` or `2`, indicating whether the `real` vector `X` should be overwritten
by `A * X`  or `A' * X` at the next call.
On the final return, `KASE` will again be `0` and `EST` is an estimate (a lower bound)
for the 1-norm of `A`. `V` is a real work vector, `ISGN` is an integer work
vector and `ISAVE` is a 3-dimensional integer vector used to save information
between the calls.
Interface to the LAPACK subroutines DLACN2/SLACN2.
"""
lacn2!(V::AbstractVector{Union{Float32,Float64}}, X::AbstractVector{Union{Float32,Float64}}, ISGN::AbstractVector{Int64}, EST::Union{Float32,Float64}, KASE::Int64, ISAVE::AbstractVector{Int64})

for (fn, elty, relty) in ((:zlacn2_, :ComplexF64, :Float64),
                          (:clacn2_, :ComplexF32, :Float32))
    @eval begin
        function lacn2!(V::AbstractVector{$elty}, X::AbstractVector{$elty},
                        EST::$relty, KASE::BlasInt, ISAVE::AbstractVector{BlasInt})
            @assert !has_offset_axes(V, X)
            chkstride1(V,X)
            n = length(V)
            n == length(X) || throw(DimensionMismatch("dimensions of V and X must be equal"))
            """
            *       SUBROUTINE ZLACN2( N, V, X, EST, KASE, ISAVE )
            *
            *       .. Scalar Arguments ..
            *       INTEGER            KASE, N
            *       DOUBLE PRECISION   EST
            *       ..
            *       .. Array Arguments ..
            *       INTEGER            ISGN( * ), ISAVE( 3 )
            *       COMPLEX*16         V( * ), X( * )
            """
            KASE1 = Array{BlasInt,1}(undef,1)
            KASE1[1] = KASE
            EST1 = Vector{$relty}(undef, 1)
            EST1[1] = EST
            ccall((@blasfunc($fn), liblapack), Cvoid,
                (Ref{BlasInt}, Ptr{$elty}, Ptr{$elty},
                Ptr{$relty}, Ptr{BlasInt}, Ptr{BlasInt}),
                n, V, X, EST1, KASE1, ISAVE)
            return EST1[1], KASE1[1]
        end
    end
end
"""
    lacn2!(V, X, EST, KASE, ISAVE ) -> (EST, KASE )
Estimate the 1-norm of a `complex` linear operator `A`, using reverse communication
by applying the operator or its adjoint to a vector.
`KASE` is a parameter to control the norm evaluation process as follows.
On the initial call, `KASE` should be `0`. On an intermediate return,
`KASE` will be `1` or `2`, indicating whether the `complex` vector `X` should be overwritten
by `A * X`  or `A' * X` at the next call.
On the final return, `KASE` will again be `0` and `EST` is an estimate (a lower bound)
for the 1-norm of `A`. `V` is a complex work vector and `ISAVE` is a 3-dimensional
integer vector used to save information between the calls.
Interface to the LAPACK subroutines ZLACN2/CLACN2.
"""
lacn2!(V::AbstractVector{Union{ComplexF32,ComplexF64}}, X::AbstractVector{Union{ComplexF32,ComplexF64}}, EST::Union{Float32,Float64}, KASE::Int64, ISAVE::AbstractVector{Int64})

end


text:
using SchumakerSpline
from = 0.5
to = 10
x1 = collect(range(from, stop=to, length=40))
y1 = (x1).^2
x2 = [0.5,0.75,0.8,0.93,0.9755,1.0,1.1,1.4,2.0]
y2 = sqrt.(x2)
left_spline = Schumaker(x1,y1)
right_spline = Schumaker(x2,y2)

crossover_point = get_intersection_points(left_spline,right_spline)

splice_point = crossover_point[1]
spliced = splice_splines(left_spline, right_spline, splice_point)
# Testing at splice point
abs(evaluate(spliced, splice_point) - evaluate(right_spline, splice_point))          < 100*eps()
# As this was a continuous split
abs(evaluate(spliced, splice_point- 100*eps()) - evaluate(spliced, splice_point + 100*eps()))          < 10000*eps()
# Testing in left spline territory.
abs(evaluate(spliced, splice_point-0.5) - evaluate(left_spline, splice_point-0.5))   < 100*eps()
abs(evaluate(spliced, splice_point-0.5) - evaluate(right_spline, splice_point-0.5))  > 0.1
# Testing solidly into right spline territory.
abs(evaluate(spliced, splice_point+0.5) - evaluate(left_spline, splice_point+0.5))   > 100*eps()
abs(evaluate(spliced, splice_point+0.5) - evaluate(right_spline, splice_point+0.5))  < 0.1

splice_point = 1.7
spliced = splice_splines(left_spline, right_spline, splice_point)
# Testing at splice point
abs(evaluate(spliced, splice_point) - evaluate(right_spline, splice_point))          < 100*eps()
# As this was NOT a continuous split
abs(evaluate(spliced, splice_point- 100*eps()) - evaluate(spliced, splice_point + 100*eps()))  > 0.1
# Testing in left spline territory.
abs(evaluate(spliced, splice_point-0.5) - evaluate(left_spline, splice_point-0.5))   < 100*eps()
abs(evaluate(spliced, splice_point-0.5) - evaluate(right_spline, splice_point-0.5))  > 0.1
# Testing solidly into right spline territory.
abs(evaluate(spliced, splice_point+0.5) - evaluate(left_spline, splice_point+0.5))   > 100*eps()
abs(evaluate(spliced, splice_point+0.5) - evaluate(right_spline, splice_point+0.5))  < 0.1


text:
function aetherlang_eval!(tokens::Vector{Token}, namespace::Namespace, line::Ref{UInt})::AetherlangObject
    """Evaluates a line of code (as tokens), returns the result as AetherlangObject{T}"""
    # one token evaluation
    if length(tokens) == 1
        if typeof(tokens[1]) != String
            return tokens[1]
        end
        token::String = tokens[1]
        # ritual/endrit evaluation
        if token == "ritual"
            ritname::String = "RIT#"*string(line[])
            endrit::UInt = line[]
            while !("endrit" == lines[endrit][1])
                endrit += 1
                if (endrit > length(lines))
                    throw(AetherlangError("Line `endrit` not found"))
                end
            end
            push!(rituals, Ritual(line[]+1, UInt(endrit-1), setindex!(copy(namespace), AetherlangObject(ritname), "dimname"), ritname))
            line[] = endrit
        elseif token[1] == '(' && token[end] == ')'
            # if one token and parentheses, remove them and re-eval
            return aetherlang_eval!(token[2:end-1], namespace, line)
        elseif token[1] == ':'
            return AetherlangObject(token[2:end], true)
        elseif token[1] == '\"' || isdigit(token[1]) || (token[1] == '-' && length(token)>1 && isdigit(token[2]))
            # try evaluating as a literal
            return AetherlangObject(eval(Meta.parse(token)))
        else
            # get from the namespace
            try
                return length(token)>7&&"aether."==token[1:7] ? ether[token[8:end]] : namespace[token]
            catch
                throw(AetherlangError("Name `$token` is not defined"))
            end
        end
    # A = X expression evaluation
    elseif tokens[2] == "="
        o::AetherlangObject = AetherlangObject()
        if length(tokens[1])>7&&"aether."==tokens[1][1:7]
            # aether.NAME = EXP
            o = ether[tokens[1][8:end]] = aetherlang_eval!(tokens[3:end], namespace, line)
            if o.type[] == AETH_BUILTIN_TYPES.timepoint
                throw(AetherlangError("Definition of timepoints in aether is forbidden"))
            end
        else
            # NAME = EXP
            if haskey(namespace, tokens[1])
                throw(AetherlangError("Trying to redefine a local dimensional constant `$(tokens[1])` of type "*namespace[tokens[1]].type[]))
            end
            o = namespace[tokens[1]] = aetherlang_eval!(tokens[3:end], namespace, line)
        end
        return o
    # dim X/enddim X evaluation
    elseif length(tokens) == 2 && tokens[1] == "dim"
        enddim::UInt = line[]
        while !("enddim" == lines[enddim][1] && tokens[2] == lines[enddim][2])
            enddim += 1
            if (enddim > length(lines))
                throw(AetherlangError("Line `enddim $(tokens[2])` not found"))
            end
        end
        push!(dimensions, Dimension(line[]+1, UInt(enddim-1), setindex!(copy(namespace), AetherlangObject(tokens[2]), "dimname"), tokens[2]))
        line[] = enddim
    # declaring a function
    elseif tokens[1] == "function"
        # try
        fname::String = tokens[2]
        eqsign::UInt = findfirst(t->t=="=", tokens)
        body::Vector{Union{String, AetherlangObject}} = Union{String, AetherlangObject}[]
        for t in tokens[eqsign+1:end]
            try
                push!(body, aetherlang_eval!(t, namespace, Ref{UInt}(0)))
            catch
                push!(body, t)
            end
        end
        f::AetherlangObject = AetherlangObject(AetherlangFunction(tokens[3:eqsign-1], body))
        if length(fname)>7&&"aether."==fname[1:7]
            ether[fname[8:end]] = f
        else
            if haskey(namespace, fname)
                throw(AetherlangError("Trying to redefine a local dimensional constant `$fname` of type "*namespace[fname].type[]))
            end
            namespace[fname] = f
        end
        return f
    # calling a function evaluation
    else
        callable::AetherlangObject = aetherlang_eval!(tokens[1], namespace, line)
        if callable.type[] != AETH_BUILTIN_TYPES.callable
            throw(AetherlangError("Object `$(tokens[1])` is not of type callable"))
        end
        return aetherlang_call(callable.dataref, namespace, Vector{AetherlangObject}([aetherlang_eval!(t, namespace, line) for t in tokens[2:end]]), line)
    end
    AetherlangObject()
end

function aetherlang_eval!(token::AetherlangObject, namespace::Namespace, line::Ref{UInt})::AetherlangObject
    """Evaluates an existing token, returns the same AetherlangObject as given"""
    token
end

function aetherlang_eval!(code::Vector{String}, namespace::Namespace, line::Ref{UInt})::AetherlangObject
    """Evaluates a line of code as string tokens, returns the result as AetherlangObject{T}"""
    aetherlang_eval!(Token[code...], namespace, line)
end

function aetherlang_eval!(code::String, namespace::Namespace, line::Ref{UInt})::AetherlangObject
    """Evaluates a line of code, returns the result as AetherlangObject{T}"""
    aetherlang_eval!(aetherlang_tokenize(code), namespace, line)
end

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
overwrite_configs = Dict()

configs = Dict(
    :substeps    => 8,                 # This controls how many steps will occur for each CESM coupling. Example: ocean couple to atmosphere every 24 hours but itself steps every 3 hours. This means we would expect `Δt` = 86400, and we set `substeps` = 8.
    :daily_record              => [],
    :monthly_record            => :ESSENTIAL,
    :enable_archive            => true,
    :archive_list              => "archive_list.txt",
    :rpointer_file             => "rpointer.hoom",
    :timeout                   => 60.0 * 20,
)

#=
configs = Dict(
    :casename    => "casename",
    :substeps    => 8,                 # This controls how many steps will occur for each CESM coupling. Example: ocean couple to atmosphere every 24 hours but itself steps every 3 hours. This means we would expect `Δt` = 86400, and we set `substeps` = 8.
    :caseroot                  => pwd(),
    :caserun                   => pwd(),
    :domain_file               => "/home/tienyiah/cesm_inputdata/cesm1/share/domains/domain.ocn.gx3v7.120323.nc",
    :daily_record              => [],
    :monthly_record            => :ALL,
    :enable_archive            => true,
    :archive_root              => pwd(),
    :archive_list              => "archive_list.txt",
    :rpointer_file             => "rpointer.hoom",
    :timeout                   => 60.0 * 20,
)
=#

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--config"
            help = "Configuration file."
            arg_type = String

        "--core"
            help = "Core of the model."
            arg_type = String

    end

    return parse_args(s)
end

parsed_args = parse_commandline()
#for (arg,val) in parsed_args
#    println("$arg  =>  ", repr(val))
#end

config_file = parsed_args["config"]
if config_file != nothing
    config_file = normpath( (isabspath(config_file)) ? config_file : joinpath(pwd(), config_file) )
    println("===== Load config file: ", config_file, " =====")
    include(config_file)
end

if !(@isdefined OMMODULE)
    core_name = parsed_args["core"]
    if core_name == nothing
        throw(ErrorException("Core ocean module is not provided. Please set --core option, or define OMMODULE in configuration file."))
    else
        module_name = "CESMCORE_" * core_name
        println("BASE: " * @__FILE__)
        core_file = normpath(joinpath(dirname(@__FILE__), "..", "models", core_name, module_name * ".jl"))
        println("Selected core: ", core_name, " => ", core_file )

        include(core_file)
        OMMODULE = getfield(Main, Symbol(module_name))
    end
end

println("===== Defining variables BEGIN =====")
for (k, v) in overwrite_configs
    if k in keys(configs)
        println("Overwrite config ", k, "...")
    else
        println("Add config ", k, "...")
    end

    configs[k] = v
end

if ! ( "tmp_folder" in keys(configs) )

    configs[:tmp_folder] = joinpath(configs[:caserun], "x_tmp")

end

print(json(configs, 4))
println("===== Defining variables END =====")

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
""" Julia serialization format layout

Formats a `Event` object in its Julia serialized form.
"""
type SerializedLayout <: Layout
end

SerializedLayout(conf::Dict) = SerializedLayout()

# Interface implementation
header(lyt::SerializedLayout) = UInt8[]

footer(lyt::SerializedLayout) = UInt8[]

contenttype(lyt::SerializedLayout) = "application/octet-stream"

function serialize(lyt::SerializedLayout, evnt::Event)
    iob = IOBuffer()
    serialize(io, evnt)
    return takebuf_array(iob)
end


text:
mutable struct zmp_data_t <: LCMType
    timestamp::Int64
    A::SMatrix{4, 4, Float64, 16}
    B::SMatrix{4, 2, Float64, 8}
    C::SMatrix{2, 4, Float64, 8}
    D::SMatrix{2, 2, Float64, 4}
    x0::SMatrix{4, 1, Float64, 4}
    y0::SMatrix{2, 1, Float64, 2}
    u0::SMatrix{2, 1, Float64, 2}
    R::SMatrix{2, 2, Float64, 4}
    Qy::SMatrix{2, 2, Float64, 4}
    S::SMatrix{4, 4, Float64, 16}
    s1::SMatrix{4, 1, Float64, 4}
    s1dot::SMatrix{4, 1, Float64, 4}
    s2::Float64
    s2dot::Float64

    com::SMatrix{4, 1, Float64, 4}
end

@lcmtypesetup(zmp_data_t)