text:
"""
    TriMesh(p::PointCloud, res::Int = 32; algo = :MarchingCubes)

Initialize TriMesh from PointCloud structures having specified `resolution`.

`algo` is the mode to be used to convert binary voxels. Available `algo` are
[:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets].

See also: [`PointCloudToTriMesh`](@ref)
"""
function TriMesh(p::PointCloud, res::Int = 32; algo = :MarchingCubes)
    voxel = pointcloud_to_voxel(p, res)
    v = VoxelGrid(voxel)
    return TriMesh(v; algo = algo)
end

"""
    TriMesh(v::VoxelGrid; thresh::Number = 0.5f0, algo = :MarchingCubes)

Initialize TriMesh from the VoxelGrid structures.

`threshold` is the threshold from which to make binary voxels, and `algo`
is the mode to be used to convert binary voxels. Available `algo` are
[:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets].

See also: [`VoxelGridToTriMesh`](@ref)
"""
function TriMesh(v::VoxelGrid; thresh::Number = 0.5f0, algo = :MarchingCubes)
    verts, faces = voxel_to_trimesh(v, Float32(thresh), algo)
    return TriMesh(verts, faces)
end

"""
    PointCloud(m::TriMesh, npoints::Int = 1000)

Initialize PointCloud, having `npoints`, from the TriMesh structure.

See also: [`TriMeshToPointCloud`](@ref)
"""
function PointCloud(m::TriMesh, npoints::Int = 1000)
    p = sample_points(m, npoints)
    return PointCloud(p)
end

"""
    PointCloud(v::VoxelGrid, npoints::Int = 1000; thresh::Number = 0.5f0, algo = :MarchingCubes)

Initialize PointCloud from the VoxelGrid Structures having `npoints`.

`thresh` is the threshold from which to make binary voxels, and `algo`
is the mode to be used to convert binary voxels. Available `algo` are
[:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets].

See also: [`VoxelGridToPointCloud`](@ref)
"""
function PointCloud(
    v::VoxelGrid,
    npoints::Int = 1000;
    thresh::Number = 0.5f0,
    algo = :MarchingCubes,
)
    m = TriMesh(v, thresh = thresh, algo = algo)
    points = sample_points(m, npoints)
    return PointCloud(points)
end

"""
    VoxelGrid(m::TriMesh, res::Int = 32)

Initialize VoxelGrid, having specified `resolution`, from TriMesh Structure.

See also: [`TriMeshToVoxelGrid`](@ref)
"""
function VoxelGrid(m::TriMesh, res::Int = 32)
    v = trimesh_to_voxel(m, res)
    return VoxelGrid(v)
end

"""
    VoxelGrid(p::PointCloud, res::Int = 32)

Initialize VoxelGrid, having specified `resolution`, from PointCloud Structure.

See also: [`PointCloudToVoxelGrid`](@ref)
"""
function VoxelGrid(p::PointCloud, res::Int = 32)
    v = pointcloud_to_voxel(p, res)
    return VoxelGrid(v)
end

function pointcloud_to_voxel(pcloud::PointCloud, resolution::Int = 32)
    p = cpu(pcloud.points)
    _, _N, _B = size(p)
    verts_max = maximum(maximum(p, dims = 1), dims = 2)
    verts_min = minimum(minimum(p, dims = 1), dims = 2)
    cloud = ((p .- verts_min) ./ (verts_max - verts_min))

    ind1 = Int[]
    ind2 = Int[]
    ind3 = Int[]

    for x = 1:resolution, y = 1:resolution, z = 1:resolution
        push!(ind1, x)
        push!(ind2, y)
        push!(ind3, z)
    end
    ind1 = reshape(ind1, 1, :)
    ind2 = reshape(ind2, 1, :)
    ind3 = reshape(ind3, 1, :)

    grid_points = cat(ind1, ind2, ind3, dims = 1)
    grid_points = ((grid_points .+ 0.5) ./ resolution)

    nn_idx = cat(
        [
            CartesianIndex.(
                reduce(vcat, knn(KDTree(cloud[:, :, i]), grid_points, 1)[1]),
                i,
            ) for i = 1:_B
        ]...,
        dims = 2,
    )

    dists_vec = (grid_points .- cloud[:, nn_idx])
    dists = sum((dists_vec .^ 2), dims = 1)
    dists = reshape(dists, resolution, resolution, resolution, _B)
    voxels = typeof(similar(pcloud.points, 1, 1, 1, 1))(
        dists .<= (0.6 / (resolution * resolution)),
    )
    return voxels
end

function trimesh_to_voxel(m::TriMesh{T,R,S}, res::Int = 32) where {T,R,S}
    verts_list = get_verts_list(m)
    faces_list = get_faces_list(m)

    _B = length(verts_list)
    voxels = S{T,4}(undef, res, res, res, _B)
    # print(typeof(verts_list))
    for (i, v, f) in zip(1:_B, verts_list, faces_list)
        voxels[:, :, :, i] = _voxelize(cpu(v), f, res)
    end
    return voxels
end

function _voxelize(v::Array{<:AbstractFloat,2}, f::Array{<:Integer,2}, resolution::Int = 32)
    verts_max = maximum(v)
    verts_min = minimum(v)
    verts = ((v .- verts_min) ./ (verts_max - verts_min))

    points = verts
    smallest_side = (1.0 / resolution)^2

    v1 = verts[:, f[1, :]]
    v2 = verts[:, f[2, :]]
    v3 = verts[:, f[3, :]]

    while true
        side_1 = sum((v1 .- v2) .^ 2, dims = 1)
        side_2 = sum((v2 .- v3) .^ 2, dims = 1)
        side_3 = sum((v3 .- v1) .^ 2, dims = 1)
        sides = cat(side_1, side_2, side_3, dims = 1)
        sides = reshape(maximum(sides, dims = 1), :)

        keep = sides .> smallest_side
        !(any(keep)) && break

        v1 = v1[:, keep]
        v2 = v2[:, keep]
        v3 = v3[:, keep]

        v4 = (v1 + v3) / 2
        v5 = (v1 + v2) / 2
        v6 = (v2 + v3) / 2

        points = cat(points, v4, v5, v6, dims = 2)
        vertex_set = [v1, v2, v3, v4, v5, v6]
        new_traingles = [
            1 4 5
            5 2 6
            5 4 6
            4 3 6
        ]
        new_verts = []
        for i = 1:4
            for j = 1:3
                if i == 1
                    push!(new_verts, vertex_set[new_traingles[i, j]])
                else
                    new_verts[j] =
                        cat(new_verts[j], vertex_set[new_traingles[i, j]], dims = 2)
                end
            end
        end

        v1, v2, v3 = new_verts
    end

    voxels = zeros(resolution, resolution, resolution)
    idx =
        (
            (round.(Int, points .* (resolution - 1), RoundToZero) .+ resolution) .%
            resolution
        ) .+ 1
    voxels[CartesianIndex.(idx[1, :], idx[2, :], idx[3, :])] .= 1
    return voxels
end

function voxel_to_trimesh(v::VoxelGrid, thresh, algo)
    algo in [:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets] ||
        error("given algo: $(algo) is not supported. Accepted algo are
              {:Exact,:MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets}.")

    _assert_voxel(v) ||
        error("invalid VoxelGrid, found element which is not between [0,1].")
    voxel = v.voxels
    T = typeof(similar(voxel, 1, 1))
    R = Array{UInt32,2}
    verts = Vector{T}(undef, 0)
    faces = Vector{R}(undef, 0)
    method = algo == :Exact ? _voxel_exact : _voxel_algo

    for i = 1:size(voxel, 4)
        v, f = method(cpu(voxel[:, :, :, i]), thresh, algo)
        v = T(v)
        v = v ./ maximum(v)
        f = R(f)
        push!(verts, v)
        push!(faces, f)
    end
    return verts, faces
end

function _voxel_algo(v, thresh, algo)
    res = size(v, 1)
    voxel = zeros(eltype(v), res + 2, res + 2, res + 2)
    voxel[2:end-1, 2:end-1, 2:end-1] .= v
    algo = eval(algo)
    voxel = voxel .>= thresh
    v, f = isosurface(voxel, algo(iso = 0, insidepositive = true))
    verts = reduce(hcat, v)
    faces = reduce(hcat, f)
    return verts, faces
end

function _voxel_exact(voxel, thresh, algo)
    voxel = voxel .>= thresh
    res = size(voxel, 1)
    if res >= 3
        voxel[2:res-1, 2:res-1, 2:res-1] =
            voxel[2:res-1, 2:res-1, 2:res-1] .& .!(
                voxel[1:res-2, 2:res-1, 2:res-1] .&
                voxel[3:res, 2:res-1, 2:res-1] .&
                voxel[2:res-1, 1:res-2, 2:res-1] .&
                voxel[2:res-1, 3:res, 2:res-1] .&
                voxel[2:res-1, 2:res-1, 1:res-2] .&
                voxel[2:res-1, 2:res-1, 3:res],
            )
    end

    v = Float32[]
    f = UInt32[]

    for i in CartesianIndices(voxel)
        if voxel[i]
            _add_face!(i, v, f)
        end
    end

    verts = reshape(v, 3, :)
    faces = reshape(f, 3, :)

    return verts, faces
end

function _add_face!(i, v, f)

    x, y, z = i[1], i[2], i[3]

    cube_verts = [
        x - 1,
        y - 1,
        z - 1,
        x - 1,
        y - 1,
        z,
        x - 1,
        y,
        z - 1,
        x - 1,
        y,
        z,
        x,
        y - 1,
        z - 1,
        x,
        y - 1,
        z,
        x,
        y,
        z - 1,
        x,
        y,
        z,
    ]

    cube_faces = [
        1,
        7,
        5,
        1,
        3,
        7,
        1,
        4,
        3,
        1,
        2,
        4,
        3,
        8,
        7,
        3,
        4,
        8,
        5,
        7,
        8,
        5,
        8,
        6,
        1,
        5,
        6,
        1,
        6,
        2,
        2,
        6,
        8,
        2,
        8,
        4,
    ]

    offset = length(v) / 3
    append!(v, cube_verts)
    append!(f, cube_faces .+ offset)
end


text:
push!(LOAD_PATH,"../")
using Documenter, QXZoo, Random

makedocs(;
    modules=[QXZoo],
    authors="QuantEx team",
    repo="https://github.com/JuliaQX/QXZoo.jl/blob/{commit}{path}#L{line}",
    sitename="QXZoo.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaQX.github.io/QXZoo.jl",
        assets=String[],
    ),
    pages=Any[
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Manual" => Any[
            "Gates" => Any[ "gates_circuits/GateOps.md",
                            "gates_circuits/DefaultGates.md",
                            "gates_circuits/CompositeGates.md",
                            "gates_circuits/GateMap.md",
                        ],
            "Circuits" => "gates_circuits/circuits.md",
            "Algorithms" => Any[    "NCU" => "algo/ncu.md",
                                    "Grover" => "algo/grover.md",
                                    "QFT" => "algo/qft.md",
                                    "RQC" => "algo/rqc.md"
                                ],
        ],
    ]
)
deploydocs(;
    repo="github.com/JuliaQX/QXZoo.jl",
)

ext: jl


text:
using LinearAlgebra

@testset "Hyper Cubes" begin
    @testset "Core" begin
        x = centered(HyperCube)
        @test origin(x) == Vec3f0(-0.5)
        @test width(x) == 1.0f0
        @test widths(x) == Vec3f0(1.0)
        @test maximum(x) == Vec3f0(0.5)
        @test minimum(x) == Vec3f0(-0.5)
    end
end
@testset "Hyper Rectangles" begin

@testset "constructors and containment" begin
    # HyperRectangle(vals...)
    a = HyperRectangle(0,0,1,1)
    @test a == HyperRectangle{2,Int}(Vec(0,0),Vec(1,1))

    b = HyperRectangle(0,0,1,1,1,2)
    @test b == HyperRectangle{3,Int}(Vec(0,0,1),Vec(1,1,2))

    a = AABB(0,0,0,1,1,1)
    b = AABB(Vec(0,0,0), Vec(1,1,1))
    @test a == HyperRectangle(0,0,0,1,1,1)
    @test a == b

    a = HyperRectangle{4, Float64}()
    @test a == HyperRectangle{4, Float64}(Vec(Inf,Inf,Inf,Inf), Vec(-Inf,-Inf,-Inf,-Inf))

    a = update(a, Vec(1,2,3,4))

    @test a == HyperRectangle{4, Float64}(Vec(1.0,2.0,3.0,4.0),Vec4f0(0.0))
    @test !( a == HyperRectangle{4, Float64}(Vec(3.0,2.0,3.0,4.0),Vec(1.0,2.0,3.0,4.0)) )

    a = update(a, Vec(5,6,7,8))

    b = HyperRectangle{4, Float64}(Vec(1.0,2.0,3.0,4.0),Vec4f0(4.0))
    @test widths(b) == Vec(4., 4., 4., 4.)
    @test a == b
    @test isequal(a,b)

    @test maximum(a) == Vec(5.0,6.0,7.0,8.0)
    @test minimum(a) == Vec(1.0,2.0,3.0,4.0)
    @test origin(a) == Vec(1.0,2.0,3.0,4.0)

    @test_throws MethodError HyperRectangle(Vec(1.0,2.0,3.0), Vec(1.0,2.0,3.0,4.0))

    @test (in(a,b) && in(b,a) && contains(a,b) && contains(b,a))

    c = HyperRectangle(Vec(1.1,2.1,3.1,4.1),Vec4f0(3.8))

    @test (!in(a,c) && in(c,a) && contains(a,c) && !contains(c,a))

    #conversion
    h = HyperRectangle(0,0,1,1)
    @test HyperRectangle{2,Float64}(h) == HyperRectangle(0.,0.,1.,1.)

    # AABB
    a = AABB(0,0,1,1)
    @test a == HyperRectangle{2,Int}(Vec(0,0),Vec(1,1))

    centered_rect = centered(HyperRectangle)
    @test centered_rect == HyperRectangle{3,Float32}(Vec3f0(-0.5),Vec3f0(1))
    centered_rect = centered(HyperRectangle{2})
    @test centered_rect == HyperRectangle{2,Float32}(Vec2f0(-0.5),Vec2f0(1))

    centered_rect = centered(HyperRectangle{2, Float64})
    @test centered_rect == HyperRectangle{2,Float64}(Vec(-0.5, -0.5),Vec(1.,1.))

    centered_rect = centered(HyperRectangle{3, Float32})
    @test centered_rect == HyperRectangle{3,Float64}(Vec(-0.5,-0.5,-0.5),Vec(1.,1.,1.))

