text:
module Checks_Test

using Compat
using Compat.Test
using IntervalArithmetic
using EAGOBranchBound

S = BnBSolver()
B = BnBModel([Interval(1.0,2.0)])

B.LBD = Float64[]
B.first_num = 100
@test EAGOBranchBound.Term_Check(S,B,10) == false

B.first_num = 0
@test EAGOBranchBound.Term_Check(S,B,10) == false

S.max_nodes = -1
B.LBD = Float64[1.0]
@test EAGOBranchBound.Term_Check(S,B,10) == false

S.max_nodes = 100
S.iter_lim = true
S.max_iter = 5
@test EAGOBranchBound.Term_Check(S,B,10) == false
end

ext: jl


text:
# This is a cache of [Val(1), Val(2), ...]
# Hard-coded for now to only handle tensors up to order 100
const ValCache = Val[Val(n) for n in 0:100]
# Faster conversions of collection to tuple than `Tuple(::AbstractVector)`
_NTuple(::Val{N}, v::Vector{T}) where {N,T} = ntuple(n -> v[n], Val(N))
_Tuple(v::Vector{T}) where {T} = _NTuple(ValCache[length(v) + 1], v)
_Tuple(t::Tuple) = t

"""
    ValLength(::Type{NTuple{N}}) = Val{N}
"""
ValLength(::Type{NTuple{N,T}}) where {N,T} = Val{N}

"""
    ValLength(::NTuple{N}) = Val(N)
"""
ValLength(::NTuple{N}) where {N} = Val(N)

# Only to help with backwards compatibility, this
# is not type stable and therefore not efficient.
ValLength(v::Vector) = Val(length(v))

ValLength(::Tuple{Vararg{<:Any,N}}) where {N} = Val(N)

ValLength(::Type{<:Tuple{Vararg{<:Any,N}}}) where {N} = Val{N}

ValLength(::CartesianIndex{N}) where {N} = Val(N)
ValLength(::Type{CartesianIndex{N}}) where {N} = Val{N}

push(s::Tuple, val) = (s..., val)

pushfirst(s::Tuple, val) = (val, s...)

pop(s::NTuple{N}) where {N} = ntuple(i -> s[i], Val(N - 1))

popfirst(s::NTuple{N}) where {N} = ntuple(i -> s[i + 1], Val(N - 1))

# Permute some other type by perm
# (for example, tuple, MVector, etc.)
# as long as the constructor accepts a tuple
@inline function _permute(s, perm)
  return ntuple(i -> s[perm[i]], ValLength(perm))
end

permute(s::Tuple, perm) = _permute(s, perm)

# TODO: is this needed?
function permute(s::T, perm) where {T<:NTuple}
  return T(_permute(s, perm))
end

function permute(s::T, perm) where {T}
  return T(_permute(Tuple(s), perm))
end

# TODO: This is to handle Vector, is this correct?
permute(s::AbstractVector, perm) = _permute(s, perm)

sim(s::NTuple) = s

# type stable findfirst
@inline _findfirst(args...) = (i = findfirst(args...); i === nothing ? 0 : i)

"""
    getperm(col1,col2)

Get the permutation that takes collection 2 to collection 1,
such that col2[p].==col1
"""
@inline function getperm(s1, s2)
  return ntuple(i -> _findfirst(==(@inbounds s1[i]), s2), length(s1))
end

"""
    getperms(col1,col2,col3)

Get the permutations that takes collections 2 and 3 to collection 1.
"""
function getperms(s, s1, s2)
  N = length(s)
  N1 = length(s1)
  N2 = length(s2)
  N1 + N2 ≠ N && error("Size of partial sets don't match with total set")
  perm1 = ntuple(i -> findfirst(==(s1[i]), s), Val(N1))
  perm2 = ntuple(i -> findfirst(==(s2[i]), s), Val(N2))
  isperm((perm1..., perm2...)) ||
    error("Combined permutations are $((perm1...,perm2...)), not a valid permutation")
  return perm1, perm2
end

function invperm!(permres, perm)
  for i in 1:length(perm)
    permres[perm[i]] = i
  end
  return permres
end

function invperm(perm::NTuple{N,Int}) where {N}
  mpermres = MVector{N,Int}(undef)
  invperm!(mpermres, perm)
  return Tuple(mpermres)
end

function invperm(perm)
  permres = similar(perm)
  invperm!(permres, perm)
  return permres
end

# Override TupleTools.isperm to speed up
# Strided.permutedims a bit (see:
# https://github.com/Jutho/Strided.jl/issues/15)
function isperm(p::NTuple{N}) where {N}
  N < 6 && return Base.isperm(p)
  used = @MVector zeros(Bool, N)
  for a in p
    (0 < a <= N) && (used[a] ⊻= true) || return false
  end
  return true
end

"""
    is_trivial_permutation(P)

Determine if P is a trivial permutation.
"""
function is_trivial_permutation(P)
  #isperm(P) || error("Input is not a permutation")
  # TODO: use `all(n->P[n]==n,1:length(P))`?
  N = length(P)
  for n in 1:N
    @inbounds P[n] != n && return false
  end
  return true
end

# Combine a bunch of tuples
@inline flatten(x) = x
@inline flatten(x, y) = (x..., y...)
@inline flatten(x, y, z...) = (x..., flatten(y, z...)...)

function _deleteat(t, pos, i)
  i < pos && return t[i]
  return t[i + 1]
end

function deleteat(t::NTuple{N}, pos::Integer) where {N}
  return ntuple(i -> _deleteat(t, pos, i), Val(N - 1))
end

deleteat(t::Tuple, I::Tuple{Int}) = deleteat(t, I[1])
function deleteat(t::Tuple, I::Tuple{Int,Int,Vararg{Int}})
  return deleteat_sorted(t, sort(I; rev=true))
end

deleteat_sorted(t::Tuple, pos::Int64) = deleteat(t, pos[1])
deleteat_sorted(t::Tuple, pos::Tuple{Int}) = deleteat(t, pos[1])
function deleteat_sorted(t::Tuple, pos::NTuple{N,Int}) where {N}
  return deleteat_sorted(deleteat_sorted(t, pos[1]), Base.tail(pos))
end

# Make a slice of the block on the specified dimensions
# Make this a generic tupletools function (TupleTools.jl calls it getindices)
function getindices(t::Tuple, I::NTuple{N,Int}) where {N}
  return ntuple(i -> t[I[i]], Val(N))
end

# Taken from TupleTools.jl
"""
    sort(t::Tuple; lt=isless, by=identity, rev::Bool=false) -> ::Tuple
Sorts the tuple `t`.
"""
Base.sort(t::Tuple; lt=isless, by=identity, rev::Bool=false) = _sort(t, lt, by, rev)
@inline function _sort(t::Tuple, lt=isless, by=identity, rev::Bool=false)
  t1, t2 = _split(t)
  t1s = _sort(t1, lt, by, rev)
  t2s = _sort(t2, lt, by, rev)
  return _merge(t1s, t2s, lt, by, rev)
end
_sort(t::Tuple{Any}, lt=isless, by=identity, rev::Bool=false) = t
_sort(t::Tuple{}, lt=isless, by=identity, rev::Bool=false) = t

function _split(t::NTuple{N}) where {N}
  M = N >> 1
  return ntuple(i -> t[i], M), ntuple(i -> t[i + M], N - M)
end

function _merge(t1::Tuple, t2::Tuple, lt, by, rev)
  if lt(by(first(t1)), by(first(t2))) != rev
    return (first(t1), _merge(Base.tail(t1), t2, lt, by, rev)...)
  else
    return (first(t2), _merge(t1, Base.tail(t2), lt, by, rev)...)
  end
end
_merge(t1::Tuple{}, t2::Tuple, lt, by, rev) = t2
_merge(t1::Tuple, t2::Tuple{}, lt, by, rev) = t1

#function tail(t::NTuple{N})
#  return ntuple(i -> t[i+1],Val(N-1))
#end

function _insertat(t, pos, n_insert, val, i)
  if i < pos
    return t[i]
  elseif i > pos + n_insert - 1
    return t[i - n_insert + 1]
  end
  return val[i - pos + 1]
end

"""
    insertat

Remove the value at pos and insert the elements in val
"""
function insertat(t::NTuple{N}, val::NTuple{M}, pos::Integer) where {N,M}
  return ntuple(i -> _insertat(t, pos, M, val, i), Val(N + M - 1))
end

function insertat(t::NTuple{N}, val, pos::Integer) where {N}
  return insertat(t, tuple(val), pos)
end

function _insertafter(t, pos, n_insert, val, i)
  if i <= pos
    return t[i]
  elseif i > pos + n_insert
    return t[i - n_insert]
  end
  return val[i - pos]
end

"""
    insertafter(t, val, pos)

Insert the elements in val after the position pos
"""
function insertafter(t::NTuple{N}, val::NTuple{M}, pos::Integer) where {N,M}
  return ntuple(i -> _insertafter(t, pos, M, val, i), Val(N + M))
end

function insertafter(t::NTuple{N}, val, pos::Integer) where {N}
  return insertafter(t, tuple(val), pos)
end

"""
    isdisjoint(s1, s2)

Determine if s1 and s2 have no overlapping elements.
"""
function isdisjoint(s1, s2)
  for i1 in 1:length(s1)
    for i2 in 1:length(s2)
      s1[i1] == s2[i2] && return false
    end
  end
  return true
end

"""
    diff(t::Tuple)

For a tuple of length N, return a tuple of length N-1
where element i is t[i+1] - t[i].
"""
diff(t::NTuple{N}) where {N} = ntuple(i -> t[i + 1] - t[i], Val(N - 1))

function count_unique(labelsT1, labelsT2)
  count = 0
  for l1 in labelsT1
    l1 ∉ labelsT2 && (count += 1)
  end
  return count
end

function count_common(labelsT1, labelsT2)
  count = 0
  for l1 in labelsT1
    l1 ∈ labelsT2 && (count += 1)
  end
  return count
end

function intersect_positions(labelsT1, labelsT2)
  for i1 in 1:length(labelsT1)
    for i2 in 1:length(labelsT2)
      if labelsT1[i1] == labelsT2[i2]
        return i1, i2
      end
    end
  end
  return nothing
end

function is_replacement(labelsT1, labelsT2)
  return count_unique(labelsT1, labelsT2) == 1 && count_common(labelsT1, labelsT2) == 1
end

function is_combiner(labelsT1, labelsT2)
  return count_unique(labelsT1, labelsT2) == 1 && count_common(labelsT1, labelsT2) > 1
end

function is_uncombiner(labelsT1, labelsT2)
  return count_unique(labelsT1, labelsT2) > 1 && count_common(labelsT1, labelsT2) == 1
end


text:
using DelimitedFiles
using Printf
using Combinatorics

# Import Data
nums = readdlm("./input_a.dat")

# Function to check all possible sums of `arr` in groups of `n` that equal `target`
function check_sums(arr, n, target)
    pairs = collect(combinations(arr, n))
    for p in pairs
        if sum(p) == target
            @printf("Combination of %d values that results in %d\n", n, target)
            println(p)
            @printf("Sum: %i\tProduct: %i\n\n", sum(p), prod(p))
            return p
        end
    end
end

check_sums(nums, 2, 2020) # Check 2 values that sum to 2020
check_sums(nums, 3, 2020) # Check 3 values that sum to 2020


text:
function LeastSquaresBuilder()
    return LeastSquaresBuilder(())
end

function build(obj::LeastSquaresBuilder)
    return jcall(obj, "build", LeastSquaresProblem, ())
end

function checker(obj::LeastSquaresBuilder, arg0::ConvergenceChecker)
    return jcall(obj, "checker", LeastSquaresBuilder, (ConvergenceChecker,), arg0)
end

function checker_pair(obj::LeastSquaresBuilder, arg0::ConvergenceChecker)
    return jcall(obj, "checkerPair", LeastSquaresBuilder, (ConvergenceChecker,), arg0)
end

function lazy_evaluation(obj::LeastSquaresBuilder, arg0::jboolean)
    return jcall(obj, "lazyEvaluation", LeastSquaresBuilder, (jboolean,), arg0)
end

function max_evaluations(obj::LeastSquaresBuilder, arg0::jint)
    return jcall(obj, "maxEvaluations", LeastSquaresBuilder, (jint,), arg0)
end

function max_iterations(obj::LeastSquaresBuilder, arg0::jint)
    return jcall(obj, "maxIterations", LeastSquaresBuilder, (jint,), arg0)
end

function model(obj::LeastSquaresBuilder, arg0::MultivariateJacobianFunction)
    return jcall(obj, "model", LeastSquaresBuilder, (MultivariateJacobianFunction,), arg0)
end

function model(obj::LeastSquaresBuilder, arg0::MultivariateVectorFunction, arg1::MultivariateMatrixFunction)
    return jcall(obj, "model", LeastSquaresBuilder, (MultivariateVectorFunction, MultivariateMatrixFunction), arg0, arg1)
end

function parameter_validator(obj::LeastSquaresBuilder, arg0::ParameterValidator)
    return jcall(obj, "parameterValidator", LeastSquaresBuilder, (ParameterValidator,), arg0)
end

function start(obj::LeastSquaresBuilder, arg0::RealVector)
    return jcall(obj, "start", LeastSquaresBuilder, (RealVector,), arg0)
end

function start(obj::LeastSquaresBuilder, arg0::Vector{jdouble})
    return jcall(obj, "start", LeastSquaresBuilder, (Vector{jdouble},), arg0)
end

function target(obj::LeastSquaresBuilder, arg0::RealVector)
    return jcall(obj, "target", LeastSquaresBuilder, (RealVector,), arg0)
end

function target(obj::LeastSquaresBuilder, arg0::Vector{jdouble})
    return jcall(obj, "target", LeastSquaresBuilder, (Vector{jdouble},), arg0)
end

function weight(obj::LeastSquaresBuilder, arg0::RealMatrix)
    return jcall(obj, "weight", LeastSquaresBuilder, (RealMatrix,), arg0)
