text:
using Bridge, StaticArrays, Distributions
using Test, Statistics, Random, LinearAlgebra
using Bridge.Models
using DelimitedFiles
using RCall
using DataFrames

T = 2.0
dt = 1/200
tt = 0.:dt:T

# specify target process
struct IntegratedDiffusion <: ContinuousTimeProcess{ℝ{2}}
    γ::Float64
end

Bridge.b(t, x, P::IntegratedDiffusion) = ℝ{2}(x[2], -(x[2]+sin(x[2]*x[1]) +1/2))
Bridge.σ(t, x, P::IntegratedDiffusion) = ℝ{2}(0.0, P.γ*(1+0.5*atan(x[1])))
Bridge.constdiff(::IntegratedDiffusion) = false

# Generate Data
Random.seed!(1)
P = IntegratedDiffusion(0.7)
W = sample(tt, Wiener())
x0 = ℝ{2}(2.0, 0.0)
X = solve(Euler(), x0, W, P)

# Specify observation scheme
L = @SMatrix [1. 0.]
Σ = @SMatrix [0.0]
v = ℝ{1}(16.0)
#v = ℝ{1}(X.yy[end][1] + √Σ[1,1]*randn())
#v = ℝ{2}(X.yy[end])

# L = @SMatrix [1. 0.; 0. 1.]
# Σ = @SMatrix [0.0 0.0 ; 0.0 0.0]
# v = ℝ{2}(-5.0, 2)    #ℝ{1}(X.yy[end][1] + √Σ[1,1]*randn())

# specify auxiliary process
struct IntegratedDiffusionAux <: ContinuousTimeProcess{ℝ{2}}
    γ::Float64
end

Bridge.B(t, P::IntegratedDiffusionAux) = @SMatrix [0.0 1.0; 0.0 -1.0]
Bridge.β(t, P::IntegratedDiffusionAux) = ℝ{2}(0, 1/2)
Bridge.σ(t, x, P::IntegratedDiffusionAux) = ℝ{2}(0.0, P.γ)
Bridge.constdiff(::IntegratedDiffusionAux) = false

Bridge.b(t, x, P::IntegratedDiffusionAux) = Bridge.B(t,P) * x + Bridge.β(t,P)
Bridge.a(t, P::IntegratedDiffusionAux) = Bridge.σ(t,0,P) * Bridge.σ(t, 0, P)'

#Pt = IntegratedDiffusionAux(0.7)
Pt = IntegratedDiffusionAux(Bridge.σ(T, [v 0], P)[2])

# Solve Backward Recursion
Po = Bridge.PartialBridge(tt, P, Pt, L, v, Σ)

# initalisation
W = sample(tt, Wiener())
x0 = ℝ{2}(2.0, 1.0)
Xo = copy(X)
bridge!(Xo, x0, W, Po)
sample!(W, Wiener())
bridge!(X, x0, W, Po)
ll = llikelihood(Bridge.LeftRule(), X, Po)

# settings sampler
iterations = 50000
subsamples = 0:100:iterations
ρ = 0.9

# further initialisation
Wo = copy(W)
W2 = copy(W)
XX = Any[]
if 0 in subsamples
    push!(XX, copy(X))
end

acc = 0

for iter in 1:iterations
    # Proposal
    sample!(W2, Wiener())
    Wo.yy .= ρ*W.yy + sqrt(1-ρ^2)*W2.yy
    bridge!(Xo, x0, Wo, Po)

    llo = llikelihood(Bridge.LeftRule(), Xo, Po)
    print("ll $ll $llo ")#, X[10], " ", Xo[10])
    if log(rand()) <= llo - ll
        X.yy .= Xo.yy
        W.yy .= Wo.yy
        ll = llo
        print("✓")
        acc +=1
    end
    println()
    if iter in subsamples
        push!(XX, copy(X))
    end
end

# write mcmc iterates to csv file
f = open("output/integrated_diff/iterates.csv","w")
head = "iteration, time, component, value \n"
write(f, head)
iterates = [Any[s, tt[j], d, XX[i].yy[j][d]] for d in 1:2, j in 1:length(X), (i,s) in enumerate(subsamples) ][:]
writecsv(f,iterates)
close(f)

println("Average acceptance percentage: ",100*round(acc/iterations,2),"\n")

# Plotting with ggplot
dd = DataFrame(iteration=map(x->x[1],iterates), time=map(x->x[2],iterates),
        component=map(x->x[3],iterates),value=map(x->x[4],iterates))

ind =  map(x-> x in [0:200:1000;  40000:2000:50000],dd[:iteration])
ddsub = dd[ind,:]

ggplot(ddsub,aes(x=:time,y=:value,group=:iteration,colour=:iteration)) +
        geom_path()+facet_wrap(R"~component",ncol=1,scales="free_y")+
         scale_colour_gradient(low="green",high="blue")

# Calling R, somehow the pdf device does not work
@rlibrary ggplot2
@rlibrary tidyverse
@rlibrary grDevices

R"""
library(ggplot2)
library(tidyverse)
library(grDevices)

theme_set(theme_minimal())
setwd("output/integrated_diff")

d <- read.csv("output/integrated_diff/iterates.csv")
dsub <- d[d$iteration %in% c(seq(0,1000,by=200),seq(40000,50000,by=2000)),]

pdf("iterates.pdf",width=8,height=7)
dsub %>% ggplot(aes(x=time,y=value,colour=iteration)) +
  geom_path(aes(group=iteration)) +
  facet_wrap(~component,ncol=1,scales='free_y')+
  scale_colour_gradient(low='yellow',high='Darkred')
dev.off()
"""

# import variables in Julia
@rget d


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   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.

"""
AbstractBlade.jl defines the AbstractBlade type, interface, and core methods
"""

# --- Exports

# Types
export AbstractBlade

# Properties
import Base.sign
export basis, grade, volume

# Functions
import Base.inv

# --- Types

"""
    AbstractBlade{<:AbstractFloat}

Supertype for all blade types.

Interface
=========

Properties
----------

    grade(B::AbstractBlade)::Int

    basis(B::AbstractBlade)::Matrix

    volume(B::AbstractBlade)::AbstractFloat

    sign(B::AbstractBlade)::Int8

Operations
----------

    dual(B::AbstractBlade, C::AbstractBlade)::AbstractBlade

    inv(B::AbstractBlade)::AbstractBlade

Implementation
==============

* For the AbstractBlade type, the `volume` of a blade should encode both norm
  and orientation information.

* The return value of all methods should preserve the precision of its
  AbstractBlade arguments (when possible).

"""
abstract type AbstractBlade{T<:AbstractFloat} <: AbstractMultivector{T} end

# --- Public functions/methods

"""
    grade(B::AbstractBlade)::Int

Return grade of `B`.
"""
function grade end

"""
    basis(B::AbstractBlade)::Matrix

Return an orthonormal basis for the subspace represented by `B`.

When `B` is an AbstractScalar, 1 (at the precision of `B`) is returned. When
`B` is a Pseudoscalar, LinearAlgebra.I is returned.
"""
function basis end

"""
    volume(B::AbstractBlade)::AbstractFloat

Return the signed volume of `B`.

When `B` is a Blade, `volume(B)` is the signed norm of the blade _relative_ to
its unit basis.
"""
function volume end

"""
    sign(B::AbstractBlade)::Int8

Return the sign of `B` relative to its unit basis.
"""
sign(B::AbstractBlade)::Int8 = sign(volume(B))

"""
    dual(B::AbstractBlade, C::AbstractBlade)::AbstractBlade

Compute the dual `B` relative to the subspace represented by `C`.

Notes
=====

* `dual(B, C)` is only defined if (1) `B` and `C` are extended from real
  vector spaces of the same dimension and (2) the subspace represented by `B`
  is contained in subspace represented by `C`.

* The volume of `C` is ignored.
"""
dual(B::AbstractBlade, C::AbstractBlade)::AbstractBlade =
    nothing  # no-op method to attach docstring to because an empty generic function for
             # `dual` already exists in AbstractMultivector.jl

"""
    inv(B::AbstractBlade)::AbstractBlade

Compute the multiplicative inverse of blade `B`.
"""
function inv end

# --- Method definitions for AbstractMultivector interface functions

grades(B::AbstractBlade) = [grade(B)]

blades(B::AbstractBlade) = Vector{AbstractBlade}([B])

norm(B::AbstractBlade) = abs(volume(B))

Base.getindex(B::AbstractBlade, k::Int) =
    k == grade(B) ? Vector{AbstractBlade}([B]) : Vector{AbstractBlade}()

ext: jl


text:
@testset "DynamicAxis" begin
    let r = DynamicAxis(-5)
        @test isempty(r)
        @test length(r) == 0
        @test size(r) == (0,)
    end
    let r = DynamicAxis(3)
        @test !isempty(r)
        @test last(r) == r.stop
        @test getindex(r, DynamicAxis(2)) == DynamicAxis(2)

        @test getindex(r, OneTo(2)) === OneTo(2)

        @test getindex(r, DynamicAxis(2)) == DynamicAxis(2)
        @test_throws ErrorException r.notfield
        @test length(r) == 3
        @test size(r) == (3,)
        @test first(r) == 1
        @test last(r) == 3
        @test minimum(r) == 1
        @test maximum(r) == 3
        @test r[2] == 2
        @test r[2:3] == 2:3
        @test_throws BoundsError r[4]
        @test_throws BoundsError r[0]
        @test broadcast(+, r, 1) == 2:4
        @test 2*r == mrange(2, step=2, stop=6)
        @test r + r == mrange(2, step=2, stop=6)
        k = 0
        for i in r
            @test i == (k += 1)
        end
        @test intersect(r, DynamicAxis(2)) == DynamicAxis(2)
        @test intersect(r, OneTo(2)) == DynamicAxis(2)
        @test intersect(OneTo(2), r) == DynamicAxis(2)
        @test intersect(r, 0:5) == 1:3
        @test intersect(r, 2) == intersect(2, r) == mrange(2, 2)
        @test findall(in(r), r) == findall(in(mrange(1, length(r))), r) ==
            findall(in(r), mrange(1, length(r))) == mrange(1, length(r))
        @test in(1, r) == true
        @test in(1.0, r) == true
    end

    @test DynamicAxis(DynamicAxis(10)) == DynamicAxis(10)
    @test DynamicAxis(DynamicAxis(10)) == DynamicAxis(UInt(10))
    @test intersect(DynamicAxis(2), OneTo(2)) == OneTo(2)
    @test intersect(OneTo(2), DynamicAxis(2)) == DynamicAxis(2)

    @test issubset(DynamicAxis(2), OneTo(2)) == true
    @test issubset(OneTo(2), DynamicAxis(2)) == true

    let r = DynamicAxis(7)
        @test findall(in(mrange(2, (length(r) - 1))), r) == mrange(2, (length(r) - 1))
        @test findall(in(r), mrange(2, (length(r) - 1))) == mrange(1, (length(r) - 2))
    end
    @test convert(DynamicAxis, 1:2) == DynamicAxis(2)
    @test_throws ArgumentError("first element must be 1, got 2") convert(DynamicAxis, 2:3)
    @test_throws ArgumentError("step must be 1, got 2") convert(DynamicAxis, 1:2:5)
    @test DynamicAxis(1:2) == DynamicAxis(2)
    @test DynamicAxis(1:1:2) == DynamicAxis(2)
    @test DynamicAxis(1:2) == DynamicAxis(2)
    @test DynamicAxis(Int32(1):Int32(2)) == DynamicAxis(2)
    @test DynamicAxis(3.0) == DynamicAxis(3)
    @test_throws InexactError DynamicAxis(3.2)

    x = DynamicAxis(10)
    StaticRanges.grow_end!(x, 1)
    @test length(x) == 11
    StaticRanges.shrink_end!(x, 2)
    @test length(x) == 9
end

if VERSION > v"1.2"
    @test mod(22, DynamicAxis(10)) == mod(22, Base.OneTo(10))
end

r = DynamicAxis(10)
empty!(r)
@test last(r) == 0

@test AbstractUnitRange{Int}(DynamicAxis(10)) isa DynamicAxis

@test_throws ErrorException r.throw_error_please
@testset "show" begin
    r = DynamicAxis(10)
    io = IOBuffer()
    show(io, r)
    str = String(take!(io))
    @test str == "DynamicAxis(10)"
end


text:
"""
StrateGrid(x::NamedTuple)

Let x = (k1 = v1, k2 = v2, ...., kN = vN) be a NamedTuple of AbstractVectors for state space,
StateGrid(x) returns an AbstractArray M such that
M[i1, i2..., iN] = (k1 = v1[i1], k2 = v2[i2], ...., kN = vN[iN]
"""
struct StateGrid{T, N, C <: NamedTuple} <: AbstractArray{T, N}
    x::C
end
function StateGrid(x::NamedTuple{Names, <: NTuple{N, <: AbstractVector{T}}}) where {Names, N, T}
    StateGrid{T, N, typeof(x)}(x)
end
function Base.eltype(stategrid::StateGrid{T, N, <: NamedTuple{Names, V}}) where {T, N, Names, V}
    NamedTuple{Names, NTuple{N, T}}
end
Base.ndims(stategrid::StateGrid{T, N}) where {T, N} = N
Base.size(stategrid::StateGrid{T, N}) where {T, N} = ntuple(i -> length(stategrid.x[i]), N)
Base.eachindex(stategrid::StateGrid) = CartesianIndices(size(stategrid))
function Base.getindex(stategrid::StateGrid, args::CartesianIndex)
    eltype(stategrid)(ntuple(i -> stategrid.x[i][args[i]], ndims(stategrid)))
end
Base.getindex(stategrid::StateGrid, x::Symbol) = stategrid.x[x]

#========================================================================================

Derive

