text:
# Copyright 2021 Gravifer
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Test
using BenchmarkTools
using OffsetArrays

@testset "transfinite_interpolate" begin
    include(normpath(joinpath(@__DIR__, "../../src/misc-util.jl")))
    using .__CFD2021__misc_util__: tuplejoin

    include("../src/transfinite-interpolate.jl")
    import .transfinite_interpolate: transfinite_interpolate_2d!, transfinite_interpolate_2d
    @testset "transfinite_interpolate_2d" begin
        Mx, My = 3, 3

        @testset "on functions" begin
            interp_func = transfinite_interpolate_2d( (x->2x, y->3y), (x->4x+3, y->2+5y) )
            @test interp_func((0.5,0.5)) == 3.0
        end

        @testset "on arrays" begin
            x_lo = ([LinRange(0, 1, Mx)...], [LinRange(0, 0, My)...]); x_hi = ([LinRange(0, 1, Mx)...], [LinRange(1, 1, My)...]);
            y_lo = ([LinRange(0, 0, Mx)...], [LinRange(0, 1, My)...]); y_hi = ([LinRange(1, 1, Mx)...], [LinRange(0, 1, My)...]);
            in_situ_interp_x = Array{Float64}(undef, (length(v) for v in x_lo)...);
            in_situ_interp_y = Array{Float64}(undef, (length(v) for v in y_lo)...);
            transfinite_interpolate_2d!(in_situ_interp_x, x_lo, x_hi)
            transfinite_interpolate_2d!(in_situ_interp_y, y_lo, y_hi)
            ab_situ_interp_x = transfinite_interpolate_2d(x_lo, x_hi)
            ab_situ_interp_y = transfinite_interpolate_2d(y_lo, y_hi)

            @testset "on ordinary arrays" begin
                @test in_situ_interp_x == ab_situ_interp_x == [0.0 0.0 0.0; 0.5 0.5 0.5; 1.0 1.0 1.0]
                @test in_situ_interp_y == ab_situ_interp_y == [0.0 0.5 1.0; 0.0 0.5 1.0; 0.0 0.5 1.0]
            end

            x_lo = map(x->OffsetArray(x, OffsetArrays.Origin(0)), x_lo)
            y_lo = map(x->OffsetArray(x, OffsetArrays.Origin(0)), y_lo)
            x_hi = map(x->OffsetArray(x, OffsetArrays.Origin(0)), x_hi)
            y_hi = map(x->OffsetArray(x, OffsetArrays.Origin(0)), y_hi)
            in_situ_interp_x = OffsetArray(in_situ_interp_x, tuplejoin((axes(v) for v in x_lo)...))
            in_situ_interp_y = OffsetArray(in_situ_interp_y, tuplejoin((axes(v) for v in y_lo)...))
            transfinite_interpolate_2d!(in_situ_interp_x, x_lo, x_hi)
            transfinite_interpolate_2d!(in_situ_interp_y, y_lo, y_hi)
            ab_situ_interp_x = transfinite_interpolate_2d(x_lo, x_hi)
            ab_situ_interp_y = transfinite_interpolate_2d(y_lo, y_hi)

            @testset "on OffsetArrays" begin
                @test in_situ_interp_x == ab_situ_interp_x != [0.0 0.0 0.0; 0.5 0.5 0.5; 1.0 1.0 1.0]
                @test in_situ_interp_y == ab_situ_interp_y != [0.0 0.5 1.0; 0.0 0.5 1.0; 0.0 0.5 1.0]
                @test in_situ_interp_x == ab_situ_interp_x == OffsetArray([0.0 0.0 0.0; 0.5 0.5 0.5; 1.0 1.0 1.0], OffsetArrays.Origin(0))
                @test in_situ_interp_y == ab_situ_interp_y == OffsetArray([0.0 0.5 1.0; 0.0 0.5 1.0; 0.0 0.5 1.0], OffsetArrays.Origin(0))
            end

            @testset "exception handling" begin
                x_lo = ([LinRange(0, 1, Mx)...], [LinRange(0, 0, My)...]); x_hi = ([LinRange(0, 1, 2Mx+1)...], [LinRange(1, 1, 2My+1)...]);
                @test_throws DimensionMismatch transfinite_interpolate_2d(x_lo, x_hi)
            end

        end # @testset "on arrays"
    end # @testset "transfinite_interpolate_2d"
end # @testset "transfinite_interpolate"

@testset "elliptic_grid_gen" begin
    include("../src/elliptic-grid-gen.jl")
    import .elliptic_grid_gen

    @testset "inverted_poisson_jacobi" begin
        @test_skip elliptic_grid_gen.inverted_poisson_jacobi_step
        @test_skip elliptic_grid_gen.inverted_poisson_jacobi_step!
        @test_skip elliptic_grid_gen.inverted_poisson_jacobi_iterate
        @test_skip elliptic_grid_gen.inverted_poisson_jacobi_iterate!
    end # @testset "inverted_poisson_jacobi"
end # @testset "elliptic_grid_gen"


text:
using GalerkinSparseGrids
using Base.Test
using Cubature
using ODE

#--------------------------------------
# Testing PDE solvers
#--------------------------------------

# The wave equation, testing conservation of energy

# In the position basis:

import GalerkinSparseGrids.energy_func_1D
import GalerkinSparseGrids.energy_func

print("Testing wave equation solver 1-D position DG basis... ")

k = 4
level = 4
f0 = x->sin(2*pi*x[1])
v0 = x->2*pi*cos(2*pi*x[1])
pos_soln = wave_evolve_1D(k, level, f0, v0, 0, 1; basis="pos", order="45")
energy_soln = energy_func_1D(k, level, pos_soln; basis="pos")
for i in 1:length(energy_soln[2])
    @test sqrt(energy_soln[2][i]) ≈ 2*pi atol=1.0e-7
end

println("Test Passed.")

# In the hierarchical basis:

print("Testing wave equation solver 1-D hierarchical DG basis... ")

hier_soln  = wave_evolve_1D(k, level, f0, v0, 0, 1; basis="hier", order="45")
energy_soln = energy_func_1D(k, level, hier_soln; basis="hier")
for i in 1:length(energy_soln[2])
    @test sqrt(energy_soln[2][i]) ≈ 2*pi atol=1.0e-7
end

hier_soln  = wave_evolve(1, k, level, f0, v0, 0, 1; order="78")
energy_soln = energy_func(1, k, level, hier_soln)
for i in 1:length(energy_soln[2])
    @test sqrt(energy_soln[2][i]) ≈ 2*pi atol=1.0e-7
end

println("Test Passed.")

# 2D Sparse case:

print("Testing wave equation solver on 2-D sparse DG basis... ")

D = 2
k_used = 3
n_used = 5
f0 = x->sin(2*pi*x[1])*sin(2*pi*x[2])
v0 = x->0
sparse_soln = wave_evolve(D, k_used, n_used, f0, v0, 0, 1; order="78", scheme="sparse");
senergy = energy_func(D, k_used, n_used, sparse_soln)
@test senergy[2][1]-senergy[2][end] > 0
@test abs(senergy[2][1]-senergy[2][end]) < 1.0e-8
for i in 1:length(senergy[1])
    @test sqrt(senergy[2][i]) ≈ sqrt(2)*pi atol=1.0e-4
end

sparse_soln = wave_evolve(D, k_used, n_used, f0, v0, 0,1; order="45", scheme="sparse");
senergy = energy_func(D, k_used, n_used, sparse_soln)
@test senergy[2][1] > senergy[2][end]
@test senergy[2][1] ≈ senergy[2][end] atol=1.0e-8
for i in 1:length(senergy[1])
    @test sqrt(senergy[2][i]) ≈ sqrt(2)*pi atol=1.0e-4
end

println("Test Passed.")

print("Testing traveling wave example in 2-D... ")

m = [1, 2]
truesoln = x -> cos(2*pi*(vecdot(m,x) - sqrt(vecdot(m,m))*0.54))

D = length(m)
k_used = 3
n_used = 6

soln = traveling_wave_solver(k_used, n_used, m, 0, 0.54; order="45")
dict = V2D(D, k_used, n_used, soln[2][end]; scheme="sparse")
@test mcerr(x->reconstruct_DG(dict, [x...]), truesoln, D) < 0.05

soln = traveling_wave_solver(k_used, n_used, m, 0, 0.54; order="78")
dict = V2D(D, k_used, n_used, soln[2][end]; scheme="sparse")
@test mcerr(x->reconstruct_DG(dict, [x...]), truesoln, D) < 0.05

println("Test passed.")

ext: jl


text:
using Test, DataFrames, aoc2019

@testset "base" begin
    @test length(aoc2019.answerfuncs) > 0
    @test length(aoc2019.answerfuncs) == length(aoc2019.inputreads)
    @test length(aoc2019.answerfuncs) == length(aoc2019.inputpaths)
    @test aoc2019.get_answers() isa DataFrames.DataFrame
end

include("d01.jl")
include("d02.jl")
include("d03.jl")


text:
@pp function model()
    local g = mem(_->flip(0.7))
    flip(0.6)
    g(0)
end

for algo in [enum, enumdepthfirst, enumlikelyfirst]
    dist = @pp algo(model)
    @test_approx_eq(exp(score(dist, true)), 0.7)
    @test_approx_eq(exp(score(dist, false)), 0.3)
end

comps = [
() -> begin
    @pp pmcmc(1,100) do
        local f = mem(_->flip())
        f(0)
    end
end,
() -> begin
    @pp begin
        local f = mem(_->flip())
        pmcmc(1,100) do
            f(0)
        end
    end
end,
() -> begin
    f = @pp mem(_->flip())
    @pp pmcmc(1,100) do
        f(0)
    end
end]

for comp in comps
    supp = support(comp())
    @test length(supp) == 2
    @test true in supp
    @test false in supp
end

# This is a hack used to ensure that all the first particle has all of
# the probability mass at the first factor statement. Therefore, after
# resampling, all particles are descendants of the same particle. When
# the store isn't correctly handled the memo cache ends up been shared
# resulting in incorrect results.

# I think this made more sense when continuations were copied during
# resampling. I'm not sure how this could fail now, but it's worth
# keeping around.

c = 0
function count(store,k,address)
    global c
    c += 1
    k(store,c)
end

dist = @pp pmcmc(1,100) do
    local f = mem(_->flip())
    factor(count() == 1 ? 0 : -Inf)
    f(0)
end
supp = support(dist)

@test length(supp) == 2
@test true in supp
@test false in supp


text:
module HigherOrderKernels

export GaussianKernel,
       PolynomialKernel,
       UniformKernel,
       EpanechnikovKernel,
       BiweightKernel,
       TriweightKernel,
       density_bandwidth,
       distribution_bandwidth,
       kpdf,
       kcdf,
       order,
       poly_order

using GaussQuadrature,
      Polynomials,
      Statistics,
      SpecialFunctions

abstract type AbstractKernel{ν} end
struct GaussianKernel{ν} <: AbstractKernel{ν} end
struct PolynomialKernel{s,ν} <: AbstractKernel{ν} end
const UniformKernel = PolynomialKernel{0}
const EpanechnikovKernel = PolynomialKernel{1}
const BiweightKernel = PolynomialKernel{2}
const TriweightKernel = PolynomialKernel{3}

poly_order(k::Type{T}) where {s,ν,T<:PolynomialKernel{s,ν}} =
  degree(_u_poly(k))

order(::Type{T}) where {ν,T<:AbstractKernel{ν}} = ν

function roughness(k::Type{T}) where {s,ν,T<:PolynomialKernel{s,ν}}
  # N = 2s+2(div(ν,2)-1)+1
  # x, w = legendre(N)
  # R = sum(w[i] * density_kernel(k, x[i])^2 for i=1:N)
  Float64(polyint(_u_poly(k)^2, -1, 1))
end

function distribution_roughness(k::Type{T}) where {s,ν,T<:PolynomialKernel{s,ν}}
  I = polyint(_u_poly(k))
  K = I - I(-1)
  Float64(1 - polyint(K^2, -1, 1))
end

function firstnzmoment(k::Type{T}) where {s,ν,T<:PolynomialKernel{s,ν}}
  # N = s+2*div(ν,2)
  # x, w = legendre(N)
  # κ = sum(w[i] * x[i]^ν * density_kernel(k, x[i]) for i=1:N)
  polyint(Poly([zeros(ν); 1]) * _u_poly(k), -1, 1)
