text:
using ManifoldsBase

using LinearAlgebra
using DoubleFloats
using ForwardDiff
using ReverseDiff
using StaticArrays
using Test

struct CustomDefinedRetraction <: ManifoldsBase.AbstractRetractionMethod end
struct CustomUndefinedRetraction <: ManifoldsBase.AbstractRetractionMethod end

function ManifoldsBase.injectivity_radius(
    ::ManifoldsBase.DefaultManifold,
    ::CustomDefinedRetraction,
)
    return 10.0
end

struct MatrixVectorTransport{T} <: AbstractVector{T}
    m::Matrix{T}
end

Base.getindex(x::MatrixVectorTransport, i) = x.m[:, i]

Base.size(x::MatrixVectorTransport) = (size(x.m, 2),)

@testset "Testing Default (Euclidean)" begin
    M = ManifoldsBase.DefaultManifold(3)
    types = [
        Vector{Float64},
        SizedVector{3,Float64},
        MVector{3,Float64},
        Vector{Float32},
        SizedVector{3,Float32},
        MVector{3,Float32},
        Vector{Double64},
        MVector{3,Double64},
        SizedVector{3,Double64},
    ]

    @test repr(M) == "DefaultManifold(3; field = ℝ)"
    @test isa(manifold_dimension(M), Integer)
    @test manifold_dimension(M) ≥ 0
    @test base_manifold(M) == M
    @test number_system(M) == ℝ
    @test ManifoldsBase.representation_size(M) == (3,)

    @test injectivity_radius(M) == Inf

    rm = ManifoldsBase.ExponentialRetraction()
    irm = ManifoldsBase.LogarithmicInverseRetraction()

    rm2 = CustomDefinedRetraction()
    rm3 = CustomUndefinedRetraction()

    for T in types
        @testset "Type $T" begin
            pts = convert.(Ref(T), [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])

            @test injectivity_radius(M, pts[1]) == Inf
            @test injectivity_radius(M, pts[1], rm) == Inf
            @test injectivity_radius(M, rm) == Inf
            @test injectivity_radius(M, rm2) == 10
            @test injectivity_radius(M, pts[1], rm2) == 10
            @test_throws ErrorException injectivity_radius(M, rm3)
            @test_throws ErrorException injectivity_radius(M, pts[1], rm3)

            tv1 = log(M, pts[1], pts[2])

            for pt in pts
                @test is_point(M, pt)
            end
            @test is_vector(M, pts[1], tv1; atol = eps(eltype(pts[1])))

            tv2 = log(M, pts[2], pts[1])
            @test isapprox(M, pts[2], exp(M, pts[1], tv1))
            @test isapprox(M, pts[1], exp(M, pts[1], tv1, 0))
            @test isapprox(M, pts[2], exp(M, pts[1], tv1, 1))
            @test isapprox(M, pts[1], exp(M, pts[2], tv2))
            @test is_point(M, retract(M, pts[1], tv1))
            @test isapprox(M, pts[1], retract(M, pts[1], tv1, 0))

            @test is_point(M, retract(M, pts[1], tv1, rm))
            @test isapprox(M, pts[1], retract(M, pts[1], tv1, 0, rm))

            new_pt = exp(M, pts[1], tv1)
            retract!(M, new_pt, pts[1], tv1)
            @test is_point(M, new_pt)
            for x in pts
                @test isapprox(M, zero_vector(M, x), log(M, x, x); atol = eps(eltype(x)))
                @test isapprox(
                    M,
                    zero_vector(M, x),
                    inverse_retract(M, x, x);
                    atol = eps(eltype(x)),
                )
                @test isapprox(
                    M,
                    zero_vector(M, x),
                    inverse_retract(M, x, x, irm);
                    atol = eps(eltype(x)),
                )
            end
            zero_vector!(M, tv1, pts[1])
            @test isapprox(M, pts[1], tv1, zero_vector(M, pts[1]))
            log!(M, tv1, pts[1], pts[2])
            @test norm(M, pts[1], tv1) ≈ sqrt(inner(M, pts[1], tv1, tv1))

            @test isapprox(M, exp(M, pts[1], tv1, 1), pts[2])
            @test isapprox(M, exp(M, pts[1], tv1, 0), pts[1])

            @test distance(M, pts[1], pts[2]) ≈ norm(M, pts[1], tv1)

            @test mid_point(M, pts[1], pts[2]) == [0.5, 0.5, 0.0]
            midp = allocate(pts[1])
            @test mid_point!(M, midp, pts[1], pts[2]) === midp
            @test midp == [0.5, 0.5, 0.0]

            @testset "Geodesic interface test" begin
                @test isapprox(M, geodesic(M, pts[1], tv1)(0.0), pts[1])
                @test isapprox(M, geodesic(M, pts[1], tv1)(1.0), pts[2])
                @test isapprox(M, geodesic(M, pts[1], tv1, 1.0), pts[2])
                @test isapprox(M, geodesic(M, pts[1], tv1, 1.0 / 2), (pts[1] + pts[2]) / 2)
                @test isapprox(M, shortest_geodesic(M, pts[1], pts[2])(0.0), pts[1])
                @test isapprox(M, shortest_geodesic(M, pts[1], pts[2])(1.0), pts[2])
                @test isapprox(M, shortest_geodesic(M, pts[1], pts[2], 0.0), pts[1])
                @test isapprox(M, shortest_geodesic(M, pts[1], pts[2], 1.0), pts[2])
                @test all(
                    isapprox.(
                        Ref(M),
                        geodesic(M, pts[1], tv1, [0.0, 1.0 / 2, 1.0]),
                        [pts[1], (pts[1] + pts[2]) / 2, pts[2]],
                    ),
                )
                @test all(
                    isapprox.(
                        Ref(M),
                        shortest_geodesic(M, pts[1], pts[2], [0.0, 1.0 / 2, 1.0]),
                        [pts[1], (pts[1] + pts[2]) / 2, pts[2]],
                    ),
                )
            end

            @testset "basic linear algebra in tangent space" begin
                @test isapprox(
                    M,
                    pts[1],
                    0 * tv1,
                    zero_vector(M, pts[1]);
                    atol = eps(eltype(pts[1])),
                )
                @test isapprox(M, pts[1], 2 * tv1, tv1 + tv1)
                @test isapprox(M, pts[1], 0 * tv1, tv1 - tv1)
                @test isapprox(M, pts[1], (-1) * tv1, -tv1)
            end

            @testset "Hat and vee in the tangent space" begin
                X = log(M, pts[1], pts[2])
                a = vee(M, pts[1], X)
                b = similar(a)
                vee!(M, b, pts[1], X)
                Y = hat(M, pts[1], a)
                Z = similar(Y)
                hat!(M, Z, pts[1], a)
                @test a == b
                @test X == Y
                @test Z == X
                @test a == vec(X)
            end

            @testset "broadcasted linear algebra in tangent space" begin
                @test isapprox(M, pts[1], 3 * tv1, 2 .* tv1 .+ tv1)
                @test isapprox(M, pts[1], -tv1, tv1 .- 2 .* tv1)
                @test isapprox(M, pts[1], -tv1, .-tv1)
                v = similar(tv1)
                v .= 2 .* tv1 .+ tv1
                @test v ≈ 3 * tv1
            end

            @testset "project test" begin
                # point
                @test isapprox(M, pts[1], project(M, pts[1]))
                pt = similar(pts[1])
                project!(M, pt, pts[1])
                @test isapprox(M, pt, pts[1])

                @test isapprox(M, pts[1], embed(M, pts[1]))
                pt = similar(pts[1])
                embed!(M, pt, pts[1])
                @test isapprox(M, pt, pts[1])

                # tangents
                @test isapprox(M, pts[1], tv1, project(M, pts[1], tv1))
                tv = similar(tv1)
                project!(M, tv, pts[1], tv1)
                @test isapprox(M, pts[1], tv, tv1)

                @test isapprox(M, pts[1], tv1, embed(M, pts[1], tv1))
                tv = similar(tv1)
                embed!(M, tv, pts[1], tv1)
                @test isapprox(M, pts[1], tv, tv1)
            end

            @testset "vector transport" begin
                # test constructor and alias
                @test DifferentiatedRetractionVectorTransport(ExponentialRetraction()) ==
                      ParallelTransport()
                v1 = log(M, pts[1], pts[2])
                v2 = log(M, pts[1], pts[3])
                v1t1 = vector_transport_to(M, pts[1], v1, pts[3])
                v1t2 = zero(v1t1)
                vector_transport_to!(M, v1t2, pts[1], v1, v2, ProjectionTransport())
                v1t3 = vector_transport_direction(M, pts[1], v1, v2)
                @test is_vector(M, pts[3], v1t1)
                @test is_vector(M, pts[3], v1t3)
                @test isapprox(M, pts[3], v1t1, v1t3)
                # along a `Vector` of points
                c = [pts[1]]
                v1t4 = vector_transport_along(M, pts[1], v1, c)
                @test isapprox(M, pts[1], v1, v1t4)
                v1t5 = allocate(v1)
                vector_transport_along!(M, v1t5, pts[1], v1, c)
                @test isapprox(M, pts[1], v1, v1t5)
                # along a custom type of points
                T = eltype(pts[1])
                c2 = MatrixVectorTransport{T}(reshape(pts[1], length(pts[1]), 1))
                v1t4c2 = vector_transport_along(M, pts[1], v1, c2)
                @test isapprox(M, pts[1], v1, v1t4c2)
                v1t5c2 = allocate(v1)
                vector_transport_along!(M, v1t5c2, pts[1], v1, c2)
                @test isapprox(M, pts[1], v1, v1t5c2)
                # On Euclidean Space Schild & Pole are identity
                @test vector_transport_to(
                    M,
                    pts[1],
                    v2,
                    pts[2],
                    SchildsLadderTransport(),
                ) == v2
                @test vector_transport_to(M, pts[1], v2, pts[2], PoleLadderTransport()) ==
                      v2
                @test vector_transport_to(
                    M,
                    pts[1],
                    v2,
                    pts[2],
                    ScaledVectorTransport(ParallelTransport()),
                ) == v2

                # along is also the identity
                c = [0.5 * (pts[1] + pts[2]), pts[2], 0.5 * (pts[2] + pts[3]), pts[3]]
                @test vector_transport_along(M, pts[1], v2, c, SchildsLadderTransport()) ==
                      v2
                @test vector_transport_along(M, pts[1], v2, c, PoleLadderTransport()) == v2
                @test vector_transport_along(M, pts[1], v2, c, ParallelTransport()) == v2
                # check mutating ones with defaults
                p = allocate(pts[1])
                ManifoldsBase.pole_ladder!(M, p, pts[1], pts[2], pts[3])
                # -log_p3 p == log_p1 p2
                @test isapprox(M, -log(M, pts[3], p), log(M, pts[1], pts[2]))
                ManifoldsBase.schilds_ladder!(M, p, pts[1], pts[2], pts[3])
                @test isapprox(M, log(M, pts[3], p), log(M, pts[1], pts[2]))

                @test repr(ParallelTransport()) == "ParallelTransport()"
                @test repr(ScaledVectorTransport(ParallelTransport())) ==
                      "ScaledVectorTransport(ParallelTransport())"
            end

            @testset "ForwardDiff support" begin
                exp_f(t) = distance(M, pts[1], exp(M, pts[1], t * tv1))
                d12 = distance(M, pts[1], pts[2])
                for t in 0.1:0.1:0.9
                    @test d12 ≈ ForwardDiff.derivative(exp_f, t)
                end

                retract_f(t) = distance(M, pts[1], retract(M, pts[1], t * tv1))
                for t in 0.1:0.1:0.9
                    @test ForwardDiff.derivative(retract_f, t) ≥ 0
                end
            end

            isa(pts[1], Union{Vector,SizedVector}) && @testset "ReverseDiff support" begin
                exp_f(t) = distance(M, pts[1], exp(M, pts[1], t[1] * tv1))
                d12 = distance(M, pts[1], pts[2])
                for t in 0.1:0.1:0.9
                    @test d12 ≈ ReverseDiff.gradient(exp_f, [t])[1]
                end

                retract_f(t) = distance(M, pts[1], retract(M, pts[1], t[1] * tv1))
                for t in 0.1:0.1:0.9
                    @test ReverseDiff.gradient(retract_f, [t])[1] ≥ 0
                end
            end
        end
    end

    @testset "mid_point on 0-index arrays" begin
        M = ManifoldsBase.DefaultManifold(1)
        p1 = fill(0.0)
        p2 = fill(1.0)
        @test isapprox(M, fill(0.5), mid_point(M, p1, p2))
    end

    @testset "Retracion" begin
        a = NLsolveInverseRetraction(ExponentialRetraction())
        @test a.retraction isa ExponentialRetraction
    end

    @testset "copy of points and vectors" begin
        M = DefaultManifold(2)
        p = [2.0, 3.0]
        q = similar(p)
        copyto!(M, q, p)
        @test p == q
        r = copy(M, p)
        @test r == p
        X = [4.0, 5.0]
        Y = similar(X)
        copyto!(M, Y, p, X)
        @test Y == X
        Z = copy(M, p, X)
        @test Z == X
    end