end


text:
using oneAPI, Test

function vadd(a, b, c)
    i = get_global_id()
    @inbounds c[i] = a[i] + b[i]
    return
end

dims = (2,)
a = round.(rand(Float32, dims) * 100)
b = round.(rand(Float32, dims) * 100)
c = similar(a)

d_a = oneArray(a)
d_b = oneArray(b)
d_c = oneArray(c)

len = prod(dims)
@oneapi items=len vadd(d_a, d_b, d_c)
c = Array(d_c)
@test a+b ≈ c


text:
include("exp_mnist.jl")

include("exp_cifar.jl")

# julia --project main.jl

max_stars_repo_name: lihebi/AdvAE

max_issues_repo_name: lihebi/AdvAE

max_forks_repo_name: lihebi/AdvAE


text:
include("partition_algorithms.jl")
include("kd_tree.jl")
include("partition_sample.jl")


text:
include("../utils.jl")
include("group_utils.jl")

@testset "Translation group" begin
    @testset "real" begin
        G = TranslationGroup(2, 3)
        @test repr(G) == "TranslationGroup(2, 3; field = ℝ)"
        @test repr(TranslationGroup(2, 3; field=ℂ)) == "TranslationGroup(2, 3; field = ℂ)"

        @test (@inferred invariant_metric_dispatch(G, LeftAction())) === Val(true)
        @test (@inferred invariant_metric_dispatch(G, RightAction())) === Val(true)
        @test (@inferred Manifolds.biinvariant_metric_dispatch(G)) === Val(true)
        @test is_default_metric(MetricManifold(G, EuclideanMetric())) === true
        @test Manifolds.default_metric_dispatch(MetricManifold(G, EuclideanMetric())) ===
              Val{true}()
        types = [Matrix{Float64}]
        @test base_manifold(G) === Euclidean(2, 3)

        pts = [reshape(i:(i + 5), (2, 3)) for i in 1:3]
        vpts = [reshape(-2:3, (2, 3))]
        for T in types
            gpts = convert.(T, pts)
            vgpts = convert.(T, vpts)
            @test compose(G, gpts[1], gpts[2]) ≈ gpts[1] + gpts[2]
            test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true)
        end
    end

    @testset "complex" begin
        G = TranslationGroup(2, 3; field=ℂ)
        @test repr(G) == "TranslationGroup(2, 3; field = ℂ)"

        types = [Matrix{ComplexF64}]
        @test base_manifold(G) === Euclidean(2, 3; field=ℂ)

        pts = [reshape(complex.(i:(i + 5), (i + 1):(i + 6)), (2, 3)) for i in 1:3]
        vpts = [reshape(complex.(-2:3, -1:4), (2, 3))]
        for T in types
            gpts = convert.(T, pts)
            vgpts = convert.(T, vpts)
            @test compose(G, gpts[1], gpts[2]) ≈ gpts[1] + gpts[2]
            @test translate_diff(G, gpts[2], gpts[1], vgpts[1]) ≈ vgpts[1]
            test_group(G, gpts, vgpts, vgpts; test_diff=true, test_invariance=true)
        end
    end
end


text:
using LinearAlgebra

mutable struct Broyden{R<:Real,C<:Union{R,Complex{R}},T<:AbstractArray{C}}
    H
    theta_bar::R
    function Broyden{R,C,T}(x::T, H, theta_bar) where {R,C,T}
        new(H, theta_bar)
    end
end

Broyden(
    x::T;
    H = I,
    theta_bar = R(0.2),
) where {R<:Real,C<:Union{R,Complex{R}},T<:AbstractArray{C}} =
    Broyden{R,C,T}(x, H, theta_bar)

_sign(x::R) where {R} = x == 0 ? R(1) : sign(x)

function update!(L::Broyden{R,C,T}, s, y) where {R,C,T}
    Hy = L.H * y
    sH = s' * L.H
    delta = dot(Hy, s) / norm(s)^2
    theta = if abs(delta) >= L.theta_bar
        R(1)
    else
        (1 - _sign(delta) * L.theta_bar) / (1 - delta)
    end
    L.H += (s - Hy) / dot(s, (1 / theta - 1) * s + Hy) * sH
end

import Base: *

function (*)(L::Broyden, v)
    w = similar(v)
    mul!(w, L, v)
end

import LinearAlgebra: mul!

mul!(d::T, L::Broyden, v::T) where {T} = mul!(d, L.H, v)


text:
# Tests for distances
function t_lr()

N = 10	# variables
M = 25  # samples

vectordata = rand(M)
matrixdata = rand(N,M)

labels = 3*vectordata+0.1*randn(M)
multilabels = [collect(1:10) 3*collect(1:10)]'*matrixdata+0.1*randn(2,M)

A01 = j4pr.datacell(vectordata)
A02 = j4pr.datacell(vectordata, labels)
A03 = j4pr.datacell(vectordata, multilabels)
A04 = j4pr.datacell(matrixdata)
A05 = j4pr.datacell(matrixdata, labels)
A06 = j4pr.datacell(matrixdata, multilabels)

Test.@test try
	for A in [A01, A02, A03, A04, A05, A06]
		w = A |> j4pr.lr(1)
		w2 = A |> j4pr.lr(1,bias=false)
		A |> w
		A |> w2
	end
	true
catch
	false
end

end


text:
using Test
using Printf
using Random
using Statistics
using LinearAlgebra
using Logging

using JLD2
using FFTW
using OffsetArrays

using Oceananigans.Architectures: @hascuda
@hascuda begin
    import CUDAdrv
    using CuArrays
    using CUDAnative
end

using Oceananigans
using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Operators
using Oceananigans.BoundaryConditions
using Oceananigans.Fields
using Oceananigans.Coriolis
using Oceananigans.Buoyancy
using Oceananigans.Forcing
using Oceananigans.Solvers
using Oceananigans.Models
using Oceananigans.Simulations
using Oceananigans.Diagnostics
using Oceananigans.OutputWriters
using Oceananigans.TurbulenceClosures
using Oceananigans.AbstractOperations
using Oceananigans.Logger
using Oceananigans.Utils

using Statistics: mean
using LinearAlgebra: norm
using GPUifyLoops: @launch, @loop
using NCDatasets: Dataset

import Oceananigans.Fields: interior
import Oceananigans.Utils: datatuple

using Oceananigans.Diagnostics: run_diagnostic, velocity_div!
using Oceananigans.TimeSteppers: _compute_w_from_continuity!
using Oceananigans.AbstractOperations: Computation, compute!

#####
##### On CI servers select the GPU with the most available memory or with the
##### highest capability if testing needs to be thorough).
##### Source credit: https://github.com/JuliaGPU/CuArrays.jl/pull/526
#####

@hascuda begin
    gpu_candidates = [(dev=dev, cap=CUDAdrv.capability(dev),
                       mem=CUDAdrv.CuContext(ctx->CUDAdrv.available_memory(), dev))
                       for dev in CUDAdrv.devices()]

    thorough = parse(Bool, get(ENV, "CI_THOROUGH", "false"))
    if thorough
        sort!(gpu_candidates, by=x->(x.cap, x.mem))
    else
        sort!(gpu_candidates, by=x->x.mem)
    end

    pick = last(gpu_candidates)
    device!(pick.dev)
end

#####
##### Useful utilities
#####

function get_model_field(field_name, model)
    if field_name ∈ (:u, :v, :w)
        return getfield(model.velocities, field_name)
    else
        return getfield(model.tracers, field_name)
    end
end

datatuple(A) = NamedTuple{propertynames(A)}(Array(data(a)) for a in A)

function get_output_tuple(output, iter, tuplename)
    file = jldopen(output.filepath, "r")
    output_tuple = file["timeseries/$tuplename/$iter"]
    close(file)
    return output_tuple
end

#####
##### Testing parameters
#####

float_types = (Float32, Float64)

         archs = (CPU(),)
@hascuda archs = (CPU(), GPU())

closures = (
    :ConstantIsotropicDiffusivity,
    :ConstantAnisotropicDiffusivity,
    :AnisotropicBiharmonicDiffusivity,
    :TwoDimensionalLeith,
    :SmagorinskyLilly,
    :BlasiusSmagorinsky,
    :RozemaAnisotropicMinimumDissipation,
    :VerstappenAnisotropicMinimumDissipation
)

#####
##### Run tests!
#####

with_logger(ModelLogger()) do
    @testset "Oceananigans" begin
        include("test_grids.jl")
        include("test_operators.jl")
        include("test_boundary_conditions.jl")
        include("test_fields.jl")
        include("test_halo_regions.jl")
        include("test_solvers.jl")
        include("test_pressure_solvers.jl")
        include("test_coriolis.jl")
        include("test_buoyancy.jl")
        include("test_surface_waves.jl")
        include("test_models.jl")
        include("test_simulations.jl")
        include("test_time_stepping.jl")
        include("test_time_stepping_bcs.jl")
        include("test_forcings.jl")
        include("test_turbulence_closures.jl")
        include("test_dynamics.jl")
        include("test_diagnostics.jl")
        include("test_output_writers.jl")
        include("test_abstract_operations.jl")
        include("test_regression.jl")
        include("test_examples.jl")
        include("test_verification.jl")
    end
end


text:
#### thermodynamics

export new_thermo_state, recover_thermo_state

"""
    new_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars)

Create a new thermodynamic state, based on the `state`, and _not_
the `aux` state.

!!! note
    This method calls the iterative saturation adjustment
    procedure for EquilMoist models.
"""
function new_thermo_state end

# First dispatch on compressibility
new_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars) =
    new_thermo_state(compressibility_model(atmos), atmos, state, aux)

new_thermo_state(::Compressible, atmos::AtmosModel, state::Vars, aux::Vars) =
    new_thermo_state(
        atmos,
        energy_model(atmos),
        moisture_model(atmos),
        state,
        aux,
    )
new_thermo_state(::Anelastic1D, atmos::AtmosModel, state::Vars, aux::Vars) =
    new_thermo_state_anelastic(atmos, state, aux)

"""
    recover_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars)

An atmospheric thermodynamic state.

!!! warn
    For now, we are directly calling new_thermo_state to avoid
    inconsistent aux states in kernels where the aux states are
    out of sync with the boundary state.

# TODO: Define/call `recover_thermo_state` when it's safely implemented
  (see https://github.com/CliMA/ClimateMachine.jl/issues/1648)
"""
function recover_thermo_state end

# First dispatch on compressibility
recover_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars) =
    new_thermo_state(compressibility_model(atmos), atmos, state, aux)

recover_thermo_state(
    ::Compressible,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
) = new_thermo_state(
    atmos,
    energy_model(atmos),
    moisture_model(atmos),
    state,
    aux,
)

recover_thermo_state(::Anelastic1D, atmos::AtmosModel, state::Vars, aux::Vars) =
    new_thermo_state_anelastic(atmos, state, aux)

function new_thermo_state(
    atmos::AtmosModel,
    energy::TotalEnergyModel,
    moist::DryModel,
    state::Vars,
    aux::Vars,
)
    param_set = parameter_set(atmos)
    e_int = internal_energy(atmos, state, aux)
    return PhaseDry(param_set, e_int, state.ρ)
end