end

# Testing split function
@testset "Testing split function" begin
    d = HyperRectangle{4, Float64}(Vec(1.0,2.0,3.0,4.0),Vec4f0(1.0))
    d1, d2 = split(d, 3, 3.5)

    @test maximum(d1)[3] == 3.5
    @test minimum(d1)[3] == 3.0
    @test maximum(d2)[3] == 4.0
    @test minimum(d2)[3] == 3.5

    d = HyperRectangle(0,0,2,2)
    d1, d2 = split(d, 1, 1)
    @test d1 == HyperRectangle(0,0,1,2)
    @test d2 == HyperRectangle(1,0,1,2)
end

@testset "Test distance functions" begin
    a = HyperRectangle(Vec(0.0,0.0), Vec(1.0, 1.0))
    b = HyperRectangle(Vec(2.0,3.0), Vec(1.0, 1.0))
    p = Vec(2.5, 1.5)

    # Rect - Rect
    @test min_dist_dim(a, b, 1) == 1.0
    @test min_dist_dim(a, b, 2) == 2.0
    @test max_dist_dim(a, b, 1) == 3.0
    @test max_dist_dim(a, b, 2) == 4.0

    @test min_euclideansq(a, b) == 5.0
    @test max_euclideansq(a, b) == 25.0
    @test minmax_euclideansq(a, b) == (5.0, 25.0)

    @test min_euclidean(a, b) == sqrt(5.0)
    @test max_euclidean(a, b) == sqrt(25.0)
    @test minmax_euclidean(a, b) == (sqrt(5.0), sqrt(25.0))

    # Rect - Point
    @test min_dist_dim(a, p, 1) == 1.5
    @test max_dist_dim(a, p, 1) == 2.5
    @test minmax_dist_dim(a, p, 1) == (1.5, 2.5)

    @test min_dist_dim(a, p, 2) == 0.5
    @test max_dist_dim(a, p, 2) == 1.5
    @test minmax_dist_dim(a, p, 2) == (0.5, 1.5)

    @test min_euclideansq(a, p) == 2.5
    @test max_euclideansq(a, p) == 8.5
    @test minmax_euclideansq(a, p) == (2.5, 8.5)

    @test min_euclidean(a, p) == sqrt(2.5)
    @test max_euclidean(a, p) == sqrt(8.5)
    @test minmax_euclidean(a, p) == (sqrt(2.5), sqrt(8.5))

    p2 = Vec(0.75, 0.75)
    @test min_dist_dim(a, p2, 1) == 0.0
    @test min_dist_dim(a, p2, 2) == 0.0

    b2  = HyperRectangle(Vec(0.25, 0.25), Vec(0.75,0.75))
    @test min_dist_dim(a, b2, 1) == 0
    @test min_dist_dim(a, b2, 2) == 0
end

@testset "set operations" begin
    a = HyperRectangle(Vec(0,0),Vec(1,1))
    b = HyperRectangle(Vec(1,1),Vec(1,1))

    @test union(a,b) == union(b,a)
    @test union(a,b) == HyperRectangle(Vec(0,0), Vec(2,2))
    @test intersect(a,b) == intersect(b,a)
    @test intersect(a,b) == HyperRectangle(Vec(1,1), Vec(0,0))
    @test diff(a,b) == a
    @test diff(b,a) == b

    c = HyperRectangle(Vec(0,0),Vec(2,2))
    d = HyperRectangle(Vec(1,1),Vec(2,2))
    @test union(c,d) == union(d,c)
    @test union(c,d) == HyperRectangle(Vec(0,0), Vec(3,3))
    @test intersect(c,d) == intersect(d,c)
    @test intersect(c,d) == HyperRectangle(Vec(1,1), Vec(1,1))
    @test diff(c,d) == c
    @test diff(d,c) == d

end

# fact relations
@testset "relations" begin
    a = HyperRectangle(Vec(0,0),Vec(1,1))
    b = HyperRectangle(Vec(1,1),Vec(1,1))
    c = HyperRectangle(Vec(0,0),Vec(2,2))
    d = HyperRectangle(Vec(1,1),Vec(2,2))
    e = HyperRectangle(Vec(3,3),Vec(1,1))
    d = HyperRectangle(Vec(0.25,0.25),Vec(0.5,0.5))
    f = HyperRectangle(Vec(0.9,0.9), Vec(1.1,1.1))

    @test (finishes(b,c) && !finishes(c,b))
    @test (!finishes(a,b))
    @test (meets(a,b) && !meets(b,a))
    @test (before(a,e) && !before(e,a))
    @test (during(d,a) && !during(a,d))
    @test (starts(a,c) && !starts(c,a))
    @test (!starts(a,b))
    @test (!overlaps(a,b) && !overlaps(b,a))
    @test (overlaps(a,f) && !overlaps(f,a))
end

@testset "point membership" begin
    a = HyperRectangle{4, Float64}(Vec(1.0,2.0,3.0,4.0),Vec4f0(4.0))
    @test (in(Vec(4,5,6,7),a) && contains(a,Vec(4,5,6,7)))
    @test (in(Vec(5,6,7,8),a) && contains(a,Vec(5,6,7,8)))
    @test (!in(Vec(6,7,8,9),a) && !contains(a,Vec(6,7,8,9)))
end

@testset "from Points" begin
    a = HyperRectangle([Point(1,1), Point(2,3), Point(4,5), Point(0,-1)])
    @test a == HyperRectangle(0, -1, 4, 6)
    a = HyperRectangle{3,Int}([Point(1,1), Point(2,3), Point(4,5), Point(0,-1)])
    @test a == HyperRectangle(0,-1,0,4,6,0)
end

@testset "transforms" begin
    t = Mat3(
        1, 0, 0,
        0, 1, 0,
        1, 2, 0
    )

    h = t * HyperRectangle(0, 0, 1, 1)
    @test h == HyperRectangle(1, 2, 1, 1)
    t = Mat{2, 2}(0, 1, 1, 0)
    h = t * HyperRectangle(0, 0, 1, 2)

    @test h == HyperRectangle(0, 0, 2, 1)
    m = Mat4f0(1.0I)
    h = centered(HyperRectangle{3, Float32})
    @test h == h

end

@testset "boundingboxes" begin
    s = Sphere(Point3f0(0), 1f0)
    f = decompose(Face{2, Int}, s, 3)
    v = decompose(Point3f0, s, 3)
    x = view(v, f)
    bb1 = HyperRectangle(x)
    bb2 = HyperRectangle(v)
    @test bb1 == bb2
end

@testset "face-orientation" begin
    cube = HyperRectangle(Vec3f0(-0.5), Vec3f0(1))

    cube_faces = decompose(Face{3,Int32}, cube)
    cube_vertices = decompose(Point{3,Float32}, cube)

    cube_tris = [cube_vertices[f] for f in cube_faces]

    normals = [cross(t[2] - t[1], t[3] - t[1]) for t in cube_tris]
    centroids = sum.(cube_tris) ./ 3f0

    for (p, n) in zip(centroids, normals)
        @test (p ⋅ n) > 0f0
    end
end

end


text:
using Particles: candidates, instantiate

# original, one-step fit method
function old_fit(p::InfiniteParticle, y, x::Int)
    # first calculate log-prior
    Δlogweight = log_prior(p.stateprior, x)
    # then update sufficient stats and convert x to an index
    stateprior, x = add(p.stateprior, x)
    0 < x ≤ length(p.components)+1 ||
        throw(ArgumentError("can't fit component $x: must be between 0 and " *
                            "$(length(p.components)+1)"))

    components = copy(p.components)
    if x ≤ length(components)
        # likelihood adjustment for old observations
        Δlogweight -= marginal_log_lhood(components[x])
    else
        push!(components, p.prior)
    end
    components[x] = add(components[x], y)
    # likelihood of new observation
    Δlogweight += marginal_log_lhood(components[x])
    weight = exp(log(p.weight) + Δlogweight)

    return InfiniteParticle(components, p, x, weight, p.prior, stateprior)
end

@testset "Updating particles" begin

    @testset "putatives vs. direct fit" begin
        p = InfiniteParticle(NormalInverseChisq(0., 1., 0.1, 0.1), 2.0)
        ys = rand(10)

        global ps_put = [p]
        global ps_oldfit = [p]
        global ps_fit = [p]
        for y in ys
            ps_put = reduce(vcat, Particles.instantiate.(putatives(p, y)) for p in ps_put)
            ps_oldfit = reduce(vcat, old_fit.(p, y, candidates(p.stateprior)) for p in ps_oldfit)
            ps_fit = reduce(vcat, fit.(p, y, candidates(p.stateprior)) for p in ps_fit)
            @test components.(ps_put) == components.(ps_fit) == components.(ps_oldfit)
        end
    end

    @testset "Weight update" begin
        # weights should remain proportional to the conditional posterior p(z|y,x)
        p = InfiniteParticle(NormalInverseChisq(0., 1., 0.1, 0.1), 2.0)
        ps = [p]
        for x in 1.:5.
            ps = Particles.instantiate.(vcat(collect.(putatives.(ps, x))...))
            @test weight.(ps) ≈ marginal_posterior.(ps)
        end
    end

end


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

@doc raw"""SubjectAccessReviewSpec is a description of the access request.  Exactly one of ResourceAuthorizationAttributes and NonResourceAuthorizationAttributes must be set

    IoK8sApiAuthorizationV1SubjectAccessReviewSpec(;
        extra=nothing,
        groups=nothing,
        nonResourceAttributes=nothing,
        resourceAttributes=nothing,
        uid=nothing,
        user=nothing,
    )

    - extra::Dict{String, Vector{String}} : Extra corresponds to the user.Info.GetExtra() method from the authenticator.  Since that is input to the authorizer it needs a reflection here.
    - groups::Vector{String} : Groups is the groups you're testing for.
    - nonResourceAttributes::IoK8sApiAuthorizationV1NonResourceAttributes : NonResourceAttributes describes information for a non-resource access request
    - resourceAttributes::IoK8sApiAuthorizationV1ResourceAttributes : ResourceAuthorizationAttributes describes information for a resource access request
    - uid::String : UID information about the requesting user.
    - user::String : User is the user you're testing for. If you specify \"User\" but not \"Groups\", then is it interpreted as \"What if User were not a member of any groups
"""
mutable struct IoK8sApiAuthorizationV1SubjectAccessReviewSpec <: SwaggerModel
    extra::Any # spec type: Union{ Nothing, Dict{String, Vector{String}} } # spec name: extra
    groups::Any # spec type: Union{ Nothing, Vector{String} } # spec name: groups
    nonResourceAttributes::Any # spec type: Union{ Nothing, IoK8sApiAuthorizationV1NonResourceAttributes } # spec name: nonResourceAttributes
    resourceAttributes::Any # spec type: Union{ Nothing, IoK8sApiAuthorizationV1ResourceAttributes } # spec name: resourceAttributes
    uid::Any # spec type: Union{ Nothing, String } # spec name: uid
    user::Any # spec type: Union{ Nothing, String } # spec name: user

    function IoK8sApiAuthorizationV1SubjectAccessReviewSpec(;extra=nothing, groups=nothing, nonResourceAttributes=nothing, resourceAttributes=nothing, uid=nothing, user=nothing)
        o = new()
        validate_property(IoK8sApiAuthorizationV1SubjectAccessReviewSpec, Symbol("extra"), extra)
        setfield!(o, Symbol("extra"), extra)
        validate_property(IoK8sApiAuthorizationV1SubjectAccessReviewSpec, Symbol("groups"), groups)
        setfield!(o, Symbol("groups"), groups)
        validate_property(IoK8sApiAuthorizationV1SubjectAccessReviewSpec, Symbol("nonResourceAttributes"), nonResourceAttributes)
        setfield!(o, Symbol("nonResourceAttributes"), nonResourceAttributes)
        validate_property(IoK8sApiAuthorizationV1SubjectAccessReviewSpec, Symbol("resourceAttributes"), resourceAttributes)
        setfield!(o, Symbol("resourceAttributes"), resourceAttributes)
        validate_property(IoK8sApiAuthorizationV1SubjectAccessReviewSpec, Symbol("uid"), uid)
        setfield!(o, Symbol("uid"), uid)
        validate_property(IoK8sApiAuthorizationV1SubjectAccessReviewSpec, Symbol("user"), user)
        setfield!(o, Symbol("user"), user)
        o
    end
end # type IoK8sApiAuthorizationV1SubjectAccessReviewSpec