end


text:
using ChebyshevExpansions
using Test

@testset "ChebyshevExpansions.jl" begin

    # using Quantica:hamiltonian, LatticePresets, RegionPresets, LatticePresets.honeycomb,
    #     RegionPresets.circle, hopping, unitcell

    # using LinearAlgebra
    # h = hamiltonian(unitcell(LatticePresets.honeycomb(a0 = .2),
    #     region = RegionPresets.circle(10)), hopping(I))

    # hmat = h.harmonics[1].h
    # ket = ChebyshevExpansions.randomkets(eltype(hmat), eltype(hmat), size(hmat,1), 1)
    # momentaKPM(hmat)
    # momentaKPM(hmat, hmat, ket = 2)
    # densityKPM(hmat, hmat, ket = 2)
    # e,d = dosKPM(hmat, kets = randomkets(h,I,10), order = 1000, bandrange = (-3,3))

end

ext: jl

max_stars_repo_name: fernandopenaranda/ChebyshevExpansions

max_issues_repo_name: fernandopenaranda/ChebyshevExpansions

max_forks_repo_name: fernandopenaranda/ChebyshevExpansions


text:
# tutorial on conventional 2D SLAM example

# 0. Load required packages
using Revise
using RoME, Distributions
using DistributedFactorGraphs

# 1. Setting up our DistributedFactorGraph layers...
# Configure the mirror here.
# include(joinpath(ENV["HOME"],"Documents","blandauthlocal.jl"))
# backend_config, user_config = standardcloudgraphsetup(addrdict=addrdict)
# user_config["session"] = "SLAM2D_TUTORIAL"
# addrdict["robotId"] = "TestRobot"

# Configuring DFG with a local Graphs.jl instance and our mirror.
dfg = DistributedFactorGraph()
include("CloudDFGAPI.jl")
push!(dfg.mirrors, cloudDFG)
info(dfg)

####### ---- MAIN Example

# ... And away we go!
# fg = initfg()

# Add a first pose variable - x0
addV!(dfg, :x0, Pose2()) # , labels=["POSE"]
# Add a PriorPose2 to pin the first pose at a fixed location
addF!(dfg, [:x0], PriorPose2(MvNormal(zeros(3), 0.01*eye(3))) )

# ls(fg, :x0)
# isInitialized(fg, :x0)
# # to debug the zero lengths issue
# ensureAllInitialized!(fg)
# getVal(fg, :x0)
# Prior(MvNormal([0.0;0.0;0], diagm([1.0;1.0;0.01].^2)))

# Drive around in a hexagon
for i in 0:5
  psym = Symbol("x$i")
  nsym = Symbol("x$(i+1)")
  addVariable!(fg, nsym, Pose2) # , labels=["VARIABLE";"POSE"]
  addFactor!(fg, [psym;nsym], Pose2Pose2(reshape([10.0;0;pi/3],3,1), 0.01*eye(3), [1.0]), autoinit=true )
  # Pose2Pose2_NEW(MvNormal([10.0;0;pi/3], diagm([0.1;0.1;0.1].^2)))
en

# Graphs.plot(fg.g)
isInitialized(fg, :x5)
isInitialized(fg, :x6)

ensureAllInitialized!(fg)

using RoMEPlotting

drawPoses(fg)

getVal(fg, :x0)
v = getVert(fg, :x0)
getVal(v)

importall CloudGraphs

exvid = fg.IDs[:x0]
neoID = fg.cgIDs[exvid]
cvr = CloudGraphs.get_vertex(fg.cg, neoID, false)
exv = CloudGraphs.cloudVertex2ExVertex(cvr)
getVal(exv)
# getData(exv)

getData(fg.g.vertices[1])

tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree)
inferOverTreeR!(fg, tree)

using RoMEPlotting, Gadfly

pl = plotKDE(fg, [:x0; :x1; :x2; :x3; :x4; :x5; :x6]);

Gadfly.draw(PDF("tmpX0123456.pdf", 15cm, 15cm), pl)

@async run(`evince tmpX0123456.pdf`)

# pl = drawPoses(fg)
pl = drawPosesLandms(fg)
Gadfly.draw(PDF("tmpPosesFg.pdf", 15cm, 15cm), pl)
@async run(`evince tmpPosesFg.pdf`)

tree = wipeBuildNewTree!(fg)

@async Graphs.plot(tree.bt)

@time inferOverTree!(fg, tree)

# Graphs.plot(tree.bt)

@async Graphs.plot(fg.g)


text:
using Random

# generate a random number between 0 and 1
println(rand())

# select a 'random' element from a collection
for i in 1:5
    println(rand([1,2,3,4,5,6]))
end

# populate a collection with random values - specify 'type' and the 'number' of
# elements you want
a1 = rand(UInt16, 5)
println(typeof(a1))
println(a1)

# generate a 'random' set of characters
s1 = "World"
println(randstring(s1)) # note - default length of the output is 8
println(randstring(s1, 4)) # set the length of the output string
println(randstring('a':'z', 6)) # give it a range ...

# shuffle a collection
a2 = [1,2,3,4,5,6,7]
Random.shuffle!(a2)
println(a2)

# use a seed function
Random.seed!(100)
println(rand())
println(rand())

max_stars_repo_name: rpillar/JuliaExamples

max_issues_repo_name: rpillar/JuliaExamples

max_forks_repo_name: rpillar/JuliaExamples


text:
A = -5:0.2:5
B = -5:5
@testset "Indicator" begin
    @testset "above, nonstrict $y, $b" for y in A, b in B
        if y ≥ b
            @test indicator(y, b) == 1
        else
            @test indicator(y, b) == 0
        end
    end
    @testset "above, strict $y, $b" for y in A, b in B
        if y > b
            @test indicator(y, b, ineq=:strict) == 1
        else
            @test indicator(y, b, ineq=:strict) == 0
        end
    end
    @testset "below, nonstrict $y, $b" for y in A, b in B
        if y ≤ b
            @test indicator(y, b, indregion=:below) == 1
        else
            @test indicator(y, b, indregion=:below) == 0
        end
    end
    @testset "below, strict $y, $b" for y in A, b in B
        if y < b
            @test indicator(y, b, indregion=:below, ineq=:strict) == 1
        else
            @test indicator(y, b, indregion=:below, ineq=:strict) == 0
        end
    end
    @testset "Unknown input error checks" begin
        @test_throws ErrorException("Unknown indicator region: test") indicator(1, 0, indregion=:test)
        @test_throws ErrorException("Unknown inequality type: test") indicator(1, 0, ineq=:test)
    end
end


text:
###########################
#### Evapotranspiration ###
###########################

#' Potential Evapotranspiration
#'
#' Potential evapotranspiration according to Priestley & Taylor 1972 or
#'              the Penman-Monteith equation with a prescribed surface conductance.
#'
#' - data      DataFrame or matrix containing all required variables; optional
#' - Tair      Air temperature (degC)
#' - pressure  Atmospheric pressure (kPa)
#' - Rn        Net radiation (W m-2)
#' - G         Ground heat flux (W m-2); optional
#' - S         Sum of all storage fluxes (W m-2); optional
#' - VPD       Vapor pressure deficit (kPa); only used if `approach = Val(:PenmanMonteith)`.
#' - Ga        Aerodynamic conductance to heat/water vapor (m s-1); only used if `approach = Val(:PenmanMonteith)`.
#' - approach  Approach used. Either `Val(:PriestleyTaylor)` (default), or `Val(:PenmanMonteith)`.
#' - alpha     Priestley-Taylor coefficient; only used if `approach = Val(:PriestleyTaylor)`.
#' - Gs_pot    Potential/maximum surface conductance (mol m-2 s-1); defaults to 0.6 mol m-2 s-1;
#'                  only used if `approach = Val(:PenmanMonteith)`.
#' - missing_G_as_NA  if `true`, missing G are treated as `missing`s, otherwise set to 0.
#' - missing_S_as_NA  if `true`, missing S are treated as `missing`s, otherwise set to 0.
#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
#'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
#'                      See [`Esat_slope`](@ref).
#' - constants cp - specific heat of air for constant pressure (J K-1 kg-1)
#'                  eps - ratio of the molecular weight of water vapor to dry air
#'                  Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
#'                  Rd - gas constant of dry air (J kg-1 K-1) (only used if `approach = Val(:PenmanMonteith)`)
#'                  Rgas - universal gas constant (J mol-1 K-1) (only used if `approach = Val(:PenmanMonteith)`)
#'                  Kelvin - conversion degree Celsius to Kelvin (only used if `approach = Val(:PenmanMonteith)`)
#'
#' # Details
#'  Potential evapotranspiration is calculated according to Priestley & Taylor, 1972
#'          if `approach = Val(:PriestleyTaylor)` (the default):
#'
#'            ``LE_pot,PT = (\\alpha * \\Delta * (Rn - G - S)) / (\\Delta + \\gamma)``
#'
#'          ``\\alpha`` is the Priestley-Taylor coefficient, ``\\Delta`` is the slope
#'          of the saturation vapor pressure curve (kPa K-1), and ``\\gamma`` is the
#'          psychrometric constant (kPa K-1).
#'          if `approach = Val(:PenmanMonteith)`, potential evapotranspiration is calculated according
#'          to the Penman-Monteith equation:
#'
#'          ``LE_pot,PM = (\\Delta * (Rn - G - S) + \\rho * cp * VPD * Ga) / (\\Delta + \\gamma * (1 + Ga/Gs_pot)``
#'
#'          where ``\\Delta`` is the slope of the saturation vapor pressure curve (kPa K-1),
#'          ``\\rho`` is the air density (kg m-3), and ``\\gamma`` is the psychrometric constant (kPa K-1).
#'          The value of `Gs_pot` is typically a maximum value of Gs observed at the site, e.g. the 90th
#'          percentile of Gs within the growing season.
#'
#' # Value
#'  a DataFrame with the following columns:
#'         - ET_pot: Potential evapotranspiration (kg m-2 s-1)
#'         - LE_pot: Potential latent heat flux (W m-2)
#'
#' # Note
#' If the first argument `data` is provided (either a matrix or a DataFrame),
#'       the following variables can be provided as character (in which case they are interpreted as
#'       the column name of `data`) or as numeric vectors, in which case they are taken
#'       directly for the calculations. If `data` is not provided, all input variables have to be
#'       numeric vectors.
#'
#' #References
#' Priestley, C_H_B., Taylor, R_J., 1972: On the assessment of surface heat flux
#'             and evaporation using large-scale parameters. Monthly Weather Review 100, 81-92.
#'
#'             Allen, R_G., Pereira L_S., Raes D., Smith M., 1998: Crop evapotranspiration -
#'             Guidelines for computing crop water requirements - FAO Irrigation and drainage
#'             paper 56.
#'
#'             Novick, K_A., et al. 2016: The increasing importance of atmospheric demand
#'             for ecosystem water and carbon fluxes. Nature Climate Change 6, 1023 - 1027.
#'
#' #See also
#' [`surface_conductance`](@ref)
#'
#' ```@example; output = false
#' ```
#' # Calculate potential ET of a surface that receives a net radiation of 500 Wm-2
#' # using Priestley-Taylor:
#' potential_ET(Tair=30,pressure=100,Rn=500,alpha=1.26,approach=Val(:PriestleyTaylor))
#'
#' # Calculate potential ET for a surface with known Gs (0.5 mol m-2 s-1) and Ga (0.1 m s-1)
#' # using Penman-Monteith:
#' LE_pot_PM = potential_ET(Gs_pot=0.5,Tair=20,pressure=100,VPD=2,Ga=0.1,Rn=400,
#'                           approach=Val(:PenmanMonteith))[,"LE_pot"]
#' LE_pot_PM
#'
#' # now cross-check with the inverted equation
#' surface_conductance(Tair=20,pressure=100,VPD=2,Ga=0.1,Rn=400,LE=LE_pot_PM)
"""
"""
function potential_ET(data,Tair="Tair",pressure="pressure",Rn="Rn",G=nothing,S=nothing,
                         VPD="VPD",Ga="Ga_h",approach=c(Val(:PriestleyTaylor),Val(:PenmanMonteith)),
                         alpha=1.26,Gs_pot=0.6,missing_G_as_NA=false,missing_S_as_NA=false,
                         Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                         constants=BigleafConstants())

  approach = match_arg(approach)

  check_input(data,list(Tair,pressure,Rn,G,S))

  if(!isnothing(G))
    if (!missing_G_as_NA)
      G[ismissing(G)] = 0
