text:
using OrdinaryDiffEq, Plots

function ball!(du,u,p,t)
    du[1] = u[2]
    du[2] = 0.0
    du[3] = u[4]
    du[4] = -p[1]
end

ground_condition(u,t,integrator) = u[3]
ground_affect!(integrator) = integrator.u[4] = -integrator.p[2] * integrator.u[4]
ground_cb = ContinuousCallback(ground_condition, ground_affect!)

u0 = [0.0,2.0,50.0,0.0]
tspan = (0.0,50.0)
p = [9.807, 0.9]

prob = ODEProblem(ball!,u0,tspan,p)
sol = solve(prob,Tsit5(),callback=ground_cb)
plot(sol, vars=(1,3), label = nothing, xlabel="x", ylabel="y")

stop_condition(u,t,integrator) = u[1] - 25.0
stop_cb = ContinuousCallback(stop_condition, terminate!)
cbs = CallbackSet(ground_cb, stop_cb)

tspan = (0.0, 1500.0)
prob = ODEProblem(ball!,u0,tspan,p)
sol = solve(prob,Tsit5(),callback=cbs)

rectangle(xc, yc, w, h) = Shape(xc .+ [-w,w,w,-w]./2.0, yc .+ [-h,-h,h,h]./2.0)

begin
    plot(sol, vars=(1,3), label=nothing, lw = 3, c=:black)
    xlabel!("x [m]")
    ylabel!("y [m]")
    plot!(rectangle(27.5, 25, 5, 50), c=:red, label = nothing)
    scatter!([25],[25],marker=:star, ms=10, label = nothing,c=:green)
    ylims!(0.0,50.0)
end

using Distributions

cor_dist = truncated(Normal(0.9, 0.02), 0.9-3*0.02, 1.0)
trajectories = 100000

prob_func(prob,i,repeat) = remake(prob, p = [p[1], rand(cor_dist)])
ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
ensemblesol = solve(ensemble_prob,Tsit5(),EnsembleThreads(),trajectories=trajectories, callback=cbs)

begin # plot
    plot(ensemblesol, vars = (1,3), lw=1,alpha=0.2, label=nothing, idxs = 1:350)
    xlabel!("x [m]")
    ylabel!("y [m]")
    plot!(rectangle(27.5, 25, 5, 50), c=:red, label = nothing)
    scatter!([25],[25],marker=:star, ms=10, label = nothing, c=:green)
    plot!(sol, vars=(1,3), label=nothing, lw = 3, c=:black, ls=:dash)
    xlims!(0.0,27.5)
end

obs(sol) = abs2(sol[3,end]-25)

mean_ensemble = mean([obs(sol) for sol in ensemblesol])

using DiffEqUncertainty

p_uncertain = [9.807, cor_dist]
expectation(obs, prob, u0, p_uncertain, Koopman(), Tsit5();
            ireltol = 1e-5, callback=cbs)

using NLopt, DiffEqSensitivity, ForwardDiff

make_u0(θ) = [θ[1],θ[2],θ[3], 0.0]

function 𝔼_loss(θ)   # \bbE
    u0 = make_u0(θ)
    expectation(obs, prob, u0, p_uncertain, Koopman(), Tsit5();
                 ireltol = 1e-5, callback=cbs)[1]
end

function 𝔼_loss_nlopt(x,∇)
    length(∇) > 0 ? ForwardDiff.gradient!(∇, 𝔼_loss,x) : nothing
    𝔼_loss(x)
end

opt = Opt(:LD_MMA, 3)
opt.lower_bounds = [-100.0,1.0, 10.0]
opt.upper_bounds = [0.0,3.0, 50.0]
opt.xtol_rel = 1e-3
opt.min_objective = 𝔼_loss_nlopt
(minf,minx,ret) = NLopt.optimize(opt, [-1.0, 2.0, 50.0])

ensembleprob = EnsembleProblem(remake(prob,u0 = make_u0(minx)),prob_func=prob_func)
ensemblesol = solve(ensembleprob,Tsit5(),EnsembleThreads(), trajectories=100_000, callback=cbs)

begin
    plot(ensemblesol, vars = (1,3), lw=1,alpha=0.1, label=nothing, idxs = 1:350)
    plot!(solve(remake(prob, u0=make_u0(minx)),Tsit5(), callback=cbs),
            vars=(1,3),label = nothing, c=:black, lw=3,ls=:dash)
    xlabel!("x [m]")
    ylabel!("y [m]")
    plot!(rectangle(27.5, 25, 5, 50), c=:red, label = nothing)
    scatter!([25],[25],marker=:star, ms=10, label = nothing,c=:green)
    ylims!(0.0,50.0)
    xlims!(minx[1], 27.5)
end

using BenchmarkTools

@btime NLopt.optimize($opt, $[-1.0, 2.0, 50.0])

constraint = [20.0, 25.0]
begin
    plot(rectangle(27.5, 25, 5, 50), c=:red, label = nothing)
    xlabel!("x [m]")
    ylabel!("y [m]")
    plot!([constraint[1], constraint[1]],[0.0,constraint[2]], lw=5, c=:black, label=nothing)
    scatter!([25],[25],marker=:star, ms=10, label = nothing,c=:green)
    ylims!(0.0,50.0)
    xlims!(minx[1], 27.5)
end

constraint_condition(u,t,integrator) = u[1] - constraint[1]
constraint_affect!(integrator) = integrator.u[3] < constraint[2] ? terminate!(integrator) : nothing
constraint_cb = ContinuousCallback(constraint_condition, constraint_affect!, save_positions=(true,false));
constraint_cbs = CallbackSet(ground_cb, stop_cb, constraint_cb)

ensemblesol = solve(ensembleprob,Tsit5(),EnsembleThreads(), trajectories=350, callback=constraint_cbs, maxstep=0.1)

begin
    plot(ensemblesol, vars = (1,3), lw=1,alpha=0.1, label=nothing)
    plot!(solve(remake(prob, u0=make_u0(minx)),Tsit5(), callback=constraint_cbs),
            vars=(1,3),label = nothing, c=:black, lw=3, ls=:dash)

    xlabel!("x [m]")
    ylabel!("y [m]")
    plot!(rectangle(27.5, 25, 5, 50), c=:red, label = nothing)
    plot!([constraint[1], constraint[1]],[0.0,constraint[2]], lw=5, c=:black)
    scatter!([25],[25],marker=:star, ms=10, label = nothing,c=:green)
    ylims!(0.0,50.0)
    xlims!(minx[1], 27.5)
end

constraint_obs(sol) = sol[1,end] ≈ constraint[1] ? one(sol[1,end]) : zero(sol[1,end])

expectation(constraint_obs, prob, make_u0(minx), p_uncertain, Koopman(), Tsit5();
            ireltol= 1e-9, iabstol = 1e-9, callback=constraint_cbs)[1]

function 𝔼_constraint(θ)
    u0 = [θ[1],θ[2],θ[3], 0.0]
    expectation(constraint_obs, prob, u0, p_uncertain, Koopman(), Tsit5(),
                ireltol= 1e-9, iabstol = 1e-9,callback=constraint_cbs)[1]
end

function 𝔼_constraint_nlopt(x,∇)
    length(∇) > 0 ? ForwardDiff.gradient!(∇, 𝔼_constraint,x) : nothing
    𝔼_constraint(x) - 0.01
end

opt = Opt(:LD_MMA, 3)
opt.lower_bounds = [-100.0, 1.0, 10.0]
opt.upper_bounds = [0.0, 3.0, 50.0]
opt.xtol_rel = 1e-3
opt.min_objective = 𝔼_loss_nlopt
inequality_constraint!(opt,𝔼_constraint_nlopt, 1e-5)
(minf2,minx2,ret2) = NLopt.optimize(opt, [-1.0, 2.0, 50.0])

λ = 𝔼_constraint(minx2)

λ - 0.01 <= 1e-5

ensembleprob = EnsembleProblem(remake(prob,u0 = make_u0(minx2)),prob_func=prob_func)
ensemblesol = solve(ensembleprob,Tsit5(),EnsembleThreads(),
                    trajectories=350, callback=constraint_cbs)

begin
    plot(ensemblesol, vars = (1,3), lw=1,alpha=0.1, label=nothing)
    plot!(solve(remake(prob, u0=make_u0(minx2)),Tsit5(), callback=constraint_cbs),
            vars=(1,3),label = nothing, c=:black, lw=3, ls=:dash)
    plot!([constraint[1], constraint[1]],[0.0,constraint[2]], lw=5, c=:black)

    xlabel!("x [m]")
    ylabel!("y [m]")
    plot!(rectangle(27.5, 25, 5, 50), c=:red, label = nothing)
    scatter!([25],[25],marker=:star, ms=10, label = nothing,c=:green)
    ylims!(0.0,50.0)
    xlims!(minx[1], 27.5)
end

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
"""
$(TYPEDEF)

Adjacency describing edges per grid cell
"""
abstract type CellEdges  <: AbstractGridAdjacency end

"""
$(TYPEDEF)

Adjacency describing cells per grid edge
"""
abstract type EdgeCells  <: AbstractGridAdjacency end

"""
$(TYPEDEF)

Adjacency describing nodes per grid edge
"""
abstract type EdgeNodes <: AbstractGridAdjacency end

"""
$(TYPEDEF)

Adjacency describing cells per boundary or interior face
"""
abstract type BFaceCells <: AbstractGridAdjacency end

"""
$(TYPEDEF)

Adjacency describing outer normals to boundary faces
"""
abstract type BFaceNormals <: AbstractGridComponent end

"""
$(SIGNATURES)

Prepare edge adjacencies (celledges, edgecells, edgenodes)

Currently depends on ExtendableSparse, we may want to remove this
adjacency.
"""
function prepare_edges!(grid::ExtendableGrid)
    Ti=eltype(grid[CellNodes])
    cellnodes=grid[CellNodes]
    geom=grid[CellGeometries][1]
    # Create cell-node incidence matrix
    cellnode_adj=asparse(cellnodes)

    # Create node-node incidence matrix for neigboring
    # nodes.
    nodenode_adj=cellnode_adj*transpose(cellnode_adj)

    # To get unique edges, we set the lower triangular part
    # including the diagonal to 0
    for icol=1:length(nodenode_adj.colptr)-1
        for irow=nodenode_adj.colptr[icol]:nodenode_adj.colptr[icol+1]-1
            if nodenode_adj.rowval[irow]>=icol
                nodenode_adj.nzval[irow]=0
            end
        end
    end
    dropzeros!(nodenode_adj)

    # Now we know the number of edges and
    nedges=length(nodenode_adj.nzval)

    if dim_space(grid)==2
        # Let us do the Euler test (assuming no holes in the domain)
        v=num_nodes(grid)
        e=nedges
        f=num_cells(grid)+1
        @assert v-e+f==2
    end

    if dim_space(grid)==1
        @assert nedges==num_cells(grid)
    end

    # Calculate edge nodes and celledges
    edgenodes=zeros(Ti,2,nedges)
    celledges=zeros(Ti,num_edges(geom),num_cells(grid))
    cen=local_celledgenodes(geom)

    for icell=1:num_cells(grid)
        for iedge=1:num_edges(geom)
            n1=cellnodes[cen[1,iedge],icell]
            n2=cellnodes[cen[2,iedge],icell]

            # We need to look in nodenod_adj for upper triangular part entries
            # therefore, we need to swap accordingly before looking
	    if (n1<n2)
		n0=n1
		n1=n2
		n2=n0;
	    end

            for irow=nodenode_adj.colptr[n1]:nodenode_adj.colptr[n1+1]-1
                if nodenode_adj.rowval[irow]==n2
                    # If the coresponding entry has been found, set its
                    # value. Note that this introduces a different edge orientation
                    # compared to the one found locally from cell data
                    celledges[iedge,icell]=irow
                    edgenodes[1,irow]=n1
                    edgenodes[2,irow]=n2
                end
            end
        end
    end

    # Create sparse incidence matrix for the cell-edge adjacency
    celledge_adj=asparse(celledges)

    # The edge cell matrix is the transpose
    edgecell_adj=SparseMatrixCSC(transpose(celledge_adj))

    # Get the adjaency array from the matrix
    edgecells=zeros(Ti,2,nedges) ## for 3D we need more here!
    for icol=1:length(edgecell_adj.colptr)-1
        ii=1
        for irow=edgecell_adj.colptr[icol]:edgecell_adj.colptr[icol+1]-1
            edgecells[ii,icol]=edgecell_adj.rowval[irow]
            ii+=1
        end
    end

    grid[EdgeCells]=edgecells
    grid[CellEdges]=celledges
    grid[EdgeNodes]=edgenodes
    true
end

ExtendableGrids.instantiate(grid, ::Type{CellEdges})=prepare_edges!(grid) && grid[CellEdges]
ExtendableGrids.instantiate(grid, ::Type{EdgeCells})=prepare_edges!(grid) && grid[EdgeCells]
ExtendableGrids.instantiate(grid, ::Type{EdgeNodes})=prepare_edges!(grid) && grid[EdgeNodes]