const _property_map_IoK8sApiAuthorizationV1SubjectAccessReviewSpec = Dict{Symbol,Symbol}(Symbol("extra")=>Symbol("extra"), Symbol("groups")=>Symbol("groups"), Symbol("nonResourceAttributes")=>Symbol("nonResourceAttributes"), Symbol("resourceAttributes")=>Symbol("resourceAttributes"), Symbol("uid")=>Symbol("uid"), Symbol("user")=>Symbol("user"))
const _property_types_IoK8sApiAuthorizationV1SubjectAccessReviewSpec = Dict{Symbol,String}(Symbol("extra")=>"Dict{String, Vector{String}}", Symbol("groups")=>"Vector{String}", Symbol("nonResourceAttributes")=>"IoK8sApiAuthorizationV1NonResourceAttributes", Symbol("resourceAttributes")=>"IoK8sApiAuthorizationV1ResourceAttributes", Symbol("uid")=>"String", Symbol("user")=>"String")
Base.propertynames(::Type{ IoK8sApiAuthorizationV1SubjectAccessReviewSpec }) = collect(keys(_property_map_IoK8sApiAuthorizationV1SubjectAccessReviewSpec))
Swagger.property_type(::Type{ IoK8sApiAuthorizationV1SubjectAccessReviewSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAuthorizationV1SubjectAccessReviewSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiAuthorizationV1SubjectAccessReviewSpec }, property_name::Symbol) =  _property_map_IoK8sApiAuthorizationV1SubjectAccessReviewSpec[property_name]

function check_required(o::IoK8sApiAuthorizationV1SubjectAccessReviewSpec)
    true
end

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


text:
# test_convergence.jl

function get_problem(ns)
    n = m = ns
    L = laplace2d(n, m)
    A = V_cycle(L, (n, m))
    pts = 1/n:1/n:1-1/n
    b = Float64[h(x,y) for x in pts, y in pts][:]
    return (A,b)
end

@testset "Convergence for example from Briggs et. al." begin
    for n in [16 32 64 128 256 512 1024]
        (A, b) = get_problem(n)
        A.grids[1].b .= b # copy rhs
        for item in Base.Iterators.take(A,15) end # iterate
        @test A.resnorm[end] < 1/n^2
        log(A,2)
    end
end


text:
"""
# Black-box combinatorial optimization with the Cakewalk method.

## Usage
`cakewalk(objective_function::Function, M::Int, N::Int, args...)`\n
objective_function - a julia function that returns a value that represents
the objective value. This function is called by cakewalk in the following fashion:\n
`objective_value = objective_function(seq, args...)`\n
where `seq::Vector{Int}` is a sequence of integers representing some input.\n

M - each dimension can be a value between 1 and M (for example, M = 2 for binary strings).\n
N - number of dimensions.\n
args - any additional information that might be needed by objective_function.\n

## Parallelism

First, add workers by calling addprocs (for example, `addprocs(number_of_workers)`).\n
Then, make Combo available to all processes: `using Combo`.\n
Lastly, define an objective that is available to all workers. An easy way to achieve
this is by using the `@everywhere` macro.\n\n
Any additional data that is passed to cakewalk will be passed automatically to the workers by Combo.\n

## Additional important keyword parameters

minimize::Bool - if true minimize, otherwise maximize. Default: true.\n
parallel::Bool - if true, use all worker processes. Default: true.\n

### Convergence

Convergence is determined by comparing two exponentially running averages, or by
prespefying a number of samples.\n

tot_samps::Int - total number of samples. Default: 0, and is ignored.\n
min_samps::Int - minimal number of samples. \n
max_samps::Int - maximal number of samples. \n

tau_short::Float64 - time constant for the short average. \n
tau_long::Float64 - time constant for the long average. \n
tau_ratio::Float64 - ratio between the former two. Default: 2. If convergence occurs
too soon increase this value. \n
epsi::Float64 - minimal difference ratio between the two averages. When the
                     difference ratio is lower than this value, cakewalk stops.

### Verbosity
verbose::Bool - print last evaluation. Default: true. \n
verbose_fid - where to print. Default: STDOUT. \n
verbose_mod::Int - perform a print whenever number of samples mod verbose_mod is 0. Default: 1. \n
"""
function cakewalk(objective_function::Function, M::Int, N::Int, args...;
  minimize::Bool = true,

  dist::Dist = Dist(),
  comp_num::Int = max(compnum(dist), 1),

  exhaustive::Bool = false,

  tot_samps::Int = 0,
  min_samps::Int = (tot_samps == 0) ? max(M*N*comp_num, 1000) : tot_samps,
  max_samps::Int = (tot_samps == 0) ? round(Int, 1e6)*M*N : tot_samps,
  min_weight::Float = Float(.99),
  tau_short::Float = tau(min_samps,min_weight),
  tau_ratio::Float = Float(2),
  tau_long::Float = tau(min_samps*tau_ratio,min_weight),
  epsi::Float = Float(1e-2),

  opt_obj::Float = Float(Inf),

  transform_seq::Bool = false,

  max_init::Float = (comp_num == 1) ? Float(0) : Float(3),

  verbose::Bool = true,
  verbose_fid::Union{Base.TTY, IOStream, Base.PipeEndpoint} = STDOUT,
  verbose_mod::Int = Int(1),
  init_file::String = "",
  save_dir::String = "",
  save_mod::Int = Int(1e3),

  parallel::Bool = true,

  throw_exceptions::Bool = true,

  kwargs...)

  # initializations
  if isempty(init_file)
    opt_seq = zeros(Int, N)
    opt_val = (Float(Inf), Float(Inf))

    if compnum(dist) != comp_num
      init!(dist, M, N, comp_num, max_init = max_init)
    end

    avg_long = 0.
    weight_long = 0.

    avg_short = 0.
    weight_short = 0.

    seq_iter = product([1:M for i in 1:N]...)
    seq_state = start(seq_iter)

    cakewalk_log = []
  else
    opt_seq, opt_val,
    dist, avg_short, weight_short, avg_long, weight_long,
    seq_iter, seq_state,
    cakewalk_log =
    open(init_file, "r") do fid
      deserialize(fid)
    end
  end
  # for the shared array bug when using init_file
  dist = deepcopy(dist)

  function save_data(file_name::String)
    try
      open(file_name, "w") do fid
        serialize(fid,
            (opt_seq, opt_val,
            dist, avg_short, weight_short, avg_long, weight_long,
            seq_iter, seq_state,
            cakewalk_log)
        )
      end
    catch
      println("error saving $file_name")
    end
  end

  # looping
  loop = true
  rr_loop = RemoteChannel()
  set!(rr_loop, true)

  last_accepted = length(cakewalk_log)
  loop_cond = Condition()
  input_lock = ReentrantLock()
  result_lock = ReentrantLock()
  exception_lock = ReentrantLock()
  total_samps = length(cakewalk_log)
  has_savedir = !isempty(save_dir)

  worker_tasks = Tuple{Int, RemoteChannel}[]

  dist_obj = -log(eps(Float))
  diff_ratio = Float(1)

  is_better(a, b, ::Val{true}) = a < b
  is_better(a, b, ::Val{false}) = a > b
  min_type = Val{minimize}()
  min_sign = Float(-1)^minimize
  opt_obj *= min_sign

  exception = Nullable{Exception}()

  local seq, opt_samp

  function get_input()
    try
      lock(input_lock)

      if exhaustive
        seq_tuple, seq_state = next(seq_iter, seq_state)
        seq = Int[seq_tuple...]
        comp = -1
      else
        seq, comp = randseq(dist)
      end

      return (seq, comp)

    catch e
      if throw_exceptions
        throw(e)
      end
    finally
      unlock(input_lock)
    end
  end

  function send_result(res_seq, res_seq_tag, res_comp, res_val)
    try
      lock(result_lock)

      if loop
        total_samps += 1

        if is_better(res_val[2], opt_val[2], min_type) || ((res_val[2] == opt_val[2]) && is_better(res_val[1], opt_val[1], min_type))
          opt_seq = copy(res_seq_tag)
          opt_val = res_val
          opt_samp = total_samps
        end

        if !exhaustive
          dist_obj = update!(dist, transform_seq ? res_seq_tag : res_seq, res_comp, min_sign*Float(res_val[1]); kwargs...)

          avg_stat = res_val[1]

          avg_short = (1 - tau_short)*avg_short + tau_short*avg_stat
          weight_short = (1 - tau_short)*weight_short + tau_short
          avg_short_bar = avg_short/weight_short

          avg_long = (1 - tau_long)*avg_long + tau_long*avg_stat
          weight_long = (1 - tau_long)*weight_long + tau_long
          avg_long_bar = avg_long/weight_long

          diff_ratio = abs(fixnan((avg_long_bar - avg_short_bar)/avg_long_bar))
        end

        file_name = String(has_savedir && (mod(total_samps, save_mod) == 0) ? joinpath(save_dir, string("cakewalk_data_", randstring(), ".jls")) : "")

        if !isempty(file_name)
          save_data(file_name)
        end

        if verbose && (mod(total_samps, verbose_mod) == 0)
          @printf(verbose_fid, "%-14d %-14.6f %-14.6f %-14.6f %-14.6f %-14.6f %-14.6f %-s\n", total_samps, res_val[1], opt_val[1], opt_val[2], dist_obj, min(weight_short, weight_long), diff_ratio, file_name);
        end

        if !exhaustive
          if (weight_short > min_weight) && (weight_long > min_weight)
            loop = diff_ratio > epsi
          end
          if total_samps >= max_samps
            loop = false
          end
        else
          loop = total_samps < length(seq_iter)
        end

        if !loop
          set!(rr_loop, false)
          println()
          notify(loop_cond)
        end

      end

      return :ok
    catch e
      if throw_exceptions
        throw(e)
      end
    finally
      unlock(result_lock)
    end
  end

  function send_exception(res_exception)
    try
      lock(exception_lock)
      if throw_exceptions
        loop = false
        exception = Nullable{Exception}(res_exception)
        set!(rr_loop, false)
        println()
        notify(loop_cond)
      end
    finally
      unlock(exception_lock)
    end
  end
  # function references
  rr_get_input = RemoteChannel()
  rr_send_result = RemoteChannel()
  rr_send_exception = RemoteChannel()
  put!(rr_get_input, get_input)
  put!(rr_send_result, send_result)
  put!(rr_send_exception, send_exception)

  # start execution
  if verbose
    @printf(verbose_fid, "%-14s %-14s %-14s %-14s %-14s %-14s %-14s %-14s\n", "sample num", "objective", "opt obj", "opt true obj", "dist obj", "weight", "diff ratio", "file");
  end

  workers_set = parallel ? workers() : Int[myid()]

  try
    @sync for w in workers_set
      try
        @async rr_worker_task = remotecall_fetch(cakewalk_worker, w, myid(),
          throw_exceptions,
          rr_loop, rr_get_input, rr_send_result, rr_send_exception,
          objective_function, args...)
        push!(worker_tasks, (w, rr_worker_task))
      catch
      end
    end

    wait(loop_cond)

    if !isnull(exception)
      throw(get(exception))
    end

  catch e
    if throw_exceptions
      throw(e)
    end
  finally

    # exit workers
    @sync for (w, rr_worker_task) in worker_tasks
      if w in workers_set
        @async remotecall((rr) -> (fetch(rr).exception = InterruptException(); nothing), w, rr_worker_task)
      end
    end
  end

  if !isempty(save_dir)
    file_name = String(joinpath(save_dir, string("cakewalk_result_", randstring(), ".jls")));
    save_data(file_name)
    @printf("\nSaved - %s\n", file_name);
  end

  return (opt_seq, opt_val)
end

function cakewalk_worker(
    manager_id::Int,
    throw_exceptions::Bool,
    rr_loop::RemoteChannel,
    rr_get_input::RemoteChannel,
    rr_send_result::RemoteChannel,
    rr_send_exception::RemoteChannel,
    objective_function::Function,
    args...)

  exception = Nullable{Exception}()
  task =
  @schedule begin
    loop = true
    while loop && fetch(rr_loop)
      try

        # get input data
        seq, comp = remotecall_fetch(fetchcall, manager_id, rr_get_input)

        # perform evalution
        call_result = objective_function(seq, args...)

        # for functions that mutate the sequence
        if isa(call_result, Tuple) && isa(call_result[1], Vector{Int})
          seq_tag = call_result[1]
          if length(call_result) == 3
            call_val = call_result[2:3]
          else
            call_val = call_result[2]
          end
        else
          seq_tag = seq
          call_val = call_result
        end

        # always return a tuple of two values
        val = isa(call_val, Tuple) ? call_val : (call_val, call_val)

        # send back results
        remotecall_fetch(fetchcall, manager_id, rr_send_result, seq, seq_tag, comp, val)

      catch exception
        is_interrupt = isa(exception, InterruptException)
        loop = !throw_exceptions && !is_interrupt
        if !is_interrupt
          remotecall_fetch(fetchcall, manager_id, rr_send_exception, exception)
        end
      end
    end
  end

  rr_task = RemoteChannel()
  put!(rr_task, task)

  return rr_task
end


text:
# export bvecs_read, fvecs_read, ivecs_read
export fvecs_write, ivecs_write

# TODO bvecs_write
# TODO refactor these with metaprogrammming

"Write a dataset in fvecs format"
function fvecs_write(X::Matrix{Float32}, filename::AbstractString)
  d, n = size(X)
  Xd = Matrix{Float32}(undef, d+1, n)
  Xd[1, :] .= reinterpret(Float32, Int32(d))
  Xd[2:end, :] = X
  write(filename, Xd)
end

"Write a dataset in ivecs format"
function ivecs_write(X::Matrix{Int32}, filename::AbstractString)
  d, n = size(X)
  Xd = Matrix{Int32}(undef, d+1, n)
  Xd[1, :] .= Int32(d)
  Xd[2:end, :] = X
  write(filename, Xd)
end


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

mutable struct DomainStruct

    boxlo::Array{Float64,2};
    boxhi::Array{Float64,2};
    boxtilt::Array{Float64,2};
    pbc::Array{Int64,2};
    bcs::Array{Int64,2};
    triclinic::Int64;

    DomainStruct() = new();
end

function setdomain(boxhi, boxtilt, pbc=[1 1 1], bcs=[1 1 1 1 1 1])

dom = DomainStruct();