function new_thermo_state(
    atmos::AtmosModel,
    energy::TotalEnergyModel,
    moist::EquilMoist,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    param_set = parameter_set(atmos)
    return PhaseEquil(
        param_set,
        e_int,
        state.ρ,
        state.moisture.ρq_tot / state.ρ,
        moist.maxiter,
        moist.tolerance,
    )
end

function new_thermo_state(
    atmos::AtmosModel,
    energy::TotalEnergyModel,
    moist::NonEquilMoist,
    state::Vars,
    aux::Vars,
)
    param_set = parameter_set(atmos)
    e_int = internal_energy(atmos, state, aux)
    q = PhasePartition(
        state.moisture.ρq_tot / state.ρ,
        state.moisture.ρq_liq / state.ρ,
        state.moisture.ρq_ice / state.ρ,
    )

    return PhaseNonEquil{eltype(state), typeof(param_set)}(
        param_set,
        e_int,
        state.ρ,
        q,
    )
end

function new_thermo_state(
    atmos::AtmosModel,
    energy::θModel,
    moist::DryModel,
    state::Vars,
    aux::Vars,
)
    param_set = parameter_set(atmos)
    θ_liq_ice = state.energy.ρθ_liq_ice / state.ρ
    return PhaseDry_ρθ(param_set, state.ρ, θ_liq_ice)
end

function new_thermo_state(
    atmos::AtmosModel,
    energy::θModel,
    moist::EquilMoist,
    state::Vars,
    aux::Vars,
)
    param_set = parameter_set(atmos)
    θ_liq_ice = state.energy.ρθ_liq_ice / state.ρ
    return PhaseEquil_ρθq(
        param_set,
        state.ρ,
        θ_liq_ice,
        state.moisture.ρq_tot / state.ρ,
        moist.maxiter,
        moist.tolerance,
    )
end

function new_thermo_state(
    atmos::AtmosModel,
    energy::θModel,
    moist::NonEquilMoist,
    state::Vars,
    aux::Vars,
)
    param_set = parameter_set(atmos)
    θ_liq_ice = state.energy.ρθ_liq_ice / state.ρ
    q = PhasePartition(
        state.moisture.ρq_tot / state.ρ,
        state.moisture.ρq_liq / state.ρ,
        state.moisture.ρq_ice / state.ρ,
    )

    return PhaseNonEquil{eltype(state), typeof(param_set)}(
        param_set,
        state.ρ,
        θ_liq_ice,
        q,
    )
end


text:
@testset "Core" begin
    @test sprint(showerror, Stratum.StratumError(-32700, "FOO", "BAR")) == "ParseError: FOO (BAR)"
    @test sprint(showerror, Stratum.StratumError(-32600, "FOO", "BAR")) == "InvalidRequest: FOO (BAR)"
    @test sprint(showerror, Stratum.StratumError(-32601, "FOO", "BAR")) == "MethodNotFound: FOO (BAR)"
    @test sprint(showerror, Stratum.StratumError(-32602, "FOO", "BAR")) == "InvalidParams: FOO (BAR)"
    @test sprint(showerror, Stratum.StratumError(-32603, "FOO", "BAR")) == "InternalError: FOO (BAR)"
    @test sprint(showerror, Stratum.StratumError(-32099, "FOO", "BAR")) == "serverErrorStart: FOO (BAR)"
    @test sprint(showerror, Stratum.StratumError(-32000, "FOO", "BAR")) == "serverErrorEnd: FOO (BAR)"
    @test sprint(showerror, Stratum.StratumError(-32002, "FOO", "BAR")) == "ServerNotInitialized: FOO (BAR)"
    @test sprint(showerror, Stratum.StratumError(-32001, "FOO", "BAR")) == "UnknownErrorCode: FOO (BAR)"
    @test sprint(showerror, Stratum.StratumError(-32800, "FOO", "BAR")) == "RequestCancelled: FOO (BAR)"
    @test sprint(showerror, Stratum.StratumError(-32801, "FOO", "BAR")) == "ContentModified: FOO (BAR)"
    @test sprint(showerror, Stratum.StratumError(1, "FOO", "BAR")) == "Unkonwn: FOO (BAR)"
end


text:
module WaveFD

using Base.Threads, CvxCompress, DSP, Distributed, DistributedArrays, FFTW, LinearAlgebra, NearestNeighbors, Random, SpecialFunctions, StaticArrays, Statistics, WaveFD_jll

import
Base.convert,
Base.copy,
Base.get,
Base.min,
Base.max,
Base.maximum,
Base.show,
Base.size

abstract type Language end
struct LangC <: Language end
struct LangJulia <: Language end
show(io::IO, l::LangC) = write(io, "C")
show(io::IO, l::LangJulia) = write(io, "Julia")

abstract type ImagingCondition end
struct ImagingConditionStandard <: ImagingCondition end
struct ImagingConditionWaveFieldSeparationFWI <: ImagingCondition end
struct ImagingConditionWaveFieldSeparationRTM <: ImagingCondition end

include("stencil.jl")
include("spacetime.jl")
include("compressedio.jl")
include("wavelet.jl")
include("absorb.jl")
include("illumination.jl")

include("prop2DAcoIsoDenQ_DEO2_FDTD.jl")
include("prop2DAcoTTIDenQ_DEO2_FDTD.jl")
include("prop2DAcoVTIDenQ_DEO2_FDTD.jl")

include("prop3DAcoIsoDenQ_DEO2_FDTD.jl")
include("prop3DAcoTTIDenQ_DEO2_FDTD.jl")
include("prop3DAcoVTIDenQ_DEO2_FDTD.jl")

export
dtmod,
Prop2DAcoIsoDenQ_DEO2_FDTD,
Prop2DAcoTTIDenQ_DEO2_FDTD,
Prop2DAcoVTIDenQ_DEO2_FDTD,
Prop3DAcoIsoDenQ_DEO2_FDTD,
Prop3DAcoTTIDenQ_DEO2_FDTD,
Prop3DAcoVTIDenQ_DEO2_FDTD,
fieldfile!,
free,
Ginsu,
linearadjoint,
linearforward,
nonlinearforward,
ntmod,
reportfreq!,
sourceillum!,
sourceillum,
traces,
Wavelet,
WaveletDerivRicker,
WaveletMinPhaseRicker,
WaveletCausalRicker,
WaveletOrmsby,
WaveletMinPhaseOrmsby,
WaveletRicker,
WaveletSine

end

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
#!/usr/bin/env julia

using Luxor, Colors

include("julia-logo.jl") # the julia logo coordinates

function draw_logo_clip(x, y)
    foregroundcolors = diverging_palette(rand(0:360), rand(0:360), 200, s = 0.99, b=0.8)
    gsave()
    translate(x-100, y)
    julialogo(false, true) # use julia logo as clipping mask
    clip()
    for i in 1:500
        sethue(foregroundcolors[rand(1:end)])
        circle(rand(-50:350), rand(0:300), 15, :fill)
    end
    clipreset()
    grestore()
end

function clipping_test(fname)
    Drawing(500, 500, fname)
    origin()
    background("white")
    setopacity(.4)
    foregroundcolors = diverging_palette(rand(0:360), rand(0:360), 200, s = 0.99, b=0.8)
    draw_logo_clip(0, 0)
    finish()
    println("finished test: output in $(fname)")
end

clipping_test("/tmp/clipping-tests.png")


text:
using NNlib, Test, Statistics, Random
using ChainRulesCore, ChainRulesTestUtils
using Base.Broadcast: broadcasted
import FiniteDifferences
using FiniteDifferences: FiniteDifferenceMethod, central_fdm
import Zygote
using Zygote: gradient
using StableRNGs

const rng = StableRNG(123)

include("test_utils.jl")

@testset "Activation Functions" begin
    include("activations.jl")
end

@testset "Batched Multiplication" begin
    include("batchedmul.jl")
end

@testset "Convolution" begin
    include("conv.jl")
    include("conv_bias_act.jl")
end

@testset "Inference" begin
    include("inference.jl")
end

@testset "Pooling" begin
    include("pooling.jl")
end

@testset "Softmax" begin
    include("softmax.jl")
end

@testset "Misc Stuff" begin
    include("misc.jl")
end


text:
"Concatenate posterior draws for different variables."
function vcat(pds::PosteriorDraws...)
    isempty(pds) && error("Need at least one argument.")
    PosteriorDraws(pds[1].len,
                   vcat([pd.keys for pd in pds]...),
                   vcat([pd.vars for pd in pds]...))
end

"Add variables to `pd`, names by the given keys in `key => var` pairs."
function addvars{T}(pd::PosteriorDraws, key_var_pairs::Pair{Symbol, T}...)
    vcat(pd, PosteriorDraws(key_var_pairs...))
end

"Return `pd` without the variables names by `keys`."
function dropvars(pd::PosteriorDraws, keys...)
    keep = setdiff(1:size(pd, 1), _key2index(pd, [keys...]))
    PosteriorDraws(pd.len, pd.keys[keep], pd.vars[keep])
end

"""
Map variables for `keys` by `f`, add the result with `result_key` to `pd`.
"""
function map(pd::PosteriorDraws, result_key::Symbol, f, keys::Symbol...)
    result_key ∈ pd.keys &&
        error(ArgumentError("Can't overwrite an existing variable."))
    vars = pd.vars[_key2index(pd, [keys...])]
    addvars(pd, result_key => map(f, vars...))
end

"""
Transform expression by replacing each `@v(symbol)` with a gensym, and
collecting `symbol => gensym` in `captured_names` as a side effect.
"""
_transform_expr(expr, captured_names::Dict{Symbol, Symbol}) = expr

function _transform_expr(expr::Expr, captured_names::Dict{Symbol, Symbol})
    @match expr begin
        @v(var_) => get!(()->gensym(var), captured_names, var)
        @v(var__) => error("@v takes only a single variable name.")
        :e_ => expr             # pass through quoted
        e_ => Expr(e.head, [_transform_expr(arg, captured_names)
                            for arg in e.args]...)
    end
end

"""
When called without the second argument, return a tuple with the
captured names as the second element.
"""
function _transform_expr(expr)
    captured_names = Dict{Symbol,Symbol}()
    (_transform_expr(expr, captured_names), captured_names)
end

"""
Helper function that creates the body of the @pdmap macro.
"""
function _map_helper(key_and_form_pair)
    (result_key, form) = @match key_and_form_pair begin
        (key_ => form_) => (key, form)
        _ => error("expecting key => form")
    end
    (transformed_form, captured_names) = _transform_expr(form)
    arguments = :($([gensym_ for (_, gensym_) in captured_names]...),)
    closure = :($arguments -> $(transformed_form))
    keys = [Meta.quot(varname) for (varname, _) in captured_names]
    (Meta.quot(result_key), esc(closure), keys...)
end

"""
Wrapper macro for `map(pd::PosteriorDraws, ...)`. Uses the syntax
```julia
@pdmap posteriordraws newvar => form(@v(somevar), @v(someothervar), ...)
```
where `newvar` will contain the result, while the `@v` forms are
parsed by a codewalker which generates the appropriate syntax for map.
"""
macro pdmap(pd, key_and_form_pair)
    :(map($(esc(pd)), $(_map_helper(key_and_form_pair)...)))
end


text:
function TestInputGrid(OrientFace)
filename='kiteMeshOsw.rtf';
filename='kite_mesh_small.rtf';
Grid=InputGrid(filename,OrientFace);
end


text:
#! /usr/bin/julia

# Rosetta Code, Primality by trial division

function isprime_trialdivision{T<:Integer}(n::T)
    1 < n || return false
    n != 2 || return true
    isodd(n) || return false
    for i in 3:isqrt(n)
        n%i != 0 || return false
    end
    return true
end

n = 100
a = filter(isprime_trialdivision, [1:n])

if all(a .== primes(n))
    println("The primes <= ", n, " are:\n    ", a)
else
    println("The function does not accurately calculate primes.")
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
points = map(str -> [parse(Int64,str[l]) for l in findall(r"[[:digit:]]+", str)], readlines("inputs/05.txt"))

grid = zeros(Int64, 1000,1000);

smart_range(x1, x2) = x1:(x2 > x1 ? 1 : -1):x2

function smart_2d_range(x1, y1, x2, y2)
  x,y = smart_range(x1, x2), smart_range(y1, y2)
  if x1 == x2
    x = Iterators.repeated(x.start, length(y))
  elseif y1 == y2
    y = Iterators.repeated(y.start, length(x))
  end
  return zip(x, y)
end

function day5p2(points)
  for p in points
    for (x,y) in smart_2d_range(p...)
      grid[x,y] += 1
    end
  end

  count(i -> i > 1, grid)
end

using BenchmarkTools
@btime day5p2(points)


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ActiveRenewableControllerAB <: ActivePowerControl
        bus_control::Int
        from_branch_control::Int
        to_branch_control::Int
        branch_id_control::String
        Freq_Flag::Int
        K_pg::Float64
        K_ig::Float64
        T_p::Float64
        fdbd_pnts::Tuple{Float64, Float64}
        fe_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        P_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        T_g::Float64
        D_dn::Float64
        D_up::Float64
        dP_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        P_lim_inner::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        T_pord::Float64
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters of Active Power Controller including REPCA1 and REECB1

# Arguments
- `bus_control::Int`: Bus Number for voltage control; , validation range: `(0, nothing)`
- `from_branch_control::Int`: Monitored branch FROM bus number for line drop compensation (if 0 generator power will be used), validation range: `(0, nothing)`
- `to_branch_control::Int`: Monitored branch TO bus number for line drop compensation (if 0 generator power will be used), validation range: `(0, nothing)`
- `branch_id_control::String`: Branch circuit id for line drop compensation (as a string). If 0 generator power will be used
- `Freq_Flag::Int`: Frequency Flag for REPCA1: 0: disable, 1:enable, validation range: `(0, 1)`
- `K_pg::Float64`: Active power PI control proportional gain, validation range: `(0, nothing)`
- `K_ig::Float64`: Active power PI control integral gain, validation range: `(0, nothing)`
- `T_p::Float64`: Real power measurement filter time constant (s), validation range: `(0, nothing)`
- `fdbd_pnts::Tuple{Float64, Float64}`: Frequency error dead band thresholds `(fdbd1, fdbd2)`
- `fe_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Upper/Lower limit on frequency error `(fe_min, fe_max)`
- `P_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Upper/Lower limit on power reference `(P_min, P_max)`
- `T_g::Float64`: Power Controller lag time constant, validation range: `(0, nothing)`
- `D_dn::Float64`: Droop for over-frequency conditions, validation range: `(nothing, 0)`
- `D_up::Float64`: Droop for under-frequency conditions, validation range: `(0, nothing)`
- `dP_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Upper/Lower limit on power reference ramp rates`(dP_min, dP_max)`
- `P_lim_inner::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Upper/Lower limit on power reference for REECB`(P_min_inner, P_max_inner)`
- `T_pord::Float64`: Power filter time constant REECB time constant, validation range: `(0, nothing)`
- `P_ref::Float64`: Reference Power Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the ActiveRenewableControllerAB model depends on the Flag
- `n_states::Int`: The states of the ActiveRenewableControllerAB model depends on the Flag
"""
mutable struct ActiveRenewableControllerAB <: ActivePowerControl
    "Bus Number for voltage control; "
    bus_control::Int
    "Monitored branch FROM bus number for line drop compensation (if 0 generator power will be used)"
    from_branch_control::Int
    "Monitored branch TO bus number for line drop compensation (if 0 generator power will be used)"
    to_branch_control::Int
    "Branch circuit id for line drop compensation (as a string). If 0 generator power will be used"
    branch_id_control::String
    "Frequency Flag for REPCA1: 0: disable, 1:enable"
    Freq_Flag::Int
    "Active power PI control proportional gain"
    K_pg::Float64
    "Active power PI control integral gain"
    K_ig::Float64
    "Real power measurement filter time constant (s)"
    T_p::Float64
    "Frequency error dead band thresholds `(fdbd1, fdbd2)`"
    fdbd_pnts::Tuple{Float64, Float64}
    "Upper/Lower limit on frequency error `(fe_min, fe_max)`"
    fe_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Upper/Lower limit on power reference `(P_min, P_max)`"
    P_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Power Controller lag time constant"
    T_g::Float64
    "Droop for over-frequency conditions"
    D_dn::Float64
    "Droop for under-frequency conditions"
    D_up::Float64
    "Upper/Lower limit on power reference ramp rates`(dP_min, dP_max)`"
    dP_lim::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Upper/Lower limit on power reference for REECB`(P_min_inner, P_max_inner)`"
    P_lim_inner::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Power filter time constant REECB time constant"
    T_pord::Float64
    "Reference Power Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    "The states of the ActiveRenewableControllerAB model depends on the Flag"
    states::Vector{Symbol}
    "The states of the ActiveRenewableControllerAB model depends on the Flag"
    n_states::Int
end

function ActiveRenewableControllerAB(bus_control, from_branch_control, to_branch_control, branch_id_control, Freq_Flag, K_pg, K_ig, T_p, fdbd_pnts, fe_lim, P_lim, T_g, D_dn, D_up, dP_lim, P_lim_inner, T_pord, P_ref=1.0, ext=Dict{String, Any}(), )
    ActiveRenewableControllerAB(bus_control, from_branch_control, to_branch_control, branch_id_control, Freq_Flag, K_pg, K_ig, T_p, fdbd_pnts, fe_lim, P_lim, T_g, D_dn, D_up, dP_lim, P_lim_inner, T_pord, P_ref, ext, PowerSystems.get_activeRETypeAB_states(Freq_Flag)[1], PowerSystems.get_activeRETypeAB_states(Freq_Flag)[2], )
end

function ActiveRenewableControllerAB(; bus_control, from_branch_control, to_branch_control, branch_id_control, Freq_Flag, K_pg, K_ig, T_p, fdbd_pnts, fe_lim, P_lim, T_g, D_dn, D_up, dP_lim, P_lim_inner, T_pord, P_ref=1.0, ext=Dict{String, Any}(), states=PowerSystems.get_activeRETypeAB_states(Freq_Flag)[1], n_states=PowerSystems.get_activeRETypeAB_states(Freq_Flag)[2], )
    ActiveRenewableControllerAB(bus_control, from_branch_control, to_branch_control, branch_id_control, Freq_Flag, K_pg, K_ig, T_p, fdbd_pnts, fe_lim, P_lim, T_g, D_dn, D_up, dP_lim, P_lim_inner, T_pord, P_ref, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function ActiveRenewableControllerAB(::Nothing)
    ActiveRenewableControllerAB(;
        bus_control=0,
        from_branch_control=0,
        to_branch_control=0,
        branch_id_control="0",
        Freq_Flag=0,
        K_pg=0,
        K_ig=0,
        T_p=0,
        fdbd_pnts=(0.0, 0.0),
        fe_lim=(min=0.0, max=0.0),
        P_lim=(min=0.0, max=0.0),
        T_g=0,
        D_dn=0,
        D_up=0,
        dP_lim=(min=0.0, max=0.0),
        P_lim_inner=(min=0.0, max=0.0),
        T_pord=0,
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`ActiveRenewableControllerAB`](@ref) `bus_control`."""
get_bus_control(value::ActiveRenewableControllerAB) = value.bus_control
"""Get [`ActiveRenewableControllerAB`](@ref) `from_branch_control`."""
get_from_branch_control(value::ActiveRenewableControllerAB) = value.from_branch_control
"""Get [`ActiveRenewableControllerAB`](@ref) `to_branch_control`."""
get_to_branch_control(value::ActiveRenewableControllerAB) = value.to_branch_control
"""Get [`ActiveRenewableControllerAB`](@ref) `branch_id_control`."""
get_branch_id_control(value::ActiveRenewableControllerAB) = value.branch_id_control
"""Get [`ActiveRenewableControllerAB`](@ref) `Freq_Flag`."""
get_Freq_Flag(value::ActiveRenewableControllerAB) = value.Freq_Flag
"""Get [`ActiveRenewableControllerAB`](@ref) `K_pg`."""
get_K_pg(value::ActiveRenewableControllerAB) = value.K_pg
"""Get [`ActiveRenewableControllerAB`](@ref) `K_ig`."""
get_K_ig(value::ActiveRenewableControllerAB) = value.K_ig
"""Get [`ActiveRenewableControllerAB`](@ref) `T_p`."""
get_T_p(value::ActiveRenewableControllerAB) = value.T_p
"""Get [`ActiveRenewableControllerAB`](@ref) `fdbd_pnts`."""
get_fdbd_pnts(value::ActiveRenewableControllerAB) = value.fdbd_pnts
"""Get [`ActiveRenewableControllerAB`](@ref) `fe_lim`."""
get_fe_lim(value::ActiveRenewableControllerAB) = value.fe_lim
"""Get [`ActiveRenewableControllerAB`](@ref) `P_lim`."""
get_P_lim(value::ActiveRenewableControllerAB) = value.P_lim
"""Get [`ActiveRenewableControllerAB`](@ref) `T_g`."""
get_T_g(value::ActiveRenewableControllerAB) = value.T_g
"""Get [`ActiveRenewableControllerAB`](@ref) `D_dn`."""
get_D_dn(value::ActiveRenewableControllerAB) = value.D_dn
"""Get [`ActiveRenewableControllerAB`](@ref) `D_up`."""
get_D_up(value::ActiveRenewableControllerAB) = value.D_up
"""Get [`ActiveRenewableControllerAB`](@ref) `dP_lim`."""
get_dP_lim(value::ActiveRenewableControllerAB) = value.dP_lim
"""Get [`ActiveRenewableControllerAB`](@ref) `P_lim_inner`."""
get_P_lim_inner(value::ActiveRenewableControllerAB) = value.P_lim_inner
"""Get [`ActiveRenewableControllerAB`](@ref) `T_pord`."""
get_T_pord(value::ActiveRenewableControllerAB) = value.T_pord
"""Get [`ActiveRenewableControllerAB`](@ref) `P_ref`."""
get_P_ref(value::ActiveRenewableControllerAB) = value.P_ref
"""Get [`ActiveRenewableControllerAB`](@ref) `ext`."""
get_ext(value::ActiveRenewableControllerAB) = value.ext
"""Get [`ActiveRenewableControllerAB`](@ref) `states`."""
get_states(value::ActiveRenewableControllerAB) = value.states
"""Get [`ActiveRenewableControllerAB`](@ref) `n_states`."""
get_n_states(value::ActiveRenewableControllerAB) = value.n_states

"""Set [`ActiveRenewableControllerAB`](@ref) `bus_control`."""
set_bus_control!(value::ActiveRenewableControllerAB, val) = value.bus_control = val
"""Set [`ActiveRenewableControllerAB`](@ref) `from_branch_control`."""
set_from_branch_control!(value::ActiveRenewableControllerAB, val) = value.from_branch_control = val
"""Set [`ActiveRenewableControllerAB`](@ref) `to_branch_control`."""
set_to_branch_control!(value::ActiveRenewableControllerAB, val) = value.to_branch_control = val
"""Set [`ActiveRenewableControllerAB`](@ref) `branch_id_control`."""
set_branch_id_control!(value::ActiveRenewableControllerAB, val) = value.branch_id_control = val
"""Set [`ActiveRenewableControllerAB`](@ref) `Freq_Flag`."""
set_Freq_Flag!(value::ActiveRenewableControllerAB, val) = value.Freq_Flag = val
"""Set [`ActiveRenewableControllerAB`](@ref) `K_pg`."""
set_K_pg!(value::ActiveRenewableControllerAB, val) = value.K_pg = val
"""Set [`ActiveRenewableControllerAB`](@ref) `K_ig`."""
set_K_ig!(value::ActiveRenewableControllerAB, val) = value.K_ig = val
"""Set [`ActiveRenewableControllerAB`](@ref) `T_p`."""
set_T_p!(value::ActiveRenewableControllerAB, val) = value.T_p = val
"""Set [`ActiveRenewableControllerAB`](@ref) `fdbd_pnts`."""
set_fdbd_pnts!(value::ActiveRenewableControllerAB, val) = value.fdbd_pnts = val
"""Set [`ActiveRenewableControllerAB`](@ref) `fe_lim`."""
set_fe_lim!(value::ActiveRenewableControllerAB, val) = value.fe_lim = val
"""Set [`ActiveRenewableControllerAB`](@ref) `P_lim`."""
set_P_lim!(value::ActiveRenewableControllerAB, val) = value.P_lim = val
"""Set [`ActiveRenewableControllerAB`](@ref) `T_g`."""
set_T_g!(value::ActiveRenewableControllerAB, val) = value.T_g = val
"""Set [`ActiveRenewableControllerAB`](@ref) `D_dn`."""
set_D_dn!(value::ActiveRenewableControllerAB, val) = value.D_dn = val
"""Set [`ActiveRenewableControllerAB`](@ref) `D_up`."""
set_D_up!(value::ActiveRenewableControllerAB, val) = value.D_up = val
"""Set [`ActiveRenewableControllerAB`](@ref) `dP_lim`."""
set_dP_lim!(value::ActiveRenewableControllerAB, val) = value.dP_lim = val
"""Set [`ActiveRenewableControllerAB`](@ref) `P_lim_inner`."""
set_P_lim_inner!(value::ActiveRenewableControllerAB, val) = value.P_lim_inner = val
"""Set [`ActiveRenewableControllerAB`](@ref) `T_pord`."""
set_T_pord!(value::ActiveRenewableControllerAB, val) = value.T_pord = val
"""Set [`ActiveRenewableControllerAB`](@ref) `P_ref`."""
set_P_ref!(value::ActiveRenewableControllerAB, val) = value.P_ref = val
"""Set [`ActiveRenewableControllerAB`](@ref) `ext`."""
set_ext!(value::ActiveRenewableControllerAB, val) = value.ext = val


text:
@polly function kernel_fdtd_2d(ex, ey, hz, _fict_)
    nx,ny = size(ex)
    tmax = size(_fict_,1)

    for t = 1:tmax
        for j = 1:ny
            ey[1,j] = _fict_[t]
        end
        for i = 2:nx, j = 1:ny
            ey[i,j] = ey[i,j] - 0.5 * (hz[i,j] - hz[i-1,j])
        end
        for i = 1:nx, j = 2:ny
            ex[i,j] = ex[i,j] - 0.5 * (hz[i,j] - hz[i,j-1])
        end
        for i = 1:(nx - 1), j = 1:(ny - 1)
            hz[i,j] = hz[i,j] - 0.7 * (ex[i,j+1] - ex[i,j] + ey[i+1,j] - ey[i,j])
        end
    end
end

let
    tmax = 500
    nx = 1000
    ny = 1200

    ex = zeros(Float32,nx,ny)
    ey = zeros(Float32,nx,ny)
    hz = zeros(Float32,nx,ny)
    _fict_ = zeros(Float32,tmax)

    for i = 1:tmax
        _fict_[i] = i
    end
    for i = 1:nx, j = 1:ny
        ex[i,j] = (i*(j+1)) / nx
        ey[i,j] = (i*(j+2)) / ny
        hz[i,j] = (i*(j+3)) / nx
    end

    SUITE["fdtd-2d"] = @benchmarkable kernel_fdtd_2d(ex, ey, hz, _fict_) setup = (ex = copy($ex); ey = copy($ey); hz = copy($hz); _fict_ = copy($_fict_))
end


text:
function flatten(data::Vector)
    res = Int[]
    for x in data
        append!(res, x)
    end
    res
end

function encode(dict::IdDict, data::Vector)
    chars = Int[]
    tags = Int[]
    for sent in data
        for items in sent
            char, tag = items[1], items[2]
            charid = push!(dict, string(char))
            push!(chars, charid)
            push!(tags, parse(Int,tag))
        end
    end
    x = Vector{Int}[]
    y = Vector{Int}[]
    for i = 1:1000:length(chars)-1000
        push!(x, chars[i:i+1000])
        push!(y, tags[i:i+1000])
    end
    x, y
end

function train(t::Tokenizer, nepochs::Int, traindata::Vector, testdata::Vector)
    train_x, train_y = encode(t.dict, traindata)
    test_x, test_y = encode(t.dict, testdata)
    info("#dict: $(length(t.dict))")

    #opt = AdaGrad(0.01)
    opt = SGD(0.0001, momentum=0.9)
    for epoch = 1:nepochs
        println("epoch: $(epoch)")
        loss = fit(t.model, crossentropy, opt, train_x, train_y)
        println("loss: $(loss)")

        test_z = map(test_x) do x
            argmax(t.model(x).data, 1)
        end
        #acc = accuracy(flatten(test_y), flatten(test_z))
        accuracy_sent(flatten(test_y), flatten(test_z))

        #println("test acc.: $(acc)")
        println("")
    end
end

function accuracy_sent(golds::Vector{Int}, preds::Vector{Int})
    @assert length(golds) == length(preds)
    correct_g, correct_p = 0, 0
    total_g, total_p = 0, 0
    for i = 1:length(golds)
        if golds[i] == 4
            golds[i] == preds[i] && (correct_g += 1)
            total_g += 1
        end
        if preds[i] == 4
            golds[i] == preds[i] && (correct_p += 1)
            total_p += 1
        end
    end
    prec = correct_g / total_g
    recall = correct_p / total_p
    println("precision: $(prec)")
    println("recall: $(recall)")
    println("f-value: $(2*recall*prec/(recall+prec))")
end

function accuracy(golds::Vector{Int}, preds::Vector{Int})
    @assert length(golds) == length(preds)
    correct = 0
    total = 0
    for i = 1:length(golds)
        golds[i] == preds[i] && (correct += 1)
        total += 1
    end
    correct / total
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
#Generic runtime policy for Metropolis-Hastings runner
immutable MHRuntimePolicy{N<:Number,T<:AbstractFloat} <: AbstractPolicy
    runner::MHRunnerType
    initialize::InitializeFrom
    propose::ProposeFrom
    indicator::IndicatorType
    jobsegments::JobSegments
    chain::ChainType
    store::StoreDuring
    sampletype::Type{N}
    calculationtype::Type{T}
    MHRuntimePolicy(t,i,p,d,j,c,s) = new(t,i,p,d,j,c,s,N,T)
end

#Local function mapping number of proposals to the ProposeFrom type
_num2runner(n::Integer) =  (@assert n > 0 ; n>1?trait(:mhrunner,:generalized):trait(:mhrunner,:standard))
_num2propose(n::Integer) = (@assert n > 0 ; n>1?trait(:propose,:auxiliary):trait(:propose,:indicator))
_num2segments(n::Integer,s::Symbol) = (@assert n > 0 ; n>1?trait(:jobsegments,s):trait(:jobsegments,:none))

_policy(::Type{Val{:mh}},nprops::Int;initialize = :prior,indicator = :stationary,jobsegments = :workers,
        chain = :standard,store = :main,sampletype =Float64,calculationtype =Float64) =
    MHRuntimePolicy{sampletype,calculationtype}(_num2runner(nprops),trait(:initialize,initialize),_num2propose(nprops),trait(:indicator,indicator),
                                                _num2segments(nprops,jobsegments),trait(:chain,chain),trait(:store,store))

function show(io::IO,p::MHRuntimePolicy)
    println(io,"MHRuntimePolicy with traits:")
    println(io,"  runner = ",traitvalue(p.runner))
    println(io,"  initialize = ",traitvalue(p.initialize))
    println(io,"  propose = ",traitvalue(p.propose))
    println(io,"  indicator = ",traitvalue(p.indicator))
    println(io,"  jobsegments = ",traitvalue(p.jobsegments))
    println(io,"  chain = ",traitvalue(p.chain))
    println(io,"  store = ",traitvalue(p.store))
    println(io,"  sampletype = ",p.sampletype)
    println(io,"  calculationtype = ",p.calculationtype)
    println(io)
    nothing
end


text:
# make your scripts automatically re-activate your project
cd(@__DIR__)
using Pkg; Pkg.activate("../..")

# load packages
using ml_for_pt
using DelimitedFiles
using LaTeXStrings
using Plots
using Random
using Flux
ENV["GKSwstype"]="nul"

# set linear system size and corresponding data folder
L = 10
# L=60
data_folder = "../../data/ising/L="*string(L)*"/"

# set path to save folder
save_folder = "./results/L="*string(L)*"/"

# critical parameter value for the Ising model
p_crit = Float32(2/log(1+sqrt(2)))

# define parameter ranges
p_min = 0.05f0
p_max = 10.0f0
dp = 0.05f0
p_range = collect(p_min:dp:p_max)
p_range_LBC = collect(p_min-dp/2:dp:p_max+dp/2)

# import Monte Carlo data

# unique energies at p=p_1
energy = readdlm(data_folder*"E_T="*string(0)*".txt", eltype(dp))

# number of times each energy is drawn at T=T_1
num = readdlm(data_folder*"num_T="*string(0)*".txt", eltype(dp))

# iterate over entire parameter range
energies = [energy]
numbers = [num]
for index in 1:length(p_range)-1
  energyy = readdlm(data_folder*"E_T="*string(index)*".txt", eltype(dp))
  numm = readdlm(data_folder*"num_T="*string(index)*".txt", eltype(dp))
  push!(energies, energyy)
  push!(numbers, numm)
end
unique_energies = unique(vcat(energies...))
n_samples = sum(numbers[1])

# construct probability distribution from Monte Carlo data
# distr(x,p) gives probability to draw the sample x at parameter value p
# samples correspond to enumeration of unique energies with integers
distr,samples = MLP.distr_approx(energies, unique_energies, numbers, p_range)

# optionally run methods using the analytical expressions to compare with neural network results

# supervised learning using analytical expression
# returns optimal predictions, indicator, and loss
pred_opt_SL, indicator_opt_SL, loss_opt_SL = MLP.get_indicators_SL_analytical(samples, distr, p_range, dp, p_min, p_max)

# prediction-based method using analytical expression
# returns optimal predictions, indicator, and loss
pred_opt_PBM, indicator_opt_PBM, loss_opt_PBM = MLP.get_indicators_PBM_analytical(samples, distr, p_range, dp)

# learning by confusion using analytical expression
# returns optimal indicator and loss
indicator_opt_LBC, loss_opt_LBC = MLP.get_indicators_LBC_analytical(samples, distr, p_range, p_range_LBC)

# construct dataset for training neural networks
dataset = MLP.get_dataset(p_range, distr, samples)
dataset_train_SL = MLP.get_modified_dataset_train_SL(distr, samples, p_range, n_samples)

# standardize inputs
inputs_one_hot = reshape(unique_energies, 1, length(unique_energies))
mean_train,std_train = MLP.get_dataset_stats(dataset, inputs_one_hot, length(p_range))
inputs_one_hot_stand = (inputs_one_hot.-mean_train)./std_train
mean_train_SL,std_train_SL = MLP.get_dataset_stats(dataset_train_SL, inputs_one_hot, length(p_range))
inputs_one_hot_stand_SL = (inputs_one_hot.-mean_train_SL)./std_train_SL

# supervised learning using neural networks

# initialize neural network
n_nodes_SL = [size(inputs_one_hot)[1], 64, 64, 64, 2]
NN_SL = Chain(
  Dense(n_nodes_SL[1], n_nodes_SL[2], relu),
  Dense(n_nodes_SL[2], n_nodes_SL[3], relu),
  Dense(n_nodes_SL[3], n_nodes_SL[4], relu),
  Dense(n_nodes_SL[4], n_nodes_SL[end], identity))
pnn_SL, re_SL = Flux.destructure(NN_SL)

# set hyperparameters
lr_SL = 0.001f0
epochs_SL = 10000
saveat_SL = 100
opt_SL = ADAM(lr_SL)
verbose = true

# train neural network
# returns predictions, indicators, and loss saved at epochs specified by saveat_SL variable
pred_logger_SL, losses_SL, NN_logger_SL = MLP.get_indicators_SL_numerical(pnn_SL,  re_SL, dataset_train_SL, dataset, epochs_SL, p_range, dp, p_min, p_max, opt_SL, inputs_one_hot_stand, verbose=verbose, saveat=saveat_SL)

# extract results at last save point
pred_NN_SL = pred_logger_SL[end][1]
indicator_NN_SL = pred_logger_SL[end][2]
loss_NN_SL = pred_logger_SL[end][3]

# save results of supervised learning in save_folder
open(save_folder*"hyperparams_SL.txt", "w") do io
  writedlm(io, [saveat_SL, lr_SL, epochs_SL])
end

open(save_folder*"NN_logger_SL.txt", "w") do io
  writedlm(io, NN_logger_SL)
end

for i in 1:length(pred_logger_SL)
  open(save_folder*"pred_logger_SL_"*string(i)*".txt", "w") do io
  writedlm(io, pred_logger_SL[i][1])
  end

  open(save_folder*"indicator_logger_SL_"*string(i)*".txt", "w") do io
  writedlm(io, pred_logger_SL[i][2])
  end

  open(save_folder*"loss_logger_SL_"*string(i)*".txt", "w") do io
  writedlm(io, pred_logger_SL[i][3])
  end

end

open(save_folder*"pred_NN_SL.txt", "w") do io
  writedlm(io, pred_NN_SL)
end

open(save_folder*"indicator_NN_SL.txt", "w") do io
  writedlm(io, indicator_NN_SL)
end

open(save_folder*"loss_NN_SL.txt", "w") do io
  writedlm(io, loss_NN_SL)
end

open(save_folder*"NN_structure_SL.txt", "w") do io
  writedlm(io, n_nodes_SL)
end

open(save_folder*"losses_SL.txt", "w") do io
  writedlm(io, losses_SL)
end

# plot and save results of supervised learning in save_folder
plt = plot(p_range, pred_NN_SL, dpi=300, c="black", label="NN", ylims=(0.0,1.0))
plot!(p_range, pred_opt_SL, dpi=300, c="blue", label="analytical", ylims=(0.0,1.0))
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$\hat{y}_{\mathrm{SL}}$")
savefig(save_folder*"pred_NN_SL.png")

plt = plot(p_range[2:end-1], indicator_NN_SL, dpi=300, c="black", label="NN")
plot!(p_range[2:end-1], indicator_opt_SL, dpi=300, c="blue", label="analytical")
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$I_{\mathrm{SL}}$")
savefig(save_folder*"indicator_NN_SL.png")

plt = plot(collect(1:length(losses_SL)), losses_SL, dpi=300, c="black", label="NN")
hline!([loss_opt_SL], label="analytical", c="blue")
xlabel!("Epochs")
ylabel!(L"$\mathcal{L}_{\mathrm{SL}}$")
savefig(save_folder*"loss_NN_SL.png")

# prediction-based method using neural networks

# initialize neural network
n_nodes_PBM = [size(inputs_one_hot)[1], 64, 64, 64, 1]
NN_PBM = Chain(
  Dense(n_nodes_PBM[1], n_nodes_PBM[2], relu),
  Dense(n_nodes_PBM[2], n_nodes_PBM[3], relu),
  Dense(n_nodes_PBM[3], n_nodes_PBM[4], relu),
  Dense(n_nodes_PBM[4], n_nodes_PBM[end], identity))
pnn_PBM, re_PBM = Flux.destructure(NN_PBM)

# set hyperparameters
lr_PBM = 0.001f0
epochs_PBM = 1000
saveat_PBM = 100
opt_PBM = ADAM(lr_PBM)
verbose = true

# train neural network
# returns predictions, indicators, and loss saved at epochs specified by saveat_PBM variable
pred_logger_PBM, losses_PBM, NN_logger_PBM =MLP.get_indicators_PBM_numerical(pnn_PBM, re_PBM, dataset, epochs_PBM, p_range, dp, opt_PBM, inputs_one_hot_stand, verbose=verbose, saveat=saveat_PBM)

# extract results at last save point
pred_NN_PBM = pred_logger_PBM[end][1]
indicator_NN_PBM = pred_logger_PBM[end][2]
loss_NN_PBM = pred_logger_PBM[end][3]

# save results of prediction-based method in save_folder
open(save_folder*"hyperparams_PBM.txt", "w") do io
  writedlm(io, [saveat_PBM, lr_PBM, epochs_PBM])
end

open(save_folder*"NN_logger_PBM.txt", "w") do io
  writedlm(io, NN_logger_PBM)
end

for i in 1:length(pred_logger_PBM)
  open(save_folder*"pred_logger_PBM_"*string(i)*".txt", "w") do io
  writedlm(io, pred_logger_PBM[i][1])
  end

  open(save_folder*"indicator_logger_PBM_"*string(i)*".txt", "w") do io
  writedlm(io, pred_logger_PBM[i][2])
  end

  open(save_folder*"loss_logger_PBM_"*string(i)*".txt", "w") do io
  writedlm(io, pred_logger_PBM[i][3])
  end

end

open(save_folder*"pred_NN_PBM.txt", "w") do io
  writedlm(io, pred_NN_PBM)
end

open(save_folder*"indicator_NN_PBM.txt", "w") do io
  writedlm(io, indicator_NN_PBM)
end

open(save_folder*"loss_NN_PBM.txt", "w") do io
  writedlm(io, loss_NN_PBM)
end

open(save_folder*"NN_structure_PBM.txt", "w") do io
  writedlm(io, n_nodes_PBM)
end

open(save_folder*"losses_PBM.txt", "w") do io
  writedlm(io, losses_PBM)
end

# plot and save results of prediction-based method in save_folder
plt = plot(p_range, pred_NN_PBM, dpi=300, c="black", label="NN", ylims=(p_range[1],p_range[end]))
plot!(p_range, pred_opt_PBM, dpi=300, c="blue", label="analytical", ylims=(p_range[1],p_range[end]))
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$\hat{y}_{\mathrm{PBM}}$")
savefig(save_folder*"pred_NN_PBM.png")

plt = plot(p_range[2:end-1], indicator_NN_PBM, dpi=300, c="black", label="NN")
plot!(p_range[2:end-1], indicator_opt_PBM, dpi=300, c="blue", label="analytical")
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$I^{\mathrm{NN}}_{\mathrm{PBM}}$")
savefig(save_folder*"indicator_NN_PBM.png")

plt = plot(collect(1:length(losses_PBM)), losses_PBM, dpi=300, c="black", label="NN")
hline!([loss_opt_PBM], label="analytical", c="blue")
xlabel!("Epochs")
ylabel!(L"$\mathcal{L}_{\mathrm{PBM}}$")
savefig(save_folder*"loss_NN_PBM.png")

# learning by confusion using neural networks

# initialize neural network
n_nodes_LBC = [size(inputs_one_hot)[1], 64, 64, 64, 2]
NN_LBC = Chain(
  Dense(n_nodes_LBC[1], n_nodes_LBC[2], relu),
  Dense(n_nodes_LBC[2], n_nodes_LBC[3], relu),
  Dense(n_nodes_LBC[3], n_nodes_LBC[4], relu),
  Dense(n_nodes_LBC[4], n_nodes_LBC[end], identity))
pnn_LBC, re_LBC = Flux.destructure(NN_LBC)

# set hyperparameters
lr_LBC = 0.001f0
epochs_LBC = 100
saveat_LBC = epochs_LBC
opt_LBC = ADAM(lr_LBC)
verbose = true

indicator_NN_LBC = []
loss_NN_LBC = []
for indx in collect(1:length(p_range_LBC))
  pred_logger_LBC, losses_LBC, NN_logger_LBC = MLP.get_indicators_LBC_numerical_fixed_p(deepcopy(pnn_LBC), re_LBC, dataset, epochs_LBC, p_range, dp, opt_LBC, p_range_LBC, indx, inputs_one_hot_stand, saveat=saveat_LBC, verbose=verbose)

  indicator_NN_LBC_p = pred_logger_LBC[end][1][1]
  push!(indicator_NN_LBC, indicator_NN_LBC_p)
  loss_NN_LBC_p = pred_logger_LBC[end][2][1]
  push!(loss_NN_LBC, loss_NN_LBC_p)
end

# save results of learning by confusion in save_folder
open(save_folder*"hyperparams_LBC.txt", "w") do io
  writedlm(io, [saveat_LBC, lr_LBC, epochs_LBC])
end

open(save_folder*"indicator_NN_LBC.txt", "w") do io
  writedlm(io, indicator_NN_LBC)
end

open(save_folder*"loss_NN_LBC.txt", "w") do io
  writedlm(io, loss_NN_LBC)
end

open(save_folder*"NN_structure_LBC.txt", "w") do io
  writedlm(io, n_nodes_LBC)
end

# plot and save results of learning by confusion in save_folder
plt = plot(p_range_LBC, indicator_NN_LBC, dpi=300, c="black", label="NN")
plot!(p_range_LBC, indicator_opt_LBC, dpi=300, c="blue", label="analytical")
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$I_{\mathrm{LBC}}$")
savefig(save_folder*"indicator_NN_LBC.png")

plt = plot(p_range_LBC, loss_NN_LBC, dpi=300, c="black", label="NN")
plot!(p_range_LBC, loss_opt_LBC, dpi=300, c="blue", label="analytical")
vline!([p_crit], label=L"p_{c}", c="red")
xlabel!(L"$p$")
ylabel!(L"$\mathcal{L}_{\mathrm{LBC}}$")
savefig(save_folder*"loss_NN_LBC.png")


text:
include("transit_poly.jl")

# Randomizer seed
using Random
Random.seed!(42)

nd = 2+ceil(Int64,rand()*20)
r=rand(); b=rand()*(1+r); d = rand(nd)
#r=1.2*rand(); b=0.0; d = rand(nd)
#r = rand(); b = rand(); d = [1.0,1.0,10.0]
#r=0.7821623254303613; b= 0.5; d = [0.0,1.0]

flux = transit_poly_d(r,b,d)
@time flux = transit_poly_d(r,b,d)

# Now, integrate by hand:
function transit_poly_int_d(r,b,d)
s_1 = maximum([0.0,b-r]); s_2=minimum([1.0,b+r])
ns = 500; ds=(s_2-s_1)/ns
s = linspace(s_1+.5*ds,s_2-.5*ds,ns)
fobs = zero(r)
for j=1:ns
  if s[j] < r-b
    dphi = 2pi
  else
    dphi = 2*acos((s[j]^2+b^2-r^2)/(2*b*s[j]))
  end
  z = sqrt(1-s[j]^2)
  imu = d[1] + d[2]*z
  for n=2:length(d)-1
    imu += d[n+1]*(z^n*(n+2)-n*z^(n-2))
  end
  fobs += s[j]*ds*dphi*imu
end
fobs /= pi*(d[1]+2/3*d[2])
fobs = 1-fobs
#println("r: ",r," b: ",b,"f_an: ",flux," flux: ",flux0," f_num: ",fobs," f_exp: ",pi*(d[1]+2/3*d[2]))
return fobs
end
f_num = transit_poly_int_d(r,b,d)
@time f_num = transit_poly_int_d(r,b,d)
println("r: ",r," b: ",b," f_an: ",flux," f_num: ",f_num)

#nphi = 500; ns = 500; dphi = 2pi/nphi; ds=1/ns
#phi = linspace(.5*dphi,2pi-.5*dphi,nphi); s = linspace(.5*ds,1-.5*ds,ns)
#flux0 = zero(r); fobs = zero(r)
#for i=1:nphi, j=1:ns
#  x=s[j]*cos(phi[i]); y = s[j]*sin(phi[i])
#  z = sqrt(1-s[j]^2)
#  imu = d[1]+d[2]*z
#  for k=2:length(d)-1
#    imu += d[k+1]*(z^k*(k+2)-k*z^(k-2))
#  end
#  flux0 += s[j]*ds*dphi*imu
#  if (y-b)^2+x^2 > r^2
#    fobs += s[j]*ds*dphi*imu
#  end
#end
##println("r: ",r," b: ",b,"f_an: ",flux," flux: ",flux0," f_num: ",fobs/flux0," f_exp: ",pi*(d[1]+2/3*d[2]))
#println("r: ",r," b: ",b," f_an: ",flux," f_num: ",fobs/flux0)


text:
abstract type AbstractScene end

struct PlanarScene <: AbstractScene

end


text:
module TestCBSOData4

    using Test
    using CBSOData4

    include("01_tables_interface.jl")
    include("02_get_meta.jl")
    include("03_get_tables.jl")
end


text:
function transform(x, q, shape)
    scale_transform = LinearMap(diagm(shape.scale))
    x_transform = Translation(x + vrotate(shape.xoffset, q))
    q_transform = LinearMap(q * shape.qoffset)

    return compose(x_transform, q_transform, scale_transform)
end

MeshCat.js_scaling(s::AbstractVector) = s
MeshCat.js_position(p::AbstractVector) = p

function preparevis!(storage::Storage{T,N}, id, shape, animation, shapevisualizer, framevisualizer, showshape, showframes) where {T,N}
    if showshape
        for i=1:N
            x = storage.x[id][i]
            q = storage.q[id][i]
            # TODO currently setting props directly because MeshCat/Rotations doesn't convert scaled rotation properly.
            # If this changes, do similarily to origin
            atframe(animation, i) do
                setprop!(shapevisualizer, "scale", js_scaling(shape.scale))
                setprop!(shapevisualizer, "position", js_position(x + vrotate(shape.xoffset, q)))
                setprop!(shapevisualizer, "quaternion", js_quaternion(q * shape.qoffset))
            end
        end
    end

    if showframes
        for i=1:N
            x = storage.x[id][i]
            q = storage.q[id][i]
            atframe(animation, i) do
                setprop!(framevisualizer, "scale", js_scaling(SA{T}[1;1;1]))
                setprop!(framevisualizer, "position", js_position(x))
                setprop!(framevisualizer, "quaternion", js_quaternion(q))
            end
        end
    end

    return
end

function MeshCat.setobject!(subvisshape, visshape, shape::Shape)
    setobject!(subvisshape, visshape, MeshPhongMaterial(color=shape.color))
end

function MeshCat.setobject!(subvisshape, visshape, shape::ConstrainedDynamics.Mesh)
    if visshape.mtl_library == ""
        visshape = MeshFileGeometry(visshape.contents, visshape.format)
        setobject!(subvisshape, visshape, MeshPhongMaterial(color=shape.color))
    else
        setobject!(subvisshape, visshape)
    end
end

"""
    visualize(mechanism, storage; env, showframes)

Visualize a `mechanism` with a trajectory stored in `storage`.

# Available kwargs
* `showframes`: Display the coordinate frames of the bodies.
* `env`:        Choose the visualization environment ("blink", "browser", "editor").
"""
function visualize(mechanism::AbstractMechanism, storage::Storage{T,N}; env::String = "blink", showframes::Bool = false) where {T,N}
    storage = deepcopy(storage)
    bodies = mechanism.bodies
    origin = mechanism.origin
    if showframes
        triads = [Triad(0.33) for i=1:length(bodies)]
    end

    vis = Visualizer()
    setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0))
    setprop!(vis["/Background"], "bottom_color", RGBA(0.5, 0.6, 0.7))
    # Somehow delete! doesn't work in a function call, so set axes to not visible for now
    setvisible!(vis["/Axes"],false)
    # delete!(vis["/Axes"])

    if env == "blink"
        open(vis, Blink.Window())
    elseif env == "browser"
        open(vis)
    elseif env == "editor"
        # Needs to be returned to show. Is there a better way of doing this?
    else
        @error "No valid visualization environment option selected."
    end

    framerate = Int64(round(1/mechanism.Δt))
    animation = Animation(Dict{MeshCat.SceneTrees.Path,MeshCat.AnimationClip}(), framerate)

    for (id,body) in enumerate(bodies)
        shape = body.shape
        visshape = convertshape(shape)
        subvisshape = nothing
        subvisframe = nothing
        showshape = false
        if visshape !== nothing
            subvisshape = vis["bodies/body:"*string(id)]
            setobject!(subvisshape,visshape,shape)
            showshape = true
        end
        if showframes
            subvisframe = vis["frames/body:"*string(id)]
            setobject!(subvisframe, triads[id])
        end
        preparevis!(storage, id, shape, animation, subvisshape, subvisframe, showshape, showframes)
    end

    id = origin.id
    shape = origin.shape
    visshape = convertshape(shape)
    if visshape !== nothing
        subvisshape = vis["bodies/origin:"*string(id)]
        setobject!(subvisshape,visshape,shape)
        shapetransform = transform(szeros(T,3), one(UnitQuaternion{T}), shape)
        settransform!(subvisshape, shapetransform)
    end
    if showframes
        subvisframe = vis["frames/origin:"*string(id)]
        setobject!(subvisframe, Triad(0.5))
    end

    setanimation!(vis, animation)
    env == "editor" ? (return render(vis)) : (return vis)