end
else
    cat("Ground heat flux G is not provided and set to 0.",fill=true)
    G = 0
end

  if(!isnothing(S))
    if(!missing_S_as_NA)
      S[ismissing(S)] = 0
end
else
    cat("Energy storage fluxes S are not provided and set to 0.",fill=true)
    S = 0
end

  gamma  = psychrometric_constant(Tair,pressure,constants)
  Delta  = Esat_slope(Tair,Esat_formula,constants)[,"Delta"]

  if (approach == Val(:PriestleyTaylor))

    LE_pot = (alpha * Delta * (Rn - G - S)) / (Delta + gamma)
    ET_pot = LE_to_ET(LE_pot,Tair)

elseif (approach == Val(:PenmanMonteith))

    check_input(data,list(Gs_pot,VPD,Ga))

    Gs_pot = mol_to_ms(Gs_pot,Tair=Tair,pressure=pressure,constants=constants)
    rho    = air_density(Tair,pressure,constants)

    LE_pot = (Delta * (Rn - G - S) + rho * constants.cp * VPD * Ga) /
      (Delta + gamma * (1 + Ga / Gs_pot))
    ET_pot = LE_to_ET(LE_pot,Tair)
end

  return(DataFrame(ET_pot,LE_pot))

end

#' Reference Evapotranspiration
#'
#' Reference evapotranspiration calculated from the Penman-Monteith
#'              equation with a prescribed surface conductance.
#'              This function is deprecated. Use potential_ET(...,approach=Val(:PenmanMonteith)) instead.
#'
#' - data      DataFrame or matrix containing all required variables; optional
#' - Gs_ref    Reference surface conductance (m s-1); defaults to 0.0143 m s-1.
#' - Tair      Air temperature (degC)
#' - pressure  Atmospheric pressure (kPa)
#' - VPD       Vapor pressure deficit (kPa)
#' - Ga        Aerodynamic conductance to heat/water vapor (m s-1)
#' - Rn        Net radiation (W m-2)
#' - G         Ground heat flux (W m-2); optional
#' - S         Sum of all storage fluxes (W m-2); optional
#' - missing_G_as_NA  if `true`, missing G are treated as `missing`s, otherwise set to 0.
#' - missing_S_as_NA  if `true`, missing S are treated as `missing`s, otherwise set to 0.
#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
#'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
#'                      See [`Esat_slope`](@ref).
#' - constants cp - specific heat of air for constant pressure (J K-1 kg-1)
#'                  eps - ratio of the molecular weight of water vapor to dry air
#'                  Rd - gas constant of dry air (J kg-1 K-1) (only if `approach = Val(:PenmanMonteith)`)
#'                  Rgas - universal gas constant (J mol-1 K-1) (only if `approach = Val(:PenmanMonteith)`)
#'                  Kelvin - conversion degree Celsius to Kelvin (only if `approach = Val(:PenmanMonteith)`)
#'
#' @export
function reference_ET(data,Gs_ref=0.0143,Tair="Tair",pressure="pressure",VPD="VPD",Rn="Rn",Ga="Ga_h",
                         G=nothing,S=nothing,missing_G_as_NA=false,missing_S_as_NA=false,
                         Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                         constants=BigleafConstants())

  stop("this function is deprecated (since bigleaf version 0.6.0). For the calculation of potential ET from the Penman-Monteith equation (as formerly calculated with this function), use function potential_ET() with the argument approach='Penman-Monteith'. Note that the default value for argument 'Gs_pot' is expressed now in mol m-2 s-1 for simplicity (0.6 mol m-2 s-1).")

end

#' Equilibrium and Imposed Evapotranspiration
#'
#' Evapotranspiration (ET) split up into imposed ET and equilibrium ET.
#'
#' - data      DataFrame or matrix containing all required input variables
#' - Tair      Air temperature (deg C)
#' - pressure  Atmospheric pressure (kPa)
#' - VPD       Air vapor pressure deficit (kPa)
#' - Gs        surface conductance to water vapor (m s-1)
#' - Rn        Net radiation (W m-2)
#' - G         Ground heat flux (W m-2); optional
#' - S         Sum of all storage fluxes (W m-2); optional
#' - missing_G_as_NA  if `true`, missing G are treated as `missing`s, otherwise set to 0.
#' - missing_S_as_NA  if `true`, missing S are treated as `missing`s, otherwise set to 0.
#' - Esat_formula  Optional: Esat_formula to be used for the calculation of esat and the slope of esat.
#'                      One of `"Sonntag_1990"` (Default), `"Alduchov_1996"`, or `"Allen_1998"`.
#'                      See [`Esat_slope`](@ref).
#' - constants cp - specific heat of air for constant pressure (J K-1 kg-1)
#'                  eps - ratio of the molecular weight of water vapor to dry air (-)
#'                  Pa2kPa - conversion pascal (Pa) to kilopascal (kPa)
#'
#' # Details
 Total evapotranspiration can be written in the form (Jarvis & McNaughton 1986):
#'
#'            ``ET = \\Omega ET_eq + (1 - \\Omega)ET_imp``
#'
#'          where ``\\Omega`` is the decoupling coefficient as calculated from
#'          [`decoupling`](@ref). `ET_eq` is the equilibrium evapotranspiration rate,
#'          the ET rate that would occur under uncoupled conditions, where the heat budget
#'          is dominated by radiation (when Ga -> 0):
#'
#'            ``ET_eq = (\\Delta * (Rn - G - S) * \\lambda) / (\\Delta + \\gamma)``
#'
#'          where ``\\Delta`` is the slope of the saturation vapor pressure curve (kPa K-1),
#'          ``\\lambda`` is the latent heat of vaporization (J kg-1), and ``\\gamma``
#'          is the psychrometric constant (kPa K-1).
#'          `ET_imp` is the imposed evapotranspiration rate, the ET rate
#'          that would occur under fully coupled conditions (when Ga -> inf):
#'
#'            ``ET_imp = (\\rho * cp * VPD * Gs * \\lambda) / \\gamma``
#'
#'          where ``\\rho`` is the air density (kg m-3).
#'
#' # Note
#' Surface conductance (Gs) can be calculated with [`surface_conductance`](@ref).
#'       Aerodynamic conductance (Ga) can be calculated using [`aerodynamic_conductance!`](@ref).
#'
#' # Value
 A DataFrame with the following columns:
#'         - ET_eq: Equilibrium ET (kg m-2 s-1)
#'         - ET_imp: Imposed ET (kg m-2 s-1)
#'         - LE_eq: Equilibrium LE (W m-2)
#'         - LE_imp: Imposed LE (W m-2)
#'
#' #References
#' Jarvis, P_G., McNaughton, K_G., 1986: Stomatal control of transpiration:
#'             scaling up from leaf to region. Advances in Ecological Research 15, 1-49.
#'
#'             Monteith, J_L., Unsworth, M_H., 2008: Principles of Environmental Physics.
#'             3rd edition. Academic Press, London.
#'
#' #See also
#' [`decoupling`](@ref)
#'
#' ```@example; output = false
#' ```
#' df = DataFrame(Tair=20,pressure=100,VPD=seq(0.5,4,0.5),
#'                  Gs_ms=seq(0.01,0.002,length_out=8),Rn=seq(50,400,50))
#' equilibrium_imposed_ET(df)
#'
"""
"""
function equilibrium_imposed_ET(data,Tair="Tair",pressure="pressure",VPD="VPD",Gs="Gs_ms",
                                   Rn="Rn",G=nothing,S=nothing,missing_G_as_NA=false,missing_S_as_NA=false,
                                   Esat_formula=c("Sonntag_1990","Alduchov_1996","Allen_1998"),
                                   constants=BigleafConstants())

  check_input(data,list(Tair,pressure,VPD,Rn,Gs,G,S))

  if(!isnothing(G))
    if (!missing_G_as_NA){G[ismissing(G)] = 0}
else
    cat("Ground heat flux G is not provided and set to 0.",fill=true)
    G = 0
end

  if(!isnothing(S))
    if(!missing_S_as_NA){S[ismissing(S)] = 0 }
else
    cat("Energy storage fluxes S are not provided and set to 0.",fill=true)
    S = 0
end

  rho    = air_density(Tair,pressure,constants)
  gamma  = psychrometric_constant(Tair,pressure,constants)
  Delta  = Esat_slope(Tair,Esat_formula,constants)[,"Delta"]

  LE_eq  = (Delta * (Rn - G - S)) / (gamma + Delta)
  LE_imp = (rho * constants.cp * Gs * VPD) / gamma

  ET_imp = LE_to_ET(LE_imp,Tair)
  ET_eq  = LE_to_ET(LE_eq,Tair)

  return(DataFrame(ET_eq,ET_imp,LE_eq,LE_imp))
end


text:
using Test
using CartesianGeneticProgramming
import YAML

@testset "CGPInd construction" begin
    cfg = get_config("test.yaml")
    ind = CGPInd(cfg)

    @test length(ind.nodes) == 3 * cfg.columns + cfg.n_in
    for node in ind.nodes
        if node.active
            @test node.x >= 1
            @test node.x <= length(ind.nodes)
            @test node.y >= 1
            @test node.y <= length(ind.nodes)
        end
        # test that stringifying works
        @test typeof(string(node)) == String
    end

    @test size(ind.genes) == (cfg.rows, cfg.columns, 3)
end

"""
using random values, sort individuals that are different
"""
function select_random(pop::Array{CGPInd}, elite::Int; n_in=113, n_sample=100)
    actions = zeros(Int, length(pop))
    dists = zeros(n_sample, length(pop))
    inputs = rand(n_in, n_sample)

    for i in 1:n_sample
        for j in eachindex(pop)
            actions[j] = argmax(process(pop[j], inputs[:, i]))
        end
        for j in eachindex(pop)
            dists[i, j] = sum(actions[j] .!= actions)
        end
    end
    d = sum(dists, dims=1)[:]
    ds = sortperm(d)[1:elite]
    pop[ds]
end

@testset "Node genes" begin
    cfg = get_config("test.yaml")
    ind = CGPInd(cfg)
    ind2 = CGPInd(cfg)
    @test any(ind.chromosome .!= ind2.chromosome)

    genes = get_genes(ind, 1)
    @test all(genes .== 0)

    genes = get_genes(ind, ind.n_in+1)
    @test all(genes .>= 0.0)
    @test all(genes .<= 1.0)

    set_genes!(ind2, ind.n_in+1, genes)
    @test all(ind.chromosome[1:3] .== ind2.chromosome[1:3])

    for i in 1:length(ind.nodes)
        genes = get_genes(ind, i)
        set_genes!(ind2, i, genes)
    end
    o = length(ind.chromosome) - cfg.n_out
    @test all(ind.chromosome[1:o] .== ind2.chromosome[1:o])

    all_genes = get_genes(ind, collect((ind.n_in+1):length(ind.nodes)))
    @test all(ind.chromosome[1:o] .== all_genes)
end

@testset "Processing" begin
    cfg = get_config("test.yaml"; functions=["f_abs", "f_add", "f_mult"])
    ind = CGPInd(cfg)

    # test that f(0, 0, 0, 0) = 0
    inputs = zeros(cfg.n_in)
    set_inputs(ind, inputs)
    for i in 1:cfg.n_in
        @test ind.buffer[i] == 0.0
    end
    output = process(ind)
    @test output[1] == 0.0
    for i in eachindex(ind.nodes)
        if ind.nodes[i].active
            @test ind.buffer[i] == 0.0
        end
    end

    # test that f(1, 1, 1, 1) = 1
    for i in eachindex(ind.nodes)
        ind.buffer[i] = 1.0 # requires that buffer is 1
    end
    output = process(ind, ones(cfg.n_in))
    @test output[1] == 1.0
    for i in eachindex(ind.nodes)
        if ind.nodes[i].active
            @test ind.buffer[i] == 1.0
        end
    end

    cfg = get_config("test.yaml"; functions=["f_abs", "f_add", "f_mult"], recur=1.0)
    ind = CGPInd(cfg)
    output = process(ind, rand(cfg.n_in))
    @test output[1] <= 1.0 && output[1] >= -1.0
    for i in eachindex(ind.nodes)
        if ind.nodes[i].active
            @test ind.buffer[i] <= 1.0 && ind.buffer[i] >= -1.0
        end
    end

    pop = [CGPInd(cfg) for i in 1:10]
    sp = select_random(pop, 2; n_in=cfg.n_in, n_sample=5)
    @test length(sp) == 2
    @test sp[1].buffer[1] != 0.0

    reset!(ind)
    @test all(ind.buffer .== 0.0)

    f_conns = forward_connections(ind)
    @test any(map(x->length(x)>1, f_conns))

    ot = get_output_trace(ind, 1)
    @test length(ot) > 0

    all_traces = get_output_trace(ind)
    @test issubset(ot, all_traces)