reg.boxlo = [0.0 0.0 0.0];
dom.boxhi = boxhi;
dom.boxtilt = boxtilt;
dom.pbc = pbc;
dom.bcs = bcs;

if (dom.boxtilt[1] == 0 && dom.boxtilt[2] == 0 && dom.boxtilt[3] == 0)
    dom.triclinic = 0;
else
    dom.triclinic = 1;
end

return dom

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# only need to be run once to install packages
#Pkg.clone("https://github.com/JuliaData/SplitApplyCombine.jl.git")
#Pkg.clone("https://github.com/xiaodaigh/FastGroupBy.jl.git")

using FastGroupBy, PooledArrays, Compat, BenchmarkTools

const N = Int(2e9/8);
const K = UInt(100);

const id3 = rand(Int32(1):Int32(K), N);
const id6 = rand(Int32(1):Int32(round(N/K)), N);
const v1 =  rand(Int32(1):Int32(5), N);

@elapsed sumby(id6[1:2],v1[1:2])

function b()
    id6c = copy(id6)
    v1c = copy(v1)
    @elapsed sumby(id6c,v1c)
end

res = [b() for i = 1:5]
mean(res)

function a()
    id3c = copy(id3)
    v1c = copy(v1)
    @elapsed sumby(id3c,v1c)
end

res = [a() for i = 1:10]
mean(res)

# generate string ids
function randstrarray1(pool, N)
    K = length(pool)
    PooledArray(PooledArrays.RefArray(rand(1:K, N)), pool)
enduisng

pool1 = [@sprintf "id%010d" k for k in 1:(N/K)]
id3 = randstrarray1(pool1, N)

# treat it as Pooledarray
@time sumby(id3, v1)

# treat by as strings and use dictionary method
id3_str = rand(pool1, N)
@time sumby_dict(id3_str, v1)

# parallelized sum
@time addprocs() # create Julia workers
@time using FastGroupBy
@everywhere using FastGroupBy
@everywhere using SplitApplyCombine
@time psumby(id6,v1) # 35 seconds


text:
using DiffEqFlux, Flux
import KitML

X = rand(4, 10)
Y = X.^2

# dense layers
nn = FastChain(FastDense(4, 4, tanh), FastDense(4, 4))
res = KitML.sci_train(nn, (X, Y))
res = KitML.sci_train(nn, (X, Y), res.minimizer, ADAM(); maxiters = 10000)

# fast affines
nn = FastChain(
    KitML.FastAffine(4, 4, tanh; precision = Float64),
    KitML.FastAffine(4, 4; precision = Float64),
)
res = KitML.sci_train(nn, (X, Y))
res = KitML.sci_train(nn, (X, Y), res.minimizer, ADAM(); maxiters = 10000)

# fast icnn
icnn = KitML.FastICNN(4, 4, [4, 4], tanh)
res2 = KitML.sci_train(icnn, (X, Y))
res2 = KitML.sci_train(icnn, (X, Y), res2.minimizer, ADAM(); maxiters = 10000)

# icnn
icnn2 = KitML.ICNNChain(4, 4, [4, 4], tanh)
KitML.sci_train!(icnn2, (X, Y), ADAM(), 10000, 10)


text:
"""

# Method `update_model_file`

Update Stan language model file if necessary

### Method
```julia
CmdStan.update_model_file(
  file::String,
  str::String
)
```
### Required arguments
```julia
* `file::AbstractString`                : File holding existing Stan model
* `str::AbstractString`                 : Stan model string
```

### Related help
```julia
?CmdStan.Stanmodel                 : Create a StanModel
```
"""
function update_model_file(file::AbstractString, model::AbstractString)

  model1 = parse_and_interpolate(model)
  model2 = ""
  if isfile(file)
    model2 = parse_and_interpolate(read(file, String))
    model1 != model2 && rm(file)
  end
  if model1 != model2
    println("\nFile $(file) will be updated.\n")
    write(file, model1)
  end

end

function parse_and_interpolate(model)
  newmodel = ""
  lines = split(model, "\n")
  for l in lines
    ls = String(strip(l))
    if VERSION.minor >= 3
      replace_strings = findall("#include", ls)
    else
      replace_strings = CmdStan._findall("#include", ls)
    end
    if length(replace_strings) == 1 &&
        # handle the case the include line is commented out
        length(ls) > 2 && !(ls[1:2] == "//")
      for r in replace_strings
        ls = split(strip(ls[r[end]+1:end]), " ")[1]
        func = open(f -> read(f, String), strip(ls))
        newmodel *= "   "*func*"\n"
      end
    else
      if length(replace_strings) > 1
        error("Improper number of includes in line `$l`")
      else
        newmodel *= l*"\n"
      end
    end
  end
  newmodel
end


text:
import Base.convert
using OpenMDAO
import OpenMDAO: detect_apply_nonlinear, detect_guess_nonlinear, detect_apply_linear
using CCBlade: solve

struct CCBladeResidualComp{TAF} <: OpenMDAO.AbstractImplicitComp
    num_nodes::Int
    num_radial::Int
    af::Vector{TAF}
    B::Int
    turbine::Bool
end

function CCBladeResidualComp(; num_nodes, num_radial, af, B, turbine)
    # Check if the airfoil interpolation passed is a num_radial-length array.
    try
        num_af = length(af)
        if num_af != num_radial
            throw(DomainError("af has length $num_af, but should have length $num_radial"))
        end
    catch e
        if isa(e, MethodError)
            # af is not an array of stuff, so assume it's just a single
            # function, and make it have shape (num_radial,).
            af = fill(af, num_radial)
        else
            # Some other error happened, so rethrow it.
            rethrow(e)
        end
    end

    TAF = eltype(af)
    return CCBladeResidualComp{TAF}(num_nodes, num_radial, af, B, turbine)
end

function OpenMDAO.setup(self::CCBladeResidualComp)
    num_nodes = self.num_nodes
    num_radial = self.num_radial

    input_data = [
        VarData("r", shape=(num_nodes, num_radial), val=1., units="m"),
        VarData("chord", shape=(num_nodes, num_radial), val=1., units="m"),
        VarData("theta", shape=(num_nodes, num_radial), val=1., units="rad"),
        VarData("Vx", shape=(num_nodes, num_radial), val=1., units="m/s"),
        VarData("Vy", shape=(num_nodes, num_radial), val=1., units="m/s"),
        VarData("rho", shape=(num_nodes, 1), val=1., units="kg/m**3"),
        VarData("mu", shape=(num_nodes, 1), val=1., units="N/m**2*s"),
        VarData("asound", shape=(num_nodes, 1), val=1., units="m/s"),
        VarData("Rhub", shape=(num_nodes, 1), val=1., units="m"),
        VarData("Rtip", shape=(num_nodes, 1), val=1., units="m"),
        VarData("pitch", shape=(num_nodes, 1), val=1., units="rad"),
        VarData("precone", shape=(num_nodes, 1), val=1., units="rad")]

    output_data = [
        VarData("phi", shape=(num_nodes, num_radial), val=1., units="rad"),
        VarData("Np", shape=(num_nodes, num_radial), val=1., units="N/m"),
        VarData("Tp", shape=(num_nodes, num_radial), val=1., units="N/m"),
        VarData("a", shape=(num_nodes, num_radial), val=1.),
        VarData("ap", shape=(num_nodes, num_radial), val=1.),
        VarData("u", shape=(num_nodes, num_radial), val=1., units="m/s"),
        VarData("v", shape=(num_nodes, num_radial), val=1., units="m/s"),
        VarData("W", shape=(num_nodes, num_radial), val=1., units="m/s"),
        VarData("cl", shape=(num_nodes, num_radial), val=1.),
        VarData("cd", shape=(num_nodes, num_radial), val=1.),
        VarData("cn", shape=(num_nodes, num_radial), val=1.),
        VarData("ct", shape=(num_nodes, num_radial), val=1.),
        VarData("F", shape=(num_nodes, num_radial), val=1.),
        VarData("G", shape=(num_nodes, num_radial), val=1.)]

    partials_data = Vector{PartialsData}()

    of_names = ["phi", "Np", "Tp", "a", "ap", "u", "v", "W", "cl", "cd", "cn", "ct", "F", "G"]

    ss_sizes = Dict(:i=>num_nodes, :j=>num_radial)
    rows, cols = get_rows_cols(ss_sizes=ss_sizes, of_ss=[:i, :j], wrt_ss=[:i])
    for name in of_names
        push!(partials_data, PartialsData(name, "Rhub", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "Rtip", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "pitch", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "precone", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "rho", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "mu", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "asound", rows=rows, cols=cols))
    end

    rows, cols = get_rows_cols(ss_sizes=ss_sizes, of_ss=[:i, :j], wrt_ss=[:i, :j])
    for name in of_names
        push!(partials_data, PartialsData(name, "r", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "chord", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "theta", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "Vx", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "Vy", rows=rows, cols=cols))
        push!(partials_data, PartialsData(name, "phi", rows=rows, cols=cols))
    end

    for name in ["Np", "Tp", "a", "ap", "u", "v", "W", "cl", "cd", "cn", "ct", "F", "G"]
        push!(partials_data, PartialsData(name, name, rows=rows, cols=cols, val=1.0))
    end

    return input_data, output_data, partials_data
end

function OpenMDAO.linearize!(self::CCBladeResidualComp, inputs, outputs, partials)
    num_nodes = self.num_nodes
    num_radial = self.num_radial
    af = self.af
    B = self.B
    turbine = self.turbine

    r = inputs["r"]
    chord = inputs["chord"]
    theta = inputs["theta"]
    Vx = inputs["Vx"]
    Vy = inputs["Vy"]
    rho = inputs["rho"]
    mu = inputs["mu"]
    asound = inputs["asound"]
    Rhub = inputs["Rhub"]
    Rtip = inputs["Rtip"]
    pitch = inputs["pitch"]
    precone = inputs["precone"]

    phi = outputs["phi"]

    for i in 1:num_nodes
        for j in 1:num_radial
            # Create the input structs.
            rotor = Rotor(Rhub[i, 1], Rtip[i, 1], B, turbine, pitch[i, 1], precone[i, 1])
            section = Section(r[i, j], chord[i, j], theta[i, j], af[j])
            op = OperatingPoint(Vx[i, j], Vy[i, j], rho[i, 1], mu[i, 1], asound[i, 1])

            # Find the partials of the phi residual.
            residual_deriv = residual_partials(phi[i, j], rotor, section, op)

            # Store the derivatives of the phi residual wrt phi.
            str = "phi"
            sym = Symbol(str)
            # reshape does not copy data: https://github.com/JuliaLang/julia/issues/112
            deriv = transpose(reshape(partials["phi", str], (num_radial, num_nodes)))
            deriv[i, j] = getfield(residual_deriv, sym)

            # Store the derivatives of the phi residual wrt the inputs.
            for str in keys(inputs)
                sym = Symbol(str)
                deriv = transpose(reshape(partials["phi", str], (num_radial, num_nodes)))
                deriv[i, j] = getfield(residual_deriv, sym)
            end

            # Get the derivatives of the explicit outputs.
            output_derivs = output_partials(phi[i, j], rotor, section, op)

            # Copy the derivatives of the explicit outputs into the partials dict.
            for of_str in keys(outputs)
                if of_str == "phi"
                    continue
                else
                    of_sym = Symbol(of_str)
                    for wrt_str in keys(inputs)
                        wrt_sym = Symbol(wrt_str)
                        # reshape does not copy data: https://github.com/JuliaLang/julia/issues/112
                        deriv = transpose(reshape(partials[of_str, wrt_str], (num_radial, num_nodes)))
                        deriv[i, j] = -getfield(getfield(output_derivs, of_sym), wrt_sym)
                    end
                    # Also need the derivative of each output with respect to phi, but
                    # phi is an output, not an input, so we'll need to handle that
                    # seperately.
                    wrt_str = "phi"
                    wrt_sym = Symbol(wrt_str)
                    # reshape does not copy data: https://github.com/JuliaLang/julia/issues/112
                    deriv = transpose(reshape(partials[of_str, wrt_str], (num_radial, num_nodes)))
                    deriv[i, j] = -getfield(getfield(output_derivs, of_sym), wrt_sym)
                end
            end
        end
    end

    return nothing
end

function OpenMDAO.solve_nonlinear!(self::CCBladeResidualComp, inputs, outputs)
    num_nodes = self.num_nodes
    num_radial = self.num_radial
    af = self.af
    B = self.B
    turbine = self.turbine

    r = inputs["r"]
    chord = inputs["chord"]
    theta = inputs["theta"]
    Vx = inputs["Vx"]
    Vy = inputs["Vy"]
    rho = inputs["rho"]
    mu = inputs["mu"]
    asound = inputs["asound"]
    Rhub = inputs["Rhub"]
    Rtip = inputs["Rtip"]
    pitch = inputs["pitch"]
    precone = inputs["precone"]

    for i in 1:num_nodes
        for j in 1:num_radial
            # Create the input structs.
            rotor = Rotor(Rhub[i, 1], Rtip[i, 1], B, turbine, pitch[i, 1], precone[i, 1])
            section = Section(r[i, j], chord[i, j], theta[i, j], af[j])
            op = OperatingPoint(Vx[i, j], Vy[i, j], rho[i, 1], mu[i, 1], asound[i, 1])

            # Solve the BEMT equation.
            out = solve(rotor, section, op)

            # Save the outputs.
            for str in keys(outputs)
                sym = Symbol(str)
                outputs[str][i, j] = getfield(out, sym)
            end
        end
    end

    return nothing
end

detect_apply_nonlinear(::Type{<:CCBladeResidualComp}) = false
detect_guess_nonlinear(::Type{<:CCBladeResidualComp}) = false
detect_apply_linear(::Type{<:CCBladeResidualComp}) = false


text:
struct TrialFunction{dim,T}
    fs::AbstractDiscreteFunctionSpace{dim,T}
end