end


text:
# You should set up a Git repo with an example Julia package for the integration tests.
# You should create a bot user GitHub account, and give the bot user GitHub account push access to the testing repo.
# You should generate a personal access token (PAT) for the bot user GitHub account.
# You should add the PAT as an encrypted Travis CI environment variable named `BCBI_TEST_USER_GITHUB_TOKEN`.
# You should generate an SSH deploy key. Upload the public key as a deploy key for the testing repo.
# You should make the private key an encrypted Travis CI environment variable named `COMPATHELPER_PRIV`. Probably you will want to Base64-encode it.

COMPATHELPER_INTEGRATION_TEST_REPO = ENV["COMPATHELPER_INTEGRATION_TEST_REPO"]::String
TEST_USER_GITHUB_TOKEN = ENV["BCBI_TEST_USER_GITHUB_TOKEN"]::String
sleep(5)
auth = CompatHelper.my_retry(() -> GitHub.authenticate(TEST_USER_GITHUB_TOKEN))
sleep(5)
whoami = username(auth)
repo_url_without_auth = "https://github.com/$(COMPATHELPER_INTEGRATION_TEST_REPO)"
repo_url_with_auth = "https://$(whoami):$(TEST_USER_GITHUB_TOKEN)@github.com/$(COMPATHELPER_INTEGRATION_TEST_REPO)"
sleep(5)
repo = GitHub.repo(COMPATHELPER_INTEGRATION_TEST_REPO; auth = auth)
sleep(5)
api = GitHub.GitHubWebAPI(HTTP.URI("https://api.github.com"))
sleep(5)
Test.@test success(`git --version`)
@info("Authenticated to GitHub as \"$(whoami)\"")