end


text:
#!/usr/bin/env julia
# A very thin wrapper to start a GeMM simulation. If you need something fancier,
# have a look at `rungemmparallel.jl`, or import the GeMM module into your own
# wrapper script.

thisDir = joinpath(pwd(), "src")
any(path -> path == thisDir, LOAD_PATH) || push!(LOAD_PATH, thisDir)
using Pkg
Pkg.activate(".")
using GeMM

using Profile
Profile.clear()

rm("results/taita_test", recursive=true, force=true)

@profile rungemm("zosterops.config")

open("profile_flat.txt", "w") do s
    Profile.print(IOContext(s, :displaysize=>(300,145)), format=:flat, mincount=10, sortedby=:count)
end

open("profile_tree.txt", "w") do s
    Profile.print(IOContext(s, :displaysize=>(300,300)), mincount=10)
end

max_stars_repo_name: CCTB-Ecomods/gemm

max_issues_repo_name: CCTB-Ecomods/gemm

max_forks_repo_name: CCTB-Ecomods/gemm


text:
using Behavior
using CoffeeMachine
using Test

@test runspec(pkgdir(CoffeeMachine))


text:
############################################################################################
# Models to be used in construction
objectives = [
    Objective(ModelWrapper(MyBaseModel(), myparameter1, FlattenDefault()), data_uv),
    Objective(ModelWrapper(MyBaseModel(), myparameter1, FlattenDefault(; output = Float32)), data_uv)
]
Nchains = 4
tempermethods = [
    IterationTempering(Float64, UpdateFalse(), 1.0, 1000),
    IterationTempering(Float64, UpdateTrue(), 1.0, 1000),
    JointTempering(Float64, UpdateFalse(), .5, Float64(Nchains), Nchains),
    JointTempering(Float64, UpdateTrue(), .5, Float64(Nchains), Nchains)
]

#=
iter = 2
tempermethod = tempermethods[iter]
=#
############################################################################################
@testset "Sampling, type conversion" begin
    for tempermethod in tempermethods
        for iter in eachindex(objectives)
                sampledefault = SampleDefault(;
                    dataformat=Batch(),
                    tempering=deepcopy(tempermethod), #IterationTempering(Float64, UpdateFalse(), 1.0, 1000),
                    chains=4,
                    iterations=100,
                    burnin=max(1, Int64(floor(10/10))),
                    thinning = 1,
                    safeoutput=false,
                    printoutput=false,
                    printdefault=PrintDefault(),
                    report=ProgressReport(;
                        bar=false,
                        log=SilentLog()
                    ),
                )
                temperupdate = sampledefault.tempering.adaption
                _obj = deepcopy(objectives[iter])
                _flattentype = _obj.model.info.reconstruct.default.output
        #MCMC
                ## Assign kernels
                mcmc = MCMC(NUTS,(:μ, :σ,); stepsize = ConfigStepsize(;stepsizeadaption = UpdateFalse()))
                trace, algorithms = sample(_rng, _obj.model, _obj.data, mcmc ; default = deepcopy(sampledefault))
                ## If single mcmc kernel assigned, can capture previous results
                @test isa(trace.info.sampling.captured, typeof(temperupdate))
                ## Continue sampling
                trace2, algorithms2 = sample!(100, _rng, _obj.model, _obj.data, trace, algorithms)
                @test isa(trace2.info.sampling.captured, typeof(temperupdate))
                postmean = trace_to_posteriormean(trace, _obj.model, _obj.tagged, 0, 1)
        #SMC
                ibis = SMCConstructor(mcmc, SMCDefault(jitterthreshold=0.99, resamplingthreshold=1.0))
                trace, algorithms = sample(_rng, _obj.model, _obj.data, ibis; default = deepcopy(sampledefault))
                ## If single mcmc kernel assigned, can capture previous results
                @test isa(trace.info.sampling.captured, UpdateFalse)
                ## Continue sampling
                newdat = randn(_rng, length(_obj.data)+100)
                trace2, algorithms2 = sample!(100, _rng, _obj.model, newdat, trace, algorithms)
                @test isa(trace2.info.sampling.captured, UpdateFalse)
        # Combinations
                trace, algorithms = sample(_rng, _obj.model, _obj.data, mcmc, ibis; default = deepcopy(sampledefault))
                ## If single mcmc kernel assigned, can capture previous results
                @test isa(trace.info.sampling.captured, UpdateTrue)
                ## Continue sampling
                newdat = randn(_rng, length(_obj.data)+100)
                trace2, algorithms2 = sample!(100, _rng, _obj.model, newdat, trace, algorithms)
                @test isa(trace2.info.sampling.captured, UpdateTrue)
        end
    end
end

############################################################################################
chains = [1, Nchains, 1, Nchains]
temperchainmethods = [
    IterationTempering(Float64, UpdateTrue(), float(chains[1]), 1000),
    IterationTempering(Float64, UpdateTrue(), float(chains[2]), 1000),
    JointTempering(Float64, UpdateTrue(), .5, float(chains[1]), chains[1]),
    JointTempering(Float64, UpdateTrue(), .5, float(chains[2]), chains[2])
]

@testset "Sampling, Chain Management" begin
    for (iter, tempermethod) in enumerate(temperchainmethods)
            ## Set SampleDefault
            sampledefault = SampleDefault(;
                dataformat=Batch(),
                tempering=deepcopy(tempermethod),
                chains=chains[iter],
                iterations=500,
                burnin=0,
                thinning = 1,
                safeoutput=false,
                printoutput=true,
                printdefault=PrintDefault(),
                report=ProgressReport(;
                    bar=false,
                    log=SilentLog()
                ),
            )
            temperupdate = sampledefault.tempering.adaption
            _obj = deepcopy(objectives[1])
    #MCMC
            ## Assign kernels
            mcmc = MCMC(NUTS,(:μ, :σ,); stepsize = ConfigStepsize(;stepsizeadaption = UpdateFalse()))
            trace, algorithms = sample(_rng, _obj.model, _obj.data, mcmc ; default = deepcopy(sampledefault))
            ## If single mcmc kernel assigned, can capture previous results
            @test isa(trace.info.sampling.captured, typeof(temperupdate))
            ## Continue sampling
            trace2, algorithms2 = sample!(100, _rng, _obj.model, _obj.data, trace, algorithms)
            @test isa(trace2.info.sampling.captured, typeof(temperupdate))
    # Combinations
            mcmc = MCMC(NUTS,(:μ, :σ,); stepsize = ConfigStepsize(;stepsizeadaption = UpdateFalse()))
            trace, algorithms = sample(_rng, _obj.model, _obj.data, mcmc, mcmc ; default = deepcopy(sampledefault))
            ## If single mcmc kernel assigned, can capture previous results
            @test isa(trace.info.sampling.captured, typeof(temperupdate))
            ## Continue sampling
            trace2, algorithms2 = sample!(100, _rng, _obj.model, _obj.data, trace, algorithms)
            @test isa(trace2.info.sampling.captured, typeof(temperupdate))
    end
end

############################################################################################
## Make model for several parameter types
@testset "Sampling, printing and ConsoleLog" begin
    for (iter, tempermethod) in enumerate(tempermethods)
            ## Set SampleDefault
            sampledefault = SampleDefault(;
                dataformat=Batch(),
                tempering=deepcopy(tempermethod),
                chains=4,
                iterations=100,
                burnin=max(1, Int64(floor(100/10))),
                thinning = 1,
                safeoutput=false,
                printoutput=true,
                printdefault=PrintDefault(),
                report=ProgressReport(;
                    bar=true,
                    log=ConsoleLog()
                ),
            )
            temperupdate = sampledefault.tempering.adaption
            _obj = deepcopy(objectives[1])
    #MCMC
            ## Assign kernels
            mcmc = MCMC(NUTS,(:μ, :σ,); stepsize = ConfigStepsize(;stepsizeadaption = UpdateFalse()))
            trace, algorithms = sample(_rng, _obj.model, _obj.data, mcmc ; default = deepcopy(sampledefault))
            ## If single mcmc kernel assigned, can capture previous results
            @test isa(trace.info.sampling.captured, typeof(temperupdate))
            ## Continue sampling
            trace2, algorithms2 = sample!(100, _rng, _obj.model, _obj.data, trace, algorithms)
            @test isa(trace2.info.sampling.captured, typeof(temperupdate))
    #SMC
            ibis = SMCConstructor(mcmc, SMCDefault(jitterthreshold=0.99, resamplingthreshold=1.0))
            trace, algorithms = sample(_rng, _obj.model, _obj.data, ibis; default = deepcopy(sampledefault))
            ## If single mcmc kernel assigned, can capture previous results
            @test isa(trace.info.sampling.captured, UpdateFalse)
            ## Continue sampling
            newdat = randn(_rng, length(_obj.data)+100)
            trace2, algorithms2 = sample!(100, _rng, _obj.model, newdat, trace, algorithms)
            @test isa(trace2.info.sampling.captured, UpdateFalse)
    # Combinations
            trace, algorithms = sample(_rng, _obj.model, _obj.data, mcmc, ibis; default = deepcopy(sampledefault))
            ## If single mcmc kernel assigned, can capture previous results
            @test isa(trace.info.sampling.captured, UpdateTrue)
            ## Continue sampling
            newdat = randn(_rng, length(_obj.data)+100)
            trace2, algorithms2 = sample!(100, _rng, _obj.model, newdat, trace, algorithms)
            @test isa(trace2.info.sampling.captured, UpdateTrue)
    end
end

############################################################################################
#SMC
@testset "Sampling, Sequential Estimation" begin
    for (iter, tempermethod) in enumerate(tempermethods)
            ## Set SampleDefault
            sampledefault = SampleDefault(;
                dataformat=Expanding(100),
                tempering= deepcopy(tempermethod),
                chains=4,
                iterations=500,
                burnin=0,
                thinning = 1,
                safeoutput=false,
                printoutput=true,
                printdefault=PrintDefault(),
                report=ProgressReport(;
                    bar=false,
                    log=SilentLog()
                ),
            )
            temperupdate = sampledefault.tempering.adaption
            _obj = deepcopy(objectives[1])
    #MCMC
            ## Assign kernels
            mcmc = MCMC(NUTS,(:μ, :σ,); stepsize = ConfigStepsize(;stepsizeadaption = UpdateFalse()))
            trace, algorithms = sample(_rng, _obj.model, _obj.data, mcmc ; default = deepcopy(sampledefault))
            ## If single mcmc kernel assigned, can capture previous results
            #!NOTE: If Expanding/Increasing data, update always true
            @test isa(trace.info.sampling.captured, UpdateTrue)
            ## Continue sampling
            trace2, algorithms2 = sample!(500, _rng, _obj.model, _obj.data, trace, algorithms)
            @test isa(trace2.info.sampling.captured, UpdateTrue)
    #SMC
            ibis = SMCConstructor(mcmc, SMCDefault(jitterthreshold=0.99, resamplingthreshold=1.0))
            trace, algorithms = sample(_rng, _obj.model, _obj.data, ibis; default = deepcopy(sampledefault))
            ## If single mcmc kernel assigned, can capture previous results
            @test isa(trace.info.sampling.captured, UpdateFalse)
            ## Continue sampling
            newdat = randn(_rng, length(_obj.data)+500)
            trace2, algorithms2 = sample!(500, _rng, _obj.model, newdat, trace, algorithms)
            @test isa(trace2.info.sampling.captured, UpdateFalse)
    # SMC2
        _obj = deepcopy(myobjective_mcmc)
        _tagged_pf = myobjective_pf.tagged
        _tagged_mcmc = myobjective_mcmc.tagged
        pfdefault_propagate = ParticleFilterDefault(referencing = Marginal(),)
        pfdefault_pmcmc = ParticleFilterDefault(referencing = Ancestral(),)
        mcmcdefault = MCMCDefault(; stepsize = ConfigStepsize(;stepsizeadaption = UpdateFalse()))
        pmcmc = ParticleGibbs(
            ParticleFilterConstructor(keys(_tagged_pf.parameter), pfdefault_pmcmc),
            MCMCConstructor(NUTS, keys(_tagged_mcmc.parameter), mcmcdefault)
        )
        pf = ParticleFilterConstructor(keys(_tagged_pf.parameter), pfdefault_propagate)
        smc2 = SMCConstructor(SMC2Constructor(pf, pmcmc), SMCDefault(jitterthreshold=0.75, resamplingthreshold=0.75))
        trace, algorithms = sample(_rng, _obj.model, data, smc2; default = deepcopy(sampledefault))
        ## If single mcmc kernel assigned, can capture previous results
        @test isa(trace.info.sampling.captured, UpdateFalse)
        #Check if correct parameter are printed
        allparam, printparam = Baytes.showparam(_obj.model, trace.info.datatune, smc2)
        @test allparam == keys(_obj.model.val)
        @test printparam == (:μ, :σ, :p)
        ## Continue sampling
        newdat = randn(_rng, length(data)+500)
        trace2, algorithms2 = sample!(500, _rng, _obj.model, newdat, trace, algorithms)
        @test isa(trace2.info.sampling.captured, UpdateFalse)
    end