@inline getfunctionspace(u::TrialFunction) = u.fs

# Test Functions

struct TestFunction{dim,T}
    fs::AbstractDiscreteFunctionSpace{dim,T}
end

@inline getfunctionspace(u::TestFunction) = u.fs

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
println("# Test Estimate of Posterior Similarity Matrix")

println(" * generate test data")

Z = hcat(vec([1 1 2 2]),
	vec([1 1 2 2]),
	vec([1 2 2 2]),
	vec([2 2 1 1]))

println(" * compute psm")

psm = compute_psm(Z)

truePSM = [1.00 0.75 0.00 0.00;
						0.75 1.00 0.25 0.25;
						0.00 0.25 1.00 1.00;
						0.00 0.25 1.00 1.00]

@test all(psm .== truePSM)

println(" * generate test data 2")

Z = hcat(vec([1 1 2 2]),
	vec([1 1 2 2]),
	vec([1 2 2 2]),
	vec([2 2 1 1]),
	vec([1 1 2 1]))

println(" * compute psm")

psm = compute_psm(Z)

truePSM = [1.0  0.8  0.0  0.2;
						0.8  1.0  0.2  0.4;
						0.0  0.2  1.0  0.8;
						0.2  0.4  0.8  1.0]

@test all(psm .== truePSM)

println("# FINISHED")


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    euclidsq(x, y)

Return the Euclidean distance between two indexable objects.
"""
euclidsq(x, y) = sum((x[i] - y[i])^2 for i in eachindex(x))

# Support matrix data in a GPU and AD compatible way

"""
    euclidsq(X::T, Y::T) where {T<:AbstractMatrix}

Return the Euclidean distance between columns in two matrices.
"""
function euclidsq(X::T, Y::T) where {T<:AbstractMatrix}
    XiXj = transpose(X) * Y
    x² = sum(X .^ 2; dims=1)
    y² = sum(Y .^ 2; dims=1)
    transpose(x²) .+ y² - 2XiXj
end

"""
    euclidsq(X::T) where {T<:AbstractMatrix}

Effective version of `euclidsq(X, X)`.
"""
function euclidsq(X::T) where {T<:AbstractMatrix}
    XiXj = transpose(X) * X
    x² = sum(X .^ 2; dims=1)
    transpose(x²) .+ x² - 2XiXj
end

"""
    gaussian_gramian(xs, ys, σ=1)

Gramian matrix for samples `xs` and `ys` using a Gaussian kernel
kernel with bandwidth `σ`.
"""
gaussian_gramian(xs; kwargs...) = gaussian_gramian(xs, xs; kwargs...)
gaussian_gramian(xs, ys; σ=1) =
  [exp(-euclidsq(x, y) / 2σ^2) for x in xs, y in ys]

function gaussian_gramian(X::T, Y::T; σ=1) where {T<:AbstractMatrix}
    gaussian_gramian(euclidsq(X, Y), σ)
end

gaussian_gramian(esq, σ::AbstractFloat) = exp.(-esq ./ 2σ^2)

"""
    safe_diagm(mat, a)

Generate a squared matrix whose diagonal is `a` that is
compatible to perform addition on `mat`. It hebaves
differently based on `mat` is on CPU or GPU.

It is compatible with
- CuArrays.jl (see lib/cuarrays.jl)
- Zygote.jl (see lib/zygote.jl)
"""
safe_diagm(mat, a) = a * I


text:
using Test
using UnivariateGaussianProcess

global atol = 1e-8

@testset "Acquisition function optimization" begin
    @test UnivariateGaussianProcess.find_local_opt(x -> x^2, -1., 1.)[1] ≈ 0. atol=atol
    @test UnivariateGaussianProcess.find_local_opt(x -> -sin(x), 0., π)[1] ≈ π/2 atol=atol
end

@testset "Kernel computation" begin
    cov_func = (x1,x2) -> abs(x2 - x1)
    data_x = [1., 2., 3.]
    noise = 0.
    @test UnivariateGaussianProcess.create_kernel(cov_func, data_x; noise) ≈ [
        0. 1. 2.
        1. 0. 1.
        2. 1. 0.
    ] atol=atol

    cov_func = (x1,x2) -> abs(x2 - x1)
    data_x = [1., 2., 3.]
    noise = 0.1
    @test UnivariateGaussianProcess.create_kernel(cov_func, data_x; noise) ≈ [
        0.01 1. 2.
        1. 0.01 1.
        2. 1. 0.01
    ] atol=atol

    cov_func = (x1,x2) -> abs(x2 - x1)
    K = [0.01 1. 2.
         1. 0.01 1.
         2. 1. 0.01]
    old_data_x = [1., 2., 3.]
    new_x = 5.
    noise = 0.1
    @test UnivariateGaussianProcess.update_kernel(cov_func, K, old_data_x, new_x; noise) ≈ [
        0.01 1. 2. 4.
        1. 0.01 1. 3.
        2. 1. 0.01 2.
        4. 3. 2. 0.01
    ] atol=atol
end


text:
module Mixtures

using Roots, SpecialFunctions, StatsFuns
using Random, Distributions, Statistics

export fitmixture, EM

function MLE(x; γ=missing)
    if ismissing(γ)
        γ = ones(size(x))
    end
    N = sum(γ)
    μ = sum((γ.*x))/N
    σ = log(μ) - sum(γ.*log.(x))/N

    # initial guess
    α = ((3 - σ) + sqrt((σ-3)^2 + 24σ)) / (12σ)
    if isinf(α)
        α = 10
    end

    try
        F(α) = log(α) - digamma(α) - σ
        α = find_zero(F, α, Order1())
    catch
        println("ERROR:")
        println("bounds=($(minimum(x)), $(maximum(x)))")
        println("N=$(N)")
        println("σ=$(σ)")
        println("μ=$(μ)")
        println(((3 - σ) + sqrt((σ-3)^2 + 24σ)) / (12σ))
    end
    β = α/μ

    return α, β
end

function weights(γ)
    return sum(γ, dims=1)'/sum(γ)
end

logρ(x,α,β) = (α*log(β) - loggamma(α)) .+ (α-1).*log.(x) .- β.*x
ρ(x,α,β)    = (x != 0) ? exp(logρ(x,α,β)) : β^α / gamma(α)

function posterior(x, Π, α, β; dropout=false)
    γ = Array{Float64, 2}(undef, length(x), length(α))

    for k = 1:length(α)
        γ[:,k] = Π[k]*ρ.(x, α[k], β[k])
    end

    if dropout
        χ = minimum(x[x.>0])
        p = vec([gammacdf(α[k], 1/β[k], χ) for k in 1:length(α)])
        γ[findall(x.==0),:] .= p'
    end

    if any(isinf.(γ))
        γ[isinf.(γ)] .= 1
    end

    γ = γ ./ sum(γ, dims=2)

    return γ
end

loglikelihood(x, Π, α, β) = sum(log.(sum(Π' .* hcat([ρ.(x,α[k],β[k]) for k in 1:length(α)]...), dims=2)))

function sample(Π, α, β, N; k=missing)
    if ismissing(k)
        MM = MixtureModel(Gamma, [(α[k], 1/β[k]) for k in 1:length(α)], vec(Π))
    else
        MM = Gamma(α[k], 1/β[k])
    end
    return rand(MM, N)
end

function kmeans(x, k::Integer; verbose = false, maxᵢ::Integer = 1000)
    q  = [quantile(x, q) for q in range(0, 1, length=k+1)]
    μ  = [mean(x[(q[i-1] .<= x) .& (x .<= q[i])]) for i in 2:length(q)]
    S₀ = zeros(length(x), 1)
    for n in 1:maxᵢ
        S = getindex.(argmin(abs.(x .- μ'), dims=2), 2)
        if !any(x->x!=0, S-S₀)
            break
        end
        if verbose
            println("iteration $(n): μ = $(μ)")
        end
        μ  = [mean(x[findall(S.==i)]) for i in 1:k]
        S₀ = S
    end

    σ = [std(x[findall(S₀.==i)]) for i in 1:k]

    return μ, σ, S₀
end

function EM(x, K::Integer; ϵ::Real = 1e-5, maxᵢ::Integer = 10000, verbose::Bool = false)
    μ, σ, S₀ = kmeans(x, K)

    β = μ./σ.^2
    α = μ.*β
    Π = [sum(S₀.==k)/length(S₀) for k in 1:K]

    δ  = Inf
    L₀ = loglikelihood(x, Π, α, β)
    i  = 0
    while abs(δ) > ϵ && i < maxᵢ
        # expectation + maximization
        γ = posterior(x, Π, α, β)

        Π = weights(γ)
        Θ = [MLE(x, γ=γ[:,k]) for k in 1:K]
        α = [θ[1] for θ in Θ]
        β = [θ[2] for θ in Θ]

        L = loglikelihood(x, Π, α, β)
        if verbose
            println("iteration $(i). loglikelihood = $(L)")
        end
        δ  = L - L₀
        L₀ = L
        i += 1
    end

    ι = sortperm(α./β)

    return L₀, (Π[ι], α[ι], β[ι])
end

function cluster(x, Π, α, β; dropout=false)
    println("computing posterior")
    γ = posterior(x, Π, α, β, dropout=dropout)
    println("getting index")
    χ = getindex.(argmax(γ, dims=2), 2)
    println("returning")
    return [x[findall(χ .== k)] for k in 1:length(α)]
end

# XXX: right now just returns the posteriors for component with largest mean
function fitmixture(data)
    _, params = EM(data.+1e-6, 2)
    ρ = posterior(data, params...)
    return ρ[:,2]
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
export MembershipTestResult, MembershipTestOptions, membershiptest

"""
    MembershipTestResult{T, N}

**Fields**
* `successfull::Bool`
* `tries::Int` The number of times Newton's method was started
* `converged_iterations::Int` Number iterations needed in a successfull iteration. If the test was not successfull `imax` is the maximal number of iterations.
* `startvalue::SVector{N, T}` The start value of the successfull iteration. Otherwise a random value.
* `solution::SVector{N, T}` The solution value of the successfull iteration. Otherwise a random value.
"""
struct MembershipTestResult{T, N}
    successfull::Bool
    tries::Int
    converged_iterations::Int
    startvalue::SVector{N, T}
    solution::SVector{N, T}
end

"""
    MembershipTestOptions(iterations=50, ntries=30, tol=1e-7)

# Details
* `iterations` The maximal number of iterations in one try.
* `ntries` The number of times Newton's method will be tried
* `tol` The desired accuracy when Newton's method converges.
"""
struct MembershipTestOptions
    maxiters::Int
    tol::Float64
    ntries::Int
end
function MembershipTestOptions(;maxiters=50, ntries=30, tol=1e-7)
    MembershipTestOptions(maxiters, tol, ntries)
end

"""
    membershiptest(F::AbstractFiber{T}, w:, options=MembershipTestOptions()) where T

Check wether the point `w` is contained in the fiber `F`.
The tests uses Newton's method with configuration provided by `options`.
Returns a [`MembershipTestResult`](@ref).
"""
function membershiptest(F, w, gen, o=MembershipTestOptions())
    membershiptest(F, w, gen, (startvalue(gen),), o)
end
function membershiptest(F::AbstractFiber{T}, w, gen::StartValueGenerator, initial_startvalues, o::MembershipTestOptions) where T
    update_fiber!(F, w)
    k = 0
    i = 1
    if isempty(initial_startvalues)
        θ = startvalue(gen)
    else
        θ = initial_startvalues[i]
    end

    maxiters = o.maxiters
    while true
        k += 1
        # try
        res = newton(F, θ, maxiters, o.tol)
        if res.converged
            return MembershipTestResult(true, k, res.iterations, θ, res.solution)
        # If it seems that we could get a zero we iterate a little bit more
        # We should only need a small amount of iterations due to the
        # locally quadratically convergence of the newton iteration
        elseif res.residual < 1e-2
            refined = newton(F, res.solution, 10, o.tol)
            if refined.converged
                return MembershipTestResult(true, k, maxiters + refined.iterations, θ, res.solution)
            end
        end
        # catch err
        #     @show err
        # end
        if k > o.ntries
            break
        end
        if i < length(initial_startvalues)
            i += 1
            θ = initial_startvalues[i]
        else
            maxiters = o.maxiters
            θ = startvalue(gen)
        end
    end
    return MembershipTestResult(false, o.ntries, maxiters, θ, θ)
end


text:
module JavisTimer
import Base

export timer, finterval

framerate = 30
currentframe = 1

"""
module contains some functions for manipulating time directly
not that in general manipulating time is tricky because of
discrete nature of frames.