end

function roughness(k::Type{T}) where {ν,T<:GaussianKernel{ν}}
  N = 2(div(ν,2)-1)+1
  x, w = hermite(N)
  R = sum(w[i] * exp(x[i]^2) * density_kernel(k, x[i])^2 for i=1:N)
end

function firstnzmoment(k::Type{T}) where {ν,T<:GaussianKernel{ν}}
  N = 2*div(ν,2)
  x, w = hermite(N)
  κ = sum(w[i] * (sqrt(2) * x[i])^ν * exp(x[i]^2) * density_kernel(k, sqrt(2) * x[i]) * sqrt(2) for i=1:N)
end

function _std(x)
  iqr = diff(quantile(x, [0.25, 0.75]))[1] / 1.3489795003921636 # diff(quantile(Normal, [0.25, 0.75]))
  σ = std(x)
  min(iqr, σ)
end

const density_bandwidth_constant_lookup = Dict{Type{T} where T<:AbstractKernel,Float64}()

function density_bandwidth_constant(k::Type{T}) where {ν,T<:AbstractKernel{ν}}
  if haskey(density_bandwidth_constant_lookup, k)
    return density_bandwidth_constant_lookup[k]
  end
  R = roughness(k)
  κ = firstnzmoment(k)
  C = Float64(2*((sqrt(pi)*factorial(big(ν))^3*R)/(2*ν*factorial(big(2ν))*κ^2))^(1/(2ν+1)))
  density_bandwidth_constant_lookup[k] = C
end

density_bandwidth(k::Type{T}, data) where {ν,T<:AbstractKernel{ν}} =
  density_bandwidth_constant(k) * _std(data) * size(data, 1)^(-1/(2ν+1))

const distribution_bandwidth_constant_lookup = Dict{Type{T} where T<:AbstractKernel,Float64}()

function distribution_bandwidth_constant(k::Type{T}) where {ν,T<:AbstractKernel{ν}}
  if haskey(distribution_bandwidth_constant_lookup, k)
    return distribution_bandwidth_constant_lookup[k]
  end
  R = distribution_roughness(k)
  κ = firstnzmoment(k)
  C = Float64(2*((sqrt(pi)*factorial(big(ν))^3*R*(2ν-1))/(factorial(big(2ν))*κ^2*ν))^(1/(2ν-1)))
  distribution_bandwidth_constant_lookup[k] = C
end

distribution_bandwidth(k::Type{T}, data) where {ν,T<:AbstractKernel{ν}} =
  distribution_bandwidth_constant(k) * _std(data) * size(data, 1)^(-1/(2ν-1))

pochhammer(x, n) = gamma(x + n) / gamma(x)

# efficient evaluation of polynomial using Horner's rule
function polyval_expr(p::AbstractArray{T,1}, x::S) where {T,S}
  lenp = length(p)
  if lenp == 0
    return :($(zero(T)) * x)
  else
    y = convert(T, p[end])
    for i = lenp-1:-1:1
      y = :(muladd($x, $y, $(p[i])))
    end
    return y
  end
end

function _u2_poly(::Type{PolynomialKernel{s,ν}}) where {s,ν}
  r = div(ν, 2)
  c = pochhammer(big(1/2), s+1) * pochhammer(big(3/2+s), r-1) * pochhammer(big(3/2), r-1) / factorial(big(s+r-1))
  Poly([c * (-1)^k * pochhammer(big(1/2+s+r), k) / (factorial(big(k)) * factorial(big(r-1-k)) * pochhammer(big(3/2), k)) for k=0:r-1], :u2)
end

function _u_poly(k::Type{PolynomialKernel{s,ν}}) where {s,ν}
  p0 = Poly([1, -1], :u2)^s
  p1 = _u2_poly(k)
  p = p0 * p1
  u2 = Poly([0, 0, 1])
  p(u2)
end

@generated function density_kernel(k::Type{PolynomialKernel{s,ν}}, u) where {s,ν}
  e0 = (s == 0) ? :(abs(u) <= 1) :
       (s == 1) ? :((1-u2)*(abs(u) <= 1)) :
                  :((1-u2)^$s*(abs(u) <= 1))
  p1 = _u2_poly(PolynomialKernel{s,ν})
  e1 = polyval_expr(Float64.(coeffs(p1)), :u2)
  quote
    $(Expr(:meta, :inline))
    u2 = u^2
    $e1 * $e0
  end
end

@generated function distribution_kernel(k::Type{PolynomialKernel{s,ν}}, u) where {s,ν}
  p = _u_poly(PolynomialKernel{s,ν})
  P = polyint(p)
  P -= P(-1)
  e = polyval_expr(Float64.(coeffs(P)), :u)
  quote
    $(Expr(:meta, :inline))
    $e * (abs(u) <= 1) + (u > 1)
  end
end

@generated function density_kernel(::Type{GaussianKernel{ν}}, u) where {ν}
  r = div(ν, 2)
  ϕ = :(exp(-0.5u2))
  p = Poly(Float64.([(-1)^i / sqrt(2*pi) * factorial(big(2r))/(big(2)^(2r-i-1)*factorial(big(r))*factorial(big(2i+1))*factorial(big(r-i-1))) for i=0:r-1]))
  e = polyval_expr(coeffs(p), :u2)
  quote
    $(Expr(:meta, :inline))
    u2 = u^2
    $e * $ϕ
  end
end

kpdf(k, x, data) = kpdf(k, x, data, density_bandwidth(k, data))

function kpdf(::Type{T}, x, data, h) where {T<:AbstractKernel}
  sum = 0.0
  @fastmath @simd for y in data
    sum += density_kernel(T, (y - x) / h)
  end
  sum /= length(data) * h
  return sum
end

function kpdf(::Type{T}, x, data, h, lb, ub) where {T<:AbstractKernel}
  sum = 0.0
  if !(lb <= x <= ub)
    return sum
  end
  @simd for y in data
    sum += density_kernel(T, (y - x) / h)
  end
  if isfinite(lb)
    @simd for y in data
      sum += density_kernel(T, (y - (2lb - x)) / h)
    end
  end
  if isfinite(ub)
    @simd for y in data
      sum += density_kernel(T, (y - (2ub - x)) / h)
    end
  end
  sum /= length(data) * h
  return sum
end

kcdf(k, x, data) = kpdf(k, x, data, distribution_bandwidth(k, data))

function kcdf(::Type{T}, x, data, h) where {T<:AbstractKernel}
  sum = 0.0
  @fastmath @simd for y in data
    sum += distribution_kernel(T, (x - y) / h)
  end
  sum /= length(data)
  return sum
end

end # module


text:
module Simulate_PendulumWithFixedJoint

using  Modia3D
import Modia3D.ModiaMath

# Solids
groundMaterial = Modia3D.Material(color="DarkGreen", transparency=0.5)

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)
vmat2 = Modia3D.Material(color="Red")

@assembly Pendulum(;Lx = 1.0, Ly=0.2*Lx, Lz=0.2*Lx, m=1.0) begin
   world  = Modia3D.Object3D(Modia3D.CoordinateSystem(0.5*Lx))

   # Pendulum
   body1  = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBeam(Lx,Ly,Lz), Modia3D.MassProperties(m=m), vmat1))
   frame1 = Modia3D.Object3D(body1; r=[-Lx/2, 0.0, 0.0])
   cyl    = Modia3D.Object3D(frame1,Modia3D.Cylinder(Ly/2,1.2*Ly; material=vmat2))
   body2  = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), Modia3D.MassProperties(m=m), vmat1))
   frame2 = Modia3D.Object3D(body2; r=[-Lx/2, 0.0, 0.0])

   # Connect body1 to world and body2 to body1
   rev    = Modia3D.Revolute(world, frame1)
   fixed  = Modia3D.Fixed(body1, frame2, r=[Lx/2-Ly/2, Ly/2, 0.0], R=ModiaMath.rot3(pi/2))
end

const Lx = 1.6
const m  = 0.5
# Modia3D.visualizeAssembly!( Pendulum(Lx=Lx, m=m) )

pendulum = Pendulum(Lx=Lx, m=m, sceneOptions=Modia3D.SceneOptions(visualizeFrames=true, defaultFrameLength=0.3))
model    = Modia3D.SimulationModel( pendulum )
result   = ModiaMath.simulate!(model, stopTime=4.5, tolerance=1e-6,log=true)

ModiaMath.plot(result, ["rev.phi", "rev.w", "rev.a"])

println("... success of Simulate_PendulumWithFixedJoint.jl!")
end


text:
hdwf = HDWF(0)
szError = zeros(Cchar, 512)

@info "Automatically opening the first available device"
if FDwfDeviceOpen(-1, Ref(hdwf)) != 1
    FDwfGetLastErrorMsg(szError)
    error("Device open failed $szError")
end

@info "Generating waveform for 5 seconds..."
# enable first channel
FDwfAnalogOutNodeEnableSet(hdwf, 0, AnalogOutNodeCarrier, true)
# set sine function
FDwfAnalogOutNodeFunctionSet(hdwf, 0, AnalogOutNodeCarrier, funcSine)
# 10kHz
FDwfAnalogOutNodeFrequencySet(hdwf, 0, AnalogOutNodeCarrier, 10000.0)
# 1.41V amplitude (1Vrms), 2.82V pk2pk
FDwfAnalogOutNodeAmplitudeSet(hdwf, 0, AnalogOutNodeCarrier, 1.41)
# 1.41V offset
FDwfAnalogOutNodeOffsetSet(hdwf, 0, AnalogOutNodeCarrier, 1.41)
# start signal generation
FDwfAnalogOutConfigure(hdwf, 0, true)
# it will run until stopped, reset, parameter changed or device closed
sleep(5)
@info "done"

    # on close device is stopped and configuration lost
FDwfDeviceClose(hdwf)