function prepare_bfacecells!(grid)
    cn=grid[CellNodes]
    bn=grid[BFaceNodes]
    dim=dim_space(grid)
    abc=asparse(cn)'*asparse(bn)
    abcx=dropzeros!(SparseMatrixCSC(abc.m,abc.n, abc.colptr, abc.rowval, map( i-> i==dim,  abc.nzval)))
    grid[BFaceCells]=VariableTargetAdjacency(abcx)
    true
end

ExtendableGrids.instantiate(grid, ::Type{BFaceCells})=prepare_bfacecells!(grid) && grid[BFaceCells]

normal!(normal,nodes,coord,::Type{Val{1}}) = normal[1]=1.0

function normal!(normal,nodes,coord,::Type{Val{2}})
    normal[1] = -(coord[2,nodes[1]]-coord[2,nodes[2]])
    normal[2] =   coord[1,nodes[1]]-coord[1,nodes[2]]
    d=norm(normal)
    normal[1]/=d
    normal[2]/=d
end

function normal!(normal,nodes,coord,::Type{Val{3}})
    ax=coord[1,nodes[1]]-coord[1,nodes[2]]
    ay=coord[2,nodes[1]]-coord[2,nodes[2]]
    az=coord[3,nodes[1]]-coord[3,nodes[2]]

    bx=coord[1,nodes[1]]-coord[1,nodes[3]]
    by=coord[2,nodes[1]]-coord[2,nodes[3]]
    bz=coord[3,nodes[1]]-coord[3,nodes[3]]

    normal[1]=(ay*bz-by*az)
    normal[2]=(az*bx-bz*ax)
    normal[3]=(ax*by-bx*ay)

    d=norm(normal)
    normal[1]/=d
    normal[2]/=d
    normal[3]/=d
end

function midpoint!(mid,nodes,coord)
    dim=size(coord,1)
    nn=size(nodes,1)
    for i=1:dim
        mid[i]=sum(coord[i,nodes])/nn
    end
end

function adjust!(normal,cmid,bmid)
    d=dot(normal,bmid-cmid)
    if d<0.0
        normal.*= -1
    end
end

function prepare_bfacenormals!(grid)
    bfnodes=grid[BFaceNodes]
    nbf=size(bfnodes,2)
    cellnodes=grid[CellNodes]
    dim=dim_space(grid)
    bfcells=grid[BFaceCells]
    bfnormals=zeros(dim,nbf)
    coord=grid[Coordinates]
    cmid=zeros(dim)
    bmid=zeros(dim)
    for ibf=1:nbf
        icell=bfcells[1,ibf]
        @views normal!(bfnormals[:,ibf],bfnodes[:,ibf],coord,Val{dim})
        @views midpoint!(cmid,cellnodes[:,icell],coord)
        @views midpoint!(bmid,bfnodes[:,ibf],coord)
        @views adjust!(bfnormals[:,ibf],cmid,bmid)
    end
    grid[BFaceNormals]=bfnormals
    true
end

ExtendableGrids.instantiate(grid, ::Type{BFaceNormals})=prepare_bfacenormals!(grid) && grid[BFaceNormals]

# The following methods are uses in VoronoiFVM.
# Do we need a more systematic approach here ?
"""
$(SIGNATURES)

Number of edges in grid.
"""
num_edges(grid::ExtendableGrid)=haskey(grid,EdgeNodes) ?  num_sources(grid[EdgeNodes]) : 0

"""
$(SIGNATURES)

Number of nodes of 0D vertex
"""
num_nodes(::Type{Vertex0D})=1

"""
$(SIGNATURES)

Number of edges of 0D vertex
"""
num_edges(::Type{Vertex0D})=0

const cen_Edge1D=reshape([1 2],:,1)
"""
$(SIGNATURES)

Cell-edege node numbering for 1D edge
"""
local_celledgenodes(::Type{Edge1D})=cen_Edge1D

"""
$(SIGNATURES)

Number of nodes for 1D edge
"""
num_nodes(::Type{Edge1D})=2

"""
$(SIGNATURES)

Number of edges of 1D edge
"""
num_edges(::Type{Edge1D})=1

const cen_Triangle2D=[ 2 3 1; 3 1 2]

"""
$(SIGNATURES)

Cell-edege node numbering for 2D triangle
"""
local_celledgenodes(::Type{Triangle2D})=cen_Triangle2D

const cen_Tetrahedron3D=[ 3 4 2  1 1 1; 4 2 3  4 3 2]

"""
$(SIGNATURES)

Cell-edege node numbering for 2D triangle
"""
local_celledgenodes(::Type{Tetrahedron3D})=cen_Tetrahedron3D

"""
$(SIGNATURES)

Number of nodes in 2D triangle
"""
num_nodes(::Type{Triangle2D})=3

"""
$(SIGNATURES)

Number of nodes in 3D tetrahedron
"""
num_nodes(::Type{Tetrahedron3D})=4

"""
$(SIGNATURES)

Number of edges in 2D triangle
"""
num_edges(::Type{Triangle2D})=3

"""
$(SIGNATURES)

Number of edges in 3D tetrahedron
"""
num_edges(::Type{Tetrahedron3D})=6

ext: jl


text:
r = HTTP.request("GET", "http://$localIP:$port/api/runIBox?res=oxy115829.dat&IBox_port=8888&IBox_path=$pathToIBox&IBox_host=$localIP&config=IBOpen&resName=111&arg=-open")
r = HTTP.request("GET", "http://$localIP:$port/api/getDataTree")
tree = JSON.parse(String(r.body))[1]

r = HTTP.request("GET", "http://127.0.0.1:8888/api/getData?dataName=QPoint&all")
QPoint = reinterpret(Int32, base64decode(r.body)) |> collect

r = HTTP.request("GET", "http://127.0.0.1:8888/api/getData?dataName=ClassQRS&all")
Class = reinterpret(Int32, base64decode(r.body)) |> collect
unCl = unique(Class)
QPoint = QPoint[1:length(Class)]

r = HTTP.request("GET", "http://127.0.0.1:8888/api/getData?dataName=SubClassQRS&all")
SubClass = reinterpret(Int16, base64decode(r.body)) |> collect
unSubCl = unique(SubClass)

treeClasses = Dict{String,Any}()
treeClasses["nodes"] = []
treeClasses["name"] = "allClasses"
treeClasses["size"] = length(Class)

for i = 1:length(unCl)
    isCl = Class.==unCl[i]
    node =  Dict{String,Any}()
    node["nodes"] = []
    node["name"] = unCl[i]
    node["size"] = sum(isCl)
    node["QPoint"] = QPoint[isCl]
    for j = 1:length(unSubCl) #теперь добавляем подклассы
        isSubCl= (SubClass.== unSubCl[j]) .& isCl
        subnode =  Dict{String,Any}()
        subnode["nodes"] = []
        subnode["name"] = unSubCl[j]
        subnode["size"] = sum(isSubCl)
        subnode["QPoint"] = QPoint[isSubCl]
        push!(node["nodes"],subnode)
    end
    push!(treeClasses["nodes"],node)
end

result = JSON.json(treeClasses)


text:
module GLAbstraction

using ModernGL
using FixedPointNumbers

import FileIO: load, save

import FixedPointNumbers: N0f8, N0f16, N0f8, Normed

import Base: merge, resize!,  similar, length, getindex, setindex!

include("context.jl")

include("AbstractGPUArray.jl")

#Methods which get overloaded by GLExtendedFunctions.jl:
import ModernGL.glShaderSource
import ModernGL.glGetAttachedShaders
import ModernGL.glGetActiveUniform
import ModernGL.glGetActiveAttrib
import ModernGL.glGetProgramiv
import ModernGL.glGetIntegerv
import ModernGL.glGenBuffers
import ModernGL.glGetProgramiv
import ModernGL.glGenVertexArrays
import ModernGL.glGenTextures
import ModernGL.glGenFramebuffers
import ModernGL.glGetTexLevelParameteriv
import ModernGL.glGenRenderbuffers
import ModernGL.glDeleteTextures
import ModernGL.glDeleteVertexArrays
import ModernGL.glDeleteBuffers
import ModernGL.glGetShaderiv
import ModernGL.glViewport
import ModernGL.glScissor

include("utils.jl")
include("buffer.jl")
include("texture.jl")
include("framebuffer.jl")
include("vertexarray.jl")
include("uniformbuffer.jl")

include("shader/uniforms.jl")
include("shader/shader.jl")
include("shader/program.jl")
include("shader/glsl_typenames.jl")
include("conversions.jl")
include("for_moderngl.jl")
dir(dirs...) = joinpath(dirname(@__FILE__), "..", dirs...)
end # module


text:
module GenieDiscordBot

using Discord

c = Client("")

end


text:
module SceneObjects
  using CG;

  export SceneObject, Sphere, AABB, Scene;

  abstract SceneObject;

  type Sphere <: SceneObject
    center::Vec4f;
    radius::Float32;
  end

  type AABB <: SceneObject
    center::Vec4f;
    # positive half length from center to face of box
    hx::Float32;
    hy::Float32;
    hz::Float32;
  end

  type Scene
    sceneObjects::Vector{SceneObject};
  end
end

max_stars_repo_name: fpoppinga/cg

max_issues_repo_name: fpoppinga/cg

max_forks_repo_name: fpoppinga/cg


text:
#SignalProcessing: Time/Frequency domain signals
#-------------------------------------------------------------------------------

import MDDatasets: _operators1, _operators2, _basefn1

#==Helper functions
===============================================================================#
missing(s::Symbol)=throw(ArgumentError("Missing argument: $s"))

Base.isvalid(::DS{:time}, d::DataTF) = d.tvalid
Base.isvalid(::DS{:freq}, d::DataTF) = d.fvalid
#invalidate(::DS{:time}, d::DataTF) = (d.tvalid = false)
#invalidate(::DS{:freq}, d::DataTF) = (d.fvalid = false)
#Implement validate(DataTF)?

isteven(d::DataTF) = evenlength(d.t)
isteven(d::DataFreq) = isteven(d.data)
isteven(d::DataTime) = isteven(d.data)

#==Basic DataTime/DataFreq functionality
===============================================================================#
Base.copy(d::DataTF) = DataTF(d.tperiodic, d.tvalid, d.fvalid, t, copy(d.xt), f, copy(d.Xf))
Base.copy(d::DataTime) = validatelengths(DataTime(copy(d.data)))
Base.copy(d::DataFreq) = validatelengths(DataFreq(copy(d.data)))

Base.length(d::DataTime) = length(validatelengths(d).data.t)
Base.length(d::DataFreq) = length(validatelengths(d).data.f)

#==Accessors/Data converters
===============================================================================#
DataF1(d::DataTime) = DataF1(collect(d.data.t), copy(d.data.xt))
DataF1(d::DataFreq) = DataF1(collect(d.data.f), copy(d.data.Xf))

datavec(::DS{:time}, d::DataTF) = d.t
datavec(::DS{:time}, d::DataTime) = datavec(TIME, d.data)
datavec(::DS{:time}, d::DataFreq) = datavec(TIME, d.data)

datavec(::DS{:freq}, d::DataTF) = d.f
datavec(::DS{:freq}, d::DataTime) = datavec(FREQ, d.data)
datavec(::DS{:freq}, d::DataFreq) = datavec(FREQ, d.data)

datavec(::DS{:sig}, d::DataTime) = d.data.xt
datavec(::DS{:sig}, d::DataFreq) = d.data.Xf

#Returns sampled frequency spectrum (non-periodic signals)
#Normalizes FFT results -> |X(f)| samples or {ak} coefficients
function fspectrum(d::DataFreq)
	ensure(!d.data.tperiodic,
		ArgumentError("Frequency spectrum is only valid for non-periodic signals"))
	return DataF1(collect(d.data.f), d.data.Xf)
end

#Returns Fourier-series coefficients (t-periodic signals)
#Normalizes FFT results -> |X(f)| samples or {ak} coefficients
function fcoeff(d::DataFreq)
	ensure(d.data.tperiodic,
		ArgumentError("Fourier-series coefficients are only valid for time-periodic signals"))
	return DataF1(collect(d.data.f), d.data.Xf./length(d.data.t))
end

#==Support basic math operations
===============================================================================#

for op in _operators1; @eval begin #CODEGEN--------------------------------------

#op DataTime:
Base.$op(d::DataTime) = DataTime(d.data, ($op).(d.data.xt))

#op DataFreq:
Base.$op(d::DataFreq) = DataFreq(d.data, ($op).(d.data.Xf))

end; end #CODEGEN---------------------------------------------------------------

for op in _operators2; @eval begin #CODEGEN--------------------------------------

#DataTime op Number:
Base.$op(d::DataTime, n::Number) = DataTime(d.data, ($op).(d.data.xt, n))