sleep(5)

_delete_branches_older_than = Dates.Hour(3)
delete_old_pull_request_branches(
    repo_url_with_auth,
    _delete_branches_older_than
)

sleep(5)

const GLOBAL_PR_TITLE_PREFIX = Random.randstring(8)

with_master_branch(templates("master_1"), "master"; repo_url = repo_url_with_auth) do master_1
    withenv("GITHUB_REPOSITORY" => COMPATHELPER_INTEGRATION_TEST_REPO,
            "GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN) do
        ci_cfg = CompatHelper.GitHubActions(whoami, "41898282+github-actions[bot]@users.noreply.github.com")
        sleep(30)
        CompatHelper.main(ENV, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-1c] ",
                          master_branch = master_1,
                          keep_existing_compat = true,
                          drop_existing_compat = true)
        sleep(5)
    end
end

with_master_branch(templates("master_1"), "master"; repo_url = repo_url_with_auth) do master_1
    withenv("GITHUB_REPOSITORY" => COMPATHELPER_INTEGRATION_TEST_REPO,
            "GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN) do
        env = deepcopy(Dict(ENV))
        delete!(env, "COMPATHELPER_PRIV")
        ci_cfg = CompatHelper.GitHubActions(whoami, "41898282+github-actions[bot]@users.noreply.github.com")
        sleep(30)
        CompatHelper.main(env, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-1c] ",
                          master_branch = master_1,
                          keep_existing_compat = true,
                          drop_existing_compat = true)
        sleep(5)
    end