text:
@testset "Events/Conditional Section" begin
    evaluate(s) = collect(E.events(L.State(IOBuffer(s))))

    E = Events
    L = Events.Lexical

    CSStart = E.ConditionalSectionStart
    CSEnd   = E.ConditionalSectionEnd

    DC = E.DataContent
    ME = E.MarkupError
    PE = E.EntityReferenceParameter

    @testset "Events/Conditional Section (Positive ... trivial content)" begin
        @test evaluate("<![INCLUDE[]]>") == [ CSStart("INCLUDE", L.Location("a buffer", -1)),
                                              CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![INCLUDE []]>") == [ CSStart("INCLUDE", L.Location("a buffer", -1)),
                                               CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![ INCLUDE[]]>") == [ CSStart("INCLUDE", L.Location("a buffer", -1)),
                                               CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![ INCLUDE []]>") == [ CSStart("INCLUDE", L.Location("a buffer", -1)),
                                                CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![\nINCLUDE\t[]]>") == [ CSStart("INCLUDE", L.Location("a buffer", -1)),
                                                  CSEnd(L.Location("a buffer", -1)) ]

        @test evaluate("<![IGNORE[]]>") == [ CSStart("IGNORE", L.Location("a buffer", -1)),
                                             CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![IGNORE []]>") == [ CSStart("IGNORE", L.Location("a buffer", -1)),
                                              CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![ IGNORE[]]>") == [ CSStart("IGNORE", L.Location("a buffer", -1)),
                                              CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![ IGNORE []]>") == [ CSStart("IGNORE", L.Location("a buffer", -1)),
                                               CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![\nIGNORE\t[]]>") == [ CSStart("IGNORE", L.Location("a buffer", -1)),
                                                 CSEnd(L.Location("a buffer", -1)) ]

        @test evaluate("<![%e;[]]>") == [ CSStart(PE("e", L.Location("a buffer", -1)), L.Location("a buffer", -1)),
                                          CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![%e; []]>") == [ CSStart(PE("e", L.Location("a buffer", -1)), L.Location("a buffer", -1)),
                                           CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![ %e;[]]>") == [ CSStart(PE("e", L.Location("a buffer", -1)), L.Location("a buffer", -1)),
                                           CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![ %e; []]>") == [ CSStart(PE("e", L.Location("a buffer", -1)), L.Location("a buffer", -1)),
                                            CSEnd(L.Location("a buffer", -1)) ]
        @test evaluate("<![\n%e;\t[]]>") == [ CSStart(PE("e", L.Location("a buffer", -1)), L.Location("a buffer", -1)),
                                              CSEnd(L.Location("a buffer", -1)) ]
    end

    @testset "Events/Conditional Section (Negative ... mixed-case keywords)" begin
        @test evaluate("<![include[") == [ ME("ERROR: The keyword 'include' must be uppercased.", L.Location("a buffer", -1)),
                                           CSStart(true, "include", L.Location("a buffer", -1)) ]
        @test evaluate("<![incLUDE[") == [ ME("ERROR: The keyword 'incLUDE' must be uppercased.", L.Location("a buffer", -1)),
                                           CSStart(true, "incLUDE", L.Location("a buffer", -1)) ]

        @test evaluate("<![ignore[") == [ ME("ERROR: The keyword 'ignore' must be uppercased.", L.Location("a buffer", -1)),
                                          CSStart(true, "ignore", L.Location("a buffer", -1)) ]
        @test evaluate("<![ignORE[") == [ ME("ERROR: The keyword 'ignORE' must be uppercased.", L.Location("a buffer", -1)),
                                          CSStart(true, "ignORE", L.Location("a buffer", -1)) ]
    end

    @testset "Events/Conditional Section (Negative ... no DSO)" begin
        @test evaluate("<![INCLUDE") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                             L.Location("a buffer", -1)),
                                          CSStart(true, "INCLUDE", L.Location("a buffer", -1)) ]
        @test evaluate("<![IGNORE") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                            L.Location("a buffer", -1)),
                                         CSStart(true, "IGNORE", L.Location("a buffer", -1)) ]
        @test evaluate("<![%e;") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                         L.Location("a buffer", -1)),
                                      CSStart(true, PE("e", L.Location("a buffer", -1)), L.Location("a buffer", -1)) ]

        # Make sure that leading and trailing white space doesn't throw off the error-reporting.
        #
        @test evaluate("<![INCLUDE ") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                              L.Location("a buffer", -1)),
                                           CSStart(true, "INCLUDE", L.Location("a buffer", -1)) ]
        @test evaluate("<![IGNORE ") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                             L.Location("a buffer", -1)),
                                          CSStart(true, "IGNORE", L.Location("a buffer", -1)) ]
        @test evaluate("<![%e; ") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                          L.Location("a buffer", -1)),
                                       CSStart(true, PE("e", L.Location("a buffer", -1)), L.Location("a buffer", -1)) ]

        @test evaluate("<![ INCLUDE") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                              L.Location("a buffer", -1)),
                                           CSStart(true, "INCLUDE", L.Location("a buffer", -1)) ]
        @test evaluate("<![ IGNORE") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                             L.Location("a buffer", -1)),
                                          CSStart(true, "IGNORE", L.Location("a buffer", -1)) ]
        @test evaluate("<![ %e;") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                          L.Location("a buffer", -1)),
                                       CSStart(true, PE("e", L.Location("a buffer", -1)), L.Location("a buffer", -1)) ]

        @test evaluate("<![ INCLUDE ") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                               L.Location("a buffer", -1)),
                                            CSStart(true, "INCLUDE", L.Location("a buffer", -1)) ]
        @test evaluate("<![ IGNORE ") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                              L.Location("a buffer", -1)),
                                           CSStart(true, "IGNORE", L.Location("a buffer", -1)) ]
        @test evaluate("<![ %e; ") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                           L.Location("a buffer", -1)),
                                        CSStart(true, PE("e", L.Location("a buffer", -1)), L.Location("a buffer", -1)) ]

        @test evaluate("<![\nINCLUDE\t") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                                 L.Location("a buffer", -1)),
                                              CSStart(true, "INCLUDE", L.Location("a buffer", -1)) ]
        @test evaluate("<![\nIGNORE\t") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                                L.Location("a buffer", -1)),
                                             CSStart(true, "IGNORE", L.Location("a buffer", -1)) ]
        @test evaluate("<![\n%e;\t") == [ ME("ERROR: Expecting '[' to open a conditional marked section.",
                                             L.Location("a buffer", -1)),
                                          CSStart(true, PE("e", L.Location("a buffer", -1)), L.Location("a buffer", -1)) ]
    end
end

max_stars_repo_name: jlegare/NativeXML

max_issues_repo_name: jlegare/NativeXML

max_forks_repo_name: jlegare/NativeXML


text:
rng = MemorizingRNG(MersenneTwister(12))

a = rand(rng)
rand(rng)
rand(rng)
rand(rng)

rand(rng, Float64)

rng.idx=0
@test rand(rng) == a

randn(rng)

rng.idx=0
b = randn(rng, 2)
randn(rng, 2)
randn(rng, 2)
randn(rng, 2)
rng.idx=0
@test b == randn(rng, 2)

src = MemorizingSource(1, 2, MersenneTwister(12), grow_reserve=false)

rng = ARDESPOT.get_rng(src, 1, 1)
r1 = rand(rng)
rng2 = ARDESPOT.get_rng(src, 1, 2)
r2 = rand(rng2)
rng = ARDESPOT.get_rng(src, 1, 1)
@test r1 == rand(rng)
r3 = rand(rng)
@test src.move_count == 1
rng2 = ARDESPOT.get_rng(src, 1, 2)
@test r2 == rand(rng2)
rng = ARDESPOT.get_rng(src, 1, 1)
@test r1 == rand(rng)
@test r3 == rand(rng)

@test src.move_count == 1

println("A warning should be thrown below.")
ARDESPOT.check_consistency(src)

# grow reserve
src = MemorizingSource(1, 4, MersenneTwister(12))

rng = ARDESPOT.get_rng(src, 1, 1)
r1 = rand(rng)
rng2 = ARDESPOT.get_rng(src, 1, 2)
r2 = rand(rng2)
rng = ARDESPOT.get_rng(src, 1, 1)
@test r1 == rand(rng)
r3 = rand(rng)
@test src.move_count == 1
@test src.min_reserve == 2
rng2 = ARDESPOT.get_rng(src, 1, 2)
@test r2 == rand(rng2)
rng = ARDESPOT.get_rng(src, 1, 1)
@test r1 == rand(rng)
@test r3 == rand(rng)

# create a new one, but don't get any numbers until creating the next - should have automatically reserved enough
rng3 = ARDESPOT.get_rng(src, 1, 3)
rng4 = ARDESPOT.get_rng(src, 1, 4)
r4 = rand(rng4)
rng3 = ARDESPOT.get_rng(src, 1, 3)
r5 = rand(rng3)
r6 = rand(rng3)

# check that this didn't trigger a move
@test src.move_count == 1


text:
for Arch in Archs
    @testset "Models [$Arch]" begin
        @info "Testing models [$Arch]..."

        @testset "Energy thermodynamic variable" begin
            @info "  Testing model construction with energy thermodynamic variable..."

            grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
            model = CompressibleModel(architecture=Arch(), grid = grid, gases = DryEarth(),
                                      thermodynamic_variable = Energy())

            @test model isa CompressibleModel
        end

        @testset "Entropy thermodynamic variable" begin
            @info "  Testing model construction with entropy thermodynamic variable..."

            grid = RegularCartesianGrid(size=(16, 16, 16), extent=(1, 1, 1))
            model = CompressibleModel(architecture=Arch(), grid = grid, gases = DryEarth(),
                                      thermodynamic_variable = Entropy())

            @test model isa CompressibleModel
        end
    end
end


text:
#
# This functions compute electrostatic interaction between two selections
#

#
# Computes the electrostatic interaction of a pair of atoms, given the distance
#

function qpair(d,q1,q2)
  Coulomb = 332.05382e0
  qpair = Coulomb*q1*q2 / d
  return qpair
end

#
# For a trajectory: x, y, z coordinates read from a DCD are given
#

# First, the general function with all parameters

function coulomb(atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int},
                 sides :: Vector{Float64},
                 x :: Vector{Float32}, y :: Vector{Float32}, z :: Vector{Float32},
                 usecutoff :: Bool, cutoff :: Float64, pbc :: Bool)
  n1 = length(sel1)
  n2 = length(sel2)
  coulomb = 0.
  for i in 1:n1
    center = [ x[sel1[i]], y[sel1[i]], z[sel1[i]] ]
    if pbc
      wrap!(sides,x,y,z,center=center,sel=sel2)
    end
    for j in 1:n2
      d = sqrt((x[sel1[i]] - x[sel2[j]])^2 + (y[sel1[i]] - y[sel2[j]])^2 + (z[sel1[i]] - z[sel2[j]])^2)
      coulomb = coulomb + qpair(d,atoms[sel1[i]].charge,atoms[sel2[j]].charge)
      if usecutoff
        if d < cutoff
          coulomb = coulomb - qpair(cutoff,atoms[sel1[i]].charge,atoms[sel2[j]].charge)
        end
      end
    end
  end
  return coulomb
end

# For a lazy user which does not want to setup the cutoff and the pbc, use defaults

coulomb(atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int},
        sides :: Vector{Float64},
        x :: Vector{Float32}, y :: Vector{Float32}, z :: Vector{Float32} ) =
  coulomb(atoms, sel1, sel2, sides, x, y, z, true, 12., true)

# If the user does not provide the sides, than we will not use PBCs

coulomb(atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int},
        x :: Vector{Float32}, y :: Vector{Float32}, z :: Vector{Float32} ) =
  coulomb(atoms, sel1, sel2, [0., 0., 0.], x, y, z, true, 12., false)

#
# For computing interaction energies from the coordinates available in the atoms array
#

function coulomb2( atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int}, usecutoff :: Bool, cutoff :: Float64 )
  n1 = length(sel1)
  n2 = length(sel2)
  coulomb = 0.
  for isel1 in 1:n1
    i = sel1[isel1]
    for isel2 in 1:n2
      j = sel2[isel2]
      d = sqrt( (atoms[i].coor[1] - atoms[j].coor[1])^2 +
                (atoms[i].coor[2] - atoms[j].coor[2])^2 +
                (atoms[i].coor[3] - atoms[j].coor[3])^2 )
      coulomb = coulomb + qpair(d,atoms[i].charge,atoms[j].charge)
      if usecutoff
        if d < cutoff
          coulomb = coulomb - qpair(cutoff,atoms[i].charge,atoms[j].charge)
        end
      end
    end
  end
  return coulomb
end

# Mapping to use a single function name

coulomb(  atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int}, usecutoff :: Bool, cutoff :: Float64 ) =
  coulomb2( atoms, sel1, sel2, usecutoff, cutoff )

# If the cutoff is not set, use default

coulomb( atoms :: Vector{Atom}, sel1 :: Vector{Int}, sel2 :: Vector{Int} ) =
  coulomb( atoms, sel1, sel2, true, 12. )


text:
using SchwarzChristoffel
using Test

using StaticArrays

@testset "integral" begin
    f = SchwarzChristoffelDerivate(SVector(0),SVector(0)); @test integral(f,1,1) ≈ 1
    f = SchwarzChristoffelDerivate(SVector(0),SVector(0)); @test integral(f,1,2) ≈ 2
    f = SchwarzChristoffelDerivate(SVector(1),SVector(0)); @test integral(f,1,2) ≈ 1

    f = SchwarzChristoffelDerivate(SVector(0),SVector(-1)); @test integral(f,1,1) ≈ 0.5
    f = SchwarzChristoffelDerivate(SVector(0),SVector(-1)); @test integral(f,1,2) ≈ 2
    f = SchwarzChristoffelDerivate(SVector(1),SVector(-1)); @test integral(f,1,2) ≈ 0.5

    f = SchwarzChristoffelDerivate(SVector(0),SVector(0.5)); @test integral(f,1,1) ≈ 2
    f = SchwarzChristoffelDerivate(SVector(0),SVector(0.5)); @test integral(f,1,2) ≈ 2*sqrt(2)
    f = SchwarzChristoffelDerivate(SVector(1),SVector(0.5)); @test integral(f,1,2) ≈ 2

    f = SchwarzChristoffelDerivate(SVector(-1,0),SVector(-1,-1)); @test integral(f,2,1) ≈ 5/6
    f = SchwarzChristoffelDerivate(SVector(-1,0),SVector(-1,-1)); @test integral(f,2,2) ≈ 14/3

    f = SchwarzChristoffelDerivate(SVector(-1e-8,0),SVector(0.5,0)); @test integral(f,2,1) ≈ (sqrt(1e8+1)-1)/5000

    f = SchwarzChristoffelDerivate(SVector(0,1),SVector(0,0)); @test segment(f,1) ≈ 1