"""

struct finterval #frame interval
  """
  an interval of frames
  iterable.

  finterval(b::Int,e::Int)
    makes an interval starting at the b'th frame
    and ending at the e'th frame.

  finterval(e::Int)
    makes an interval starting at `currentframe`
    and ending  `currentframe+e` frame
  """
  beginframe::Int
  endframe::Int
  iterable::Any
end

finterval(b::Int, e::Int) = finterval(b, e, b:e)
finterval(e) = finterval(currentframe, e, currentframe:e)

"""implement iteration on finterval"""
Base.iterate(f1::finterval) = Base.iterate(f1.iterable)
function Base.iterate(f1::finterval, state)
  return Base.iterate(f1.iterable, state)
end

function Base.:+(f1::finterval, f2::finterval)
  """adds two intervals, addition is defined as an finterval from
  f1's first frame to f2's last frame
  """
  return finterval(f1.beginframe, f2.endframe)
end

function Base.:*(x::Real, f1::finterval)
  """scales the interval by x , keeps the beginframe the same"""
  return finterval(f1.beginframe, floor(Int, x * f1.endframe))
end
#another other ideas for * on finterval is to scale around midpoint frame

function Base.length(f::finterval)
  return length(f.iterable)
end

function timer(t::Real)
  """advance currentframe by t seconds of time,
  note that the currentframe is obtained by
  flooring to an `Int`,
  do not expect running `timer(x);timer(y)` to add `x+y` seconds
  to the current timestamp."""
  global currentframe
  beginframe = currentframe
  endframe = currentframe + floor(Int, t * framerate) - 1
  currentframe = endframe + 1
  return finterval(beginframe, endframe)
end

function timer(f::finterval)
  """advance currentframe by the duration same as f,
  and return apropriate finterval"""
  global currentframe, beginframe = currentframe
  Δframe = f.endframe - f.beginframe + 1
  endframe = currentframe + Δframe
  return finterval(beginframe, endframe)
end

function getframe()
  """ returns currentframe """
  global currentframe
  return currentframe
end

function setframe(f::Int)
  """ sets current frame to f """
  global currentframe
  currentframe = f
end

function gettime()
  """gets timestamp, time is tricky to work with.
  """
  global currentframe
  return currentframe * framerate
end

function settime(t::Real)
  global currentframe
  currentframe = floor(Int, s * framerate)
end

function settime(f::finterval)
  """set currentframe to end frame of finterval,
  """
  global currentframe
  currentframe = f1.endframe
end

end

max_stars_repo_name: JuliaAnimators/JavisModes

max_issues_repo_name: JuliaAnimators/JavisModes

max_forks_repo_name: JuliaAnimators/JavisModes


text:
"""
    WaveEquationNonperiodicSemidiscretization(D, left_bc, right_bc)

A semidiscretization of the linear wave equation
    \$\\partial_t^2 u(t,x) = \\partial_x^2 u(t,x)\$.