end

with_master_branch(templates("master_2"), "master"; repo_url = repo_url_with_auth) do master_2
    withenv("GITHUB_REPOSITORY" => COMPATHELPER_INTEGRATION_TEST_REPO,
            "GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN) do
        ci_cfg = CompatHelper.GitHubActions(whoami, "41898282+github-actions[bot]@users.noreply.github.com")
        sleep(30)
        CompatHelper.main(ENV, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-2c] ",
                          master_branch = master_2,
                          keep_existing_compat = true,
                          drop_existing_compat = true)
        sleep(5)
    end
end

with_master_branch(templates("master_3"), "master"; repo_url = repo_url_with_auth) do master_3
    withenv("GITHUB_REPOSITORY" => COMPATHELPER_INTEGRATION_TEST_REPO,
            "GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN) do
        ci_cfg = CompatHelper.GitHubActions(whoami, "41898282+github-actions[bot]@users.noreply.github.com")
        sleep(30)
        CompatHelper.main(ENV, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-3a] ",
                          master_branch = master_3,
                          keep_existing_compat = false,
                          drop_existing_compat = true)
        sleep(30)
        CompatHelper.main(ENV, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-3b] ",
                          master_branch = master_3,
                          keep_existing_compat = true,
                          drop_existing_compat = false)
        sleep(30)
        CompatHelper.main(ENV, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-3c] ",
                          master_branch = master_3,
                          keep_existing_compat = true,
                          drop_existing_compat = true)
        sleep(30)
        CompatHelper.main(ENV, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-3c] ",
                          master_branch = master_3,
                          keep_existing_compat = true,
                          drop_existing_compat = true)
        sleep(5)
    end
end

with_master_branch(templates("master_4"), "master"; repo_url = repo_url_with_auth) do master_4
    withenv("GITHUB_REPOSITORY" => COMPATHELPER_INTEGRATION_TEST_REPO,
            "GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN) do
        ci_cfg = CompatHelper.GitHubActions(whoami, "41898282+github-actions[bot]@users.noreply.github.com")
        sleep(30)
        CompatHelper.main(ENV, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-4c] ",
                          master_branch = master_4,
                          keep_existing_compat = true,
                          drop_existing_compat = true)
        sleep(5)
    end
end