========================================================================================#
# 1 state variable
@generated function differentiate(::Type{Tsolution}, grid::StateGrid{T1, 1, <: NamedTuple{N}}, y::AbstractArray{T}, icar, bc) where {Tsolution, T1, N, T}
    statename = N[1]
    expr = Expr[]
    for k in 1:length(Tsolution.parameters[1])
        solname = Tsolution.parameters[1][k]
        push!(expr, Expr(:(=), solname, :(y[i, $k])))
        push!(expr, Expr(:(=), Symbol(solname, statename, :_up), :((i < size(y, 1)) ? (y[i+1, $k] - y[i, $k]) / Δxp : convert($T, bc[end, $k]))))
        push!(expr, Expr(:(=), Symbol(solname, statename, :_down), :((i > 1) ? (y[i, $k] - y[i-1, $k]) / Δxm : convert($T, bc[1, $k]))))
        push!(expr, Expr(:(=), Symbol(solname, statename, statename), :((1 < i < size(y, 1)) ? (y[i + 1, $k] / (Δxp * Δx) + y[i - 1, $k] / (Δxm * Δx) - 2 * y[i, $k] / (Δxp * Δxm)) : ((i == 1) ? (y[2, $k] / (Δxp * Δx) + (y[1, $k] - bc[1, $k] * Δxm) / (Δxm * Δx) - 2 * y[1, $k] / (Δxp * Δxm)) : ((y[end, $k] + bc[end, $k] * Δxp) / (Δxp * Δx) + y[end - 1, $k] / (Δxm * Δx) - 2 * y[end, $k] / (Δxp * Δxm))))))
    end
    quote
        $(Expr(:meta, :inline))
        i = icar[1]
        grida = grid.x[1]
        @inbounds Δxm = grida[max(i, 2)] - grida[max(i-1, 1)]
        @inbounds Δxp = grida[min(i+1, size(y, 1))] - grida[min(i, size(y, 1) - 1)]
        Δx = (Δxm + Δxp) / 2
        @inbounds $(Expr(:tuple, expr...))
    end
end

# 2 state variables
@generated function differentiate(::Type{Tsolution}, grid::StateGrid{T1, 2, <: NamedTuple{N}}, y::AbstractArray{T}, icar, bc) where {Tsolution, T1, N, T}
    statename1 = N[1]
    statename2 = N[2]
    expr = Expr[]
    for k in 1:length(Tsolution.parameters[1])
        solname = Tsolution.parameters[1][k]
        push!(expr, Expr(:(=), solname, :(y[i1, i2, $k])))
        push!(expr, Expr(:(=), Symbol(solname, statename1, :_up), :((i1 < size(y, 1)) ? (y[i1+1, i2, $k] - y[i1, i2, $k]) / Δx1p : convert($T, bc[end, i2, $k]))))
        push!(expr, Expr(:(=), Symbol(solname, statename1, :_down), :((i1 > 1) ? (y[i1, i2, $k] - y[i1-1, i2, $k]) / Δx1m : convert($T, bc[1, i2, $k]))))
        push!(expr, Expr(:(=), Symbol(solname, statename2, :_up), :((i2 < size(y, 2)) ? (y[i1, i2+1, $k] - y[i1, i2, $k]) / Δx2p : convert($T, bc[i1, end, $k]))))
        push!(expr, Expr(:(=), Symbol(solname, statename2, :_down),  :((i2 > 1) ? (y[i1, i2, $k] - y[i1, i2-1, $k]) / Δx2m : convert($T, bc[i1, 1, $k]))))
        push!(expr, Expr(:(=), Symbol(solname, statename1, statename1), :((1 < i1 < size(y, 1)) ? (y[i1 + 1, i2, $k] / (Δx1p * Δx1) + y[i1 - 1, i2, $k] / (Δx1m * Δx1) - 2 * y[i1, i2, $k] / (Δx1p * Δx1m)) : ((i1 == 1) ? (y[2, i2, $k] / (Δx1p * Δx1) + (y[1, i2, $k] - bc[1, i2, $k] * Δx1m) / (Δx1m * Δx1) - 2 * y[1, i2, $k] / (Δx1p * Δx1m)) : ((y[end, i2, $k] + bc[end, i2, $k] * Δx1p) / (Δx1p * Δx1) + y[end - 1, i2, $k] / (Δx1m * Δx1) - 2 * y[end, i2, $k] / (Δx1p * Δx1m))))))
        push!(expr, Expr(:(=), Symbol(solname, statename2, statename2), :((1 < i2 < size(y, 2)) ? (y[i1, i2 + 1, $k] / (Δx2p * Δx2) + y[i1, i2 - 1, $k] / (Δx2m * Δx2) - 2 * y[i1, i2, $k] / (Δx2p * Δx2m)) : ((i2 == 1) ? (y[i1, 2, $k] / (Δx2p * Δx2) + (y[i1, 1, $k] - bc[i1, 1, $k] * Δx2m) / (Δx2m * Δx2) - 2 * y[i1, 1, $k] / (Δx2p * Δx2m)) : ((y[i1, end, $k] + bc[i1, end, $k] * Δx2p) / (Δx2p * Δx2) + y[i1, end - 1, $k] / (Δx2m * Δx2) - 2 * y[i1, end, $k] / (Δx2p * Δx2m))))))
        push!(expr, Expr(:(=), Symbol(solname, statename1, statename2), :((y[min(i1 + 1, size(y, 1)), min(i2 + 1, size(y, 2)), $k] - y[min(i1 + 1, size(y, 1)), max(i2 - 1, 1), $k] - y[max(i1 - 1, 1), min(i2 + 1, size(y, 2)), $k] + y[max(i1 - 1, 1), max(i2 - 1, 1), $k]) / (4 * Δx1 * Δx2))))
    end
    quote
        $(Expr(:meta, :inline))
        i1, i2 = icar[1], icar[2]
        grid1, grid2 = grid.x[1], grid.x[2]
        @inbounds Δx1m = grid1[max(i1, 2)] - grid1[max(i1-1, 1)]
        @inbounds Δx1p = grid1[min(i1+1, size(y, 1))] - grid1[min(i1, size(y, 1) - 1)]
        Δx1 = (Δx1m + Δx1p) / 2
        @inbounds Δx2m = grid2[max(i2, 2)] - grid2[max(i2-1, 1)]
        @inbounds Δx2p = grid2[min(i2+1, size(y, 2))] - grid2[min(i2, size(y, 2) - 1)]
        Δx2 = (Δx2m + Δx2p) / 2
        @inbounds $(Expr(:tuple, expr...))
    end
end


text:
module Contractors
using ValidatedNumerics
using DualNumbers

export N_rig, root, range_estimate, ShootingMethod

function N_rig(f, f′, x::Interval{T}) where {T}
	x_mid = Interval{T}(mid(x))
	return intersect(x, x_mid-f(x_mid)/f′(x))
end

root(f, x::Interval, ϵ; max_iter = 100) = root(f, x->f(Dual(x, 1)).epsilon, x, ϵ; max_iter = max_iter)

function root(f, f′, x::Interval, ϵ; max_iter = 100)
	for i in 1:max_iter
		x = N_rig(f, f′, x)
		if diam(x)<ϵ
			return x
		end
		if isempty(x)
			return x
		end
	end
	@info "Maximum iterates reached" max_iter
	return x
end

# superseded by IntervalOptimisation.jl
function range_estimate(f, domain, recstep = 5)
	if recstep == 1
		return f(domain)
	else
		a, b = bisect(domain)
		Iₐ = range_estimate(f, a, recstep-1)
		Iᵦ = range_estimate(f, b, recstep-1)
		return union(Iₐ, Iᵦ)
	end
end

# a shooting method to find the k-th preimage of a point

#function that gives me the type of the endpoints of the interval

using LinearAlgebra

# this function generates the Jacobian for
# f(x_0)=x_1, f(x_1)=x_2, ..., f(x_{n-1})=y

coeff_interval(x::Array{Interval{T}, 1}) where {T} = T

function Jac(fprime, v::Vector{T}) where {T}
    dv = fprime.(v)
    ev = -ones(T, length(v)-1)
    return Bidiagonal{T}(dv, ev, :U)
end

function ShootingMethod(f, fprime, n, x, y, rigstep = 10)
	F = x->(f.(x)-[x[2:end]; y])

	for i in 1:rigstep
		x_mid = Interval{coeff_interval(x)}.(mid.(x))
		x = intersect.(x, x_mid-Jac(fprime, x)\F(x_mid))
	end
	return x
end

end

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
using PyPlot
using HDF5
using Printf
import Statistics

function load_recall_curves(ntrials, toplot, dataset_name, RPATH, fname)
  recalls = zeros(ntrials, 1000)
  for trial = 1:ntrials
    recalls[trial, :] = h5read(joinpath(RPATH, dataset_name, fname), "$trial/recall")
  end
  meanrecall, stdrecall = Statistics.mean(recalls, dims=1)[toplot], Statistics.std(recalls, dims=1)[toplot]
end

function load_recall_curves(ntrials, toplot, dataset_name, RPATH, fname, nilsiters)
  recalls = zeros(ntrials, 1000)
  for trial = 1:ntrials
    # @show RPATH, dataset_name, fname
    recalls[trial, :] = h5read(joinpath(RPATH, dataset_name, fname), "$trial/recall_$nilsiters")
  end
  meanrecall, stdrecall = Statistics.mean(recalls, dims=1)[toplot], Statistics.std(recalls, dims=1)[toplot]
end

function get_title(dataset_name)
  if dataset_name == "sift1m"
    return "SIFT1M"
  elseif dataset_name == "convnet1m"
    return "Convnet1M"
  elseif dataset_name == "mnist"
    return "MNIST"
  elseif dataset_name == "labelme"
    return "LabelMe"
  end
end

function print_recalls(meanrecall, stdrecall)
  for j = [1, 2, 5, 10, 20, 50, 100]
    @printf("r@%d %.2f \\pm %.2f\n", j, 100 * meanrecall[j], 100 * stdrecall[j]);
  end
  println()
end

function make_plots(RPATH, dnames, sfile)
  ms      = [7]
  ilsit   = [32]

  sz = 13
  linew = 2
  toplot = 1:1000
  iterations = 25

  # Set the overall size here
  fig = figure("recall_plot", figsize=( (5+0.6)*length(dnames), 5))
  ntrials = 10

  subplotidx = 1;
  for dataset_name = dnames
    ax = subplot(1,length(dnames),subplotidx)
    for midx = 1:length( ms )

      ax[:set_prop_cycle](nothing)  # Reset colour cycle
      m = ms[ midx ]

      # Load srd
      for i=length(ilsit):-1:1
        meanrecall, stdrecall = load_recall_curves(ntrials, toplot, dataset_name, RPATH, "srd_m$(m)_it$(iterations).h5")
        println("$dataset_name srd m=$m"); print_recalls(meanrecall, stdrecall)
        plot(toplot, meanrecall[toplot], label="SRD-$(ilsit[i]) $((m+1)*8) bits", m == 7 ? "--" : "", lw=linew)
      end

      # Load lsq
      for i=length(ilsit):-1:1
        meanrecall, stdrecall = load_recall_curves(ntrials, toplot, dataset_name, RPATH, "lsq_m$(m)_it$(iterations).h5")
        println("$dataset_name lsq m=$m"); print_recalls(meanrecall, stdrecall)
        plot(toplot, meanrecall[toplot], label="LSQ-$(ilsit[i]) $((m+1)*8) bits", m == 7 ? "--" : "", lw=linew)
      end

      # if dataset_name == "convnet1m"
      #   # Skip the color for consistency
      #   plot([], [])
      # else
      #   bpath = joinpath("/home/julieta/Desktop/CQ/build/results/", lowercase(dataset_name) * "_$(m+1)", "trial_3")
      #   recall = h5read(joinpath(bpath, "recall.h5"), "recall")
      #   println("$dataset_name cq m=$(m+1)"); print_recalls(recall, zeros(size(recall)))
      #   plot(toplot, recall[toplot], label="CQ $((m+1)*8) bits", m == 7 ? "--" : "", lw=linew)
      # end

      # Load rvq and ervq
      for baseline = ["ervq", "rvq"]
        meanrecall, stdrecall = load_recall_curves(ntrials, toplot, dataset_name, RPATH, "$(baseline)_m$(m)_it$(iterations).h5")
        println("$dataset_name $baseline m=$m"); print_recalls(meanrecall, stdrecall)
        plot(toplot, meanrecall[toplot], label="$(uppercase(baseline)) $((m+1)*8) bits", m == 7 ? "--" : "", lw=linew)
      end

      # Load pq and opq
      for baseline = ["opq", "pq"]
        meanrecall, stdrecall = load_recall_curves(ntrials, toplot, dataset_name, RPATH, "$(baseline)_m$(m+1)_it$(iterations).h5")
        println("$dataset_name $baseline m=$m"); print_recalls(meanrecall, stdrecall)
        plot(toplot, meanrecall[toplot], label="$(uppercase(baseline)) $((m+1)*8) bits", m == 7 ? "--" : "", lw=linew)
      end

      grid(true)
      title(get_title(dataset_name), size=sz)
      ax[:set_xscale]("log")

      ymin = dataset_name == "sift1m" ? 0.2 : 0
      ymin = dataset_name == "mnist" ? 0.3 : ymin
      ymin = dataset_name == "labelme" ? 0.2 : ymin
      ax[:set_ylim]([ymin, 1])
      ax[:set_xlim]([1, dataset_name == "sift1m" || dataset_name == "convnet1m" ? 1000 : 100])

      ylabel("Recall@N", size=sz)
      xlabel("N", size=sz)
      legend(loc="lower right", fontsize=sz-2, ncol=2)

    end
    subplotidx += 1
  end

  PyPlot.tight_layout()
  savefig( sfile )
end

RPATH = "./results/"
dataset = "sift1m"
savename = "./results/sift_plot.pdf"
make_plots(RPATH, [dataset], savename)


text:
using Flux

"""
Get `m` random samples from a dataset.
"""
function Base.rand(data::Flux.Data.DataLoader, m::Int=1)
    idx = rand(data.indices, m)

    x = selectrand(data.data[1], idx)
    y = selectrand(data.data[2], idx)
    return (x,y)
end