`D` is assumed to be a second-derivative SBP operator and the boundary conditions
can be `Val(:HomogeneousNeumann)`, `Val(:HomogeneousDirichlet)`, or
`Val(:NonReflecting)`.
"""
struct WaveEquationNonperiodicSemidiscretization{Derivative<:AbstractDerivativeOperator,
                                                 LeftBC, RightBC} <: AbstractSemidiscretization
    derivative::Derivative
    left_bc::LeftBC
    right_bc::RightBC
end

function semidiscretize(du0func, u0func, semi::WaveEquationNonperiodicSemidiscretization, tspan)
  du0 = compute_coefficients(du0func, semi.derivative)
  u0 = compute_coefficients(u0func, semi.derivative)
  ode = SecondOrderODEProblem(semi, du0, u0, tspan)
end

function Base.show(io::IO, semi::WaveEquationNonperiodicSemidiscretization)
    if get(io, :compact, false)
        print(io, "Semidiscretization of the linear wave equation (non-periodic)")
    else
        print(io, "Semidiscretization of the linear wave equation\n")
        print(io, "  \$ \\partial_t^2 u(t,x) = \\partial_x^2 u(t,x) \$ \n")
        print(io, "with nonperiodic boundaries using")
        print(io, semi.derivative)
    end
end

function (disc::WaveEquationNonperiodicSemidiscretization)(ddu, du, u, p, t)
    @unpack derivative, left_bc, right_bc = disc
    @boundscheck begin
        @argcheck length(u) == length(du)
        @argcheck length(u) == length(ddu)
    end

    # volume terms of the second derivative
    mul!(ddu, derivative, u)

    # boundary conditions using SATs
    if left_bc == Val(:HomogeneousNeumann)
        @inbounds ddu[1] += derivative_left(derivative, u, Val(1)) / left_boundary_weight(derivative)
    elseif left_bc == Val(:HomogeneousDirichlet)
        mul_transpose_derivative_left!(ddu, derivative, Val(1), -u[1] / left_boundary_weight(derivative), true)
        # only a rough estimate, should probably be improved
        # see Mattsson, Ham, Iaccarino (2009) Stable boundary ...
        # parameter α in Table 1 instead of accuracy_order(derivative)
        @inbounds ddu[1] -= accuracy_order(derivative) * u[1] / left_boundary_weight(derivative)^2
    elseif left_bc == Val(:NonReflecting)
        @inbounds ddu[1] -= (du[1] - derivative_left(derivative, u, Val(1))) / left_boundary_weight(derivative)
    else
        throw(ArgumentError("Boundary condition $left_bc not implemented."))
    end

    if right_bc == Val(:HomogeneousNeumann)
       @inbounds ddu[end] -= derivative_right(derivative, u, Val(1)) / right_boundary_weight(derivative)
    elseif right_bc == Val(:HomogeneousDirichlet)
        mul_transpose_derivative_right!(ddu, derivative, Val(1), u[end] / right_boundary_weight(derivative), true)
        # only a rough estimate, should probably be improved, see above
        @inbounds ddu[end] -= accuracy_order(derivative) * u[end] / right_boundary_weight(derivative)^2
    elseif right_bc == Val(:NonReflecting)
        @inbounds ddu[end] -= (du[end] + derivative_right(derivative, u, Val(1))) / right_boundary_weight(derivative)
    else
        throw(ArgumentError("Boundary condition $right_bc not implemented."))
    end

    nothing
end


text:
module Day1

using Advent: save_input

function solve_a(filename, target)
    data = [parse(Int, i) for i in eachline(filename)]
    for i in Iterators.product(data, data)
        sum(i) == target && return prod(i)
    end
end

function solve_b(filename, target)
    data = [parse(Int, i) for i in eachline(filename)]
    for i in Iterators.product(data, data, data)
        sum(i) == target && return prod(i)
    end
end

solve(filename = save_input(2020, 1), t = 2020) = solve_a(filename, t), solve_b(filename, t)
solve_big() = solve("day1_big.in", 99920044)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
# This file includes the MonteCarlo simulation of ASK modulation scheme and compares
# the numerical results with the theoretical results.

using DigitalCommunications
using Plots

# Simulation parameters
k = 3
M = 2^k
nsymbols = Int(1e6)
nbits = k * nsymbols
Ep = 100.                       # Modulating pulse energy.
ebno = collect(0 : 10)          # Snr per bit
esno = ebno .+ 10 * log10(k)    # Snr ber symbol

# Communcation system components
gen = BitGenerator(nbits)
coding = GrayCoding(M)
modulator = Modulator(PAM(M))
channel = AWGNChannel()
detector = Detector(alphabet(modulator))

# Monte Carlo simulation
message = coding(gen.bits)  # Message signal
symerr = zeros(length(esno))
for i in 1 : length(symerr)
    channel.esno = esno[i]  # Update channel snr
    mbar = gen.bits |> coding |> modulator |> channel |> detector  # Extracted message signal
    symerr[i] = sum(mbar .!= message) / length(message)  # Symbol error rate
end

# Plots
plt = plot(title="$M-PAM", xlabel="ebno [dB]", ylabel="Pe")
plot!(ebno, berpam.(ebno, M), marker=:circle, yscale=:log10, label="theoretical")
plot!(ebno, symerr, marker=:circle, yscale=:log10, label="montecarlo")


text:
module HiOscSolver

include("interface.jl")
export HiOscODEProblem, HiOscInterpolation, AbstractHiOscSolution, HiOscODESolution
export solve, getexactsol
end # module


text:
import Graphs: inclist, ExVertex, ExEdge, add_vertex!, add_edge!, AttributeDict, num_vertices, num_edges, GenericIncidenceList
import Base: getindex, setindex!

export indexof, vertexof, edgeof
export for_vertices, for_edges, weights

import Graphs: to_dot
export         to_dot

import Base: +
export       +

typealias Label AbstractString

type GraphStore
  graph::GenericIncidenceList
  attributes::AttributeDict
  GraphStore(graph) = new(graph, AttributeDict())
end

getindex(g::GraphStore, label::Label) = vertexof(g, label)
function vertexof(g::GraphStore, label::Label)
  for v in g.graph.vertices
    v.label == label && return v
  end
  nothing
end

function edgeof(g::GraphStore, from::Label, to::Label)
  for el in g.graph.inclist
    for e::ExEdge{ExVertex} in el
      e.source.label == from && e.target.label == to && return e
    end
  end
  nothing
end

function indexof(g::GraphStore, label::Label)
  vtx = vertexof(g, label)
  isa(vtx, Void) ? 0 : vtx.index
end

for_vertices(fn::Function, g::GraphStore) = [fn(v) for v in g.graph.vertices]
function for_edges(fn::Function, g::GraphStore)
  results = []
  for edges in g.graph.inclist
    for edge in edges
      push!(results, fn(edge))
    end
  end
  results
end

function weights(g::GraphStore)
  for_edges(g) do e
    haskey(e.attributes, "weight") ? e.attributes["weight"] : 0
  end
end

function add_vertex!(g::GraphStore, label::Label)
  vtx = vertexof(g, label)
  if isa(vtx, Void)
    cnt = num_vertices(g.graph)+1
    vtx = ExVertex(cnt, label)
    vtx.attributes["label"] = label
    add_vertex!(g.graph, vtx)
  else
    vtx
  end
end

function +(a::GraphStore, b::GraphStore)
  for el in b.graph.inclist
    for e::ExEdge{ExVertex} in el
      vtxfrom, vtxto = add_vertex!(a, e.source.label), add_vertex!(a, e.target.label)
      edge = ExEdge{ExVertex}(num_edges(a.graph), vtxfrom, vtxto, e.attributes)
      add_edge!(a.graph, edge)
    end
  end
  if !issubset(map(x->x.label, b.graph.vertices), map(x->x.label, a.graph.vertices))
    for v in b.graph.vertices
      vtx = vertexof(a, v.label)
      isa(vtx, Void) && add_vertex!(a, v.label)
    end
  end
  if a.graph.is_directed != b.graph.is_directed
    a.graph.is_directed = true
  end
  a
end

function setindex!(ex::Union{ExVertex,ExEdge}, value, attr)
  ex.attributes[attr] = value
end

include("parse.jl")
macro graph(args...)
  list = []
  for el in args
    list = vcat(list, isa(el, Expr) ? adjacency(el.args...) : isa(el, Symbol) ? string(el) : el)
  end
  parse(list)
end

include("draw.jl")


text:
# lnprob evaluation
# for a single process only, designed for profiling purposes
# that means the starting point is receiving the pars::Parameters object
# and the end point is returning a lnprop
using Pkg; Pkg.activate("DiskJockey")

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--run_index", "-r"
    help = "Output run index"
    arg_type = Int
    "config"
    help = "a YAML configuration file"
    required = true
end

parsed_args = parse_args(ARGS, s)

import YAML
cfg = YAML.load(open(parsed_args["config"]))
outfmt(run_index::Int) = cfg["out_base"] * @sprintf("run%02d/", run_index)
basefmt(id::Int) = cfg["base_dir"] * @sprintf("run%02d/", id)

# This code is necessary for multiple simultaneous runs on odyssey
# so that different runs do not write into the same output directory
if parsed_args["run_index"] == nothing
    run_index = 0
    outdir = outfmt(run_index)
    while ispath(outdir)
        println(outdir, " exists")
        run_index += 1
        outdir = outfmt(run_index)
    end
else
    run_index = parsed_args["run_index"]
    outdir = outfmt(run_index)
    println("Deleting old $outdir")
    run(`rm -rf $outdir`)
end

# make the directories
println("Creating ", outdir)
mkdir(outdir)

@everywhere basefmt(id::Int) = cfg["base_dir"] * @sprintf("run%02d/", id)

@everywhere const global basedir = basefmt(run_index)

# make the internal Judith directory, if it doesn't exist
if !ispath(basedir)
    println("Creating ", basedir)
    mkdir(basedir)
end

using constants
using visibilities
using image
using gridding
using model

# This is the likelihood function called by each individual process
function f(data, keys::Vector{Int}, p::Parameters)
    # split into datasets and interpolation functions
    dvarr, int_arr = data

    nkeys = length(keys)

    # We are using the following convention: inclination ranges from
    # 0 to 180 degrees. 0 means face on, angular momentum vector pointing
    # at observer; 90 means edge on; and 180 means face on, angular momentum
    # vector pointing away from observer.
    # These are also the RADMC conventions.
    incl = p.incl # [deg]

    # We also adopt the RADMC convention for position angle, which defines position angle
    # by the angular momentum vector.
    # A positive PA angle means the disk angular momentum vector will be
    # rotated counter clockwise (from North towards East).
    PA = p.PA # [deg]

    vel = p.vel # [km/s]
    npix = cfg["npix"] # number of pixels, can alternatively specify x and y separately

    # Doppler shift the dataset wavelengths to rest-frame wavelength
    beta = vel/c_kms # relativistic Doppler formula
    lams = Array{Float64}(length(keys))
    for i=1:nkeys
        lams[i] =  dvarr[i].lam * sqrt((1. - beta) / (1. + beta)) # [microns]
    end

    # key = keys[1]
    # keydir = basedir * "jud$key/"
    # println("Writing cameralambda into $keydir")
    write_lambda(lams, "") # write into current directory

    # Run RADMC3D, redirect output to /dev/null
    run(`radmc3d image incl $incl posang $PA npix $npix loadlambda` |> DevNull)

    # Read the RADMC3D images from disk (we should already be in sub-directory)
    im = imread()

    # Convert raw images to the appropriate distance
    skim = imToSky(im, p.dpc)

    # Apply the gridding correction function before doing the FFT
    # shifts necessary as if the image were already offset
    corrfun!(skim, p.mu_RA, p.mu_DEC) # alpha = 1.0

    lnprobs = Array{Float64}(nkeys)
    # Do the Fourier domain stuff per channel
    for i=1:nkeys
        dv = dvarr[i]
        # FFT the appropriate image channel
        vis_fft = transform(skim, i)

        # Apply the phase correction here, before we do the interpolation
        phase_shift!(vis_fft, p.mu_RA, p.mu_DEC)

        # Interpolate the `vis_fft` to the same locations as the DataSet
        # mvis = ModelVis(dv, vis_fft)

        mvis = int_arr[i](dv, vis_fft)

        # Calculate chi^2 between these two
        lnprobs[i] = lnprob(dv, mvis)
    end

    # Sum them all together and feed back to the master process
    return sum(lnprobs)
end

# Create the model grid
grd = cfg["grid"]
global const grid = Grid(grd["nr"], grd["ntheta"], grd["r_in"], grd["r_out"], true)

# Regenerate all of the static files (e.g., amr_grid.inp)
# so that they may be later copied
write_grid(basedir, grid)

# Essentially, `mach_three.initfunc` is replicated here, since we'll just be profiling a single channel

keys = [11, 12, 13, 14, 15, 16]

@everywhere function initfunc(keys)
    # In contrast to `burma_shave.jl`, now each dataset receives several
    # channels, specified by keys

    # Load the relevant chunk of the dataset, return array
    dvarr = DataVis(cfg["data_file"], keys)
    for dset in dvarr
        # Conjugation is necessary for the SMA, methinks
        visibilities.conj!(dset) # Swap UV convention
    end

    # Create a directory where all RADMC files will reside and be driven from
    # using the first key as the directory index
    key = keys[1]
    keydir = basedir * "jud$key"
    println("Creating $keydir")
    mkdir(keydir)

    # Copy all relevant configuration scripts to this subdirectory
    # these are mainly setup files which will not change throughout the run
    run(`cp radmc3d.inp $keydir`)
    ag = basedir * "amr_grid.inp"
    run(`cp $ag $keydir`)
    run(`cp lines.inp $keydir`)
    run(`cp molecule_co.inp $keydir`)
    run(`cp wavelength_micron.inp $keydir`)

    # change the subprocess to reside in this directory for the remainder of the run
    # where it will drive its own independent RADMC3D process
    cd(keydir)

    # For each channel, also calculate the interpolation closures
    npix = cfg["npix"]
    pixsize = cfg["pixsize"] # [cm]
    pix_AU = pixsize/AU # [AU]

    dl = sin(pix_AU/cfg["parameters"]["dpc"][1] * arcsec)

    uu = fftshift(fftfreq(npix, dl)) * 1e-3 # [kλ]
    vv = fftshift(fftfreq(npix, dl)) * 1e-3 # [kλ]

    int_arr = Array{Function}(length(keys))
    for (i, dset) in enumerate(dvarr)
        int_arr[i] = plan_interpolate(dset, uu, vv)
    end

    # return the array of datasets
    return (dvarr, int_arr)
end

# Load the relevant chunk of the dataset, return array
data = initfunc(keys)

# READ THE PARAMETERS FROM DISK
pp = cfg["parameters"]
# The parameters we'll be using
params = ["M_star", "r_c", "T_10", "q", "logM_CO", "ksi", "incl", "PA", "vel", "mu_RA", "mu_DEC"]
nparam = length(params)
p = Array{Float64}(nparam)
jumps = Array{Float64}(nparam)

for i=1:nparam
    p[i], jumps[i] = pp[params[i]]
end

# FPROB START
gamma = 1.0 # surface temperature gradient exponent
dpc = cfg["parameters"]["dpc"][1] # [pc] distance

M_star, r_c, T_10, q, logM_CO, ksi, incl, PA, vel, mu_RA, mu_DEC = p

# Enforce hard priors on physical parameters
# Short circuit evaluation if we know the RADMC won't be valid.
if ksi <= 0. || T_10 <= 0. || r_c <= 0.0 || M_star <= 0.0
    return -Inf
end

if incl < 0. || incl > 180.
    return -Inf
end

M_CO = 10^logM_CO

# If we are going to fit with some parameters dropped out, here's the place to do it
pars = Parameters(M_star, r_c, T_10, q, gamma, M_CO, ksi, dpc, incl, PA, vel, mu_RA, mu_DEC)

# Compute parameter file using model.jl, write to disk
write_model(pars, basedir, grid)

# Time it again
tic()
write_model(pars, basedir, grid)
toc()

nd = basedir * "numberdens_co.inp"
gv = basedir * "gas_velocity.inp"
gt = basedir * "gas_temperature.inp"
mt = basedir * "microturbulence.inp"

# Copy new parameter files to all subdirectories

key = keys[1]
keydir = basedir * "jud$key"
run(`cp $nd $keydir`)
run(`cp $gv $keydir`)
run(`cp $gt $keydir`)
run(`cp $mt $keydir`)

# FPROB DONE

# given these made up parameters, call fprob
println(f(data, keys, pars))

println("Now for the profiling.")

clear_malloc_data()
f(data, keys, pars)

@profile f(data, keys, pars)
Profile.print()

println("f Timing")
tic()
f(data, keys, pars)
toc()


text:
using TerminalUserInterfaces: Block, Buffer, Rect, draw, inner, set, height

Base.@kwdef struct Text
    block::Block
    text::String
end

function TerminalUserInterfaces.draw(t::Text, rectouter::Rect, buf::Buffer)
    draw(t.block, rectouter, buf)
    rect = inner(t.block, rectouter)
    height(rect) < 1 && return

    col, row = rect.x, rect.y

    for char in t.text
        if char == '\n'
            row += 1
            col = rect.x
            continue
        end
        if col >= 200#(rect.x + rect.width - 1)
            row += 1
            col = rect.x
        end
        if row > (rect.y + rect.height)
            break
        end
        set(buf, col, row, char)
        col += 1
    end
end


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

mutable struct ApplicationGatewayUrlConfiguration <: SwaggerModel
    modifiedPath::Any # spec type: Union{ Nothing, String } # spec name: modifiedPath
    modifiedQueryString::Any # spec type: Union{ Nothing, String } # spec name: modifiedQueryString
    reroute::Any # spec type: Union{ Nothing, Bool } # spec name: reroute

    function ApplicationGatewayUrlConfiguration(;modifiedPath=nothing, modifiedQueryString=nothing, reroute=nothing)
        o = new()
        validate_property(ApplicationGatewayUrlConfiguration, Symbol("modifiedPath"), modifiedPath)
        setfield!(o, Symbol("modifiedPath"), modifiedPath)
        validate_property(ApplicationGatewayUrlConfiguration, Symbol("modifiedQueryString"), modifiedQueryString)
        setfield!(o, Symbol("modifiedQueryString"), modifiedQueryString)
        validate_property(ApplicationGatewayUrlConfiguration, Symbol("reroute"), reroute)
        setfield!(o, Symbol("reroute"), reroute)
        o
    end
end # type ApplicationGatewayUrlConfiguration

const _property_map_ApplicationGatewayUrlConfiguration = Dict{Symbol,Symbol}(Symbol("modifiedPath")=>Symbol("modifiedPath"), Symbol("modifiedQueryString")=>Symbol("modifiedQueryString"), Symbol("reroute")=>Symbol("reroute"))
const _property_types_ApplicationGatewayUrlConfiguration = Dict{Symbol,String}(Symbol("modifiedPath")=>"String", Symbol("modifiedQueryString")=>"String", Symbol("reroute")=>"Bool")
Base.propertynames(::Type{ ApplicationGatewayUrlConfiguration }) = collect(keys(_property_map_ApplicationGatewayUrlConfiguration))
Swagger.property_type(::Type{ ApplicationGatewayUrlConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayUrlConfiguration[name]))}
Swagger.field_name(::Type{ ApplicationGatewayUrlConfiguration }, property_name::Symbol) =  _property_map_ApplicationGatewayUrlConfiguration[property_name]

function check_required(o::ApplicationGatewayUrlConfiguration)
    true
end

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


text:
include(joinpath(AS.topdirectory, "test", "EpemaData.jl"))

wing = EpemaData.wing
rotor = EpemaData.rotor
setup = EpemaData.setup
results = EpemaData.results

Vinf = setup[:"Vinf"]
Vref = Vinf
vinfs = [Vinf]
plotstepi = [1]
alphas = [setup[:"wing_aoa"]] .* pi/180
wing_b = wing[:"span"]
wing_TR = wing["chord_tip"] / wing[:"chord_root"]
wing_AR = wing["span"]^2 / wing[:"area"]
wing_θroot = 0.0
wing_θtip = 0.0
rotor_omegas = [AS.get_omega(Vinf=setup[:"Vinf"], J=setup[:"J"], D=rotor[:"diameter"])]
nblades = [rotor[:"Nblades"]]
rhub = [rotor[:"radius_hub"]]
rtip = [rotor[:"radius"]]
rotor_chord_data = rotor["r/R vs chord/R"]
rotor_twist_data = rotor["r/R vs twist"]
radii = radii_list
rotor_chords = chords_list
rotor_twists = twists_list
# radii = [rotor_chord_data[:,1] * rtip[1]]
# rotor_chords = [rotor_chord_data[:,2] * rtip[1]]
# rotor_twists = [FM.linear(rotor_twist_data[:,1]*rtip[1], rotor_twist_data[:,2], radii[1])] * pi/180

airfoilcontour = joinpath(contour_directory, "epema_interpolated_bspline_n30_0700.dat")
airfoilcontours = fill(airfoilcontour, length(radii))
# airfoilcontours = [fill(airfoilcontour, length(radii[1]))]
index = [1]
rotor_positions = [[0.0; rotor[:"y"]; 0.0]]
rotor_orientations = [[-1.0; 0.0; 0.0]] # positive x downstream
spin_directions = [true]

epema_chords = EpemaData.wing[:"chords"]
xle = [[0.0, 0.0, 0.0]]
yle = [epema_chords[:,1]]
zle = [[0.0, 0.0, 0.0]]
wing_chord = [epema_chords[:,2]]
wing_twist = [[0.0, 0.0, 0.0]]
wing_phi = [[0.0, 0.0, 0.0]]
wing_le_sweep = 0.0
wing_ϕ = 0.0

surfacenames = ["epema wing"]

wing_npanels = 256
wing_n_chordwise_panels = 1
Res = [5e4, 1e5, 5e5, 1e6, 1e7]
Machs = [0.0, 0.1, 0.2, 0.3]

# Res_list = fill(fill(Res, length(radii[1])), length(nblades))
# Ms_list = fill(fill(Machs, length(radii[1])), length(nblades))
Res_list = fill(Res, length(nblades))
Ms_list = fill(Machs, length(nblades))

args = AS.vlm_bem_template(vinfs, plotstepi, alphas,
    xle, yle, zle, wing_chord, wing_twist, wing_phi, wing_npanels, wing_n_chordwise_panels,
    rotor_omegas, nblades, rhub, rtip, radii, rotor_chords, rotor_twists,
    airfoilcontours, index,
    rotor_positions, rotor_orientations, spin_directions, Res_list, Ms_list;
        Vref, symmetric = true, iref = 1, static_margin = 0.10, liftingline_x_over_c = 0.25,
        cambers = [fill((xc) -> 0, length(lex)) for lex in xle],
        spanwise_spacings = fill(VL.Uniform(), length(xle)),
        chordwise_spacings = fill(VL.Uniform(), length(xle)),
        wake_developement_factor = 1.0, # fully developed by default
        swirl_recovery_factor = 0.5, # as described in Veldhuis' paper
        surfacenames,
        polar_directory,
        rotor_names = ["rotor"],
        mass = 0.0,
        inertia_x = 0.0,
        inertia_y = 0.0,
        inertia_z = 0.0,
        plot_directory = joinpath(AS.topdirectory, "data","plots",TODAY),
        plot_base_name = "Epema_props_on",
        plot_extension = ".png",
        step_symbol = LS.L"\alpha ",
)

outs = AS.runsimulation!(args...)

aircraft = args[1]
parameters = args[2]
panels = aircraft.wing_system.system.surfaces[1]

span_plot = aircraft.wing_system.lifting_line_rs[1][2,2:end] / (wing[:"span"]/2) # normalized
cls_plot = parameters.cfs[1][1][3,:]
wing_chords = FM.linear(epema_chords[:,1], epema_chords[:,2], range(epema_chords[1,1], stop=epema_chords[3,1], length=length(span_plot)))
# cls_plot = cls_plot .* wing_chords / wing[:"mac"] #VortexLattice already normalizes it for now

gammas = [panel.gamma * Vref for panel in aircraft.wing_system.system.properties[1]]
velocities = [panel.velocity * Vref for panel in aircraft.wing_system.system.properties[1]]
velocities_minus_rotor = velocities .- parameters.wake_function[1].(VL.top_center.(panels)) #w/o rotor-induced velocities
# cr = CartesianIndices(panels)
Δs = VL.top_vector.(panels)

# back of the envelope approach
cls_new = 2 * gammas' ./  (Vinf * wing[:"mac"])

# including wing induced velocities
v_perp = LA.cross.(velocities_minus_rotor, Δs) ./ LA.norm.(Δs)
v_perp = LA.norm.([[v[1]; 0.0; v[3]] for v in v_perp])
cls_new2 = 2 * gammas' .* v_perp' ./ (wing[:"mac"] * Vinf^2)

# including wing and rotor induced velocities
v_perp3 = LA.cross.(velocities, Δs) ./ LA.norm.(Δs)
v_perp3 = LA.norm.([[v[1]; 0.0; v[3]] for v in v_perp3])
cls_new3 = 2 * gammas' .* v_perp3' ./ (wing[:"mac"] * Vinf^2)

epema_fig = plt.figure("epema_blown_wing_cl_epema_polar")
epema_fig.clear()
plt.plot(span_plot, cls_plot, label="VortexLattice output")
plt.plot(results[:"lift distribution VLM rotors on"][:,1], results[:"lift distribution VLM rotors on"][:,2], "--", label="Epema VLM model")
plt.scatter(results[:"normalized cl experimental"][:,1], results[:"normalized cl experimental"][:,2], label="Epema experimental")
plt.plot(span_plot, cls_new, "--", label="2 * Γ / V_{infty}")
plt.plot(span_plot, cls_new2, "--", label="2 * Γ (V_{infty} + V_w)/ V_{infty}^2")
plt.plot(span_plot, cls_new3, "--", label="2 * Γ (V_{infty} + V_w + V_r)/ V_{infty}^2")
plt.xlabel("2y/b")
plt.ylabel("normalized c_l (c/c_{mac})")
plt.legend()
if save_figs
    plt.savefig(joinpath(plot_directory,"epema_blown_wing_cl_epema_polar" * ".pdf"), bbox_inches="tight")
end

# no tests for now - eventually would like to add some tests on cl values, as well as add results to the plot

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()

app.layout = html_div() do
    html_p(id="tcid", children = "Hello Props Check"),
    dcc_graph(id="broken", animate = 3)
end

run_server(app, debug = true, dev_tools_hot_reload = false, dev_tools_props_check = false)


text:
using ChaosTools
using Test, StaticArrays
using OrdinaryDiffEq: Vern9
using DynamicalSystemsBase.Systems: lorenz

ν_thresh_lower, ν_thresh_upper = 0.1, 0.9
C_thresh_lower, C_thresh_upper = 0.1, 0.9

println("\nTesting Partially predictable chaos...")

@testset "Predictability Continuous" begin
    @testset "Lorenz map" begin
        @testset "Lorenz map - strongly chaotic" begin
            lz = lorenz(ρ=180.70)
            chaos_type, ν, C = predictability(lz; alg=Vern9(), maxiters=1e9)
            @test chaos_type == :SC
            @test ν < ν_thresh_lower
            @test C < C_thresh_lower
        end
        @testset "Lorenz map - PPC 1" begin
            lz = lorenz(ρ=180.78)
            chaos_type, ν, C = predictability(lz; alg=Vern9(), maxiters=1e9)
            @test chaos_type == :PPC
            @test ν < ν_thresh_lower
            @test C > C_thresh_upper
        end
        @testset "Lorenz map - PPC 2" begin
            lz = lorenz(ρ=180.95)
            chaos_type, ν, C = predictability(lz; alg=Vern9(), maxiters=1e9)
            @test chaos_type == :PPC
            @test ν < ν_thresh_lower
            @test C > C_thresh_upper
        end
        @testset "Lorenz map - laminar" begin
            lz = lorenz(ρ=181.10)
            chaos_type, ν, C = predictability(lz; T_max = 400, alg=Vern9(), maxiters=1e9)
            @test chaos_type == :REG
            @test ν > ν_thresh_upper
            @test C > C_thresh_upper
        end
    end
end

@testset "Predictability Discrete" begin
    @testset "Henon map" begin
        ds = Systems.henon()
        a_reg = 0.8; a_ppc = 1.11; a_reg2 = 1.0; a_cha = 1.4
        res = [:REG, :REG, :PPC, :SC]

        for (i, a) in enumerate((a_reg, a_reg2, a_ppc, a_cha))
            set_parameter!(ds, 1, a)
            chaos_type, ν, C = predictability(ds; T_max = 400000)
            @test chaos_type == res[i]
        end
    end
end


text:
### Inputs ###

"""
Supertype to any concrete input type such as [`Bosonic`](@ref), [`PartDist`](@ref), [`Distinguishable`](@ref)
and [`Undef`](@ref).
"""
abstract type InputType end

"""
Type used to notify that the input is made of FockState indistiguishable photons.
"""
struct Bosonic <: InputType
end

"""
Type used to notify that the input is made of FockState partially distinguishable
photons.
"""
abstract type PartDist <: InputType
end

"""
One parameter model of partial distinguishability interpolating between indistinguishable
and fully distinguishable photons FockState.

