text:
"""
`Manopt.jl` – Optimization on Manifolds in Julia.
"""
module Manopt
using Colors
using ColorSchemes
using ColorTypes
using Markdown
using LinearAlgebra
using Dates: Period, Nanosecond, value
using Random: shuffle!
using DataStructures: CircularBuffer, capacity, length, size, push!
using StaticArrays
import Random: rand, randperm
import Base: copy, identity
import ManifoldsBase:
    ℝ,
    ℂ,
    ×,
    ^,
    AbstractBasis,
    AbstractPowerManifold,
    AbstractVectorTransportMethod,
    AbstractRetractionMethod,
    AbstractInverseRetractionMethod,
    CachedBasis,
    DefaultOrthonormalBasis,
    ExponentialRetraction,
    LogarithmicInverseRetraction,
    ParallelTransport,
    PowerManifold,
    Manifold,
    allocate,
    allocate_result,
    allocate_result_type,
    distance,
    exp,
    exp!,
    log,
    log!,
    injectivity_radius,
    inner,
    geodesic,
    manifold_dimension,
    mid_point,
    mid_point!,
    norm,
    project,
    project!,
    retract,
    retract!,
    inverse_retract,
    inverse_retract!,
    shortest_geodesic,
    vector_transport_to,
    vector_transport_to!,
    zero_tangent_vector,
    zero_tangent_vector!,
    DiagonalizingOrthonormalBasis,
    get_basis,
    get_coordinates,
    get_vector,
    get_vectors,
    representation_size
using Manifolds: #temporary for random
    Circle,
    Euclidean,
    Grassmann,
    Hyperbolic,
    ProductManifold,
    Rotations,
    SymmetricPositiveDefinite,
    Stiefel,
    Sphere
import Manifolds: mid_point, mid_point!
using Manifolds: # Wishlist for Base
    NestedPowerRepresentation,
    ArrayPowerRepresentation,
    mean,
    median,
    get_iterator,
    _read,
    _write,
    power_dimensions,
    ArrayReshaper,
    prod_point,
    ShapeSpecification,
    ProductRepr,
    submanifold_components,
    submanifold_component,
    get_component,
    set_component!,
    getindex,
    setindex!

include("plans/plan.jl")
# Functions
include("functions/manifold.jl")
include("functions/bezier_curves.jl")
include("functions/adjoint_differentials.jl")
include("functions/costs.jl")
include("functions/differentials.jl")
include("functions/gradients.jl")
include("functions/Jacobi_fields.jl")
include("functions/proximal_maps.jl")
# solvers general framework
include("solvers/solver.jl")
# specific solvers
include("solvers/ChambollePock.jl")
include("solvers/conjugate_gradient_descent.jl")
include("solvers/cyclic_proximal_point.jl")
include("solvers/DouglasRachford.jl")
include("solvers/NelderMead.jl")
include("solvers/gradient_descent.jl")
include("solvers/particle_swarm.jl")
include("solvers/quasi_Newton.jl")
include("solvers/truncated_conjugate_gradient_descent.jl")
include("solvers/trust_regions.jl")
include("solvers/stochastic_gradient_descent.jl")
include("solvers/subgradient.jl")
include("solvers/debug_solver.jl")
include("solvers/record_solver.jl")
include("helpers/errorMeasures.jl")
include("helpers/exports/Asymptote.jl")
include("data/artificialDataFunctions.jl")

include("random.jl")

export ×, ^, ℝ, ℂ

export Problem,
    ProximalProblem,
    CostProblem,
    SubGradientProblem,
    GradientProblem,
    HessianProblem,
    PrimalDualProblem,
    StochasticGradientProblem

#
# Options
export Options,
    AbstractGradientDescentOptions,
    ChambollePockOptions,
    ConjugateGradientDescentOptions,
    CyclicProximalPointOptions,
    DouglasRachfordOptions,
    GradientDescentOptions,
    HessianOptions,
    NelderMeadOptions,
    ParticleSwarmOptions,
    PrimalDualOptions,
    StochasticGradientDescentOptions,
    SubGradientMethodOptions,
    TruncatedConjugateGradientOptions,
    TrustRegionsOptions
#
# Accessors and helpers for Options
export linesearch_backtrack
export get_cost,
    get_gradient,
    get_subgradient,
    get_proximal_map,
    get_options,
    get_initial_stepsize,
    get_gradients
export getHessian, approxHessianFD
export is_options_decorator, dispatch_options_decorator
export primal_residual, dual_residual

export QuasiNewtonOptions, QuasiNewtonLimitedMemoryDirectionUpdate
export QuasiNewtonCautiousDirectionUpdate,
    BFGS, InverseBFGS, DFP, InverseDFP, SR1, InverseSR1
export InverseBroyden, Broyden
export AbstractQuasiNewtonDirectionUpdate, AbstractQuasiNewtonUpdateRule
export WolfePowellLineseach,
    StrongWolfePowellLineseach,
    operator_to_matrix,
    square_matrix_vector_product,
    WolfePowellBinaryLinesearch

export ConjugateGradientDescentOptions,
    GradientDescentOptions,
    HessianOptions,
    SubGradientMethodOptions,
    NelderMeadOptions,
    TruncatedConjugateGradientOptions,
    TrustRegionsOptions,
    ParticleSwarmOptions
export AbstractOptionsAction, StoreOptionsAction
export has_storage, get_storage, update_storage!

#
# Direction Update Rules
export DirectionUpdateRule,
    IdentityUpdateRule, StochasticGradient, AverageGradient, MomentumGradient, Nesterov
export DirectionUpdateRule,
    SteepestDirectionUpdateRule,
    HeestenesStiefelCoefficient,
    FletcherReevesCoefficient,
    PolakRibiereCoefficient,
    ConjugateDescentCoefficient,
    LiuStoreyCoefficient,
    DaiYuanCoefficient,
    HagerZhangCoefficient
#
# Solvers
export ChambollePock,
    ChambollePock!,
    conjugate_gradient_descent,
    conjugate_gradient_descent!,
    cyclic_proximal_point,
    cyclic_proximal_point!,
    DouglasRachford,
    DouglasRachford!,
    gradient_descent,
    gradient_descent!,
    NelderMead,
    NelderMead!,
    particle_swarm,
    particle_swarm!,
    quasi_Newton,
    quasi_Newton!,
    stochastic_gradient_descent,
    stochastic_gradient_descent!,
    subgradient_method,
    subgradient_method!,
    truncated_conjugate_gradient_descent,
    truncated_conjugate_gradient_descent!,
    trust_regions,
    trust_regions!
# Solver helpers
export decorate_options
export initialize_solver!, step_solver!, get_solver_result, stop_solver!
export solve
#
# Stepsize
export ConstantStepsize, DecreasingStepsize
export Linesearch, ArmijoLinesearch, NonmonotoneLinesearch
export get_stepsize, get_initial_stepsize, get_last_stepsize
#
# Stopping Criteria
export StopIfResidualIsReducedByFactor,
    StopIfResidualIsReducedByPower,
    StopWhenCurvatureIsNegative,
    StopWhenTrustRegionIsExceeded
export StopAfterIteration, StopWhenChangeLess, StopWhenGradientNormLess, StopWhenCostLess
export StopAfter, StopWhenAll, StopWhenAny
export get_active_stopping_criteria, get_stopping_criteria, get_reason
export StoppingCriterion, StoppingCriterionSet, Stepsize

export random_point, random_tangent, mid_point, mid_point!, reflect, sym_rem
#
# Data functions
export artificial_S1_signal, artificial_S1_slope_signal, artificialIn_SAR_image
export artificial_SPD_image, artificial_SPD_image2
export artificial_S2_whirl_image, artificial_S2_whirl_patch
export artificial_S2_rotation_image
export artificial_S2_whirl_patch, artificial_S2_lemniscate
export artificial_S2_composite_bezier_curve
#
# Exports
export asymptote_export_S2_signals, asymptote_export_S2_data, asymptote_export_SPD
export render_asymptote
#
# Coeffs & Helpers for differentials
export βdifferential_geodesic_startpoint, βdifferential_exp_basepoint
export βdifferential_exp_argument, βdifferential_log_basepoint, βdifferential_log_argument
export jacobi_field, adjoint_Jacobi_field
#
# Adjoint differentials
export adjoint_differential_geodesic_startpoint, adjoint_differential_geodesic_endpoint
export adjoint_differential_exp_basepoint, adjoint_differential_exp_argument
export adjoint_differential_log_basepoint, adjoint_differential_log_argument
export adjoint_differential_forward_logs, adjoint_differential_bezier_control
#
# Differentials
export differential_geodesic_startpoint, differential_geodesic_endpoint
export differential_exp_basepoint, differential_exp_argument
export differential_log_basepoint, differential_log_argument, differential_forward_logs
export differential_bezier_control
#
# Functions
export costL2TV, costL2TVTV2, costL2TV2, costTV, costTV2, costIntrICTV12
export cost_L2_acceleration_bezier, cost_acceleration_bezier
# Gradients
export ∇TV, ∇TV2, ∇intrinsic_infimal_convolution_TV12, forward_logs, ∇distance
export ∇acceleration_bezier, ∇L2_acceleration_bezier
# Proximal maps
export prox_distance, prox_TV, prox_parallel_TV, prox_TV2, project_collaborative_TV
# Error measures
export meanSquaredError, meanAverageError
#
# Bézier
export BezierSegment,
    de_casteljau,
    get_bezier_degrees,
    get_bezier_degree,
    get_bezier_inner_points,
    get_bezier_junction_tangent_vectors,
    get_bezier_junctions,
    get_bezier_points,
    get_bezier_segments
#
# Debugs
export DebugOptions, DebugAction, DebugGroup, DebugEntry, DebugEntryChange, DebugEvery
export DebugChange, DebugIterate, DebugIteration, DebugDivider
export DebugCost, DebugStoppingCriterion, DebugFactory, DebugActionFactory
export DebugGradient, DebugGradientNorm, DebugStepsize
export DebugPrimalBaseChange, DebugPrimalBaseIterate, DebugPrimalChange, DebugPrimalIterate
export DebugDualBaseChange, DebugDualBaseIterate, DebugDualChange, DebugDualIterate
export DebugDualResidual, DebugPrimalDualResidual, DebugPrimalResidual
export DebugProximalParameter
export DebugGradient, DebugGradientNorm, DebugStepsize
#
# Records
export RecordGradient, RecordGradientNorm, RecordStepsize
export RecordPrimalBaseChange,
    RecordPrimalBaseIterate, RecordPrimalChange, RecordPrimalIterate
export RecordDualBaseChange, RecordDualBaseIterate, RecordDualChange, RecordDualIterate
export RecordProximalParameter
end


text:
using PublicSuffix
using Base.Test

psl = PublicSuffixList()

test_tld_valid = ["com", "co.uk", "in", "org"]
test_tld_invalid = ["uk", "sdfdsfsdfdsf.in"]

for tld in test_tld_valid
    println("testing existing tld: $tld")
    @test tld_exists(tld)
end

for tld in test_tld_invalid
    println("testing non existing tld: $tld")
    @test !tld_exists(tld)
end

d = Domain("www.google.com")
@test d.typ == Domain
@test d.full == "www.google.com"
@test d.sub_domain == "www"
@test d.public_suffix == "google.com"
@test d.top_domain == "com"

d = Domain("192.168.1.1")
@test d.typ == IPv4
@test d.full == "192.168.1.1"
@test isempty(d.sub_domain)
@test isempty(d.public_suffix)
@test isempty(d.top_domain)

d = Domain("2001:0db8:85a3:0000:0000:8a2e:0370:7334")
@test d.typ == IPv6
@test d.full == "2001:0db8:85a3:0000:0000:8a2e:0370:7334"
@test isempty(d.sub_domain)
@test isempty(d.public_suffix)
@test isempty(d.top_domain)

ext: jl


text:
module AStarSearch
export astar, AbstractAStarSearch, search

include("a_star.jl")

end # module


text:
function ClasspathCrawler(arg0::ClassLoader, arg1::Vector{JString})
    return ClasspathCrawler((ClassLoader, Vector{JString}), arg0, arg1)
end

function ClasspathCrawler(arg0::Vector{JString})
    return ClasspathCrawler((Vector{JString},), arg0)
end

function feed(obj::ClasspathCrawler, arg0::Pattern, arg1::DataLoader, arg2::DataProvidersManager)
    return jcall(obj, "feed", jboolean, (Pattern, DataLoader, DataProvidersManager), arg0, arg1, arg2)
end


text:
#s = open("/home/neomatrix/Escritorio/datos.txt", "a+")
# a+ es creacion, lectura y escritura
#for i=1:1000
#          write(s, AbstractString(rand(1701:254595)))
#          write(s, "\n")
#       end
#close(s)

include("persona.jl")
include("direccion.jl")
include("trabajo.jl")
include("datatime.jl")
include("phone_number.jl")
include("color.jl")
include("base.jl")
include("company.jl")
include("file.jl")
include("internet.jl")
include("lorem.jl")
include("misc.jl")
include("profile.jl")

t= time("%H-%M:%S")
d = date("%Y-%m:%d")
p = bothify("## ??")

for i=1:2
	println(" \n $(prefixe()) $(first_name()) $(last_name()) vive en $(state()) trabaja en $(job()) $(zip_code("OAX")) $(timezone()) $(am_pm())  $(day_of_week()) $(months()) $(month_name()) $(time()) $(t) $(date()) $(d) $(iso8601())$(date_time()) $(unix_time()) \n numero telefonico $(phone_number()) \n colores: $(color_name())  $(safe_color_name()) $(hex_color())  $(safe_hex_color())  $(rgb_color_list())  $(rgb_css_color()) \n$(street_name() ) \n $(p) \n $(company()) \n $(mime_type()) \n $(uri()) \n $(text(10)) \n $(password()) \n $(profile())")