#Number op DataTime:
Base.$op(n::Number, d::DataTime) = DataTime(d.data, ($op).(n, d.data.xt))

#DataFreq op Number:
Base.$op(d::DataFreq, n::Number) = DataFreq(d.data, ($op).(d.data.Xf, n))

#Number op DataFreq:
Base.$op(n::Number, d::DataFreq) = DataFreq(d.data, ($op).(n, d.data.Xf))

end; end #CODEGEN---------------------------------------------------------------

#==Support for 1-argument functions from Base
===============================================================================#

for fn in _basefn1; @eval begin #CODEGEN----------------------------------------

#fn(DataTime)
$fn(d::DataTime) = DataTime(d.data, ($fn).(d.data.xt))

#fn(DataFreq)
$fn(d::DataFreq) = DataFreq(d.data, ($fn).(d.data.Xf))

end; end #CODEGEN---------------------------------------------------------------

#==Fourier transforms/series
===============================================================================#
#TODO: use \mscrF (ℱ)??
#function Fourier(x::DataF1)

function timedomain(d::DataFreq)
	if !d.data.tvalid
		d.data.xt = irfft(d.data.Xf, length(d.data.t))
		d.data.tvalid = true
	end
	return DataTime(d.data) #Reference same data
end

function freqdomain(d::DataTime)
	if !d.data.fvalid
		d.data.Xf = rfft(d.data.xt)
		d.data.fvalid = true
	end
	return DataFreq(d.data) #Reference same data
end

#Last line


text:
struct Stardate
  sd::Float64
end

function Stardate(zdt::ZonedDateTime)
  tx = astimezone(zdt, tz"UTC")
  y0 = year(tx)
  t0 = Stardates.getstartofyear(y0).zoneddatetime
  t1 = Stardates.getstartofyear(y0 + 1).zoneddatetime
  Stardate(y0 + (tx - t0) / (t1 - t0))
end

function Stardate(dt::Dates.DateTime)
  # Assume dt is in UTC
  y0 = Dates.year(dt)
  t0 = Stardates.getstartofyear(y0).unix
  t1 = Stardates.getstartofyear(y0 + 1).unix
  Stardate(y0 + (Dates.datetime2unix(dt) - t0) / (t1 - t0))
end

function Stardate(d::Dates.Date, hr, mi, ss, tz)
  Stardate(ZonedDateTime(
    Dates.year(d),
    Dates.month(d),
    Dates.day(d),
    hr,
    mi,
    ss,
    tz,
  ))
end

function stardate(args...)
  Stardate(args...).sd
end

"Return the stardate of the closing bell of the date"
function nyse_closing_stardate(d::Dates.Date)
  stardate(d, 16, 0, 0, tz"America/New_York")
end

"Return the stardate of the closing bell of the date represented by iso8601"
function nyse_closing_stardate(iso8601)
  nyse_closing_stardate(parse(Dates.Date, iso8601))
end

function Stardate(st::Base.StatStruct)
  Stardate(Dates.unix2datetime(st.mtime))
end

function Stardate()
  Stardate(now(tz"UTC"))
end

function canonical(sd::Stardate)
  Printf.@sprintf("%0.15f", sd.sd)
end

function short(sd::Stardate)
  Printf.@sprintf("%0.3f", sd.sd)
end


text:
# Continuous Moving Bed Reconstruction using Gridding
#
# All data and results are in single-precision (32-bit) float, to
# reduce memory requirements.

using NFFT
using MAT
include("correct.jl")
require("Julia Recon Scripts/coords.jl")
require("Julia Recon Scripts/pnfft.jl")
require("Julia Recon Scripts/coilcomb.jl")

data_path = "../data_input/";
base_name = ARGS[1]
infile = string(data_path,base_name ,".mat")
out_dir = "../data_output/";

println(ARGS)

# ADJUSTABLE PARAMETERS
golden = true # this should correspond to the sequence
npe_per_z = int(ARGS[4])
sliceΔz = int(float(ARGS[3])/(20 * float(ARGS[2])))

nvirt=2  # number of virtual coils to compress to before the recon
fcc_phase = 1.4580  # from SIN file
fcc_scale = 1.1784

# load some COMBI data
println("reading data from $infile")
mat = matread(infile)
data = mat["Data"]
#info = mat["info"]  # in case we want to pull params from INFO
#noise = info["FRC_NOISE_DATA"][1]
#npe = size(data,3)  # reduce data for testing
#data = data[:,:,int(0.4npe):int(0.5npe)]

function recon_combi_ms2d(data, npe_per_z, sliceΔz; kernwidth=3,
                          oversamp=1.5f0, golden=false)

  nchan, nro, npe = size(data)
  #n = int(nro / 2)  # reconstructed image dimension
   n = int(nro)  # reconstructed image dimension
  pstart = [1:sliceΔz:npe-npe_per_z]
  pend = pstart + npe_per_z - 1
  nz = length(pstart)

  pcenter = int(pstart+npe_per_z/2-1)
  matwrite(string(out_dir,"pcenter.mat"),["pcenter"=>pcenter])

  # combine QBC channels data with Fast Channel Combination
  if nchan == 2
    println("performing FCC with scale=$fcc_scale and phase $fcc_phase")
    data = data[1,:,:] + fcc_scale*exp(-1im*fcc_phase)*data[2,:,:]
  else
    #prewhiten!(data, noise)
    data = coilcompress(data; nvirt=4)
  end

  # other corrections?
  tic()
  data, slope, intercepts = unshift(data; golden=golden)
  toc()

  # Generate coordinates for all profiles up front
  # TODO figure out why I need offset = 1 here
  println("generating radial coordinates")
  K = radialcoords(nro, npe; golden=golden, offset=0)

  # Generate the sample density compensation
  # Should be the same for each slice as long as they have
  # equal numbers of profiles.
  #println("computing SDC")
  #weights = radialsdc(nro, npe_per_z; golden=golden)

  println("reconstructing ... $(size(data))")
  fComb = pnfft_combi(data, K, pstart, pend, kernwidth, oversamp; verbose=true)

  return fComb
end

tic()
fComb = recon_combi_ms2d(data, npe_per_z, sliceΔz; golden=golden)
toc()

#println("saving to $outfile")
#matwrite(outfile, ["data"=>fComb])
matwrite(string(out_dir,base_name,"_JULIA_RECON.mat"),["Image"=>fComb])

#using PyPlot

#nx, ny, nz = size(fComb)
#figure(1)
#clf()
#imshow(abs(squeeze(fComb[:,:,int(0.3nz)],3)), interpolation="nearest", cmap="gray")
#title("axial")
#colorbar()

#figure(2)
#clf()
#aspect_ratio = (0.4/1.5)*nz/nx
#imshow(abs(squeeze(fComb[:,110,:],2)),aspect=aspect_ratio, cmap="gray", interpolation="nearest")
#title("sag or cor")
#colorbar()

#figure(7)
#clf()
#plot(p.windowLUT[1])
#title("kernel")

max_stars_repo_path: code/Julia Recon Scripts/recon.jl

max_issues_repo_path: code/Julia Recon Scripts/recon.jl

max_forks_repo_path: code/Julia Recon Scripts/recon.jl


text:
using ConfParser, Compat
using Compat.Test

conf = ConfParse("confs/config.ini")
parse_conf!(conf)
loc = ConfParse("confs/local.ini")
parse_conf!(loc)

@test retrieve(conf, "database", "user") == "dbuser"
@test retrieve(conf, "database", "password") == "abc123"
@test retrieve(loc, "database", "password") == "xyz789"
@test retrieve(loc, "foobar", "foo") == ["foo","bar"]
@test merge!(conf, loc) == true
save!(conf, "confs/out.conf")

# validate new conf write
conf = ConfParse("confs/out.conf")
parse_conf!(conf)
@test retrieve(conf, "database", "user") == "dbuser"
@test retrieve(conf, "database", "password") == "xyz789"
@test retrieve(conf, "foobar", "foo") == ["foo","bar"]


text:
using Documenter, Turing
using Plots
using LibGit2: clone
using Weave

# Get paths.
markdown_path = joinpath(@__DIR__, "markdown")

function polish_latex(path::String)
    txt = open(f -> read(f, String), path)
    open(path, "w+") do f
        txt = replace(txt, raw"$$" => raw"\$\$")

        beg_align = r"\\\$\\\$\n?\\begin{align}"
        end_align = r"\\end{align}\n?\\\$\\\$?"

        txt = replace(txt, beg_align => "\$\$\n\\begin{align}")
        txt = replace(txt, end_align => "\\end{align}\n\$\$")

        # Make inline math work.
        inline = r"(?<=[^\\])\$(?!\$)"
        txt = replace(txt, inline => raw"$$")

        # If there's more than three dollar signs together, use only two.
        extra_dollars = r"\${3,}"
        txt = replace(txt, extra_dollars => raw"$$")

        write(f, txt)
    end
end

function add_yaml(path)
    lines = readlines(path, keep=true)

    title = ""
    permalink = "/:collection/:name/"

    # Find the ID sections.
    for i = 1:length(lines)
        line = lines[i]

        if startswith(line, "# ")
            title = replace(line, "# " => "")
            title = replace(title, "\n" => "")

            lines[i] = ""
            break
        end
    end

    yaml_lines = ["---\n", "title: $title\n", "permalink: $permalink\n", "---\n"]

    # Write lines back.
    open(path, "w+") do f
        # Prepend YAML.
        for line in yaml_lines
            write(f, line)
        end

        # Add remaining content.
        for line in lines
            write(f, replace(line, "![](figure" => "![](/tutorials/figure"))
        end
    end
end

# Weave all examples
try
    for file in readdir(@__DIR__)
        if endswith(file, ".ipynb") || endswith(file, ".jmd")
            out_name = split(file, ".")[1] * ".md"
            out_path = joinpath(markdown_path, out_name)

            full_path = joinpath(@__DIR__, file)

            if mtime(out_path) < mtime(full_path)
                @warn "Weaving $full_path as it has been updated since the least weave."
                Weave.weave(
                    full_path,
                    doctype = "github",
                    out_path = out_path,
                    mod = Main,
                    throw_errors = true
                )

                polish_latex(out_path)
                add_yaml(out_path)
            else
                @warn "Skipping $full_path as it has not been updated."
            end
        end
    end
catch e
    println("Weaving error: $e")
    rethrow(e)
end

max_stars_repo_name: sivapvarma/TuringTutorials

max_issues_repo_name: sivapvarma/TuringTutorials

max_forks_repo_name: sivapvarma/TuringTutorials


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestMPS

import MathOptInterface
using Test

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MPS = MOI.FileFormats.MPS
const MPS_TEST_FILE = "test.mps"

function _test_model_equality(model_string, variables, constraints; kwargs...)
    model = MPS.Model(; kwargs...)
    MOIU.loadfromstring!(model, model_string)
    MOI.write_to_file(model, MPS_TEST_FILE)
    model_2 = MPS.Model()
    MOI.read_from_file(model_2, MPS_TEST_FILE)
    return MOI.Test.util_test_models_equal(
        model,
        model_2,
        variables,
        constraints,
    )
end

function test_show()
    @test sprint(show, MPS.Model()) ==
          "A Mathematical Programming System (MPS) model"
end

function test_quadratic()
    model = MPS.Model()
    @test_throws(
        MOI.UnsupportedAttribute,
        MOIU.loadfromstring!(
            model,
            """
variables: x
minobjective: 1.0*x*x
""",
        )
    )
end

function test_nonempty()
    model = MPS.Model()
    @test MOI.is_empty(model)
    MOI.add_variable(model)
    @test !MOI.is_empty(model)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    MOI.add_variable(model)
    @test_throws Exception MOI.read_from_file(
        model,
        joinpath(@__DIR__, "failing_models", "bad_name.mps"),
    )
end

function test_failing_models()
    @testset "$(filename)" for filename in filter(
        f -> endswith(f, ".mps"),
        readdir(joinpath(@__DIR__, "failing_models")),
    )
        @test_throws Exception MOI.read_from_file(
            MPS.Model(),
            joinpath(@__DIR__, "failing_models", filename),
        )
    end
end

function test_empty_row_name()
    model = MPS.Model()
    x = MOI.add_variable(model)
    MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
        MOI.LessThan(1.0),
    )
    @test_throws Exception sprint(MPS.write_rows, model)
end

function test_sos()
    model = MPS.Model()
    x = MOI.add_variables(model, 3)
    names = Dict{MOI.VariableIndex,String}()
    for i in 1:3
        MOI.set(model, MOI.VariableName(), x[i], "x$(i)")
        names[x[i]] = "x$(i)"
    end
    MOI.add_constraint(
        model,
        MOI.VectorOfVariables(x),
        MOI.SOS1([1.5, 2.5, 3.5]),
    )
    MOI.add_constraint(
        model,
        MOI.VectorOfVariables(x),
        MOI.SOS2([1.25, 2.25, 3.25]),
    )
    @test sprint(MPS.write_sos, model, names) ==
          "SOS\n" *
          " S1 SOS1\n" *
          "    x1        1.5\n" *
          "    x2        2.5\n" *
          "    x3        3.5\n" *
          " S2 SOS2\n" *
          "    x1        1.25\n" *
          "    x2        2.25\n" *
          "    x3        3.25\n"