selectrand(data::Union{Matrix, Flux.OneHotMatrix}, idx) = data[:, idx]
selectrand(data::Union{BitArray, Vector}, idx) = data[idx]

function getmapping(data)
    mapping = Dict()
    for (i, (x,y)) in enumerate(data)
        mapping[i] = (x, y)
    end
    return mapping
end


text:
function assess_convergence(x::Array,
                            x_previous::Array,
                            f_x::Real,
                            f_x_previous::Real,
                            g::Array,
                            x_tol::Real,
                            f_tol::Real,
                            g_tol::Real)
    x_converged, f_converged, g_converged = false, false, false

    if maxdiff(x, x_previous) < x_tol
        x_converged = true
    end

    # Absolute Tolerance
    # if abs(f_x - f_x_previous) < f_tol
    # Relative Tolerance
    if abs(f_x - f_x_previous) / (abs(f_x) + f_tol) < f_tol || nextfloat(f_x) >= f_x_previous
        f_converged = true
    end

    if vecnorm(g, Inf) < g_tol
        g_converged = true
    end

    converged = x_converged || f_converged || g_converged

    return x_converged, f_converged, g_converged, converged
end

function assess_convergence(state, options)
    x_converged, f_converged, g_converged = false, false, false

    if maxdiff(state.x, state.x_previous) < options.x_tol
        x_converged = true
    end

    # Absolute Tolerance
    # if abs(f_x - f_x_previous) < f_tol
    # Relative Tolerance
    if abs(state.f_x - state.f_x_previous) / (abs(state.f_x) + options.f_tol) < options.f_tol || nextfloat(state.f_x) >= state.f_x_previous
        f_converged = true
    end

    if vecnorm(state.g, Inf) < options.g_tol
        g_converged = true
    end

    converged = x_converged || f_converged || g_converged

    return x_converged, f_converged, g_converged, converged
end

function assess_convergence(state::NelderMeadState, options)
    g_converged = state.nm_x <= options.g_tol # Hijact g_converged for NM stopping criterior
    return false, false, g_converged, g_converged
end

function assess_convergence(state::Union{ParticleSwarmState, SimulatedAnnealingState}, options)
    false, false, false, false
end

function assess_convergence(state::NewtonTrustRegionState, options)
    x_converged, f_converged, g_converged, converged = false, false, false, false
    if state.rho > state.eta
        # Accept the point and check convergence
        x_converged,
        f_converged,
        g_converged,
        converged = assess_convergence(state.x,
                                       state.x_previous,
                                       state.f_x,
                                       state.f_x_previous,
                                       state.g,
                                       options.x_tol,
                                       options.f_tol,
                                       options.g_tol)
    end
    x_converged, f_converged, g_converged, converged
end


text:
@testset "Test Circle" begin
    M = Circle()
    x = S1Point(0.)
    y = opposite(M,x)
    z = randomMPoint(M)
    η = randomTVector(M,x)
    @test getValue(x) == 0.
    ξ = S1TVector(2*π/3)
    @test getValue(ξ) == 2*π/3
    @test distance(M, x, exp(M,x,ξ)) ≈ 2*π/3
    @test dot(M,x,ξ, S1TVector(0.5)) ≈ π/3
    @test embed(M,x) == SnPoint([1., 0.])
    @test exp(M,x,ξ) == S1Point(2*π/3)
    @test distance(M, exp(M,x,2*ξ), S1Point(-2*π/3) ) ≈ 0 atol=10^(-15)
    @test log(M,x,exp(M,x,ξ)) == ξ
    @test manifoldDimension(x) == 1
    @test manifoldDimension(M) == 1
    @test norm(M,x,ξ) ≈ 2*π/3
    @test distance(M, opposite(M,S1Point(π/2)), S1Point(-π/2) ) ≈ 0 atol=10^(-15)
    @test parallelTransport(M,x,S1Point(π/2),ξ) == ξ
    @test validateMPoint(M,randomMPoint(M))
    @test validateTVector(M,x,randomTVector(M,x))
    @test typeofTVector(typeof(x)) == S1TVector
    @test typeofMPoint(typeof(ξ)) == S1Point
    @test typeofTVector(x) == S1TVector
    @test typeofMPoint(ξ) == S1Point
    @test typicalDistance(M) == π/2
    @test zeroTVector(M,x) == S1TVector(0.)
    @test "$M" == "The manifold S1 consisting of angles"
    @test "$x" == "S1(0.0)"
    @test "$(S1TVector(0.))" == "S1T(0.0)"
    @test_throws ErrorException validateMPoint(M,S1Point(2*π))
    @test_throws ErrorException validateMPoint(M,S1Point(-2*π))

    # Test manifoldDimension
    @test manifoldDimension(M) == 1
    @test manifoldDimension(x) == 1
    # Test norm
    @test norm(M,x,ξ) == 2*π/3
    # Test opposite
    @test getValue(y) - pi ≈ 0 atol = 10.0^(-16)
    # Test parallelTransport
    @test getValue(parallelTransport(M,x,y,ξ)) - getValue(ξ) ≈ 0 atol = 10.0^(-16)
    # Test typicalDistance
    @test typicalDistance(M) == pi/2
    # Test validateMPoint and validateTVector
    xnot = S1Point(3.15)
    @test_throws ErrorException validateMPoint(M,xnot)
end


text:
function test_boundary_function(B, X1, X2, func)
    boundary_function = BoundaryFunction{B, X1, X2}(func)
    return true
end

function test_z_boundary_condition_simple(arch, FT, fldname, bctype, bc, Nx, Ny)
    Nz = 16
    bc = BoundaryCondition(bctype, bc)
    fieldbcs = HorizontallyPeriodicBCs(top=bc)
    modelbcs = BoundaryConditions(; Dict(fldname=>fieldbcs)...)

    model = Model(grid=RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(0.1, 0.2, 0.3)), architecture=arch,
                  float_type=FT, boundary_conditions=modelbcs)

    time_step!(model, 1, 1e-16)

    return typeof(model) <: Model
end

function test_z_boundary_condition_top_bottom_alias(arch, FT, fldname)
    N, val = 16, 1.0
    top_bc = BoundaryCondition(Value, val)
    bottom_bc = BoundaryCondition(Value, -val)
    fieldbcs = HorizontallyPeriodicBCs(top=top_bc, bottom=bottom_bc)
    modelbcs = BoundaryConditions(; Dict(fldname=>fieldbcs)...)

    model = Model(grid=RegularCartesianGrid(FT; size=(N, N, N), length=(0.1, 0.2, 0.3)), architecture=arch,
                  float_type=FT, boundary_conditions=modelbcs)

    bcs = getfield(model.boundary_conditions.solution, fldname)

    time_step!(model, 1, 1e-16)

    return getbc(bcs.z.top) == val && getbc(bcs.z.bottom) == -val
end

function test_z_boundary_condition_array(arch, FT, fldname)
    Nx = Ny = Nz = 16

    bcarray = rand(FT, Nx, Ny)

    if arch == GPU()
        bcarray = CuArray(bcarray)
    end

    value_bc = BoundaryCondition(Value, bcarray)
    fieldbcs = HorizontallyPeriodicBCs(top=value_bc)
    modelbcs = BoundaryConditions(; Dict(fldname=>fieldbcs)...)

    model = Model(grid=RegularCartesianGrid(FT; size=(Nx, Ny, Nz), length=(0.1, 0.2, 0.3)), architecture=arch,
                  float_type=FT, boundary_conditions=modelbcs)

    bcs = getfield(model.boundary_conditions.solution, fldname)

    time_step!(model, 1, 1e-16)

    return bcs.z.top[1, 2] == bcarray[1, 2]
end

function test_flux_budget(arch, FT, fldname)
    N, κ, Lz = 16, 1, 0.7

    bottom_flux = FT(0.3)
    flux_bc = BoundaryCondition(Flux, bottom_flux)
    fieldbcs = HorizontallyPeriodicBCs(bottom=flux_bc)
    modelbcs = BoundaryConditions(; Dict(fldname=>fieldbcs)...)

    grid = RegularCartesianGrid(FT; size=(N, N, N), length=(1, 1, Lz))
    closure = ConstantIsotropicDiffusivity(FT; ν=κ, κ=κ)
    model = Model(grid=grid, closure=closure, architecture=arch,
                  float_type=FT, buoyancy=nothing, boundary_conditions=modelbcs)

    if fldname ∈ (:u, :v, :w)
        field = getfield(model.velocities, fldname)
    else
        field = getfield(model.tracers, fldname)
    end

    @. field.data = 0

    bcs = getfield(model.boundary_conditions.solution, fldname)
    mean_init = mean(interior(field))

    τκ = Lz^2 / κ   # Diffusion time-scale
    Δt = 1e-6 * τκ  # Time step much less than diffusion time-scale
    Nt = 100        # Number of time steps

    time_step!(model, Nt, Δt)

    # budget: Lz*∂<ϕ>/∂t = -Δflux = -top_flux/Lz (left) + bottom_flux/Lz (right)
    # therefore <ϕ> = bottom_flux * t / Lz
    return isapprox(mean(interior(field)) - mean_init, bottom_flux * model.clock.time / Lz)
end

@testset "Boundary conditions" begin
    println("Testing boundary conditions...")

    @testset "Boundary functions" begin
        simple_bc(ξ, η, t) = exp(ξ) * cos(η) * sin(t)
        for B in (:x, :y, :z)
            for X1 in (:Face, :Cell)
                @test test_boundary_function(B, X1, Cell, simple_bc)
            end
        end
    end

    funbc(args...) = π
    boundaryfunbc = BoundaryFunction{:z, Face, Cell}((ξ, η, t) -> exp(ξ) * cos(η) * sin(t))

    @testset "Boundary condition instatiation and time-stepping" begin
        Nx = Ny = 16
        for arch in archs
            for FT in float_types
                for fld in (:u, :v, :T, :S)
                    for bctype in (Gradient, Flux, Value)

                        arraybc = rand(FT, Nx, Ny)
                        if arch == GPU()
                            arraybc = CuArray(arraybc)
                        end

                        for bc in (FT(0.6), arraybc, funbc, boundaryfunbc)
                            @test test_z_boundary_condition_simple(arch, FT, fld, bctype, bc, Nx, Ny)
                        end
                    end

                    @test test_z_boundary_condition_top_bottom_alias(arch, FT, fld)
                    @test test_z_boundary_condition_array(arch, FT, fld)
                    @test test_flux_budget(arch, FT, fld)
                end
            end
        end
    end
end


text:
# ---
# title: 520. Detect Capital
# id: problem520
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/detect-capital/description/>
# hidden: true
# ---
#
# Given a word, you need to judge whether the usage of capitals in it is right
# or not.
#
# We define the usage of capitals in a word to be right when one of the
# following cases holds:
#
#   1. All letters in this word are capitals, like "USA".
#   2. All letters in this word are not capitals, like "leetcode".
#   3. Only the first letter in this word is capital, like "Google".
#
# Otherwise, we define that this word doesn't use capitals in a right way.
#
#
#
# **Example 1:**
#
#
#
#     Input: "USA"
#     Output: True
#
#
#
#
# **Example 2:**
#
#
#
#     Input: "FlaG"
#     Output: False
#
#
#
#
# **Note:** The input will be a non-empty word consisting of uppercase and
# lowercase latin letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test
using EFIT

g = readg(EFIT.test_gfile)

lt, ut = triangularity(g)
@test lt == 0.059614676027860296
@test ut == 0.05822145848512557

@test ellipticity(g) == 1.475784591289634
@test major_radius(g) == 1.648852555
@test minor_radius(g) == 0.633397135
@test aspect_ratio(g) == 2.6031891587258285


text:
isdefined(Base, :__precompile__) && __precompile__()

module CTM

include("generated/includes.jl")

export Bus
export Generator

end # module

max_stars_repo_name: NREL-SIIP/CommonTransmissionModel

max_issues_repo_name: NREL-SIIP/CommonTransmissionModel

max_forks_repo_name: NREL-SIIP/CommonTransmissionModel


text:
# Fast sincos function for Julia 0.5
# (see https://github.com/JuliaLang/julia/pull/21589)
@inline function sincos(v::Float64)
    return Base.llvmcall("""
        %f = bitcast i8 *%1 to void (double, double *, double *)*
        %ps = alloca double
        %pc = alloca double
        call void %f(double %0, double *%ps, double *%pc)
        %s = load double, double* %ps
        %c = load double, double* %pc
        %res0 = insertvalue [2 x double] undef, double %s, 0
        %res = insertvalue [2 x double] %res0, double %c, 1
        ret [2 x double] %res
        """, Tuple{Float64,Float64}, Tuple{Float64,Ptr{Void}}, v,
        cglobal((:sincos, Base.Math.libm)))
end

sincos(v::Real) = ( (r, i) = reim(cis(v)); (i, r) )


text:
module SSLTestMatching
using SimpleSDMLayers
using Test

M = rand(Bool, 6, 9)
for rep in 1:10
    l = rand()
    r = 4l
    b = rand()
    t = 12b
    S = SimpleSDMPredictor(M, l, r, b, t)

    @test SimpleSDMLayers._match_latitude(S, S.bottom) == 1
    @test SimpleSDMLayers._match_latitude(S, S.top) == 6
    @test SimpleSDMLayers._match_longitude(S, S.left) == 1
    @test SimpleSDMLayers._match_longitude(S, S.right) == 9
end

S = SimpleSDMResponse(rand(Bool, 4, 4), 0.3, 0.7, 0.5, 0.9)

@test SimpleSDMLayers._match_latitude(S, 0.2) === nothing
@test SimpleSDMLayers._match_latitude(S, 1.8) === nothing
@test SimpleSDMLayers._match_latitude(S, 0.5) == 1
for (i,l) in enumerate(latitudes(S))
    @test SimpleSDMLayers._match_latitude(S, l) == i
end
for (i,l) in enumerate(longitudes(S))
    @test SimpleSDMLayers._match_longitude(S, l) == i
end

end