end


text:
#### Precipitation component in atmosphere model
abstract type PrecipitationModel end

export NoPrecipitation, Rain

using ..Microphysics

vars_state_conservative(::PrecipitationModel, FT) = @vars()
vars_state_gradient(::PrecipitationModel, FT) = @vars()
vars_state_gradient_flux(::PrecipitationModel, FT) = @vars()
vars_state_auxiliary(::PrecipitationModel, FT) = @vars()

function atmos_nodal_update_auxiliary_state!(
    ::PrecipitationModel,
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function flux_precipitation!(
    ::PrecipitationModel,
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_flux!(
    ::PrecipitationModel,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
    ρD_t,
) end
function flux_second_order!(
    ::PrecipitationModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
) end
function flux_first_order!(
    ::PrecipitationModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_argument!(
    ::PrecipitationModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end

"""
    NoPrecipitation <: PrecipitationModel

No precipitation.
"""
struct NoPrecipitation <: PrecipitationModel end

"""
    Rain <: PrecipitationModel

Precipitation model with rain only.
"""
struct Rain <: PrecipitationModel end

vars_state_conservative(::Rain, FT) = @vars(ρq_rain::FT)
vars_state_gradient(::Rain, FT) = @vars(q_rain::FT)
vars_state_gradient_flux(::Rain, FT) = @vars(ρd_q_rain::SVector{3, FT})
vars_state_auxiliary(::Rain, FT) =
    @vars(terminal_velocity::FT, src_q_rai_tot::FT)

function atmos_nodal_update_auxiliary_state!(
    rain::Rain,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    q_rain = state.precipitation.ρq_rain / state.ρ
    if q_rain > FT(0) # TODO - need a way to prevent negative values
        aux.precipitation.terminal_velocity = terminal_velocity(q_rain, state.ρ)
    else
        aux.precipitation.terminal_velocity = FT(0)
    end

    ρ = state.ρ
    ρinv = 1 / state.ρ
    p = aux.pressure
    # TODO - ensure positive definite
    q_tot = max(FT(0), state.ρq_tot * ρinv)
    q_rai = max(FT(0), state.ρq_rain * ρinv)

    # current state
    ts = thermo_state(rain, state, aux)
    # q     = PhasePartition(q_tot, q_liq, q_ice)
    q = PhasePartition(ts)
    # T     = air_temperature(e_int, q)
    T = air_temperature(ts)
    # equilibrium state at current T
    q_eq = PhasePartition_equil(ts)

    # tendency from cloud water condensation/evaporation
    # src_q_liq = conv_q_vap_to_q_liq(q_eq, q)# TODO - temporary handling ice

    # tendencies from rain
    src_q_rai_evap = conv_q_rai_to_q_vap(atmos.param_set, q_rai, q, T, p, ρ)
    src_q_rai_acnv = conv_q_liq_to_q_rai_acnv(atmos.param_set, q.liq)
    src_q_rai_accr = conv_q_liq_to_q_rai_accr(atmos.param_set, q.liq, q_rai, ρ)

    aux.precipitation.src_q_rai_tot =
        src_q_rai_acnv + src_q_rai_accr + src_q_rai_evap
end

function compute_gradient_argument!(
    rain::Rain,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.precipitation.q_rain = state.precipitation.ρq_rain / state.ρ
end

function compute_gradient_flux!(
    rain::Rain,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    ρν::Union{Real, AbstractMatrix},
    D_T,
)
    # diffusive flux
    diffusive.precipitation.ρd_q_rain =
        state.ρ .* (-D_T) .* ∇transform.precipitation.q_rain
end

function flux_precipitation!(
    rain::Rain,
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    u = state.ρu / state.ρ
    flux.precipitation.ρq_rain += state.precipitation.ρq_rain * u
end

function flux_second_order!(
    rain::Rain,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    D_T,
)
    flux.precipitation.ρq_rain = diffusive.precipitation.ρd_q_rain
end

function boundarycondition_precipitation!(
    rain::Rain,
    stateP::Vars,
    diffP::Vars,
    auxP::Vars,
    nM,
    stateM::Vars,
    diffM::Vars,
    auxM::Vars,
    bctype::BC,
    t,
) where {BC}
    stateP.precipitation.ρq_rain = eltype(stateP)(0)
end


text:
if !(:ModelClockSystem in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "..", "share", "ModelClockSystem.jl")))
end

if !(:ConfigCheck in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "..", "share", "ConfigCheck.jl")))
end

if !(:CyclicData in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "..", "share", "CyclicData.jl")))
end

if !(:LogSystem in names(Main))
    include(normpath(joinpath(dirname(@__FILE__), "..", "..", "share", "LogSystem.jl")))
end

if ! ( :DataManager in names(Main) )
    include(joinpath(@__DIR__, "..", "..", "share", "DataManager.jl"))
end

macro hinclude(path)
    return :(include(normpath(joinpath(@__DIR__, $path))))
end

module EMOM

    using LinearAlgebra
    using MPI
    using Dates
    using Printf
    using Formatting
    using SharedArrays
    using Distributed
    using SparseArrays
    using NCDatasets
    using JLD2
    using DataStructures

    using ..ModelClockSystem
    using ..ConfigCheck
    using ..CyclicData
    using ..LogSystem
    using ..DataManager

    macro hinclude(path)
        return :(include(normpath(joinpath(@__DIR__, $path))))
    end

    @hinclude("../../share/constants.jl")
    @hinclude("../../share/ocean_state_function.jl")

    # classes
    @hinclude("../../share/GridFile.jl")
    @hinclude("../../share/PolelikeCoordinate.jl")
    @hinclude("../../share/BasicMatrixOperators.jl")
    @hinclude("../../share/AdvancedMatrixOperators.jl")

    @hinclude("entry_list.jl")

    @hinclude("VerticalDiffusion.jl")

    @hinclude("Topography.jl")
    @hinclude("Workspace.jl")
    @hinclude("Env.jl")
    @hinclude("TempField.jl")
    @hinclude("Field.jl")
    @hinclude("Core.jl")

    @hinclude("ModelBlock.jl")

    @hinclude("Leonard1979.jl")
    @hinclude("setupForcing.jl")
    @hinclude("stepAdvection.jl")
    @hinclude("stepColumn.jl")
    @hinclude("checkBudget.jl")

    @hinclude("var_list.jl")
    @hinclude("var_desc.jl")

    @hinclude("snapshot_funcs.jl")

    @hinclude("updateDatastream.jl")
    @hinclude("updateBuoyancy.jl")

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
include(joinpath(@__DIR__, "associations.jl"))
include(joinpath(@__DIR__, "glade_maker.jl"))
include(joinpath(@__DIR__, "util.jl"))

function wire_poi_gui(vfs, folder, add_poi)
    files = OrderedSet([vf.file for vf in vfs])
    points = strip.(vec(readcsv(joinpath(folder, "metadata", "poi.csv"), String)))

    poi_in = Signal(POI(points[1], Point(files[1], Second(0)), Point(files[1], Second(0)), "", ""))
    name_in = map(x -> x.name, poi_in)
    f1_in = map(x -> x.start.file, poi_in)
    f2_in = map(x -> x.stop.file, poi_in)
    start_time_in = map(x -> Time(0) + x.start.time, poi_in)
    stop_time_in = map(x -> Time(0) + x.stop.time, poi_in)
    poi_label_in = map(x -> x.label, poi_in)
    poi_comment_in = map(x -> x.comment, poi_in)

    poi_builder = Builder(filename=joinpath(@__DIR__, "poi.glade"))

    # widgets
    name = dropdown(points; widget=poi_builder["names"], signal=name_in)
    f1 = dropdown(files; widget=poi_builder["start.file"], signal=f1_in)
    f2 = dropdown(files; widget=poi_builder["stop.file"], signal=f2_in)
    start_time = timewidget(Time(0); widget=poi_builder["start.time"], signal=start_time_in)
    stop_time = timewidget(Time(0); widget=poi_builder["stop.time"], signal=stop_time_in)
    start = map(t -> Dates.Second(t - Dates.Time(0,0,0)), start_time)
    stop = map(t -> Dates.Second(t - Dates.Time(0,0,0)), stop_time)
    poi_label = textarea(;widget=poi_builder["label"], signal=poi_label_in)
    poi_comment = textarea(;widget=poi_builder["comment"], signal=poi_comment_in)
    playstart = button(;widget=poi_builder["start.play"])
    playstop = button(;widget=poi_builder["stop.play"])
    done_poi = togglebutton(true, widget=poi_builder["done"])
    cancel_poi = button(widget=poi_builder["cancel"])
    # functions
    foreach(playstart, init=nothing) do _
        @spawn openit(joinpath(folder, value(f1)))
        return nothing
    end
    foreach(playstop, init=nothing) do _
        @spawn openit(joinpath(folder, value(f2)))
        return nothing
    end
    #=tsksstrt, rsltsstrt = async_map(nothing, signal(playstart)) do _
        openit(joinpath(folder, value(f1)))
        return nothing
    end
    tsksstp, rsltsstp = async_map(nothing, signal(playstop)) do _
        openit(joinpath(folder, value(f2)))
        return nothing
    end=#
    start_point = map(Point, f1, start)
    stop_point = map(Point, f2, stop)
    foreach(start_point) do p
        if p.file == value(f2) && p.time > value(stop)
            push!(stop_time, Dates.Time(0,0,0) + p.time)
        end
    end
    foreach(stop_point) do p
        if p.file == value(f1) && p.time < value(start)
            push!(start_time, Dates.Time(0,0,0) + p.time)
        end
    end
    poi_temp = map(POI, name, start_point, stop_point, poi_label, poi_comment)

    poi_new = map(_ -> value(poi_temp), done_poi, init = value(poi_temp))
    poi_out = droprepeats(poi_new)

    foreach(cancel_poi) do _
        push!(poi_in, value(poi_out))
        push!(done_poi, true)
    end

    showall(poi_builder["window"])
    foreach(x -> visible(poi_builder["window"], !x), done_poi)

    bindmap!(signal(done_poi), !, signal(add_poi), !)

    return (poi_in, poi_out)
end

function wire_run_gui(folder, add_run)

    # prepare glade
    glade_widgets = OrderedDict{Symbol, Tuple{Symbol, String}}()
    drops = Dict{Symbol, Vector{String}}()
    tmp = readcsv(joinpath(folder, "metadata", "run.csv"))
    for i = 1:size(tmp,1)
        x = filter(!isempty, strip.(tmp[i,:]))
        l = x[1]
        id = Symbol(l)
        widget = length(x) == 1 ? :textbox : :dropdown
        glade_widgets[id] = (widget, l)
        if widget == :dropdown
            drops[id] = x[2:end]
        end
    end
    parse2glade(glade_widgets)
    run_builder = Builder(filename=joinpath(@__DIR__, "run.glade"))
    run_in = Signal(Run(Dict(k => haskey(drops, k) ? first(drops[k]) : "" for k in keys(glade_widgets)), ""))

    widgets = Dict{Symbol, Union{GtkReactive.Textarea, GtkReactive.Dropdown}}()
    for (id, (widget, _)) in glade_widgets
        if widget == :textbox
            widgets[id] = textarea(widget=run_builder[String(id)], signal = map(x -> x.metadata[id], run_in))
        else
            widgets[id] = dropdown(drops[id], widget=run_builder[String(id)], signal = map(x -> x.metadata[id], run_in))
        end
    end

    run_comment = textarea(widget=run_builder["comment.run.wJqRk"], signal = map(x -> x.comment, run_in))
    done_run = togglebutton(true, widget=run_builder["done.run.wJqRk"])
    cancel_run = button(widget=run_builder["cancel.run.wJqRk"])

    # function
    run_update = merge(map(signal, values(widgets))..., signal(run_comment))
    run_temp = map(_ -> Run(Dict(k => value(v) for (k,v) in widgets), value(run_comment)), run_update)

    new_run = map(_ -> value(run_temp), done_run, init = value(run_temp))

    counter = foldp(+, 1, signal(done_run))
    odd = map(isodd, counter)

    run_out = filterwhen(odd, value(new_run), new_run)

    foreach(cancel_run, init=nothing) do _
        r = value(run_out)
        for (k, v) in widgets
            push!(v, r.metadata[k])
        end
        push!(run_comment, r.comment)
        push!(done_run, true)
        a = value(association)
        delete!(a, last(a.runs))
        push!(association, a)
    end

    showall(run_builder["window.run.wJqRk"])
    foreach(x -> visible(run_builder["window.run.wJqRk"], !x), done_run)

    bindmap!(signal(done_run), !, signal(add_run), !)

    return (run_in, run_out)
end

function log_gui(folder)
    vfs = getVideoFiles(folder)
    # files = OrderedSet([vf.file for vf in vfs])
    # points = strip.(vec(readcsv(joinpath(folder, "metadata", "poi.csv"), String)))

    ##################################### LOG ######################################
    log_builder = Builder(filename=joinpath(@__DIR__, "log.glade"))

    # widgets
    add_poi = togglebutton(false, widget=log_builder["add.poi"])
    add_run = togglebutton(false, widget=log_builder["add.run"])
    clear = button(widget=log_builder["clear"])
    cancel = button(widget=log_builder["cancel"])
    done = button(widget=log_builder["done"])
    g = log_builder["poi.run.grid"]

    # functions
    foreach(x -> visible(log_builder["window"], !x), add_poi)
    foreach(x -> visible(log_builder["window"], !x), add_run)
    foreach(_ -> destroy(log_builder["window"]), cancel, init=nothing)

    showall(log_builder["window"])

    ##################################### POI ######################################

    poi_in, poi = wire_poi_gui(vfs, folder, add_poi)

    ##################################### RUN ######################################

    run_in, run = wire_run_gui(folder, add_run)

    ################################## ASSOCIATIONS ################################

    added = merge(poi, run)
    association = foldp(push!, loadAssociation(folder), added)

    merged_in = merge(poi_in, run_in)
    new_in = map(_ -> true, merged_in, init=false)
    replaceit = filterwhen(new_in, value(added), added)
    foreach(replaceit, init=nothing) do x
        a = value(association)
        replace!(a, value(merged_in), x)
        push!(new_in, false)
        push!(association, a)
    end

    assdone = map(association) do a
        empty!(g)
        g[0,0] = log_builder["poi.diagonal.run"]
        for (x, p) in enumerate(a.pois)
            file = MenuItem("_$(p.name) $(p.label)")
            filemenu = Menu(file)
            check_ = MenuItem("Check")
            checkh = signal_connect(check_, :activate) do _
                for r in a.runs
                    push!(a, (p, r))
                end
                push!(association, a)
            end
            push!(filemenu, check_)
            uncheck_ = MenuItem("Uncheck")
            uncheckh = signal_connect(uncheck_, :activate) do _
                for r in a.runs
                    delete!(a, (p, r))
                end
                push!(association, a)
            end
            push!(filemenu, uncheck_)
            #=hide_ = MenuItem("Hide")
            hideh = signal_connect(hide_, :activate) do _
                p.visible = false
                push!(association, a)
            end
            push!(filemenu, hide_)=#
            edit_ = MenuItem("Edit")
            edith = signal_connect(edit_, :activate) do _
                push!(add_poi, true)
                push!(poi_in, p)
                # foreach(add_poi, init=nothing) do _
                #=poi_ = poi_gui(p, points, files, folder)
                poi_new = droprepeats(poi_)
                foreach(poi_new, init = nothing) do n=#
                    # replace!(a, p, n)
                    # push!(association, a)
                    # nothing
                # end
            end
            push!(filemenu, edit_)
            push!(filemenu, SeparatorMenuItem())
            delete = MenuItem("Delete")
            deleteh = signal_connect(delete, :activate) do _
                delete!(a, p)
                push!(association, a)
            end
            push!(filemenu, delete)
            mb = MenuBar()
            push!(mb, file)
            g[x,0] = mb
        end
        g[length(a.pois) + 1,0] = togglebutton(false, widget=log_builder["add.poi"])
        for (y, r) in enumerate(a.runs)
            file = MenuItem(string("_", shorten(string(join(values(r.run.metadata), ":")..., ":", r.repetition), 30)))
            filemenu = Menu(file)
            check_ = MenuItem("Check")
            checkh = signal_connect(check_, :activate) do _
                for p in a.pois
                    push!(a, (p, r))
                end
                push!(association, a)
            end
            push!(filemenu, check_)
            uncheck_ = MenuItem("Uncheck")
            uncheckh = signal_connect(uncheck_, :activate) do _
                for p in a.pois
                    delete!(a, (p, r))
                end
                push!(association, a)
            end
            push!(filemenu, uncheck_)
            #=hide_ = MenuItem("Hide")
            hideh = signal_connect(hide_, :activate) do _
                r.visible = false
                push!(association, a)
            end
            push!(filemenu, hide_)=#
            edit_ = MenuItem("Edit")
            edith = signal_connect(edit_, :activate) do _
                push!(add_run, true)
                push!(run_in, r.run)
                #=run_ = run_gui(r.run, metadata)
                run_new = droprepeats(run_)
                foreach(run_new, init = nothing) do n
                    replace!(a, r, n)
                    push!(association, a)
                    nothing
                end=#
            end
            push!(filemenu, edit_)
            push!(filemenu, SeparatorMenuItem())
            delete = MenuItem("Delete")
            deleteh = signal_connect(delete, :activate) do _
                delete!(a, r)
                push!(association, a)
            end
            push!(filemenu, delete)
            mb = MenuBar()
            push!(mb, file)
            g[0,y] = mb
        end
        g[0, length(a.runs) + 1] = togglebutton(false, widget=log_builder["add.run"])
        for (x, p) in enumerate(a.pois), (y, r) in enumerate(a.runs)
            key = (p, r)
            cb = checkbox(key in a)
            foreach(cb) do tf
                tf ? push!(a, key) : delete!(a, key)
            end
            g[x,y] = cb
        end
        showall(log_builder["window"])
    end

    # clear the log window grid
    foreach(clear, init=nothing) do _
        a = value(association)
        empty!(a)
        push!(association, a)
        nothing
    end

    # done and move on to assesing the videos
    foreach(done, init=nothing) do _
        a = value(association)
        destroy(log_builder["window"])
        if !isempty(a)
            save(folder, a)
            checkvideos(a, folder, vfs)
        end
        nothing
    end
end

function return_selected_videos(a::Association, vfs::OrderedSet{VideoFile})
    uvfs = Set{String}(vf for poi in a.pois for vf in [poi.start.file, poi.stop.file])
    ft = Dict{String, VideoFile}()
    for vf in vfs
        if vf.file in uvfs
            ft[vf.file] = vf
        end
    end
    return collect(values(ft))
end

function checkvideos(a::Association, folder::String, vfs::OrderedSet{VideoFile})
    builder = Builder(filename=joinpath(@__DIR__, "video.glade"))
    # data
    ft = return_selected_videos(a, vfs)
    n = length(ft)
    # widgets
    done = button(widget=builder["done"])
    previous = button(widget=builder["previous"])
    next = button(widget=builder["next"])
    play = button(widget=builder["play"])
    # functions
    down = map(_ -> -1, previous)
    up = map(_ -> +1, next)
    step = merge(down, up)
    _state = foldp(1, step) do x,y
        clamp(x + y, 1, n)
    end
    state = droprepeats(_state)
    pb = progressbar(n, widget=builder["progressbar"], signal=state)
    file = map(state) do i
        ft[i].file
    end
    foreach(play, init=nothing) do _
        @spawn openit(joinpath(folder, value(file)))
        nothing
    end
    #=tsk, rslt = async_map(nothing, signal(play)) do _
        openit(joinpath(folder, value(file)))
        return nothing
    end=#
    datetime = map(state) do i
        ft[i].datetime
    end
    dt = datetimewidget(ft[1].datetime, widget=builder["datetime"], signal=datetime)
    label(ft[1].file, widget=builder["file.name"], signal=file)

    foreach(dt) do x
        i = value(state)
        ft[i] = VideoFile(ft[i].file, x)
    end
    foreach(done,  init = nothing) do _
        save(folder, OrderedSet{VideoFile}(ft))
        destroy(builder["window"])
    end
    showall(builder["window"])
end


text:
using ForwardDiff

"Parameters for right-hand side function of formal Lagrangian Runge-Kutta methods."
mutable struct ParametersFLRK{DT, TT, D, S, ET <: NamedTuple} <: Parameters{DT,TT}
    equs::ET
    tab::Tableau{TT}
    Δt::TT

    t::TT
    q::Vector{DT}
    p::Vector{DT}

    function ParametersFLRK{DT,D}(equs::ET, tab::Tableau{TT}, Δt::TT) where {DT, TT, D, ET <: NamedTuple}
        new{DT,TT,D,tab.s,ET}(equs, tab, Δt, zero(TT), zeros(DT,D), zeros(DT,D))
    end
end

struct IntegratorCacheFLRK{DT,D,S} <: ODEIntegratorCache{DT,D}
    q̃::Vector{DT}
    ṽ::Vector{DT}
    s̃::Vector{DT}

    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    Y::Vector{Vector{DT}}

    function IntegratorCacheFLRK{DT,D,S}() where {DT,D,S}
        # create internal stage vectors
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        Y = create_internal_stage_vector(DT, D, S)

        new(zeros(DT,D), zeros(DT,D), zeros(DT,D), Q, V, Y)
    end
end

function IntegratorCache{ST}(params::ParametersFLRK{DT,TT,D,S}; kwargs...) where {ST,DT,TT,D,S}
    IntegratorCacheFLRK{ST,D,S}(; kwargs...)
end

@inline CacheType(ST, params::ParametersFLRK{DT,TT,D,S}) where {DT,TT,D,S} = IntegratorCacheFLRK{ST,D,S}

"Formal Lagrangian Runge-Kutta integrator."
struct IntegratorFLRK{DT, TT, D, S, PT <: ParametersFLRK{DT,TT},
                                    ST <: NonlinearSolver{DT},
                                    IT <: InitialGuessODE{TT}} <: AbstractIntegratorPRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    ϑ::Vector{Vector{DT}}
    P::Vector{Vector{DT}}
    F::Vector{Vector{DT}}
    G::Vector{Vector{DT}}
    Z::Vector{Vector{DT}}
    J::Vector{Matrix{DT}}
    A::Array{DT,2}

    function IntegratorFLRK(params::ParametersFLRK{DT,TT,D,S}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,ST,IT}
        # create internal stage vectors
        ϑ = create_internal_stage_vector(DT, D, S)
        P = create_internal_stage_vector(DT, D, S)
        F = create_internal_stage_vector(DT, D, S)
        G = create_internal_stage_vector(DT, D, S)
        Z = create_internal_stage_vector(DT, D, S)

        J = [zeros(DT,D,D) for i in 1:S]
        A = zeros(DT,D*S,D*S)

        new{DT, TT, D, S, typeof(params), ST, IT}(params, solver, iguess, caches, ϑ, P, F, G, Z, J, A)
    end

    function IntegratorFLRK{DT,D}(equations::NamedTuple, tableau::Tableau{TT}, Δt::TT) where {DT,TT,D}
        # get number of stages
        S = tableau.s

        # check if tableau is fully implicit
        if get_config(:verbosity) ≥ 1
            if isexplicit(tableau)
                @warn "Initializing IntegratorFIRK with explicit tableau $(q.name).\nYou might want to use IntegratorERK instead."
            elseif isdiagnonallyimplicit(tableau)
                @warn "Initializing IntegratorFIRK with diagonally implicit tableau $(q.name).\nYou might want to use IntegratorDIRK instead."
            end
        end

        # create params
        params = ParametersFLRK{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, D*S, params, caches)

        # create initial guess
        iguess = InitialGuessODE(get_config(:ig_extrapolation), equations[:v̄], Δt)

        # create integrator
        IntegratorFLRK(params, solver, iguess, caches)
    end

    function IntegratorFLRK(equation::LODE{DT}, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT}
        IntegratorFLRK{DT, equation.d}(get_functions(equation), tableau, Δt; kwargs...)
    end

end

@inline Base.ndims(::IntegratorFLRK{DT,TT,D,S}) where {DT,TT,D,S} = D

function initialize!(int::IntegratorFLRK, sol::AtomicSolutionODE)
    sol.t̄ = sol.t - timestep(int)

    equations(int)[:v̄](sol.t, sol.q, sol.v)

    initialize!(int.iguess, sol.t, sol.q, sol.v,
                            sol.t̄, sol.q̄, sol.v̄)
end

function update_params!(int::IntegratorFLRK, sol::AtomicSolutionODE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
end

function update_params!(int::IntegratorFLRK, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
    int.params.p .= sol.p
end

function initial_guess!(int::IntegratorFLRK{DT}, sol::Union{AtomicSolutionODE{DT},AtomicSolutionPODE{DT}},
                        cache::IntegratorCacheFLRK{DT}=int.caches[DT]) where {DT}
    local offset::Int

    # compute initial guess for internal stages
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q̄, sol.v̄, sol.q, sol.v, cache.Q[i], cache.V[i], tableau(int).c[i])
    end
    for i in eachstage(int)
        offset = ndims(int)*(i-1)
        for k in eachdim(int)
            int.solver.x[offset+k] = 0
            for j in eachstage(int)
                int.solver.x[offset+k] += timestep(int) * tableau(int).a[i,j] * cache.V[j][k]
            end
        end
    end
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Y::Vector{Vector{ST}},
                         params::ParametersFLRK{DT,TT,D,S}) where {ST,DT,TT,D,S}

    local tᵢ::TT

    @assert S == length(Q) == length(V) == length(Y)

    # copy x to Y and compute Q = q + Δt Y
    for i in 1:S
        @assert D == length(Q[i]) == length(V[i]) == length(Y[i])
        for k in 1:D
            Y[i][k] = x[D*(i-1)+k]
            Q[i][k] = params.q[k] + Y[i][k]
        end
    end

    # compute V = v(Q)
    for i in 1:S
        tᵢ = params.t + params.Δt * params.tab.c[i]
        params.equs[:v̄](tᵢ, Q[i], V[i])
    end
end

# Compute stages of formal Lagrangian Runge-Kutta methods.
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersFLRK{DT,TT,D,S},
                          caches::CacheDict) where {ST,DT,TT,D,S}
    # temporary variables
    local y1::ST
    local y2::ST

    # get cache for internal stages
    cache = caches[ST]

    # compute stages from nonlinear solver solution x
    compute_stages!(x, cache.Q, cache.V, cache.Y, params)

    # compute b = - (Y-AV)
    for i in 1:S
        for k in 1:D
            y1 = 0
            y2 = 0
            for j in 1:S
                y1 += params.tab.a[i,j] * cache.V[j][k]
                y2 += params.tab.â[i,j] * cache.V[j][k]
            end
            b[D*(i-1)+k] = - cache.Y[i][k] + params.Δt * (y1 + y2)
        end
    end