end

function test_maximization()
    model = MPS.Model()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    @test sprint(MPS.write_columns, model, true, ["x"], Dict(x => "x")) ==
          "COLUMNS\n    x         OBJ       -1\n"
end

function test_maximization_objsense_false()
    model = MPS.Model(; print_objsense = true)
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    sprint(MPS.write, model)
    @test sprint(MPS.write_columns, model, false, ["x"], Dict(x => "x")) ==
          "COLUMNS\n    x         OBJ       1\n"
end

function test_stacked_data()
    model = MPS.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "stacked_data.mps"))
    model_2 = MPS.Model()
    MOIU.loadfromstring!(
        model_2,
        """
variables: x, y, z
maxobjective: x + y + z + 2.5
con1: 1.0 * x in Interval(1.0, 5.0)
con2: 1.0 * x in Interval(2.0, 6.0)
con3: 1.0 * x in Interval(3.0, 7.0)
con4: 2.0 * x in Interval(4.0, 8.0)
y in Integer()
y in Interval(1.0, 4.0)
z in ZeroOne()
""",
    )
    MOI.set(model_2, MOI.Name(), "stacked_data")
    return MOI.Test.util_test_models_equal(
        model,
        model_2,
        ["x", "y", "z"],
        ["con1", "con2", "con3", "con4"],
        [
            ("y", MOI.Integer()),
            ("y", MOI.Interval{Float64}(1.0, 4.0)),
            ("z", MOI.ZeroOne()),
        ],
    )
end

function test_free_integer()
    model = MPS.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "free_integer.mps"))
    model_2 = MPS.Model()
    MOIU.loadfromstring!(
        model_2,
        """
variables: x
minobjective: x
con1: 1.0 * x >= 1.0
x in Integer()
""",
    )
    return MOI.Test.util_test_models_equal(
        model,
        model_2,
        ["x"],
        ["con1"],
        [("x", MOI.Integer())],
    )
end

function test_min_objective()
    return _test_model_equality(
        """
    variables: x
    minobjective: x
""",
        ["x"],
        String[],
    )
end

function test_objconst()
    return _test_model_equality(
        """
    variables: x
    minobjective: 1.1 * x + 1.2
""",
        ["x"],
        String[],
    )
end

function test_default_rhs_greater()
    return _test_model_equality(
        """
variables: x
minobjective: x
c1: 2.0 * x >= 0.0
""",
        ["x"],
        ["c1"],
    )
end

function test_default_rhs_less()
    return _test_model_equality(
        """
    variables: x
    minobjective: x
    c1: 2.0 * x <= 0.0
""",
        ["x"],
        ["c1"],
    )
end

function test_default_rhs_equal()
    return _test_model_equality(
        """
variables: x
minobjective: x
c1: 2.0 * x == 0.0
""",
        ["x"],
        ["c1"],
    )
end

function test_min_scalaraffine()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x
""",
        ["x"],
        String[],
    )
end

function test_scalaraffine_greaterthan()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x
c1: 1.1 * x >= 2.0
""",
        ["x"],
        ["c1"],
    )
end

function test_scalaraffine_lessthan()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x
c1: 1.1 * x <= 2.0
""",
        ["x"],
        ["c1"],
    )
end

function test_scalaraffine_equalto()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x
c1: 1.1 * x == 2.0
""",
        ["x"],
        ["c1"],
    )
end

function test_scalaraffine_interval()
    return _test_model_equality(
        """
variables: x
minobjective: 1.2x
c1: 1.1 * x in Interval(1.0, 2.0)
""",
        ["x"],
        ["c1"],
    )
end

function test_objsense_max()
    _test_model_equality(
        """
variables: x
maxobjective: 1.2x
c1: 1.0 * x >= 0.0
""",
        ["x"],
        ["c1"],
        print_objsense = true,
    )
    _test_model_equality(
        """
variables: x
minobjective: 1.2x
c1: 1.0 * x >= 0.0
""",
        ["x"],
        ["c1"],
        print_objsense = true,
    )
    return
end

function test_MARKER_INT()
    model = MPS.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: x, y, z
minobjective: x + y + z
x in Integer()
c2: 2 * x + -1.0 * z <= 1.0
z in ZeroOne()
x >= 1.0
""",
    )
    MOI.write_to_file(model, MPS_TEST_FILE)
    model_2 = MPS.Model()
    MOI.read_from_file(model_2, MPS_TEST_FILE)
    return MOI.Test.util_test_models_equal(
        model,
        model_2,
        ["x", "y", "z"],
        ["c2"],
        [
            ("x", MOI.Integer()),
            ("z", MOI.ZeroOne()),
            ("x", MOI.GreaterThan(1.0)),
        ],
    )
end

function test_zero_variable_bounds()
    model = MPS.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: x, y, z
minobjective: x + y + z
x >= 0.0
y <= 0.0
""",
    )
    MOI.write_to_file(model, MPS_TEST_FILE)
    model_2 = MPS.Model()
    MOI.read_from_file(model_2, MPS_TEST_FILE)
    return MOI.Test.util_test_models_equal(
        model,
        model_2,
        ["x", "y", "z"],
        String[],
        [
            ("x", MOI.GreaterThan{Float64}(0.0)),
            ("y", MOI.LessThan{Float64}(0.0)),
        ],
    )
end

function test_nonzero_variable_bounds()
    model = MPS.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: w, x, y, z
minobjective: w + x + y + z
x == 1.0
y >= 2.0
z <= 3.0
w in Interval(4.0, 5.0)
""",
    )
    MOI.write_to_file(model, MPS_TEST_FILE)
    model_2 = MPS.Model()
    MOI.read_from_file(model_2, MPS_TEST_FILE)
    return MOI.Test.util_test_models_equal(
        model,
        model_2,
        ["w", "x", "y", "z"],
        String[],
        [
            ("x", MOI.EqualTo{Float64}(1.0)),
            ("y", MOI.GreaterThan{Float64}(2.0)),
            ("z", MOI.LessThan{Float64}(3.0)),
            ("w", MOI.Interval{Float64}(4.0, 5.0)),
        ],
    )
end

function test_multiple_variable_bounds()
    model = MPS.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: a_really_long_name
minobjective: a_really_long_name
a_really_long_name >= 1.0
a_really_long_name <= 2.0
""",
    )
    MOI.write_to_file(model, MPS_TEST_FILE)
    @test read(MPS_TEST_FILE, String) ==
          "NAME          \n" *
          "ROWS\n" *
          " N  OBJ\n" *
          "COLUMNS\n" *
          "    a_really_long_name OBJ       1\n" *
          "RHS\n" *
          "RANGES\n" *
          "BOUNDS\n" *
          " LO bounds    a_really_long_name 1\n" *
          " UP bounds    a_really_long_name 2\n" *
          "ENDATA\n"
end

function test_unused_variable()
    # In this test, `x` will be written to the file with a 0 objective
    # coefficient since it does not appear in the objective or in the
    # constraints.
    model = MPS.Model()
    MOIU.loadfromstring!(
        model,
        """
variables: x, y
minobjective: y
c1: 2.0 * y >= 1.0
x >= 0.0
""",
    )
    MOI.write_to_file(model, MPS_TEST_FILE)
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    model2 = MPS.Model()
    MOI.read_from_file(model2, MPS_TEST_FILE)
    @test MOI.get(model2, MOI.NumberOfVariables()) == 2
end

function test_names_with_spaces()
    model = MPS.Model()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x[1, 2]")
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
        MOI.EqualTo(1.0),
    )
    MOI.set(model, MOI.ConstraintName(), c, "c c")
    @test sprint(write, model) ==
          "NAME          \n" *
          "ROWS\n" *
          " N  OBJ\n" *
          " E  c_c\n" *
          "COLUMNS\n" *
          "    x[1,_2]   c_c       1\n" *
          "RHS\n" *
          "    rhs       c_c       1\n" *
          "RANGES\n" *
          "BOUNDS\n" *
          " FR bounds    x[1,_2]\n" *
          "ENDATA\n"
end

function test_objsense_default()
    model = MPS.Model()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    @test sprint(write, model) ==
          "NAME          \n" *
          "ROWS\n" *
          " N  OBJ\n" *
          "COLUMNS\n" *
          "    x         OBJ       -1\n" *
          "RHS\n" *
          "RANGES\n" *
          "BOUNDS\n" *
          " FR bounds    x\n" *
          "ENDATA\n"
end

function test_objsense_true()
    model = MPS.Model(; print_objsense = true)
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), x)
    @test sprint(write, model) ==
          "NAME          \n" *
          "OBJSENSE MAX\n" *
          "ROWS\n" *
          " N  OBJ\n" *
          "COLUMNS\n" *
          "    x         OBJ       1\n" *
          "RHS\n" *
          "RANGES\n" *
          "BOUNDS\n" *
          " FR bounds    x\n" *
          "ENDATA\n"
end

function test_sos_constraints()
    model = MPS.Model()
    x = MOI.add_variables(model, 3)
    MOI.set.(model, MOI.VariableName(), x, ["x1", "x2", "x3"])
    MOI.add_constraint(
        model,
        MOI.VectorOfVariables(x),
        MOI.SOS1([1.0, 2.0, 3.0]),
    )
    MOI.add_constraint(
        model,
        MOI.VectorOfVariables(reverse(x)),
        MOI.SOS2([1.2, 2.3, 3.4]),
    )
    @test sprint(write, model) ==
          "NAME          \n" *
          "ROWS\n" *
          " N  OBJ\n" *
          "COLUMNS\n" *
          "    x1        OBJ       0\n" *
          "    x2        OBJ       0\n" *
          "    x3        OBJ       0\n" *
          "RHS\n" *
          "RANGES\n" *
          "BOUNDS\n" *
          " FR bounds    x1\n" *
          " FR bounds    x2\n" *
          " FR bounds    x3\n" *
          "SOS\n" *
          " S1 SOS1\n" *
          "    x1        1\n" *
          "    x2        2\n" *
          "    x3        3\n" *
          " S2 SOS2\n" *
          "    x3        1.2\n" *
          "    x2        2.3\n" *
          "    x1        3.4\n" *
          "ENDATA\n"
    io = IOBuffer()
    write(io, model)
    seekstart(io)
    new_model = MPS.Model()
    Base.read!(io, new_model)
    constraints = MOI.get(new_model, MOI.ListOfConstraintTypesPresent())
    @test (MOI.VectorOfVariables, MOI.SOS1{Float64}) in constraints
    @test (MOI.VectorOfVariables, MOI.SOS2{Float64}) in constraints
    return
end

function test_generic_names()
    model = MPS.Model(; generic_names = true)
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [x, y]), 0.0),
        MOI.EqualTo(1.0),
    )
    MOI.add_constraint(model, y, MOI.GreaterThan(2.0))
    @test sprint(write, model) ==
          "NAME          \n" *
          "ROWS\n" *
          " N  OBJ\n" *
          " E  R1\n" *
          "COLUMNS\n" *
          "    C1        R1        1\n" *
          "    C2        R1        1\n" *
          "RHS\n" *
          "    rhs       R1        1\n" *
          "RANGES\n" *
          "BOUNDS\n" *
          " FR bounds    C1\n" *
          " LO bounds    C2        2\n" *
          " PL bounds    C2\n" *
          "ENDATA\n"
end

function test_rew_filename()
    model = MOI.FileFormats.Model(; filename = "test.rew")
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [x, y]), 0.0),
        MOI.EqualTo(1.0),
    )
    MOI.add_constraint(model, y, MOI.GreaterThan(2.0))
    @test sprint(write, model) ==
          "NAME          \n" *
          "ROWS\n" *
          " N  OBJ\n" *
          " E  R1\n" *
          "COLUMNS\n" *
          "    C1        R1        1\n" *
          "    C2        R1        1\n" *
          "RHS\n" *
          "    rhs       R1        1\n" *
          "RANGES\n" *
          "BOUNDS\n" *
          " FR bounds    C1\n" *
          " LO bounds    C2        2\n" *
          " PL bounds    C2\n" *
          "ENDATA\n"
end