text:
"""
    RegressionGSA(rank) <: GSAMethod

RegressionGSA methods for global sensitivity analysis. Providing this to `gsa` results
in a calculation of the following statistics, provided as a `RegressionGSAResult`. If
the function f to be analyzed is of dimensionality f: R^n -> R^m, then these coefficients
are returned as a matrix, with the corresponding statistic in the (i, j) entry.

- `pearson`: This is equivalent to the correlation coefficient matrix between input and
output. The rank version is known as the Spearman coefficient.
- `standard_regression`: Standard regression coefficients, also known as sigma-normalized
derivatives
- `partial_correlation`: Partial correlation coefficients, related to the precision matrix
and a measure of the correlation of linear models of the

# Arguments
- `rank::Bool = false`: Flag determining whether to also run a rank regression analysis
"""
@with_kw mutable struct RegressionGSA <: GSAMethod
    rank::Bool = false
end

struct RegressionGSAResult{T, TR}
    pearson::T
    standard_regression::T
    partial_correlation::T
    pearson_rank::TR
    standard_rank_regression::TR
    partial_rank_correlation::TR
end

function gsa(f, method::RegressionGSA, p_range::AbstractVector;samples::Int = 1000, batch::Bool = false, kwargs...)
    lb = [i[1] for i in p_range]
    ub = [i[2] for i in p_range]
    X = QuasiMonteCarlo.sample(samples, lb, ub, QuasiMonteCarlo.SobolSample())
    desol = false

    if batch
        _y = f(X)
        multioutput = _y isa AbstractMatrix
        Y = multioutput ? _y : reshape(_y, 1, length(_y))
    else
        _y = [f(X[:, j]) for j in axes(X, 2)]
        multioutput = !(eltype(_y) <: Number)
        if eltype(_y) <: RecursiveArrayTools.AbstractVectorOfArray
            y_size = size(_y[1])
            _y = vec.(_y)
            desol = true
        end
        Y = multioutput ? reduce(hcat,_y) : reshape(_y, 1, length(_y))
    end

    srcs = _calculate_standard_regression_coefficients(X, Y)
    corr = _calculate_correlation_matrix(X, Y)
    partials = _calculate_partial_correlation_coefficients(X, Y)

    if method.rank
        X_rank = vcat((sortperm(view(X, i, :))' for i in axes(X, 1))...)
        Y_rank = vcat((sortperm(view(Y, i, :))' for i in axes(Y, 1))...)

        srcs_rank = _calculate_standard_regression_coefficients(X_rank, Y_rank)
        corr_rank = _calculate_correlation_matrix(X_rank, Y_rank)
        partials_rank = _calculate_partial_correlation_coefficients(X_rank, Y_rank)

        return RegressionGSAResult(
            corr,
            srcs,
            partials,
            corr_rank,
            srcs_rank,
            partials_rank
        )
    end

    return RegressionGSAResult(
        corr,
        srcs,
        partials,
        nothing, nothing, nothing
    )
end

function _calculate_standard_regression_coefficients(X, Y)
    β̂ = X' \ Y'
    srcs = β̂ .* std(X, dims = 2) ./ std(Y, dims = 2)'
    return srcs
end

function _calculate_correlation_matrix(X, Y)
    corr = cov(X, Y, dims = 2) ./ (std(X, dims = 2) .* std(Y, dims = 2)')
    return corr
end

function _calculate_partial_correlation_coefficients(X, Y)
    XY = vcat(X, Y)
    corr = cov(XY, dims = 2) ./ (std(XY, dims = 2) .* std(XY, dims = 2)')
    prec = pinv(corr) # precision matrix
    pcc_XY = -prec ./ sqrt.(diag(prec) .* diag(prec)')
    # return partial correlation matrix relating f: X -> Y model values
    return pcc_XY[axes(X, 1), lastindex(X, 1) .+ axes(Y, 1)]
end


text:
using LinearAlgebra
using SparseArrays

using MeshCat
using GeometryBasics
using Colors
using CoordinateTransformations
using FiniteDiff
using Rotations

# visualize
function set_background!(vis::Visualizer; top_color=RGBA(1,1,1.0), bottom_color=RGBA(1,1,1.0))
    setprop!(vis["/Background"], "top_color", top_color)
    setprop!(vis["/Background"], "bottom_color", bottom_color)
end
vis = Visualizer()
open(vis)
# set_background!(vis)

## setup
xa = [0.0; 0.0; 0.0]
qa = [1.0; 0.0; 0.0; 0.0]
ra = 0.25

xb = [1.0; 1.0e-13; 0.0]
qb = [1.0; 0.0; 0.0; 0.0]
rb = 0.25
collision = SphereSphereCollision{Float64,3,3,9}(
        szeros(3),
        szeros(3),
        ra,
        rb)

pa = contact_point(:parent, collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false))
pb = contact_point(:child,  collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false))

n = contact_normal(collision,  xa, Quaternion(qa..., false), xb, Quaternion(qb..., false))
t = contact_tangent(collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false))

# sphere a
sa = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), ra)
color1 = Colors.RGBA(0.7, 0.7, 0.7, 0.5);
setobject!(vis[:spherea], sa, MeshPhongMaterial(color=color1))
settransform!(vis[:spherea], Translation(xa...))

# sphere b
sb = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), rb)
color2 = Colors.RGBA(0.7, 0.7, 0.7, 0.5);
setobject!(vis[:sphereb], sb, MeshPhongMaterial(color=color2))
settransform!(vis[:sphereb], Translation(xb...))

# closest points
color_cp = Colors.RGBA(0.0, 0.0, 0.0, 1.0);
cs1 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), 0.025)
cs2 = GeometryBasics.Sphere(Point(0.0, 0.0, 0.0), 0.025)
setobject!(vis[:cs1], cs1, MeshPhongMaterial(color=color_cp))
setobject!(vis[:cs2], cs2, MeshPhongMaterial(color=color_cp))
settransform!(vis[:cs1], Translation(pa...))
settransform!(vis[:cs2], Translation(pb...))

# normals (NOTE: visualize opposite directions)
set_arrow!(vis, pa, -vec(n), name=:parent_normal)
set_arrow!(vis, pb, vec(n), name=:child_normal)

set_arrow!(vis, pa, vec(t[1, :]), name=:parent_tangent1, color=RGBA(0.0, 1.0, 0.0, 1.0))
set_arrow!(vis, pa, vec(t[2, :]), name=:parent_tangent2, color=RGBA(0.0, 0.0, 1.0, 1.0))

set_arrow!(vis, pb, -1.0 * vec(t[1, :]), name=:child_tangent1, color=RGBA(0.0, 1.0, 0.0, 1.0))
set_arrow!(vis, pb, -1.0 * vec(t[2, :]), name=:child_tangent2, color=RGBA(0.0, 0.0, 1.0, 1.0))

## Jacobians contact_normal(collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false))

# FiniteDiff.finite_difference_jacobian(x -> vec(contact_normal(collision, x, Quaternion(qa..., false), xb, Quaternion(qb..., false))), xa)

λ = rand(3)
λ' * FiniteDiff.finite_difference_jacobian(x -> contact_normal(collision, x, Quaternion(qa..., false), xb, Quaternion(qb..., false)), xa)
λ' * FiniteDiff.finite_difference_jacobian(q -> contact_normal(collision, xa, Quaternion(q..., false), xb, Quaternion(qb..., false)), vector(qa))

∂contact_normal_jvp∂x(:parent, collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false), λ)
∂contact_normal_jvp∂q(:parent, collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false), λ)

contact_tangent(collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false))

λ' * FiniteDiff.finite_difference_jacobian(x -> contact_tangent(collision, x, Quaternion(qa..., false), xb, Quaternion(qb..., false))[1, :]', xa)
λ' * FiniteDiff.finite_difference_jacobian(x -> contact_tangent(collision, x, Quaternion(qa..., false), xb, Quaternion(qb..., false))[2, :]', xa)

λ' * FiniteDiff.finite_difference_jacobian(q -> contact_tangent(collision, xa, Quaternion(q..., false), xb, Quaternion(qb..., false))[1, :], vector(qa))
λ' * FiniteDiff.finite_difference_jacobian(q -> contact_tangent(collision, xa, Quaternion(q..., false), xb, Quaternion(qb..., false))[2, :], vector(qa))

∂contact_tangent_jvp∂x(:parent, collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false), λ)
∂contact_tangent_jvp∂q(:parent, collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false), λ)

∂contact_tangent_jvp∂x(:child, collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false), λ)
∂contact_tangent_jvp∂q(:child, collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false), λ)

λ = rand(2)
∂contact_tangent_vjp∂x(:parent, collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false), λ)
∂contact_tangent_vjp∂q(:parent, collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false), λ)

∂contact_tangent_vjp∂x(:child, collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false), λ)
∂contact_tangent_vjp∂q(:child, collision, xa, Quaternion(qa..., false), xb, Quaternion(qb..., false), λ)


text: n = 9223372036854775809

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
struct vcovHetero <: vcov
    ssc::String
    function vcovHetero(; ssc::String = "HC1")
        new(ssc)
    end
end