end

############################################################################################
#Utility
@testset "Utility, maxiterations" begin
    Nstart = 10^2
    Niter1 = 10^4
    Niter2 = 10^1
    Ndata = length(data) - Nstart

    datatune = DataTune(data, Expanding(Nstart))
    @test Baytes.maxiterations(datatune, Niter1) == Ndata
    @test Baytes.maxiterations(datatune, Niter2) == Niter2
end


text:
include(joinpath(splitdir(@__FILE__())[1],"EvoqusPreamble.jl"))

#Run the script
println("==========================================================================")
println("Beginning execution of script GMHPhotoReceptor.jl/evoqus/PhotoReceptor1.jl")
println("==========================================================================")

println("Number of parallel processes running: ",nprocs())

@everywhere using GeneralizedMetropolisHastings
@everywhere using GMHPhotoReceptor

###Print a message indicating that the GMH package has loaded correctly
print_gmh_module_loaded()

println("================================")
println("Initialize Simulation Parameters")
println("================================")

#Standard M-H for nproposals == 1
#Generalized M-H for nproposals > 1
nproposals = 100000

#MCMC iteration specifications
nburnin = 0
niterations = 20
ntunerperiod = 5

###Values of the model
numvilli1 = 30000

#specify the values that determine the priors on the parameters
latencylocation = (2.0,3.5) #uniform distribution with (low,high) values
latencyscale = (0.2,0.7) #uniform distribution with (low,high) values
refractorylocation = (4.0,6.0) #uniform distribution with (low,high) values
refractoryscale = (1.5,2.5) #uniform distribution with (low,high) values
bumpamplitude = (3.0,5.0) #uniform distribution with (low,high) values
bumpshape = (log(3.0),0.3) #lognormal distribution with (location,scale) values; variable can vary roughly between 2.0 and 4.0, but becomes increasingly penalized outside
bumpscale = (log(2.5),0.3) #lognormal distribution with (location,scale) values; variable can vary roughly between 1.5 and 3.5, but becomes increasingly penalized outside

photonfilename = joinpath(splitdir(@__FILE__())[1],"../data/naturallight.jld")
photons1 = photonsequence(photonfilename)
current1 = lightinducedcurrent(photonfilename)

modelpolicy1 = policy(:photoreceptor) #4-parameter model with stochastic lognormal latency and refractory parameters and fixed bump parameters
params1 = parameters(:photoreceptor,modelpolicy1,latencylocation,latencyscale,refractorylocation,refractoryscale)

####Variance for the noise model, estimated from previous runs
variance1 = [3600.0]

println("==========================================")
println("Simulation parameters defined successfully")
println("==========================================")

###Create a PhotoReceptor model
model1 = model(:photoreceptor,params1,photons1,current1,variance1,numvilli1,modelpolicy1)

###Show the model
println("==========================")
println("Model defined successfully")
println("==========================")
show(model1)

###Create a Metropolis sampler with a Normal proposal density
propcov = [0.1,0.01,0.1,0.01]
sampler1 = sampler(:mh,:normal,0.1,propcov)
println("============================")
println("Sampler defined successfully")
println("============================")
show(sampler1)

###Create a tuner that scales the proposal density
#tuner1 = tuner(:scale,ntunerperiod,0.5,:erf)
tuner1 = tuner(:scale,ntunerperiod,0.5,:erf)
println("==========================")
println("Tuner defined successfully")
println("==========================")
show(tuner1)

###Create a Generalized Metropolis-Hastings runner (which will default to Standard MH when nproposals=1)
runnerpolicy1 = policy(:mh,nproposals;initialize=:prior)
runner1 = runner(runnerpolicy1,niterations,nproposals;numburnin = nburnin)
println("===========================")
println("Runner defined successfully")
println("===========================")
show(runner1)

###Run the MCMC (can take quite a bit of time)
println("=======================")
println("Run the MCMC simulation")
println("=======================")
@time chain1 = run!(runner1,model1,sampler1,tuner1)
println("=========================")
println("Completed MCMC simulation")
println("=========================")

###Show the result of the simulations
show(chain1)

nparas = numparas(model1)
meanparamvals = mean(samples(chain1),2)
stdparamvals = std(samples(chain1),2)

println("Results of the MCMC simulation:")
for i=1:nparas
    println("mean $(parameters(model1)[i].key): $(meanparamvals[i])")
    println("std $(parameters(model1)[i].key): $(stdparamvals[i])")
end


text:
function test_exact_forecast(metric; n = 1000, seed::Int = 1)
    Random.seed!(seed)
    a = rand(n)
    @test metric(a, a) == 0
    return
end


text:
using ForwardDiff, Test
using MixedModeBroadcastAD: unwrt

include("kernels.jl")

@testset "hmlstm kernels" begin
    dims = 128
    for usegpu in (false, true)
        println("testing hmlstm kernels for usegpu=", usegpu)
        tfkernel!, _, tfderivs, tfbuffers = get_hmlstm_kernel(true, usegpu, dims)
        kernel!, inputs, derivs, buffers = get_hmlstm_kernel(false, usegpu, dims)
        kernel!(inputs, derivs, buffers)
        tfkernel!(unwrt.(inputs), tfderivs, tfbuffers)
        @test Array(derivs[1]) ≈ Array(tfderivs[1])
        @test Array(derivs[2]) ≈ Array(tfderivs[2])
        @test Array(derivs[3]) ≈ Array(tfderivs[3])
        @test Array(derivs[4]) ≈ Array(tfderivs[4])
        cpu_inputs = Array.(unwrt.(inputs))
        deriv_i = 1
        for i in 1:length(inputs)
            if isa(inputs[i], Wrt)
                println("\t...checking gradient for input $i")
                cpu_kernel_i = x -> begin
                    before = cpu_inputs[1:(i - 1)]
                    after = cpu_inputs[(i + 1):end]
                    return sum(broadcast(cpu_hmlstm_update_c_scalar, before..., x, after...))
                end
                @test Array(derivs[deriv_i]) ≈ ForwardDiff.gradient(cpu_kernel_i, cpu_inputs[i])
                deriv_i += 1
            end
        end
    end
end

@testset "arity scaling kernels" begin
    dims = 128
    for arity in 1:3, usegpu in (false, true)
        println("testing arity scaling kernel for usegpu=", usegpu)
        kernel!, inputs, derivs, buffers = get_arity_scaling_kernel(usegpu, dims, arity)
        kernel!(inputs, derivs, buffers)
        cpu_inputs = Array.(unwrt.(inputs))
        for (i, cpu_input) in enumerate(cpu_inputs)
            println("\t...checking gradient for input $i")
            cpu_kernel_i = x -> begin
                before = cpu_inputs[1:(i - 1)]
                after = cpu_inputs[(i + 1):end]
                return sum(broadcast(cpu_arity_scaling, before..., x, after...))
            end
            @test Array(derivs[i]) ≈ ForwardDiff.gradient(cpu_kernel_i, cpu_input)
        end
    end
end


text:
using PushVectors, Test, BenchmarkTools

@testset "implementation sanity checks" begin

    # create
    v = PushVector{Int}()
    @test v.len == 0
    @test v.parent isa Vector{Int}
    @test v == Int[]

    # add one element
    push!(v, 1)
    @test @inferred(v[1]) == 1
    @test v.len == 1 == length(v)
    @test size(v) == (1, )
    @test v == Int[1]

    # empty
    @test empty!(v) == Int[] == v

    # make it double
    for i in 1:5
        push!(v, i)
    end
    @test v == 1:5
    @test v.len == 5 == length(v)
    @test size(v) == (5, )
    @test length(v.parent) > 4
    @test @inferred setindex!(v, 9, 3) == 9
    w = [1, 2, 9, 4, 5]

    # sizehint!
    sizehint!(v, 20)            # up
    @test length(v.parent) == 20
    @test v == w

    sizehint!(v, 10)            # down
    @test v == w
    @test length(v.parent) == 10

    sizehint!(v, 3)             # ignored
    @test v == w
    @test length(v.parent) == 10

    @test @inferred(finish!(v)) == w
end

@testset "append!" begin
    v = Vector{Float64}()
    pv = PushVector{Float64}()
    for _ in 1:1000
        z = randn(rand(1:10))
        @test append!(v, z) == @inferred append!(pv, z)
        @test v == pv
    end
end

function pushit!(v)
    for i in 1:10^4
        push!(v, i)
    end
end

function appendit!(v, cycled)
    n = length(cycled)
    for i in 1:1000
        append!(v, cycled[(i % n) + 1])
    end
end

@testset "relative benchmarking" begin
    T_PushVector = @belapsed begin
        p = PushVector{Int64}()
        pushit!(p)
        finish!(p)
    end

    T_Vector = @belapsed begin
        p = Vector{Int64}()
        pushit!(p)
    end

    cycled = [randn(i) for i in 1:5:(5*17)]

    A_PushVector = @belapsed begin
        p = PushVector{Float64}()
        appendit!(p, $cycled)
        finish!(p)
    end

    A_Vector = @belapsed begin
        p = Vector{Float64}()
        appendit!(p, $cycled)
    end

    @info "benchmarks" T_PushVector T_Vector push!_ratio = round(T_PushVector / T_Vector; digits = 2)
    @test T_PushVector ≤ T_Vector
    @info "benchmarks" A_PushVector A_Vector append!_ratio = round(A_PushVector / A_Vector; digits = 2)
    # here just ensure that it is not much worse, testing is noisy
    # @test A_PushVector ≤ A_Vector * 1.5
end


text:
function _verify_indexorder(network, indexorder)
    allindices = Vector{Any}()
    for ind in network
        append!(allindices, ind)
    end
    while length(allindices) > 0
        i = pop!(allindices)
        if i ∈ indexorder # contraction index
            k = findfirst(isequal(i), allindices)
            k === nothing && return false
            l = findnext(isequal(i), allindices, k+1)
            l !== nothing && return false
            deleteat!(allindices, k)
        else
            findfirst(isequal(i), allindices) === nothing || return false
        end
    end
    return true
end

function indexordertree(network, indexorder)
    _verify_indexorder(network, indexorder) ||
        throw(ArgumentError("not a valid contraction tree with contraction indices $indexorder"))
    contractionindices = Vector{Vector{Any}}(undef, length(network))
    for k = 1:length(network)
        indices = network[k]
        # trace indices have already been removed, remove open indices by filtering on positive values
        contractionindices[k] = filter(in(indexorder), indices)
    end
    partialtrees = collect(Any, 1:length(network))
    _indexordertree!(partialtrees, contractionindices, indexorder)
end

function _indexordertree!(partialtrees, contractionindices, indexorder)
    while length(indexorder) > 0
        c = indexorder[1]
        i1 = findfirst(x->in(c,x), contractionindices)
        if i1 === nothing
            indexorder = indexorder[2:end]
        else
            i2 = findnext(x->in(c,x), contractionindices, i1+1)
            @assert i2 !== nothing
            cs = intersect(contractionindices[i1], contractionindices[i2])
            indexorder = setdiff(indexorder, cs)
            newindices = unique2(vcat(contractionindices[i1], contractionindices[i2]))
            newtree = Any[partialtrees[i1], partialtrees[i2]]
            partialtrees[i1] = newtree
            deleteat!(partialtrees, i2)
            contractionindices[i1] = newindices
            deleteat!(contractionindices, i2)
        end
    end
    if length(partialtrees) > 1
        @assert all(isempty, contractionindices) # disconnected network
        while length(partialtrees) > 1
            partialtrees[end-1] = Any[partialtrees[end-1], partialtrees[end]]
            pop!(partialtrees)
            pop!(contractionindices)
        end
    end
    return partialtrees[1]
end


text:
__precompile__(true)

module MathematicalPredicates

export Predicate, Atom, Negation, Conjunction, Disjunction,
       evaluate

include("Predicate.jl")
include("Atom.jl")
include("Negation.jl")
include("Conjunction.jl")
include("Disjunction.jl")

# external functionality using 'Requires'
using Requires
include("init.jl")

end  # module