function test_rew_format()
    model = MOI.FileFormats.Model(; format = MOI.FileFormats.FORMAT_REW)
    x = MOI.add_variable(model)
    y = MOI.add_variable(model)
    c = MOI.add_constraint(
        model,
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, [x, y]), 0.0),
        MOI.EqualTo(1.0),
    )
    MOI.add_constraint(model, y, MOI.GreaterThan(2.0))
    @test sprint(write, model) ==
          "NAME          \n" *
          "ROWS\n" *
          " N  OBJ\n" *
          " E  R1\n" *
          "COLUMNS\n" *
          "    C1        R1        1\n" *
          "    C2        R1        1\n" *
          "RHS\n" *
          "    rhs       R1        1\n" *
          "RANGES\n" *
          "BOUNDS\n" *
          " FR bounds    C1\n" *
          " LO bounds    C2        2\n" *
          " PL bounds    C2\n" *
          "ENDATA\n"
end

function test_infinite_interval()
    model = MPS.Model()
    x = MOI.add_variable(model)
    MOI.add_constraint(model, 1.0 * x, MOI.Interval(-Inf, Inf))
    MOI.add_constraint(model, 1.0 * x, MOI.Interval(-Inf, 1.0))
    MOI.add_constraint(model, 1.0 * x, MOI.Interval(2.0, Inf))
    MOI.add_constraint(model, 1.0 * x, MOI.Interval(3.0, 4.0))
    @test sprint(write, model) ==
          "NAME          \n" *
          "ROWS\n" *
          " N  OBJ\n" *
          " N  c1\n" *
          " L  c2\n" *
          " G  c3\n" *
          " L  c4\n" *
          "COLUMNS\n" *
          "    x1        c1        1\n" *
          "    x1        c2        1\n" *
          "    x1        c3        1\n" *
          "    x1        c4        1\n" *
          "RHS\n" *
          "    rhs       c2        1\n" *
          "    rhs       c3        2\n" *
          "    rhs       c4        4\n" *
          "RANGES\n" *
          "    rhs       c4        1\n" *
          "BOUNDS\n" *
          " FR bounds    x1\n" *
          "ENDATA\n"
    return
end

function runtests()
    for name in names(@__MODULE__, all = true)
        if startswith("$(name)", "test_")
            @testset "$name" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    sleep(1.0)  # Allow time for unlink to happen.
    rm(MPS_TEST_FILE, force = true)
    return
end

end