with_master_branch(templates("master_5"), "master"; repo_url = repo_url_with_auth) do master_5
    withenv("GITHUB_REPOSITORY" => COMPATHELPER_INTEGRATION_TEST_REPO,
            "GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN) do
        ci_cfg = CompatHelper.GitHubActions(whoami, "41898282+github-actions[bot]@users.noreply.github.com")
        sleep(30)
        CompatHelper.main(ENV, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-5a] ",
                          master_branch = master_5,
                          keep_existing_compat = false,
                          drop_existing_compat = true)
        sleep(5)
    end
end

# Same as master_1 and master_5, but in sub-directories
with_master_branch(templates("master_6"), "master"; repo_url = repo_url_with_auth) do master_6
    withenv("GITHUB_REPOSITORY" => COMPATHELPER_INTEGRATION_TEST_REPO,
            "GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN) do
        ci_cfg = CompatHelper.GitHubActions(whoami, "41898282+github-actions[bot]@users.noreply.github.com")
        sleep(30)
        CompatHelper.main(ENV, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-1c] ",
                          master_branch = master_6,
                          keep_existing_compat = true,
                          drop_existing_compat = true,
                          subdirs = ["subdir_1", "subdir_2"])
        sleep(5)
    end
end

with_master_branch(templates("master_7"), "master"; repo_url = repo_url_with_auth) do master_7
    withenv("GITHUB_REPOSITORY" => COMPATHELPER_INTEGRATION_TEST_REPO,
            "GITHUB_TOKEN" => TEST_USER_GITHUB_TOKEN) do
        ci_cfg = CompatHelper.GitHubActions(whoami, "41898282+github-actions[bot]@users.noreply.github.com")
        sleep(30)
        CompatHelper.main(ENV, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-7a] ",
                          master_branch = master_7,
                          keep_existing_compat = false,
                          drop_existing_compat = true,
                          bump_compat_containing_equality_specifier = false)
        sleep(30)
        CompatHelper.main(ENV, ci_cfg;
                          pr_title_prefix = "$(GLOBAL_PR_TITLE_PREFIX) [test-7b] ",
                          master_branch = master_7,
                          keep_existing_compat = false,
                          drop_existing_compat = true,
                          bump_compat_containing_equality_specifier = true)
        sleep(5)
    end
end

sleep(5)

all_prs = CompatHelper.get_all_pull_requests(api, repo, "open";
                                             auth = auth,
                                             per_page = 3,
                                             page_limit = 3)

sleep(5)


text:
function leaf_bytes(seen::Set{UInt64}, inode::INode{K,V}) where {K,V}
    id = objectid(inode)
    in(id, seen) && return 0
    push!(seen, id)
    n = 0
    for child in inode.entries
        n += leaf_bytes(seen, child)
    end
    return n
end

function leaf_bytes(seen::Set{UInt64}, leaf::Leaf{K,V}) where {K,V}
    id = objectid(leaf)
    in(id, seen) && return 0
    push!(seen, id)
    n = Base.summarysize(leaf)
    if length(leaf.big_slots) == 0
        n -= Base.summarysize(leaf.big_slots)
    end
    return n
end

function dump_tree(seen::Set{UInt64}, prefix::String, inode::INode{K,V}) where {K,V}
    id = objectid(inode)
    if in(id,seen)
        sflag="*"
    else
        push!(seen, id)
        sflag=" "
    end

    count=0
    for i=1:length(inode.entries)
        child_prefix = "$(prefix)$(i-1)|"
        count += dump_tree(seen, child_prefix, inode.entries[i])
    end
    println("$(sflag) $(prefix) \t\t$(count)")
    sflag=="*" && return 0
    return count
end

function dump_tree(seen::Set{UInt64}, prefix::String, leaf::Leaf{K,V}) where {K,V}
    id = objectid(leaf)
    if in(id,seen)
        sflag="*"
    else
        push!(seen, id)
        sflag=" "
    end

    count=0
    num_slots = length(leaf.entries)
    for i=1:num_slots
        if _slot_occupied(leaf, UInt64(i), num_slots)
            count += 1
        end
    end
    println("$(sflag) $(prefix) $(count) in $(num_slots) slots")
    sflag=="*" && return 0
    return count
end


text:
include("kern_proc.jl")

d = 10
n = [50, 100, 500, 1000]

srand(1)

cov_procs = Array(Proc, 0)
push!(cov_procs, KernelTest(LinIso(0.0), d, GaussianProcesses.cov))
push!(cov_procs, KernelTest(LinArd(rand(d)), d, GaussianProcesses.cov))
push!(cov_procs, KernelTest(Noise(0.0), d, GaussianProcesses.cov))
push!(cov_procs, KernelTest(SEIso(1.0, 1.0) + Noise(0.0), d, GaussianProcesses.cov))
push!(cov_procs, KernelTest(SEArd(rand(d), 1.0) + Noise(0.0), d, GaussianProcesses.cov))
cov_table = run(cov_procs, n)

grad_stack_procs = Array(Proc, 0)
push!(grad_stack_procs, KernelTest(LinIso(0.0), d, GaussianProcesses.grad_stack))
push!(grad_stack_procs, KernelTest(LinArd(rand(d)), d, GaussianProcesses.grad_stack))
push!(grad_stack_procs, KernelTest(Noise(0.0), d, GaussianProcesses.grad_stack))
push!(grad_stack_procs, KernelTest(SEIso(1.0, 1.0) + Noise(0.0), d, GaussianProcesses.grad_stack))
push!(grad_stack_procs, KernelTest(SEIso(1.0, 1.0) + Noise(0.0), d, GaussianProcesses.grad_stack))
push!(grad_stack_procs, KernelTest(SEArd(rand(d), 1.0) + Noise(0.0), d, GaussianProcesses.grad_stack))
grad_stack_table = run(grad_stack_procs, n)

println("\ncov calculations")
show(cov_table)
println("\grad_stack calculations")
show(grad_stack_table)


text:
importall ProfileAligner

println("testing AlignmentMatrix...")
empty_matrix  = AlignmentMatrix{Float64}(3, 7)
@test size(empty_matrix.matrix) == size(empty_matrix.path) == (3, 7)

println("testing Profile...")
profile1 = Profile{Float64}("ACQ")
@test typeof(profile1) == Profile{Float64}
@test score(profile1, profile1, 1, 1) == 1.0

dd2 = Profile{Float64}(['A' 'C'; 'Q' 'A'; 'B' 'Q'], ["str1", "str2"])
#println(dd2)
# TODO: add alignment test: when align profile to itself, frequency matrix shouldn't change
#println(ProfileAligner.align(dd2, dd2))
@test length(getstrings(dd2)) == 2
#for str in ProfileAligner.getstrings(dd2)
#  println(str)
#end

println("testing Profile alignments...")

str1 = "CAPTIAN"
str2 = "CAPTA"

p1 = Profile{Float64}(str1)
p2 = Profile{Float64}(str2)
aligned = align(p1, p2)
@test length(getstrings(aligned)) == 2

aligned_expected = [
  FastaRecord("","CAPTIAN"),
  FastaRecord("","CAPT-A-")
]
aligned_actual = getstrings(aligned)
println("retrieving aligned strings for profile...")
for i in 1:2
  @test aligned_expected[i] == aligned_actual[i]
end

aligned2 = ProfileAligner.align(aligned, ProfileAligner.Profile{Float64}("CAPTAIN"))
aligned2_expected = [
  FastaRecord("","CAPTIA-N"),
  FastaRecord("","CAPT-A--"),
  FastaRecord("","CAPT-AIN")
]
aligned2_actual = getstrings(aligned2)

println("retrieving strings for profile2...")
for i in 1:3
  @test aligned2_expected[i] == aligned2_actual[i]
end

aligned3 = ProfileAligner.align(aligned2, ProfileAligner.Profile{Float64}("APTAN"))
aligned3_expected = [
  FastaRecord("","CAPTIA-N"),
  FastaRecord("","CAPT-A--"),
  FastaRecord("","CAPT-AIN"),
  FastaRecord("","-APT-A-N")
]
aligned3_actual = getstrings(aligned3)
println("retrieving strings for profile3...")
for i in 1:4
  @test aligned3_expected[i] == aligned3_actual[i]
end


text:
# for Julia v0.6
# Pkg.add("ScikitLearn")
using DecisionTree
using ScikitLearn
using PyPlot
using Random

# Create a random dataset
Random.seed!(42)
X = sort(5 * rand(80))
XX = reshape(X, 80, 1)
y = sin.(X)
y[1:5:end] += 3 * (0.5 - rand(16))

# Fit regression model
regr_1 = DecisionTreeRegressor()
regr_2 = DecisionTreeRegressor(pruning_purity_threshold=0.05)
regr_3 = RandomForestRegressor(ntrees=20)
fit!(regr_1, XX, y)
fit!(regr_2, XX, y)
fit!(regr_3, XX, y)

# Predict
X_test = 0:0.01:5.0
y_1 = predict(regr_1, hcat(X_test))
y_2 = predict(regr_2, hcat(X_test))
y_3 = predict(regr_3, hcat(X_test))

# Plot the results
scatter(X, y, c="k", label="data")
plot(X_test, y_1, c="g", label="no pruning", linewidth=2)
plot(X_test, y_2, c="r", label="pruning_purity_threshold=0.05", linewidth=2)
plot(X_test, y_3, c="b", label="RandomForestClassifier", linewidth=2)
xlabel("data")
ylabel("target")
title("Decision Tree Regression")
legend(prop=Dict("size" => 10));

updateSettingsMod!(sett, crit="difference")
strs, resmT = dtm(dtmtable, sett)
updateSettingsMod!(sett, crit="maxabsvalue")
strs, resmT = dtm(dtmtable, sett)
updateSettingsMod!(sett, crit="maxminusvalue")
strs, resmT = dtm(dtmtable, sett)
updateSettingsMod!(sett, crit="maxvalue")
strs, resmT = dtm(dtmtable, sett)


text:
module TransientFEOperatorsTests

using Gridap
using ForwardDiff
using LinearAlgebra
using Test
using GridapODEs.ODETools
using GridapODEs.TransientFETools
using Gridap.FESpaces: get_algebraic_operator

# using GridapODEs.ODETools: ThetaMethodLinear

import Gridap: ∇
import GridapODEs.TransientFETools: ∂t

θ = 0.4

# Analytical functions
u(x,t) = (1.0-x[1])*x[1]*(1.0-x[2])*x[2]*(t+3.0)
u(t::Real) = x -> u(x,t)
v(x) = t -> u(x,t)

f(t) = x -> ∂t(u)(x,t)-Δ(u(t))(x)

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

order = 2

reffe = ReferenceFE(lagrangian,Float64,order)
V0 = FESpace(
  model,
  reffe,
  conformity=:H1,
  dirichlet_tags="boundary")

U = TransientTrialFESpace(V0,u)

Ω = Triangulation(model)
degree = 2*order
dΩ = Measure(Ω,degree)

a(u,v) = ∫(∇(v)⊙∇(u))dΩ
m(u,v) = ∫(v*u)dΩ
b(v,t) = ∫(v*f(t))dΩ

res(t,u,v) = a(u,v) + m(∂t(u),v) - b(v,t)
jac(t,u,du,v) = a(du,v)
jac_t(t,u,dut,v) = m(dut,v)

op = TransientFEOperator(res,jac,jac_t,U,V0)

t0 = 0.0
tF = 1.0
dt = 0.1

U0 = U(0.0)
uh0 = interpolate_everywhere(u(0.0),U0)

ls = LUSolver()
ode_solver = ThetaMethod(ls,dt,θ)

sol_t = solve(ode_solver,op,uh0,t0,tF)

l2(w) = w*w

tol = 1.0e-6
_t_n = t0

for (uh_tn, tn) in sol_t
  global _t_n
  _t_n += dt
  @test tn≈_t_n
  e = u(tn) - uh_tn
  el2 = sqrt(sum( ∫(l2(e))dΩ ))
  @test el2 < tol
end

#

u0 = get_free_dof_values(uh0)
uf = get_free_dof_values(uh0)

odeop = get_algebraic_operator(op)

ode_cache = allocate_cache(odeop)
vθ = similar(u0)
nl_cache = nothing

# tf = t0+dt

ode_solver.θ == 0.0 ? dtθ = dt : dtθ = dt*ode_solver.θ
tθ = t0+dtθ
ode_cache = update_cache!(ode_cache,odeop,tθ)

using GridapODEs.ODETools: ThetaMethodNonlinearOperator
nlop = ThetaMethodNonlinearOperator(odeop,tθ,dtθ,u0,ode_cache,vθ)

nl_cache = solve!(uf,ode_solver.nls,nlop,nl_cache)

K = nl_cache.A
h = nl_cache.b

# Steady version of the problem to extract the Laplacian and mass matrices
# tf = 0.1
tf = tθ
Utf = U(tf)
# fst(x) = -Δ(u(tf))(x)
fst(x) = f(tf)(x)
a(u,v) = ∫(∇(v)⊙∇(u))dΩ

function extract_matrix_vector(a,fst)
  btf(v) = ∫(v*fst)dΩ
  op = AffineFEOperator(a,btf,Utf,V0)
  ls = LUSolver()
  solver = LinearFESolver(ls)
  uh = solve(solver,op)

  tol = 1.0e-6
  e = uh-u(tf)
  l2(e) = e*e
  l2e = sqrt(sum( ∫(l2(e))dΩ ))
  # @test l2e < tol

  Ast = op.op.matrix
  bst = op.op.vector

  @test uh.free_values ≈ Ast \ bst

  return Ast, bst
end

A,rhs = extract_matrix_vector(a,fst)

gst(x) = u(tf)(x)
m(u,v) = ∫(u*v)dΩ

M,_ = extract_matrix_vector(m,gst)

@test rhs ≈ h
@test A+M/(θ*dt) ≈ K

rhs
h

end #module


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