!!! note "Reference"
    [Sampling of partially distinguishable bosons and the relation to the
    multidimensional permanent](https://arxiv.org/pdf/1410.7687.pdf)
"""
struct OneParameterInterpolation <: PartDist
end

"""
Model of partially distinguishable photons FockState described by a randomly generated [`GramMatrix`](@ref).
"""
struct RandomGramMatrix <: PartDist
end

"""
Model of partially distinguishable photons FockState described by a provided [`GramMatrix`](@ref).
"""
struct UserDefinedGramMatrix <: PartDist
end

"""
Model of distinguishable photons FockState.
"""
struct Distinguishable <: InputType
end

"""
Model of photons FockState with undefined [`GramMatrix`](@ref).
"""
struct Undef <: InputType
end

"""
    OrthonormalBasis(vector_matrix::Union{Matrix, Nothing})

Basis of vectors ``v_1,...,v_n`` stored as columns in a ``n``-by-``r`` matrix
possibly empty.

    Fields:
        vectors_matrix::Union{Matrix,Nothing}
"""
mutable struct OrthonormalBasis

    vectors_matrix::Union{Matrix,Nothing}
    function OrthonormalBasis(vectors_matrix = nothing)
        if vectors_matrix == nothing
            new(nothing)
        else
            is_orthonormal(vectors_matrix, atol = ATOL) ? new(vectors_matrix) : error("invalid orthonormal basis")
        end
    end

end

"""
    GramMatrix{T}(n::Int) where {T<:InputType}
    GramMatrix{T}(n::Int, distinguishability_param::Real) where {T<:InputType}
    GramMatrix{T}(n::Int, S::Matrix) where {T<:InputType}

Matrix of partial distinguishability. Will automatically generate the proper
matrix related to the provided [`InputType`](@ref).

    Fields:
        - n::Int: photons number
        - S::Matrix: Gram matrix
        - rank::Union{Int, Nothing}
        - distinguishability_param::Union{Real, Nothing}
        - generating_vectors::OrthonormalBasis
"""
struct GramMatrix{T<:InputType}

    n::Int
    S::Matrix
    rank::Union{Int,Nothing}
    distinguishability_param::Union{Real,Nothing}
    generating_vectors::OrthonormalBasis

    function GramMatrix{T}(n::Int) where {T<:InputType}
        if T == Bosonic
            return new{T}(n, ones(ComplexF64,n,n), nothing, nothing, OrthonormalBasis())
        elseif T == Distinguishable
            return new{T}(n, Matrix{ComplexF64}(I,n,n), nothing, nothing, OrthonormalBasis())
        elseif T == RandomGramMatrix
            return new{T}(n, rand_gram_matrix(n), nothing, nothing, OrthonormalBasis())
        elseif T == Undef
            return new{T}(n, Matrix{ComplexF64}(undef,n,n), nothing, nothing, OrthonormalBasis())
        else
            error("type ", T, " not implemented")
        end
    end

    function GramMatrix{T}(n::Int, distinguishability_param::Real) where {T<:InputType}
        if T == OneParameterInterpolation
            return new{T}(n, gram_matrix_one_param(n,distinguishability_param), nothing, distinguishability_param, OrthonormalBasis())
        else
            T in [Bosonic,Distinguishable,RandomGramMatrix,Undef] ? error("S matrix should not be specified for type ", T) : error("type ", T, " not implemented")
        end
    end

    function GramMatrix{T}(n::Int, S::Matrix) where {T<:InputType}
        if T == UserDefinedGramMatrix
            return new{T}(n, S, nothing, nothing, OrthonormalBasis())
        else
            T in [Bosonic,Distinguishable,RandomGramMatrix,Undef] ? error("S matrix should not be specified for type ", T) : error("type ", T, " not implemented")
        end
    end

end

"""

    Input{T<:InputType}
    Input{T}(r::ModeOccupation) where {T<:InputType}
    Input{T}(r::ModeOccupation, G::GramMatrix) where {T<:InputType}

Input state at the entrance of the interferometer.

    Fields:
        - r::ModeOccupation
        - n::Int
        - m::Int: modes numbers
        - G::GramMatrix
        - distinguishability_param::Union{Real, Nothing}
"""
struct Input{T<:InputType}

    r::ModeOccupation
    n::Int
    m::Int
    G::GramMatrix
    distinguishability_param::Union{Real,Nothing}

    function Input{T}(r::ModeOccupation) where {T<:InputType}
        if T in [Bosonic, Distinguishable, Undef, RandomGramMatrix]
            return new{T}(r, r.n, r.m, GramMatrix{T}(r.n), nothing)
        else
            error("type ", T, " not implemented")
        end
    end

    function Input{T}(r::ModeOccupation, distinguishability_param::Real) where {T<:InputType}
        if T == OneParameterInterpolation
            return new{T}(r, r.n, r.m, GramMatrix{T}(r.n,distinguishability_param), distinguishability_param)
        else
            error("type ", T, " not implemented")
        end
    end

    function Input{T}(r::ModeOccupation, S::Matrix) where {T<:InputType}
        if T == UserDefinedGramMatrix
            return new{T}(r, r.n, r.m, GramMatrix{T}(r.n,S), nothing)
        else
            error("type ", T, " not implemented")
        end
    end

end

at_most_one_photon_per_bin(input::Input) = at_most_one_photon_per_bin(input.r)


text:
module YaVocabulary

import Base: size, length

export Vocabulary, get_id, get_word, size, length, one_hot_encode,
  add!, save_vocab, load_vocab

include("./tensor_dt.jl")
using ..TensorDT: Tensor

const JSON_DATA_DIR = "./vocab_data"
const D_SI = Dict{String, Integer}
const D_IS = Dict{Integer, String}
const S_N = Union{String, Nothing}
const I_N = Union{Integer, Nothing}
const T_N = Union{Tensor, Nothing}

struct Vocabulary
  word2ix::D_SI
  ix2word::D_IS

  Vocabulary() = new(D_SI(), D_IS())

  Vocabulary(w2ix::D_SI, ix2w::D_IS) = new(w2ix, ix2w)

  function Vocabulary(words::Vector{String})
    @assert length(words) > 0
    new(_add(words)...)
  end
end

##
## Public API
##

function get_id(self::Vocabulary, w::String)::I_N
  !haskey(self.word2ix, w) && (return nothing)
  self.word2ix[w]
end

function get_word(self::Vocabulary, ix::Integer)::S_N
  !(1 ≤ ix ≤ length(self.word2ix)) && (return nothing)
  self.ix2word[ix]
end

size(self::Vocabulary) = length(self.word2ix)
length(self::Vocabulary) = length(self.word2ix)

function one_hot_encode(self::Vocabulary, w::String)::T_N
  !haskey(self.word2ix, w) && (return nothing)
  w_ix = self.word2ix[w]
  [ix == w_ix ? 1. : 0. for ix ∈ 1:length(self.word2ix)]
  # FIXME: use a BitArray?
end

function add!(self::Vocabulary, word::String)
  if !haskey(self.word2ix, word)
    w_ix = length(self.word2ix) + 1
    self.word2ix[word] = w_ix
    self.ix2word[w_ix] = word
  end
end

function save_vocab(self::Vocabulary, filename::String)
  open(string(JSON_DATA_DIR, "/", filename), "w") do f
    JSON.print(f, self.word2ix)
  end
end

function load_vocab(filename::String)
  word2ix = open(string(JSON_DATA_DIR, "/", filename), "r") do f
    JSON.parse(f)
  end |> h -> D_SI(h)
  ix2word = D_IS(ix => w for (w, ix) ∈ word2ix)
  Vocabulary(word2ix, ix2word)
end

##
## Internal
##

function _add(words::Vector{String})
  w2ix = D_SI()
  ix2w = D_IS()
  for w ∈ words
    haskey(w2ix, w) && continue
    w_ix = length(w2ix) + 1
    w2ix[w] = w_ix
    ix2w[w_ix] = w
  end
  (w2ix, ix2w)
end

end ## Module

max_stars_repo_path: Data Science From Scratch/src/YaVocabulary.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
module TestGroupGenerics

using Base.Test
using RLEVectors

@testset begin

# compare group
vec = [1,1,2,2,4,4]
x = RLEVector(vec)
@test (x .< 3) == RLEVector([true,false],[4,6])

# math group
vec = [1,1,2,2,4,4]
x = RLEVector(vec)
@test x + 5 == 5 + x
@test x + 4 == RLEVector([5,5,6,6,8,8])
@test div(vec,2) == collect(div(x,2))
@test median(RLEVector([1,2,3,2,1,5,4])) == median([1,2,3,2,1,5,4])
@test median(RLEVector([1,2,3,2,1,5])) == median([1,2,3,2,1,5])
@test median(RLEVector([2, 2, 2, 3, 3, 3])) == 2.5 # median is average of end of run and next value
@test sum(RLEVector([4, 4, 5, 5, 6, 6])) == 30
@test mean(RLEVector([4, 4, 5, 5, 6, 6])) == 5.0

# findmax, findmin
@test findmin(RLEVector([1,2,3,4,1,1])) == findmin([1,2,3,4,1,1])
@test findmax(RLEVector([1,2,3,4,1,1])) == findmax([1,2,3,4,1,1])

# indexin
foo = IntegerRle( collect(1:1000), collect(5:5:5000))
x = RLEVector([2,2,4,4,3,3])
y = RLEVector([0,0,0,3,3,3,4,4])
@test indexin(x,y) == RLEVector([0,8,6],[2,4,6])
@test indexin(x,collect(3:11)) == RLEVector([0,2,1],[2,4,6])
@test indexin([200,200,1,1,5,5],foo) == [1000,1000,5,5,25,25]

# findin
@test findin(RLEVector([1,1,2,2,3,3]), RLEVector(collect(3:10))) == collect(5:6)
@test findin(RLEVector([1,1,2,2,3,3]), 3:10) == collect(5:6)
@test findin(RLEVector([1,1,2,2,3,3]), collect(3:10)) == collect(5:6)
@test findin([3,4,5],RLEVector(collect(1:4))) == [1,2]

# in
@test in(3, RLEVector( [ 1,2,2,3 ] )) == true
@test in(4, RLEVector( [ 1,2,2,3 ] )) == false
@test setdiff( Set([1,2,3,4,5]), RLEVector([1,1,2,2,4,4,5,5]) ) == [3]

end # testset

end # module