text:
function log_ticks(loglims)
    a = round(loglims[1], RoundDown)
    b = round(loglims[2], RoundUp)
    v = log10.(1:9)
    ticks = vcat([i .+ a .+ v for i in 0:b-a]..., [b])
    labels = [tick == round(tick) ? to_latex("10^{$(round(Int, tick))}") : " " for tick in ticks]
    return ticks, labels
end

function find_limits(vals; start=[first(vals), first(vals)], margins=0.05)
    vmin = start[1]
    vmax = start[2]
    for v in vals
        if v < vmin
            vmin = v
        elseif v > vmax
            vmax = v
        end
    end
    diff = vmax - vmin
    vmax += margins * diff
    vmin -= margins * diff
    return vmin, vmax
end

function print_tf(roots, gain)
    # Convert the numerator and denominator to strings
    num = []
    den = []

    for root in roots
        s = printroot(root.pos)
        if root.pole
            push!(den, s)
        else
            push!(num, s)
        end
    end

    if length(num) == 0
        numstr = "1.0"
    elseif length(num) == 1
        numstr = num[1]
    else
        numstr = "(" * join(num, ")(") * ")"
    end
    if length(den) == 0
        denstr = "1.0"
    elseif length(den) == 1
        denstr = den[1]
    else
        denstr = "(" * join(den, ")(") * ")"
    end

    gainstr = "$(round(gain, sigdigits=3)) "
    gainpre = "$(repeat(" ", length(gainstr)))"

    # Figure out the length of the separating line
    len_num = length(numstr)
    len_den = length(denstr)
    dashcount = max(len_num, len_den)

    # Center the numerator or denominator
    if len_num < dashcount
        numstr = "$(repeat(" ", div(dashcount - len_num, 2)))$numstr"
    else
        denstr = "$(repeat(" ", div(dashcount - len_den, 2)))$denstr"
    end
    return gainpre * to_latex(numstr) * "\n " * gainstr * repeat("-", dashcount) * "\n " * gainpre * to_latex(denstr)
end

function printroot(z)
    a = z[1]
    b = z[2]
    if b == 0
        return "s $(a < 0 ? "+" : "-") $(round(abs(a), sigdigits=3))"
    else
        return "s^2 $(a < 0 ? "+" : "-") $(round(2abs(a), sigdigits=3))s + $(round(a^2 + b^2, sigdigits=3))"
    end
end

max_stars_repo_name: albheim/PoleExploration

max_issues_repo_name: albheim/PoleExploration

max_forks_repo_name: albheim/PoleExploration


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

# RUN: julia --startup-file=no %s %t
# RUN: cat %t/* | FileCheck %s

include(joinpath("..", "testhelpers", "llvmpasses.jl"))

struct Foo
    x::Int32
    y::Int32
end

@generated foo(x)=:(ccall("extern foo", llvmcall, $x, ($x,), x))
bar(x) = ntuple(i -> VecElement{Float16}(x[i]), 2)

# CHECK: call half @foo(half zeroext %{{[0-9]+}})
emit(foo, Float16)

# CHECK: call [2 x half] @foo([2 x half] %{{[0-9]+}})
emit(foo, NTuple{2, Float16})

# CHECK: call <2 x half> @foo(<2 x half> %{{[0-9]+}})
emit(foo, NTuple{2, VecElement{Float16}})

# CHECK: call i8 addrspace(3)* @foo(i8 addrspace(3)* %{{[0-9]+}})
emit(foo, Core.LLVMPtr{Float32, 3})

# CHECK: call { i32, i32 } @foo({ i32, i32 } %{{[0-9]+}})
emit(foo, Foo)

# CHECK: define <2 x i16> @julia_bar_{{[0-9]+}}([2 x i16]
emit(bar, NTuple{2, Float16})

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

module SimilaritySearch
abstract type Index end
export Index

include("distances/bits.jl")
include("distances/sets.jl")
include("distances/strings.jl")
include("distances/vectors.jl")
include("distances/cos.jl")
include("res/knn.jl")
include("nns/performance.jl")
include("indexes/pivotselection.jl")
include("indexes/seq.jl")
include("indexes/laesa.jl")
include("indexes/pivotselectiontables.jl")
include("knr/SimilarReferences.jl")
include("graph/graph.jl")
include("utils/aknn.jl")
end


text:
#using SymmetrySnap
using Test
include("../src/SymmetrySnap.jl")

@testset "SymmetrySnap.jl" begin
    # Make sure that we get the same number of lattice points no matter what
    @test length(SymmetrySnap.getLatticePoints([1.234324,3.654,6.213], [5.12436,4.712346,9.2778], [1.23436346,2.678323,9.7])) == length(SymmetrySnap.getLatticePoints([2.435,6.54,6.213], [5.12436,7.12346,2.778], [2.3436346,6.78323,97]))

    # Test cubic symmetry, we should get 48 symmetry operations
    a = [1,2,3]
    b = sqrt(14/3)*[-1,-1,1]
    c = [5,-4,1]/sqrt(3)
    newLattice = Array{Float64, 2}(undef, 3, 3)
    newLattice[1:end, 1] = a
    newLattice[1:end, 2] = b
    newLattice[1:end, 3] = c
    inverse = inv(newLattice)
    @test length(SymmetrySnap.getSymmetryOps(a,b,c,inverse)) == 48

    # Test orthorhombic symmetry, we should get 8 symmetry operations
    a = [1,2,3]
    b = [-1,-1,1]
    c = [5,-4,1]
    newLattice = Array{Float64, 2}(undef, 3, 3)
    newLattice[1:end, 1] = a
    newLattice[1:end, 2] = b
    newLattice[1:end, 3] = c
    inverse = inv(newLattice)
    @test length(SymmetrySnap.getSymmetryOps(a,b,c,inverse)) == 8

    # Test hexagonal symmetry, we should get 24 symmetry operations
    a = [1,0,0]
    b = [-1/2,sqrt(3)/2,0]
    c = [0,0,1.212624]
    newLattice = Array{Float64, 2}(undef, 3, 3)
    newLattice[1:end, 1] = a
    newLattice[1:end, 2] = b
    newLattice[1:end, 3] = c
    inverse = inv(newLattice)
    @test length(SymmetrySnap.getSymmetryOps(a,b,c,inverse)) == 24

    # Make sure we can map atoms to a cell, we should get [0.478362, 0.76541, 0.0030137]
    a = [1,0,0]
    b = [0,1,0]
    c = [0,0,1]
    newLattice = Array{Float64, 2}(undef, 3, 3)
    newLattice[1:end, 1] = a
    newLattice[1:end, 2] = b
    newLattice[1:end, 3] = c
    inverse = inv(newLattice)
    atomPosition = [113.478362, -2093.234590, 1.0030137]
    @test isapprox(SymmetrySnap.mapAtomsToCell(a, b, c, inverse, [atomPosition]), [[0.478362, 0.76541, 0.0030137]])
    """
    a = [1,0,0]
    b = [0,1,0]
    c = [0,0,1]
    newLattice = Array{Float64, 2}(undef, 3, 3)
    newLattice[1:end, 1] = a
    newLattice[1:end, 2] = b
    newLattice[1:end, 3] = c
    inverse = inv(newLattice)
    atom1 = [0.25,0.25,0.25]
    atom2 = [0.75,0.25,0.25]
    operator = [-1 0 0; 0 1 0; 0 0 1]
    @test isapprox(SymmetrySnap.testAtom(a, b, c, inverse, operator, atom1, 1, atom2, 2, 2), [-1 0 0 -1 0 0 -1; 0 1 0 0 -1 0 0; 0 0 1 0 0 -1 0], [0,0,0])
    """
    a = [1,0,0]
    b = [0,1,0]
    c = [0,0,1]
    atomPos = [[0.25001,0.25003,0.24998], [0.2499992,0.25,0.751], [0.2519203,0.749210,0.25777], [0.751359,0.252468,0.248642], [0.75,0.74867530999,0.2525600], [0.753141592653,0.2510293,0.74987654321], [0.2500012387593,0.75130958,0.75], [0.74989132,0.750019284,0.751837592985409287374987810293892834750917823984723894574012983975]]
    newPos = vcat(atomPos...)
    solutionSpace = [0.25,0.25,0.25, 0.25,0.25,0.75, 0.25,0.75,0.25, 0.75,0.25,0.25, 0.75,0.75,0.25, 0.75,0.25,0.75, 0.25,0.75,0.75, 0.75,0.75,0.75]
    answer = SymmetrySnap.proj(solutionSpace, newPos)
    display(answer)
    display(SymmetrySnap.symmetrize(a, b, c, atomPos, [8]))
end


text:
# Modelo de separador de corrente de 4 saidas
type Splitter4
	Splitter4()=begin
		new(
			fill(Fracao(Dict{Symbol,Any}(
				:Brief=>"Fracao de massa"
			)),(3)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in}",
				:PosX=>1,
				:PosY=>0.5
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out1}",
				:PosX=>0,
				:PosY=>0.19
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out2}",
				:PosX=>0,
				:PosY=>0.4
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out3}",
				:PosX=>0,
				:PosY=>0.62
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out4}",
				:PosX=>0,
				:PosY=>0.85
			)),
			[
				:(Fout1.P = Fin.P),
				:(Fout1.T = Fin.T),
				:(Fout1.S = Fin.S),
				:(Fout1.H = Fin.H),
				:(Fout2.P = Fin.P),
				:(Fout2.T = Fin.T),
				:(Fout2.S = Fin.S),
				:(Fout2.H = Fin.H),
				:(Fout3.P = Fin.P),
				:(Fout3.T = Fin.T),
				:(Fout3.S = Fin.S),
				:(Fout3.H = Fin.H),
				:(Fout4.P = Fin.P),
				:(Fout4.T = Fin.T),
				:(Fout4.S = Fin.S),
				:(Fout4.H = Fin.H),
				:(Fout1.F = Fin.F * y(1)),
				:(Fout2.F = Fin.F * y(2)),
				:(Fout3.F = Fin.F * y(3)),
				:(Fout4.F = Fin.F - Fout1.F - Fout2.F - Fout3.F),
			],
			[
				"","","","","","","","","","","","","","","","","","","","",
			],
			[:y,:Fin,:Fout1,:Fout2,:Fout3,:Fout4,]
		)
	end
	y::Array{Fracao}
	Fin::Corrente
	Fout1::Corrente
	Fout2::Corrente
	Fout3::Corrente
	Fout4::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Splitter4
function setEquationFlow(in::Splitter4)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	#Fout4.F = Fin.F * y(4);
	addEquation(20)