end


text:
# exponential kernel with total variation distance
struct TVExponentialKernel <: KernelFunctions.SimpleKernel end

KernelFunctions.kappa(kernel::TVExponentialKernel, d) = exp(-d)
KernelFunctions.metric(::TVExponentialKernel) = TotalVariation()

# Simplifies computations of distances but is type piracy... :(
function KernelFunctions._scale(t::ScaleTransform, metric::TotalVariation, x, y)
    return first(t.s) * evaluate(metric, x, y)
end


text:
for ST in LazySets.subtypes(LazySet, true)
    # check that certain methods are implemented
    isconvextype(ST)
    isoperationtype(ST)
end


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

mutable struct BlobRestoreParameters <: SwaggerModel
    timeToRestore::Any # spec type: Union{ Nothing, DateTime } # spec name: timeToRestore
    blobRanges::Any # spec type: Union{ Nothing, Vector{BlobRestoreRange} } # spec name: blobRanges

    function BlobRestoreParameters(;timeToRestore=nothing, blobRanges=nothing)
        o = new()
        validate_property(BlobRestoreParameters, Symbol("timeToRestore"), timeToRestore)
        setfield!(o, Symbol("timeToRestore"), timeToRestore)
        validate_property(BlobRestoreParameters, Symbol("blobRanges"), blobRanges)
        setfield!(o, Symbol("blobRanges"), blobRanges)
        o
    end
end # type BlobRestoreParameters

const _property_map_BlobRestoreParameters = Dict{Symbol,Symbol}(Symbol("timeToRestore")=>Symbol("timeToRestore"), Symbol("blobRanges")=>Symbol("blobRanges"))
const _property_types_BlobRestoreParameters = Dict{Symbol,String}(Symbol("timeToRestore")=>"DateTime", Symbol("blobRanges")=>"Vector{BlobRestoreRange}")
Base.propertynames(::Type{ BlobRestoreParameters }) = collect(keys(_property_map_BlobRestoreParameters))
Swagger.property_type(::Type{ BlobRestoreParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_BlobRestoreParameters[name]))}
Swagger.field_name(::Type{ BlobRestoreParameters }, property_name::Symbol) =  _property_map_BlobRestoreParameters[property_name]

function check_required(o::BlobRestoreParameters)
    (getproperty(o, Symbol("timeToRestore")) === nothing) && (return false)
    (getproperty(o, Symbol("blobRanges")) === nothing) && (return false)
    true
end

function validate_property(::Type{ BlobRestoreParameters }, name::Symbol, val)
end


text:
module SecondFinDiffPlans

import ..PoincareInvariants: compute!, getpoints, getpointnum, getpointspec
import ..PoincareInvariants: SecondPoincareInvariant

using ..PoincareInvariants: @argcheck

using LinearAlgebra: dot

struct SecondFinDiffPlan{T, D} end

function SecondFinDiffPlan{T, D}(ω, ps::NTuple{2, Int}) where {T, D}
    SecondFinDiffPlan{T, D}()
end

function compute!(
    pinv::SecondPoincareInvariant{T, D, ωT, <:Any, P}, t, p
) where {T, D, ωT, P <: SecondFinDiffPlan}
    nx, ny = pinv.pointspec
    points = pinv.points

    @argcheck size(points) == (nx * ny, D) "Expected points mtarix to have size $((nx * ny, D))"

    I = zero(T)

    ∂xi = Vector{T}(undef, D)
    ∂yi = Vector{T}(undef, D)

    colviews = [view(points, :, d) for d in 1:D]

    i = 1
    for ix in 1:nx
        for iy in 1:ny
            for d in 1:D
                ∂xi[d], ∂yi[d] = differentiate(colviews[d], ix, iy, (nx, ny))
            end

            # This if statement should hopefully get optimised away by the compiler
            if ωT <: AbstractMatrix
                ωi = pinv.ω
            else
                pnti = view(points, i, :)
                ωi = pinv.ω(pnti, t, p)
            end

            w = getsimpweight(T, ix, iy, (nx, ny))

            I += w * dot(∂yi, ωi, ∂xi)

            i += 1
        end
    end

    return I
end

## getpoints, getpointnum and getpointspec ##

getpointnum((nx, ny)::NTuple{2, Integer}) = nx * ny

nextodd(n::Integer)::Int = n <= 3 ? 3 : 2 * fld(n, 2) + 1

function getpointspec(dims::NTuple{2, Integer}, ::Type{<:SecondFinDiffPlan})
    return (nextodd(dims[1]), nextodd(dims[2]))
end

function getpointspec(N::Integer, ::Type{<:SecondFinDiffPlan})
    n = nextodd(ceil(Int, sqrt(Int(N))))
    return (n, n)
end

function getpoints(f, ::Type{T}, dims::NTuple{2, Integer}, ::Type{<:SecondFinDiffPlan}) where T
    D = length(f(zero(T), zero(T)))
    nx, ny = dims
    N = nx * ny
    out = Matrix{T}(undef, N, D)

    i = 1
    for x in range(0, 1, length=nx)
        for y in range(0, 1, length=ny)
            out[i, :] .= f(x, y)
            i += 1
        end
    end

    # Should return vector instead of matrix in 1D case
    return D == 1 ? vec(out) : out
end

## Differentiation ##

function _getmid(ix, iy, nx, ny)
    x0, y0 = ix, iy
    edge = false

    if ix == 1
        x0 = 2
        edge = true
    end

    if iy == 1
        y0 = 2
        edge = true
    end

    if ix == nx
        x0 = nx - 1
        edge = true
    end

    if iy == ny
        y0 = ny - 1
        edge = true
    end

    return x0, y0, edge
end

# evaluate derivative of quadratic approxmiation around (x0, y0) at (x0 + x, y0 + y)
# We have approximation f(x, y) = c + cx*x + cy*y + cxx*x^2 + cxy*x*y + cyy*y^2
function differentiate(
    vals::AbstractVector{T}, ix::Integer, iy::Integer, (nx, ny)::NTuple{2, Integer}
) where T
    x0, y0, edge = _getmid(ix, iy, nx, ny)

    mat = reshape(vals, ny, nx)
    f(x, y) = @inbounds mat[y0 + y, x0 + x]

    # All coefficients are multiplied by 2*Δx or 2*Δy
    bx = f(1, 0); ax = f(-1, 0)
    fx = bx - ax  # cx

    by = f(0, 1); ay = f(0, -1)
    fy = by - ay  # cy

    if edge
        f00 = f(0, 0)
        cxx = bx + ax - 2*f00
        cyy = by + ay - 2*f00
        cxy = (f(1, 1) + f(-1, -1) - f(1, -1) - f(-1, 1)) / 2

        dx, dy = ix - x0, iy - y0
        fx += 2*cxx * dx + cxy * dy
        fy += 2*cyy * dy + cxy * dx
    end

    invΔx = nx - 1; invΔy = ny - 1
    return fx * invΔx / 2, fy * invΔy / 2
end

## Integration ##

function _sw(i, n)
    if i == 1 || i == n
        return 1
    elseif iseven(i)
        return 4
    else
        @assert isodd(i)
        return 2
    end
end

function getsimpweight(::Type{T}, x, y, (nx, ny)) where T
    return T(_sw(x, nx) * _sw(y, ny)) / (9 * (nx - 1) * (ny - 1))
end

end  # SecondFinDiffPlans


text:
"""
Matrix of a 2D piecewise linear FEM discretization
=======================

*Input options:*

+ n: size of the matrix
"""
function fem2d(n::Int; nnzrow=4)
    A = ExtendableSparseMatrix(n, n)
    sprand_sdd!(A, nnzrow=nnzrow)
    flush!(A)
    return A.cscmatrix
end

"""
Matrix for a mock finite difference operator for a diffusion
problem with random coefficients on a unit hypercube
=======================

*Input options:*

+ n: number of unknowns in x, y, z direction
"""
function fdm3d(n::Int)
	nx = floor(Int, n^(1/3))
	ny = floor(Int, n^(1/3))
	nz = floor(Int, n^(1/3))
    A = ExtendableSparseMatrix(nx*ny*nz, nx*ny*nz)
    fdrand!(A, nx, ny, nz)
    flush!(A)
    return  A.cscmatrix
end

max_stars_repo_name: jkrch/JuliaParGMRES-BenchCLI

max_issues_repo_name: jkrch/JuliaParGMRES-BenchCLI

max_forks_repo_name: jkrch/JuliaParGMRES-BenchCLI


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