function inference(N::Int,
                   K::Int,
                   resid::Vector,
                   β::Vector,
                   XX_inv::Matrix,
                   X::Matrix,
                   vcov::vcovHetero)

    resid_sq = resid.^2
    vcov_matrix = ((N-1)/(N-K)) * XX_inv * (X' * diagm(resid_sq) * X) * XX_inv

    se = sqrt.(diag(vcov_matrix))
    pval = 2 .* cdf.(TDist(N - K), -abs.(β ./ se))
    σ_sq = sum(resid.^2) / (N - K)
    return se, pval, σ_sq, vcov_matrix
end

max_stars_repo_name: EdJeeOnGitHub/jeeves

max_issues_repo_name: EdJeeOnGitHub/jeeves

max_forks_repo_name: EdJeeOnGitHub/jeeves


text:
using Test

@time @testset "ProbInts tests" begin include("probints.jl") end
@time @testset "Koopman Tests" begin include("koopman.jl") end


text:
# ---
# title: 1663. Smallest String With A Given Numeric Value
# id: problem1663
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/smallest-string-with-a-given-numeric-value/description/>
# hidden: true
# ---
#
# The **numeric value** of a **lowercase character** is defined as its position
# `(1-indexed)` in the alphabet, so the numeric value of `a` is `1`, the numeric
# value of `b` is `2`, the numeric value of `c` is `3`, and so on.
#
# The **numeric value** of a **string** consisting of lowercase characters is
# defined as the sum of its characters' numeric values. For example, the numeric
# value of the string `"abe"` is equal to `1 + 2 + 5 = 8`.
#
# You are given two integers `n` and `k`. Return _the **lexicographically
# smallest string** with **length** equal to `n` and **numeric value** equal to
# `k`._
#
# Note that a string `x` is lexicographically smaller than string `y` if `x`
# comes before `y` in dictionary order, that is, either `x` is a prefix of `y`,
# or if `i` is the first position such that `x[i] != y[i]`, then `x[i]` comes
# before `y[i]` in alphabetic order.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 3, k = 27
#     Output: "aay"
#     Explanation: The numeric value of the string is 1 + 1 + 25 = 27, and it is the smallest string with such a value and length equal to 3.
#
#
# **Example 2:**
#
#
#
#     Input: n = 5, k = 73
#     Output: "aaszz"
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 105`
#   * `n <= k <= 26 * n`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using PyPlot
using Convex
using PyCall
using HDF5, JLD
@pyimport matplotlib2tikz

save_plots = true
do_belief = true

plt[:close]("all")
pygui(!save_plots)

include("CenterOfMass.jl")
include("setup_parameters.jl")

type TrialData
  belief::Histogram
  estimate
  error
  applied_f
  f_hat
  csqmis
  robot_measurement
end

folder = "greedy_vs_random"

@load "$(folder)/data" data_array thetas

forces = map(x->x.applied_f, data_array)
errors = map(x->x.error, data_array)
n_measurement = size(data_array, 3)

colors = ["b" "g" "m"]

mean_error = mean(errors, 1)
for ii = 1:size(mean_error, 2)
  plot(1:n_measurement, mean_error[:,ii,:][:], color = colors[ii], linewidth = 3.0)
end
legend(["greedy", "random", "cyclic"])

indices = collect(1:n_measurement)
for ii = 1:size(errors,2)
  stds = std(errors[:,ii,:],1)[:]/sqrt(size(errors,1))
  es = mean_error[:,ii,:][:]
  fill([indices;reverse(indices)], [es+stds;reverse(es-stds)], color = colors[ii],
    alpha=0.2, linewidth=0.0)
end

if false
for ii = 1:size(errors, 2)
  for jj = 1:size(errors, 1)
    plot(1:n_measurement, errors[jj, ii,:][:], color = colors[ii], alpha = 0.1,
      linewidth = 1.0)
  end
end
end

plot([0;20],[0.1;0.1],linestyle="--", color="k")
plot([7;7],[0.0;0.1],linestyle="--", color="b")
plot([15;15],[0.0;0.1],linestyle="--", color="g")
plot([11;11],[0.0;0.1],linestyle="--", color="m")

xlabel("Iteration")
ylabel("Normalized error")

if save_plots
  matplotlib2tikz.save("$(folder)/convergence.tex", figureheight="\\figureheight",
    figurewidth="\\figurewidth")
end

if save_plots && do_belief
  for ii = 1:min(10,size(errors,1)) # num robots

    for jj = 1:size(errors,2) # function
      trial_folder = "$(folder)/fun_$(jj)_$(ii)"
      mkpath(trial_folder)
      for kk = 1:size(errors,3) # trial
        plt[:close]("all")

        println("belief: trial=$(ii), fun=$(jj), iter=$(kk)")

        data = data_array[ii,jj,kk]
        belief = data.belief
        csqmis = data.csqmis
        robot_measurement = data.robot_measurement
        theta = thetas[:,ii]
        com_p = theta[1:2]
        mass = theta[3]
        applied_p = attachment_ps[data.robot_measurement]

        figure()
        cloud = to_cloud(belief, interior_q)

        plot_attachment_csqmis(circle_ps, attachment_ps, csqmis)
        fill3d(hcat(circle_ps...), alpha=0.2, color="k")

        plot_solution(applied_p)

        scaled_ps = (15*cloud[4,:]'/maximum(cloud[4,:]))

        # plot probabilities
        scatter3D(cloud[1,:]', cloud[2,:]', cloud[3,:]', "z", scaled_ps.^2, "purple", alpha = 0.5)

        # plot CoM point
        scatter3D([com_p[1]], [com_p[2]], [mass], "z", 200, "red", marker="*", alpha=1)

        xlabel("X")
        ylabel("Y")
        zlabel("M")
        savefig("$(trial_folder)/belief_$(kk).png", pad_inches=0.01, bbox_inches="tight")
      end
    end
  end
end

pygui(true)

max_stars_repo_name: mcorah/CenterOfMass

max_issues_repo_name: mcorah/CenterOfMass

max_forks_repo_name: mcorah/CenterOfMass


text:
# This file contains the definition for the Storages struct.

mutable struct Storages
	unit
	node
	technology
	fuel
	power
	efficiency
	storage
	varcost

	function Storages(storages_df::DataFrame)

		S = Symbol.(storages_df[!, 1])
		storages_dict = df_to_dict_with_id(storages_df)

		return new(S,
				   storages_dict[:node],
				   storages_dict[:technology],
				   storages_dict[:fuel],
				   storages_dict[:power],
   				   storages_dict[:efficiency],
				   storages_dict[:storage],
				   storages_dict[:varcost])
	end

end

max_issues_repo_name: bobbyxiong/redispatch-ptg

max_forks_repo_name: bobbyxiong/redispatch-ptg


text:
module Chemistry

using JuLIP.FIO: load_dict
import JuLIP.FIO: read_dict, write_dict

export atomic_number,
       chemical_symbol,
       atomic_mass,
       symmetry,
       lattice_constant,
       rnn,
       AtomicNumber

"""
`struct AtomicNumber` : wrapper type to encode atomic numbers. Makes it
easier to dispatch on `AtomicNumber` and not confuse with a species index in
a list e.g.
"""
struct AtomicNumber <: Number
   z::Int16
end

AtomicNumber(z::AtomicNumber) = z
AtomicNumber(s::Symbol) = atomic_number(s)
(TI::Type{<: Integer})(z::AtomicNumber) = TI(z.z)
Base.convert(::Type{AtomicNumber}, z::Integer) = AtomicNumber(z)
Base.convert(::Type{AtomicNumber}, z::AtomicNumber) = z
Base.convert(TI::Type{<: Integer}, z::AtomicNumber) = convert(TI, z.z)
Base.hash(z::AtomicNumber, h::UInt) = hash(z.z, h)
Base.show(io::IO, z::AtomicNumber) = print(io, "<$(z.z)>")
Base.promote(x::AtomicNumber, y::Number) = promote(x.z, y)
Base.promote(x::Number, y::AtomicNumber) = promote(x, y.z)
Base.isless(x::AtomicNumber, y::AtomicNumber) = (x.z < y.z)
Base.isequal(x::AtomicNumber, y::AtomicNumber) = (x.z == y.z)

write_dict(z::AtomicNumber) =
      Dict("__id__" => "JuLIP_AtomicNumber", "z" => z.z)

read_dict(::Val{:JuLIP_AtomicNumber}, D::Dict) =
      AtomicNumber(D["z"])

data = load_dict(@__DIR__()[1:end-3] * "data/asedata.json")

const _rnn = [Float64(d) for d in data["rnn"]]
const _masses = [Float64(m) for m in data["masses"]]
const _refstates = Dict{String, Any}[ d == nothing ? Dict{String, Any}() : d
                                      for d in data["refstates"]]

# const _symbols = [Symbol(s) for s in data["symbols"]]
const _symbols = Dict{Int16, Symbol}()
const _numbers = Dict{Symbol, Int16}()
for (n, sym) in enumerate(Symbol.(data["symbols"]))
   _numbers[sym] = Int16(n - 1)
   _symbols[n-1] = sym
end

# special symbols
# :X -> 0 is already included in data["symbols"]
const __zAny__ = Int16(1000)
_numbers[:Any] = __zAny__
_symbols[__zAny__] = :Any

atomic_number(sym::Symbol) = AtomicNumber(_numbers[sym])

chemical_symbol(z::AtomicNumber) = _symbols[z.z]

atomic_mass(z::AtomicNumber) = _masses[z.z+1]
atomic_mass(sym::Symbol) = atomic_mass(atomic_number(sym))

element_name(z::AtomicNumber) = _names[z.z+1]
element_name(sym::Symbol) = element_name(atomic_number(sym))

rnn(z::AtomicNumber) = _rnn[z.z+1]
rnn(sym::Symbol) = rnn(atomic_number(sym))

symmetry(z::AtomicNumber) = Symbol(_refstates[z.z+1]["symmetry"])
symmetry(sym::Symbol) = symmetry(atomic_number(sym))

lattice_constant(z::AtomicNumber) = Float64( _refstates[z.z+1]["a"] )
lattice_constant(sym::Symbol) = lattice_constant(atomic_number(sym))

refstate(z::AtomicNumber) = _refstates[z.z+1]
refstate(sym::Symbol) = refstate(atomic_number(sym))

end


text:
import LinearAlgebra

using RungeKutta.Tableaus: get_lobatto_nodes, get_lobatto_weights,
                  get_lobatto_a_coefficients, get_lobatto_b_coefficients,
                  get_lobatto_c_coefficients, get_lobatto_c̄_coefficients,
                  get_lobatto_d_coefficients, get_lobatto_e_coefficients,
                  get_lobatto_f_coefficients, get_lobatto_g_coefficients

@testset "$(rpad("Lobatto Tableaus",80))" begin

    @test_throws ErrorException get_lobatto_nodes(1)
    @test_throws ErrorException get_lobatto_weights(1)
    @test_throws ErrorException get_lobatto_nullvector(1)
    @test_throws ErrorException get_lobatto_a_coefficients(1)
    @test_throws ErrorException get_lobatto_b_coefficients(1)
    @test_throws ErrorException get_lobatto_c_coefficients(1)
    @test_throws ErrorException get_lobatto_c̄_coefficients(1)
    @test_throws ErrorException get_lobatto_d_coefficients(1)
    @test_throws ErrorException get_lobatto_e_coefficients(1)
    @test_throws ErrorException get_lobatto_f_coefficients(1)
    @test_throws ErrorException get_lobatto_g_coefficients(1)

    function _get_lobatto_nullvector(s)
        if s == 2
            d = [+1.0, -1.0]
        elseif s == 3
            d = [+1.0, -2.0, +1.0]
        elseif s == 4
            d = [+1.0, -√5, +√5, -1.0]
        elseif s == 5
            d = [+3.0, -7.0, +8.0, -7.0, +3.0]
        else
            @error("We don't have a d vector for s=$(s) stages.")
        end
        return LinearAlgebra.normalize(d) * sign(d[begin])
    end

    @test get_lobatto_nullvector(2; normalize=true) ≈ _get_lobatto_nullvector(2)
    @test get_lobatto_nullvector(3; normalize=true) ≈ _get_lobatto_nullvector(3)
    @test get_lobatto_nullvector(4; normalize=true) ≈ _get_lobatto_nullvector(4)
    @test get_lobatto_nullvector(5; normalize=true) ≈ _get_lobatto_nullvector(5)

    function _getTableauLobattoIIIA2(T=Float64)
        a = BigFloat[
                [0     0   ]
                [1/2   1/2 ]
            ]

        Tableau{T}(:LobattoIIIA2, 2, a, get_lobatto_weights(2), get_lobatto_nodes(2))
    end

    function _getTableauLobattoIIIA3(T=Float64)
        a = BigFloat[
                [0      0     0    ]
                [5/24   1/3  -1/24 ]
                [1/6    2/3   1/6  ]
            ]

        Tableau{T}(:LobattoIIIA3, 4, a, get_lobatto_weights(3), get_lobatto_nodes(3))
    end

    function _getTableauLobattoIIIA4(T=Float64)
        a = BigFloat[
                [0            0               0               0          ]
                [(11+√5)/120  (25-   √5)/120  (25-13*√5)/120  (-1+√5)/120]
                [(11-√5)/120  (25+13*√5)/120  (25+   √5)/120  (-1-√5)/120]
                [      1/12            5/12            5/12         1/12 ]
            ]

        Tableau{T}(:LobattoIIIA4, 6, a, get_lobatto_weights(4), get_lobatto_nodes(4))
    end

    function _getTableauLobattoIIIA5(T=Float64)
        a = BigFloat[
                [0                 0                   0                  0                   0               ]
                [(119+3*√21)/1960  (343-  9*√21)/2520  (392-96*√21)/2205  (343- 69*√21)/2520  (-21+3*√21)/1960]
                [         13/320   (392+105*√21)/2880             8/45    (392-105*√21)/2880            3/320 ]
                [(119-3*√21)/1960  (343+ 69*√21)/2520  (392+96*√21)/2205  (343+  9*√21)/2520  (-21-3*√21)/1960]
                [          1/20               49/180             16/45               49/180             1/20  ]
            ]

        Tableau{T}(:LobattoIIIA5, 8, a, get_lobatto_weights(5), get_lobatto_nodes(5))
    end

    function _getTableauLobattoIIIB2(T=Float64)
        a = BigFloat[
                [1/2  0]
                [1/2  0]
            ]

        Tableau{T}(:LobattoIIIB2, 2, a, get_lobatto_weights(2), get_lobatto_nodes(2))
    end

    function _getTableauLobattoIIIB3(T=Float64)
        a = BigFloat[
                [1/6  -1/6   0   ]
                [1/6   1/3   0   ]
                [1/6   5/6   0   ]
            ]

        Tableau{T}(:LobattoIIIB3, 4, a, get_lobatto_weights(3), get_lobatto_nodes(3))
    end

    function _getTableauLobattoIIIB4(T=Float64)
        a = BigFloat[
                [ 1/12  (-1-   √5)/24   (-1+   √5)/24    0 ]
                [ 1/12  (25+   √5)/120  (25-13*√5)/120   0 ]
                [ 1/12  (25+13*√5)/120  (25-   √5)/120   0 ]
                [ 1/12  (11-   √5)/24   (11+   √5)/24    0 ]
            ]

        Tableau{T}(:LobattoIIIB4, 6, a, get_lobatto_weights(4), get_lobatto_nodes(4))
    end

    function _getTableauLobattoIIIB5(T=Float64)
        a = BigFloat[
                [ 1/20  ( -7-   √21)/120             1/15   ( -7+   √21)/120    0 ]
                [ 1/20  (343+ 9*√21)/2520  (56-15*√21)/315  (343-69*√21)/2520   0 ]
                [ 1/20  ( 49+12*√21)/360             8/45   ( 49-12*√21)/360    0 ]
                [ 1/20  (343+69*√21)/2520  (56+15*√21)/315  (343- 9*√21)/2520   0 ]
                [ 1/20  (119- 3*√21)/360            13/45   (119+ 3*√21)/360    0 ]
            ]

        Tableau{T}(:LobattoIIIB5, 8, a, get_lobatto_weights(5), get_lobatto_nodes(5))
    end

    function _getTableauLobattoIIIC2(T=Float64)
        a = BigFloat[
                [1/2  -1/2 ]
                [1/2   1/2 ]
            ]

        Tableau{T}(:LobattoIIIC2, 2, a, get_lobatto_weights(2), get_lobatto_nodes(2))
    end

    function _getTableauLobattoIIIC3(T=Float64)
        a = BigFloat[
                [1/6  -1/3    1/6  ]
                [1/6   5/12  -1/12 ]
                [1/6   2/3    1/6  ]
            ]

        Tableau{T}(:LobattoIIIC3, 4, a, get_lobatto_weights(3), get_lobatto_nodes(3))
    end

    function _getTableauLobattoIIIC4(T=Float64)
        a = BigFloat[
                [ 1/12        -√5/12         √5/12   -1/12 ]
                [ 1/12          1/4   (10-7*√5)/60   √5/60 ]
                [ 1/12  (10+7*√5)/60          1/4   -√5/60 ]
                [ 1/12          5/12          5/12    1/12 ]
            ]

        Tableau{T}(:LobattoIIIC4, 6, a, get_lobatto_weights(4), get_lobatto_nodes(4))
    end

    function _getTableauLobattoIIIC5(T=Float64)
        a = BigFloat[
                [ 1/20            -21/180             2/15             -21/180     1/20  ]
                [ 1/20             29/180   (47-15*√21)/315  (203- 30*√21)/1260   -3/140 ]
                [ 1/20  (329+105*√21)/2880           73/360  (329-105*√21)/2880    3/160 ]
                [ 1/20  (203+ 30*√21)/1260  (47+15*√21)/315             29/180    -3/140 ]
                [ 1/20             49/180            16/45              49/180     1/20  ]
            ]

        Tableau{T}(:LobattoIIIC5, 8, a, get_lobatto_weights(5), get_lobatto_nodes(5))
    end

    function _getTableauLobattoIIIC̄2(T=Float64)
        a = BigFloat[
                [0  0]
                [1  0]
            ]

        Tableau{T}(:LobattoIIIC̄2, 2, a, get_lobatto_weights(2), get_lobatto_nodes(2))
    end

    function _getTableauLobattoIIIC̄3(T=Float64)
        a = BigFloat[
                [0    0    0 ]
                [1/4  1/4  0 ]
                [0    1    0 ]
            ]

        Tableau{T}(:LobattoIIIC̄3, 4, a, get_lobatto_weights(3), get_lobatto_nodes(3))
    end

    function _getTableauLobattoIIIC̄4(T=Float64)
        a = BigFloat[
                [      0             0             0     0 ]
                [ (5+√5)/60          1/6   (15-7*√5)/60  0 ]
                [ (5-√5)/60  (15+7*√5)/60          1/6   0 ]
                [      1/6      (5-√5)/12     (5+√5)/12  0 ]
            ]

        Tableau{T}(:LobattoIIIC̄4, 6, a, get_lobatto_weights(4), get_lobatto_nodes(4))
    end

    function _getTableauLobattoIIIC̄5(T=Float64)
        a = BigFloat[
                [ 0     0                0              0                0 ]
                [ 1/14            1/9    (13-3*√21)/63  (14- 3*√21)/126  0 ]
                [ 1/32  (91+21*√21)/576          11/72  (91-21*√21)/576  0 ]
                [ 1/14  (14+ 3*√21)/126  (13+3*√21)/63            1/9    0 ]
                [ 0               7/18            2/9             7/18   0 ]
            ]

        Tableau{T}(:LobattoIIIC̄5, 8, a, get_lobatto_weights(5), get_lobatto_nodes(5))
    end

    function _getTableauLobattoIIID2(T=Float64)
        lobC = _getTableauLobattoIIIC2(BigFloat)
        lobC̄ = _getTableauLobattoIIIC̄2(BigFloat)
        Tableau{T}(:LobattoIIID2, lobC.o, (lobC̄.a + lobC.a)/2, lobC.b, lobC.c)
    end

    function _getTableauLobattoIIID3(T=Float64)
        lobC = _getTableauLobattoIIIC3(BigFloat)
        lobC̄ = _getTableauLobattoIIIC̄3(BigFloat)
        Tableau{T}(:LobattoIIID3, lobC.o, (lobC̄.a + lobC.a)/2, lobC.b, lobC.c)
    end

    function _getTableauLobattoIIID4(T=Float64)
        lobC = _getTableauLobattoIIIC4(BigFloat)
        lobC̄ = _getTableauLobattoIIIC̄4(BigFloat)
        Tableau{T}(:LobattoIIID4, lobC.o, (lobC̄.a + lobC.a)/2, lobC.b, lobC.c)
    end

    function _getTableauLobattoIIID5(T=Float64)
        lobC = _getTableauLobattoIIIC5(BigFloat)
        lobC̄ = _getTableauLobattoIIIC̄5(BigFloat)
        Tableau{T}(:LobattoIIID5, lobC.o, (lobC̄.a + lobC.a)/2, lobC.b, lobC.c)
    end

    function _getTableauLobattoIIIE2(T=Float64)
        lobA = _getTableauLobattoIIIA2(BigFloat)
        lobB = _getTableauLobattoIIIB2(BigFloat)
        Tableau{T}(:LobattoIIIE2, lobA.o, (lobA.a + lobB.a)/2, lobA.b, lobA.c)
    end

    function _getTableauLobattoIIIE3(T=Float64)
        lobA = _getTableauLobattoIIIA3(BigFloat)
        lobB = _getTableauLobattoIIIB3(BigFloat)
        Tableau{T}(:LobattoIIIE3, lobA.o, (lobA.a + lobB.a)/2, lobA.b, lobA.c)
    end

    function _getTableauLobattoIIIE4(T=Float64)
        lobA = _getTableauLobattoIIIA4(BigFloat)
        lobB = _getTableauLobattoIIIB4(BigFloat)
        Tableau{T}(:LobattoIIIE4, lobA.o, (lobA.a + lobB.a)/2, lobA.b, lobA.c)
    end

    function _getTableauLobattoIIIE5(T=Float64)
        lobA = _getTableauLobattoIIIA5(BigFloat)
        lobB = _getTableauLobattoIIIB5(BigFloat)
        Tableau{T}(:LobattoIIIE5, lobA.o, (lobA.a + lobB.a)/2, lobA.b, lobA.c)
    end

    function _getTableauLobattoIIIF2(T=Float64)
        a = BigFloat[
                [1/12  -1/12 ]
                [7/12   5/12 ]
            ]

        Tableau{T}(:LobattoIIIF2, 4, a, get_lobatto_weights(2), get_lobatto_nodes(2))
    end

    function _getTableauLobattoIIIF3(T=Float64)
        a = BigFloat[
                [1/30   -1/15    1/30 ]
                [5/24    1/3    -1/24 ]
                [2/15   11/15    2/15 ]
            ]

        Tableau{T}(:LobattoIIIF3, 6, a, get_lobatto_weights(3), get_lobatto_nodes(3))
    end

    function _getTableauLobattoIIIF4(T=Float64)
        a = BigFloat[
                [  1/56                 -√5/56           √5/56   -1/56         ]
                [ 37/420+√5/120  5/24-   √5/210  5/24-47*√5/420  -1/210+√5/120 ]
                [ 37/420-√5/120  5/24+47*√5/420  5/24+   √5/210  -1/210-√5/120 ]
                [ 17/168         5/12-   √5/56   5/12+   √5/56   11/168        ]
            ]

        Tableau{T}(:LobattoIIIF4, 8, a, get_lobatto_weights(4), get_lobatto_nodes(4))
    end

    function _getTableauLobattoIIIG2(T=Float64)
        symplecticize(_getTableauLobattoIIIF2(BigFloat); name=:LobattoIIIG2, T=T)
    end

    function _getTableauLobattoIIIG3(T=Float64)
        symplecticize(_getTableauLobattoIIIF3(BigFloat); name=:LobattoIIIG3, T=T)
    end

    function _getTableauLobattoIIIG4(T=Float64)
        symplecticize(_getTableauLobattoIIIF4(BigFloat); name=:LobattoIIIG4, T=T)
    end

    @test_throws ErrorException TableauLobattoIIIA(1)
    @test_throws ErrorException TableauLobattoIIIB(1)
    @test_throws ErrorException TableauLobattoIIIC(1)
    @test_throws ErrorException TableauLobattoIIIC̄(1)
    @test_throws ErrorException TableauLobattoIIID(1)
    @test_throws ErrorException TableauLobattoIIIE(1)
    @test_throws ErrorException TableauLobattoIIIF(1)
    @test_throws ErrorException TableauLobattoIIIG(1)

    @test TableauLobattoIIIA(2) ≈ _getTableauLobattoIIIA2()
    @test TableauLobattoIIIA(3) ≈ _getTableauLobattoIIIA3()
    @test TableauLobattoIIIA(4) ≈ _getTableauLobattoIIIA4()
    @test TableauLobattoIIIA(5) ≈ _getTableauLobattoIIIA5()

    @test TableauLobattoIIIB(2) ≈ _getTableauLobattoIIIB2()
    @test TableauLobattoIIIB(3) ≈ _getTableauLobattoIIIB3()
    @test TableauLobattoIIIB(4) ≈ _getTableauLobattoIIIB4()
    @test TableauLobattoIIIB(5) ≈ _getTableauLobattoIIIB5()

    @test TableauLobattoIIIC(2) ≈ _getTableauLobattoIIIC2()
    @test TableauLobattoIIIC(3) ≈ _getTableauLobattoIIIC3()
    @test TableauLobattoIIIC(4) ≈ _getTableauLobattoIIIC4()
    @test TableauLobattoIIIC(5) ≈ _getTableauLobattoIIIC5()

    @test TableauLobattoIIIC̄(2) ≈ _getTableauLobattoIIIC̄2()
    @test TableauLobattoIIIC̄(3) ≈ _getTableauLobattoIIIC̄3()
    @test TableauLobattoIIIC̄(4) ≈ _getTableauLobattoIIIC̄4()
    @test TableauLobattoIIIC̄(5) ≈ _getTableauLobattoIIIC̄5()

    @test TableauLobattoIIID(2) ≈ _getTableauLobattoIIID2()
    @test TableauLobattoIIID(3) ≈ _getTableauLobattoIIID3()
    @test TableauLobattoIIID(4) ≈ _getTableauLobattoIIID4()
    @test TableauLobattoIIID(5) ≈ _getTableauLobattoIIID5()

    @test TableauLobattoIIIE(2) ≈ _getTableauLobattoIIIE2()
    @test TableauLobattoIIIE(3) ≈ _getTableauLobattoIIIE3()
    @test TableauLobattoIIIE(4) ≈ _getTableauLobattoIIIE4()
    @test TableauLobattoIIIE(5) ≈ _getTableauLobattoIIIE5()

    @test TableauLobattoIIIF(2) ≈ _getTableauLobattoIIIF2()
    @test TableauLobattoIIIF(3) ≈ _getTableauLobattoIIIF3()
    @test TableauLobattoIIIF(4) ≈ _getTableauLobattoIIIF4()

    @test TableauLobattoIIIG(2) ≈ _getTableauLobattoIIIG2()
    @test TableauLobattoIIIG(3) ≈ _getTableauLobattoIIIG3()
    @test TableauLobattoIIIG(4) ≈ _getTableauLobattoIIIG4()

end


text:
using Test
using SparseArrays: sparsevec
using ..ChemistryFeaturization.Utils.OrbitalFeatureUtils:
    _indexorbital, _orbitalindex, _orbitalsparse, _econf_to_name, valenceshell_conf_df

@testset "OrbitalFeatureUtils" begin
    @testset "Index - Orbital" begin
        orbitals =
            SubString.([
                "1s",
                "2s",
                "2p",
                "3s",
                "3p",
                "4s",
                "3d",
                "4p",
                "5s",
                "4d",
                "5p",
                "6s",
                "4f",
                "5d",
                "6p",
                "7s",
                "5f",
                "6d",
                "7p",
                "8s",
                "5g",
                "6f",
                "7d",
                "8p",
                "9s",
            ])
        count = [i for i in Int16.(1:25)]

        @test _indexorbital.(count) == orbitals
        @test _orbitalindex.(orbitals) == count
    end

    @testset "_orbitalsparse" begin
        # in the expected result, the first vector the positions and the second is the values.
        @test _orbitalsparse("1s2") == ([1], [2]) # 1s2
        @test _orbitalsparse("3d5.4s1") == ([7, 6], [5, 1]) # 3d5, 4s1.
    end

    @testset "Element Name - Electronic Configuration" begin

        get_config(el) = _orbitalsparse(
            getproperty(
                valenceshell_conf_df[valenceshell_conf_df.Symbol.==el, :][1, :],
                Symbol("Electronic Structure"),
            ),
        )

        function test_econf_to_name(I, V, element)
            @test _econf_to_name(sparsevec(Int16.(I), Int16.(V))) == element
        end

        I_H, V_H = get_config("H")
        @test (I_H, V_H) == ([1], [1])
        test_econf_to_name(I_H, V_H, "H")

        I_He, V_He = get_config("He")
        @test (I_He, V_He) == ([1], [2])
        test_econf_to_name(I_He, V_He, "He")

        I_Li, V_Li = get_config("Li")
        @test (I_Li, V_Li) == ([2], [1])
        test_econf_to_name(I_He, V_He, "He")
    end
end


text:
using Libdl

using BinaryProvider

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS

issource_build = "ARBNUMERICS_SOURCE_BUILD" in keys(ENV) && ENV["ARBNUMERICS_SOURCE_BUILD"] == "1"

const prefixpath = joinpath(@__DIR__, "usr")

const wdir = joinpath(@__DIR__)

if !issource_build
  # Dependencies that must be installed before this package can be built
  dependencies = [
    # This has to be in sync with the corresponding commit in the source build below (for flint, arb, antic)
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/GMP-v6.1.2-1/build_GMP.v6.1.2.jl",
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/MPFR-v4.0.2-1/build_MPFR.v4.0.2.jl",
    "https://github.com/thofma/Flint2Builder/releases/download/f46562/build_libflint.v0.0.0-f465622699d5c4c22bb3617596f8ae86e4570652.jl",
    "https://github.com/thofma/ArbBuilder/releases/download/6c3738-v2/build_libarb.v0.0.0-6c3738555d00b8b8b24a1f5e0065ef787432513c.jl",
    "https://github.com/thofma/AnticBuilder/releases/download/364f97-v2/build_libantic.v0.0.0-364f97edd9b6af537787113cf910f16d7bbc48a3.jl"
   ]

  const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))

  products = []

  for url in dependencies
      build_file = joinpath(@__DIR__, basename(url))
      if !isfile(build_file)
          download(url, build_file)
      end
  end

  # Execute the build scripts for the dependencies in an isolated module to avoid overwriting
  # any variables/constants here
  for url in dependencies
      build_file = joinpath(@__DIR__, basename(url))
      m = @eval module $(gensym()); include($build_file); end
      append!(products, m.products)
  end

  filenames = ["libgmp.la", "libgmpxx.la", "libmpfr.la"]
  for filename in filenames
    fpath = joinpath(prefixpath, "lib", filename)
    txt = read(fpath, String)
    open(fpath, "w") do f
      write(f, replace(txt, "/workspace/destdir" => prefixpath))
    end
  end

else
  println("Doing a source build for C dependencies...")
  if "ARBNUMERICS_BUILD_THREADS" in keys(ENV)
     build_threads = ENV["ARBNUMERICS_BUILD_THREADS"]
     println("Using $build_threads threads for building as specified by ARBNUMERICS_BUILD_THREADS.")
  else
     build_threads = Sys.CPU_THREADS
     println("Using $build_threads threads (detected that many CPU threads).")
  end

  @show M4_VERSION = "1.4.17"
  @show YASM_VERSION = "1.3.0"
  @show MPIR_VERSION = "3.0.0-90740d8fdf03b941b55723b449831c52fd7f51ca"
  @show MPFR_VERSION = "4.0.0"
  @show FLINT_VERSION = "f465622699d5c4c22bb3617596f8ae86e4570652"
  @show ARB_VERSION = "6c3738555d00b8b8b24a1f5e0065ef787432513c"

  if Sys.iswindows()
    error("Source build not available on Windows")
  end

  println("Removing old binaries ...")

  rm(prefixpath, force = true, recursive = true)
  mkdir(prefixpath)
  mkdir(joinpath(prefixpath, "lib"))

  if Sys.isapple() && !("CC" in keys(ENV))
   ENV["CC"] = "clang"
   ENV["CXX"] = "clang++"
  end

  LDFLAGS = "-Wl,-rpath,$prefixpath/lib -Wl,-rpath,\$\$ORIGIN/../share/julia/site/v$(VERSION.major).$(VERSION.minor)/Nemo/local/lib"
  DLCFLAGS = "-fPIC -fno-common"

  cd(wdir)

  try
    run(`m4 --version`)
  catch
    println("Building m4 ... ")
    M4_FILE = "m4-" * M4_VERSION * ".tar.bz2"
    download("http://ftp.gnu.org/gnu/m4/$M4_FILE", joinpath(wdir, "$M4_FILE"))

      run(`tar -xvf $M4_FILE`)
      run(`rm $M4_FILE`)
      cd(joinpath("$wdir", "m4-$M4_VERSION"))
      run(`./configure --prefix=$prefixpath`)
      run(`make`)
      run(`make install`)
      println("DONE")
   end

  cd(wdir)

  # install yasm

  if !ispath(joinpath(wdir, "yasm-$YASM_VERSION"))
     println("Building yasm ... ")
     YASM_FILE = "yasm-" * YASM_VERSION * ".tar.gz"
     download("https://github.com/yasm/yasm/archive/v$(YASM_VERSION).tar.gz", YASM_FILE)
     run(`tar -xvf $YASM_FILE`)
     run(`rm $YASM_FILE`)
     cd(joinpath("$wdir","yasm-$YASM_VERSION"))
     run(`./autogen.sh`)
     run(`./configure`)
     run(`make`)
     println("DONE")
  end

  cd(wdir)

  # install GMP/MPIR

  MPIR_FILE = "mpir-" * MPIR_VERSION * ".tar.bz2"

  if !ispath(joinpath(wdir, "mpir-$MPIR_VERSION"))
    println("Downloading MPIR sources ... ")
    download("http://nemocas.org/binaries/$MPIR_FILE", joinpath(wdir, MPIR_FILE))
    println("DONE")
  end

  println("Building MPIR ... ")
  if isfile(joinpath(wdir, MPIR_FILE))
     run(`tar -xvf $MPIR_FILE`)
     run(`rm $MPIR_FILE`)
  end
  cd("$wdir/mpir-$MPIR_VERSION")
  try
     run(`m4 --version`)
     run(`./configure --with-yasm=$wdir/yasm-$YASM_VERSION/yasm --prefix=$prefixpath --enable-gmpcompat --disable-static --enable-shared`)
  catch
     run(`./configure --with-yasm=$wdir/yasm-$YASM_VERSION/yasm --prefix=$prefixpath M4=$prefixpath/bin/m4 --enable-gmpcompat --disable-static --enable-shared`)
  end
  run(`make -j$build_threads`)
  run(`make install`)
  cd(wdir)
  run(`rm -rf bin`)
  println("DONE")

  cd(wdir)

  # install MPFR

  MPFR_FILE = "mpfr-" * MPFR_VERSION * ".tar.bz2"

  if !ispath(joinpath(wdir, "mpfr-$MPFR_VERSION"))
     println("Downloading MPFR sources ... ")
     download("http://ftp.vim.org/ftp/gnu/mpfr/$MPFR_FILE", joinpath(wdir, MPFR_FILE))

     println("DONE")
  end

  println("Building MPFR ... ")
  if isfile(joinpath(wdir, MPFR_FILE))
    run(`tar -xvf $MPFR_FILE`)
    run(`rm $MPFR_FILE`)
  end
  cd("$wdir/mpfr-$MPFR_VERSION")
  withenv("LD_LIBRARY_PATH"=>"$prefixpath/lib", "LDFLAGS"=>LDFLAGS) do
    run(`./configure --prefix=$prefixpath --with-gmp=$prefixpath --disable-static --enable-shared`)
    run(`make -j$build_threads`)
    run(`make install`)
  end
  println("DONE")

  cd(wdir)

  # INSTALL FLINT

  try
    println("Cloning flint2 ... ")
    run(`git clone https://github.com/wbhart/flint2.git`)
    cd(joinpath("$wdir", "flint2"))
    run(`git checkout $FLINT_VERSION`)
    cd(wdir)
  catch
    if ispath(joinpath("$wdir", "flint2"))
      open(`patch -R --forward -d flint2 -r -`, "r", open("../deps-PIE-ftbfs.patch"))
      cd(joinpath("$wdir", "flint2"))
      run(`git fetch`)
      run(`git checkout $FLINT_VERSION`)
      cd(wdir)
    end
  end
  open(`patch --forward -d flint2 -r -`, "r", open("../deps-PIE-ftbfs.patch"))
  println("DONE")

  println("Building flint ... ")
  cd(joinpath("$wdir", "flint2"))
  withenv("LD_LIBRARY_PATH"=>"$prefixpath/lib", "LDFLAGS"=>LDFLAGS) do
    run(`./configure --prefix=$prefixpath --disable-static --enable-shared --with-mpir=$prefixpath --with-mpfr=$prefixpath`)
    run(`make -j$build_threads`)
    run(`make install`)
  end

  println("DONE")

  cd(wdir)

  # INSTALL ARB

  println("Cloning arb ... ")
  try
    run(`git clone https://github.com/fredrik-johansson/arb.git`)
    cd(joinpath("$wdir", "arb"))
    run(`git checkout $ARB_VERSION`)
    cd(wdir)
  catch
    if ispath(joinpath("$wdir", "arb"))
      #open(`patch -R --forward -d arb -r -`, "r", open("../deps-PIE-ftbfs.patch"))
      cd(joinpath("$wdir", "arb"))
      run(`git fetch`)
      run(`git checkout $ARB_VERSION`)
      cd(wdir)
    end
  end
  println("DONE")

  println("Building arb ... ")
  cd(joinpath("$wdir", "arb"))
  withenv("LD_LIBRARY_PATH"=>"$prefixpath/lib", "LDFLAGS"=>LDFLAGS) do
    run(`./configure --prefix=$prefixpath --disable-static --enable-shared --with-mpir=$prefixpath --with-mpfr=$prefixpath --with-flint=$prefixpath`)
    run(`make -j$build_threads`)
    run(`make install`)
  end
  println("DONE")

  cd(wdir)
end

push!(Libdl.DL_LOAD_PATH, joinpath(prefixpath, "lib"), joinpath(prefixpath, "bin"))


text:
import JSON3, HTTP
import UUIDs: UUID
import DataFrames: DataFrame
export list_notebooks, get_notebook_html, get_notebook_state, get_notebook_code, publish_notebook, update_notebook, delete_notebook

UUIDish = Union{UUID, AbstractString}

function endpoint(path::AbstractString; url="https://plutopublish.com")
    return "$(url)$(path)"
end

function query_string(properties::Dict)
    parts = []

    for (prop_key, prop_val) ∈ properties
        push!(parts, HTTP.escapeuri(prop_key) * "=" * HTTP.escapeuri(prop_val))
    end

    "?" * join(parts, "&")
end

function list_notebooks(; kw...)
    req = HTTP.get(endpoint("/notebooks"; kw...))
    JSON3.read(req.body) |> DataFrame
end
function get_notebook_html(id::UUIDish; kw...)
    req = HTTP.get(endpoint("/notebooks/$(id)/html"; kw...))
    String(req.body)
end
function get_notebook_state(id::UUIDish; kw...)
    req = HTTP.get(endpoint("/notebooks/$(id)/state"; kw...))
    req.body
end
function get_notebook_code(id::UUIDish; kw...)
    req = HTTP.get(endpoint("/notebooks/$(id)/code"; kw...))
    String(req.body)
end
function publish_notebook(notebook_path::AbstractString, notebook_name::AbstractString; apikey=nothing, kw...)
    open(notebook_path) do notebook_io
        body_data = Dict(
            "notebook" => HTTP.Multipart(notebook_name, notebook_io, "application/x-pluto-statefile")
        )

        body = HTTP.Form(body_data)
        headers = []
        if !isnothing(apikey)
            push!(headers, "X-Api-Key" => apikey)
        end

        req = HTTP.post(endpoint("/notebooks" * (isnothing(notebook_name) ? "" : "?name=$(HTTP.escapeuri(notebook_name))"); kw...), headers, body)
        JSON3.read(req.body) |> DataFrame
    end
end
function update_notebook(id::UUIDish, patch_data::Dict; apikey=nothing, kw...)
    notebook_file = if haskey(patch_data, "notebook")
        nbfile = patch_data["notebook"]
        delete!(patch_data, "notebook")
        nbfile
    else
        nothing
    end

    notebook_body = if !isnothing(notebook_file)
        open(notebook_file, "r") do io
            read(io)
        end
    else
        nothing
    end

    headers = []
    if !isnothing(apikey)
        push!(headers, "X-Api-Key" => apikey)
    end

    req = HTTP.patch(endpoint("/notebooks/$(id)" * query_string(patch_data); kw...), headers, isnothing(notebook_body) ? UInt8[] : notebook_body)
    JSON3.read(req.body) |> DataFrame
end
function delete_notebook(id::UUIDish; apikey=nothing, kw...)
    headers = []
    if !isnothing(apikey)
        push!(headers, "X-Api-Key" => apikey)
    end

    req = HTTP.delete(endpoint("/notebooks/$(id)"; kw...), headers)
    req.status
end


text:
module InterfaceTest

using Test
using ForneyLab: Interface, FactorNode

# Integration helper
mutable struct MockNode <: FactorNode end

@testset "Interface" begin
    # Interface should construct
    node = MockNode()
    iface = Interface(node)
    @test ===(iface.node, node)
    @test iface.edge == nothing
    @test iface.partner == nothing
end

end # module


text:
# eachslice, eachcol, eachrow exactly as in Julia 1.1
# https://github.com/JuliaLang/julia/pull/29749/files
# https://github.com/JuliaLang/julia/blob/master/HISTORY.md#new-library-functions

export eachrow, eachcol, eachslice

"""
    eachrow(A::AbstractVecOrMat)
Create a generator that iterates over the first dimension of vector or matrix `A`,
returning the rows as views.
See also [`eachcol`](@ref) and [`eachslice`](@ref).
"""
eachrow(A::AbstractVecOrMat) = (view(A, i, :) for i in axes(A, 1))

"""
    eachcol(A::AbstractVecOrMat)
Create a generator that iterates over the second dimension of matrix `A`, returning the
columns as views.
See also [`eachrow`](@ref) and [`eachslice`](@ref).
"""
eachcol(A::AbstractVecOrMat) = (view(A, :, i) for i in axes(A, 2))

"""
    eachslice(A::AbstractArray; dims)
Create a generator that iterates over dimensions `dims` of `A`, returning views that select all
the data from the other dimensions in `A`.
Only a single dimension in `dims` is currently supported. Equivalent to `(view(A,:,:,...,i,:,:
...)) for i in axes(A, dims))`, where `i` is in position `dims`.
See also [`eachrow`](@ref), [`eachcol`](@ref), and [`selectdim`](@ref).
"""
@inline function eachslice(A::AbstractArray; dims)
    length(dims) == 1 || throw(ArgumentError("only single dimensions are supported"))
    dim = first(dims)
    dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
    idx1, idx2 = ntuple(d->(:), dim-1), ntuple(d->(:), ndims(A)-dim)
    return (view(A, idx1..., i, idx2...) for i in axes(A, dim))
end

# also from 1.1, like === not ==
isnothing(::Any) = false
isnothing(::Nothing) = true


text:
function toZXY(a::AbstractArray{T, 3}, old_arrange::Symbol) where T
    if old_arrange == :zxy
        return a
    elseif old_arrange == :xyz
        return PermutedDimsArray(a, (3, 1, 2))
    else
        throw(ErrorException("Unknown arrangement: " * old_arrange))
    end
end

function toXYZ(a::AbstractArray{T, 3}, old_arrange::Symbol) where T
    if old_arrange == :xyz
        return a
    elseif old_arrange == :zxy
        return PermutedDimsArray(a, (2, 3, 1))
    else
        throw(ErrorException("Unknown arrangement: " * old_arrange))
    end
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
function get_tensorinfo(tensor_shape)
    tensor_info = ""
    tensor_info_for_title = ""
    for s in tensor_shape
        tensor_info = tensor_info*"_"*string(lpad(s, n_digits, "0"))
    end
    return tensor_info
end

function get_dkl(P, Q ; eps_thred = 1.0e-10)
    tensor_shape = size(P)
    dkl = 0.0
    for idx in CartesianIndices(tensor_shape)
        if P[idx] < eps_thred
            continue
        end
        if Q[idx] < eps_thred
            continue
        end
        dkl += P[idx] * log( P[idx] / Q[idx] )
    end
    return dkl - sum(P) + sum(Q)
    #return sum( P .* log.( P ./ Q ) ) - sum(P) + sum(Q)
end

max_stars_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_issues_repo_name: gkazunii/Legendre-tucker-rank-reduction

max_forks_repo_name: gkazunii/Legendre-tucker-rank-reduction


text:
using CUDAnative
using CuArrays
CuArrays.allowscalar(false)
try
  isdefined(NumericalFluxes, :rusanov!)
catch
  include("../../src/DGmethods/NumericalFluxes.jl")
end

include("DGBalanceLawProfileHarness.jl")

const _nstate = 5
const _ρ, _U, _V, _W, _E = 1:_nstate
const stateid = (ρid = _ρ, Uid = _U, Vid = _V, Wid = _W, Eid = _E)
const statenames = ("ρ", "U", "V", "W", "E")
const γ_exact = 7 // 5

# preflux computation
@inline function preflux(Q, _...)
  γ::eltype(Q) = γ_exact
  @inbounds ρ, U, V, W, E = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E]
  ρinv = 1 / ρ
  u, v, w = ρinv * U, ρinv * V, ρinv * W
  ((γ-1)*(E - ρinv * (U^2 + V^2 + W^2) / 2), u, v, w, ρinv)