mutable struct EncryptionService <: SwaggerModel
    enabled::Any # spec type: Union{ Nothing, Bool } # spec name: enabled
    lastEnabledTime::Any # spec type: Union{ Nothing, DateTime } # spec name: lastEnabledTime
    keyType::Any # spec type: Union{ Nothing, String } # spec name: keyType

    function EncryptionService(;enabled=nothing, lastEnabledTime=nothing, keyType=nothing)
        o = new()
        validate_property(EncryptionService, Symbol("enabled"), enabled)
        setfield!(o, Symbol("enabled"), enabled)
        validate_property(EncryptionService, Symbol("lastEnabledTime"), lastEnabledTime)
        setfield!(o, Symbol("lastEnabledTime"), lastEnabledTime)
        validate_property(EncryptionService, Symbol("keyType"), keyType)
        setfield!(o, Symbol("keyType"), keyType)
        o
    end
end # type EncryptionService

const _property_map_EncryptionService = Dict{Symbol,Symbol}(Symbol("enabled")=>Symbol("enabled"), Symbol("lastEnabledTime")=>Symbol("lastEnabledTime"), Symbol("keyType")=>Symbol("keyType"))
const _property_types_EncryptionService = Dict{Symbol,String}(Symbol("enabled")=>"Bool", Symbol("lastEnabledTime")=>"DateTime", Symbol("keyType")=>"String")
Base.propertynames(::Type{ EncryptionService }) = collect(keys(_property_map_EncryptionService))
Swagger.property_type(::Type{ EncryptionService }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_EncryptionService[name]))}
Swagger.field_name(::Type{ EncryptionService }, property_name::Symbol) =  _property_map_EncryptionService[property_name]

const _allowed_EncryptionService_keyType = ["Service", "Account"]

function check_required(o::EncryptionService)
    true
end

function validate_property(::Type{ EncryptionService }, name::Symbol, val)
    if name === Symbol("keyType")
        Swagger.validate_param(name, "EncryptionService", :enum, val, _allowed_EncryptionService_keyType)
    end
end


text:
# This file contains functions for the Reverse DDF DEA model
"""
    ReverseDDFDEAModel
An data structure representing a Reverse DDF DEA model.
"""
struct ReverseDDFDEAModel <: AbstractTechnicalDEAModel
    n::Int64
    m::Int64
    s::Int64
    measure::Symbol
    orient::Symbol
    rts::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
    Xtarget::Matrix
    Ytarget::Matrix
    Gxrddf::Union{Vector,Matrix}
    Gyrddf::Union{Vector,Matrix}
end

"""
    dearddf(X, Y, measure)
Compute data envelopment analysis reverse directional distance function (RDDF) model for inputs
`X`, outputs `Y`, and efficiency measure  `measure`.

# Measure specification:

- `:ERG`: Enhanced Russell Graph (or Slack Based Measure (SBM)).
- `:MDDF`: Modified Directional Distance Function.

# Direction specification:

For the Modified Directional Distance Function, the directions `Gx` and `Gy` can be one of the following symbols.
- `:Ones`: use ones.
- `:Observed`: use observed values.
- `:Mean`: use column means.

# Optional Arguments
- `orient=:Graph`: choose between graph oriented `:Graph`, input oriented `:Input`, or output oriented model `:Output`.
- `rts=:CRS`: choose between constant returns to scale `:CRS` or variable returns to scale `:VRS`.
- `atol=1e-6`: tolerance for DMU to be considered efficient.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [2; 4; 8; 12; 6; 14; 14; 9.412];

julia> Y = [1; 5; 8; 9; 3; 7; 9; 2.353] ;

julia> dearddf(X, Y, :ERG, orient = :Graph, rts = :VRS)
Reverse DDF DEA Model
DMUs = 8; Inputs = 1; Outputs = 1
Orientation = Graph; Returns to Scale = VRS
Associated efficiency measure = ERG
─────────────
   efficiency
─────────────
1    0.0
2    0.0
3    0.0
4    0.0
5    0.6
6    0.52381
7    0.142857
8    0.8
─────────────
```
"""
function dearddf(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector}, measure::Symbol;
    Gx::Union{Symbol,Matrix,Vector,Nothing} = nothing, Gy::Union{Symbol,Matrix,Vector,Nothing} = nothing,
    orient::Symbol = :Graph, rts::Symbol = :CRS, slack::Bool = false, atol::Float64 = 1e-6,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix,Vector,Nothing} = nothing,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::ReverseDDFDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nrefx != nrefy
        throw(DimensionMismatch("number of rows in Xref and Yref ($nrefx, $nrefy) are not equal"));
    end
    if m != mref
        throw(DimensionMismatch("number of columns in X and Xref ($m, $mref) are not equal"));
    end
    if s != sref
        throw(DimensionMismatch("number of columns in Y and Yref ($s, $sref) are not equal"));
    end

    n = nx

    # Get directions based on efficiency measure
    Gxrddf, Gyrddf = dearddfdirections(X, Y, measure, Gx = Gx, Gy = Gy, orient = orient, rts = rts, atol = atol, Xref = Xref, Yref = Yref, optimizer = optimizer)

    # Directional model using calculated directions
    rddfmodel = deaddf(X, Y, Gx = Gxrddf, Gy = Gyrddf, rts = rts, Xref = Xref, Yref = Yref, slack = slack, optimizer = optimizer)

    effi = efficiency(rddfmodel)
    lambdaeff = peersmatrix(rddfmodel)
    slackX = slacks(rddfmodel, :X)
    slackY = slacks(rddfmodel, :Y)
    Xtarget = targets(rddfmodel, :X)
    Ytarget = targets(rddfmodel, :Y)

    return ReverseDDFDEAModel(n, m, s, measure, orient, rts, names, effi, slackX, slackY, lambdaeff, Xtarget, Ytarget, Gxrddf, Gyrddf)

end

"""
    dearddfdirections(X, Y, measure)
Compute corresponding directions for data envelopment analysis Reverse DDF models
for inputs `X`, outputs `Y`, and efficiency `measure`.

# Measure specification:
- `:ERG`: Enhanced Russell Graph (or Slack Based Measure (SBM)).
- `:MDDF`: Modified Directional Distance Function.

# Direction specification:

For the Modified Directional Distance Function, the directions `Gx` and `Gy` can be one of the following symbols.
- `:Ones`: use ones.
- `:Observed`: use observed values.
- `:Mean`: use column means.

Alternatively, a vector or matrix with the desired directions can be supplied.
# Optional Arguments
- `orient=:Graph`: choose between graph oriented `:Graph`, input oriented `:Input`, or output oriented model `:Output`.
- `rts=:CRS`: choose between constant returns to scale `:CRS` or variable returns to scale `:VRS`.
- `atol=1e-6`: tolerance for DMU to be considered efficient.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
"""
function dearddfdirections(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector}, measure::Symbol;
    Gx::Union{Symbol,Matrix,Vector,Nothing} = nothing, Gy::Union{Symbol,Matrix,Vector,Nothing} = nothing,
    orient::Symbol = :Graph, rts::Symbol = :CRS, atol::Float64 = 1e-6,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix,Vector,Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)

    # Calculate direction based on efficiency measure
    Gxrddf = zeros(size(X))
    Gyrddf = zeros(size(Y))

    n = size(X, 1)

    if measure == :ERG
        # Show error if user specify directions
        if (Gx !== nothing) | (Gy !== nothing)
            throw(ArgumentError("`Gx` and `Gy` should not be specified for ERG measure."));
        end

        # Enhanced Russels Graph  (ERG) (or Slack Based Measure (SBM))
        if orient == :Graph
            ergmodel = deaerg(X, Y, rts = rts, Xref = Xref, Yref = Yref, optimizer = optimizer)
        elseif orient == :Input
            ergmodel = dearussell(X, Y, orient = :Input, rts = rts, Xref = Xref, Yref = Yref, slack = false, optimizer = optimizer)
        elseif orient == :Output
            ergmodel = dearussell(X, Y, orient = :Output, rts = rts, Xref = Xref, Yref = Yref, slack = false, optimizer = optimizer)
        else
            throw(ArgumentError("`orient` must be :Input or :Output"));
        end

        eff = efficiency(ergmodel)
        if orient == :Graph || orient == :Input
            ineff = 1 .- eff
        elseif orient == :Output
            #ineff = eff .- 1
            ineff = 1 .- (1 ./ eff)
        end

        Xtarget = targets(ergmodel, :X)
        Ytarget = targets(ergmodel, :Y)

        for i = 1:n
            if abs(ineff[i]) <= atol
                Gxrddf[i,:] .= 1
                Gyrddf[i,:] .= 1
            else
                Gxrddf[i,:] = (X[i,:] - Xtarget[i,:]) ./ ineff[i]
                Gyrddf[i,:] = (Ytarget[i,:] - Y[i,:]) ./ ineff[i]
            end
        end
    elseif measure == :MDDF
        # Modified DDF
        if orient != :Graph
            throw(ArgumentError("Measure MDDF only available for Graph orientation"))
        end

        # Build or get user directions
        if typeof(Gx) == Symbol
            if Gx == :Ones
                Gx = ones(size(X))
            elseif Gx == :Observed
                Gx = X
            elseif Gx == :Mean
                Gx = repeat(mean(X, dims = 1), size(X, 1))
            else
                throw(ArgumentError("Invalid `Gx`"));
            end
        end

        if typeof(Gy) == Symbol
            if Gy == :Ones
                Gy = ones(size(Y))
            elseif Gy == :Observed
                Gy = Y
            elseif Gy == :Mean
                Gy = repeat(mean(Y, dims = 1), size(Y, 1))
            else
                throw(ArgumentError("Invalid `Gy`"));
            end
        end

        mddfmodel = deamddf(X, Y, Gx = Gx, Gy = Gy, rts = rts, Xref = Xref, Yref = Yref, slack = false, optimizer = optimizer)

        betax = mddfmodel.betax
        betay = mddfmodel.betay

        for i = 1:n
            if betax[i] != 0
                Gxrddf[i,:] = (betax[i] / (betax[i] + betay[i])) .* Gx[i,:]
            else
                Gxrddf[i,:] .= 0
            end

            if betay[i] != 0
                Gyrddf[i,:] = (betay[i] / (betax[i] + betay[i])) .* Gy[i,:]
            else
                Gyrddf[i,:] .= 0
            end
        end

    else
        throw(ArgumentError("Invalid efficiency `measure`"));
    end

    # Set directions to 0 if Input or Output oriented
    if orient == :Input
        Gyrddf = zeros(size(Y))
    elseif orient == :Output
        Gxrddf = zeros(size(X))
    end

    return Gxrddf, Gyrddf

end

function Base.show(io::IO, x::ReverseDDFDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    eff = efficiency(x)
    dmunames = names(x)

    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)
    hasslacks = ! isempty(slackX)

    if !compact
        print(io, "Reverse DDF DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = ", string(x.orient))
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        print(io, "Associated efficiency measure = ", string(x.measure))
        print(io, "\n")

        if hasslacks == true
            show(io, CoefTable(hcat(eff, slackX, slackY), ["efficiency"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
        else
            show(io, CoefTable(hcat(eff), ["efficiency"], dmunames))
        end
    end

end


text:
mutable struct joint_state_t <: LCMType
    utime::Int64
    num_joints::Int16
    joint_name::Vector{String}
    joint_position::Vector{Float32}
    joint_velocity::Vector{Float32}
    joint_effort::Vector{Float32}
end

@lcmtypesetup(joint_state_t,
    joint_name => (num_joints,),
    joint_position => (num_joints,),
    joint_velocity => (num_joints,),
    joint_effort => (num_joints,),
)


text:
## Deprecated

# deprecated at 0.13
@deprecate kmpp(X, k) initseeds(:kmpp, X, k)
@deprecate kmpp_by_costs(costs, k) initseeds_by_costs(:kmpp, costs, k)

# deprecated at 0.13.1
@deprecate copyseeds(X, iseeds) copyseeds!(Matrix{eltype(X)}(undef, size(X, 1), length(iseeds)), X, iseeds)

# deprecated as of 0.13.2
@deprecate varinfo(k1::Int, a1::AbstractVector{Int},
                   k2::Int, a2::AbstractVector{Int}) varinfo(a1, a2)
@deprecate varinfo(R::ClusteringResult, k0::Int, a0::AbstractVector{Int}) varinfo(R, a0)

# FIXME remove after deprecation period for merge/labels/height/method
Base.propertynames(hclu::Hclust, private::Bool = false) =
    (:merges, :heights, :order, :linkage,
     #= deprecated as of 0.12 =# :height, :labels, :merge, :method)

# FIXME remove after deprecation period for merge/labels/height/method
@inline function Base.getproperty(hclu::Hclust, prop::Symbol)
    if prop === :height
        Base.depwarn("Hclust::height is deprecated, use Hclust::heights", Symbol("Hclust::height"))
        return getfield(hclu, :heights)
    elseif prop === :labels
        Base.depwarn("Hclust::labels is deprecated and will be removed in future versions", Symbol("Hclust::labels"))
        return 1:nnodes(hclu)
    elseif prop === :merge
        Base.depwarn("Hclust::merge is deprecated, use Hclust::merges", Symbol("Hclust::merge"))
        return getfield(hclu, :merges)
    elseif prop === :method
        Base.depwarn("Hclust::method is deprecated, use Hclust::linkage", Symbol("Hclust::method"))
        return getfield(hclu, :linkage)
    else
        return getfield(hclu, prop)
    end
end

# FIXME remove after deprecation period for cweights
Base.propertynames(clu::KmeansResult, private::Bool = false) =
    (:centers, :assignments, :costs, :counts, :wcounts,
     :totalcost, :iterations, :converged,
     #= deprecated as of 0.13.2 =# :cweights)

# FIXME remove after deprecation period for cweights
@inline function Base.getproperty(clu::KmeansResult, prop::Symbol)
    if prop === :cweights
        Base.depwarn("KmeansResult::cweights is deprecated, use wcounts(clu::KmeansResult)",
                     Symbol("KmeansResult::cweights"))
        return clu.wcounts
    else
        return getfield(clu, prop)
    end
end