end
function atributes(in::Splitter4,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/splitter4"
	drive!(fields,_)
	return fields
end
Splitter4(_::Dict{Symbol,Any})=begin
	newModel=Splitter4()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
#!/usr/bin/env julia
# del2z <delta.z@aliyun.com>

module ClassifyController

using Genie.Renderer
using Genie.Router

mutable struct Book
    title::String
    author::String

    Book(title::String, author::String) = new(title, author)
    Book(; title::String = "", author::String = "") = new(title, author)
end

const BillGatesBooks = Book[
    Book("The Best We Could Do", "Thi Bui"),
    Book("Evicted: Poverty and Profit in the American City", "Matthew Desmond"),
    Book("Believe Me: A Memoir of Love, Death, and Jazz Chickens", "Eddie Izzard"),
    Book("The Sympathizer", "Viet Thanh Nguyen"),
    Book("Energy and Civilization, A History", "Vaclav Smil")
]

function billgatesbooks()
    html!(:classify, :classify, books = BillGatesBooks)
end

function new()
    html!(:classify, :anno, books = BillGatesBooks)
end

function create()
    @params(:label) |> println
    redirect_to(:get_bgbooks)
end

end


text:
import Base: permutedims

function gemm_sp_dense_small(alpha, M::SparseMatrixCSC, B::Matrix, result::Matrix)
    if isone(alpha)
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                row = M.rowval[i]
                val = M.nzval[i]
                @inbounds for j=1:size(B, 2)
                    result[row, j] += val*B[colindex, j]
                end
            end
        end
    else
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                row = M.rowval[i]
                val = alpha*M.nzval[i]
                @inbounds for j=1:size(B, 2)
                    result[row, j] += val*B[colindex, j]
                end
            end
        end
    end
end

function gemm_sp_dense_big(alpha, M::SparseMatrixCSC, B::Matrix, result::Matrix)
    if isone(alpha)
        @inbounds for j=1:size(B, 2)
            @inbounds for colindex = 1:M.n
                m2 = B[colindex, j]
                @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                    row = M.rowval[i]
                    result[row, j] += M.nzval[i]*m2
                end
            end
        end
    else
        @inbounds for j=1:size(B, 2)
            @inbounds for colindex = 1:M.n
                m2 = alpha*B[colindex, j]
                @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                    row = M.rowval[i]
                    result[row, j] += M.nzval[i]*m2
                end
            end
        end
    end
end

function gemm_dense_adj_sp(alpha, B::Matrix, M::SparseMatrixCSC, result::Matrix)
    if isone(alpha)
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                row = M.rowval[i]
                val = conj(M.nzval[i])
                @inbounds for j=1:size(B, 1)
                    result[j, row] += val*B[j, colindex]
                end
            end
        end
    else
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                row = M.rowval[i]
                val = alpha*conj(M.nzval[i])
                @inbounds for j=1:size(B, 1)
                    result[j, row] += val*B[j, colindex]
                end
            end
        end
    end
end

function gemm_adj_sp_dense_small(alpha, M::SparseMatrixCSC, B::Matrix, result::Matrix)
    dimB = size(result,2)
    if isone(alpha)
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                mi = conj(M.nzval[i])
                mrowvali = M.rowval[i]
                @inbounds for j=1:dimB
                    result[colindex, j] += mi*B[mrowvali, j]
                end
            end
        end
    else
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                mi = conj(M.nzval[i])*alpha
                mrowvali = M.rowval[i]
                @inbounds for j=1:dimB
                    result[colindex, j] += mi*B[mrowvali, j]
                end
            end
        end
    end
end

function gemm!(alpha, M::SparseMatrixCSC, B::Matrix, beta, result::Matrix)
    if iszero(beta)
        fill!(result, beta)
    elseif !isone(beta)
        rmul!(result, beta)
    end
    if nnz(M) > 550
        gemm_sp_dense_big(alpha, M, B, result)
    else
        gemm_sp_dense_small(alpha, M, B, result)
    end
end

function gemm!(alpha, B::Matrix, M::SparseMatrixCSC, beta, result::Matrix)
    if iszero(beta)
        fill!(result, beta)
    elseif !isone(beta)
        rmul!(result, beta)
    end
    dimB = size(result,1)
    if isone(alpha)
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                mi = M.nzval[i]
                mrowvali = M.rowval[i]
                @inbounds for j=1:dimB
                    result[j, colindex] += mi*B[j, mrowvali]
                end
            end
        end
    else
        @inbounds for colindex = 1:M.n
            @inbounds for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                mi = M.nzval[i]*alpha
                mrowvali = M.rowval[i]
                @inbounds for j=1:dimB
                    result[j, colindex] += mi*B[j, mrowvali]
                end
            end
        end
    end
end

function gemm!(alpha, M_::Adjoint{T,<:SparseMatrixCSC{T}}, B::Matrix, beta, result::Matrix) where T
    M = M_.parent
    if nnz(M) > 550
        LinearAlgebra.mul!(result, M_, B, alpha, beta)
    else
        if iszero(beta)
            fill!(result, beta)
        elseif !isone(beta)
            rmul!(result, beta)
        end
        gemm_adj_sp_dense_small(alpha, M, B, result)
    end
end

function gemm!(alpha, B::Matrix, M::Adjoint{T,<:SparseMatrixCSC{T}}, beta, result::Matrix) where T
    if iszero(beta)
        fill!(result, beta)
    elseif !isone(beta)
        rmul!(result, beta)
    end
    gemm_dense_adj_sp(alpha, B, M.parent, result)
end

function gemv!(alpha, M::SparseMatrixCSC, v::Vector, beta, result::Vector)
    if iszero(beta)
        fill!(result, beta)
    elseif !isone(beta)
        rmul!(result, beta)
    end
    if isone(alpha)
        @inbounds for colindex = 1:M.n
            vj = v[colindex]
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[M.rowval[i]] += M.nzval[i]*vj
            end
        end
    else
        @inbounds for colindex = 1:M.n
            vj = alpha*v[colindex]
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[M.rowval[i]] += M.nzval[i]*vj
            end
        end
    end
end

function gemv!(alpha, v::Vector, M::SparseMatrixCSC, beta, result::Vector)
    if iszero(beta)
        fill!(result, beta)
    elseif !isone(beta)
        rmul!(result, beta)
    end
    if isone(alpha)
        @inbounds for colindex=1:M.n
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[colindex] += M.nzval[i]*v[M.rowval[i]]
            end
        end
    else
        @inbounds for colindex=1:M.n
            for i=M.colptr[colindex]:M.colptr[colindex+1]-1
                result[colindex] += M.nzval[i]*alpha*v[M.rowval[i]]
            end
        end
    end
end

function sub2sub(shape1::NTuple{N, Int}, shape2::NTuple{M, Int}, I::CartesianIndex{N}) where {N, M}
    linearindex = LinearIndices(shape1)[I.I...]
    CartesianIndices(shape2)[linearindex]
end

function ptrace(x, shape_nd::Vector{Int}, indices::Vector{Int})
    shape_nd = (shape_nd...,)
    N = div(length(shape_nd), 2)
    shape_2d = (x.m, x.n)
    shape_nd_after = ([i ∈ indices || i-N ∈ indices ? 1 : shape_nd[i] for i=1:2*N]...,)
    shape_2d_after = (prod(shape_nd_after[1:N]), prod(shape_nd_after[N+1:end]))
    I_nd_after_max = CartesianIndex(shape_nd_after...)
    y = spzeros(eltype(x), shape_2d_after...)
    for I in eachindex(x)
        I_nd = sub2sub(shape_2d, shape_nd, I)
        if I_nd.I[indices] != I_nd.I[indices .+ N]
            continue
        end
        I_after = sub2sub(shape_nd_after, shape_2d_after, min(I_nd, I_nd_after_max))
        y[I_after] += x[I]
    end
    y
end

function permutedims(x, shape, perm)
    shape = (shape...,)
    shape_perm = ([shape[i] for i in perm]...,)
    y = spzeros(eltype(x), x.m, x.n)
    for I in eachindex(x)
        linear_index = LinearIndices((x.m, x.n))[I.I...]
        cartesian_index = CartesianIndices(shape)[linear_index]
        cartesian_index_perm = [cartesian_index[i] for i=perm]
        linear_index_perm = LinearIndices(shape_perm)[cartesian_index_perm...]
        J = Tuple(CartesianIndices((x.m, x.n))[linear_index_perm])
        y[J...] = x[I.I...]
    end
    y
end


text:
using PyCall

@testset "Space" begin

@testset "BoxSpace" begin
    @test occursin(0.5, BoxSpace(0, 1)) == true
    @test occursin(0.0, BoxSpace(0, 1)) == true
    @test occursin(1.0, BoxSpace(0, 1)) == true
    @test occursin(-1.0, BoxSpace(0, 1)) == false
    @test occursin(-Inf, BoxSpace(0, 1)) == false
    @test occursin([0.5], BoxSpace(0, 1)) == true

    @test occursin([0, 0], BoxSpace([-1, -2], [1, 2])) == true
    @test occursin([0, 3], BoxSpace([-1, -2], [1, 2])) == false
    @test occursin([0, 0], BoxSpace([-1, -2], [1, 2])) == true
end

@testset "DiscreteSpace" begin
    @test occursin(0, DiscreteSpace(10)) == true
    @test occursin(5, DiscreteSpace(10)) == true
    @test occursin(10, DiscreteSpace(10)) == false
end

@testset "MultiBinarySpace" begin
    @test occursin([true false; true false], MultiBinarySpace(2,2)) == true
    @test occursin([true false], MultiBinarySpace(2,2)) == false
end

@testset "MultiDiscreteSpace" begin
    @test occursin(0, MultiDiscreteSpace([2,3,2])) == true
    @test occursin(1, MultiDiscreteSpace([2,3,2])) == true
    @test occursin(2, MultiDiscreteSpace([2,3,2])) == false
    @test occursin([1,1,1], MultiDiscreteSpace([2,3,2])) == true
    @test occursin([0,0,0], MultiDiscreteSpace([2,3,2])) == true
    @test occursin([3,3,3], MultiDiscreteSpace([2,3,2])) == false
end

@testset "Space Tuple" begin
    @test occursin(([0.5], 5, [true true; true true], [1, 1]),
                   (BoxSpace(0,1), DiscreteSpace(5), MultiBinarySpace(2,2), MultiDiscreteSpace([2,2]))) == false
    @test occursin(([0.5], 0, [true true; true true], [1, 1]),
                   (BoxSpace(0,1), DiscreteSpace(5), MultiBinarySpace(2,2), MultiDiscreteSpace([2,2]))) == true
    @test contains((),
                   (BoxSpace(0,1), DiscreteSpace(5), MultiBinarySpace(2,2), MultiDiscreteSpace([2,2]))) == false
end

@testset "Space Dict" begin
    @test occursin(
        Dict(
            "sensors" => Dict(
                "position" => [-10, 0, 10],
                "velocity" => [0.1, 0.2, 0.3],
                "front_cam" => (rand(10, 10, 3), rand(10, 10, 3)),
                "rear_cam" => rand(10,10,3)),
            "ext_controller" => [2, 1, 1],
            "inner_state" => Dict(
                "charge" => 35,
                "system_checks" => rand(Bool, 10),
                "job_status" => Dict(
                    "task" => 3,
                    "progress" => 23))),
        Dict(
            "sensors" =>  Dict(
                "position"=> BoxSpace(-100, 100, (3,)),
                "velocity"=> BoxSpace(-1, 1, (3,)),
                "front_cam"=> (BoxSpace(0, 1, (10, 10, 3)),
                               BoxSpace(0, 1, (10, 10, 3))),
                "rear_cam" => BoxSpace(0, 1, (10, 10, 3))),
            "ext_controller" => MultiDiscreteSpace([5, 2, 2]),
            "inner_state" => Dict(
                "charge" => DiscreteSpace(100),
                "system_checks" => MultiBinarySpace(10),
                "job_status" => Dict(
                    "task" => DiscreteSpace(5),
                    "progress" => BoxSpace(0, 100))))) == true
end

@testset "gym space transformation" begin
    @pyimport gym
    env = gym.make("CartPole-v0")
    @test gymspace2jlspace(env[:action_space]) == DiscreteSpace(2)
    @test gymspace2jlspace(env[:observation_space]) == BoxSpace(Array{Float32}([-4.8, -3.4028235e38, -0.41887903, -3.4028235e38]),
                                                                 Array{Float32}([4.8, 3.4028235e38, 0.41887903, 3.4028235e38]))
    env = gym.make("Blackjack-v0")
    @test gymspace2jlspace(env[:observation_space]) == (DiscreteSpace(32), DiscreteSpace(11), DiscreteSpace(2))

end
end


text:
using DataStructures

function main()
    N, M = parse.(Int, split(readline()))
    G = [Vector{Int}() for _ = 1:N]
    for _ = 1:M
        u, v = parse.(Int, split(readline()))
        push!(G[u], v)
    end
    S, T = parse.(Int, split(readline()))

    # Solve
    dists = -ones(Int, N, 3)
    dists[S, 1] = 0
    que = Queue{Tuple{Int, Int}}()
    enqueue!(que, (S, 1))
    while !isempty(que)
        v, s = dequeue!(que)
        for v′ ∈ G[v]
            s′ = s % 3 + 1
            if dists[v′, s′] == -1
                dists[v′, s′] = dists[v, s] + 1
                enqueue!(que, (v′, s′))
            end
        end
    end

    if dists[T, 1] == -1
        println("-1")
    else
        println("$(dists[T, 1] ÷ 3)")
    end
end

main()


text:
export trtest, trcauchy, trdogleg

# test if step is acceptable and update size of trust region
function trtest(f::Function,fc,df,H,xc,st,Delta::Real;
    rhom=1e-4, rhoL=.25, rhoH=0.75, wdown=0.5, wup=1.5,C=1e2)

    ared = fc - f(xc+st)
    pred = - dot(df,st) - dot(st,H*st)/2

    if ared/pred < rhom # reject and decrease Delta
        return xc,wdown*Delta,-2,ared
    elseif (rhom <= ared/pred) && (ared/pred < rhoL)
        return xc+st, wdown*Delta,-1,ared
    elseif (rhoL <= ared/pred) && (ared/pred < rhoH)
        return xc+st, Delta,0,ared
    elseif (rhoH < ared/pred) && (norm(st) <Delta)
        return xc+st, Delta,0,ared
    elseif (rhoH < ared/pred)
        return xc+st, min(wup*Delta,C*norm(df)),1,ared
    end
end

# implementation of cauchy method
function trcauchy(f::Function,J::Function,H::Function,x::Vector;
	Delta=2.0,maxIter=20,atol=1e-8,out::Int=0,storeInterm::Bool=false,
	trt::Function=(f,fc,df,d2f,x,st,Delta)->trtest(f,fc,df,d2f,x,st,Delta))

    his = zeros(maxIter,4)
    X   = (storeInterm) ? zeros(length(x),maxIter) : []
    i   = 1; flag = -1; alph=0.0
    while i<=maxIter
        fc = f(x)
        df = J(x)
        his[i,1:2] = [norm(fc) norm(df)]
        if storeInterm; X[:,i] = x; end;

        if(norm(df)<atol)
            flag = 0
            his = his[1:i,:]
            break
        end

        # get step size
        d2f  = H(x)
        curv = dot(df,d2f*df)
        if curv <= 0
            alph = Delta/norm(df)
        else
            alph = min(Delta/norm(df), norm(df)^2/curv)
        end

        # check for acceptance
        xn,Deltan,trflag,ared = trt(f,fc,df,d2f,x,-alph*df,Delta)
        his[i,3:4] = [trflag,Delta]
        if out>0
            @printf "iter=%04d\t|f|=%1.4e\t |df|=%1.2e\tTR=%d\tDelta=%1.2f\n" i his[i,1] his[i,2] his[i,3] his[i,4]
        end
        # update
        x     = copy(xn)
        Delta = Deltan
        i+=1
    end
    i = min(maxIter,i)

    if out>=0
        if flag==-1
            println(@sprintf("trcauchy iterated maxIter (=%d) times but reached only atol of %1.2e instead of tol=%1.2e",i,his[i,2],atol))
        elseif out>1
            println(@sprintf("trcauchy achieved desired atol of %1.2e at iteration %d.",atol,i))
        end
    end
    if storeInterm; X = X[:,1:i]; end
    return x,flag,his,X
end

# dogleg method
function trdogleg(f::Function,J::Function,H::Function,x::Vector;
	Delta=2.0,maxIter=20,atol=1e-8,out::Int=0,storeInterm::Bool=false,
	trt::Function=(f,fc,df,d2f,x,st,Delta)->trtest(f,fc,df,d2f,x,st,Delta))

    his = zeros(maxIter,5)
    X   = (storeInterm) ? zeros(length(x),maxIter) : []
    i   = 1; flag = -1; alph=0.0
    while i<=maxIter
        fc = f(x)
        df = J(x)
        his[i,1:2] = [norm(fc) norm(df)]
        if storeInterm; X[:,i] = x; end;

        if(norm(df)<atol)
            flag = 0
            his = his[1:i,:]
            break
        end

        # get cauchy point
        d2f  = H(x)
        curv = dot(df,d2f*df)
        if curv <= 0
            alph = Delta/norm(df)
        else
            alph = min(Delta/norm(df), norm(df)^2/curv)
        end
        st = -alph*df
        xcp = x + st
        flagnt = 0

        if norm(xcp-x) < Delta
            # get Newton point
            try
               sn = - d2f \ df
               xn = x + sn
               # check angle condition
               angl = dot(xn-xcp,xcp-x)
               if (angl>0) && (norm(xcp-x) < norm(sn)) && norm(sn) <= Delta
                  flagnt=1
                  st = sn # take newton point
               elseif (angl>0) # build dogleg
                  flagnt = 2
                  step = Delta/10; cnt = 1;

                  while (norm(xcp + cnt*step*(xn-xcp) - x) < Delta) && (cnt < 100)
                       cnt+=1
                  end
                  xt = xcp + (cnt-1)*step*(xn-xcp)
                  st = xt-x
               end
            catch err
                if ~isa(err,Base.LinAlg.SingularException)
                    throw(err)
                end
            end
        end

        # check for acceptance
        xn,Deltan,trflag,ared = trt(f,fc,df,d2f,x,st,Delta)
        his[i,3:5] = [trflag,flagnt,Delta]
        if out>0
            @printf "iter=%04d\t|f|=%1.2e\t|df|=%1.2e\tTR=%d\tnt?=%d\tDelta=%1.2f\n" i his[i,1] his[i,2] his[i,3] his[i,4] his[i,5]
        end
        # update
        x     = xn
        Delta = Deltan
        i+=1
    end
    i = min(maxIter,i)

    if out>=0
        if flag==-1
            println(@sprintf("trdogleg iterated maxIter (=%d) times but reached only atol of %1.2e instead of tol=%1.2e",i,his[i,2],atol))
        elseif out>1
            println(@sprintf("trdogleg achieved desired atol of %1.2e at iteration %d.",atol,i))
        end
    end
    if storeInterm; X = X[:,1:i]; end
    return x,flag,his,X
end


text:
#####################
# Save lines
#####################

function save_lines(name::String, lines::CircuitFormatLines)
    open(name, "w") do f
        for line in lines
            println(f, line)
        end
    end
end

#####################
# decompile for nodes
#####################

"Decompile for sdd circuit, used during saving of circuits to file"
decompile(n::StructLiteralNode, node2id, vtree2id)::UnweightedLiteralLine =
    UnweightedLiteralLine(node2id[n], vtree2id[n.vtree], literal(n), false)

decompile(n::StructConstantNode, node2id, vtree2id)::AnonymousConstantLine =
    AnonymousConstantLine(node2id[n], constant(n), false)

decompile(n::Struct⋁Node, node2id, vtree2id)::DecisionLine{SDDElement} =
    DecisionLine(node2id[n], vtree2id[n.vtree], UInt32(num_children(n)), map(c -> make_element(c, node2id), children(n)))

make_element(n::Struct⋀Node, node2id) =
    SDDElement(node2id[n.children[1]],  node2id[n.children[2]])

make_element(n::StructLogicalΔNode, node2id) =
    error("Given circuit is not an SDD, its decision node elements are not conjunctions.")

# TODO: decompile for logical circuit to some file format

#####################
# build maping
#####################

function get_node2id(ln::AbstractVector{X}, T::Type)where X #<: T#::Dict{T, ID}
    node2id = Dict{T, ID}()
    outnodes = filter(n -> !(GateType(n) isa ⋀Gate), ln)
    sizehint!(node2id, length(outnodes))
    index = ID(0) # node id start from 0
    for n in outnodes
        node2id[n] = index
        index += ID(1)
    end
    node2id
end

function get_vtree2id(ln::PlainVtree):: Dict{PlainVtreeNode, ID}
    vtree2id = Dict{PlainVtreeNode, ID}()
    sizehint!(vtree2id, length(ln))
    index = ID(0) # vtree id start from 0

    for n in ln
        vtree2id[n] = index
        index += ID(1)
    end
    vtree2id
end

#####################
# saver for circuits
#####################

"Returns header for SDD file format"
function sdd_header()
    """
    c ids of sdd nodes start at 0
    c sdd nodes appear bottom-up, children before parents
    c
    c file syntax:
    c sdd count-of-sdd-nodes
    c F id-of-false-sdd-node
    c T id-of-true-sdd-node
    c L id-of-literal-sdd-node id-of-vtree literal
    c D id-of-decomposition-sdd-node id-of-vtree number-of-elements {id-of-prime id-of-sub}*
    c
    c File generated by Juice.jl
    c"""
end

function save_sdd_file(name::String, circuit::DecoratorΔ, vtree::PlainVtree)
    save_sdd_file(name, origin(circuit, StructLogicalΔNode), vtree)
end

"Save a SDD circuit to file"
function save_sdd_file(name::String, circuit::StructLogicalΔ, vtree::PlainVtree)
    #TODO no need to pass the vtree, we can infer it from origin?
    @assert endswith(name, ".sdd")
    node2id = get_node2id(circuit, StructLogicalΔNode)
    vtree2id = get_vtree2id(vtree)
    formatlines = Vector{CircuitFormatLine}()
    append!(formatlines, parse_sdd_file(IOBuffer(sdd_header())))
    push!(formatlines, SddHeaderLine(num_nodes(circuit)))
    for n in filter(n -> !(GateType(n) isa ⋀Gate), circuit)
        push!(formatlines, decompile(n, node2id, vtree2id))
    end
    save_lines(name, formatlines)
end

"Save a circuit to file"
save_circuit(name::String, circuit::StructLogicalΔ, vtree::PlainVtree) = save_sdd_file(name, circuit, vtree)


text: Dict{String, Any}("black" => Dict{String, Any}("allow_prereleases" => Dict{String, Any}("value" => "true", "type" => "bool"), "python" => Dict{String, Any}("value" => ">3.6", "type" => "string"), "version" => Dict{String, Any}("value" => ">=18.9b0", "type" => "string")))


text:
# This file was generated, do not modify it. # hide
rng = StableRNG(5123) # for reproducibility
m = machine(forest, X, y)
r = range(forest, :n, lower=10, upper=1000)
curves = learning_curve!(m, resampling=Holdout(fraction_train=0.8, rng=rng),
                         range=r, measure=rms);


text:
using Test
using LoopVectorization

using LinearAlgebra
import InteractiveUtils

InteractiveUtils.versioninfo(stdout; verbose = true)

const START_TIME = time()

function clenshaw(x, coeff)
    len_c = length(coeff)
    tmp = zero(x)
    ret = zero(x)
    for i in len_c:-1:2
        ret     = muladd(x,2tmp,coeff[i]-ret)
        ret,tmp = tmp,ret
    end
    ret = muladd(x,tmp,coeff[1]-ret)
    return ret
end

"""
Causes `check_args` to fail.
"""
struct FallbackArrayWrapper{T,N} <: AbstractArray{T,N}
    data::Array{T,N}
end
Base.size(A::FallbackArrayWrapper) = size(A.data)
Base.@propagate_inbounds Base.getindex(A::FallbackArrayWrapper, i::Vararg{Int, N}) where {N} = getindex(A.data, i...)
Base.@propagate_inbounds Base.setindex!(A::FallbackArrayWrapper, v, i::Vararg{Int, N}) where {N} = setindex!(A.data, v, i...)
Base.IndexStyle(::Type{<:FallbackArrayWrapper}) = IndexLinear()

@show LoopVectorization.REGISTER_COUNT

const RUN_SLOW_TESTS = LoopVectorization.REGISTER_COUNT ≤ 16 || !parse(Bool, get(ENV, "GITHUB_ACTIONS", "false"))
@show RUN_SLOW_TESTS

@time @testset "LoopVectorization.jl" begin

    @test isempty(detect_unbound_args(LoopVectorization))

    @time include("printmethods.jl")

    @time include("can_avx.jl")

    @time include("fallback.jl")

    @time include("utils.jl")

    @time include("arraywrappers.jl")

    @time include("check_empty.jl")

    @time include("offsetarrays.jl")

    @time include("tensors.jl")

    @time include("map.jl")

    @time include("filter.jl")

    @time include("mapreduce.jl")

    @time include("ifelsemasks.jl")

    @time include("dot.jl")

    @time include("special.jl")

    @time include("gemv.jl")

    @time include("miscellaneous.jl")

    @time include("copy.jl")

    @time include("broadcast.jl")

    @time include("gemm.jl")
end

const ELAPSED_MINUTES = (time() - START_TIME)/60
@test ELAPSED_MINUTES < 180


text:
module ParticleFilters

using POMDPs
import POMDPs: pdf, mode, update, initialize_belief, support
import POMDPs: statetype, isterminal, observation
import POMDPs: generate_s
import POMDPs: action, value
import POMDPs: implemented
import POMDPs: sampletype
import POMDPModelTools: obs_weight

using StatsBase
using Random
using Statistics
using LinearAlgebra
using POMDPPolicies
using POMDPModelTools # for weighted_iterator

import Random: rand
import Statistics: mean, var, cov, std
import StatsBase: weights

export
    AbstractParticleBelief,
    ParticleCollection,
    WeightedParticleBelief,
    BasicParticleFilter,
    ImportanceResampler,
    LowVarianceResampler,
    SIRParticleFilter,
    UnweightedParticleFilter,
    ParticleFilterModel,
    PredictModel,
    ReweightModel

export
    resample,
    predict,
    predict!,
    reweight,
    reweight!,
    particle_memory

export
    weighted_posterior_belief,
    predicted_belief

export
    n_particles,
    particles,
    weighted_particles,
    weight_sum,
    weight,
    particle,
    weights,
    obs_weight,
    n_init_samples,
    runfilter

export
    pdf,
    get_probs,
    mode,
    update,
    support,
    initialize_belief,
    mean,
    cov,
    var,
    std

export
    SimpleParticleFilter

include("beliefs.jl")
include("basic.jl")
include("resamplers.jl")
include("sir.jl")
include("unweighted.jl")
include("models.jl")
include("pomdps.jl")
include("policies.jl")
include("runfilter.jl")
include("deprecated.jl")

end # module


text:
using jlfmt
using Test

@testset "jlfmt" begin
    @test jlfmt.main(["--help"]) === nothing
end

max_stars_repo_name: tkf/jlfmt

max_issues_repo_name: tkf/jlfmt

max_forks_repo_name: tkf/jlfmt