end

# max eigenvalue
@inline function wavespeed(n, Q, aux, t, P, u, v, w, ρinv)
  γ::eltype(Q) = γ_exact
  @inbounds abs(n[1] * u + n[2] * v + n[3] * w) + sqrt(ρinv * γ * P)
end

# physical flux function
eulerflux!(F, Q, QV, aux, t) =
eulerflux!(F, Q, QV, aux, t, preflux(Q)...)

@inline function eulerflux!(F, Q, QV, aux, t, P, u, v, w, ρinv)
  @inbounds begin
    ρ, U, V, W, E = Q[_ρ], Q[_U], Q[_V], Q[_W], Q[_E]

    F[1, _ρ], F[2, _ρ], F[3, _ρ] = U          , V          , W
    F[1, _U], F[2, _U], F[3, _U] = u * U  + P , v * U      , w * U
    F[1, _V], F[2, _V], F[3, _V] = u * V      , v * V + P  , w * V
    F[1, _W], F[2, _W], F[3, _W] = u * W      , v * W      , w * W + P
    F[1, _E], F[2, _E], F[3, _E] = u * (E + P), v * (E + P), w * (E + P)
  end
end

numerical_flux!(x...) = NumericalFluxes.rusanov!(x..., eulerflux!, wavespeed,
                                                 preflux)