TestMPS.runtests()


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    ExecutableProduct(prefix, "readstat", :readstat),
    LibraryProduct(prefix, String["libreadstat"], :libreadstat),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/davidanthoff/ReadStatBuilder/releases/download/v0.1.1-build.3"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/ReadStat.v0.1.1-dev.aarch64-linux-gnu.tar.gz", "3a670d3ebc6e7eb2ed4f228142adb1ab975f2d0fa9489037453e098c9a3717a5"),
    Linux(:aarch64, :musl) => ("$bin_prefix/ReadStat.v0.1.1-dev.aarch64-linux-musl.tar.gz", "03bcd9b7402739ba674881c9a9af731aa3ee8644d2050a07554c4a3ab3320ebe"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/ReadStat.v0.1.1-dev.arm-linux-gnueabihf.tar.gz", "7255a92bf6b320499061a0478705df3960bb3cd39db259e53602e7d18cfa51e2"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/ReadStat.v0.1.1-dev.arm-linux-musleabihf.tar.gz", "6a2688f0dfbc1fd8912b494642e194d09a778356989f25793f8baaf1313e3b15"),
    Linux(:i686, :glibc) => ("$bin_prefix/ReadStat.v0.1.1-dev.i686-linux-gnu.tar.gz", "c969c3cda37aa9924b98abe86ddcc72dc6c9b59744d2983e930b8d4b604f4de6"),
    Linux(:i686, :musl) => ("$bin_prefix/ReadStat.v0.1.1-dev.i686-linux-musl.tar.gz", "89d42858955bdeed036c07856b8ebdd87d73bc60ad771e246d67195454c31319"),
    Windows(:i686) => ("$bin_prefix/ReadStat.v0.1.1-dev.i686-w64-mingw32.tar.gz", "9d4919bf033d67c72a7c2c10d91b3a021cddf1d33db64db449fd63c17e73c14c"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/ReadStat.v0.1.1-dev.powerpc64le-linux-gnu.tar.gz", "b7ffcdbc6eab6cfb673f581403508722982ef935f556c25fc53b7102c138eae8"),
    MacOS(:x86_64) => ("$bin_prefix/ReadStat.v0.1.1-dev.x86_64-apple-darwin14.tar.gz", "653c06b1f5be334510bf29588c7c0a8c34ad9dcab60693858cbba05f3083944a"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/ReadStat.v0.1.1-dev.x86_64-linux-gnu.tar.gz", "87b2b96614f61ecdc63af890d3ecd926374bb7403b084c1fedc5cfd14db5a6ce"),
    Linux(:x86_64, :musl) => ("$bin_prefix/ReadStat.v0.1.1-dev.x86_64-linux-musl.tar.gz", "b95c666cdef269c6994bc009c7d09ff21b688c395e7578bbe8f6d723aa919342"),
    FreeBSD(:x86_64) => ("$bin_prefix/ReadStat.v0.1.1-dev.x86_64-unknown-freebsd11.1.tar.gz", "b41746a4845db5a4b1732fcc8f5c13b464e425c4a0ab44161c76c26dd41507bf"),
    Windows(:x86_64) => ("$bin_prefix/ReadStat.v0.1.1-dev.x86_64-w64-mingw32.tar.gz", "576faea230119c7ea6590a802c4695961d88a7f3ea515ecc5e01e740fbb2354c"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
module ChampionMasteryV4 # LOLTools

using ..LOLTools: AbstractDTO, lol_api_server, http_action, call_api, nothing_in_event
using HTTP

struct ChampionMasteryDTO <: AbstractDTO
    chestGranted
    championLevel
    championPoints
    championId
    championPointsUntilNextLevel
    lastPlayTime
    tokensEarned
    championPointsSinceLastLevel
    summonerId
    ChampionMasteryDTO(chestGranted, championLevel, championPoints, championId, championPointsUntilNextLevel, lastPlayTime, tokensEarned, championPointsSinceLastLevel, summonerId) = new(chestGranted, championLevel, championPoints, championId, championPointsUntilNextLevel, lastPlayTime, tokensEarned, championPointsSinceLastLevel, summonerId)
end

"""
         by_summoner_id(api_key::String,
                        platform::String,
                        encryptedSummonerId::String ;
                        endpoint::HTTP.URI = lol_api_server(platform),
                        action::Function = http_action,
                        event::Function = nothing_in_event)::Vector{ChampionMasteryDTO}
"""
function by_summoner_id(api_key::String,
                        platform::String,
                        encryptedSummonerId::String ;
                        endpoint::HTTP.URI = lol_api_server(platform),
                        action::Function = http_action,
                        event::Function = nothing_in_event)::Vector{ChampionMasteryDTO}
    call_api(Vector{ChampionMasteryDTO}, api_key, action, endpoint, "/lol/champion-mastery/v4/champion-masteries/by-summoner/$encryptedSummonerId", event, by_summoner_id)
end

"""
         by_summoner_id(api_key::String,
                        platform::String,
                        encryptedSummonerId::String,
                        championId::Int64 ;
                        endpoint::HTTP.URI = lol_api_server(platform),
                        action::Function = http_action,
                        event::Function = nothing_in_event)::ChampionMasteryDTO
"""
function by_summoner_id(api_key::String,
                        platform::String,
                        encryptedSummonerId::String,
                        championId::Int64 ;
                        endpoint::HTTP.URI = lol_api_server(platform),
                        action::Function = http_action,
                        event::Function = nothing_in_event)::ChampionMasteryDTO
    call_api(ChampionMasteryDTO, api_key, action, endpoint, "/lol/champion-mastery/v4/champion-masteries/by-summoner/$encryptedSummonerId/by-champion/$championId", event, by_summoner_id)
end

end # module LOLTools.ChampionMasteryV4


text:
#                                  _|
#    _|_|_|  _|    _|    _|_|_|  _|_|_|_|    _|_|    _|_|_|  _|_|
#  _|_|      _|    _|  _|_|        _|      _|_|_|_|  _|    _|    _|
#      _|_|  _|    _|      _|_|    _|      _|        _|    _|    _|
#  _|_|_|      _|_|_|  _|_|_|        _|_|    _|_|_|  _|    _|    _|
#                  _|
#              _|_|

@with_kw struct ThreeBodyMasses
    m1::Float64
    m2::Float64
    m3::Float64
    m0::Float64
    ThreeBodyMasses(m1,m2,m3,m0) =
        m0<m1+m2+m3 ? error("m₀ should be bigger than m₁+m₂+m₃") : new(m1,m2,m3,m0)
end

lims(k::Int, ms::ThreeBodyMasses) = k==1 ? lims1(ms) : ((k==2) ? lims2(ms) : lims3(ms))
lims1(ms::ThreeBodyMasses) = ((ms.m2+ms.m3)^2, (ms.m0-ms.m1)^2)
lims2(ms::ThreeBodyMasses) = ((ms.m3+ms.m1)^2, (ms.m0-ms.m2)^2)
lims3(ms::ThreeBodyMasses) = ((ms.m1+ms.m2)^2, (ms.m0-ms.m3)^2)
#
import Base: getindex, ^, length
^(ms::ThreeBodyMasses,i::Int) = (ms.m1, ms.m2, ms.m3, ms.m0).^i
#
function getindex(ms::ThreeBodyMasses, i::Int)
    (i==0 || i==4) && return ms.m0
    i==2 && return ms.m2
    i==3 && return ms.m3
    i!=1 && error("i should be equal to 0,1,2,3 or 4")
    return ms.m1
end
nt(ms::ThreeBodyMasses) = NamedTuple{(:m0,:m1,:m2,:m3)}([ms.m0,ms.m1,ms.m2,ms.m3])
#
import Base: iterate
#
iterate(ms::ThreeBodyMasses) = iterate((ms.m1,ms.m2,ms.m3,ms.m0))
iterate(ms::ThreeBodyMasses, state) = iterate((ms.m1,ms.m2,ms.m3,ms.m0),state)
#
ThreeBodyMasses(m1,m2,m3;
    m0=error("used the format ThreeBodyMasses(0.15,0.5,0.94; m0=4.5)")) =
    ThreeBodyMasses(m1=m1,m2=m2,m3=m3,m0=m0)
#
# -----------------------------------------------------

@with_kw struct ThreeBodySpins{T}
    two_h1::T
    two_h2::T
    two_h3::T
    two_h0::T
    ThreeBodySpins{T}(two_h1,two_h2,two_h3,two_h0) where {T} =
        isodd(two_h1+two_h2+two_h3+two_h0) ?
            error("baryon number is not conserved") :
            new(two_h1,two_h2,two_h3,two_h0)
end
function getindex(two_hs::ThreeBodySpins, i::Int)
    (i==0 || i==4) && return two_hs.two_h0
    i==2 && return two_hs.two_h2
    i==3 && return two_hs.two_h3
    i!=1 && error("i should be equal to 0,1,2,3 or 4")
    return two_hs.two_h1
end
length(σs::ThreeBodySpins) = 4
iterate(two_hs::ThreeBodySpins)        = iterate((two_hs.two_h1,two_hs.two_h2,two_hs.two_h3,two_hs.two_h0))
iterate(two_hs::ThreeBodySpins, state) = iterate((two_hs.two_h1,two_hs.two_h2,two_hs.two_h3,two_hs.two_h0),state)
#
import Base:getproperty
function getproperty(two_js::ThreeBodySpins, sym::Symbol)
    (sym in fieldnames(ThreeBodySpins)) && return getfield(two_js, sym)
    (sym == :two_j1 || sym == :two_λ1) && return getfield(two_js, :two_h1)
    (sym == :two_j2 || sym == :two_λ2) && return getfield(two_js, :two_h2)
    (sym == :two_j3 || sym == :two_λ3) && return getfield(two_js, :two_h3)
    (sym == :two_j0 || sym == :two_λ0) && return getfield(two_js, :two_h0)
    error("no property $(sym)")
end
#
ThreeBodySpins(two_h1,two_h2,two_h3;
    two_h0=error("used the format ThreeBodySpins(1,1,0; two_j0=2)")) =
    ThreeBodySpins(two_h1=two_h1,two_h2=two_h2,two_h3=two_h3,two_h0=two_h0)

#
#
@with_kw struct ThreeBodySystem{T,K}
    ms::T
    two_js::K = ThreeBodySpins(0,0,0,0)
end
#
# convenient constructors
ThreeBodySystem(ms::ThreeBodyMasses) = ThreeBodySystem(ms=ms)
ThreeBodySystem(m1,m2,m3; m0, two_js=ThreeBodySpins(0,0,0,0)) = ThreeBodySystem(ThreeBodyMasses(m1,m2,m3; m0=m0), two_js)
#
two_j0(tbs::ThreeBodySystem) = tbs.two_js[4]
two_j1(tbs::ThreeBodySystem) = tbs.two_js[1]
two_j2(tbs::ThreeBodySystem) = tbs.two_js[2]
two_j3(tbs::ThreeBodySystem) = tbs.two_js[3]

# -----------------------------------------------------
@with_kw struct ThreeBodyParities
    P1::Char
    P2::Char
    P3::Char
    P0::Char
end
function getindex(Ps::ThreeBodyParities, i::Int)
    (i==0 || i==4) && return Ps.P0
    i==2 && return Ps.P2
    i==3 && return Ps.P3
    i!=1 && error("i should be equal to 0,1,2,3 or 4")
    return Ps.P1
end
length(σs::ThreeBodyParities) = 4
iterate(Ps::ThreeBodyParities)        = iterate((Ps.P1,Ps.P2,Ps.P3,Ps.P0))
iterate(Ps::ThreeBodyParities, state) = iterate((Ps.P1,Ps.P2,Ps.P3,Ps.P0),state)
#
ThreeBodyParities(P1,P2,P3;
    P0=error("used the format ThreeBodyParities(1,1,0; two_j0=2)")) =
    ThreeBodyParities(P1=P1,P2=P2,P3=P3,P0=P0)

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

# Dynamic variables
@with_kw struct Invariants
    σ1::Float64
    σ2::Float64
    σ3::Float64
end
function Invariants(ms::ThreeBodyMasses;σ1=-1.0,σ2=-1.0,σ3=-1.0)
    sign(σ1)+sign(σ2)+sign(σ3)!=1 && error("the method works with TWO invariants given: $((σ1,σ2,σ3))")
    σ3 < 0 && return Invariants(;σ1=σ1,σ2=σ2,σ3=sum(ms^2)-σ1-σ2)
    σ1 < 0 && return Invariants(;σ2=σ2,σ3=σ3,σ1=sum(ms^2)-σ2-σ3)
    return Invariants(;σ3=σ3,σ1=σ1,σ2=sum(ms^2)-σ3-σ1)
end
#
iterate(σs::Invariants)        = iterate((σs.σ1,σs.σ2,σs.σ3))
iterate(σs::Invariants, state) = iterate((σs.σ1,σs.σ2,σs.σ3),state)
length(σs::Invariants) = 3
function getindex(σs::Invariants, i::Int)
    i==1 && return σs.σ1
    i==2 && return σs.σ2
    i!=3 && error("i should be equal to 1,2, or 3")
    return σs.σ3
end
nt(σs::Invariants) = NamedTuple{(:σ1,:σ2,:σ3)}([σs.σ1,σs.σ2,σs.σ3])

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

@with_kw struct DalitzPlotPoint{I,S}
    σs::I
    two_λs::S
end
#
function randomPoint(ms::ThreeBodyMasses)
    σ1 = lims1(ms)[1] + rand()* (lims1(ms)[2]-lims1(ms)[1])
    σ3 = σ3of1(2rand()-1, σ1, ms^2)
    return Invariants(ms;σ1=σ1,σ3=σ3);
end

function randomPoint(tbs::ThreeBodySystem)
    DalitzPlotPoint(σs=randomPoint(tbs.ms),
        two_λs=[rand(-j:2:j) for j in tbs.two_js])
end

#                                                                _|
#    _|_|_|    _|_|    _|_|_|      _|_|    _|  _|_|    _|_|_|  _|_|_|_|    _|_|
#  _|    _|  _|_|_|_|  _|    _|  _|_|_|_|  _|_|      _|    _|    _|      _|_|_|_|
#  _|    _|  _|        _|    _|  _|        _|        _|    _|    _|      _|
#    _|_|_|    _|_|_|  _|    _|    _|_|_|  _|          _|_|_|      _|_|    _|_|_|
#        _|
#    _|_|

function border(k, ms::ThreeBodyMasses; Nx::Int=300)
    (i,j) = ij_from_k(k)
    σiv = range(lims(i,ms)..., length=Nx)
    σkm = [σkofi(k,-1.0,σi,ms^2) for σi in σiv]
    σkp = [σkofi(k, 1.0,σi,ms^2) for σi in σiv]
    return (σiv, [σkm σkp])
end
#
border31(ms; Nx::Int=300) = border(3, ms; Nx=Nx)
border12(ms; Nx::Int=300) = border(1, ms; Nx=Nx)
border23(ms; Nx::Int=300) = border(2, ms; Nx=Nx)
#
function flatDalitzPlotSample(ms::ThreeBodyMasses; Nev::Int=10000, σbins::Int=500)
    @unpack m0,m1,m2,m3 = ms
    density = getbinned1dDensity(σ1->sqrt(Kallen(σ1,m2^2,m3^2)*Kallen(σ1,m0^2,m1^2))/σ1, lims1(ms), σbins)
    σ1v = [rand(density) for _ in 1:Nev]
    σ3v = [σ3of1(2*rand()-1,σ1,ms^2) for σ1 in σ1v]
    return [Invariants(ms; σ1=σ1, σ3=σ3) for (σ1,σ3) in zip(σ1v,σ3v)]
end

#
inrange(x,r) = r[1]<x<r[2]
inphrange(σs::Invariants, ms::ThreeBodyMasses) = Kibble(σs,ms^2) < 0 &&
    inrange(σs[1],lims1(ms)) && inrange(σs[2],lims2(ms)) && inrange(σs[3],lims3(ms))
#
#
change_basis_3from1(τ1, ms::ThreeBodyMasses) = change_basis_3from1(τ1..., ms.m1^2, ms.m2^2, ms.m3^2, ms.m0^2)
change_basis_1from2(τ2, ms::ThreeBodyMasses) = change_basis_3from1(τ2..., ms.m2^2, ms.m3^2, ms.m1^2, ms.m0^2)
change_basis_2from3(τ3, ms::ThreeBodyMasses) = change_basis_3from1(τ3..., ms.m3^2, ms.m1^2, ms.m2^2, ms.m0^2)

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
module DaprClients
using gRPCClient
using ProtoBuf:google

include("dapr.jl")
using .dapr.proto.runtime.v1

import Base: show

# begin service: dapr.proto.runtime.v1.Dapr

export DaprBlockingClient, DaprClient

struct DaprBlockingClient
    controller::gRPCController
    channel::gRPCChannel
    stub::DaprBlockingStub

    function DaprBlockingClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = DaprBlockingStub(channel)
        new(controller, channel, stub)
    end
end

struct DaprClient
    controller::gRPCController
    channel::gRPCChannel
    stub::DaprStub

    function DaprClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = DaprStub(channel)
        new(controller, channel, stub)
    end
end

show(io::IO, client::DaprBlockingClient) = print(io, "DaprBlockingClient(", client.channel.baseurl, ")")
show(io::IO, client::DaprClient) = print(io, "DaprClient(", client.channel.baseurl, ")")

import .dapr.proto.runtime.v1: InvokeService
"""
    InvokeService

- input: dapr.proto.runtime.v1.InvokeServiceRequest
- output: dapr.proto.common.v1.InvokeResponse
"""
InvokeService(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.InvokeServiceRequest) = InvokeService(client.stub, client.controller, inp)
InvokeService(client::DaprClient, inp::dapr.proto.runtime.v1.InvokeServiceRequest, done::Function) = InvokeService(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: GetState
"""
    GetState

- input: dapr.proto.runtime.v1.GetStateRequest
- output: dapr.proto.runtime.v1.GetStateResponse
"""
GetState(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.GetStateRequest) = GetState(client.stub, client.controller, inp)
GetState(client::DaprClient, inp::dapr.proto.runtime.v1.GetStateRequest, done::Function) = GetState(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: GetBulkState
"""
    GetBulkState

- input: dapr.proto.runtime.v1.GetBulkStateRequest
- output: dapr.proto.runtime.v1.GetBulkStateResponse
"""
GetBulkState(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.GetBulkStateRequest) = GetBulkState(client.stub, client.controller, inp)
GetBulkState(client::DaprClient, inp::dapr.proto.runtime.v1.GetBulkStateRequest, done::Function) = GetBulkState(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: SaveState
"""
    SaveState

- input: dapr.proto.runtime.v1.SaveStateRequest
- output: google.protobuf.Empty
"""
SaveState(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.SaveStateRequest) = SaveState(client.stub, client.controller, inp)
SaveState(client::DaprClient, inp::dapr.proto.runtime.v1.SaveStateRequest, done::Function) = SaveState(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: DeleteState
"""
    DeleteState

- input: dapr.proto.runtime.v1.DeleteStateRequest
- output: google.protobuf.Empty
"""
DeleteState(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.DeleteStateRequest) = DeleteState(client.stub, client.controller, inp)
DeleteState(client::DaprClient, inp::dapr.proto.runtime.v1.DeleteStateRequest, done::Function) = DeleteState(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: DeleteBulkState
"""
    DeleteBulkState

- input: dapr.proto.runtime.v1.DeleteBulkStateRequest
- output: google.protobuf.Empty
"""
DeleteBulkState(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.DeleteBulkStateRequest) = DeleteBulkState(client.stub, client.controller, inp)
DeleteBulkState(client::DaprClient, inp::dapr.proto.runtime.v1.DeleteBulkStateRequest, done::Function) = DeleteBulkState(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: ExecuteStateTransaction
"""
    ExecuteStateTransaction

- input: dapr.proto.runtime.v1.ExecuteStateTransactionRequest
- output: google.protobuf.Empty
"""
ExecuteStateTransaction(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.ExecuteStateTransactionRequest) = ExecuteStateTransaction(client.stub, client.controller, inp)
ExecuteStateTransaction(client::DaprClient, inp::dapr.proto.runtime.v1.ExecuteStateTransactionRequest, done::Function) = ExecuteStateTransaction(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: PublishEvent
"""
    PublishEvent

- input: dapr.proto.runtime.v1.PublishEventRequest
- output: google.protobuf.Empty
"""
PublishEvent(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.PublishEventRequest) = PublishEvent(client.stub, client.controller, inp)
PublishEvent(client::DaprClient, inp::dapr.proto.runtime.v1.PublishEventRequest, done::Function) = PublishEvent(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: InvokeBinding
"""
    InvokeBinding

- input: dapr.proto.runtime.v1.InvokeBindingRequest
- output: dapr.proto.runtime.v1.InvokeBindingResponse
"""
InvokeBinding(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.InvokeBindingRequest) = InvokeBinding(client.stub, client.controller, inp)
InvokeBinding(client::DaprClient, inp::dapr.proto.runtime.v1.InvokeBindingRequest, done::Function) = InvokeBinding(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: GetSecret
"""
    GetSecret

- input: dapr.proto.runtime.v1.GetSecretRequest
- output: dapr.proto.runtime.v1.GetSecretResponse
"""
GetSecret(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.GetSecretRequest) = GetSecret(client.stub, client.controller, inp)
GetSecret(client::DaprClient, inp::dapr.proto.runtime.v1.GetSecretRequest, done::Function) = GetSecret(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: GetBulkSecret
"""
    GetBulkSecret

- input: dapr.proto.runtime.v1.GetBulkSecretRequest
- output: dapr.proto.runtime.v1.GetBulkSecretResponse
"""
GetBulkSecret(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.GetBulkSecretRequest) = GetBulkSecret(client.stub, client.controller, inp)
GetBulkSecret(client::DaprClient, inp::dapr.proto.runtime.v1.GetBulkSecretRequest, done::Function) = GetBulkSecret(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: RegisterActorTimer
"""
    RegisterActorTimer

- input: dapr.proto.runtime.v1.RegisterActorTimerRequest
- output: google.protobuf.Empty
"""
RegisterActorTimer(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.RegisterActorTimerRequest) = RegisterActorTimer(client.stub, client.controller, inp)
RegisterActorTimer(client::DaprClient, inp::dapr.proto.runtime.v1.RegisterActorTimerRequest, done::Function) = RegisterActorTimer(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: UnregisterActorTimer
"""
    UnregisterActorTimer

- input: dapr.proto.runtime.v1.UnregisterActorTimerRequest
- output: google.protobuf.Empty
"""
UnregisterActorTimer(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.UnregisterActorTimerRequest) = UnregisterActorTimer(client.stub, client.controller, inp)
UnregisterActorTimer(client::DaprClient, inp::dapr.proto.runtime.v1.UnregisterActorTimerRequest, done::Function) = UnregisterActorTimer(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: RegisterActorReminder
"""
    RegisterActorReminder

- input: dapr.proto.runtime.v1.RegisterActorReminderRequest
- output: google.protobuf.Empty
"""
RegisterActorReminder(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.RegisterActorReminderRequest) = RegisterActorReminder(client.stub, client.controller, inp)
RegisterActorReminder(client::DaprClient, inp::dapr.proto.runtime.v1.RegisterActorReminderRequest, done::Function) = RegisterActorReminder(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: UnregisterActorReminder
"""
    UnregisterActorReminder

- input: dapr.proto.runtime.v1.UnregisterActorReminderRequest
- output: google.protobuf.Empty
"""
UnregisterActorReminder(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.UnregisterActorReminderRequest) = UnregisterActorReminder(client.stub, client.controller, inp)
UnregisterActorReminder(client::DaprClient, inp::dapr.proto.runtime.v1.UnregisterActorReminderRequest, done::Function) = UnregisterActorReminder(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: GetActorState
"""
    GetActorState

- input: dapr.proto.runtime.v1.GetActorStateRequest
- output: dapr.proto.runtime.v1.GetActorStateResponse
"""
GetActorState(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.GetActorStateRequest) = GetActorState(client.stub, client.controller, inp)
GetActorState(client::DaprClient, inp::dapr.proto.runtime.v1.GetActorStateRequest, done::Function) = GetActorState(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: ExecuteActorStateTransaction
"""
    ExecuteActorStateTransaction

- input: dapr.proto.runtime.v1.ExecuteActorStateTransactionRequest
- output: google.protobuf.Empty
"""
ExecuteActorStateTransaction(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.ExecuteActorStateTransactionRequest) = ExecuteActorStateTransaction(client.stub, client.controller, inp)
ExecuteActorStateTransaction(client::DaprClient, inp::dapr.proto.runtime.v1.ExecuteActorStateTransactionRequest, done::Function) = ExecuteActorStateTransaction(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: InvokeActor
"""
    InvokeActor

- input: dapr.proto.runtime.v1.InvokeActorRequest
- output: dapr.proto.runtime.v1.InvokeActorResponse
"""
InvokeActor(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.InvokeActorRequest) = InvokeActor(client.stub, client.controller, inp)
InvokeActor(client::DaprClient, inp::dapr.proto.runtime.v1.InvokeActorRequest, done::Function) = InvokeActor(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: GetMetadata
"""
    GetMetadata

- input: google.protobuf.Empty
- output: dapr.proto.runtime.v1.GetMetadataResponse
"""
GetMetadata(client::DaprBlockingClient, inp::google.protobuf.Empty) = GetMetadata(client.stub, client.controller, inp)
GetMetadata(client::DaprClient, inp::google.protobuf.Empty, done::Function) = GetMetadata(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: SetMetadata
"""
    SetMetadata

- input: dapr.proto.runtime.v1.SetMetadataRequest
- output: google.protobuf.Empty
"""
SetMetadata(client::DaprBlockingClient, inp::dapr.proto.runtime.v1.SetMetadataRequest) = SetMetadata(client.stub, client.controller, inp)
SetMetadata(client::DaprClient, inp::dapr.proto.runtime.v1.SetMetadataRequest, done::Function) = SetMetadata(client.stub, client.controller, inp, done)

import .dapr.proto.runtime.v1: Shutdown
"""
    Shutdown

- input: google.protobuf.Empty
- output: google.protobuf.Empty
"""
Shutdown(client::DaprBlockingClient, inp::google.protobuf.Empty) = Shutdown(client.stub, client.controller, inp)
Shutdown(client::DaprClient, inp::google.protobuf.Empty, done::Function) = Shutdown(client.stub, client.controller, inp, done)

# end service: dapr.proto.runtime.v1.Dapr

end # module DaprClients


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.SimulationEngine(ModiaMath/SimulationEngine/_module.jl)
#

# Methods that involve IDA UserData

#=
"KLU sparse matrix type"
type SlsMat
   M::Cint                        # number of rows
   N::Cint                        # number of columns
   NNZ::Cint                      # maximum number of nonzero entries in the matrix
   data::Ptr{Sundials.realtype}   # data[NNZ]    - value of nonzero entries
   rowvals::Ptr{Cint}             # rowvals[NNZ] - row indices of data (index starts at zero)
   colptrs::Ptr{Cint}             # colptrs[N+1] - index of the first column entry into the
                                  #                data and rowvals arrays (index starts at zero)
                                  #                The last entry contains the total number of nonzero values
                                  #                in the matrix and hence points one past the end
                                  #                of the active data in the data and rowvals arrays.
end

const SlsMat_Ptr = Ptr{SlsMat}

"""
    CSCtoSlsMat!(mat, matKLU)

Copy SparseMatrixCSC mat to SlsMat data structure matKLU
(storage for both objects are provided in the calling program)
"""
function CSCtoSlsMat!(mat::SparseMatrixCSC{Float64,Cint}, matKLU::SlsMat)
   for i = 1:length(mat.nzval)
      unsafe_store!(matKLU.data   , mat.nzval[i], i)
      unsafe_store!(matKLU.rowvals, mat.rowval[i]-1, i)
   end

   for i = 1:length(mat.colptr)
      unsafe_store!(matKLU.colptrs, mat.colptr[i]-1, i)
   end
end
=#

##################################################################
#
# Wrappers to IDA code
#
##################################################################

function idasol_f(t::Sundials.realtype, _y::Sundials.N_Vector, _yp::Sundials.N_Vector, _r::Sundials.N_Vector, simModel::IntegratorData)
    Sundials.IDAGetCurrentStep(simModel.ida_mem, simModel.hcur)
    Sundials.IDAGetCurrentOrder(simModel.ida_mem, simModel.order)
    stat           = simModel.statistics
    stat.hMin      = min(stat.hMin, simModel.hcur[1])
    stat.hMax      = max(stat.hMax, simModel.hcur[1])
    stat.orderMax  = max(stat.orderMax, simModel.order[1])
    stat.nResidues += 1

    sim                  = simModel.simulationState
    sim.time             = t
    simModel.last_t      = t
    simModel.last_norm_y = norm(simModel.y, Inf)

    # Check that simModel.y is still identical to _y
    if pointer(simModel.y) === Sundials.__N_VGetArrayPointer_Serial(_y)
       # If this assumption is true, no unnecessary memory is allocated via Sundials.asarray(..)
        y  = simModel.y
    else
        y  = Sundials.asarray(_y)
        println("\n!!! Info message from ModiaMath.simulate:\n",
                 "    idasol_f assumption SimModel.y === _y is not valid; using asarray(_y).")
    end

    if pointer(simModel.yp) === Sundials.__N_VGetArrayPointer_Serial(_yp)
       # If this assumption is true, no unnecessary memory is allocated via Sundials.asarray(..)
        yp = simModel.yp
    else
        yp = Sundials.asarray(_yp)
        println("\n!!! Info message from ModiaMath.simulate:\n",
                 "    idasol_f assumption SimModel.yp === _yp is not valid; using asarray(_yp).")
    end

    ModiaMath.DAE.getResidues!(simModel.model, sim, t, y, yp, simModel.r, simModel.hcur[1])

    # Copy simModel.r to _r
    @static if VERSION >= v"0.7.0-DEV.2005"
        unsafe_copyto!(Sundials.__N_VGetArrayPointer_Serial(_r), pointer(simModel.r), simModel.simulationState.nx)
    else
        unsafe_copy!(Sundials.__N_VGetArrayPointer_Serial(_r), pointer(simModel.r), simModel.simulationState.nx)
    end
    return Cint(0)   # indicates normal return
end

@static if VERSION >= v"0.7.0-DEV.2005"
    @noinline function old_cfunction(f, r, a)
        ccall(:jl_function_ptr, Ptr{Cvoid}, (Any, Any, Any), f, r, a)
    end
else
    const idasol_fc = cfunction(idasol_f, Cint, (Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Sundials.N_Vector, Ref{IntegratorData}))
end

#------- root finding
function idasol_g(t::Sundials.realtype, y::Sundials.N_Vector, yp::Sundials.N_Vector, gout::Ptr{Sundials.realtype}, simModel::IntegratorData)
    simModel.statistics.nZeroCrossings += 1
    sim      = simModel.simulationState
    sim.time = t
    ModiaMath.DAE.getEventIndicators!(simModel.model, sim, t, Sundials.asarray(y), Sundials.asarray(yp), simModel.z)

    for i = 1:simModel.simulationState.nz
        unsafe_store!(gout, simModel.z[i], i)
    end
    return Cint(0)   # indicates normal return
end

@static if VERSION < v"0.7.0-DEV.2005"
    const idasol_gc = cfunction(idasol_g, Cint, (Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Ptr{Sundials.realtype}, Ref{IntegratorData}))
end

#=
#-------- Jacobian
function idasol_sjac(t::Sundials.realtype, c_j::Sundials.realtype, y::Sundials.N_Vector, yp::Sundials.N_Vector, r::Sundials.N_Vector,
                     jacKLU::SlsMat, simModel::IntegratorData,
                     tmp1::Sundials.N_Vector, tmp2::Sundials.N_Vector, tmp3::Sundials.N_Vector)
    # Compute non-zero values of jac
    sim      = simModel.simulationState
    sim.time = t
    h        = IDAGetCurrentStep(simModel.ida_mem)
    IDAGetErrWeights(simModel.ida_mem, simModel.eweight)

    ModiaMath.SparseJacobian.computeJacobian!(ModiaMath.DAE.getResidues!, simModel.model, simModel.simulationState,
                                              t, Sundials.asarray(y), Sundials.asarray(yp), Sundials.asarray(r),
                                              sim.cg, c_j, h, Sundials.asarray(simModel.eweight), Sundials.asarray(tmp1), sim.jac)
    simModel.statistics.nResidues += sim.cg.ngroups

    # Copy userData.jac to jacKLU
    #   function call gives a Julia crash: CSCtoSlsMat!(simModel.cm.jac, jacKLU)
    #   Manually inlining the function works
    jac = sim.jac
    for i = 1:length(jac.nzval)
      unsafe_store!(jacKLU.data   , jac.nzval[i], i)
      unsafe_store!(jacKLU.rowvals, jac.rowval[i]-1, i)
    end

    for i = 1:length(jac.colptr)
      unsafe_store!(jacKLU.colptrs, jac.colptr[i]-1, i)
    end

    return Cint(0)   # indicates normal return
end

@static if VERSION < v"0.7.0-DEV.2005"
    const idasol_sjacc = @cfunction(idasol_sjac, Int32, (Sundials.realtype, Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Sundials.N_Vector, Ref{SlsMat},
                                                         Ref{IntegratorData}, Sundials.N_Vector, Sundials.N_Vector, Sundials.N_Vector))
else
    const idasol_sjacc = cfunction(idasol_sjac, Int32, (Sundials.realtype, Sundials.realtype, Sundials.N_Vector, Sundials.N_Vector, Sundials.N_Vector, Ref{SlsMat},
                                                        Ref{IntegratorData}, Sundials.N_Vector, Sundials.N_Vector, Sundials.N_Vector))
end
=#

function idasol_ErrHandlerFn(error_code::Cint, IDAmodule::Cstring, IDAfunction::Cstring,
                             message::Cstring, simModel::IntegratorData)
    modelName = simModel.model.modelName
    if error_code > 0
        # Print information text
        println("\n\n!!! Warning from ModiaMath.simulate(", modelName, ", ...): ", unsafe_string(IDAfunction),
               "(...) returned with [", unsafe_string(IDAmodule), "] error_code  = ",
               error_code, ":\n", unsafe_string(message), "\n")
    elseif error_code < 0
        # Raise error
        time   = simModel.last_t
        norm_y = simModel.last_norm_y

        if time > typemin(Float64)
            error("\n\n!!! Error from ModiaMath.simulate(", modelName, ", ...): ", unsafe_string(IDAfunction),
                "(...) returned with an [", unsafe_string(IDAmodule), "] error:\n",
                unsafe_string(message), "\n( norm( x(t=", @sprintf("%0.3g",time), " s) ) = ", @sprintf("%0.3g",norm_y),
                "; if this value is large, the model is unstable )\n")
        else
            error("\n\n!!! Error from ModiaMath.simulate(", modelName, ", ...): ", unsafe_string(IDAfunction),
                "(...) returned with an [", unsafe_string(IDAmodule), "] error:\n",
                unsafe_string(message), "\n")
        end
    end
    nothing
end

@static if VERSION < v"0.7.0-DEV.2005"
    const idasol_ErrHandlerFnc = cfunction(idasol_ErrHandlerFn, Void, (Cint, Cstring, Cstring, Cstring, Ref{IntegratorData}))
end


text:
cluster = [1 1 1 1; 2 2 3 3; 4 4 5 5]
X = Array(reshape(1:24, 2, 3, 4))

@testset "pool" begin
    @testset "GlobalPool" begin
        glb_cltr = [1 1 1 1; 1 1 1 1; 1 1 1 1]
        p = GlobalPool(:add, 3, 4)
        @test p(X) == sumpool(glb_cltr, X)
    end

    @testset "LocalPool" begin
        p = LocalPool(:add, cluster)
        @test p(X) == sumpool(cluster, X)
    end

    @testset "TopKPool" begin
        N = 10
        k, in_channel = 4, 7
        X = rand(in_channel, N)
        for T = [Bool, Float64]
            adj = rand(T, N, N)
            p = TopKPool(adj, k, in_channel)
            @test eltype(p.p) === Float32
            @test size(p.p) == (in_channel,)
            @test eltype(p.Ã) === T
            @test size(p.Ã) == (k, k)
            y = p(X)
            @test size(y) == (in_channel, k)
        end
    end

    for T = [UInt32, UInt64]
        @testset "$(T)" begin
            @testset "sumpool" begin
                y = [40, 44, 12, 14, 36, 38, 16, 18, 40, 42]
                @test sumpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:add, cluster, X) == reshape(y, 2, 5)
            end

            @testset "maxpool" begin
                y = [19, 20, 9, 10, 21, 22, 11, 12, 23, 24]
                @test maxpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:max, cluster, X) == reshape(y, 2, 5)
            end

            @testset "minpool" begin
                y = [1, 2, 3, 4, 15, 16, 5, 6, 17, 18]
                @test minpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:min, cluster, X) == reshape(y, 2, 5)
            end
        end
    end

    for T = [Int32, Int64]
        @testset "$(T)" begin
            @testset "sumpool" begin
                y = [40, 44, 12, 14, 36, 38, 16, 18, 40, 42]
                @test sumpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:add, cluster, X) == reshape(y, 2, 5)
            end

            @testset "subpool" begin
                y = [-40, -44, -12, -14, -36, -38, -16, -18, -40, -42]
                @test subpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:sub, cluster, X) == reshape(y, 2, 5)
            end

            @testset "maxpool" begin
                y = [19, 20, 9, 10, 21, 22, 11, 12, 23, 24]
                @test maxpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:max, cluster, X) == reshape(y, 2, 5)
            end

            @testset "minpool" begin
                y = [1, 2, 3, 4, 15, 16, 5, 6, 17, 18]
                @test minpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:min, cluster, X) == reshape(y, 2, 5)
            end
        end
    end

    for T = [Float16, Float32, Float64]
        @testset "$(T)" begin
            @testset "sumpool" begin
                y = [40, 44, 12, 14, 36, 38, 16, 18, 40, 42]
                @test sumpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:add, cluster, X) == reshape(y, 2, 5)
            end

            @testset "subpool" begin
                y = [-40, -44, -12, -14, -36, -38, -16, -18, -40, -42]
                @test subpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:sub, cluster, X) == reshape(y, 2, 5)
            end

            @testset "maxpool" begin
                y = [19, 20, 9, 10, 21, 22, 11, 12, 23, 24]
                @test maxpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:max, cluster, X) == reshape(y, 2, 5)
            end

            @testset "minpool" begin
                y = [1, 2, 3, 4, 15, 16, 5, 6, 17, 18]
                @test minpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:min, cluster, X) == reshape(y, 2, 5)
            end

            @testset "prodpool" begin
                y = [1729, 4480, 27, 40, 315, 352, 55, 72, 391, 432]
                @test prodpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:mul, cluster, X) == reshape(y, 2, 5)
            end

            @testset "divpool" begin
                y = 1 ./ [1729, 4480, 27, 40, 315, 352, 55, 72, 391, 432]
                @test divpool(cluster, X) ≈ reshape(y, 2, 5)
                @test pool(:div, cluster, X) ≈ reshape(y, 2, 5)
            end

            @testset "meanpool" begin
                y = [10., 11., 6., 7., 18., 19., 8., 9., 20., 21.]
                @test meanpool(cluster, X) == reshape(y, 2, 5)
                @test pool(:mean, cluster, X) == reshape(y, 2, 5)
            end
        end
    end
end


text:
# cache strategy -- the interface
struct LocalWindowCache{N} <: AbstractCacheStrategy
    window_size::NTuple{N, Int}
end
LocalWindowCache(window_size::Tuple) = LocalWindowCache{length(window_size)}(map(i->convert(Int, i), window_size))
make_cache(T, S::LocalWindowCache{N}, ax::Tuple) where N = LocalWindowCacheArray{T}(ax..., S.window_size)

# cache implementation
struct LocalWindowCacheArray{T, N, NA, AT<:AbstractArray{T, N}, IA, IB} <: CachedView{T, N}
    cache::AT
    axesA::IA
    axesB::IB
end

function LocalWindowCacheArray{T}(axesA::NTuple{NA}, axesB::NTuple{NA}, window_size::NTuple{NA, Integer}) where {T, NA}
    window_ax = map(window_size) do w
        r = w ÷ 2
        -r:r
    end
    cache_ax = (axesA..., window_ax...)
    N = length(axesA)+length(axesB)

    FT = Union{Missing, T}
    cache = OffsetArray{FT, N}(missing, cache_ax)
    LocalWindowCacheArray{FT, N, length(axesA), typeof(cache), typeof(axesA), typeof(axesB)}(cache, axesA, axesB)
end

Base.axes(A::LocalWindowCacheArray) = (A.axesA..., A.axesB...)
Base.size(A::LocalWindowCacheArray) = map(length, axes(A))


text:
@info "Installing and precompiling packages"

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
Pkg.precompile()

Pkg.add(url="https://github.com/JuliaTeachingCTU/ImageInspector.jl", rev="master")

using BSON
using BenchmarkTools
using CSV
using DataFrames
using DifferentialEquations
using Distributions
using Flux
using GLM
using GLPK
using HypothesisTests
using Ipopt
using JuMP
using LinearAlgebra
using MLDatasets
using PkgTemplates
using Plots
using ProgressMeter
using Query
using RDatasets
using Random
using SpecialFunctions
using Statistics
using StatsPlots

@info "Downloading datasets"

ENV["DATADEPS_ALWAYS_ACCEPT"] = true

dataset("datasets", "iris"); # download iris dataset
dataset("plm", "Snmesp"); # download Snmesp dataset

MNIST.traindata();  # download MNIST dataset
FashionMNIST.traindata();  # download FashionMNIST dataset
CIFAR10.traindata();  # download CIFAR10 dataset

@info "Finished"


text:
export LSR1Operator, diag  #, InverseLSR1Operator

"A data type to hold information relative to LSR1 operators."
mutable struct LSR1Data{T,I<:Integer}
  mem::I
  scaling::Bool
  scaling_factor::T
  s::Vector{Vector{T}}
  y::Vector{Vector{T}}
  ys::Vector{T}
  a::Vector{Vector{T}}
  as::Vector{T}
  insert::I
  Ax::Vector{T}
end

function LSR1Data(T::DataType, n::I; mem::I = 5, scaling::Bool = true, inverse::Bool = false) where {I<:Integer}
  inverse && @warn "inverse LSR1 operator not yet implemented"
  LSR1Data{T,I}(
    max(mem, 1),
    scaling,
    convert(T, 1),
    [zeros(T, n) for _ = 1:mem],
    [zeros(T, n) for _ = 1:mem],
    zeros(T, mem),
    [zeros(T, n) for _ = 1:mem],
    zeros(T, mem),
    1,
    Vector{T}(undef, n),
  )
end

LSR1Data(n::I; kwargs...) where {I<:Integer} = LSR1Data(Float64, n; kwargs...)

"A type for limited-memory SR1 approximations."
mutable struct LSR1Operator{T,S,I<:Integer,F,Ft,Fct} <: AbstractLinearOperator{T}
  nrow::I
  ncol::I
  symmetric::Bool
  hermitian::Bool
  prod!::F     # apply the operator to a vector
  tprod!::Ft    # apply the transpose operator to a vector
  ctprod!::Fct   # apply the transpose conjugate operator to a vector
  Mv::S # storage vector for prod!
  Mtu::S # storage vector for tprod!
  Maw::S # storage vector for ctprod!
  inverse::Bool
  data::LSR1Data{T,I}
  nprod::I
  ntprod::I
  nctprod::I
end

LSR1Operator{T}(
  nrow::I,
  ncol::I,
  symmetric::Bool,
  hermitian::Bool,
  prod!::F,
  tprod!::Ft,
  ctprod!::Fct,
  Mv::S,
  Mtu::S,
  Maw::S,
  inverse::Bool,
  data::LSR1Data{T,I},
) where {T,S,I<:Integer,F,Ft,Fct} =
  LSR1Operator{T,S,I,F,Ft,Fct}(nrow, ncol, symmetric, hermitian, prod!, tprod!, ctprod!, Mv, Mtu, Maw, inverse, data, 0, 0, 0)

"""
    LSR1Operator(T, n; [mem=5, scaling=false)
    LSR1Operator(n; [mem=5, scaling=false)
Construct a limited-memory SR1 approximation in forward form. If the type `T` is
omitted, then `Float64` is used.
"""
function LSR1Operator(T::DataType, n::I; kwargs...) where {I<:Integer}
  lsr1_data = LSR1Data(T, n; kwargs...)

  function lsr1_multiply(q::AbstractVector, data::LSR1Data, x::AbstractArray, α, β)
    # Multiply operator with a vector.

    # q = data.Ax
    q .= α .* x ./ data.scaling_factor .+ β .* q

    # data.scaling && (q ./= data.scaling_factor)  # q = B₀ * x

    for i = 1:(data.mem)
      k = mod(data.insert + i - 2, data.mem) + 1
      if data.ys[k] != 0
        ax = dot(data.a[k], x) / data.as[k]
        for j ∈ eachindex(q)
          q[j] += α * ax * data.a[k][j]
        end
      end
    end
  end

  prod! = @closure (res, x, α, β) -> lsr1_multiply(res, lsr1_data, x, α, β)
  Mv = lsr1_data.Ax
  return LSR1Operator{T}(n, n, true, true, prod!, nothing, nothing, Mv, Mv, Mv, false, lsr1_data)
end

LSR1Operator(n::I; kwargs...) where {I<:Integer} = LSR1Operator(Float64, n; kwargs...)

"""
    push!(op, s, y)
Push a new {s,y} pair into a L-SR1 operator.
"""
function push!(op::LSR1Operator, s::AbstractVector, y::AbstractVector)

  # op.counters.updates += 1
  data = op.data
  Bs = op * s
  ymBs = y - Bs
  ys = dot(y, s)
  sNorm = norm(s)
  yy = dot(y, y)

  ϵ = eps(eltype(op))
  well_defined = abs(dot(ymBs, s)) ≥ ϵ + ϵ * norm(ymBs) * sNorm

  sufficient_curvature = true
  scaling_condition = true
  if data.scaling
    yNorm = √yy
    sufficient_curvature = abs(ys) ≥ ϵ * yNorm * sNorm
    if sufficient_curvature
      scaling_factor = ys / yy
      scaling_condition = norm(y - s / scaling_factor) >= ϵ * yNorm * sNorm
    end
  end

  if !(well_defined && sufficient_curvature && scaling_condition)
    @debug "LSR1 update rejected" well_defined sufficient_curvature scaling_condition
    # op.counters.rejects += 1
    return op
  end

  data.s[data.insert] .= s
  data.y[data.insert] .= y
  data.ys[data.insert] = ys

  # update scaling factor
  data.scaling && (data.scaling_factor = ys / yy)

  # update next insertion position
  data.insert = mod(data.insert, data.mem) + 1

  # update rank-1 terms
  for i = 1:(data.mem)
    k = mod(data.insert + i - 2, data.mem) + 1
    if data.ys[k] != 0
      data.a[k] .= data.y[k] - data.s[k] / data.scaling_factor  # = y - B₀ * s
      for j = 1:(i - 1)
        l = mod(data.insert + j - 2, data.mem) + 1
        if data.ys[l] != 0
          as = dot(data.a[l], data.s[k]) / data.as[l]
          data.a[k] .-= as * data.a[l]
        end
      end
      data.as[k] = dot(data.a[k], data.s[k])
    end
  end

  return op
end

"""
    diag(op)
    diag!(op, d)
Extract the diagonal of a L-SR1 operator in forward mode.
"""
function diag(op::LSR1Operator{T}) where {T}
  d = Vector{T}(undef, op.nrow)
  diag!(op, d)
end

function diag!(op::LSR1Operator{T}, d) where {T}
  op.inverse && throw("only the diagonal of a forward L-SR1 approximation is available")
  data = op.data

  fill!(d, 1)
  data.scaling && (d ./= data.scaling_factor)

  for i = 1:(data.mem)
    k = mod(data.insert + i - 2, data.mem) + 1
    if data.ys[k] != 0.0
      for j = 1:(op.nrow)
        d[j] += data.a[k][j]^2 / data.as[k]
      end
    end
  end
  return d
end

"""
    reset!(data)
Reset the given LSR1 data.
"""
function reset!(data::LSR1Data{T,I}) where {T,I<:Integer}
  for i = 1:(data.mem)
    fill!(data.s[i], 0)
    fill!(data.y[i], 0)
    fill!(data.a[i], 0)
  end
  fill!(data.ys, 0)
  fill!(data.as, 0)
  data.scaling_factor = T(1)
  data.insert = 1
  return data
end

"""
    reset!(op)
Resets the LSR1 data of the given operator.
"""
function reset!(op::LSR1Operator)
  reset!(op.data)
  op.nprod = 0
  op.ntprod = 0
  op.nctprod = 0
  return op
end