mutable struct ApplicationGatewayPrivateLinkIpConfigurationProperties <: SwaggerModel
    privateIPAddress::Any # spec type: Union{ Nothing, String } # spec name: privateIPAddress
    privateIPAllocationMethod::Any # spec type: Union{ Nothing, IPAllocationMethod } # spec name: privateIPAllocationMethod
    subnet::Any # spec type: Union{ Nothing, SubResource } # spec name: subnet
    primary::Any # spec type: Union{ Nothing, Bool } # spec name: primary
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function ApplicationGatewayPrivateLinkIpConfigurationProperties(;privateIPAddress=nothing, privateIPAllocationMethod=nothing, subnet=nothing, primary=nothing, provisioningState=nothing)
        o = new()
        validate_property(ApplicationGatewayPrivateLinkIpConfigurationProperties, Symbol("privateIPAddress"), privateIPAddress)
        setfield!(o, Symbol("privateIPAddress"), privateIPAddress)
        validate_property(ApplicationGatewayPrivateLinkIpConfigurationProperties, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        setfield!(o, Symbol("privateIPAllocationMethod"), privateIPAllocationMethod)
        validate_property(ApplicationGatewayPrivateLinkIpConfigurationProperties, Symbol("subnet"), subnet)
        setfield!(o, Symbol("subnet"), subnet)
        validate_property(ApplicationGatewayPrivateLinkIpConfigurationProperties, Symbol("primary"), primary)
        setfield!(o, Symbol("primary"), primary)
        validate_property(ApplicationGatewayPrivateLinkIpConfigurationProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ApplicationGatewayPrivateLinkIpConfigurationProperties

const _property_map_ApplicationGatewayPrivateLinkIpConfigurationProperties = Dict{Symbol,Symbol}(Symbol("privateIPAddress")=>Symbol("privateIPAddress"), Symbol("privateIPAllocationMethod")=>Symbol("privateIPAllocationMethod"), Symbol("subnet")=>Symbol("subnet"), Symbol("primary")=>Symbol("primary"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ApplicationGatewayPrivateLinkIpConfigurationProperties = Dict{Symbol,String}(Symbol("privateIPAddress")=>"String", Symbol("privateIPAllocationMethod")=>"IPAllocationMethod", Symbol("subnet")=>"SubResource", Symbol("primary")=>"Bool", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ ApplicationGatewayPrivateLinkIpConfigurationProperties }) = collect(keys(_property_map_ApplicationGatewayPrivateLinkIpConfigurationProperties))
Swagger.property_type(::Type{ ApplicationGatewayPrivateLinkIpConfigurationProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayPrivateLinkIpConfigurationProperties[name]))}
Swagger.field_name(::Type{ ApplicationGatewayPrivateLinkIpConfigurationProperties }, property_name::Symbol) =  _property_map_ApplicationGatewayPrivateLinkIpConfigurationProperties[property_name]

function check_required(o::ApplicationGatewayPrivateLinkIpConfigurationProperties)
    true
end

function validate_property(::Type{ ApplicationGatewayPrivateLinkIpConfigurationProperties }, name::Symbol, val)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Libdl, MbedTLS_jll

@testset "MbedTLS_jll" begin
    vstr = zeros(UInt8, 32)
    ccall((:mbedtls_version_get_string, libmbedcrypto), Cvoid, (Ref{UInt8},), vstr)
    vn = VersionNumber(unsafe_string(pointer(vstr)))
    @test vn == v"2.24.0"
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
using LinearAlgebra
using Printf
using BenchmarkTools

mutable struct LinOp
    matvec
    dim::Int
end

Base.size(lop::LinOp) = return (lop.dim,lop.dim)
Base.:(*)(lop::LinOp, v::AbstractVector{T}) where {T} = return lop.matvec(v)
Base.:(*)(lop::LinOp, v::AbstractMatrix{T}) where {T} = return lop.matvec(v)

mutable struct Davidson
    op
    dim::Int
    nroots::Int
    max_iter::Int
    max_ss_vecs::Int
    tol::Float64
    converged::Bool
    status::Vector{Bool}        # converged status of each root
    iter_curr::Int
    vec_curr::Array{Float64,2}
    sig_curr::Array{Float64,2}
    vec_prev::Array{Float64,2}
    sig_prev::Array{Float64,2}
    ritz_v::Array{Float64,2}
    ritz_e::Vector{Float64}
    resid::Vector{Float64}
end

function Davidson(op; max_iter=100, max_ss_vecs=8, tol=1e-8, nroots=1, v0=nothing)
    size(op)[1] == size(op)[2] || throw(DimensionError())
    dim = size(op)[1]
    if v0==nothing
        v0 = rand(dim,nroots)
        S = v0'*v0
        v0 = v0*inv(sqrt(S))
    end
    #display(v0'*v0)
    return Davidson(op,
                    dim,
                    nroots,
                    max_iter,
                    max_ss_vecs*nroots,
                    tol,
                    false,
                    [false for i in 1:nroots],
                    0,
                    v0,
                    v0,
                    zeros(dim,0),
                    zeros(dim,0),
                    zeros(nroots,nroots),
                    zeros(nroots),
                    zeros(nroots))
end

function print_iter(solver::Davidson)
    @printf(" Iter: %3i SS: %-4i", solver.iter_curr, size(solver.vec_prev)[2])
    @printf(" E: ")
    for i in 1:solver.nroots
        if solver.status[i]
            @printf("%13.8f* ", solver.ritz_e[i])
        else
            @printf("%13.8f  ", solver.ritz_e[i])
        end
    end
    @printf(" R: ")
    for i in 1:solver.nroots
        if solver.status[i]
            @printf("%5.1e* ", solver.resid[i])
        else
            @printf("%5.1e  ", solver.resid[i])
        end
    end
    println("")
    flush(stdout)
end

function _apply_diagonal_precond!(res_s::Vector{Float64}, Adiag::Vector{Float64}, denom::Float64)
    dim = size(Adiag)[1]
    length(size(Adiag)) == 1 || throw(Exception)
    length(size(res_s)) == 1 || throw(Exception)
    @inbounds @simd for i in 1:dim
        res_s[i] = res_s[i] / (denom - Adiag[i])
    end
end

function iteration(solver::Davidson; Adiag=nothing, iprint=0)
    #@printf(" Davidson Iter: %4i SS_Curr %4i\n", iter, size(solver.vec_curr,2))
    #
    # perform sig_curr = A*vec_curr
    solver.sig_curr = Matrix(solver.op * solver.vec_curr)
    #
    # add these new vectors to previous quantities
    solver.sig_prev = hcat(solver.sig_prev, solver.sig_curr)
    solver.vec_prev = hcat(solver.vec_prev, solver.vec_curr)

    #
    # form H in current subspace
    Hss = solver.vec_prev' * solver.sig_prev
    F = eigen(Hss)
    #ritz_e = F.values[1:solver.nroots]
    #ritz_v = F.vectors[:,1:solver.nroots]
    ritz_e = F.values
    ritz_v = F.vectors
    ss_size = length(F.values)
    if ss_size >= solver.max_ss_vecs
        ritz_v = ritz_v[:,sortperm(ritz_e)][:,1:solver.nroots]
        ritz_e = ritz_e[sortperm(ritz_e)][1:1:solver.nroots]
    else
        ritz_v = ritz_v[:,sortperm(ritz_e)]
        ritz_e = ritz_e[sortperm(ritz_e)]
    end
    solver.ritz_e = ritz_e
    solver.ritz_v = ritz_v

    solver.sig_prev = solver.sig_prev * ritz_v
    solver.vec_prev = solver.vec_prev * ritz_v
    Hss1 = solver.vec_prev' * solver.sig_prev
    Hss = ritz_v' * Hss * ritz_v
    # make sure these match
    all([abs(Hss1[i]-Hss[i])<1e-12 for i in 1:Base.length(Hss)]) || throw(Exception)

    if iprint>0
        println(" Determinant of subspace overlap: ",det(solver.vec_prev'*solver.vec_prev))
    end
    #Hss = solver.vec_curr' * sigma
    if iprint>0
        [@printf(" Ritz Value: %12.8f \n",i) for i in ritz_e]
    end
    res = similar(solver.vec_prev[:,1:solver.nroots])
    for s in 1:solver.nroots
        res[:,s] .= -solver.vec_prev[:,s] * Hss[s,s]
    end
    res[:,:] .+= solver.sig_prev[:,1:solver.nroots]

    #@btime $res = - $solver.vec_prev * $Hss[1:$solver.nroots,1:$solver.nroots]
    #@btime $res .+= $solver.sig_prev
    #res = solver.sig_prev - (solver.vec_prev * Hss)
    #@btime $res = $solver.sig_prev - ($solver.vec_prev * $Hss)
    ss = deepcopy(solver.vec_prev)

    new = zeros(solver.dim,0)
    #solver.statusconv = [false for i in 1:solver.nroots]
    for s in 1:solver.nroots
        # |v'> = (I-|SS><SS|)|v>
        #      = |v> - |SS><SS|V>
        solver.resid[s] = norm(res[:,s])
        if norm(res[:,s]) <= solver.tol
            solver.status[s] = true
            continue
        end
        if Adiag != nothing
            level_shift = 1e-3

            denom = (ritz_e[s]  + level_shift)
            _apply_diagonal_precond!(res[:,s], Adiag, denom)
            #@btime $_apply_diagonal_precond!($res[:,$s], $Adiag, $denom)
            #@inbounds @simd for i in 1:solver.dim
            #    res_s[i] = res_s[i] / (denom - Adiag[i])
            #end
        end
        #scr = solver.vec_prev' * res[:,s]
        #res[:,s] = res[:,s] - (solver.vec_prev * scr)
        scr = ss' * res[:,s]
        res[:,s] = res[:,s] - (ss * scr)
        nres = norm(res[:,s])
        if nres>1e-12
            ss = hcat(ss,res[:,s]/nres)
            new = hcat(new,res[:,s]/nres)
        end
    end
    return new
end

function solve(solver::Davidson; Adiag=nothing, iprint=0)

    for iter = 1:solver.max_iter
        #@btime $solver.vec_curr = $iteration($solver)
        solver.vec_curr = iteration(solver, Adiag=Adiag, iprint=iprint)
        solver.iter_curr = iter
        print_iter(solver)
        if all(solver.status)
            solver.converged == true
            break
        end
    end
    return solver.ritz_e[1:solver.nroots], solver.vec_prev[:,1:solver.nroots]
    #return solver.ritz_e[1:solver.nroots], solver.vec_prev*solver.ritz_v[:,1:solver.nroots]
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
export check_tri_index, car2lin_unchecked, car2lin, lin2car

# Convert between 1D and 2D indexing for triangular arrays

"""
    check_tri_index(Bool, layout::TriLayout, i, j)::Bool
    check_tri_index(Bool, layout::TriLayout, idx::CartesianIndex{2})::Bool

Check whether a row/column index is within the stored region of a `TriMatrix`
with the given layout. Assumes the index is otherwise valid for the
corresponding matrix size.
"""
function check_tri_index end

@inline check_tri_index(::Type{Bool}, ::TriLower{true}, i, j) = j <= i
@inline check_tri_index(::Type{Bool}, ::TriLower{false}, i, j) = j < i
@inline check_tri_index(::Type{Bool}, ::TriUpper{true}, i, j) = i <= j
@inline check_tri_index(::Type{Bool}, ::TriUpper{false}, i, j) = i < j
@inline check_tri_index(::Type{Bool}, ::TriSymmetric{true}, i, j) = true
@inline check_tri_index(::Type{Bool}, ::TriSymmetric{false}, i, j) = i != j
check_tri_index(::Type{Bool}, layout::TriLayout, idx::CartesianIndex{2}) = check_tri_index(Bool, layout, idx[1], idx[2])

"""
    check_tri_index(layout::TriLayout, i::Integer, j::Integer)
    check_tri_index(layout::TriLayout, idx::CartesianIndex{2})

Check that a row/column index is within the stored data region of a `TriMatrix`
with the given `layout` or throw a `DomainError`. Assumes the index is otherwise
valid for the corresponding matrix size.
"""
function check_tri_index(layout::TriLayout, i::Integer, j::Integer)
    check_tri_index(Bool, layout, i, j) || throw(DomainError("Invalid index for TriLayout $(layout): ($i, $j)"))
end

check_tri_index(layout::TriLayout, idx::CartesianIndex{2}) = check_tri_index(layout, idx[1], idx[2])

# Convert cartesian index of given layout to cartesian index of TriLower{true}
# which has the same index in the data array.
# Assumes index is within the stored data region.
@inline cartesian_to_tril(::TriLower{true}, i, j) = (i, j)
@inline cartesian_to_tril(::TriLower{false}, i, j) = (i - 1, j)
@inline cartesian_to_tril(::TriUpper{true}, i, j) = (j, i)
@inline cartesian_to_tril(::TriUpper{false}, i, j) = (j - 1, i)
@inline function cartesian_to_tril(::TriSymmetric{D}, i, j) where D
	if i >= j
		return cartesian_to_tril(TriLower{D}(), i, j)
	else
		return cartesian_to_tril(TriUpper{D}(), i, j)
	end
end

# Inverse of cartesian_to_tril()
@inline cartesian_from_tril(::TriLower{true}, i, j) = (i, j)
@inline cartesian_from_tril(::TriLower{false}, i, j) = (i + 1, j)
@inline cartesian_from_tril(::TriUpper{true}, i, j) = (j, i)
@inline cartesian_from_tril(::TriUpper{false}, i, j) = (j, i + 1)
@inline cartesian_from_tril(::TriSymmetric{D}, i, j) where D = cartesian_from_tril(TriLower{D}(), i, j)

"""
	$(FUNCTIONNAME)(layout::TriLayout, i::Integer, j::Integer)
	$(FUNCTIONNAME)(layout::TriLayout, idx::CartesianIndex{2})

Like [`car2lin`](@ref), but somewhat faster because it does not check that the
index is within the data region of the matrix.
"""
@inline car2lin_unchecked(::TriLower{true}, i::Integer, j::Integer) = trinum(i - 1) + j
@inline car2lin_unchecked(layout::TriLayout, i::Integer, j::Integer) = car2lin_unchecked(TriLower(), cartesian_to_tril(layout, i, j)...)
# @inline car2lin_unchecked(::TriLower{true}, i, j) = car2lin_unchecked(i, j)
# @inline car2lin_unchecked(::TriLower{false}, i, j) = car2lin_unchecked(i - 1, j)
# @inline car2lin_unchecked(::TriUpper{D}, i, j) where D = car2lin_unchecked(TriLower{D}, j, i)
# function @inline car2lin_unchecked(::TriSymmetric{D}, i, j) where D
#     if i >= j
#         return car2lin_unchecked(TriLower{D}(), i, j)
#     else
#         return car2lin_unchecked(TriUpper{D}(), j, i)
#     end
# end
@inline car2lin_unchecked(layout::TriLayout, idx::CartesianIndex{2}) = car2lin_unchecked(layout, idx[1], idx[2])

"""
	$(FUNCTIONNAME)(layout::TriLayout, i::Integer, j::Integer)
	$(FUNCTIONNAME)(layout::TriLayout, idx::CartesianIndex{2})

Convert Cartesian row/column indices of a `TriMatrix` with the given layout to
a linear index of its data array.

This function is the inverse of [`lin2car`](@ref).
"""
function car2lin(layout::TriLayout, i::Integer, j::Integer)
	check_tri_index(layout, i, j)
	return car2lin_unchecked(layout, i, j)
end

car2lin(layout::TriLayout, idx::CartesianIndex{2}) = car2lin(layout, idx[1], idx[2])

"""
    $(FUNCTIONNAME)(layout::TriLayout, i::Integer)

Convert linear index of the data array of a `TriMatrix` with the given layout to
the corresponding Cartesian index.

This function is the inverse of [`car2lin`](@ref).
"""
function lin2car(layout::TriLayout, i::Integer)
	r, c = lin2car(i)
	return cartesian_from_tril(layout, r, c)
end

function lin2car(i::Integer)
	n, i = triinv_rem(i - 1)
	return (n + 1, i + 1)
end


text:
"""
Here we define the interface that all contexts implement and provide a simple implementation

Constructor takes compute graph and any implementation specific parameters and initialises initial
tensor and parameter storage

Each context implements the following functions to access tensor and parameter information:
gettensor(ctx, sym): get the tensor for the given symbol
settensor!(ctx, value, sym): set the tensor data for the given symbol
deletetensor!(ctx, sym): delete the tensor after it is no longer required
Base.getindex(ctx, sym): get the parameter value corresponding to given symbol
Base.setindex!(ctx, value, sym): set the parameter value corresponding to given symbol
Base.haskey(ctx, sym): Check if the parameter key exists
Base.zeros(ctx, size): create array of zeros with appropritate type for context
Base.zero(ctx): create scalar with same type used in context
Base.eltype(ctx): return the element type of numeric datastructures
set_open_bonds!(ctx, bitstring::String): Set output parameters according to provided bitstring
set_slice_vals!(ctx, slice_values::Vector{Int}): Set output parameters according to provided slice values
(c::ctx)(): Execute the compute graph in the provided context. Returns final tensor
compute_amplitude!(ctx, bitstring; num_slices=nothing): Contracts the network to compute the given bitstring. Includes reduction over slices

Also provides an implementation of each command which uses above functions to get required
tensors and parameters. For example contraction command this function is defined as

(c::ContractionCommand)(ctx): implenents the contraction described by the given contraction index with the context

Each context will also implement the following functions to be implemented by distributed contexts
ctxmap(f, ctx, items): Applies function f to items
ctxreduce(f, ctx, items): Performs reduction using the function f on items
ctxgather(ctx, items): Gathers all items
"""

using QXContexts.ComputeGraphs
using OMEinsum
using DataStructures
using CUDA

export gettensor, settensor!, deletetensor!, set_open_bonds!, set_slice_vals!
export AbstractContext, QXContext, compute_amplitude!
export ctxmap, ctxgather, ctxreduce

abstract type AbstractContext end

##################################################################################
# Provide implementation of each command
##################################################################################
"""
    (c::ContractCommand)(ctx::AbstractContext)

Execute a contraction command in the given context
"""
function (c::ContractCommand)(ctx::AbstractContext)
    output_idxs = Tuple(c.output_idxs)
    left_idxs = Tuple(c.left_idxs)
    right_idxs = Tuple(c.right_idxs)

    @debug "ncon DSL command: $(c.output_name)[$(output_idxs)] = $(c.left_name)[$(c.left_idxs)] * $(c.right_name)[$(c.right_idxs)]"
    @debug "ncon shapes: left_size=$(size(gettensor(ctx, c.left_name))), right_size=$(size(gettensor(ctx, c.right_name)))"
    @nvtx_range "NCON $(c.output_name)" begin
        settensor!(ctx, EinCode((left_idxs, right_idxs), output_idxs)(gettensor(ctx, c.left_name), gettensor(ctx, c.right_name)), c.output_name)
    end
    deletetensor!(ctx, c.left_name)
    deletetensor!(ctx, c.right_name)
    nothing
end

"""
    (c::LoadCommand)(ctx::AbstractContext)

Generic implementation of load command
"""
function (c::LoadCommand)(ctx::AbstractContext)
    @nvtx_range "Load $(c.label)" begin
        settensor!(ctx, gettensor(ctx, c.label), c.name)
    end
end

"""
    (c::SaveCommand)(ctx::AbstractContext)

Generic implementation of save command
"""
function (c::SaveCommand)(ctx::AbstractContext)
    @nvtx_range "Load $(c.label)" begin
        settensor!(ctx, gettensor(ctx, c.label), c.name)
    end
end

"""
    (c::ReshapeCommand)(ctx::AbstractContext)

Implementation of reshape command
"""
function (c::ReshapeCommand)(ctx::AbstractContext)
    tensor_dims = size(gettensor(ctx, c.input))
    new_dims = [prod([tensor_dims[y] for y in x]) for x in c.dims]
    @nvtx_range "Reshape $(c.output)" begin
        settensor!(ctx, reshape(gettensor(ctx, c.input), new_dims...), c.output)
    end
    @debug "Reshape DSL command: name=$(c.input)($(tensor_dims)) -> $(c.output)($(new_dims))"
    nothing
end

"""
    (c::ViewCommand)(ctx::AbstractContext)

Execute the given view command using provided context
"""
function (c::ViewCommand)(ctx::AbstractContext)
    bond_val = haskey(ctx, c.slice_sym) ? ctx[c.slice_sym] : nothing
    @nvtx_range "View $(c.output_sym)" begin
        if bond_val !== nothing
            dims = size(gettensor(ctx, c.input_sym))
            view_index_list = [i == c.bond_index ? UnitRange(bond_val, bond_val) : UnitRange(1, dims[i]) for i in 1:length(dims)]
            new_tensor = @view gettensor(ctx, c.input_sym)[view_index_list...]
            settensor!(ctx, new_tensor, c.output_sym)
            @debug "view DSL command: $(c.output_sym) = $(c.input_sym)[$(view_index_list)]"
        else
            settensor!(ctx, gettensor(ctx, c.input_sym), c.output_sym)
            @debug "view DSL command: $(c.output_sym) = $(c.input_sym)"
        end
    end
    nothing
end

"""
    (c::OutputCommand)(ctx::AbstractContext)

Execute the given output command using provided context
"""
function (c::OutputCommand)(ctx::AbstractContext)
    @nvtx_range "Output $(c.idx)" begin
        sym = Symbol("o$(c.idx)")
        @assert haskey(ctx, sym) "Output $sym not set in context"
        out_val = ctx[sym]
        settensor!(ctx, gettensor(ctx, Symbol("output_$out_val")), c.name)
        @debug "output DSL command: $(c.name) = o$(c.idx)"
    end
end

##########################################################################################
# Provide concreate implementation of a context
#########################################################################################

"""Data structure for context implementation"""
struct QXContext{T} <: AbstractContext
    params::Dict{Symbol, Int}
    tensors::Dict{Symbol, T}
    cg::ComputeGraph
    slice_dims::OrderedDict{Symbol, Int}
    output_dims::OrderedDict{Symbol, Int}
end

"""
    QXContext{T}(cg::ComputeGraph) where T

Constuctor which initialises instance from a compute graph
"""
function QXContext{T}(cg::ComputeGraph) where T
    tensors = Dict{Symbol, T}()
    # load tensors from compute graph into this dictionary
    for (k, v) in pairs(cg.tensors)
        tensors[k] = convert(T, v)
    end
    slice_dims = convert(OrderedDict, params(cg, ViewCommand))
    output_dims = convert(OrderedDict, params(cg, OutputCommand))
    dims = collect(values(output_dims))
    if length(dims) > 0
        @assert all(dims .== dims[1]) "Multiple output dimensions not supported"
        for d in 1:dims[1]
            t = zeros(eltype(T), dims[1])
            t[d] = 1.
            tensors[Symbol("output_$(d-1)")] = convert(T, t)
        end
    end
    sort!(slice_dims)
    sort!(output_dims)
    QXContext{T}(Dict{Symbol, Int}(), tensors, cg, slice_dims, output_dims)
end

QXContext(cg::ComputeGraph) = QXContext{Array{ComplexF32}}(cg)

"""
    gettensor(ctx::QXContext, sym)

Function to retrieve tensors by key
"""
function gettensor(ctx::QXContext, sym)
    @assert haskey(ctx.tensors, sym) "Tensor $sym does not exist in this context"
    ctx.tensors[sym]
end

"""
    settensor!(ctx::QXContext, sym)

Function to set tensors by key
"""
function settensor!(ctx::QXContext, value, sym)
    ctx.tensors[sym] = value
end

"""
    deletetensor!(ctx::QXContext, sym)
Function to delte tensors by key
"""
function deletetensor!(ctx::QXContext, sym)
    delete!(ctx.tensors, sym)
end

"""Implement has key to check if parameter by this name present"""
Base.haskey(ctx::QXContext, sym) = haskey(ctx.params, sym)

"""
    Base.getindex(ctx::QXContext, sym)::Int

Implement getindex for retrieving parameter values by key
"""
function Base.getindex(ctx::QXContext, sym)::Int
    @assert haskey(ctx, sym) "Parameter $sym does not exist in this context"
    ctx.params[sym]
end

"""
    Base.setindex!(ctx::QXContext, sym)::Int

Implement setindex! for setting parameter values
"""
function Base.setindex!(ctx::QXContext, value, sym)
    ctx.params[sym] = value
end

Base.zeros(::QXContext{T}, size) where T = convert(T, zeros(eltype(T), size))
Base.zero(::QXContext{T}) where T = zero(eltype(T))
Base.eltype(::QXContext{T}) where T = eltype(T)

"""
    set_open_bonds!(ctx::QXContext, bitstring::String)

Given a bitstring, set the open bonds to values so contracting the network will
calculate the amplitude of this bitstring
"""
function set_open_bonds!(ctx::QXContext, bitstring::String="")
    if bitstring == "" bitstring = "0"^length(ctx.output_dims) end
    @assert length(bitstring) == length(ctx.output_dims) "Bitstring length must match nubmer of outputs"
    for (i, key) in enumerate(keys(ctx.output_dims)) ctx[key] = parse(Int, bitstring[i]) end
end

"""
    set_slice_vals!(ctx::QXContext, slice_values::Vector{Int})

For each bond that is being sliced set the dimension to slice on.
"""
function set_slice_vals!(ctx::QXContext, slice_values::Vector{Int})
    # set slice values and remove any already set
    for (i, key) in enumerate(keys(ctx.slice_dims))
        if i > length(slice_values)
            delete!(ctx.params, key)
        else
            ctx[key] = slice_values[i]
        end
    end
end

"""
    (ctx::QXContext)()

Funciton to execute compute graph when struct is called. Returns final tensor
"""
function (ctx::QXContext)()
    for n in get_commands(ctx.cg) n(ctx) end
    gettensor(ctx, output(ctx.cg))
end

"""
    compute_amplitude!(ctx::QXContext, bitstring::String; max_slices=nothing)

Calculate a single amplitude with the given context and bitstring. Involves a sum over
contributions from each slice. Can optionally set the number of bonds. By default all slices
are used.
"""
function compute_amplitude!(ctx::QXContext, bitstring::String; max_slices=nothing)
    set_open_bonds!(ctx, bitstring)
    amplitude = nothing
    for p in SliceIterator(ctx.cg, max_slices=max_slices)
        set_slice_vals!(ctx, p)
        if amplitude === nothing
            amplitude = ctx()
        else
            amplitude += ctx()
        end
    end
    amplitude = convert(Array, amplitude) # if a gpu array convert back to gpu
    if ndims(amplitude) == 0
        amplitude = amplitude[]
    end
    amplitude
end

"""Map over items as placeholder for more complicated contexts"""
ctxmap(f, ctx::QXContext, items) = map(f, items)

"""Simple gather as placeholder for distributed contexts"""
ctxgather(ctx::QXContext, items) = items

"""Simple gather as placeholder for distributed contexts"""
ctxreduce(f, ctx::QXContext, items) = reduce(f, items)


text:
using GPFlux
using Test, Random
using Flux

@testset "model_gradient_check" begin
	Random.seed!(4)
	T = Float64
	X = rand(T, 4, 10)
	y = rand(T, 10)

	# kernels
	rbf_kernel = GaussKernel(4)
	isorbf_kernel = IsoGaussKernel()
	rq_kernel = RQKernel(4)
	isorq_kernel = IsoRQKernel()
	isoperiod_kernel = IsoPeriodKernel()
	linear_kernel = LinearKernel(4)
	kernels = [rbf_kernel, isorbf_kernel, rq_kernel, isorq_kernel, isoperiod_kernel, linear_kernel]

	# means
	c_mean = ConstantMean()
	simple_nn_mean = SimpleNeuralNetworkMean(4, 1)
	means = [c_mean, simple_nn_mean]

	for kernel in kernels, mean in means
		mean = reset(mean)
		kernel = reset(kernel)
		@show (mean, kernel)
		gp = GaussProcess(mean, kernel)
		ps = params(gp)
		f = () -> negloglik(gp, X, y)
		@test model_gradient_check(f, ps)
	end
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
# sta gives the Spike Triggered Average (average neuron activity preceding a spike)
# stim is the recording of voltage over time in the neuron
# ρ is a matching time series which records the neuron's spikes
# num_timesteps gives the desired time period to record prior to the spike (preset length)
function sta(stim, ρ, num_timesteps)

	# sta_val will store the Spike Triggered Average of the desired preset length
    sta_val = zeros(num_timesteps)

    # Get the indices at the start of each of the presets (ρ records spikes as non-zero values)
    # ignore spikes before num_timesteps as there is not enough data to include those in the average
    start_times = find(ρ[num_timesteps:end])

    # The number of presets
    num_presets = length(start_times)

	# Current positions in the presets (vector)
	current = start_times
	for i = 1:length(sta_val)
		# Average of the current index (there are num_presets)
		sta_val[i] = sum(stim[current]) / num_presets
		current += 1
	end

	return sta_val
end

max_stars_repo_name: voutasaurus/comp-neuro

max_issues_repo_name: voutasaurus/comp-neuro

max_forks_repo_name: voutasaurus/comp-neuro


text:
using Test

ti = time()

include("orbitdiagrams/orbitdiagram_tests.jl")
include("orbitdiagrams/poincare_tests.jl")

include("chaosdetection/lyapunov_exponents.jl")
include("chaosdetection/gali_tests.jl")
include("chaosdetection/partially_predictable_tests.jl")
include("chaosdetection/01test.jl")
include("chaosdetection/expansionentropy_tests.jl")

include("period_return/periodicity_tests.jl")
include("period_return/period_tests.jl")
include("period_return/transit_time_tests.jl")

include("dimensions/entropy_dimension.jl")
include("dimensions/correlation_dimension_takens.jl")
include("nlts_tests.jl")

ti = time() - ti
println("\nTest took total time of:")
println(round(ti, digits = 33), " seconds or ",
round(ti/60, digits = 3), " minutes")


text:
module CorrelationDimension

import Plots

using Statistics: mean, covm, varm

using Distributions: fit
using LinearAlgebra: norm
using StatsBase: Histogram

export correlation_dimension, plot_correlation

"""
    correlation_dimension(x; dist=(x1,x2)->norm(x1.-x2)) -> (corrdim, logr, logNr, intercept)

Compute an estimate of the correlation dimension for a set of points in space, with
coordinates for the ith point given as `x[:,i]`.

This is given by

δ̂ = ∂[log N(r)]/∂[log r]

where N(r) is the number of point pairs within a distance r of each other.  This is
estimated by binning the pair-distances into bins (chosen automatically)
and fitting a straight line through the log-distribution.

The distance function `dist` by default is the Cartesian 2-norm, but any function may be
supplied.

The function returns a named tuple of the estimate of the correlation dimension, `corrdim`,
the lower edges of the distance bins, `logr`, the weights in the distance bins `logNr`, and
the `intercept` of the best fitting line (of which the slope is `corrdim`).
"""
function correlation_dimension(x; dist=(a,b)->norm(a .- b))
    D, n = size(x)
    # List of log(pairwise distances)
    point1 = Vector{float(eltype(x))}(undef, D)
    point2 = similar(point1)
    logr = float(eltype(x))[]
    for i in 1:n, j in (i+1):n
        point1 .= x[:,i]
        point2 .= x[:,j]
        push!(logr, log(dist(point1, point2)))
    end
    # Construct histogram
    h = fit(Histogram, logr)
    # Linear regression
    logr = h.edges[1][1:end-1]
    logNr = log.(cumsum(h.weights./(n*(n-1)/2)))
    c, δ̂ = linear_regression(logr, logNr)
    (corrdim=δ̂, logr=logr, logNr=logNr, intercept=c)
end

"""
    correlation_dimension(x1, x2, ..., xN; kwargs...) -> δ̂

Specify coordinates with dimension N as N separate vectors.
"""
correlation_dimension(x, y, args...; kwargs...) = correlation_dimension(permutedims(hcat(x, y, args...)))

"""
    plot_correlation((logr, logNr, slope, intercept))

Create a plot of the output from `correlation_dimension`.  Input is a tuple
of the log of the distance bins (`logr`), the weights of the histogram (`logNr`),
plus the `slope` and `intercept` of the best fitting line.  The `slope` therefore
is the estimate of the ccorrelation dimension.
"""
function plot_correlation((corrdim, logr, logNr, intercept))
    Plots.scatter(logr, logNr, label="Data", legend=:topleft)
    Plots.plot!(x->corrdim*x+intercept, label="y = $(round(corrdim, sigdigits=3))x + $(round(intercept, sigdigits=3))")
end

"""
    linear_regression(x, y)

Perform simple linear regression using Ordinary Least Squares. Returns `a` and `b` such
that `a + b*x` is the closest straight line to the given points `(x, y)`, i.e., such that
the squared error between `y` and `a + b*x` is minimized.
"""
function linear_regression(x::AbstractVector, y::AbstractVector)
    size(x) == size(y) || throw(DimensionMismatch("x and y must be the same size"))
    mx, my = mean(x), mean(y)
    b = covm(x, mx, y, my)/varm(x, mx)
    a = my - b*mx
    return a, b
end

end # module

max_stars_repo_name: anowacki/assorted-julia-modules

max_issues_repo_name: anowacki/assorted-julia-modules

max_forks_repo_name: anowacki/assorted-julia-modules


text:
function matrix_free_apply2f!(
    f::AbstractVector{T},
    elementinfo::ElementFEAInfo{dim,T},
    M,
    vars,
    problem::StiffnessTopOptProblem,
    penalty,
    xmin,
    applyzero::Bool=false,
) where {dim,T}
    @unpack Kes, black, white, varind, metadata = elementinfo
    @unpack dof_cells, cell_dofs = metadata
    @unpack ch = problem
    @unpack values, prescribed_dofs = ch

    update_f!(
        f,
        values,
        prescribed_dofs,
        applyzero,
        dof_cells,
        cell_dofs,
        black,
        white,
        Kes,
        xmin,
        penalty,
        vars,
        varind,
        M,
    )

    return nothing
end

function update_f!(
    f::Vector{T},
    values,
    prescribed_dofs,
    applyzero,
    dof_cells,
    cell_dofs,
    black,
    white,
    Kes,
    xmin,
    penalty,
    vars,
    varind,
    M,
) where {T}
    for ind in 1:length(values)
        d = prescribed_dofs[ind]
        v = values[ind]
        m = size(Kes[ind], 1)

        r = dof_cells.offsets[d]:(dof_cells.offsets[d + 1] - 1)
        if !applyzero && v != 0
            for idx in r
                (i, j) = dof_cells.values[idx]
                if PENALTY_BEFORE_INTERPOLATION
                    px = ifelse(
                        black[i],
                        one(T),
                        ifelse(white[i], xmin; px=density(penalty(vars[varind[i]]), xmin)),
                    )
                else
                    px = ifelse(
                        black[i],
                        one(T),
                        ifelse(white[i], xmin; px=penalty(density(vars[varind[i]], xmin))),
                    )
                end
                if eltype(Kes) <: Symmetric
                    Ke = Kes[i].data
                else
                    Ke = Kes[i]
                end
                for row in 1:m
                    f[cell_dofs[row, i]] -= px * v * Ke[row, j]
                end
            end
        end
        f[d] = M * v
    end
    return nothing
end


text:
export calcCalibResidualAprilTags!, calcCornerProjectionsAprilTags!

function _calcCornerProjectionsAprilTag!( cimg_, tags_, idx, horid, verid, tl, tl_2, cTt_,
                                          f_width, f_height, c_width, c_height,
                                          HORI::Int, VERT::Int;
                                          dodraw::Bool=false    )
  #
  resid = 0.0

  verz = (verid - idx.I[1])*(2*tl)
  horz = (horid - idx.I[2])*(2*tl)

  # reproject the corners
  tC = [-tl_2  tl_2   tl_2 -tl_2;
        tl_2  tl_2  -tl_2 -tl_2;
          0     0     0     0;
          1     1     1     1]
  #

  # offset to the specific tag on board
  tC[1,:] .+= horz
  tC[2,:] .+= verz

  # get corners in camera frame
  cP_ = (cTt_*tC)[1:3,:]
  # get corners in image reference frame (iPc)
  # very basic pinhole camera model
  iPc_ = cP_[1:2,:]
  for c in 1:4
    # inverse camera transform? should more than just f_width
    iPc_[:,c] .*= f_width/cP_[3,c]
  end
  iPc_[1,:] .+= c_width
  iPc_[2,:] .+= c_height

  # get the right tag reference measurement
  j_ = VERT*(horid-1) + verid
  # where are the measured corners
  for c in 1:4
    resid += (tags_[j_].p[c] - iPc_[1:2,c]).^2 |> sum
  end

  # also draw the new positions if desired
  if dodraw
    # drawing requires Int
    iPc_I = round.(Int, iPc_)
    for c in 1:4
      # consolidate all drawing functions in the same prior src file
      _drawCrossOnImg!(cimg_, (iPc_I[1,c],iPc_I[2,c]), 50, RGB{N0f8}(1,0,0))
    end
  end

  return resid
end

"""
    $SIGNATURES

April grid calibration helper function to assemble the cost for a single image.
This function can also be used to draw the corner points on images
to show how the cost function is constructed as well as the calibration performance.
See [`calcCalibResidualAprilTags!`](@ref) for details.
"""
function calcCornerProjectionsAprilTags!(cimg_::AbstractMatrix{<:Colorant},
                                          tags_::AbstractVector;
                                          tagList::AbstractVector=1:length(tags_),
                                          taglength::Real=0.0315,
                                          f_width::Real = size(cimg_,1),
                                          f_height::Real = f_width,
                                          c_width::Real = size(cimg_,2)/2,
                                          c_height::Real = size(cimg_,1)/2,
                                          s::Real=0.0,
                                          VERT::Int=5,
                                          HORI::Int=8,
                                          boardPattern = reshape(1:(VERT*HORI), VERT, HORI),
                                          dodraw = false )
  #
  resid = 0.0

  tl = taglength
  tl_2 = tl/2

  # FIXME use standard CameraModels.jl package
  K_ = [f_width s c_width;
        0 f_height c_height]

  @assert length(tags_) == HORI*VERT "Number of detected tags must equal HORI*VERT=$(HORI*VERT) but only finding $(length(tags_))"
  # loop through all 40 tags
  for j in tagList
    cTt, err1 = AprilTags.tagOrthogonalIteration( tags_[j].p,
                                                  tags_[j].H,
                                                  f_width,
                                                  f_height,
                                                  c_width,
                                                  c_height,
                                                  taglength = taglength)
    #
    cTt_ = [cTt; 0 0 0 1]
    idx = findfirst(x->x==tags_[j].id, boardPattern)

    # go through all the tags on the grid, to compare against this one idx tags[j]
    for verid in 1:VERT, horid in 1:HORI
      resid += _calcCornerProjectionsAprilTag!( cimg_, tags_, idx, horid, verid, tl, tl_2, cTt_,
                                                f_width, f_height, c_width, c_height,
                                                HORI, VERT;
                                                dodraw=dodraw    )
    end
  end

  #
  return resid
end

"""
    $SIGNATURES

Grid of AprilTags calibration helper function.  This function sets up a squared cost to be minimized
by any method, including `Optim.opimize`.  The cost function is constructed by assuming a grid of 40 AprilTags
are of equal size and on a common co-planar surface like a computer screen.

The cost function is constructed by predicting from each tag detection individually where the corners of all
other 39 tags should be on the co-planar surface.  The discrepanc_height between the predicted and detected tag
positions are square accumulated.

### Notes
- Assume regular spacing grid of 40 AprilTags 36h11, ids=1:40,
  - down 1:5 by 8 columns 1:5:40 across where taglength and spacing gap are equal.
  - See grid file at AprilTags/data/CameraCalibration/AprilTagGrid1to40.png
- Current implementation requires all 40 tags to be visible and detected
  - Contributions welcome to allow a few missed detections.
- Take pictures of the grid with the camera you want to calibrate, making sure all tags are clearly visible.
  - Combinations of the grid nearer and further, centered and slanted around the edges of the field of view are best.
  - Any number of images can be used, the more the better.
  - Note that half a doen high-res images might take up to 20 mins or so to optimize.
- Julia Images.jl follows the common ``::Array` column-major---i.e. vertical-major---index convention
  - That is `img[vertical, horizontal]`
  - See https://evizero.github.io/Augmentor.jl/images/#Vertical-Major-vs-Horizontal-Major-1
- This function has the ability to draw the predicted tag corners, see keyword `dodraw=true`.
- This is not a copy of any other AprilCal or such software, this was newly written code out of pure frustration
  when getting stuff to work.  The more native Julia code the better, because Julia is much more mobile and versatile
  than any of the other calibration software out there.  See DevNotes below for roadmap of features to add,
  contributions welcome.

### Example

Also see `AprilTags/examples/AprilTagsGridClibration.jl`.

This example shows how a series of photos of the tag grid image (just use your computer screen, not a projector)
can be used to calibrate a camera.  This example only shows the basic pinhole parameters, although more are possible,
see keyword arguments for which calibration parameters are available.  The latter part shows how to draw crosses
to see before and after result.
```julia
using AprilTags
using FileIO

# where are the photos of the calibration files
filepaths = ["photo1.jpg"; "photo2.jpg";...]
# load the images into memory
imgs = load.(filepaths)

# It's imporant that you measure and specify the tag length correctly here
# 30 mm is just a guess, insert your own correct tag measurements here.
taglength = 0.03

# rough guess of what calibration parameters might be
# img[rows,columns] <==> img[height, width]
c_width = size(imgs[1],2) / 2 # columns across in Images.jl
c_height = size(imgs[1],1) / 2 # rows down in Images.jl
f_width = size(imgs[1],1)
f_height=f_width

#  detect the tags and duplicate the memory before freeing the detector
detector = AprilTagDetector()
tags = detector.(imgs) .|> deepcopy
# remember to free detector later

# setup the cost function, you can add more parameters here if you like
obj = (f_width, f_height, c_width, c_height) -> calcCalibResidualAprilTags!( imgs, tags, taglength=taglength, f_width=f_width, f_height=f_height, c_width=c_width, c_height=c_height, dodraw=false )
obj_ = (fc_wh) -> obj(fc_wh...)

# check that it works
obj_([f_width, f_height, c_width, c_height])

## Bring in the Optim.jl optimization routines
using Optim

# Run the optimization. BFGS is slower by more precise, it's okay to mix and match as coarse and fine optimization stages
result = optimize(obj_, [f_width; f_height; c_width; c_height], BFGS(), Optim.Options(x_tol=1e-8))

# see the optimized calibration parameters
@show f_width_, f_height_, c_width_, c_height_ = (result.minimizer...,)

## show the before and after images to visually confirm things are working
using ImageView

# bad calibration
img1_before = deepcopy(imgs[1])
calcCornerProjectionsAprilTags!(img1_before, taglength=taglength, f_width=f_width, f_height=f_height, c_width=c_width, c_height=c_height, dodraw=true)
imshow(img1_before)

# new calibration
img1_after = deepcopy(imgs[1])
calcCornerProjectionsAprilTags!(img1_after, taglength=taglength, f_width=f_width_, f_height=f_height_, c_width=c_width_, c_height=c_height_, dodraw=true)
imshow(img1_after)

# free the detector memory
freeDetector!(detector) # could also use a deepcopy to duplicate the memory to a secondary location and free the primary immediately
```

### DevNotes
- FIXME common JuliaRobotics/CameraModels.jl package shoudl be made
- TODO Radial distortion parameters should be added, see https://en.wikipedia.org/wiki/Distortion_(optics)
- TODO allow missed detections on grid of 40 tags
- TODO auto-detect the grid so that any grid can be used (still assuming regular spacing)

### Related

[`calcCornerProjectionsAprilTags!`](@ref)
"""
function calcCalibResidualAprilTags!( images::AbstractVector,
                                      allTags::AbstractVector;
                                      tagList::AbstractVector=1:length(allTags),
                                      taglength = 0.0315,
                                      VERT = 5,
                                      HORI = 8,
                                      f_width::Real=size(images[1],2),
                                      f_height::Real=f_width,
                                      c_width::Real = size(images[1],2) / 2,
                                      c_height::Real = size(images[1],1) / 2,
                                      s::Real=0.0,
                                      boardPattern = reshape(1:(VERT*HORI), VERT, HORI),
                                      dodraw = false  )
  #

  #
  resid = 0.0

  for count in 1:length(images)

    # detect tags in current image
    cimg_ = images[count]
    tags_ = allTags[count]

    # do the actual calculations
    resid += calcCornerProjectionsAprilTags!( cimg_, tags_,
                                              tagList=tagList,
                                              taglength=taglength,
                                              f_width=f_width,
                                              f_height=f_height,
                                              c_width=c_width,
                                              c_height=c_height,
                                              s=s,
                                              VERT=VERT, HORI=HORI,
                                              boardPattern=boardPattern,
                                              dodraw=false )
    #
  end

  return resid
end


text:
using FluxTikz
using Test

@testset "FluxTikz.jl" begin
    # Write your tests here.
end


text:
"""
    Padding <: Construct{Nothing}

Represents padding data.
"""
struct Padding <: Construct{Nothing}
    size::UInt

    Padding(size::Integer = 0) = new(convert(UInt, size))
end

function deserialize(cons::Padding, s::IO; contextkw...)
    skip(s, cons.size)
    nothing
end
serialize(cons::Padding, s::IO, ::Nothing; contextkw...) = write(s, zeros(UInt8, cons.size))
estimatesize(cons::Padding; contextkw...) = cons.size


text:
module STAMP

import Automa
import Automa.RegExp: @re_str
import Automa.Stream: @mark, @markpos, @relpos, @abspos
import BioGenerics: BioGenerics, isfilled
# import BioGenerics.Exceptions: missingerror
import BioGenerics.Automa: State
using BioSequences
import TranscodingStreams: TranscodingStreams, TranscodingStream

include("record.jl")
include("readrecord.jl")
include("reader.jl")

end # module


text:
using Symbolics, LinearAlgebra, DifferentialEquations

# ================== CUSTOM TYPES ================== #

# FINITE DIFFERENCE
abstract type Difference end
struct ComplexDifference <: Difference end
struct CentralDifference <: Difference
	points::Int64
	derivative::Int64
	function CentralDifference(points=3, derivative=1)
		return new(points, derivative)
	end
end
struct ForwardDifference <: Difference
	points::Int64
	derivative::Int64
	function ForwardDifference(points=3, derivative=1)
		return new(points, derivative)
	end
end
struct BackwardDifference <: Difference
	points::Int64
	derivative::Int64
	function BackwardDifference(points=3, derivative=1)
		return new(points, derivative)
	end
end

# ================== CONSTANTS ================== #
# CENTRAL DIFFERENCE QUERIES
cntr3pt = Dict( 1 => (f, x, δx) -> (-f(x - δx) + f(x + δx))/(2*δx),        # Accuracy: δx²
				2 => (f, x, δx) -> (f(x - δx) -2*f(x) + f(x + δx))/(δx^2)) # Accuracy: δx²

cntr5pt = Dict( 1 => (f, x, δx) -> (f(x - 2*δx) - 8*f(x - δx) + 8*f(x + δx) - f(x + 2*δx))/(12*δx),                # Accuracy: δx⁴
				2 => (f, x, δx) -> (-f(x - 2*δx) + 16*f(x - δx) - 30*f(x) + 16*f(x + δx) - f(x + 2*δx))/(12*δx^2), # Accuracy: δx⁴
				3 => (f, x, δx) -> (-f(x - 2*δx) + 2*f(x - δx) -2*f(x + δx) + f(x + 2*δx))/(2*δx^3))               # Accuracy: δx²

cntr7pt = Dict( 1 => (f, x, δx) -> (-f(x - 3*δx) + 9*f(x - 2*δx) - 45*f(x - δx) + 45*f(x + δx) - 9*f(x + 2*δx) + f(x + 3*δx))/(60*δx),                     # Accuracy: δx⁶
				2 => (f, x, δx) -> (2*f(x - 3*δx) - 27*f(x - 2*δx) + 270*f(x - δx) -490*f(x) + 270*f(x + δx) - 27*f(x + 2*δx) + 2*f(x + 3*δx))/(180*δx^2), # Accuracy: δx⁶
				3 => (f, x, δx) -> (f(x - 3*δx) - 8*f(x - 2*δx) + 13*f(x - δx) - 13*f(x + δx) + 8*f(x + 2*δx) - f(x + 3*δx))/(8*δx^3))                     # Accuracy: δx⁴

cntrDiff = Dict(3 => cntr3pt, 5 => cntr5pt, 7 => cntr7pt)

# FORWARD DIFFERENCE

# BACKWARD DIFFERENCE

# ================== FUNCTIONS ================== #

# FINITE DIFFERENCES
# function gradient(f::Function, x::Vector, h::Union{Float64, Int64}, method::Difference = ComplexDifference())
function gradient(f, x::Vector, h::Union{Float64, Int64}, method = ComplexDifference())
	return _gradient(method, f, x, h)
end

function _gradient(method::ComplexDifference, f, x, h)
	∂f = []
	L = length(x)
	for i = 1:L
		δx = zeros(Complex{Float64}, L)
		δx[i] = h*im
		push!(∂f, imag(f(x + δx))/h)
	end
	return ∂f
end

# TODO: Make work with 2nd+ order derivatives (ie: 1x1 )
function _gradient(method::CentralDifference, f, x, h)
	∂ⁿf = []
	L = length(x)
	g = cntrDiff[method.points][method.derivative]
	for i = 1:L
		δx = zeros(L)
		δx[i] = h
		push!(∂ⁿf, g(f, x, δx))
	end
	return ∂ⁿf
end

# TODO: Make Forward and Backwards Differences

# ------------- OLD CODE ------------- #

@variables r μ x y z

# SETUP
g = -μ/norm([x, y, z])^3 * [x, y, z]
jac = Array{Any}(undef, 3, 3)
hes = Array{Any}(undef, 3, 3, 3)
sub = Dict(sqrt(abs2(x) + abs2(y) + abs2(z)) => r);
D = [Differential(x), Differential(y), Differential(z)]

# CALCULATING PARTIALS
for i = 1:3, j = 1:3, k = 1:3

	# n×n×n HESSIAN
	hes[k, j, i] = substitute(
		expand_derivatives(D[i](-expand_derivatives(D[j](g[k])))),
		sub)

	# n×n JACOBIAN
	jac[k, j] = substitute(
		expand_derivatives(D[j](g[k])),
		sub)
end

# CREATING FUNCTIONS
_buildJacobian = build_function(jac, [x, y, z, μ, r], expression=false)[1]
jacobian(r⃗, μ₀) = _buildJacobian([r⃗[1], r⃗[2], r⃗[3], μ₀, norm(r⃗)])
_buildHessian  = build_function(hes, [x, y, z, μ, r], expression=false)[1]
hessian(r⃗, μ₀) = _buildHessian([r⃗[1], r⃗[2], r⃗[3], μ₀, norm(r⃗)])

export jacobian, hessian, cntrDiff, gradient


text:
function relaxed_initial_state(mdp::NoCrashProblem, steps=200,
                             rng=MersenneTwister(rand(UInt32));
                             solver=BehaviorSolver(NORMAL, true, rng))

    pp = mdp.dmodel.phys_param
    is = MLState(0.0, 0.0, CarState[CarState(pp.lane_length/2, 2.0, pp.v_med, 0.0, NORMAL, 1)], true)
    sim = HistoryRecorder(max_steps=steps, rng=rng)
    policy = solve(solver, mdp)
    #policy = RandomPolicy(mdp)
    simulate_cz(sim, mdp, policy, is)
    #println("simulate done")
    return sim.state_hist[end]
end

##########################  change it

function relaxed_initial_state_cz(mdp::NoCrashProblem, is::MLState, policy::Policy, steps=20
                             #rng=MersenneTwister(rand(UInt32));
                             )
                             #solver=BehaviorSolver(NORMAL, true, rng))
                             #solver=SimpleSolver())

    #pp = mdp.dmodel.phys_param
    #is = MLState(0.0, 0.0, CarState[CarState(pp.lane_length/2, 1.0, pp.v_med, 0.0, NORMAL, 1)])
    rng = MersenneTwister(5)
    sim = HistoryRecorder(max_steps=steps, rng=rng)
    #policy = RandomPolicy(mdp)
    #policy = solve(solver, mdp)
    #policy = Simple(mdp)
    #print("sim")
    simulate_cz(sim, mdp, policy, is)
    #return sim.state_hist[end]

    #print("simulate done")
    return sim
end

function relaxed_initial_state_cz_control(mdp::NoCrashProblem, is::MLState, policy::Policy, steps=20
                             #rng=MersenneTwister(rand(UInt32));
                             )
                             #solver=BehaviorSolver(NORMAL, true, rng))
                             #solver=SimpleSolver())

    #pp = mdp.dmodel.phys_param
    #is = MLState(0.0, 0.0, CarState[CarState(pp.lane_length/2, 1.0, pp.v_med, 0.0, NORMAL, 1)])
    rng = MersenneTwister(5)
    sim = HistoryRecorder(max_steps=steps, rng=rng)
    #policy = RandomPolicy(mdp)
    #policy = solve(solver, mdp)
    #policy = Simple(mdp)
    #print("sim")
    simulate_control(sim, mdp, policy, is)
    #return sim.state_hist[end]

    #print("simulate done")
    return sim
end

function relaxed_initial_state_cz_emergency(mdp::NoCrashProblem, is::MLState, policy::Policy, steps=20
                             #rng=MersenneTwister(rand(UInt32));
                             )
                             #solver=BehaviorSolver(NORMAL, true, rng))
                             #solver=SimpleSolver())

    #pp = mdp.dmodel.phys_param
    #is = MLState(0.0, 0.0, CarState[CarState(pp.lane_length/2, 1.0, pp.v_med, 0.0, NORMAL, 1)])
    rng = MersenneTwister(5)
    sim = HistoryRecorder(max_steps=steps, rng=rng)
    #policy = RandomPolicy(mdp)
    #policy = solve(solver, mdp)
    #policy = Simple(mdp)
    #print("sim")
    #sim, collision_num, step_before_collision = simulate_emergency(sim, mdp, policy, is)

    sim, collision_num, step_before_collision = simulate_emergency(sim, mdp, policy, is)
    #return sim.state_hist[end]

    #print("simulate done")
    return sim, collision_num, step_before_collision
end

function get_traffic_density(sim::HistoryRecorder)
    density = []
    headway = []
    for t = 1:200
        s = sim.state_hist[t]
        num = length(s.cars)
        headway1 = 0
        for i in 2:length(s.cars)
            car = s.cars[i]
            headway1 = headway1 + car.behavior.p_idm.T
        end
        d = num/100
        hw = headway1/(num-1)
        push!(density,d)
        push!(headway,hw)
    end
    return density,headway
end
###############################
#=
struct BehaviorSolver <: Solver
    b::BehaviorModel
    keep_lane::Bool
    rng::AbstractRNG
end

struct BehaviorPolicy <: Policy
    problem::NoCrashProblem
    b::BehaviorModel
    keep_lane::Bool
    rng::AbstractRNG
end
solve(s::BehaviorSolver, p::NoCrashProblem) = BehaviorPolicy(p, s.b, s.keep_lane, s.rng)
=#
############################

function simulate_cz{S,A}(sim::HistoryRecorder,
                       mdp::MDP{S,A}, policy::Policy
                      ,init_state::S=get_initial_state(sim, mdp)
                      )

                      #    init_state = get_initial_state(sim, mdp)

    max_steps = get(sim.max_steps, typemax(Int))
    if !isnull(sim.eps)
        max_steps = min(max_steps, ceil(Int,log(get(sim.eps))/log(discount(mdp))))
    end
    sizehint = get(sim.sizehint, min(max_steps, 1000))

    # aliases for the histories to make the code more concise
    sh = sim.state_hist = sizehint!(Vector{S}(0), sizehint)
    ah = sim.action_hist = sizehint!(Vector{A}(0), sizehint)
    oh = sim.observation_hist = Any[]
    bh = sim.belief_hist = Any[]
    rh = sim.reward_hist = sizehint!(Vector{Float64}(0), sizehint)

    if sim.show_progress
       prog = Progress(max_steps, "Simulating..." )
    end

    push!(sh, init_state)

    disc = 1.0
    step = 1
    try
        while !isterminal(mdp, sh[step]) && step <= max_steps
            action_calculate = action(policy, sh[step])
            action_use = control_action(mdp,sh[step])
            push!(ah, action(policy, sh[step]))
            #push!(ah, action_use)
            sp, r = generate_sr(mdp, sh[step], ah[step], sim.rng)

        #    if action_calculate == action_use
        #       sp.Control_Signal = sh[step].Control_Signal
        #    else
        #       sp.Control_Signal = false
        #    end

            push!(sh, sp)
            push!(rh, r)

            disc *= discount(mdp)
            step += 1

            if sim.show_progress
                    next!(prog)
            end
        end
    catch ex
        if sim.capture_exception
            sim.exception = ex
            sim.backtrace = catch_backtrace()
        else
            rethrow(ex)
        end
    end

    if sim.show_progress
        finish!(prog)
    end

    return MDPHistory(sh, ah, rh, discount(mdp), sim.exception, sim.backtrace)

    return 0

end

function simulate_control{S,A}(sim::HistoryRecorder,
                       mdp::MDP{S,A}, policy::Policy
                      ,init_state::S=get_initial_state(sim, mdp)
                      )

                      #    init_state = get_initial_state(sim, mdp)

    max_steps = get(sim.max_steps, typemax(Int))
    if !isnull(sim.eps)
        max_steps = min(max_steps, ceil(Int,log(get(sim.eps))/log(discount(mdp))))
    end
    sizehint = get(sim.sizehint, min(max_steps, 1000))

    # aliases for the histories to make the code more concise
    sh = sim.state_hist = sizehint!(Vector{S}(0), sizehint)
    ah = sim.action_hist = sizehint!(Vector{A}(0), sizehint)
    oh = sim.observation_hist = Any[]
    bh = sim.belief_hist = Any[]
    rh = sim.reward_hist = sizehint!(Vector{Float64}(0), sizehint)

    if sim.show_progress
       prog = Progress(max_steps, "Simulating..." )
    end

    push!(sh, init_state)

    disc = 1.0
    step = 1
    try
        while !isterminal(mdp, sh[step]) && step <= max_steps
            #action_calculate = action(policy, sh[step])
            action_use = control_action_adapt(mdp,sh[step])
            push!(ah, action_use)
            #push!(ah, action(policy, sh[step]))

            sp, r = generate_sr(mdp, sh[step], ah[step], sim.rng)
            push!(sh, sp)
            push!(rh, r)

            disc *= discount(mdp)
            step += 1

            if sim.show_progress
                    next!(prog)
            end
        end
    catch ex
        if sim.capture_exception
            sim.exception = ex
            sim.backtrace = catch_backtrace()
        else
            rethrow(ex)
        end
    end

    if sim.show_progress
        finish!(prog)
    end

    return MDPHistory(sh, ah, rh, discount(mdp), sim.exception, sim.backtrace)

    return 0

end

function judge_collision(mdp::NoCrashProblem,s::Union{MLState,MLObs})
    nb_cars = length(s.cars)
    pp = mdp.dmodel.phys_param
    car1 = s.cars[1]
    for i in 2:nb_cars
        car_k = s.cars[i]
        if abs(car_k.x-car1.x)<pp.l_car+2
            if occupation_overlap(car1.y, car_k.y)
                return true
            end
        end
    end
    return false
end

function simulate_emergency{S,A}(sim::HistoryRecorder,
                       mdp::MDP{S,A}, policy::Policy
                      ,init_state::S=get_initial_state(sim, mdp)
                      )

                      #    init_state = get_initial_state(sim, mdp)

    collision_num = 0
    max_steps = get(sim.max_steps, typemax(Int))
    if !isnull(sim.eps)
        max_steps = min(max_steps, ceil(Int,log(get(sim.eps))/log(discount(mdp))))
    end
    sizehint = get(sim.sizehint, min(max_steps, 1000))

    # aliases for the histories to make the code more concise
    sh = sim.state_hist = sizehint!(Vector{S}(0), sizehint)
    ah = sim.action_hist = sizehint!(Vector{A}(0), sizehint)
    oh = sim.observation_hist = Any[]
    bh = sim.belief_hist = Any[]
    rh = sim.reward_hist = sizehint!(Vector{Float64}(0), sizehint)

    if sim.show_progress
       prog = Progress(max_steps, "Simulating..." )
    end

    push!(sh, init_state)

    disc = 1.0
    step = 1
    try
        while !isterminal(mdp, sh[step]) && step <= max_steps
            #action_calculate = action(policy, sh[step])
            action_use = control_action_emergency(mdp,sh[step],true)
            push!(ah, action_use)
            #push!(ah, action(policy, sh[step]))

            sp, r = generate_sr(mdp, sh[step], ah[step], sim.rng)
            push!(sh, sp)
            push!(rh, r)

            disc *= discount(mdp)
            step += 1
            #######################Get Collision
            nb_cars = length(sp.cars)
            pp = mdp.dmodel.phys_param
            car1 = sp.cars[1]
            for i in 2:nb_cars
                car_k = sp.cars[i]
                if abs(car_k.x-car1.x)<pp.l_car#+5
                    if occupation_overlap(car1.y, car_k.y)
                        #collision_num = collision_num+1
                    end
                end
            end
            if judge_collision(mdp,sp)
                collision_num = collision_num+1
            end

            if collision_num > 0
                break
            end

            #######################

            if sim.show_progress
                    next!(prog)
            end
        end
    catch ex
        if sim.capture_exception
            sim.exception = ex
            sim.backtrace = catch_backtrace()
        else
            rethrow(ex)
        end
    end

    if sim.show_progress
        finish!(prog)
    end

    return MDPHistory(sh, ah, rh, discount(mdp), sim.exception, sim.backtrace), collision_num, step

    return 0

end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
"""
Constructors for covariance (kernel) functions.

      Constructor                   Description                                Isotropic/Anisotropic
    - SquaredExponentialI(γ,σ):     squared exponential covariance function    isotropic
    - ExponentialI(γ,σ):            exponential covariance function            isotropic
    - RationalQuadraticI(γ,σ,α):    rational quadratic covariance function     isotropic
    - Matern12I(γ,σ):               Matérn covariance function with ʋ = 1/2.   isotropic
    - Matern32I(γ,σ):               Matérn covariance function with ʋ = 3/2.   isotropic
    - Matern52I(γ,σ):               Matérn covariance function with ʋ = 5/2.   isotropic
    - SMP(w,μ,γ):                   Spectral mixture product cov. function

Distance metrics

    - euclidean_distance            l²-norm:  d(x,x') = || x - x' ||",
    - derivative_distance           H¹-norm:  d(x,x') = || diff(x)./diff(z) - diff(x')./diff(z) ||",
    - antiderivative_distance       H⁻¹-norm: d(x,x') = || diff(x).*diff(z) - diff(x').*diff(z) ||"
"""
abstract type Kernel end

#  *--*--*--*--*--*--*--*--*--*--*--*
#  | Isotropic covariance functions |
#  *--*--*--*--*--*--*--*--*--*--*--*

""" SquaredExponentialI(γ,σ): squared exponential covariance function, isotropic """
struct SquaredExponentialI <: Kernel
    # Hyperparameters
    # "Length scale"
    γ::Float64
    # "Signal variance"
    σ::Float64
    # "Distance metric"
    d::Function
end

# evaluates the kernel function for a given pair of inputs
function kernel_function(k::SquaredExponentialI; z=nothing)
    # k(x,x') = σ * exp( - d(x,x')² / 2γ² )
    σ = k.σ
    γ = k.γ
    d = k.d
  evaluate(a,b) = σ * exp(- d(a,b,z)^2 / 2*γ^2 )
  return evaluate
end

struct Matern12I <: Kernel
    # Hyperparameters
    γ::Float64 #length scale
    σ::Float64 #signal variance
    d::Function #distance metric
end

function kernel_function(k::Matern12I; z=nothing)
  # k(x,x') = σ * exp( - ||x-x'|| / γ )
  σ = k.σ
  γ = k.γ
  d = k.d
  evaluate(a,b) = σ * exp(- d(a,b,z) / γ )
  return evaluate
end

struct Matern32I <: Kernel
    # Hyperparameters
    γ::Float64 #length scale
    σ::Float64 #signal variance
    d::Function #distance metric
end

function kernel_function(k::Matern32I; z=nothing)
    # k(x,x') = σ * (1+c) * exp(-√(3)*||x-x'||)/γ)
    σ = k.σ
    γ = k.γ
    d = k.d
    function evaluate(a,b)
        c = sqrt(3)*d(a,b,z)/γ
        return σ * (1+c) * exp(-c)
    end
  return evaluate
end

struct Matern52I <: Kernel
    # Hyperparameters
    γ::Float64 #length scale
    σ::Float64 #signal variance
    d::Function #distance metric
end

function kernel_function(k::Matern52I; z=nothing)
    # k(x,x') = σ * ( 1 + √(5)*||x-x'||)/γ + 5*||x-x'||²/(3*γ^2) ) * exp(-√(5)*||x-x'||)/γ)
    σ = k.σ
    γ = k.γ
    d = k.d
    function evaluate(a,b)
        g = sqrt(5)*d(a,b,z)/γ
        h = 5*(d(a,b,z)^2)/(3*γ^2)
        return σ * (1+g+h) * exp(-g)
    end
  return evaluate
end

struct RationalQuadraticI <: Kernel
    # Hyperparameters
    γ::Float64 #length scale
    σ::Float64 #signal variance
    α::Float64 #shape parameter
    d::Function #distance metric
end

function kernel_function(k::RationalQuadraticI; z=nothing)
    # k(x,x') = σ * (1+(x-x')'*(x-x')/(2*α*(γ²))^(-α)
    σ = k.σ
    α = k.α
    d = k.d
    l = γ^2 # squared length scale
    function evaluate(a,b)
       return σ * (1+(a-b)'*(a-b)/(2*α*l))^(-α)
  end
  return evaluate
end

struct SpectralMixtureProductI <: Kernel
    w::Array{Float64} #mixture weights
    μ::Array{Float64} #spectral means
    γ::Array{Float64} #spectral variances
end

function SpectralMixtureProductI(hyp)
    if length(hyp) % 3 != 0
        throw(error("Inconsistent number of components. Length of parameter vector should be a multiple of three."))
    end

    Q = Int(length(hyp)/3)

    SpectralMixtureProductI(hyp[1:Q],hyp[Q+1:2Q],hyp[2Q+1:3Q])
end

function kernel_function(k::SpectralMixtureProductI; z=nothing)
    w = (k.w .^2)' # square mixture weights
    μ = k.μ
    γ = k.γ
    h(arg1, arg2) = exp.(-0.5 * arg1) .* cos.(arg2)

    function evaluate(a,b)
        τ = (a .- b) * 2*pi
        D  = length(a)

        K = 1
        for d=1:D
            # println(w)
            # println(h((τ[d] .^ 2)*k.γ, τ[d]*k.μ))
            K *= w * h((τ[d] .^ 2)*γ, τ[d]*μ)
        end
        K
    end
   return evaluate
end

struct SpectralMixtureProductA <: Kernel
    w::Array{Float64} # Mixture weights, D x Q array
    μ::Array{Float64} # Spectral means, D x Q array
    γ::Array{Float64} # Spectral variances, D x Q array
end

function SpectralMixtureProductA(hyp, D)
    Q = Int(floor(length(hyp)/(3D)))
    w = reshape(  hyp[                      1:D*Q],D,Q);   # mixture weights
    μ = reshape(  hyp[D .* Q .+           (1:D*Q)],D,Q);   # spectral means
    γ = reshape(  hyp[D .* Q .+ D .* Q .+ (1:D*Q)],D,Q);   # spectral variances
    SpectralMixtureProductA(w, μ, γ)
end

function kernel_function(k::SpectralMixtureProductA; z=nothing)
    w = (k.w .^2)' # square mixture weights
    μ = k.μ
    γ = k.γ
    h(arg1, arg2) = exp.(-0.5 * arg1).*cos.(arg2);
    D,Q  = size(w)

    function evaluate(a,b)
        τ = (a .- b) * 2*pi

        K = 1
        for d=1:D
            K = K .* (w[d,:]' * h((τ[d] .^ 2)*γ[d,:], τ[d]*μ[d,:]) )
        end
        K
    end
  return evaluate
end

##

# """
# https://github.com/alshedivat/gpml/blob/master/cov/covSM.m
#
# hyp = [ log(w(:))
#          log(m(:))
#          log(sqrt(v(:))) ]
#
# Q: number of components
# D: dimension of the input
#
# """

#

##
# function SMP(Q,hyp,a,b)
#       D = length(a);
#       n = 1                                                   # dimensionality
#       w = reshape(  hyp[                     1:D*Q],D,Q);     # mixture weights
#       m = reshape(  hyp[D .* Q .+           (1:D*Q)],D,Q);    # spectral means
#       v = reshape(2*hyp[D .* Q .+ D .* Q .+ (1:D*Q)],D,Q);    # spectral variances
#
#       T = 2*pi*bsxfun(@minus,reshape(a,n,1,D),reshape(b,1,[],D));
#
#       T = reshape(T,[],D);
#
#       K = 1;
#       w = reshape(w,Q,D)';
#       m = reshape(m,Q,D)';
#       v = reshape(v,Q,D)';
#
#       h(t2v, tm) = exp(-0.5 * t2v).*cos(tm);
#       for d=1:D
#           K = K .* ( h( (T[:,d] .* T[:,d]) * v[d,:], T(:,d)*m(d,:) )*w[d,:]' );
#       end
#       K = reshape(K.*ones(size(T,1),1),n,[]);
#
#       return K
# end

# function covSMfast(Q,hyp,x,z)
#     smp = Q<0
#     Q = abs(Q)
#     n,D = size(x); P = smp*D+(1-smp);               # dimensionality, P=D or P=1
#     w = exp(reshape(  hyp(         1:P*Q) ,P,Q));     # mixture weights
#     m = exp(reshape(  hyp(P*Q+    (1:D*Q)),D,Q));     # spectral means
#     v = exp(reshape(2*hyp(P*Q+D*Q+(1:D*Q)),D,Q));     # spectral variances
#
#     T = reshape(T,[],D);
#
#     if smp
#         h(t2v, tm) = exp(-0.5*t2v).*cos(tm);
#         K = 1; w = reshape(w,Q,P)'; m = reshape(m,Q,D)'; v = reshape(v,Q,D)';
#         for d=1:D
#             K = K .* ( h( (T(:,d).*T(:,d))*v(d,:), T(:,d)*m(d,:) )*w(d,:)' );
#     end
#     K = reshape(K.*ones(size(T,1),1),n,[]);
#     else
         # E = exp(-0.5*(T.*T)*v); H = E.*cos(T*m);
         # K = reshape(H*w',n,[]);
         # return K
#   end
# end
#
#
#
# function SM(Q,hyp,x,z)
#   Q = abs(Q)
#   n,D = size(x);                                  # dimensionality
#   w = exp(reshape(  hyp(         1:Q) ,1,Q));     # mixture weights
#   m = exp(reshape(  hyp(Q+    (1:D*Q)),D,Q));     # spectral means
#   v = exp(reshape(2*hyp(Q+D*Q+(1:D*Q)),D,Q));     # spectral variances
#
#   # w = exp(reshape(  hyp(         1) ,1,Q));     # mixture weights
#   # m = exp(reshape(  hyp(1+    (1:D)),D,Q));     # spectral means
#   # v = exp(reshape(2*hyp(1+16+ (1:D)),D,Q));     # spectral variances
#
#   T = reshape(T,[],D);
#   E = exp(-0.5*(T.*T)*v);
#   H = E.*cos(T*m);
#   K = reshape(H*w',n,[]);
#   return K
# end
#
#

# function SMP(hyp,a,b)
#
#       τ = (a .- b) * 2*pi
#
#       D = length(a);
#       Q = Int(length(hyp)/(3*D))
#
#       println(Q)
#
#       w .^= 2 # square mixture weights
#
#       h(arg1, arg2) = exp.(-0.5 * arg1).*cos.(arg2);
#
#       K = 1;
#       for d=1:D
#           K = K .* ( h((τ[d] .^ 2)*γ[d,:], τ[d]*μ[d,:] ) * w[d,:]' );
#       end
#       K
# end

# hyp = [3,3,3,4,4,4,2.095, 2.095, 2.095]
#
# γ = 3*0.0001
#
# fn = kernel_function(SMP(hyp, 3); z=nothing)
# fn2 = kernel_function(SquaredExponentialI2(0.1, 3.0, euclidean_distance))
#
# fn([1,2,3], [5,7,8])
# fn2([1,2,3], [5,7,8])
#
# a = [0.03, 0.04]
# b = [0.0, 0.0]
# hyp = [1,2,0.1,0.4,2,1]
# fn = kernel_function(SMP(hyp, 2); z=nothing)
# fn(a,b)
#
#
# SMP(zeros(18), [1,2,3],[4,5,6])
# SMP(5*ones(18), [1,2,3],[4,5,6])
# SMP(hyp, [1,2,3],[4,5,7])
#
#
# a= [1,2,3]
# b= [7,6,5]
#
# τ = (a .- b) * 2*pi
#
# D = length(a);
# Q = Int(length(hyp)/(3*D))
#
# println(Q)
# w = reshape(  hyp[                      1:D*Q],D,Q);   # mixture weights
# μ = reshape(  hyp[D .* Q .+           (1:D*Q)],D,Q);   # spectral means
# γ = reshape(2*hyp[D .* Q .+ D .* Q .+ (1:D*Q)],D,Q);   # spectral variances
#
# w .^= 2 # square mixture weights
#
# h(arg1, arg2) = exp.(-0.5 * arg1).*cos.(arg2);
#
# K = 1
# for d=1:D
#     println(h((τ[d] .^ 2)*γ[d,:], τ[d]*μ[d,:] ))
#     println(w[d,:]')
#     K .*= ( w[d,:]' * h((τ[d] .^ 2)*γ[d,:], τ[d]*μ[d,:] ));
# end
# K

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
using BibParser, BibParser.BibTeX
using Test

const PACKAGE_ROOT = pkgdir(BibParser)

@testset "BibTeX" begin
    for file in ["test.bib", "error.bib"]
        println("\nstart $file")
        parsed = parse_file(joinpath(PACKAGE_ROOT, "examples", "$file"))

        for (k, e) in pairs(parsed)
            println("|--->$(e.id)")
            for fn in fieldnames(typeof(e))
                println("|     |----$(string(fn)) -> parsed!") # TODO: better printing
            end
            for (n, v) in pairs(e.fields)
                println("|     |----$n -> $v")
            end
            println("|")
        end
    end
end

@testset "CFF" begin
    files = [
        ("CITATION.cff",        true),
        ("invalid_yaml.cff",    false),
        ("invalid_version.cff", false),
        ("invalid_schema.cff",  false)
    ]
    @testset "$file" for (file, expected_result) in files
        parsed, result = parse_file(joinpath(PACKAGE_ROOT, "examples", file), :CFF)

        @test result == expected_result
    end
end


text:
using Oceananigans.Grids: AbstractGrid
using Oceananigans.Operators: ∂xᶠᶜᵃ, ∂yᶜᶠᵃ
using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions

using Adapt

struct ExplicitFreeSurface{E, T}
    η :: E
    gravitational_acceleration :: T
end

ExplicitFreeSurface(; gravitational_acceleration=g_Earth) =
    ExplicitFreeSurface(nothing, gravitational_acceleration)

Adapt.adapt_structure(to, free_surface::ExplicitFreeSurface) =
    ExplicitFreeSurface(Adapt.adapt(to, free_surface.η), free_surface.gravitational_acceleration)

function FreeSurface(free_surface::ExplicitFreeSurface{Nothing}, velocities, arch, grid)
    η = CenterField(arch, grid, TracerBoundaryConditions(grid))
    g = convert(eltype(grid), free_surface.gravitational_acceleration)
    return ExplicitFreeSurface(η, g)
end

explicit_barotropic_pressure_x_gradient(i, j, k, grid, free_surface::ExplicitFreeSurface) =
    free_surface.gravitational_acceleration * ∂xᶠᶜᵃ(i, j, k, grid, free_surface.η)

explicit_barotropic_pressure_y_gradient(i, j, k, grid, free_surface::ExplicitFreeSurface) =
    free_surface.gravitational_acceleration * ∂yᶜᶠᵃ(i, j, k, grid, free_surface.η)


text:
## Julia wrappers of miscellaneous functions described in section 10

function wssortindicesi(m::Integer,n::Integer,ia::Vector{Cint},ja::Vector{Cint})
    info = Array(Cint,1)
    ccall((:ws_sortindices_i_,libwsmp),Void,
          (Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint}),
          &m,&n,ia,ja,info)
    info[1] == 0 || error("Error return $(info[1]) from ws_sortindices_i_")
    nothing
end

function wssortindicesd(m::Integer,n::Integer,ia::Vector{Cint},ja::Vector{Cint},
                          avals::Vector{Cdouble})
    info = Array(Cint,1)
    ccall((:ws_sortindices_d_,libwsmp),Void,
          (Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cint}),
          &m,&n,ia,ja,avals,info)
    info[1] == 0 || error("Error return $(info[1]) from ws_sortindices_i_")
    nothing
end

wssetmaxthreads(n::Integer) = ccall((:wsetmaxthrds_,libwsmp),Void,(Ptr{Cint},),&n)

wssetmaxstack(fstk::Cdouble) = ccall((:wsetmaxstack_,libwsmp),Void,(Ptr{Cdouble},),&fstk)

wssetlf(dlf::Cdouble) = ccall((:wsetlf_,libwsmp),Void,(Ptr{Cdouble},),&dlf)

function wsversion()
    V = Array(Cint,1)
    R = Array(Cint,1)
    M = Array(Cint,1)
    ccall((:wsmp_version_,libwsmp),Void,(Ptr{Cint},Ptr{Cint},Ptr{Cint}),V,R,M)
    (V[1],R[1],M[1])
end

wsclear() = ccall((:wsmp_clear_,libwsmp),Void,())

wsffree() = ccall((:wsffree_,libwsmp),Void,())

wsafree() = ccall((:wsafree_,libwsmp),Void,())

wssfree() = ccall((:wssfree_,libwsmp),Void,())

function wssmatvec(ia::Vector{Cint},ja::Vector{Cint},
                   avals::Vector{Cdouble},x::Vector{Cdouble},b::Vector{Cdouble})
    n = length(ia) - 1
    length(x) == length(b) == n || throw(DimensionMismatch(""))
    length(ja) == length(avals) == (ia[end] - 1) || throw(DimensionMismatch(""))
    ierr = Array(Cint,1)
    ccall((:wssmatvec_,libwsmp),Void,
          (Ptr{Cint},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cint}),
          &n,ia,ja,avals,x,b,ierr)
    ierr[1] == 0 || error("Error return $(ierr[1]) from wssmatvec_")
    b
end

function wssmatvec(ia::Vector{Cint},ja::Vector{Cint},avals::Vector{Cdouble},x::Vector{Cdouble})
    wssmatvec(ia,ja,avals,x,similar(x))
end

function wsstoremat(id::Integer)
    0 ≤ id < 64 || error("id = $id should be in 0:63")
    info = Array(Cint,1)
    ccall((:wstoremat_,libwsmp),Void,(Ptr{Cint},Ptr{Cint}),&id,info)
    info[1] == 0 && return id
    info[1] == -210 && error("Slot $id is occupied")
    error("Error code $(info[1]) from wstoremat_")
end

function wsrecallmat(id::Integer)
    0 ≤ id < 64 || error("id = $id should be in 0:63")
    info = Array(Cint,1)
    ccall((:wrecallmat_,libwsmp),Void,(Ptr{Cint},Ptr{Cint}),&id,info)
    info[1] == 0 && return id
    info[1] == -210 && error("Slot $id is empty")
    error("Error code $(info[1]) from wstoremat_")
end


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# ╔═╡ 13280ba4-064a-11eb-04dd-5b2c0598d84a
using Pkg, DrWatson

# ╔═╡ 134f547a-064a-11eb-0750-b590037632f5
begin
	@quickactivate "SR2TuringPluto"
	using LinearAlgebra
	using Turing
	using StatisticalRethinking
	Turing.setprogress!(false)
end

# ╔═╡ b32d14e8-0648-11eb-15d4-7d13643a1531
md"## Clip-04-26-36t.jl"

# ╔═╡ 134fd792-064a-11eb-12dc-c5fe09197eba
md"### snippet 4.26"

# ╔═╡ 13549444-064a-11eb-08db-4fec5f434e9e
begin
	df = CSV.read(sr_datadir("Howell1.csv"), DataFrame)
	df = df[df.age .>= 18, :]
end;

# ╔═╡ 1363941c-064a-11eb-29c1-abbe84e92300
md"### snippet 4.27 - 4.29"

# ╔═╡ 13647ab2-064a-11eb-191c-73b11220323a
@model function ppl4_1(heights)
    μ ~ Normal(178, 20)
    σ ~ Uniform(0, 50)
    heights .~ Normal(μ, σ)
end

# ╔═╡ 137a6b92-064a-11eb-1b02-61006c9f3d8c
md"### snippet 4.30"

# ╔═╡ 136ffa54-064a-11eb-24a2-bbadb2cd9327
m4_1t = ppl4_1(df.height)

# ╔═╡ a452e82e-12fc-11eb-3f66-e96efc73e191
precis(m4_1t)

# ╔═╡ 1370b690-064a-11eb-01ff-ad478891bb93
md"##### Evaluating above cell might occasionally fail in Optim. Below cell should work."

# ╔═╡ 137b04a8-064a-11eb-3c7f-7959d865eda3
begin
	start = [mean(Normal(178, 20)), mean(Uniform(0, 50))]
	q4_1t = quap(m4_1t, start)
end

# ╔═╡ 13878f3e-064a-11eb-3d2d-833dab5fa561
md"### snippet 4.31"

# ╔═╡ 13886080-064a-11eb-014a-4badada265d5
@model function ppl4_2(heights)
    μ ~ Normal(178, 0.1)
    σ ~ Uniform(0, 50)
    heights .~ Normal(μ, σ)
end

# ╔═╡ 1394a660-064a-11eb-3b04-9fd3f1391d7f
begin
	m4_2t = ppl4_2(df.height)
	q4_2t = quap(m4_2t, NelderMead())
end

# ╔═╡ 3933a352-12fd-11eb-20bd-5bc9941a3fdc
typeof(q4_2t)

# ╔═╡ 5c64a180-12fd-11eb-2495-b7e35230241f
q4_1t.params

# ╔═╡ 1396d066-064a-11eb-16f6-41f3d6ed887e
md"### snippets 4.32, 4.33"

# ╔═╡ 13a1b544-064a-11eb-23eb-bd4cfd891d0f
q4_1t.vcov

# ╔═╡ 13a9efde-064a-11eb-00fe-99d6a9dae8c6
diag(q4_1t.vcov)

# ╔═╡ 78bb877e-064b-11eb-19e6-0f9da5cab17c
cov2cor(Matrix(q4_1t.vcov), sqrt.(diag(q4_1t.vcov)))

# ╔═╡ 13acf544-064a-11eb-2d81-3f22922e46ac
md"### snippets 4.34 - 4.36"

# ╔═╡ 13ba0c16-064a-11eb-0c6c-f19686f83ffd
begin
	res = rand(q4_1t.distr, 10_000)
	quap4_1t = DataFrame(res', ["μ", "σ"])
	Text(precis(quap4_1t; io=String))
end

# ╔═╡ 7cd372a2-12ff-11eb-1a05-bb7f9a95c1b7
begin
	chns4_1t = sample(m4_1t, NUTS(), 1000)
end

# ╔═╡ 13bbb6a6-064a-11eb-0658-99ca806b8a22
md"## End of clip-04-26-36t.jl"

# ╔═╡ Cell order:
# ╟─b32d14e8-0648-11eb-15d4-7d13643a1531
# ╠═13280ba4-064a-11eb-04dd-5b2c0598d84a
# ╠═134f547a-064a-11eb-0750-b590037632f5
# ╟─134fd792-064a-11eb-12dc-c5fe09197eba
# ╠═13549444-064a-11eb-08db-4fec5f434e9e
# ╟─1363941c-064a-11eb-29c1-abbe84e92300
# ╠═13647ab2-064a-11eb-191c-73b11220323a
# ╟─137a6b92-064a-11eb-1b02-61006c9f3d8c
# ╠═136ffa54-064a-11eb-24a2-bbadb2cd9327
# ╠═a452e82e-12fc-11eb-3f66-e96efc73e191
# ╟─1370b690-064a-11eb-01ff-ad478891bb93
# ╠═137b04a8-064a-11eb-3c7f-7959d865eda3
# ╟─13878f3e-064a-11eb-3d2d-833dab5fa561
# ╠═13886080-064a-11eb-014a-4badada265d5
# ╠═1394a660-064a-11eb-3b04-9fd3f1391d7f
# ╠═3933a352-12fd-11eb-20bd-5bc9941a3fdc
# ╠═5c64a180-12fd-11eb-2495-b7e35230241f
# ╟─1396d066-064a-11eb-16f6-41f3d6ed887e
# ╠═13a1b544-064a-11eb-23eb-bd4cfd891d0f
# ╠═13a9efde-064a-11eb-00fe-99d6a9dae8c6
# ╠═78bb877e-064b-11eb-19e6-0f9da5cab17c
# ╟─13acf544-064a-11eb-2d81-3f22922e46ac
# ╠═13ba0c16-064a-11eb-0c6c-f19686f83ffd
# ╠═7cd372a2-12ff-11eb-1a05-bb7f9a95c1b7
# ╟─13bbb6a6-064a-11eb-0658-99ca806b8a22