let
  DFloat = Float64
  dim = 2
  nelem = 100
  N = 4
  cpu_dg = DGProfiler(Array, DFloat, dim, nelem, N, _nstate, eulerflux!,
                      numerical_flux!; stateoffset = ((_E, 20), (_ρ, 1)))
  gpu_dg = DGProfiler(CuArray, DFloat, dim, nelem, N, _nstate, eulerflux!,
                      numerical_flux!; stateoffset = ((_E, 20), (_ρ, 1)))
  volumerhs!(cpu_dg)
  volumerhs!(gpu_dg)
  facerhs!(cpu_dg)
  facerhs!(gpu_dg)
  gpu_rhs = Array(gpu_dg.rhs)
  cpu_rhs = Array(cpu_dg.rhs)
  @show gpu_rhs ≈ cpu_rhs
  nothing
end

max_stars_repo_name: jkozdon/CLIMA

max_issues_repo_name: jkozdon/CLIMA

max_forks_repo_name: jkozdon/CLIMA


text:
function append_subject end;
function append_subject(acs::AccessControlSystem, u::Int = 1)
    #acs.subject
end;

max_stars_repo_name: daviddoret/julia-cybersecurity

max_issues_repo_name: daviddoret/julia-cybersecurity

max_forks_repo_name: daviddoret/julia-cybersecurity