end

function integrate_step!(int::IntegratorFLRK, sol::AtomicSolutionODE)
    integrate_step_flrk!(int, sol)
end

function integrate_step!(int::IntegratorFLRK, sol::AtomicSolutionPODE)
    integrate_step_flrk!(int, sol)
    integrate_diag_flrk!(int, sol)
end

function integrate_step_flrk!(int::IntegratorFLRK{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                              cache::IntegratorCacheFLRK{DT}=int.caches[DT]) where {DT,TT}
    # update nonlinear solver parameters from atomic solution
    update_params!(int, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset atomic solution
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector field at internal stages
    compute_stages!(int.solver.x, cache.Q, cache.V, cache.Y, int.params)

    # compute final update
    update_solution!(sol.q, sol.q̃, cache.V, tableau(int).b, tableau(int).b̂, timestep(int))

    # copy solution to initial guess
    update_vector_fields!(int.iguess, sol.t, sol.q, sol.v)
end

function integrate_diag_flrk!(int::IntegratorFLRK{DT,TT,D,S}, sol::AtomicSolutionPODE{DT,TT},
                              cache::IntegratorCacheFLRK{DT}=int.caches[DT]) where {DT,TT,D,S}
    # create temporary arrays
    δP = zeros(DT, D*S)

    # compute ϑ = ϑ(t,Q), V(Q) = v(t, Q, V)
    # and f_0(Q, V(Q)) = f(t, Q, V, F)
    for i in 1:S
        tᵢ = int.params.t + timestep(int) * tableau(int).c[i]
        int.params.equs[:ϑ](tᵢ, cache.Q[i], cache.V[i], int.ϑ[i])
        int.params.equs[:v̄](tᵢ, cache.Q[i], cache.V[i])
        int.params.equs[:g](tᵢ, cache.Q[i], cache.V[i], int.F[i])
    end

    # compute Jacobian of v via ForwardDiff
    for i in 1:S
        tᵢ = int.params.t + timestep(int) * tableau(int).c[i]
        v_rev! = (v,q) -> int.params.equs[:v̄](tᵢ,q,v)
        ForwardDiff.jacobian!(int.J[i], v_rev!, cache.ṽ, cache.Q[i])
    end

    # contract J with ϑ and add to G
    for l in 1:S
        for i in 1:D
            int.G[l][i] = 0
            for j in 1:D
                int.G[l][i] += int.ϑ[l][j] * int.J[l][j,i]
            end
        end
    end

    # solve linear system AP=δP for P

    # compute δP
    for l in 1:S
        for i in 1:D
            # set δP = p
            δP[(l-1)*D+i] = int.params.p[i]
            # add A(F+G) to δP
            for k in 1:S
                δP[(l-1)*D+i] += timestep(int) * tableau(int).a[l,k] * (int.F[k][i] + int.G[k][i])
            end
        end
    end

    # construct A = identity(sd×sd) + A ⊗ J
    for k in 1:S
        for l in 1:S
            for i in 1:D
                for j in 1:D
                    int.A[(k-1)*D+i, (l-1)*D+j] = timestep(int) * tableau(int).a[k,l] * int.J[l][j,i]
                end
            end
        end
    end
    for i in 1:D*S
        int.A[i,i] += one(DT)
    end

    # solve AP = δP and copy result to int.P
    lu = LUSolver(int.A, δP)
    factorize!(lu)
    solve!(lu)
    tP = reshape(lu.x, (D,S))
    for l in 1:S
        int.P[l] .= tP[:,l]
    end

    # contract J with P and subtract from G, so that G = (ϑ-P)J
    for l in 1:S
        for i in 1:D
            int.G[l][i] = 0
            for j in 1:D
                int.G[l][i] += (int.ϑ[l][j] - int.P[l][j]) * int.J[l][j,i]
            end
        end
    end

    # println(int.ϑ)
    # println(int.P)
    # println(int.ϑ .- int.P)
    # println()

    # compute final update for p
    update_solution!(sol.p, sol.p̃, int.F, tableau(int).b, tableau(int).b̂, timestep(int))
    update_solution!(sol.p, sol.p̃, int.G, tableau(int).b, tableau(int).b̂, timestep(int))
end


text:
function is_leap_year(year::Int)
    (0 == year % 4) && ((0 != year % 100) || (0 == year % 400))
end

max_stars_repo_name: cjtallman/exercism

max_issues_repo_name: cjtallman/exercism

max_forks_repo_name: cjtallman/exercism


text:
using NumericalMethodsforEngineers

a = [16. 4. 8.; 4. 5. -4.; 8. -4. 22.]
b = [4., 2., 5.]

f = cholesky(a)
f.U |> display

y = f.L \ b
c = f.U \ y

maxiters = 5
x = [(i=i, cg=cg(a, b; maxiter=i)) for i in 1:maxiters]
x |> display

@test round.(x[3].cg; digits=9) == round.(c; digits=9)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetSkuCapacity <: SwaggerModel
    minimum::Any # spec type: Union{ Nothing, Int64 } # spec name: minimum
    maximum::Any # spec type: Union{ Nothing, Int64 } # spec name: maximum
    defaultCapacity::Any # spec type: Union{ Nothing, Int64 } # spec name: defaultCapacity
    scaleType::Any # spec type: Union{ Nothing, String } # spec name: scaleType

    function VirtualMachineScaleSetSkuCapacity(;minimum=nothing, maximum=nothing, defaultCapacity=nothing, scaleType=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetSkuCapacity, Symbol("minimum"), minimum)
        setfield!(o, Symbol("minimum"), minimum)
        validate_property(VirtualMachineScaleSetSkuCapacity, Symbol("maximum"), maximum)
        setfield!(o, Symbol("maximum"), maximum)
        validate_property(VirtualMachineScaleSetSkuCapacity, Symbol("defaultCapacity"), defaultCapacity)
        setfield!(o, Symbol("defaultCapacity"), defaultCapacity)
        validate_property(VirtualMachineScaleSetSkuCapacity, Symbol("scaleType"), scaleType)
        setfield!(o, Symbol("scaleType"), scaleType)
        o
    end
end # type VirtualMachineScaleSetSkuCapacity

const _property_map_VirtualMachineScaleSetSkuCapacity = Dict{Symbol,Symbol}(Symbol("minimum")=>Symbol("minimum"), Symbol("maximum")=>Symbol("maximum"), Symbol("defaultCapacity")=>Symbol("defaultCapacity"), Symbol("scaleType")=>Symbol("scaleType"))
const _property_types_VirtualMachineScaleSetSkuCapacity = Dict{Symbol,String}(Symbol("minimum")=>"Int64", Symbol("maximum")=>"Int64", Symbol("defaultCapacity")=>"Int64", Symbol("scaleType")=>"String")
Base.propertynames(::Type{ VirtualMachineScaleSetSkuCapacity }) = collect(keys(_property_map_VirtualMachineScaleSetSkuCapacity))
Swagger.property_type(::Type{ VirtualMachineScaleSetSkuCapacity }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineScaleSetSkuCapacity[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetSkuCapacity }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetSkuCapacity[property_name]

const _allowed_VirtualMachineScaleSetSkuCapacity_scaleType = ["Automatic", "None"]

function check_required(o::VirtualMachineScaleSetSkuCapacity)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetSkuCapacity }, name::Symbol, val)
    if name === Symbol("scaleType")
        Swagger.validate_param(name, "VirtualMachineScaleSetSkuCapacity", :enum, val, _allowed_VirtualMachineScaleSetSkuCapacity_scaleType)
    end
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI
import FMIZoo

# our simulation setup
t_start = 0.0
t_stop = 8.0

# we use a FMU from the FMIZoo.jl
pathToFMU = FMIZoo.get_model_filename("IO", "Dymola", "2022x")

myFMU = fmiLoad(pathToFMU)
fmiInstantiate!(myFMU; loggingOn=true)
fmiSetupExperiment(myFMU, 0.0)
fmiEnterInitializationMode(myFMU)

fmiGetString(myFMU, "p_string")

rndReal = 100 * rand()
rndInteger = round(Integer, 100 * rand())
rndBoolean = rand() > 0.5
rndString = "Not random!"

# case A: The fast way ...

fmiSet(myFMU, ["p_real", "p_integer", "p_boolean", "p_string"], [rndReal, rndInteger, rndBoolean, rndString])
fmiGet(myFMU, ["p_real", "p_integer", "p_boolean", "p_string"])

# case B: Maximum control over what happens ...

fmiSetReal(myFMU, "p_real", rndReal)
display("$rndReal == $(fmiGetReal(myFMU, "p_real"))")

fmiSetInteger(myFMU, "p_integer", rndInteger)
display("$rndInteger == $(fmiGetInteger(myFMU, "p_integer"))")

fmiSetBoolean(myFMU, "p_boolean", rndBoolean)
display("$rndBoolean == $(fmiGetBoolean(myFMU, "p_boolean"))")

fmiSetString(myFMU, "p_string", rndString)
display("$rndString == $(fmiGetString(myFMU, "p_string"))")

# clean up

fmiExitInitializationMode(myFMU)

fmiUnload(myFMU)


text:
########################################################################
#                               API                                    #
########################################################################

const DIFF_RULES = Vector{Tuple}()
const NO_DIFF_RULES = Vector{Tuple}()
const DIFF_PHS = Set([:x, :u, :v, :w, :t, :i, :j, :k,])

function resolve_old_broadcast(ex)
    # rewrite dot op symbols like .*, .+, etc. into broadcasting
    for (pat, rpat) in [
        :(.+(_xs...)) => :($broadcast($+, _xs...)),
        :(.-(_xs...)) => :($broadcast($-, _xs...)),
        :(.*(_xs...)) => :($broadcast($*, _xs...)),
        :(./(_xs...)) => :($broadcast($/, _xs...)),
    ]
        ex = Espresso.rewrite_all(ex, pat, rpat)
    end
    return ex
end

function resolve_functions_and_types!(mod::Module, ex)
    if Meta.isexpr(ex, :call)
        # replace function symbol with actual function reference and
        # recursively call on function args
        if string(ex.args[1])[1] != '.'  # .*, .+, etc.
            ex.args[1] = Core.eval(mod, ex.args[1])
        end
        resolve_functions_and_types!(mod, @view ex.args[2:end])
        # for x in ex.args[2:end]
        #     resolve_functions_and_types!(mod, x)
        # end
    elseif Meta.isexpr(ex, :(::))
        ex.args[2] = Core.eval(mod, ex.args[2])
    elseif ex isa Vector || ex isa SubArray
        for (i, x) in enumerate(ex)
            if Meta.isexpr(x, :$)
                ex[i] = Core.eval(mod, x.args[1])
            else
                resolve_functions_and_types!(mod, x)
            end
        end
    end
    return ex
end

add_diff_rule(rule) = push!(DIFF_RULES, rule)
add_primitive(primitive) = push!(PRIMITIVES, primitive)

macro diffrule(pat, var, rpat)
    esc(quote
        mod = @__MODULE__
        local pat, rpat = map($Espresso.sanitize, ($(QuoteNode(pat)), $(QuoteNode(rpat))))
        pat, rpat = map($resolve_old_broadcast, (pat, rpat))
        $resolve_functions_and_types!(mod, pat)
        $resolve_functions_and_types!(mod, rpat)
        $add_diff_rule((pat, $(QuoteNode(var)), rpat))
        $add_primitive(pat.args[1])
        nothing
        end)
end

macro diffrule_kw(pat, var, rpat)
    kw_pat = rewrite(pat, :(_fn(_args...)), :(Core.kwfunc(_fn)(_kw, _, _args...)))
    kw_rpat = rewrite(rpat, :(_fn(_args...)), :(Core.kwfunc(_fn)(_kw, _, _args...)))
    kw_rpat = subs(kw_rpat, Dict(:_ => Expr(:$, rpat.args[1])))
    return esc(
        quote
        @diffrule $pat $var $rpat
        @diffrule $kw_pat $var $kw_rpat
        @nodiff $kw_pat _kw
        @nodiff $kw_pat _
        end)
end

isparameters(a) = isa(a, Expr) && a.head == :parameters

function without_types(pat)
    rpat = copy(pat)
    for i=2:length(rpat.args)
        a = rpat.args[i]
        if !isparameters(a)  # parameters aren't modified
            rpat.args[i] = isa(a, Expr) ? a.args[1] : a
        end
    end
    return rpat
end

function get_arg_names(pat)
    return [isa(a, Expr) ? a.args[1] : a for a in pat.args[2:end] if !isparameters(a)]
end

function get_arg_types(pat)
    return [isa(a, Expr) ? Core.eval(Base, a.args[2]) : Any
            for a in pat.args[2:end] if !isparameters(a)]
end

function match_rule(rule, ex, dep_types, idx)
    tpat, vname, rpat = rule
    vidx = findfirst(isequal(vname), get_arg_names(tpat))
    if idx != vidx
        return nothing
    end
    pat_types = get_arg_types(tpat)
    if (length(dep_types) != length(pat_types) ||
        !all([t <: p for (t, p) in zip(dep_types, pat_types)]))
        return nothing
    end
    pat = without_types(tpat)
    ex_ = without_keywords(ex)
    if !matchingex(pat, ex_; phs=DIFF_PHS)
        return nothing
    else
        return pat => rpat
    end
end

function rewrite_with_keywords(ex, pat, rpat)
    if rpat isa Expr && rpat.head == :call
        op, args, kw_args = parse_call_expr(ex)
        ex_no_kw = Expr(:call, op, args...)
        rex_no_kw = rewrite(ex_no_kw, pat, rpat; phs=DIFF_PHS)
        rex = with_keywords(rex_no_kw, kw_args)
    else
        rex = rewrite(ex, pat, rpat; phs=DIFF_PHS)
    end
    return rex
end

"""
Rewrite single call expression into its derivative. Example:

```
deriv_expr(:($sin(x)), [Float64], 1)
# ==> :((cos(x) * ds))
```
"""
function deriv_expr(ex, dep_types, idx::Int)
    rex = nothing
    for rule in DIFF_RULES
        m = match_rule(rule, ex, dep_types, idx)
        if m != nothing
            pat, rpat = m
            rex = rewrite_with_keywords(ex, pat, rpat)
            break
        end
    end
    if rex == nothing
        error("Can't find differentiation rule for $ex at $idx " *
              "with types $dep_types)")
    end
    return rex
end

"""
Internal function for finding rule by function name
"""
function find_rules_for(fun)
    return [r for r in DIFF_RULES if r[1].args[1] == fun]
end

## nodiff

add_no_diff_rule(rule) = push!(NO_DIFF_RULES, rule)

macro nodiff(pat, var)
    esc(quote
        mod = @__MODULE__
        local pat = $Espresso.sanitize($(QuoteNode(pat)))
        pat, $resolve_old_broadcast(pat)
        $resolve_functions_and_types!(mod, pat)
        $add_no_diff_rule((pat, $(QuoteNode(var))))
        $add_primitive(pat.args[1])
        nothing
        end)
end

function match_nodiff_rule(rule, ex, dep_types, idx)
    tpat, vname = rule
    vidx = findfirst(isequal(vname), get_arg_names(tpat))
    if idx != vidx
        return false
    end
    pat_types = get_arg_types(tpat)
    if (length(dep_types) != length(pat_types) ||
        !all([t <: p for (t, p) in zip(dep_types, pat_types)]))
        return false
    end
    pat = without_types(tpat)
    ex_ = without_keywords(ex)
    return matchingex(pat, ex_; phs=DIFF_PHS)
end

function dont_diff(tape::Tape, op::AbstractOp, idx::Int)
    if op.fn == broadcast
        ex = to_unbroadcast_expr(tape, op)
        dep_types = [eltype(tape[arg].typ) for arg in op.args[2:end]]
        idx_ = idx - 1
    else
        ex =  to_expr(tape, op)
        dep_types = [tape[arg].typ for arg in op.args]
        idx_ = idx
    end
    for rule in NO_DIFF_RULES
        if match_nodiff_rule(rule, ex, dep_types, idx_)
            return true
        end
    end
    return false
end

include("basic.jl")


text:
using Statistics

function fill_refs!(refs::AbstractArray, X::AbstractArray,
                    breaks::AbstractVector, extend::Bool, allow_missing::Bool)
    n = length(breaks)
    lower = first(breaks)
    upper = last(breaks)

    @inbounds for i in eachindex(X)
        x = X[i]

        if extend && x == upper
            refs[i] = n-1
        elseif !extend && !(lower <= x < upper)
            throw(ArgumentError("value $x (at index $i) does not fall inside the breaks: adapt them manually, or pass extend=true"))
        else
            refs[i] = searchsortedlast(breaks, x)
        end
    end
end

function fill_refs!(refs::AbstractArray, X::AbstractArray{>: Missing},
                    breaks::AbstractVector, extend::Bool, allow_missing::Bool)
    n = length(breaks)
    lower = first(breaks)
    upper = last(breaks)

    @inbounds for i in eachindex(X)
        ismissing(X[i]) && continue

        x = X[i]

        if extend && x == upper
            refs[i] = n-1
        elseif !extend && !(lower <= x < upper)
            allow_missing || throw(ArgumentError("value $x (at index $i) does not fall inside the breaks: adapt them manually, or pass extend=true or allow_missing=true"))
            refs[i] = 0
        else
            refs[i] = searchsortedlast(breaks, x)
        end
    end
end

"""
    default_formatter(from, to, i; closed=false)

Provide the default label format for the `cut` function.
"""
default_formatter(from, to, i; closed) = string("[", from, ", ", to, closed ? "]" : ")")

@doc raw"""
    cut(x::AbstractArray, breaks::AbstractVector;
        labels::Union{AbstractVector{<:AbstractString},Function},
        extend::Bool=false, allow_missing::Bool=false)

Cut a numeric array into intervals and return an ordered `CategoricalArray` indicating
the interval into which each entry falls. Intervals are of the form `[lower, upper)`,
i.e. the lower bound is included and the upper bound is excluded.

If `x` accepts missing values (i.e. `eltype(x) >: Missing`) the returned array will
also accept them.

# Arguments
* `extend::Bool=false`: when `false`, an error is raised if some values in `x` fall
  outside of the breaks; when `true`, breaks are automatically added to include all
  values in `x`, and the upper bound is included in the last interval.
* `labels::Union{AbstractVector,Function}: a vector of strings giving the names to use for
  the intervals; or a function `f(from, to, i; closed)` that generates the labels from the
  left and right interval boundaries and the group index. Defaults to
  `"[from, to)"` (or `"[from, to]"` for the rightmost interval if `extend == true`).
* `allow_missing::Bool=true`: when `true`, values outside of breaks result in missing values.
  only supported when `x` accepts missing values.

# Examples
```jldoctest
julia> cut(-1:0.5:1, [0, 1], extend=true)
5-element CategoricalArray{String,1,UInt32}:
 "[-1.0, 0.0)"
 "[-1.0, 0.0)"
 "[0.0, 1.0]"
 "[0.0, 1.0]"
 "[0.0, 1.0]"

julia> cut(-1:0.5:1, 2)
5-element CategoricalArray{String,1,UInt32}:
 "[-1.0, 0.0)"
 "[-1.0, 0.0)"
 "[0.0, 1.0]"
 "[0.0, 1.0]"
 "[0.0, 1.0]"

julia> cut(-1:0.5:1, 2, labels=["A", "B"])
5-element CategoricalArray{String,1,UInt32}:
 "A"
 "A"
 "B"
 "B"
 "B"

julia> fmt(from, to, i; closed) = "grp $i ($from//$to)"
fmt (generic function with 1 method)

julia> cut(-1:0.5:1, 3, labels=fmt)
5-element CategoricalArray{String,1,UInt32}:
 "grp 1 (-1.0//-0.333333)"
 "grp 1 (-1.0//-0.333333)"
 "grp 2 (-0.333333//0.333333)"
 "grp 3 (0.333333//1.0)"
 "grp 3 (0.333333//1.0)"
```
"""
function cut(x::AbstractArray{T, N}, breaks::AbstractVector;
             extend::Bool=false,
             labels::Union{AbstractVector{<:AbstractString},Function}=default_formatter,
             allow_missing::Bool=false) where {T, N}
    if !issorted(breaks)
        breaks = sort(breaks)
    end

    if extend
        min_x, max_x = extrema(x)
        if !ismissing(min_x) && breaks[1] > min_x
            breaks = [min_x; breaks]
        end
        if !ismissing(max_x) && breaks[end] < max_x
            breaks = [breaks; max_x]
        end
    end

    refs = Array{DefaultRefType, N}(undef, size(x))
    try
        fill_refs!(refs, x, breaks, extend, allow_missing)
    catch err
        # So that the error appears to come from cut() itself,
        # since it refers to its keyword arguments
        if isa(err, ArgumentError)
            throw(err)
        else
            rethrow(err)
        end
    end

    n = length(breaks)
    if labels isa Function
        @static if VERSION >= v"0.7.0-DEV.4524"
            from = map(x -> sprint(show, x, context=:compact=>true), breaks[1:n-1])
            to = map(x -> sprint(show, x, context=:compact=>true), breaks[2:n])
        else
            from = map(x -> sprint(showcompact, x), breaks[1:n-1])
            to = map(x -> sprint(showcompact, x), breaks[2:n])
        end
        levs = Vector{String}(undef, n-1)
        for i in 1:n-2
            levs[i] = labels(from[i], to[i], i, closed=false)
        end
        levs[end] = labels(from[end], to[end], n-1, closed=extend)
    else
        length(labels) == n-1 || throw(ArgumentError("labels must be of length $(n-1), but got length $(length(labels))"))
        # Levels must have element type String for type stability of the result
        levs::Vector{String} = copy(labels)
    end

    pool = CategoricalPool(levs, true)
    S = T >: Missing ? Union{String, Missing} : String
    CategoricalArray{S, N}(refs, pool)
end

"""
    cut(x::AbstractArray, ngroups::Integer;
        labels::Union{AbstractVector{<:AbstractString},Function})

Cut a numeric array into `ngroups` quantiles, determined using
[`quantile`](@ref).
"""
cut(x::AbstractArray, ngroups::Integer;
    labels::Union{AbstractVector{<:AbstractString},Function}=default_formatter) =
    cut(x, Statistics.quantile(x, (1:ngroups-1)/ngroups); extend=true, labels=labels)


text:
using Stheno: CustomMean, ZeroMean, OneMean, ConstMean

@timedtestset "mean" begin
    @timedtestset "CustomMean" begin
        rng, N, D = MersenneTwister(123456), 11, 2
        x = randn(rng, N)
        foo_mean = x->sum(abs2, x)
        f = CustomMean(foo_mean)

        @test ew(f, x) == map(foo_mean, x)
        differentiable_mean_function_tests(f, randn(rng, N), x)
    end
    @timedtestset "ZeroMean" begin
        rng, P, Q, D = MersenneTwister(123456), 3, 2, 4
        X, x = ColVecs(randn(rng, D, P)), randn(rng, P)
        f = ZeroMean{Float64}()

        for x in [x, X]
            @test ew(f, x) == zeros(size(x))
            differentiable_mean_function_tests(f, randn(rng, P), x)
        end
    end
    @timedtestset "OneMean" begin
        rng, P, Q, D = MersenneTwister(123456), 3, 2, 4
        X, x = ColVecs(randn(rng, D, P)), randn(rng, P)
        f = OneMean()

        for x in [x, X]
            @test ew(f, x) == ones(size(x))
            differentiable_mean_function_tests(f, randn(rng, P), x)
        end
    end
    @timedtestset "ConstMean" begin
        rng, D, N = MersenneTwister(123456), 5, 3
        X, x, c = ColVecs(randn(rng, D, N)), randn(rng, N), randn(rng)
        m = ConstMean(c)

        for x in [x, X]
            @test ew(m, x) == fill(c, N)
            differentiable_mean_function_tests(m, randn(rng, N), x)
        end
    end
    @timedtestset "(is)zero" begin
        @test zero(ZeroMean()) == ZeroMean()
        @test zero(OneMean()) == ZeroMean()
        @test iszero(ZeroMean()) == true
        @test iszero(OneMean()) == false
    end
end


text:
using Intan, SpikeSorting, Gtk.ShortNames

myt=Task_NoTask()
myfpga=FPGA(1,[0],usb3=true)
myparams=Algorithm[DetectAbs(),ClusterTemplate(49),AlignProm(),FeatureTime(),ReductionNone(),ThresholdMeanN()];
(myrhd,ss,myfpgas)=makeRHD([myfpga],params=myparams,single_channel_mode=true);

handles = makegui(myrhd,ss,myt,myfpgas);

max_stars_repo_name: paulmthompson/WhiskerTask

max_issues_repo_name: paulmthompson/WhiskerTask

max_forks_repo_name: paulmthompson/WhiskerTask


text:
module Genia

using LightXML

type Token
    word
    cat
end

Token(xelem::XMLElement) = Token(content(xelem), attribute(xelem, "c"))

type Sentence
    tokens::Vector{Token}
end

Sentence(xelem::XMLElement) = Sentence(map(Token, get_elements_by_tagname(xelem, "w")))

type Article
    bibliomisc
    title::Sentence
    abst::Vector{Sentence}
end

Base.show(io::IO, a::Article) = println(io, "Article($(a.title))")

function Article(xnode::XMLElement)
    xbib = get_elements_by_tagname(xnode, "articleinfo")[1]
    xbib = get_elements_by_tagname(xbib, "bibliomisc")[1]
    xtitle = get_elements_by_tagname(xnode, "title")[1]
    xtitle = get_elements_by_tagname(xtitle, "sentence")[1]
    xabst = get_elements_by_tagname(xnode, "abstract")[1]
    xabst = get_elements_by_tagname(xabst, "sentence")
    Article(content(xbib), Sentence(xtitle), Sentence[])
end

function readxml(path)
    xdoc = parse_file(path)
    xroot = root(xdoc)
    doc = map(Article, get_elements_by_tagname(xroot, "article"))
    doc
end

end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
using BinaryBuilder

# Collection of sources required to build Arpack
name = "Arpack"
version = v"3.5.1" # <-- This is actually v3.5.0, but we need to build for new platforms

sources = [
    GitSource("https://github.com/opencollab/arpack-ng.git",
              "9233f7f86f063ca6ca3793cb54dec590eb146e10"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/arpack-ng*

# arpack tests require finding libgfortran when linking with C linkers,
# and gcc doesn't automatically add that search path.  So we do it for it with `rpath-link`.
EXE_LINK_FLAGS=()
if [[ ${target} != *darwin* ]]; then
    EXE_LINK_FLAGS+=("-Wl,-rpath-link,/opt/${target}/${target}/lib")
    EXE_LINK_FLAGS+=("-Wl,-rpath-link,/opt/${target}/${target}/lib64")
fi

# Symbols that have float32, float64, complexf32, and complexf64 support
SDCZ_SYMBOLS=(
    axpy copy gemv geqr2 lacpy lahqr lanhs larnv lartg
    lascl laset scal trevc trmm trsen gbmv gbtrf gbtrs
    gttrf gttrs pttrf pttrs
)

# All symbols that have float32/float64 support (including the SDCZ_SYMBOLS above)
SD_SYMBOLS=(
    ${SDCZ_SYMBOLS[@]}
    dot ger labad laev2 lamch lanst lanv2
    lapy2 larf larfg lasr nrm2 orm2r rot steqr swap
)

# All symbols that have complexf32/complexf64 support (including the SDCZ_SYMBOLS above)
CZ_SYMBOLS=(${SDCZ_SYMBOLS[@]} dotc geru unm2r)

# Add in (s|d)*_64 symbol remappings:
SYMBOL_DEFS=()
for sym in ${SD_SYMBOLS[@]}; do
    SYMBOL_DEFS+=("-Ds${sym}=s${sym}_64" "-Dd${sym}=d${sym}_64")
done

# Add in (c|z)*_64 symbol remappings:
for sym in ${CZ_SYMBOLS[@]}; do
    SYMBOL_DEFS+=("-Dc${sym}=c${sym}_64" "-Dz${sym}=z${sym}_64")
done

# Add one-off symbol mappings; things that don't fit into any other bucket:
for sym in scnrm2 dznrm2 csscal zdscal dgetrf dgetrs; do
    SYMBOL_DEFS+=("-D${sym}=${sym}_64")
done

# Set up not only lowercase symbol remappings, but uppercase as well:
SYMBOL_DEFS+=(${SYMBOL_DEFS[@]^^})

FFLAGS="${FFLAGS} -O3 -fPIE -ffixed-line-length-none -fno-optimize-sibling-calls -cpp"
LIBOPENBLAS=openblas
if [[ ${nbits} == 64 ]]; then
    LIBOPENBLAS=openblas64_
    FFLAGS="${FFLAGS} -fdefault-integer-8 ${SYMBOL_DEFS[@]}"
fi

# Work around error
#
#     Error: Rank mismatch between actual argument at (1) and actual argument at (2) (scalar and rank-1)
#
# Properly fixed upstream in v3.8.0 with https://github.com/opencollab/arpack-ng/pull/245.
# TODO: Remove this line when we upgrade to that version.
if [[ "${target}" == aarch64-apple-* ]]; then
    FFLAGS="${FFLAGS} -fallow-argument-mismatch"
fi

mkdir build
cd build
export LDFLAGS="${EXE_LINK_FLAGS[@]} -L$prefix/lib -lpthread"
cmake .. -DCMAKE_INSTALL_PREFIX="$prefix" \
    -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" -DCMAKE_BUILD_TYPE=Release \
    -DEXAMPLES=OFF \
    -DBUILD_SHARED_LIBS=ON \
    -DBLAS_LIBRARIES="-l${LIBOPENBLAS}" \
    -DLAPACK_LIBRARIES="-l${LIBOPENBLAS}" \
    -DCMAKE_Fortran_FLAGS="${FFLAGS}"

make -j${nproc} VERBOSE=1
make install VERBOSE=1

# Arpack links against a _very_ specific version of OpenBLAS on macOS by default:
if [[ ${target} == *apple* ]]; then
    # Figure out what version it probably latched on to:
    OPENBLAS_LINK=$(otool -L ${prefix}/lib/libarpack.dylib | grep libopenblas64_ | awk '{ print $1 }')
    install_name_tool -change ${OPENBLAS_LINK} @rpath/libopenblas64_.dylib ${prefix}/lib/libarpack.dylib
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We enable the full
# combinatorial explosion of GCC versions because this package most
# definitely links against libgfortran.
platforms = expand_gfortran_versions(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libarpack", :libarpack),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("OpenBLAS_jll"),
    Dependency("CompilerSupportLibraries_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"6", julia_compat="1.7")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
#
#   SVD decomposition of a matrix truncated to a rank
#
using Manopt, Manifolds
import LinearAlgebra: norm, svd, Diagonal

"""
    truncated_svd(A, p)

return a singular value decomposition of a real valued matrix A truncated to
rank p.

# Input
* `A` – a real-valued matrix A of size mxn
* `p` – an integer p in min { m, n }

# Output
* `U` – an orthonormal matrix of size mxp
* `V` – an orthonormal matrix of size nxp
* `S` – a diagonal matrix of size pxp with nonnegative and decreasing diagonal
        entries
"""
function truncated_svd(A::Array{Float64,2}=randn(42, 60), p::Int64=5)
    (m, n) = size(A)

    if p > min(m, n)
        throw(
            ErrorException(
                "The Rank p=$p must be smaller than the smallest dimension of A = $min(m, n).",
            ),
        )
    end

    M = ProductManifold(Grassmann(m, p), Grassmann(n, p))

    function cost(X::ProductRepr)
        return cost([submanifold_components(X)...])
    end
    function cost(X::Array{Matrix{Float64},1})
        return -0.5 * norm(transpose(X[1]) * A * X[2])^2
    end

    function egrad(X::Array{Matrix{Float64},1})
        U = X[1]
        V = X[2]
        AV = A * V
        AtU = transpose(A) * U
        return [-AV * (transpose(AV) * U), -AtU * (transpose(AtU) * V)]
    end

    function rgrad(M::ProductManifold, X::ProductRepr)
        eG = egrad([submanifold_components(M, X)...])
        x = [submanifold_components(M, X)...]
        return Manifolds.ProductRepr(project.(M.manifolds, x, eG)...)
    end

    function e2rHess(
        M::Grassmann, x, ξ, eGrad::Matrix{T}, Hess::Matrix{T}
    ) where {T<:Union{U,Complex{U}}} where {U<:AbstractFloat}
        pxHess = project(M, x, Hess)
        xtGrad = x' * eGrad
        ξxtGrad = ξ * xtGrad
        return pxHess - ξxtGrad
    end

    function eHess(X::Array{Matrix{Float64},1}, H::Array{Matrix{Float64},1})
        U = X[1]
        V = X[2]
        Udot = H[1]
        Vdot = H[2]
        AV = A * V
        AtU = transpose(A) * U
        AVdot = A * Vdot
        AtUdot = transpose(A) * Udot
        return [
            -(
                AVdot * transpose(AV) * U +
                AV * transpose(AVdot) * U +
                AV * transpose(AV) * Udot
            ),
            -(
                AtUdot * transpose(AtU) * V +
                AtU * transpose(AtUdot) * V +
                AtU * transpose(AtU) * Vdot
            ),
        ]
    end

    function rhess(M::ProductManifold, X::ProductRepr, H::ProductRepr)
        x = [submanifold_components(M, X)...]
        h = [submanifold_components(M, H)...]
        eG = egrad(x)
        eH = eHess(x, h)
        return Manifolds.ProductRepr(e2rHess.(M.manifolds, x, h, eG, eH)...)
    end

    x = random_point(M)
    X = trust_regions(
        M,
        cost,
        rgrad,
        x,
        rhess;
        Δ_bar=4 * sqrt(2 * p),
        debug=[:Iteration, " ", :Cost, " | ", DebugEntry(:Δ), "\n", 1, :Stop],
    )
    U = X[M, 1]
    V = X[M, 2]
    Spp = transpose(U) * A * V
    SVD = svd(Spp)
    U = U * SVD.U
    S = SVD.S
    V = V * SVD.V

    return [U, S, V]
end

A = [1.0 2.0 3.0 4.0; 5.0 6.0 7.0 8.0; 9.0 10.0 11.0 12.0; 13.0 14.0 15.0 16.0]

truncated_svd(A, 2)