text:
@testset "Sort2MILP: $(fct_type), dimension $(array_dim), $(T)" for fct_type in ["vector of variables", "vector affine function"], array_dim in [2, 3], T in [Int, Float64]
    dim = 2 * array_dim
    mock = MOIU.MockOptimizer(MILPModel{T}())
    model = COIB.Sort2MILP{T}(mock)

    if T == Int
        @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    end
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Sort,
    )

    if T == Int
        x, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:dim])
    elseif T == Float64
        x = MOI.add_variables(model, dim)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables(x)
    elseif fct_type == "vector affine function"
        MOIU.vectorize(x)
    else
        @assert false
    end

    @test_throws AssertionError MOI.add_constraint(model, fct, CP.Sort(array_dim))

    for i in 1:array_dim
        MOI.add_constraint(model, x[array_dim + i], MOI.LessThan(T(4)))
        MOI.add_constraint(model, x[array_dim + i], MOI.GreaterThan(-T(4)))
    end

    c = MOI.add_constraint(model, fct, CP.Sort(array_dim))

    for i in 1:dim
        @test MOI.is_valid(model, x[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.Sort}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.Sort) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == [(MOI.ZeroOne,)]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
            (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 2 * array_dim ^ 2
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.ZeroOne}()) == array_dim ^ 2
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == 4 * array_dim
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}()) == array_dim ^ 2
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}()) == array_dim ^ 2 + array_dim - 1

        @test Set(MOI.get(bridge, MOI.ListOfVariableIndices())) == Set(vcat(bridge.vars_flow, bridge.vars_unicity))
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.ZeroOne}()) == bridge.vars_unicity_bin
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == vcat(bridge.cons_transportation_x, bridge.cons_transportation_y, bridge.cons_unicity_x, bridge.cons_unicity_y)
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}()) == vec(bridge.cons_flow_lt)
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}()) == vcat(vec(bridge.cons_flow_gt), bridge.cons_sort)
    end

    @testset "Variables" begin
        @test length(bridge.vars_flow) == array_dim ^ 2
        @test length(bridge.vars_unicity) == array_dim ^ 2
        @test length(bridge.vars_unicity_bin) == array_dim ^ 2

        for i in 1:array_dim
            for j in 1:array_dim
                @test MOI.is_valid(model, bridge.vars_flow[i, j])
                @test MOI.is_valid(model, bridge.vars_unicity[i, j])
                @test MOI.is_valid(model, bridge.vars_unicity_bin[i, j])
            end
        end
    end

    @testset "Constraints: transportation" begin
        @test length(bridge.cons_transportation_x) == array_dim
        @test length(bridge.cons_transportation_y) == array_dim

        for i in 1:array_dim
            @test MOI.is_valid(model, bridge.cons_transportation_x[i])
            s = MOI.get(model, MOI.ConstraintSet(), bridge.cons_transportation_x[i])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_transportation_x[i])

            @test s == MOI.EqualTo(zero(T))
            @test length(f.terms) == 1 + array_dim
            @test f.constant == zero(T)

            t1 = f.terms[1]
            @test t1.coefficient === one(T)
            @test t1.variable == x[i]

            for j in 1:array_dim
                t = f.terms[1 + j]
                @test t.coefficient === -one(T)
                @test t.variable == bridge.vars_flow[i, j]
            end
        end

        for j in 1:array_dim
            @test MOI.is_valid(model, bridge.cons_transportation_y[j])
            s = MOI.get(model, MOI.ConstraintSet(), bridge.cons_transportation_y[j])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_transportation_y[j])

            @test s == MOI.EqualTo(zero(T))
            @test length(f.terms) == 1 + array_dim
            @test f.constant == zero(T)

            t1 = f.terms[1]
            @test t1.coefficient === one(T)
            @test t1.variable == x[array_dim + j]

            for i in 1:array_dim
                t = f.terms[1 + i]
                @test t.coefficient === -one(T)
                @test t.variable == bridge.vars_flow[i, j]
            end
        end
    end

    @testset "Constraints: unicity" begin
        @test length(bridge.cons_unicity_x) == array_dim
        @test length(bridge.cons_unicity_y) == array_dim

        for i in 1:array_dim
            @test MOI.is_valid(model, bridge.cons_unicity_x[i])
            s = MOI.get(model, MOI.ConstraintSet(), bridge.cons_unicity_x[i])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_unicity_x[i])

            @test s == MOI.EqualTo(one(T))
            @test length(f.terms) == array_dim
            @test f.constant == zero(T)

            for j in 1:array_dim
                t = f.terms[j]
                @test t.coefficient === one(T)
                @test t.variable == bridge.vars_unicity[i, j]
            end
        end

        for j in 1:array_dim
            @test MOI.is_valid(model, bridge.cons_unicity_y[j])
            s = MOI.get(model, MOI.ConstraintSet(), bridge.cons_unicity_y[j])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_unicity_y[j])

            @test s == MOI.EqualTo(one(T))
            @test length(f.terms) == array_dim
            @test f.constant == zero(T)

            for i in 1:array_dim
                t = f.terms[i]
                @test t.coefficient === one(T)
                @test t.variable == bridge.vars_unicity[i, j]
            end
        end
    end

    @testset "Constraints: unicity of flows" begin
        @test length(bridge.cons_flow_gt) == array_dim ^ 2
        @test length(bridge.cons_flow_lt) == array_dim ^ 2

        for i in 1:array_dim
            for j in 1:array_dim
                @test MOI.is_valid(model, bridge.cons_flow_gt[i, j])
                @test MOI.is_valid(model, bridge.cons_flow_lt[i, j])

                sgt = MOI.get(model, MOI.ConstraintSet(), bridge.cons_flow_gt[i, j])
                fgt = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_flow_gt[i, j])

                slt = MOI.get(model, MOI.ConstraintSet(), bridge.cons_flow_lt[i, j])
                flt = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_flow_lt[i, j])

                @test sgt == MOI.GreaterThan(zero(T))
                @test slt == MOI.LessThan(zero(T))

                @test length(flt.terms) == 2
                @test flt.constant == zero(T)
                @test length(fgt.terms) == 2
                @test fgt.constant == zero(T)

                t1lt = flt.terms[1]
                @test t1lt.coefficient === one(T)
                @test t1lt.variable == bridge.vars_flow[i, j]

                t1gt = fgt.terms[1]
                @test t1gt.coefficient === one(T)
                @test t1gt.variable == bridge.vars_flow[i, j]

                t2lt = flt.terms[2]
                @test t2lt.coefficient === one(T) * 4
                @test t2lt.variable == bridge.vars_unicity[i, j]

                t2gt = fgt.terms[2]
                @test t2gt.coefficient === one(T) * -4
                @test t2gt.variable == bridge.vars_unicity[i, j]
            end
        end
    end

    @testset "Constraint: sorted array" begin
        @test length(bridge.cons_sort) == array_dim - 1

        for i in 1:(array_dim - 1)
            @test MOI.is_valid(model, bridge.cons_sort[i])

            s = MOI.get(model, MOI.ConstraintSet(), bridge.cons_sort[i])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_sort[i])

            @test s == MOI.GreaterThan(zero(T))
            @test length(f.terms) == 2

            t1 = f.terms[1]
            @test t1.coefficient === one(T)
            @test t1.variable == x[i]

            t2 = f.terms[2]
            @test t2.coefficient === -one(T)
            @test t2.variable == x[i + 1]
        end
    end
end


text:
using RDatasets
using Gadfly

plot(data("car", "SLID"), x="wages", color="language",
     Geom.BarGeometry(:dodge))


text:
@testset "math" begin
    x = collect(1:100)
    y = collect(6:9)

    @test ATools.rxy(1,1) ≈ sqrt(2.)
    @test ATools.phixy(1,1) ≈ π/4

    @test ATools.dxyz([1,1,1], [2,2,2]) ≈ sqrt(3.)

    qs = ones(length(data))
    @test isapprox(ATools.wstd(data, qs), std(data), rtol=1e-4)

    m,s = ATools.mean_std(x, 5, 10)
    @test m ≈ mean(y)
    @test s ≈ std(y)

    fxy = ATools.gline2p(0, 10., 1., 2.)
    @test fxy(0) == 10.0
    @test fxy(1) == 2.0
    @test fxy(0.5) == 6.0
end

max_stars_repo_name: jjgomezcadenas/ATools

max_issues_repo_name: jjgomezcadenas/ATools

max_forks_repo_name: jjgomezcadenas/ATools


text:
# extract all sites
const LIQUID_STR = (:liquid,:LIQUID,:L,:l)

"""
    is_liquid(x::Union{Symbol,String})

Returns `true` if the symbol is in `(:liquid,:LIQUID,:L,:l)`.

If a string is passed, it is converted to symbol.
"""
is_liquid(sym::Symbol) = sym in LIQUID_STR
is_liquid(str::String) = is_liquid(Symbol(str))

const VAPOUR_STR = (:vapor,:VAPOR,:VAPOUR,:vapour,:g,:G,:v,:V,:gas,:GAS)

"""
    is_vapour(x::Union{Symbol,String})

Returns `true` if the symbol is in `(:vapor,:VAPOR,:VAPOUR,:vapour,:g,:G,:v,:V,:gas,:GAS)`.

If a string is passed, it is converted to symbol.
"""
is_vapour(sym::Symbol) = sym in VAPOUR_STR
is_vapour(str::String) = is_vapour(Symbol(str))

const SUPERCRITICAL_STR = (:sc,:SC,:supercritical,:SUPERCRITICAL)

"""
    is_supercritical(x::Union{Symbol,String})

Returns `true` if the symbol is in `(:sc,:SC,:supercritical,:SUPERCRITICAL)`.

If a string is passed, it is converted to symbol.
"""
is_supercritical(sym::Symbol) = sym in SUPERCRITICAL_STR
is_supercritical(str::String) = is_vapour(Symbol(str))

"""
    ∑(iterator)

equivalent to `sum(iterator,init=0.0)`.

"""
function ∑(iterator)
    len = Base.IteratorSize(typeof(iterator)) === Base.HasLength()
    hastype =  (Base.IteratorEltype(typeof(iterator)) === Base.HasEltype()) && (eltype(iterator) !== Any)
    local _0
    if hastype
        _0 = zero(eltype(iterator))
    else
        _0 = 0.0
    end
    len && iszero(length(iterator)) && return _0
    !len && return reduce(Base.add_sum,iterator,init=_0)
    return sum(iterator)
end

function ∑(fn,iterator)
    len = Base.IteratorSize(typeof(iterator)) === Base.HasLength()
    hastype =  (Base.IteratorEltype(typeof(iterator)) === Base.HasEltype()) && (eltype(iterator) !== Any)
    local _0
    if hastype
        _0 = zero(eltype(iterator))
    else
        _0 = 0.0
    end
    len && iszero(length(iterator)) && return _0
    !len && return mapreduce(fn,Base.add_sum,iterator,init=_0)
    return sum(fn,iterator)
end

"""
    xlogx(x::Real)
Return `x * log(x)` for `x ≥ 0`, handling ``x = 0`` by taking the downward limit.

copied from LogExpFunctions.jl
"""
function xlogx(x::Real)
    _0 = zero(x)
    iszero(x) && return _0
    ifelse(x >= _0,x*Base.log(max(_0,x)),_0/_0)
end

@inline function nan_num(V,T,z)
    _0 = zero(V+T+first(z))
    _0/_0
end

@inline function negative_vt(V,T)::Bool
    _0 = zero(V+T)
    (T <= _0) | (V <= _0)
end

function complement_index(i,ij)
    i1,i2 = ij
    ifelse(i1 == i,i2,i1)
end