text:
let
    reset_debug_file_id!()
    @test get_debug_file_id() == 1
    @test get_debug_file_id() == 2
    reset_debug_file_id!()
end
let
    spec = ProblemSpec()
    get_distance(spec,1,1)
end
let
    spec = ProjectSpec()

    TaskGraphs.set_initial_condition!(spec,ObjectID(1),OBJECT_AT(1,1))
    @test TaskGraphs.get_initial_condition(spec,ObjectID(1)) == OBJECT_AT(1,1)
    TaskGraphs.set_initial_condition!(spec,ObjectID(2),OBJECT_AT(2,1))
    @test TaskGraphs.get_initial_condition(spec,ObjectID(2)) == OBJECT_AT(2,1)
    TaskGraphs.set_final_condition!(spec,ObjectID(1),OBJECT_AT(1,2))
    @test TaskGraphs.get_final_condition(spec,ObjectID(1)) == OBJECT_AT(1,2)
    TaskGraphs.set_final_condition!(spec,ObjectID(1),OBJECT_AT(1,3))
    @test TaskGraphs.get_final_condition(spec,ObjectID(1)) == OBJECT_AT(1,3)
    TaskGraphs.set_final_condition!(spec,ObjectID(2),OBJECT_AT(2,3))
    @test TaskGraphs.get_final_condition(spec,ObjectID(2)) == OBJECT_AT(2,3)

    # add operations to a ProjectSpec and test that the dependencies are
    # stored correctly
    op1 = construct_operation(spec,-1,[],[1],0.0)
    op2 = construct_operation(spec,-1,[1],[2],0.0)
    op3 = construct_operation(spec,-1,[2],[],0.0)
    add_operation!(spec,op1)
    add_operation!(spec,op2)
    add_operation!(spec,op3)
    @test has_edge(spec, get_operation_id(op1), ObjectID(1))
    @test has_edge(spec, ObjectID(1), get_operation_id(op2))
    @test has_edge(spec, get_operation_id(op2), ObjectID(2))
    @test has_edge(spec, ObjectID(2), get_operation_id(op3))

    # now reset the spec to empty and see if the operations can be safely added
    spec = ProjectSpec()
    add_operation!(spec,op1)
    add_operation!(spec,op2)
    add_operation!(spec,op3)
    @test has_edge(spec, get_operation_id(op1), ObjectID(1))
    @test has_edge(spec, ObjectID(1), get_operation_id(op2))
    @test has_edge(spec, get_operation_id(op2), ObjectID(2))
    @test has_edge(spec, ObjectID(2), get_operation_id(op3))

end
let
    def = SimpleProblemDef(
        project_spec = ProjectSpec(),
        r0 = [1,2,3],
        s0 = [1,2,3],
        sF = [1,2,3],
    )
    @test length(def.shapes) == length(def.s0)
end
# ProjectSpec
let
    M = 3
    object_ICs = Vector{OBJECT_AT}([
        OBJECT_AT(1,1),
        OBJECT_AT(2,2),
        OBJECT_AT(3,3)
        ])
    object_FCs = Vector{OBJECT_AT}([
        OBJECT_AT(1,4),
        OBJECT_AT(2,5),
        OBJECT_AT(3,6)
        ])
    robot_ICs = [
        ROBOT_AT(1,7),
        ROBOT_AT(2,8),
        ROBOT_AT(3,9)
        ]
    # Testing root nodes
    let
        project_spec = ProjectSpec(object_ICs,object_FCs)
        add_operation!(project_spec,construct_operation(project_spec, 3, [1,2], [3], 1.0))
        add_operation!(project_spec,construct_operation(project_spec, 6, [3], [], 0.0))
    end
    let
        project_spec = ProjectSpec(object_ICs,object_FCs)
        add_operation!(project_spec,construct_operation(project_spec, 3, [1,2], [], 1.0))
        add_operation!(project_spec,construct_operation(project_spec, 6, [3], [], 0.0))
    end
    let
        project_spec = ProjectSpec(object_ICs,object_FCs)
        add_operation!(project_spec,construct_operation(project_spec, 3, [1,2], [3], 1.0))
        add_operation!(project_spec,construct_operation(project_spec, 6, [3], [], 0.0))
    end
end
let
    solver = NBSSolver()
    prob = pctapf_problem_1(solver)
    env = get_env(prob)
    sched = get_schedule(env)
    spec = get_problem_spec(env)
    n = OBJECT_AT(1,1)
    s = generate_path_spec(sched,spec,n)
    @test s.plan_path == false
    n = ROBOT_AT(1,1)
    s = generate_path_spec(sched,spec,n)
    @test s.plan_path == true
    @test s.min_duration == 0
    @test s.free         == true
    n = GO(1,1,2)
    s = generate_path_spec(sched,spec,n)
    @test s.plan_path == true
    @test s.min_duration == 1
    @test s.tight        == true
    @test s.free         == false
    n = GO(1,1,-1)
    s = generate_path_spec(sched,spec,n)
    @test s.plan_path == true
    @test s.min_duration == 0
    @test s.tight        == true
    @test s.free         == true
    n = CARRY(1,1,1,2)
    s = generate_path_spec(sched,spec,n)
    @test s.plan_path == true
    @test s.min_duration == 1
    @test s.tight == false
    @test s.free  == false
    n = COLLECT(1,1,1)
    s = generate_path_spec(sched,spec,n)
    @test s.plan_path == true
    @test s.tight == false
    @test s.free  == false
    @test s.static == true
    n = DEPOSIT(1,1,1)
    s = generate_path_spec(sched,spec,n)
    @test s.plan_path == true
    @test s.tight == false
    @test s.free  == false
    @test s.static == true
    n = Operation(Δt=4)
    s = generate_path_spec(sched,spec,n)
    @test s.plan_path == false
    @test s.min_duration == duration(n)
    n = TEAM_ACTION{DeliveryBot,GO}(instructions=[GO(1,1,2)])
    s = generate_path_spec(sched,spec,n)
    @test s.plan_path == true
    @test s.free      == false
    @test s.static    == false
    @test s.tight     == true
    @test s.min_duration == 1
    n = TEAM_ACTION{DeliveryBot,CARRY}(instructions=[CARRY(1,1,1,2)])
    s = generate_path_spec(sched,spec,n)
    @test s.plan_path == true
    @test s.free      == false
    @test s.tight     == false
    @test s.static    == false
    @test s.min_duration == 1
    n = TEAM_ACTION{DeliveryBot,COLLECT}(instructions=[COLLECT(1,1,1)])
    s = generate_path_spec(sched,spec,n)
    @test s.plan_path == true
    @test s.free      == false
    @test s.tight     == false
    @test s.static    == true
end
let
    sched = OperatingSchedule()
    node = OBJECT_AT(1,1)
    add_node!(sched,make_node(sched,node))
    replace_in_schedule!(sched,make_node(sched,node))
    get_path_spec(sched,1)
    set_path_spec!(sched,1,get_path_spec(sched,1))
    for pred in [
            OBJECT_AT(1,2),
            ROBOT_AT(1,2),
            GO(),
            COLLECT(),
            CARRY(),
            DEPOSIT(),
            TEAM_CARRY(instructions=[CARRY()])
            ]
        generate_path_spec(sched,pred)
    end

    get_graph(sched)
    get_vtx_ids(sched)

    id = get_object_id(node)
    get_node_from_id(sched,id)
    get_vtx(sched,id)
    v = 1
    get_vtx_id(sched,v)
    get_node_from_vtx(sched,v)

    get_object_ICs(sched)
    get_robot_ICs(sched)
    get_actions(sched)
    get_operations(sched)

    zero(sched)
    LightGraphs.edges(sched)
    LightGraphs.edgetype(sched)
    LightGraphs.has_edge(sched,1,1)
    LightGraphs.has_vertex(sched,1)
    LightGraphs.inneighbors(sched,1)
    LightGraphs.outneighbors(sched,1)
    LightGraphs.is_directed(sched)
    LightGraphs.ne(sched)
    LightGraphs.nv(sched)
    LightGraphs.vertices(sched)
end
let
    # test node deletion
    sched = OperatingSchedule()
    nodes = [OBJECT_AT(1,1),OBJECT_AT(2,2),OBJECT_AT(3,3)]
    for node in nodes
        add_node!(sched,make_node(sched,node))
    end
    rem_nodes!(sched,[1,2])
    @test get_node_from_vtx(sched,1) == OBJECT_AT(3,3)
end
# Test process_schedule
let
    solver = TaskGraphsMILPSolver(AssignmentMILP())
    # project_spec, problem_spec, robot_ICs, env_graph, _ = pctapf_problem_1()
    sched, problem_spec, env_graph, _ = pctapf_problem_1()

    model = formulate_milp(solver,sched,problem_spec;cost_model=MakeSpan())
    optimize!(model)
    @test termination_status(model) == MOI.OPTIMAL
    cost = Int(round(value(objective_function(model))))
    update_project_schedule!(solver,model,sched,problem_spec,get_assignment_matrix(model))

    TaskGraphs.process_schedule!(sched)
    @test get_t0(sched,RobotID(1)) == 0
    @test get_t0(sched,RobotID(2)) == 0
    # try with perturbed start times
    set_t0!(sched,RobotID(2),1)
    TaskGraphs.process_schedule!(sched)
    @test get_t0(sched,RobotID(1)) == 0
    @test get_t0(sched,RobotID(2)) == 1

    @test length(get_vtx_ids(sched)) == nv(get_graph(sched))
    for (v,id) in enumerate(sched.vtx_ids)
        @test get_vtx(sched, id) == v
    end
end
let
    r0 = [1,4]
    s0 = [5,8,14]
    sF = [13,12,15]
    project_spec, robot_ICs = TaskGraphs.empty_pctapf_problem(r0,s0,sF)
    add_operation!(project_spec,construct_operation(project_spec,-1,[1,2],[3],0))
    add_operation!(project_spec,construct_operation(project_spec,-1,[3],  [], 0))

    filename = "/tmp/project_spec.toml"
    open(filename, "w") do io
        TOML.print(io, TOML.parse(project_spec))
    end
    project_spec_mod = read_project_spec(filename)
    # @show project_spec_mod.object_id_to_idx
    # run(`rm $filename`)

    r0 = [get_id(get_initial_location_id(robot_ICs[k])) for k in sort(collect(keys(robot_ICs)))]
    s0 = map(pred->get_id(get_initial_location_id(pred)),TaskGraphs.initial_conditions_vector(project_spec))
    sF = map(pred->get_id(get_initial_location_id(pred)),TaskGraphs.final_conditions_vector(project_spec))
    problem_def = SimpleProblemDef(project_spec,r0,s0,sF)
    filename = "/tmp/problem_def.toml"
    open(filename, "w") do io
        TOML.print(io, TOML.parse(problem_def))
    end
    problem_def = read_problem_def(filename)
    # run(`rm $filename`)
end
let
    solver = NBSSolver()
    prob = pctapf_problem_1(solver)
    model = formulate_assignment_problem(assignment_solver(solver),prob)
    sched,_ = solve_assignment_problem!(assignment_solver(solver),model,prob)
    tips = robot_tip_map(sched)
    [k=>string(get_node_from_id(sched,id)) for (k,id) in tips]
end
let
    sched = OperatingSchedule()
    r1 = add_node!(sched, make_node(sched,ROBOT_AT(1,1), ))
    r2 = add_node!(sched, make_node(sched,ROBOT_AT(2,5), ))
    n1 = add_node!(sched, make_node(sched,GO(1,1,2),     ))
    o1 = add_node!(sched, make_node(sched,OBJECT_AT(1,2),))
    n2 = add_node!(sched, make_node(sched,COLLECT(1,1,2),))
    n3 = add_node!(sched, make_node(sched,CARRY(1,1,2,3),))
    add_edge!(sched,r1,n1)
    add_edge!(sched,n1,n2)
    add_edge!(sched,o1,n2)
    add_edge!(sched,n2,n3)
    for (n,vtxs) in [
        (n3,[get_vtx(sched,n2)]),
        (n2,[get_vtx(sched,n1)]),
        (n1,[get_vtx(sched,r1)]),
        (o1,Int[]),
        ]
        v = get_vtx(sched,n)
        @test vtxs == backtrack_node(sched,v)
    end
    tips = robot_tip_map(sched)
    @test tips[r1.id] == n3.id
    @test tips[r2.id] == r2.id
    # test tips for team
    n4 = add_node!(sched, make_node(sched,GO(2,5,3)))
    add_edge!(sched,r2,n4)
    n5 = add_node!(sched, make_node(sched,GO(1,3,2)))
    add_edge!(sched,n3,n5)
    n6 = add_node!(sched,make_node(sched,TEAM_COLLECT(instructions=[COLLECT(1,1,2),COLLECT(2,1,3),]))) #,ActionID(3))
    add_edge!(sched,n4,n6)
    add_edge!(sched,n5,n6)
    tips = robot_tip_map(sched)
    @test tips[r1.id] == n6.id
    @test tips[r2.id] == n6.id
end
let
    sched = OperatingSchedule()
    for (robot_ids, node) in [
        ([1],ROBOT_AT(1,1)),
        ([1],GO(1,1,2)),
        ([],OBJECT_AT(1,2)),
        ([1,2],TEAM_CARRY(instructions=[CARRY(1,1,2,3),CARRY(2,1,3,4)]))
        ]
        n = add_node!(sched,make_node(sched,node))
        @test get_valid_robot_ids(sched,n.id) == map(i->RobotID(i),robot_ids)
    end
end


text:
using Test

using DistributedArrays
using FFTW

@testset "Distributed FFT Tests" begin

@testset "brfft Tests" begin
    rtol=1e-4
    atol=1e-5
    pids = workers()
    np = length(pids)
    D = drand(Complex{Float64}, (51, 100, 100), pids, [np,1,1])
    C = copy(Array(D))
    @test isapprox(Array(brfft(D, 100)), brfft(C,100))
end

end

ext: jl


text:
function nonsym_sixvertex_mpo(a=1., b=1., c=1.)
  d = [a 0 0 0;
       0 c b 0;
       0 b c 0;
       0 0 0 a]
  return InfiniteMPO(permute(TensorMap(complex(d), ℂ^2 ⊗ ℂ^2, ℂ^2 ⊗ ℂ^2),(1,2),(4,3)))
end

function u1_sixvertex_mpo(a=1., b=1., c=1.)
  sym = Rep[U₁](-1 => 1, 1 => 1)
  mpo =  TensorMap(zeros, ComplexF64, sym * sym, sym * sym)
  blocks(mpo)[Irrep[U₁](0)] = [b c; c b]
  blocks(mpo)[Irrep[U₁](2)] = reshape([a], (1, 1))
  blocks(mpo)[Irrep[U₁](-2)] = reshape([a], (1, 1))

  return InfiniteMPO(permute(mpo,(1,2),(4,3)))
end

function cu1_sixvertex_mpo(a=1., b=1., c=1.)
  sym = Rep[CU₁](1//2 => 1)
  mpo =  TensorMap(zeros, ComplexF64, sym * sym, sym * sym)
  blocks(mpo)[Irrep[CU₁](0, 0)] = reshape([b + c], (1, 1))
  blocks(mpo)[Irrep[CU₁](0, 1)] = reshape([-b + c], (1, 1))
  blocks(mpo)[Irrep[CU₁](1, 2)] = reshape([a], (1, 1))

  return InfiniteMPO(permute(mpo,(1,2),(4,3)))
end

max_stars_repo_name: maartenvd/MPSKitModels

max_issues_repo_name: maartenvd/MPSKitModels

max_forks_repo_name: maartenvd/MPSKitModels


text:
@testset "Interface Test" begin
  agent = Agent()
  n1 = GenericNeurons(8)
  addnode!(agent, n1; name="N1")
  n2 = GenericNeurons(8)
  addnode!(agent, n2; name="N2")

  # Hooks
  interface = AgentInterface(agent)
  addhook!(interface, "Input", interface["N1"], (:state, :I))
  addhook!(interface, "Output", interface["N2"], (:state, :F))
  delhook!(interface, "Input")
  cO = gethook(interface, "Output")
  O = get(cO)
  @test typeof(O) <: GenericNeurons
  sethook!(interface, "Output", O)
end


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

module SubscriptionClient

using Random
using Dates
using Swagger
import Swagger: field_name, property_type, hasproperty, validate_property, SwaggerApi, SwaggerModel
import Base: convert, propertynames

include("modelincludes.jl")

include("api_OperationsApi.jl")
include("api_SubscriptionsApi.jl")
include("api_TenantsApi.jl")

# export models
export convert, Location
export convert, LocationListResult
export convert, LocationMetadata
export convert, ManagedByTenant
export convert, Operation
export convert, OperationDisplay
export convert, OperationListResult
export convert, PairedRegion
export convert, Subscription
export convert, SubscriptionListResult
export convert, SubscriptionPolicies
export convert, TenantIdDescription
export convert, TenantListResult

# export operations
export convert, OperationsApi, SubscriptionsApi, TenantsApi

export check_required, field_name, property_type, hasproperty, propertynames, validate_property, convert

end


text:
@testset "Minimum Eigenvalue" begin
    @testset "Split" begin
        let t = read_test_data()
            e, n = t
            s = splitting(e, n)
            # Ensure trace ordering doesn't matter
            s′ = splitting(n, e)
            splits_are_equal(s, s′)
            @test s.phi_best ≈ 40 atol=step(s.phi)
            @test s.dt_best ≈ 1.4 atol=step(s.dt)
            @test s.spol ≈ 10 atol=1.0
            # Ensure starting orientation doesn't matter
            s′ = splitting(Seis.rotate_through(n, e, rand(1:360))...)
            splits_are_equal(s, s′)
        end
    end
    @testset "Result" begin
        let t = read_test_data(), nphi=361, ndt=80, wb=10, we=35
            e, n = t
            s = splitting(e, n, 10, 35, nphi=nphi, ndt=ndt)
            @test length(s.phi) == nphi
            @test length(s.dt) == ndt
            @test size(s.lam1) == size(s.lam2) == (nphi, ndt)
            @test s.window_start == wb
            @test s.window_end == we
            @test s.trace1 == n
            @test s.trace2 == e
        end
    end
    @testset "Window" begin
        let t = read_test_data()
            e, n = t
            s = splitting(e, n, Seis.starttime(e), Seis.endtime(e))
            s′ = splitting(e, n)
            @test s.window_start == s′.window_start
            @test s.window_end == s′.window_end
        end
    end
    @testset "Errors" begin
        let t = read_test_data()
            e, n = t
            s = splitting(e, n, e.meta.SAC_user0, e.meta.SAC_user2)
            @test s.dphi ≈ 0.75 atol=0.5
            @test s.ddt ≈ 0.013 atol=0.1
            @test s.dspol ≈ 0.35 atol=0.5
            @test s.ndf ≈ 462 atol=2
        end
    end
end

max_stars_repo_name: anowacki/SeisSplit

max_issues_repo_name: anowacki/SeisSplit

max_forks_repo_name: anowacki/SeisSplit


text:
"""
    module ConsistentRNG
Provides a an array random number generators with one for each thread.

Exports [`cRand()`](@ref) and [`seedCRNG!()`](@ref). These are thread consistent.
"""
module ConsistentRNG
#__precompile__(false) # do not precompile

using RandomNumbers
import Random
# using Distributed # only for info message

export cRand, cRandn, seedCRNG!, trng, newChildRNG # threadsafe
export sync_cRandn

"""
Baseline random number generator used throughout.
Currently we are using 'Xoshiro256StarStar' from 'RandomNumbers.jl', see the
[Documentation](https://sunoru.github.io/RandomNumbers.jl/stable/man/benchmark/#Benchmark-1)
and this [Blog post](https://nullprogram.com/blog/2017/09/21/).
In order to change the random number generator, edit 'ConsistentRNG.jl'.
"""
const CRNG = RandomNumbers.Xorshifts.Xoshiro256StarStar

"""
    CRNGs[]
Defines an array of random number generators suitable for threaded code.
For MPI or distributed
runs it should be seeded separately on each process with [`seedCRNG!`](@ref).
Currently we are using 'Xoshiro256StarStar' from 'RandomNumbers.jl', see the
[Documentation](https://sunoru.github.io/RandomNumbers.jl/stable/man/benchmark/#Benchmark-1)
and this [Blog post](https://nullprogram.com/blog/2017/09/21/).
In order to change the random number generator, edit 'ConsistentRNG.jl'.

```julia
rng = CRNGs[][Threads.threadid()]
rand(rng)
```
"""
const CRNGs = Ref([CRNG(),])
# Define this as constant with the correct type for precompilation. Will be
# initialised at runtime in __init__()

"""
    trng()
Thread local random number generator.

```julia
rand(trng())
rand(trng(),UInt)
```
"""
@inline trng() = @inbounds CRNGs[][Threads.threadid()]

"""
    r = cRand(args...)
Similar to 'rand(args)' but uses consistent random number generator 'CRNGs[]'.
'cRand()' generates a single uniformly distributed random number in the interval [0,1).
Currently we are using 'Xoshiro256StarStar' from 'RandomNumbers.jl', see the
[Documentation](https://sunoru.github.io/RandomNumbers.jl/stable/man/benchmark/#Benchmark-1)
and this [Blog post](https://nullprogram.com/blog/2017/09/21/).
"""
@inline cRand(args...) = rand(trng(), args...)
# cRand(args...) = rand(CRNGs, args...)
# #

"""
    r = cRandn(args...)
Similar to 'randn(args)' but uses consistent random number generators 'CRNGs[]'.
'cRandn()' generates a single normally distributed random number.
Currently we are using 'Xoshiro256StarStar' from 'RandomNumbers.jl', see the
[Documentation](https://sunoru.github.io/RandomNumbers.jl/stable/man/benchmark/#Benchmark-1)
and this [Blog post](https://nullprogram.com/blog/2017/09/21/).
"""
@inline cRandn(args...) = randn(trng(), args...)

"""
    seedCRNG!([seed])
Seed the threaded consistent random number generators `CRNGs[]`. If a single
number is given, this will be used to seed a random sequence, which is hashed
and then used to generate seeds for each rng in
the vector [`CRNGs`](@ref). When no argument is given, each rng is seeded
randomly.
"""
function seedCRNG!(seeds::Vector)
    @assert length(seeds) == length(CRNGs[]) "Incorrect number of seeds supplied for seeding CNRGs"
    @assert length(CRNGs[]) == Threads.nthreads() "Number of CNRGs should be equal to nthreads()"
    for (i,seed) in enumerate(seeds)
        @assert seed ≠ 0 "seed must not be zero"
        Random.seed!(CRNGs[][i], seed)
    end
    return CRNGs[]
end
function seedCRNG!(seed::Number)
    @assert length(CRNGs[]) == Threads.nthreads() "Number of CNRGs should be equal to nthreads()"
    @assert seed ≠ 0 "seed must not be zero"
    rng=CRNG(seed)
    for i = 1:length(CRNGs[])
        seedi = hash(rand(rng,UInt))
        @assert seedi ≠ 0 "seed must not be zero"
        @inbounds Random.seed!(CRNGs[][i], seedi)
    end
    return CRNGs[]
end
seedCRNG!() = map(Random.seed!, CRNGs[])

"""
    ConsistentRNG.check_crng_independence(v)
Primitve test to check the random number
generators. It throws an error if some of the threaded random number
generators are equal. Returns the number of threaded RNGs.
"""
function check_crng_independence(v)
    @assert union(ConsistentRNG.CRNGs[]) == ConsistentRNG.CRNGs[] "Threaded random number generators must be independent" ConsistentRNG.CRNGs[]
    @assert length(ConsistentRNG.CRNGs[]) == Threads.nthreads() "Number of CNRGs should be equal to nthreads()"
    return length(ConsistentRNG.CRNGs[])
end

"""
    newChildRNG(parent_rng = trng())
Random number generator that is seeded deterministically from the
thread-consistent global rng [`trng()`](@ref). By scrambling with `hash()`,
a statistically independent pseudo-random sequence from the parent rng is
accessed.
"""
function newChildRNG(parent_rng = trng())
    return CRNG(hash(rand(parent_rng,UInt)))
end

"""
    sync_cRandn(v)
Generate one random number with [`cRandn()`](@ref) in a synchronous way.
Defaults to [`cRandn()`](@ref).
"""
sync_cRandn(v) = cRandn()

# generate random seeds for the RNGs when loading the module at runtime
function __init__()
    CRNGs[] = [CRNG() for i in 1:Threads.nthreads()]
end

end # module


text:
using KernelFunctions: NeuralKernelNetwork, LinearLayer, product, Primitive

@testset "neural_kernel_network" begin
    rng, N, N′, D = MersenneTwister(123456), 5, 6, 2
    x0 = collect(range(-2.0, 2.0; length=N)) .+ 1e-3 .* randn(rng, N)
    x1 = collect(range(-1.7, 2.3; length=N)) .+ 1e-3 .* randn(rng, N)
    x2 = collect(range(-1.7, 3.3; length=N′)) .+ 1e-3 .* randn(rng, N′)

    X0 = ColVecs(randn(rng, D, N))
    X1 = ColVecs(randn(rng, D, N))
    X2 = ColVecs(randn(rng, D, N′))

    # Most of the NeuralKernelNetwork tests are currently broken.
    @testset "general test" begin

        # Specify primitives.
        k1 = 0.6 * (SEKernel() ∘ ScaleTransform(0.5))
        k2 = 0.4 * (Matern32Kernel() ∘ ScaleTransform(0.1))
        primitives = Primitive(k1, k2)

        # Build NKN Kernel.
        nkn = NeuralKernelNetwork(primitives, product)

        # Apply standard test suite.
        TestUtils.test_interface(nkn, Float64)
    end
    @testset "kernel composition test" begin
        rng = MersenneTwister(123456)

        # Specify primitives.
        k1 = rand(rng) * (SEKernel() ∘ ScaleTransform(randn(rng)))
        k2 = rand(rng) * (Matern32Kernel() ∘ ScaleTransform(randn(rng)))
        primitives = Primitive(k1, k2)

        @testset "LinearLayer" begin
            # Specify linear NKN and equivalent composite kernel.
            weights = rand(rng, 1, 2)
            nkn_add_kernel = NeuralKernelNetwork(primitives, LinearLayer(weights))
            sum_k =
                LogExpFunctions.softplus(weights[1]) * k1 +
                LogExpFunctions.softplus(weights[2]) * k2

            # Vector input.
            @test kernelmatrix_diag(nkn_add_kernel, x0) ≈ kernelmatrix_diag(sum_k, x0)
            @test kernelmatrix_diag(nkn_add_kernel, x0, x1) ≈
                kernelmatrix_diag(sum_k, x0, x1)

            # ColVecs input.
            @test kernelmatrix_diag(nkn_add_kernel, X0) ≈ kernelmatrix_diag(sum_k, X0)
            @test kernelmatrix_diag(nkn_add_kernel, X0, X1) ≈
                kernelmatrix_diag(sum_k, X0, X1)
        end
        @testset "product" begin
            nkn_prod_kernel = NeuralKernelNetwork(primitives, product)
            prod_k = k1 * k2

            # Vector input.
            @test kernelmatrix(nkn_prod_kernel, x0) ≈ kernelmatrix(prod_k, x0)
            @test kernelmatrix(nkn_prod_kernel, x0, x1) ≈ kernelmatrix(prod_k, x0, x1)

            # ColVecs input.
            @test kernelmatrix(nkn_prod_kernel, X0) ≈ kernelmatrix(prod_k, X0)
            @test kernelmatrix(nkn_prod_kernel, X0, X1) ≈ kernelmatrix(prod_k, X0, X1)
        end
    end
end


text:
function freeze_order(config, trade_grid, check_key, client)
    if trade_grid[check_key]["side"] == client.SIDE_BUY
        Pinguu.freeze_buy_order(config, trade_grid, check_key, client)
    else
        Pinguu.freeze_sell_order(config, trade_grid, check_key, client)
    end
end
export freeze_order

function freeze_buy_order(config, trade_grid, check_key, client)
    price_str, info = string(check_key), trade_grid[check_key]
    test_bot = config["testing"] == "1" ? true : false
    if test_bot
        #
        ## Just overwrite it for test bots
        value = Pinguu.get_buy_value(config, trade_grid)
        qty = floor(value / parse(Float64, price_str), digits=config["precision"])
        grid_point = Dict(
            "price"  => price_str,
            "side"   => client.SIDE_BUY,
            "status" => "wait",
            "type"   => "normal",
            "qty"    => qty,
            "value"  => value
        )
        trade_grid[check_key] = grid_point
    else
        #
        ## Check order
        order = try client.get_order(symbol=info["order"]["symbol"], orderId=info["order"]["orderId"])
        catch
            return
        end

        if parse(Float64, order["executedQty"]) == 0
            #
            ## Cancel order
            if order["status"] != client.ORDER_STATUS_CANCELED
                result = try client.cancel_order(symbol=info["order"]["symbol"], orderId=info["order"]["orderId"])
                catch
                    return
                end
            end
            value = Pinguu.get_buy_value(config, trade_grid)
            qty = floor(value / parse(Float64, price_str), digits=config["precision"])
            grid_point = Dict(
                "price"  => price_str,
                "side"   => client.SIDE_BUY,
                "status" => "wait",
                "type"   => "normal",
                "qty"    => qty,
                "value"  => value
            )
            trade_grid[check_key] = grid_point
        else
            trade_grid[check_key]["status"] = "frozen"
        end
        open(joinpath(config["bot_path"], "trade_grid.json"), "w") do f
            JSON.print(f, trade_grid, 4)
        end
    end
end
export freeze_buy_order

function freeze_sell_order(config, trade_grid, check_key, client)
    trade_grid[check_key]["status"] = "frozen"
    open(joinpath(config["bot_path"], "trade_grid.json"), "w") do f
        JSON.print(f, trade_grid, 4)
    end
end
export freeze_sell_order

function cancel_frozen_buy(config, trade_grid, check_key, client)
    price_str, info = string(check_key), trade_grid[check_key]
    test_bot = config["testing"] == "1" ? true : false
    if test_bot
        #
        ## Just overwrite it for test bots
        value = Pinguu.get_buy_value(config, trade_grid)
        qty = floor(value / parse(Float64, price_str), digits=config["precision"])
        grid_point = Dict(
            "price"  => price_str,
            "side"   => client.SIDE_BUY,
            "status" => "wait",
            "type"   => "normal",
            "qty"    => qty,
            "value"  => value
        )
        trade_grid[check_key] = grid_point
        open(joinpath(config["bot_path"], "trade_grid.json"), "w") do f
            JSON.print(f, trade_grid, 4)
        end
    else
        if parse(Float64, info["order"]["executedQty"]) == 0
            #
            ## Check order
            order = try client.get_order(symbol=info["order"]["symbol"], orderId=info["order"]["orderId"])
            catch
                return
            end

            if parse(Float64, order["executedQty"]) == 0
                #
                ## Cancel order
                if order["status"] != client.ORDER_STATUS_CANCELED
                    result = try client.cancel_order(symbol=info["order"]["symbol"], orderId=info["order"]["orderId"])
                    catch
                        return
                    end
                end
                value = Pinguu.get_buy_value(config, trade_grid)
                qty = floor(value / parse(Float64, price_str), digits=config["precision"])
                grid_point = Dict(
                    "price"  => price_str,
                    "side"   => client.SIDE_BUY,
                    "status" => "wait",
                    "type"   => "normal",
                    "qty"    => qty,
                    "value"  => value
                )
                trade_grid[check_key] = grid_point
                open(joinpath(config["bot_path"], "trade_grid.json"), "w") do f
                    JSON.print(f, trade_grid, 4)
                end
            else
                @info("Partially executed")
            end
        end
    end
end
export cancel_frozen_buy


text:
@testset "catbyrow(arr)" begin
    arr1 = [1, 2, 3]
    arr2 = [4, 5, 6]
    @test catbyrow([arr1, arr2]) == [1 2 3; 4 5 6]
end


text:
export Figure, Axis, Axis3
export hist, scatter, linesegments, mesh,
       hist!, scatter!, linesegments!, mesh!

using GLMakie: Axis, Axis3, Figure, Scatter, LineSegments
using GLMakie: Mesh as GLMakieMesh
using GLMakie: current_axis, record, hist, hist!
import GLMakie: linesegments, linesegments!, mesh, mesh!, scatter, scatter!,
                convert_arguments

const plot_nonlinear_subdivisions = 5

# geometry
include("point.jl")
include("polytope.jl")
include("axisalignedbox.jl")

# mesh
include("rectilinear_grid.jl")
include("polytope_vertex_mesh.jl")

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
#This is an example of how the SetIntersectionProjection package
#can work together with another package written in an older Julia version

#We use SetIntersectionProjection in Julia 1.1 and WAVEFORM in Julia 0.6
#The idea is as follows: 1) set up experiment in julia 1.1; 2) save velocity model
#and some other info ; 3) call a Julia 0.6 script that loads the model and runs one or a few FWI iterations and saves model again;
# 4) load saved model in Julia 1.1 and project onto intersection; 5) save model....

#This is not mathematically equivalent to running (spectral) projected gradient. It is
#intended to be a practical work-around to deal with incompatible software packages/versions.
#The results obtained using this script are not entirely the same as in the corresponding paper
#because that was generated using spectral projected gradient in julia 0.6

using LinearAlgebra
using JLD
ENV["MPLBACKEND"]="qt5agg"
using PyPlot
#using Waveform

using SetIntersectionProjection

save_dir = pwd()

# Set up the model geometry
# In 2D, we organize our coordinates as (z,x) for legacy + visualization purposes

# Number of gridpoints in z-, x-
n = [151;101];

# Spacing of gridpoints in z-, x-
d = 10.0*[1;1];

# Origin point of domain in z-, x- (typically (0.,0.))
o = 0.0*[1;1];

# Constant velocity background
vel_background = 2500;

# Source wavelet time shift (in seconds)
t0 = 0.0;

# Peak frequency of Ricker Wavelet (in Hz)
f0 = 10.0;

# Units of the model
unit = "m/s";

# Model domain range
L = n.*d;

# x-coordinates of the sources
xsrc = 0.0:50.0:L[2];xsrc=convert(Vector{Float64},xsrc)

# y-coordinates of the sources (irrelevant in 2D)
ysrc = [0.0];ysrc=convert(Vector{Float64},ysrc)

# z-coordinates of the sources
zsrc = [10.0];zsrc=convert(Vector{Float64},zsrc)

# x-coordinates of the receivers
xrec = [950.0];xrec=convert(Vector{Float64},xrec)

# y-coordinates of the receivers (irrelevant in 2D)
yrec = [0.0];yrec=convert(Vector{Float64},yrec)

# z-coordinates of the receivers (this is a transmission experiment)
#zrec = 0.0:50.0:L[2];
zrec = 0.0:50.0:L[1];zrec=convert(Vector{Float64},zrec)

# Frequencies to compute
freqs = [3.0 ; 5.0 ; 7.0]# ; 8.0]#[5.0;10.0;15.0];

# Maximum wavelength
λ = vel_background/maximum(freqs);

nsrc,nrec,nfreq = length(xsrc),length(xrec),length(freqs)

# Set up computational parameters
## Size/spacing of computational domain
comp_n = n;
comp_d = d;
comp_o = o;

# Number of PML points
npml = round(Int,λ/minimum(comp_d))*[1 1; 1 1];

# Create velocity model
v0 = vel_background*ones(n...);
v = copy(v0);
#v[div(n[1],3):2*div(n[1],3),div(n[2],3):2*div(n[2],3)] = 1.25*vel_background;
v[61:80,50:70] .= 2400
v = vec(v);
v0 = vec(v0);

function plot_velocity(vs,title_str,n,d,keyword,rec_x,rec_z,src_x,src_z)
    fig = figure(figsize=(4.2, 3.2))
    FS  = 8
    LS  = 6
    TML = 2
    PD  = 2
    #ax1 = axes()
    #vplot = imshow(reshape(vs,model.n...),vmin=minimum(v)-50,vmax=maximum(v)+50,axes=ax1,cmap="jet",extent=[0.0 , model.d[2]*model.n[2] , model.d[1]*model.n[1] , 0.0])
    vplot = imshow(reshape(vs,n[1],n[2]),vmin=minimum(v)-50,vmax=maximum(v)+50,cmap="jet",extent=[0.0 , d[2]*n[2] , d[1]*n[1] , 0.0])
    tick_params(labelsize=LS,length=TML,pad=PD)
    title(title_str,FontSize=FS)
    xlabel("x [m]",FontSize=FS)
    ylabel("z [m]",FontSize=FS)
    zt = 0:20:n[1]
    xt = 0:20:n[2]
    #xticks(xt,round.(Int,xt*model.d[2]))
    #yticks(zt,round.(Int,zt*model.d[1]))
    cbar = colorbar()
    cbar[:ax][:tick_params](labelsize=LS)
    figure;plot(rec_x,rec_z,linewidth=1.0, marker="o",linestyle="",markersize=2)
    figure;plot(src_x,src_z,linewidth=1.0, marker="x",linestyle="","k",markersize=2)
    savefig(joinpath(save_dir,string("CFWI_simple_freq_m_est_",keyword,".eps")),bbox_inches="tight",dpi=300)
    return nothing
end
plot_velocity(v,"a) True velocity",n,d,"true",repeat(xrec,length(zrec)),zrec,xsrc,repeat(zsrc,length(xsrc)))
plot_velocity(v0,"b) Initial velocity",n,d,"initial",repeat(xrec,length(zrec)),zrec,xsrc,repeat(zsrc,length(xsrc)))

# Source weight matrix
Q = Matrix{Float64}(I,nsrc, nsrc)

# Frequency continuation
max_func_evals = 2

vest = v0
#opts1 = deepcopy(opts);

bound_min = minimum(v)-50.0
bound_max = maximum(v)+50.0

v_true    = v
vel_model = deepcopy(v0)
vel_model = save("vel_model.jld","vel_model",vel_model)

#Save info for FWI
FWI_info = save("FWI_info.jld","bound_min",bound_min,"bound_max",bound_max,"nsrc",nsrc,"nfreq",nfreq,"Q",Q,"max_func_evals",max_func_evals,"v_true",v_true,"n",n,"d",d,"o",o,"t0",t0,"f0",f0,"unit",unit,"freqs",freqs,"xsrc",xsrc,"ysrc",ysrc,"zsrc",zsrc,"xrec",xrec,"yrec",yrec,"zrec",zrec,"npml",npml)

# Setup constraints:
options=PARSDMM_options()
options.FL=Float64
options=default_PARSDMM_options(options,options.FL)
options.adjust_gamma           = true
options.adjust_rho             = true
options.adjust_feasibility_rho = true
options.Blas_active            = true
options.maxit                  = 1000
options.feas_tol= 0.001
options.obj_tol=0.001
options.evol_rel_tol = 0.00001

options.rho_ini=[1.0f0]

set_zero_subnormals(true)
BLAS.set_num_threads(2)
#FFTW.set_num_threads(2)
options.parallel=false
options.feasibility_only = false
options.zero_ini_guess=true

mutable struct compgrid
  d :: Tuple
  n :: Tuple
end

comp_grid=compgrid((d[1],d[2]),(n[1],n[2]))
constraint_strategy_list=[1 2 3 4 5]#
for j in constraint_strategy_list

  #reset initial model
  vel_model = deepcopy(v0)
  vel_model = save("vel_model.jld","vel_model",vel_model)

  if j==1
    keyword="bounds_only"
    title_str="c) bounds only"

    constraint = Vector{SetIntersectionProjection.set_definitions}()

    #bounds:
    m_min     = minimum(v)-50.0
    m_max     = maximum(v)+50.0
    set_type  = "bounds"
    TD_OP     = "identity"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

elseif j==2 #various types of cardinality and rank
		keyword="cardmat_cardcol_rank_bounds"
    title_str="c)"# fiber, matrix grad. card. & rank & bounds"

    constraint = Vector{SetIntersectionProjection.set_definitions}()

    #bounds:
    m_min     = minimum(v)-50.0
    m_max     = maximum(v)+50.0
    set_type  = "bounds"
    TD_OP     = "identity"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

		#cardinality on derivatives (column and row wise)
    m_min     = 0
    m_max     = 2
    set_type  = "cardinality"
    TD_OP     = "D_x"
    app_mode  = ("fiber","x")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    m_min     = 0
    m_max     = 2
    set_type  = "cardinality"
    TD_OP     = "D_z"
    app_mode  = ("fiber","z")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

		#cardinality on derivatives (matrix based)
    m_min     = 0
    m_max     = round(Integer,3*0.33*n[2])
    set_type  = "cardinality"
    TD_OP     = "D_x"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    m_min     = 0
    m_max     = round(Integer,3*0.33*n[1])
    set_type  = "cardinality"
    TD_OP     = "D_z"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

		#rank constraint
    m_min     = 0
    m_max     = 3
    set_type  = "rank"
    TD_OP     = "identity"
    app_mode  = ("matrix","")
    custom_TD_OP = ([],false)
    push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

  elseif j==3 #true tv and bounds
      keyword="trueTV_bounds"
      title_str="d) true TV & bounds"

      constraint = Vector{SetIntersectionProjection.set_definitions}()

      #bounds:
      m_min     = minimum(v)-50.0
      m_max     = maximum(v)+50.0
      set_type  = "bounds"
      TD_OP     = "identity"
      app_mode  = ("matrix","")
      custom_TD_OP = ([],false)
      push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

      #TV
      (TV,dummy1,dummy2,dummy3)=get_TD_operator(comp_grid,"TV",options.FL)
      m_min     = 0.0
      m_max     = norm(TV*vec(v),1)
      set_type  = "l1"
      TD_OP     = "TV"
      app_mode  = ("matrix","")
      custom_TD_OP = ([],false)
      push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    elseif j==4 # various types of cardinality
    			keyword="cardmat_cardcol_bounds"
          title_str="b)"#" fiber, matrix grad. card. & bounds"

    			constraint = Vector{SetIntersectionProjection.set_definitions}()

          #bounds:
          m_min     = minimum(v)-50.0
          m_max     = maximum(v)+50.0
          set_type  = "bounds"
          TD_OP     = "identity"
          app_mode  = ("matrix","")
          custom_TD_OP = ([],false)
          push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    			#cardinality on derivatives (column and row wise)
          m_min     = 0
          m_max     = 2
          set_type  = "cardinality"
          TD_OP     = "D_x"
          app_mode  = ("fiber","x")
          custom_TD_OP = ([],false)
          push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

          m_min     = 0
          m_max     = 2
          set_type  = "cardinality"
          TD_OP     = "D_z"
          app_mode  = ("fiber","z")
          custom_TD_OP = ([],false)
          push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    			#cardinality on derivatives (matrix based)
          m_min     = 0
          m_max     = round(Integer,3*0.33*n[2])
          set_type  = "cardinality"
          TD_OP     = "D_x"
          app_mode  = ("matrix","")
          custom_TD_OP = ([],false)
          push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

          m_min     = 0
          m_max     = round(Integer,3*0.33*n[1])
          set_type  = "cardinality"
          TD_OP     = "D_z"
          app_mode  = ("matrix","")
          custom_TD_OP = ([],false)
          push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

        elseif j==5
    				keyword="cardmat_bounds"
            title_str="a)"# matrix grad. card. & bounds & rank"

    				constraint = Vector{SetIntersectionProjection.set_definitions}()

            #bounds:
            m_min     = minimum(v)-50.0
            m_max     = maximum(v)+50.0
            set_type  = "bounds"
            TD_OP     = "identity"
            app_mode  = ("matrix","")
            custom_TD_OP = ([],false)
            push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

    				#cardinality on derivatives (matrix based)
            m_min     = 0
            m_max     = round(Integer,3*0.33*n[2])
            set_type  = "cardinality"
            TD_OP     = "D_x"
            app_mode  = ("matrix","")
            custom_TD_OP = ([],false)
            push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

            m_min     = 0
            m_max     = round(Integer,3*0.33*n[1])
            set_type  = "cardinality"
            TD_OP     = "D_z"
            app_mode  = ("matrix","")
            custom_TD_OP = ([],false)
            push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

            #rank constraint
            m_min     = 0
            m_max     = 3
            set_type  = "rank"
            TD_OP     = "identity"
            app_mode  = ("matrix","")
            custom_TD_OP = ([],false)
            push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

end #end if block for constraint setup choices

  #set up constraints, precompute some things and define projector
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  (TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)
  options.rho_ini        = ones(length(TD_OP))*10.0
  # for i=1:length(options.rho_ini)
  #   if set_Prop.ncvx[i]==true
  #     options.rho_ini[i]=1.0
  #   end
  # end
  proj_intersection = x-> PARSDMM(x,AtA,TD_OP,set_Prop,P_sub,comp_grid,options)
  # function prj!{T}(dst::T, src::T)
  #   (dst[:],dummy1,dummy2,dymmy3) = proj_intersection(src)
  #   return dst
  # end
   function prj!(input)
     (x,dummy1,dummy2,dymmy3) = proj_intersection(input)
     return x
   end

  for i=1:8
    #1) run software package in different julia version for a few iterations to update model
    run(`/Applications/Julia-0.6.app/Contents/Resources/julia/bin/julia GetUPDModel.jl`)

    #2) load updated model
    tmp = load("vel_model.jld")
    vel_model = tmp["vel_model"]

    #3) project new model
    vel_model = prj!(vec(vel_model))

    #4) save projected model
    save("vel_model.jld","vel_model",vel_model)
  end

  #load final model
  tmp = load("vel_model.jld")
  vel_model = tmp["vel_model"]
  plot_velocity(vel_model,title_str,n,d,keyword,repeat(xrec,length(zrec)),zrec,xsrc,repeat(zsrc,length(xsrc)))

end #end loop over FWI with various constraint combinations

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
using WeirdMatrices
using LinearAlgebra
using BandedMatrices
using SparseArrays
using Test

@testset "WeirdMatrices.jl" begin
    @testset "OneBlockMatrices" begin
        N = 7
        upper_block = Matrix(reshape(1:16, 4,4))
        A = OneBlockMatrix(upper_block, N, N)
        denseA = zeros(N, N)
        denseA[1:4,1:4] = upper_block

        for i = 1:N, j = 1:N
            @test A[i,j] == denseA[i,j]
        end
        @test A[1:5,3:7] == denseA[1:5,3:7]

        eye = Matrix{Int}(I, size(A))
        out = similar(eye)
        mul!(out, A, eye)
        @test out == Matrix(A)
        @test eye*A == Matrix(A)

        mA = -A
        @test mA isa OneBlockMatrix
        @test mA.block == -A.block

        @testset "Resultant OneBlockMatrix" begin
            o = ones(Int, N)
            D = Diagonal(o)
            T = Tridiagonal(o[2:end], 2o, 3o[2:end])
            sT = SymTridiagonal(1o, 2o[2:end])
            S = SkewTridiagonal(1:N-1)
            B = BandedMatrix(2 => 1:N-2, 3=>1:N-3)
            Bu = Bidiagonal(1:N, 1:N-1, 'U')
            Bl = Bidiagonal(1:N, 1:N-1, 'L')

            function case(A,ref)
                @test A isa OneBlockMatrix
                @test A == ref
            end

            case(A, A*D)
            case(A, D*A)

            case(T*A, Matrix(T)*Matrix(A))
            @test Matrix(T)*A == Matrix(T)*Matrix(A)

            case(A*T, Matrix(A)*Matrix(T))
            @test A*Matrix(T) == Matrix(A)*Matrix(T)

            case(T*A*T, Matrix(T)*Matrix(A)*Matrix(T))

            case(A*(T*A), Matrix(A)*Matrix(T)*Matrix(A))
            case((A*T)*A, Matrix(A)*Matrix(T)*Matrix(A))

            case((A*T)*T, Matrix(A)*Matrix(T)*Matrix(T))

            case((A*T)*A, Matrix(A)*Matrix(T)*Matrix(A))

            case(T*(T*A), Matrix(T)*Matrix(T)*Matrix(A))

            case(sT*A, Matrix(sT)*Matrix(A))
            case(A*sT, Matrix(A)*Matrix(sT))
            case((sT*A)*sT, Matrix(sT)*Matrix(A)*Matrix(sT))
            case((A*sT)*sT, Matrix(A)*Matrix(sT)*Matrix(sT))
            case((sT*A)*T, Matrix(sT)*Matrix(A)*Matrix(T))
            case((A*sT)*T, Matrix(A)*Matrix(sT)*Matrix(T))

            case(S*A, Matrix(S)*Matrix(A))

            case(B*A, Matrix(B)*Matrix(A))
            case(B*B*A, Matrix(B)*Matrix(B)*Matrix(A))
            case(B*A*B, Matrix(B)*Matrix(A)*Matrix(B))

            case(Bl*A, Matrix(Bl)*Matrix(A))
            case(Bl*(Bl*A), Matrix(Bl)*Matrix(Bl)*Matrix(A))
            case(Bl*A*Bl, Matrix(Bl)*Matrix(A)*Matrix(Bl))

            case(Bu*A, Matrix(Bu)*Matrix(A))
            case(Bu*(Bu*A), Matrix(Bu)*Matrix(Bu)*Matrix(A))
            case(Bu*A*Bu, Matrix(Bu)*Matrix(A)*Matrix(Bu))
        end

        @testset "Multiplication with arbitrary arrays" begin
            M = reshape(1.0:4N, N, 4)
            out = similar(M)
            mul!(out, A, M)
            @test out == Matrix(A)*M

            v = M[:,1]
            w = similar(v)
            mul!(w, A, v)
            @test w == Matrix(A)*v
        end
    end

    @testset "ColumnSparseMatrix" begin
        m,n = 10, 150
        u = rand(ComplexF64, m, 4)

        is = [2,4,10,6]
        A = ColumnSparseMatrix(m, n, u, is)
        @test issorted(A.col_indices)
        MA = Matrix(A)
        B = rand(ComplexF64, n, m)

        @test A*B ≈ MA*B

        C = rand(m, n)
        @test dot(A, C) ≈ dot(MA, C)
        @test dot(C, A) ≈ dot(C, MA)
        # display(@benchmark(dot(A,C)))
        # display(@benchmark(dot(MA,C)))

        @test dot(A,A) ≈ dot(MA,MA)
        # # Not as dramatic a difference as e.g. dot(A,C) vs dot(MA,C)
        # display(@benchmark(dot(A,A)))
        # display(@benchmark(dot(MA,MA)))

        CMA = ColumnSparseMatrix(MA)
        @test CMA isa ColumnSparseMatrix
        @test CMA == A

        BA = B*A
        @test BA isa ColumnSparseMatrix
        @test BA ≈ B*MA

        # display(@benchmark(B*A))
        # display(@benchmark(B*MA))AA = A'A

        AA = A'A
        @test AA isa SparseMatrixCSC
        @test AA ≈ MA'MA

        Ahalf = A/2
        @test Ahalf isa ColumnSparseMatrix
        @test Ahalf ≈ MA/2

        Atwice = A*2
        @test Atwice isa ColumnSparseMatrix
        @test Atwice ≈ MA*2
    end
end


text:
export Relaxed, Acquire, AcquRel, Release

const argtype_map = Dict(
    :hsa_signal_t => :Signal,
    :hsa_agent_t => :Agent,
    :(Ptr{hsa_queue_t}) => :Queue,
    :hsa_region_t => :Region,
)

# Memory Consistency Enumeration
const Relaxed = 0
const Acquire = 1
const AcquRel = 2
const Release = 3


text:
type NIGN
    N::Int
    sum_x::Float64
    sum_x_sq::Float64
    function NIGN()
        new(0, 0.0, 0.0)
    end
end

function incorporate!(comp::NIGN, x::Float64)
    comp.N += 1
    comp.sum_x += x
    comp.sum_x_sq += x*x
end

function unincorporate!(comp::NIGN, x::Float64)
    @assert comp.N > 0
    comp.N -= 1
    comp.sum_x -= x
    comp.sum_x_sq -= x*x
end

immutable NIGNParamsPrior
    m_prior::Normal
    r_prior::Gamma
    nu_prior::Gamma
    s_prior::Gamma
end

type NIGNParams
	m::Float64
	r::Float64
    nu::Float64
    s::Float64
    function NIGNParams(m::Float64, r::Float64, nu::Float64, s::Float64)
        new(m, r, nu, s)
    end
    function NIGNParams(prior::NIGNParamsPrior)
        m = rand(prior.m_prior)
        r = rand(prior.r_prior)
        nu = rand(prior.nu_prior)
        s = rand(prior.s_prior)
        new(m, r, nu, s)
    end
end

function log_probability_density(prior::NIGNParamsPrior, params::NIGNParams)
    l = logpdf(prior.m_prior, params.m)
    l += logpdf(prior.r_prior, params.r)
    l += logpdf(prior.nu_prior, params.nu)
    l += logpdf(prior.s_prior, params.s)
    l
end

function log_z(r::Float64, s::Float64, nu::Float64)
    lz = ((nu + 1.) / 2.) * log(2)
    lz += .5 * log(pi)
    lz -= .5 * log(r)
    lz -= (nu/2.) * log(s)
	lz += lgamma(nu/2.0)
    lz
end

function posterior_params(comp::NIGN, params::NIGNParams)
	rn = params.r + Float64(comp.N)
	nun = params.nu + Float64(comp.N)
	mn = (params.r*params.m + comp.sum_x)/rn
	sn = params.s + comp.sum_x_sq + params.r*params.m*params.m - rn*mn*mn
	if sn == Float64(0)
		sn = params.s
	end
	(mn, rn, sn, nun)
end

function log_probability_density(comp::NIGN, params::NIGNParams)
    (mn, rn, sn, nun) = posterior_params(comp, params)
    Z0 = log_z(params.r, params.s, params.nu)
    ZN = log_z(rn, sn, nun)
	-(comp.N/2.) * log(2 * pi) + ZN - Z0
end


text:
using LibSerialPort
using Test

function test_blocking_serial_loopback(sp::SerialPort)

    # This "tests" `read` indirectly, but has an important side effect of flushing
    # the input buffer, putting the test into a more well-defined state.
    println("\nRead all data currently waiting in the serial input buffer...")
    if bytesavailable(sp) > 0
        println(String(read(sp)))
    end

    for i = 1:10
        msg = "test message $i\n"
        write(sp, msg)
        print("Wrote ", msg)
        sleep(0.1)
        line = readline(sp)
        println(line)
        @test startswith(line, "test message")
    end

end

function test_nonblocking_serial_loopback(sp::SerialPort)

    println("\nnonblocking_read...")
    for i = 1:100
        print(String(nonblocking_read(sp)))
        sleep(0.001)
    end
    println()

    write(sp, "test message.\n")

    sleep(0.1)

    result = String(nonblocking_read(sp))
    show(result)
    println()
    @test occursin("test message", result)
end

"""
Check that a serial device returns the expected result under a given
setting for the number of data bits `n` in a UART frame.

The most common value for n is 8. Support for other values is device-dependent,
but only 5,6,7,8 are supported by libserialport.
"""
function test_bytestring_roundtrip(s::SerialPort, n::Int=8)

    # If 0x11 and 0x13 are missing, flow control may be enabled.
    # set_flow_control(s, xonxoff=SP_XONXOFF_DISABLED)

    msg = collect(UInt8(0):UInt8(255))

    write(s, msg)
    sp_drain(s)

    sleep(0.1)
    # rcv = nonblocking_read(s)
    rcv = read(s)

    # Expected response
    msg .&= UInt8(1 << n - 1)

    if rcv != msg
        @warn("Received data does not match transmitted data:")
        println("Bytes received: ", rcv)
        println("Bytes missing:  ", setdiff(msg, rcv));
        println("Bytes added:    ", setdiff(rcv, msg));
    end

    @test rcv == msg
end

function test_high_level_api(args...)

    if length(args) != 2
        println("Usage: $(basename(@__FILE__)) port baudrate")
        println("Available ports:")
        list_ports()
        return
    end

    if !ispath(args[1])
        println("Not found: ", args[1])
        return
    end

    sp = open(args[1], parse(Int, args[2]))

    sp.read_timeout_ms = 1000

    print_port_metadata(sp)
    print_port_settings(sp)

    sleep(2)  # MCU serial startup

    @testset "Blocking read/write" begin
        test_blocking_serial_loopback(sp)
    end

    @testset "Nonblocking read" begin
        test_nonblocking_serial_loopback(sp)
    end

    @testset "Bytestring roundtrip" begin
        test_bytestring_roundtrip(sp)
    end

    close(sp)
    return
end

test_high_level_api(ARGS...)


text:
# automatically generated by the FlatBuffers compiler, do not modify

MyGame.Example.eval(quote

FlatBuffers.@with_kw mutable struct Referrable
    id::UInt64 = 0
end
FlatBuffers.@ALIGN(Referrable, 1)
FlatBuffers.slot_offsets(::Type{T}) where {T<:Referrable} = [
    0x00000004
]

Referrable(buf::AbstractVector{UInt8}) = FlatBuffers.read(Referrable, buf)
Referrable(io::IO) = FlatBuffers.deserialize(io, Referrable)

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
"""
    struct CodeGenerator <: AbstractASTConsumer
Hold a pointer to a `clang::CodeGenerator` object.
"""
struct CodeGenerator <: AbstractASTConsumer
    ptr::CXCodeGenerator
end

Base.unsafe_convert(::Type{CXCodeGenerator}, x::CodeGenerator) = x.ptr
Base.cconvert(::Type{CXCodeGenerator}, x::CodeGenerator) = x


text:
module BesselFunctions

using AndExport
# using GSL
# using SpecialFunctions

# @xport sphj(n, x) = GSL.sf_bessel_jl(n, x)
# @xport sphn(n, x) = GSL.sf_bessel_yl(n, x)
# @xport sphh₊(n, x) = sphj(n, x) + im*sphn(n, x)
# @xport sphh₋(n, x) = sphj(n, x) - im*sphn(n, x)
# # @xport sphj(n, x) = sqrt(π/2x) * besselj(n+0.5, x)
# # @xport sphn(n, x) = sqrt(π/2x) * bessely(n+0.5, x)
# # @xport sphh₊(n, x) = sqrt(π/2x) * besselh1(n+0.5, x)
# # @xport sphh₋(n, x) = sqrt(π/2x) * besselh2(n+0.5, x)

# * Changed to use the complex versions in SpecialFunctions
using SpecialFunctions

@xport sphj(n,x) = sqrt(π/2x) * besselj(n+1/2, x)
@xport sphn(n,x) = sqrt(π/2x) * bessely(n+1/2, x)
@xport sphh₊(n, x) = sqrt(π/2x) * hankelh1(n+1/2, x)
@xport sphh₋(n, x) = sqrt(π/2x) * hankelh2(n+1/2, x)

#####

@xport dsphjdx(n, x) = n/x * sphj(n, x) - sphj(n+1,x)
@xport dsphndx(n, x) = n/x * sphn(n, x) - sphn(n+1,x)
@xport dsphh₊dx(n, x) = n/x * sphh₊(n, x) - sphh₊(n+1,x)
@xport dsphh₋dx(n, x) = n/x * sphh₋(n, x) - sphh₋(n+1,x)

@xport dsphjdr(n,k,r) = dsphjdx(n,k*r) * k
@xport dsphndr(n,k,r) = dsphndx(n,k*r) * k

@xport dsphh₊dr(n,k,r) = dsphh₊dx(n,k*r) * k
@xport dsphh₋dr(n,k,r) = dsphh₋dx(n,k*r) * k

@xport ricj(n,x) = x*sphj(n,x)
# Note: this -ve makes ricn(0,x) = cos(x)
@xport ricn(n,x) = -x*sphn(n,x)
# Note: these negatives are to recover the proper hankel behaviour
# (which is to obtain i^(-n-1) exp(iz) and i^(n+1) exp(-iz)
@xport rich₊(n,x) = ricj(n,x) - im*ricn(n,x)
@xport rich₋(n,x) = ricj(n,x) + im*ricn(n,x)

@xport dricjdx(n,x) = sphj(n,x) + x*dsphjdx(n,x)
@xport dricndx(n,x) = -sphn(n,x) + -x*dsphndx(n,x)
@xport drich₊dx(n,x) = dricjdx(n,x) - im*dricndx(n,x)
@xport drich₋dx(n,x) = dricjdx(n,x) + im*dricndx(n,x)

@xport dricjdr(n,k,r) = k * dricjdx(n,k*r)
@xport dricndr(n,k,r) = k * dricndx(n,k*r)
@xport drich₊dr(n,k,r) = k * drich₊dx(n,k*r)
@xport drich₋dr(n,k,r) = k * drich₋dx(n,k*r)

end


text:
"""
    getting(f) :: Lens

Apply a callable `f` (typically a type constructor) before getting the
value; i.e.,

```
get(obj, lens ∘ getting(f)) == f(get(obj, lens))
```

This is useful for, e.g., getting a tuple as a `StaticVector` and
converting it back to a tuple when setting.

Note that `getting` requires some properties for `f` and the values
stored in the "field."  See the details below.

# Examples
```jldoctest
julia> using Kaleido, Setfield, StaticArrays

julia> obj = (x = ((0, 1, 2), "A"), y = "B");

julia> lens = (@lens _.x[1]) ∘ getting(SVector)
(@lens _.x[1]) ∘ (←|SArray{Tuple{S},T,1,S} where T where S→)

julia> get(obj, lens) === SVector(obj.x[1])
true

julia> set(obj, lens, SVector(3, 4, 5))
(x = ((3, 4, 5), "A"), y = "B")
```

```jldoctest; filter = r"#[0-9]+"
julia> using Kaleido, Setfield, StaticArrays

julia> obj = (x = ((a = 0, b = 1, c = 2), "A"), y = "B");

julia> lens = (@lens _.x[1]) ∘ getting(Base.splat(SVector))
(@lens _.x[1]) ∘ (←|#60→)

julia> get(obj, lens) === SVector(obj.x[1]...)
true

julia> set(obj, lens, SVector(3, 4, 5))
(x = ((a = 3, b = 4, c = 5), "A"), y = "B")
```

# Details

The lens created by `getting(f)` relies on that:

* The output value `y = f(x)` can be converted back to the original
  value `x` by `C(y)` where `C` is a constructor of `x`; i.e., for any
  `x` that could be retrieved from the object through this lens,

  ```
  C(f(x)) == x
  ```

* The conversion in the reverse direction also holds; i.e., for any
  `y` that could be stored into the object through this lens,

  ```
  f(C(y)) == y
  ```

The constructor `C` can be controlled by defining
`Setfield.constructor_of` for custom types of `x`.
"""
getting(f) = ConverterLens(prefer_singleton_callable(f))

struct ConverterLens{T} <: KaleidoLens
    f::T
end

Setfield.get(obj, lens::ConverterLens) = lens.f(obj)
Setfield.set(obj, ::ConverterLens, value) = _constructor_of(typeof(obj))(value)


text:
# Xmap
# =================================

struct Xmap{F<:Transform, Tf<:Number, Ti<:Number, d} <: MapField{F,Tf,Ti,d}
    ft::F
    fd::Array{Tf,d}
    function Xmap{F,Tf,Ti,d}(ft::F, fd::AbstractArray{T,d})  where {Tf,Ti,T,d,F<:Transform{Tf,d}}
        @assert Ti == eltype_out(ft)
        @assert size(fd) == size_in(ft)
        new{F,Tf,Ti,d}(ft, fd)
    end
end

# Xfourier
# =================================

struct Xfourier{F<:Transform, Tf<:Number, Ti<:Number, d}  <: FourierField{F,Tf,Ti,d}
    ft::F
    fd::Array{Ti,d}
    function Xfourier{F,Tf,Ti,d}(ft::F, fd::AbstractArray{T,d})  where {Tf,Ti,T,d,F<:Transform{Tf,d}}
        @assert Ti == eltype_out(ft)
        @assert size(fd) == size_out(ft)
        new{F,Tf,Ti,d}(ft, fd)
    end
end

# Interface definitions required to hook into Field definitions
# =================================

FourierField(::Type{Xfourier{F,Tf,Ti,d}}) where {Tf,Ti,d, F<:Transform{Tf,d}} = Xfourier{F,Tf,Ti,d}
FourierField(::Type{Xmap{F,Tf,Ti,d}}) where {Tf,Ti,d, F<:Transform{Tf,d}}     = Xfourier{F,Tf,Ti,d}
MapField(::Type{Xmap{F,Tf,Ti,d}})     where {Tf,Ti,d, F<:Transform{Tf,d}}     = Xmap{F,Tf,Ti,d}
MapField(::Type{Xfourier{F,Tf,Ti,d}}) where {Tf,Ti,d, F<:Transform{Tf,d}}     = Xmap{F,Tf,Ti,d}

Xfield{F,Tf,Ti,d} = Union{Xfourier{F,Tf,Ti,d}, Xmap{F,Tf,Ti,d}}
@inline fieldtransform(f::Xfield{F,Tf,Ti,d})  where {Tf,Ti,d, F<:Transform{Tf,d}} = f.ft
@inline fielddata(f::Xfield{F,Tf,Ti,d}) where {Tf,Ti,d, F<:Transform{Tf,d}} = f.fd

# Extras
# =================================

function Xfourier(ft::F, fd::AbstractArray{T,d})  where {T,Tf,d,F<:Transform{Tf,d}}
    @assert size(fd) == size_out(ft)
    Ti = eltype_out(ft)
    Xfourier{F,Tf,Ti,d}(ft, fd)
end

function Xmap(ft::F, fd::AbstractArray{T,d})  where {T,Tf,d,F<:Transform{Tf,d}}
    @assert size(fd) == size_in(ft)
    Ti = eltype_out(ft)
    Xmap{F,Tf,Ti,d}(ft, fd)
end

Xmap(f::Xfield{F,Tf,Ti,d})     where {Tf,Ti,d, F<:Transform{Tf,d}} = MapField(f)
Xfourier(f::Xfield{F,Tf,Ti,d}) where {Tf,Ti,d, F<:Transform{Tf,d}} = FourierField(f)

# zero and constant constructors
Xmap(ft::F)                 where {Tf,d,F<:Transform{Tf,d}} = Xmap(ft, zeros(Tf, size_in(ft)))
Xmap(ft::F, n::Number)      where {Tf,d,F<:Transform{Tf,d}} = Xmap(ft, fill(Tf(n), size_in(ft)))
function Xfourier(ft::F) where {Tf,d,F<:Transform{Tf,d}}
    Ti = eltype_out(ft)
    Xfourier(ft, zeros(Ti, size_out(ft)))
end
function Xfourier(ft::F, n::Number) where {Tf,d,F<:Transform{Tf,d}}
    Ti = eltype_out(ft)
    Xfourier(ft, fill(Ti(n), size_out(ft)))
end


text:
using DesignAutomation

include("./functions.jl")

@testset "Random Hill Climbing" begin
    x0 = [6.0, -2.0, -2.0, 6.0]
    res = random_hill(levy3, x0; max_failed=16, verbose=true)
    res = random_hill(levy3, x0; h=.01, max_failed=16, verbose=true)

    x0 = [4.2, -1.2, 2.4, -3.8]
    res = random_hill(levy3, x0; max_failed=16, verbose=true)
    res = random_hill(levy3, x0; h=.01, max_failed=16, verbose=true)

    x0 = [0.0, 0.0, 0.0, 3.2]
    res = random_hill(levy3, x0; max_failed=16, verbose=true)
    res = random_hill(levy3, x0; h=.01, max_failed=16, verbose=true)

    x0 = [1.0, 1.0, 1.0, 1.0]
    res = random_hill(levy3, x0; max_failed=16, verbose=true)
    res = random_hill(levy3, x0; h=.01, max_failed=16, verbose=true)
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
# Affine coupling layer from Dinh et al. (2017)
# Includes 1x1 convolution from in Putzky and Welling (2019)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export LearnedCouplingLayerSLIM

"""
    CS = LearnedCouplingLayerSLIM(nx1, nx2, nx_in, ny1, ny2, ny_in, n_hidden, batchsize;
        logdet::Bool=false, permute::Bool=false, k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)

 Create an invertible SLIM coupling layer with a learned data-to-image-space map.

 *Input*:

 - `nx1`, `nx2`, `nx_in`: spatial dimensions and no. of channels of input image

 - `ny1`, `ny2`, `ny_in`: spatial dimensions and no. of channels of input data

 - `n_hidden`: number of hidden units in conditional residual block

 - `loget`: bool to indicate whether to return the logdet (default is `false`)

 - `permute`: bool to indicate whether to apply a channel permutation (default is `false`)

 - `k1`, `k2`: kernel size of convolutions in residual block. `k1` is the kernel of the first and third
    operator, `k2` is the kernel size of the second operator

 - `p1`, `p2`: padding for the first and third convolution (`p1`) and the second convolution (`p2`)

 - `s1`, `s2`: stride for the first and third convolution (`s1`) and the second convolution (`s2`)

 *Output*:

 - `CS`: Invertible SLIM coupling layer with learned data-to-image map

 *Usage:*

 - Forward mode: `Y, logdet = CS.forward(X, D, A)`    (if constructed with `logdet=true`)

 - Inverse mode: `X = CS.inverse(Y, D, A)`

 - Backward mode: `ΔX, X = CS.backward(ΔY, Y, D, A)`

 - where `A` is a linear forward modeling operator and `D` is the observed data.

 *Trainable parameters:*

 - None in `CL` itself

 - Trainable parameters in residual block `CL.RB` and 1x1 convolution layer `CL.C`

 See also: [`Conv1x1`](@ref), [`ResidualBlock`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct LearnedCouplingLayerSLIM <: NeuralNetLayer
    C::Union{Conv1x1, Nothing}
    RB::ConditionalResidualBlock
    logdet::Bool
end

@Flux.functor LearnedCouplingLayerSLIM

# Constructor from input dimensions
function LearnedCouplingLayerSLIM(nx1, nx2, nx_in, ny1, ny2, ny_in, n_hidden, batchsize;
    k1=3, k2=3, p1=1, p2=1, s1=1, s2=1, logdet::Bool=false, permute::Bool=false)

    # 1x1 Convolution and residual block for invertible layer
    permute == true ? (C = Conv1x1(nx_in)) : (C = nothing)
    RB = ConditionalResidualBlock(nx1, nx2, Int(nx_in/2), ny1, ny2, ny_in, n_hidden, batchsize; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2)

    return LearnedCouplingLayerSLIM(C, RB, logdet)
end

# Forward pass: Input X, Output Y
function forward(X::AbstractArray{Float32, 4}, D, CS::LearnedCouplingLayerSLIM)

    # Get dimensions
    nx, ny, n_s, batchsize = size(X)

    # Permute and split
    isnothing(CS.C) ? (X_ = copy(X)) : (X_ = CS.C.forward(X))
    X1_, X2_ = tensor_split(X_)

    # Coupling layer
    Y1_ = copy(X1_)
    Y2_ = X2_ + CS.RB.forward(X1_, D)[1]
    Y_ = tensor_cat(Y1_, Y2_)

    isnothing(CS.C) ? (Y = copy(Y_)) : (Y = CS.C.inverse(Y_))
    CS.logdet == true ? (return Y, 0f0) : (return Y)
end

# Inverse pass: Input Y, Output X
function inverse(Y::AbstractArray{Float32, 4}, D, CS::LearnedCouplingLayerSLIM; save=false)

    # Get dimensions
    nx, ny, n_s, batchsize = size(Y)

    # Permute and split
    isnothing(CS.C) ?  (Y_ = copy(Y)) : (Y_ = CS.C.forward(Y))
    Y1_, Y2_ = tensor_split(Y_)

    # Coupling layer
    X1_ = copy(Y1_)
    X2_ = Y2_ - CS.RB.forward(X1_, D)[1]
    X_ = tensor_cat(X1_, X2_)

    isnothing(CS.C) ? (X = copy(X_)) : (X = CS.C.inverse(X_))
    save == true ? (return X, X_) : (return X)
end

# Backward pass: Input (ΔY, Y), Output (ΔX, X)
function backward(ΔY::AbstractArray{Float32, 4}, Y::AbstractArray{Float32, 4}, D, CS::LearnedCouplingLayerSLIM)

    # Recompute forward states
    X, X_ = inverse(Y, D, CS; save=true)
    nx1, nx2, nx_in, batchsize = size(X)

    # Backpropagation
    isnothing(CS.C) ? (ΔY_ = copy(ΔY)) : (ΔY_ = CS.C.forward((ΔY, Y))[1])
    ΔY1_, ΔY2_ = tensor_split(ΔY_)
    ΔX2_ = copy(ΔY2_)
    ΔX1_, ΔD = CS.RB.backward(ΔY2_, D.*0f0, tensor_split(X_)[1], D)[1:2]
    ΔX1_ += ΔY1_

    ΔX_ = tensor_cat(ΔX1_, ΔX2_)
    isnothing(CS.C) ? (ΔX = copy(ΔX_)) : (ΔX = CS.C.inverse((ΔX_, X_))[1])

    return ΔX, ΔD, X
end

# Clear gradients
function clear_grad!(CS::LearnedCouplingLayerSLIM)
    ~isnothing(CS.C) && clear_grad!(CS.C)
    clear_grad!(CS.RB)
end

# Get parameters
function get_params(CS::LearnedCouplingLayerSLIM)
    p = get_params(CS.RB)
    ~isnothing(CS.C) && (p = cat(p, get_params(CS.C); dims=1))
    return p
end

# Put parameters
function put_params!(CS::LearnedCouplingLayerSLIM, Params::Array{Any,1})
    put_params!(CS.RB, Params[1:5])
    ~isnothing(CS.C) && put_params!(CS.C, Params[6:end])
end


text:
module test_suite_applications

using DomainSets, BasisFunctions, FrameFun, Plots, StaticArrays
# Plots loads the default backend (PyPlot unless otherwise specified)

using Test

FE = FrameFun
BA = BasisFunctions

function delimit(s::AbstractString)
    println()
    println("############")
    println("# ",s)
    println("############")
end

delimit("Plotting")

function test_plots()
    gr()
    # test plotting functionsets
    B = ChebyshevT(5)
    plot(B)

    plot(B[4])

    F = Fourier(5)
    plot(F[2],plot_complex=true)

    B = (ChebyshevT(5) → -2..0) ⊕ (Fourier(5) → 0..2)
    plot(B[2:8],plot_complex=true)

    F = Fourier(10)⊗Fourier(10)
    plot(F[82],plot_complex=false)

    G = interpolation_grid(ChebyshevT(51))
    plot(G)

    G = interpolation_grid((ChebyshevT(51) → -1..0.3) ⊗ (ChebyshevT(51) → -0.5..0.5))
    G = FrameFun.subgrid(G,mandelbrot())
    plot(G)

    B = (Fourier(21) → -1..1)^3
    D = Ball()\Rectangle(SA[-0.5,-0.5,-2], SA[0.5,0.5,2])
    Df = ExtensionFrame(D,B)
    G = sampling_grid(Df)
    plot(G,size=(400,400))

    basis = Fourier(51) → -2.0..(-0.5)
    domain = -1.7..1.0
    f(x) = cos(3*x.^2)
    F = Fun(f, basis, domain)
    # Easily combine multiple plots
    plot(sampling_grid(dictionary(F)),label="grid",markercolor=:white)
    plot!(F,label="function")
    plot!(F',title="Function and derivative",label="derivative",legend=true)

    plot(F,f,label="function")
    df(x) = -sin(3*x^2)*6*x
    plot!(F',df,label="derivative",legend=true)

    basis = (Fourier(21) → -1..1)^2
    domain = intersect(Disk(0.8),Disk(0.4))
    f(x,y) = cos(7*x-2*y^2)
    F = Fun(f, basis, domain)

    plot(domain, n=100)

    plot(F,size=(600,400))

    contourf(F,colorbar=true)

    plot(F)

    P = contour(F,f)
end

test_plots()

end


text:
@static if !isdefined(Base, Symbol("@info"))
    macro info(msg)
        return :(info($(esc(msg))))
    end
end

@static if VERSION >= v"0.7.0-DEV.3656"
    using Pkg.API: installed
end

function check_grdir()
    if "GRDIR" in keys(ENV)
        grdir = ENV["GRDIR"]
        have_dir = length(grdir) > 0
        if have_dir
            have_dir = isdir(joinpath(grdir, "fonts"))
        end
    else
        have_dir = false
        for d in (homedir(), "/opt", "/usr/local", "/usr")
            grdir = joinpath(d, "gr")
            if isdir(joinpath(grdir, "fonts"))
                have_dir = true
                break
            end
        end
    end
    if have_dir
        @info "Found exisitng GR run-time in $grdir"
    end
    have_dir
end

function get_version()
    version = v"0.39.0"
    try
@static if VERSION >= v"0.7.0-DEV.3656"
        v = installed()["GR"]
else
        v = Pkg.installed("GR")
end
        if string(v)[end:end] == "+"
            version = "latest"
        end
    catch
    end
    version
end

function get_os_release(key)
    value = try
        String(read(pipeline(`cat /etc/os-release`, `grep ^$key=`, `cut -d= -f2`)))[1:end-1]
    catch
        ""
    end
    if VERSION < v"0.7-"
        replace(value, "\"", "")
    else
        replace(value, "\"" => "")
    end
end

if !check_grdir()
  if Sys.KERNEL == :NT
    os = :Windows
  else
    os = Sys.KERNEL
  end
  arch = Sys.ARCH
  if os == :Linux && arch == :x86_64
    if isfile("/etc/redhat-release")
      rel = String(read(pipeline(`cat /etc/redhat-release`, `sed s/.\*release\ //`, `sed s/\ .\*//`)))[1:end-1]
      if rel > "7.0"
        os = "Redhat"
      end
    elseif isfile("/etc/os-release")
      id = get_os_release("ID")
      id_like = get_os_release("ID_LIKE")
      if id == "ubuntu" || id_like == "ubuntu"
        os = "Ubuntu"
      elseif id == "debian" || id_like == "debian"
        os = "Debian"
      end
    end
  elseif os == :Linux && arch in [:i386, :i686]
    arch = :i386
  end
  version = get_version()
  tarball = "gr-$version-$os-$arch.tar.gz"
  rm("downloads", force=true, recursive=true)
  @info("Downloading pre-compiled GR $version $os binary")
  mkpath("downloads")
  file = "downloads/$tarball"
  try
    url = "github.com/sciapp/gr/releases/download/v$version/$tarball"
    download("https://$url", file)
  catch
    url = "gr-framework.org/downloads/$tarball"
    try
      download("https://$url", file)
    catch
      @info("Using insecure connection")
      try
        download("http://$url", file)
      catch
        @info("Cannot download GR run-time")
      end
    end
  end
  if os == :Windows
    home = (VERSION < v"0.7-") ? JULIA_HOME : Sys.BINDIR
    success(`$home/7z x downloads/$tarball -y`)
    rm("downloads/$tarball")
    tarball = tarball[1:end-3]
    success(`$home/7z x $tarball -y -ttar`)
    rm(tarball)
  else
    run(`tar xzf downloads/$tarball`)
    rm("downloads/$tarball")
  end
  if os == :Darwin
    app = joinpath("gr", "Applications", "GKSTerm.app")
    run(`/System/Library/Frameworks/CoreServices.framework/Frameworks/LaunchServices.framework/Support/lsregister -f $app`)
    try
@static if VERSION >= v"0.7.0-DEV.3656"
      have_qml = haskey(installed(), "QML")
else
      have_qml = Pkg.installed("QML") != nothing
end
      if have_qml
        @eval import QML
        qt = QML.qt_prefix_path()
        path = joinpath(qt, "Frameworks")
        if isdir(path)
          qt5plugin = joinpath(pwd(), "gr", "lib", "qt5plugin.so")
          run(`install_name_tool -add_rpath $path $qt5plugin`)
          println("Using Qt ", splitdir(qt)[end], " at ", qt)
        end
      end
    catch
    end
  end
end


text:
module MPCCModelManager

# TODO perhaps make more generic by defining all the operations in a custom sense.  Perhaps a 2.0 type of thing.

using SparseArrays, ForwardDiff, LinearAlgebra, UnPack, RuntimeGeneratedFunctions, SymbolicUtils, Symbolics
# StaticArrays

const MPCCMM_ROOT_DIR = @__DIR__

RuntimeGeneratedFunctions.init(@__MODULE__)

@warn "Changed gradients from column to row vectors"

include("mpccmodel_common.jl")
include("mpccmodel_proc_model.jl")
include("mpccmodel_calc_forwarddiff_dense.jl")
include("mpccmodel_calc_symdiff_sparse.jl")
include("mpccmodel_calc_forwarddiff_newton_penalty.jl")
include("mpccmodel_pointeval.jl")
include("mpccmodel_test.jl")

export  mpccmodel_build_fn_from_defn,
        mpccmodel_load_defn_from_file,
        mpccmodel_build_fixed_jump_fns,
        mpccmodel_setup_forwarddiff_dense,
        mpccmodel_setup_symdiff_sparse,
        mpccmodel_setup_newton_penalty,
        mpccmodel_pointeval_basic,
        mpccmodel_test_testvectors_evalpt_basic

# More recent note: having pr and ps, then gradp, etc, is a bit werid.  cannot take derivative wrt ps, so perhaps change notation

# TODO
# - remove line numbers from generated functions, inline what need inlined, remove bounds checks where appropriate
# - deal with transpose for ce, ci.
# - remove f() in the knownsolutions eval, unncessary and creates a headache

 # - for forwarddiff processing, cast results properly

# NOTE(i): We leave the return types fairly vague and allow caller to specify whether they want sparse, dense, or even static arrays.

# NOTE reminders to self:
# - Type stability: return types should not depend on values but can depend on input types
# - Use argument types for multiple dispatch, it does not improve performance in general
# - Explicit return type annotation seems to degrate performance

# TODO should probably inline the small functions in each model definition
# TODO check ForwardDiff cfg
# TODO saving of lower order derivatives?

end


text:
#!/usr/bin/env julia

module PE13

const data=BigInt[
37107287533902102798797998220837590246510135740250,
46376937677490009712648124896970078050417018260538,
74324986199524741059474233309513058123726617309629,
91942213363574161572522430563301811072406154908250,
23067588207539346171171980310421047513778063246676,
89261670696623633820136378418383684178734361726757,
28112879812849979408065481931592621691275889832738,
44274228917432520321923589422876796487670272189318,
47451445736001306439091167216856844588711603153276,
70386486105843025439939619828917593665686757934951,
62176457141856560629502157223196586755079324193331,
64906352462741904929101432445813822663347944758178,
92575867718337217661963751590579239728245598838407,
58203565325359399008402633568948830189458628227828,
80181199384826282014278194139940567587151170094390,
35398664372827112653829987240784473053190104293586,
86515506006295864861532075273371959191420517255829,
71693888707715466499115593487603532921714970056938,
54370070576826684624621495650076471787294438377604,
53282654108756828443191190634694037855217779295145,
36123272525000296071075082563815656710885258350721,
45876576172410976447339110607218265236877223636045,
17423706905851860660448207621209813287860733969412,
81142660418086830619328460811191061556940512689692,
51934325451728388641918047049293215058642563049483,
62467221648435076201727918039944693004732956340691,
15732444386908125794514089057706229429197107928209,
55037687525678773091862540744969844508330393682126,
18336384825330154686196124348767681297534375946515,
80386287592878490201521685554828717201219257766954,
78182833757993103614740356856449095527097864797581,
16726320100436897842553539920931837441497806860984,
48403098129077791799088218795327364475675590848030,
87086987551392711854517078544161852424320693150332,
59959406895756536782107074926966537676326235447210,
69793950679652694742597709739166693763042633987085,
41052684708299085211399427365734116182760315001271,
65378607361501080857009149939512557028198746004375,
35829035317434717326932123578154982629742552737307,
94953759765105305946966067683156574377167401875275,
88902802571733229619176668713819931811048770190271,
25267680276078003013678680992525463401061632866526,
36270218540497705585629946580636237993140746255962,
24074486908231174977792365466257246923322810917141,
91430288197103288597806669760892938638285025333403,
34413065578016127815921815005561868836468420090470,
23053081172816430487623791969842487255036638784583,
11487696932154902810424020138335124462181441773470,
63783299490636259666498587618221225225512486764533,
67720186971698544312419572409913959008952310058822,
95548255300263520781532296796249481641953868218774,
76085327132285723110424803456124867697064507995236,
37774242535411291684276865538926205024910326572967,
23701913275725675285653248258265463092207058596522,
29798860272258331913126375147341994889534765745501,
18495701454879288984856827726077713721403798879715,
38298203783031473527721580348144513491373226651381,
34829543829199918180278916522431027392251122869539,
40957953066405232632538044100059654939159879593635,
29746152185502371307642255121183693803580388584903,
41698116222072977186158236678424689157993532961922,
62467957194401269043877107275048102390895523597457,
23189706772547915061505504953922979530901129967519,
86188088225875314529584099251203829009407770775672,
11306739708304724483816533873502340845647058077308,
82959174767140363198008187129011875491310547126581,
97623331044818386269515456334926366572897563400500,
42846280183517070527831839425882145521227251250327,
55121603546981200581762165212827652751691296897789,
32238195734329339946437501907836945765883352399886,
75506164965184775180738168837861091527357929701337,
62177842752192623401942399639168044983993173312731,
32924185707147349566916674687634660915035914677504,
99518671430235219628894890102423325116913619626622,
73267460800591547471830798392868535206946944540724,
76841822524674417161514036427982273348055556214818,
97142617910342598647204516893989422179826088076852,
87783646182799346313767754307809363333018982642090,
10848802521674670883215120185883543223812876952786,
71329612474782464538636993009049310363619763878039,
62184073572399794223406235393808339651327408011116,
66627891981488087797941876876144230030984490851411,
60661826293682836764744779239180335110989069790714,
85786944089552990653640447425576083659976645795096,
66024396409905389607120198219976047599490197230297,
64913982680032973156037120041377903785566085089252,
16730939319872750275468906903707539413042652315011,
94809377245048795150954100921645863754710598436791,
78639167021187492431995700641917969777599028300699,
15368713711936614952811305876380278410754449733078,
40789923115535562561142322423255033685442488917353,
44889911501440648020369068063960672322193204149535,
41503128880339536053299340368006977710650566631954,
81234880673210146739058568557934581403627822703280,
82616570773948327592232845941706525094512325230608,
22918802058777319719839450180888072429661980811197,
77158542502016545090413245809786882778948721859617,
72107838435069186155435662884062257473692284509516,
20849603980134001723930671666823555245252804609722,
53503534226472524250874054075591789781264330331690];

end

export pe13

pe13(data::Vector{<:Integer}=PE13.data)=parse(Int,string(sum(data))[1:10])

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe13()|>print
end

max_stars_repo_path: 13. Large sum/13.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 13. Large sum/13.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 13. Large sum/13.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
"""
Functions to deal with various types of norms and seminorms
"""

using FastRounding
using IntervalOptimisation
using TaylorSeries: Taylor1
using SparseArrays: getcolptr
using .DynamicDefinition

using .DynamicDefinition: derivative
using .PwDynamicDefinition: PwMap

"""
'Absolute value' definition that returns mag(I) for an interval and abs(x) for a real
"""
abs_or_mag(x::Number) = abs(x)
abs_or_mag(x::Interval) = mag(x)

"""
Certified upper bound to ||A|| (of specified NormKind)
"""
function opnormbound(::Type{L1}, A::AbstractVecOrMat{T}) where {T}
    # partly taken from JuliaLang's LinearAlgebra/src/generic.jl
    Tnorm = typeof(abs_or_mag(float(real(zero(T)))))
    Tsum = promote_type(Float64, Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for j = 1:size(A, 2)
            nrmj::Tsum = 0
            for i = 1:size(A, 1)
                nrmj = nrmj ⊕₊ abs_or_mag(A[i,j])
            end
            nrm = max(nrm,nrmj)
        end
    end
    return convert(Tnorm, nrm)
end

function opnormbound(::Type{Linf}, A::AbstractVecOrMat{T}) where {T}
    # partly taken from JuliaLang's LinearAlgebra/src/generic.jl
    Tnorm = typeof(abs_or_mag(float(real(zero(T)))))
    Tsum = promote_type(Float64, Tnorm)
    nrm::Tsum = 0
    @inbounds begin
        for i = 1:size(A, 1)
            nrmi::Tsum = 0
            for j = 1:size(A, 2)
                nrmi = nrmi ⊕₊ abs_or_mag(A[i,j])
            end
            nrm = max(nrm,nrmi)
        end
    end
    return convert(Tnorm, nrm)
end

function opnormbound(::Type{L1}, A::SparseMatrixCSC) where {T}
    # partly taken from JuliaLang's Sparsearray/src/linalg.jl
    m, n = size(A)
    Tnorm = typeof(abs_or_mag(float(real(zero(eltype(A))))))
    Tsum = promote_type(Float64,Tnorm)
    nA::Tsum = 0
    @inbounds begin
        for j=1:n
            colSum::Tsum = 0
            for i = getcolptr(A)[j]:getcolptr(A)[j+1]-1
                colSum = colSum ⊕₊ abs_or_mag(nonzeros(A)[i])
            end
            nA = max(nA, colSum)
        end
    end
    return convert(Tnorm, nA)
end

function opnormbound(::Type{Linf}, A::SparseMatrixCSC) where {T}
    # partly taken from JuliaLang's Sparsearray/src/linalg.jl
    m, n = size(A)
    Tnorm = typeof(abs_or_mag(float(real(zero(eltype(A))))))
    Tsum = promote_type(Float64,Tnorm)
    rowSum = zeros(Tsum,m)
    @inbounds begin
        for i=1:length(nonzeros(A))
            rowSum[rowvals(A)[i]] = rowSum[rowvals(A)[i]] ⊕₊ abs_or_mag(nonzeros(A)[i])
        end
    end
    return convert(Tnorm, maximum(rowSum))
end

"""
Rigorous upper bound on a vector norm. Note that Linf, L1 are the "analyst's" norms
"""
normbound(N::Type{L1}, v::AbstractVector) = opnormbound(L1, v) ⊘₊ Float64(length(v), RoundDown)
normbound(N::Type{Linf}, v::AbstractVector) = opnormbound(Linf, v)

"""
Types to compute norms iteratively by "adding a column at a time".
"""
abstract type NormCacher{T} end

mutable struct NormCacherL1 <: NormCacher{L1}
    C::Float64
    function NormCacherL1(n)
        new(0)
    end
end
"""
Create a new NormCacher to compute the normbound of the empty matrix with n rows
"""
NormCacher{L1}(n) = NormCacherL1(n)

mutable struct NormCacherLinf <: NormCacher{Linf}
    C::Vector{Float64}
    function NormCacherLinf(n)
        new(zeros(n))
    end
end
NormCacher{Linf}(n) = NormCacherLinf(n)

"""
Update a NormCacher to add one column to the matrix it is computing a norm of.
This column may be affected by an error ε (in the same norm).
"""
function add_column!(Cacher::NormCacherL1, v::AbstractVector, ε::Float64)
    Cacher.C = max(Cacher.C, opnormbound(L1, v) ⊕₊ ε)
end

function add_column!(Cacher::NormCacherLinf, v::AbstractVector, ε::Float64)
    Cacher.C = Cacher.C .⊕₊ abs.(v) .⊕₊ ε
end

"""
Return the norm of the matrix the NormCacher is working on.
"""
function get_norm(Cacher::NormCacherL1)
    return Cacher.C
end

function get_norm(Cacher::NormCacherLinf)
    return maximum(Cacher.C)
end

"""
(A, B) = dfly(strongnorm, auxnorm, dynamic)

Constants (A, B) such that ||Lf||_s ≦ A||f||_s + B||f||_aux
"""
dfly(::Type{<:NormKind}, ::Type{<:NormKind}, ::Dynamic) = @error "Not implemented"

function dfly(::Type{TotalVariation}, ::Type{L1}, D::Dynamic)
    dist = maximise(x -> distorsion(D, x), domain(D))[1]
    lam = maximise(x-> abs(1/derivative(D, x)), domain(D))[1]

    if !(abs(lam) < 1) # these are intervals, so this is *not* equal to abs(lam) >= 1.
        @error "The function is not expanding"
    end

    if is_full_branch(D)
        return lam.hi, dist.hi
    else
        if !(abs(lam) < 0.5)
            @error "Expansivity is insufficient to prove a DFLY. Try with an iterate."
        end
        # We need a way to estimate the branch widths
        @error "Not implemented"
    end
end

function dfly(::Type{TotalVariation}, ::Type{L1}, D::PwMap)
    dist = 0
    lam = 0
    disc = 0
    for i in 1:nbranches(D)
        f(x) = D.Ts[i](x)
        domain = hull(D.endpoints[i], D.endpoints[i+1])
        fprime(x) = f(Taylor1([x, 1], 1))[1]
        fsecond(x) = f(Taylor1([x, 1], 2))[2]*2
        distorsion(x)=abs(fsecond(x)/(fprime(x)^2))
        lambda(x) = abs(1/fprime(x))
        dist = max(dist, maximise(distorsion, domain)[1].hi)
        lam = max(lam, maximise(lambda, domain)[1].hi)
        low_rad = (abs(D.endpoints[i]-D.endpoints[i+1])/2).lo
        disc = max(disc, ((1/Interval(low_rad)).hi))
    end

    if is_full_branch(D)
        if !(abs(lam) < 1) # these are intervals, so this is *not* equal to abs(lam) >= 1.
            @error "The function is not expanding"
        end
        return lam, dist
    else
        if !(abs(2⊗₊lam) < 1)
            @error "Expansivity is insufficient to prove a DFLY. Try with an iterate."
        end
        return 2⊗₊lam, dist ⊕₊ disc
    end
end

function dfly(::Type{Lipschitz}, ::Type{L1}, D::Dynamic)
    # TODO: should assert that D is globally C2 instead, but we don't have that kind of infrastructure yet.
    @assert is_full_branch(D)

    dist = maximise(x -> distorsion(D, x), domain(D))[1]
    lam = maximise(x-> abs(1/derivative(D, x)), domain(D))[1]

    lam = lam.hi
    dist = dist.hi

    return lam ⊗₊ (2. ⊗₊  dist ⊕₊ 1.), dist ⊗₊ (2. ⊗₊ dist ⊕₊ 1.)
end

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
############################################################################################
module BaytesPMCMC

############################################################################################
#Import modules
using BaytesCore:
    BaytesCore,
    AbstractAlgorithm,
    AbstractTune,
    AbstractConfiguration,
    AbstractDiagnostics,
    AbstractKernel,
    AbstractConstructor,
    update,
    AcceptStatistic,
    SampleDefault

import BaytesCore:
    BaytesCore,
    update!,
    infer,
    results,
    init,
    init!,
    propose,
    propose!,
    propagate,
    propagate!,
    result!,
    get_result,
    get_sym,
    generate_showvalues,
    UpdateBool,
    UpdateTrue,
    UpdateFalse

using ModelWrappers:
    ModelWrappers,
    ModelWrapper,
    Tagged,
    Objective,
    DiffObjective,
    ℓObjectiveResult,
    dynamics,
    _checkprior

import ModelWrappers: ModelWrappers, predict, generate

using BaytesMCMC, BaytesFilters

using DocStringExtensions:
    DocStringExtensions, TYPEDEF, TYPEDFIELDS, FIELDS, SIGNATURES, FUNCTIONNAME
using ArgCheck: ArgCheck, @argcheck
using UnPack: UnPack, @unpack, @pack!

using Random: Random, AbstractRNG, GLOBAL_RNG

############################################################################################
# Define abstract types
abstract type PMCMCKernel <: AbstractKernel end

############################################################################################
# Import sub-folder
include("Core/Core.jl")
include("Kernels/Kernels.jl")

############################################################################################
export
    #BaytesCore
    UpdateBool,
    UpdateTrue,
    UpdateFalse,

    PMCMCKernel,
    update!,
    init,
    init!,
    propose,
    propose!,
    propagate,
    propagate!,
    SampleDefault,

    #ModelWrappers
    dynamics,
    predict,
    generate
end


text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

@testset "row permutation of sparse matrices" begin
  @testset "2x2" begin
    mat = sparse([1.0 0.0; 0.0 1.0])
    FirstOrderLp.row_permute_in_place(mat, [2, 1])
    @test mat == sparse([0.0 1.0; 1.0 0.0])
  end

  @testset "3x2" begin
    mat = sparse([
      1.0 0.0
      0.0 1.0
      2.0 3.0
    ])
    FirstOrderLp.row_permute_in_place(mat, [3, 1, 2])
    @test mat == sparse([
      0.0 1.0
      2.0 3.0
      1.0 0.0
    ])
  end
end


text:
using WAV
using SpeechAugment: clipWav

# degree of distortion
degree = 0.5
# read a wav file
wav = wavread("xxx.wav")[1]
# clip wav to generate distortion effect
y = clipWav(wav, degree)

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
using FinancialAssets, FactCheck

include("assets.jl")
include("id.jl")

exitstatus()


text:
@testset "Test Bias Evaluate" begin
    b = Bias(5)

    x = [1, 2, 3, 4, 5]
    Θ = [1, 2, 3, 4, 5]
    z = z_buffer(b)
    dx = dx_buffer(b)
    dΘ = dΘ_buffer(b)

    evaluate!(b, x, Θ, z, dx=dx, dΘ=dΘ)

    @test z == [2, 4, 6, 8, 10]
    @test dx == Diagonal(ones(5))
    @test dΘ == Diagonal(ones(5))
end


text:
# This is a wrapper for an int that allows the next empty index of an array to be tracked when the array is being filled by recursive functions
type Counter
	N::Int
end


text:
"""
    struct TorusMantle{T,TP,TT,D} <: AbstractCurvedSurfacePrimitive{T}

Surface primitive describing the mantle of a [`Torus`](@ref).

## Parametric types
* `T`: Precision type.
* `TP`: Type of the azimuthial angle `φ`.
    * `TP == Nothing`: Full 2π in `φ`.
* `TT`: Type of the polar angle `θ`.
    * `TT == Nothing`: Full 2π in `θ`.
* `D`: Direction in which the normal vector points (`:inwards` or `:outwards`).

## Fields
* `r_torus::T`: Distance of the center of the `TorusMantle` to the center of the tube (in m).
* `r_tube::T`: Radius of the tube of the `TorusMantle` (in m).
* `φ::TP`: Range in azimuthial angle `φ` of the `TorusMantle`.
* `θ::TT`: Range in polar angle `θ` of the `TorusMantle`.
* `origin::CartesianPoint{T}`: The position of the center of the `TorusMantle`.
* `rotation::SMatrix{3,3,T,9}`: Matrix that describes a rotation of the `TorusMantle` around its `origin`.
"""
@with_kw struct TorusMantle{T,TP,TT,D} <: AbstractCurvedSurfacePrimitive{T}
    r_torus::T = 1
    r_tube::T = 1
    φ::TP = nothing
    θ::TT = nothing

    origin::CartesianPoint{T} = zero(CartesianPoint{T})
    rotation::SMatrix{3,3,T,9} = one(SMatrix{3, 3, T, 9})
end

flip(t::TorusMantle{T,TP,TT,:inwards}) where {T,TP,TT} =
TorusMantle{T,TP,TT,:outwards}(t.r_torus, t.r_tube, t.φ, t.θ, t.origin, t.rotation )

get_φ_limits(tm::TorusMantle{T,Tuple{T,T}}) where {T} = tm.φ[1], tm.φ[2]
get_φ_limits(tm::TorusMantle{T,Nothing}) where {T} = T(0), T(2π)

get_θ_limits(tm::TorusMantle{T,<:Any,Tuple{T,T}}) where {T} = tm.θ[1], tm.θ[2]
get_θ_limits(tm::TorusMantle{T,<:Any,Nothing}) where {T} = T(0), T(2π)

function normal(tm::TorusMantle{T,TP,TT,:outwards}, pt::CartesianPoint{T}) where {T,TP,TT}
    pto = _transform_into_object_coordinate_system(pt, tm)
    cyl = CylindricalPoint(pto)
    ptt = CartesianPoint(CylindricalPoint{T}(tm.r_torus, cyl.φ, zero(T)))
    return pt - _transform_into_global_coordinate_system(ptt, tm)
end
normal(tm::TorusMantle{T,TP,TT,:inwards}, pt::CartesianPoint{T}) where {T,TP,TT} = -normal(flip(tm), pt)

function vertices(tm::TorusMantle{T}, n_arc::Int64)::Vector{CartesianPoint{T}} where {T}
    φMin, φMax = get_φ_limits(tm)
    θMin, θMax = get_θ_limits(tm)
    n_arcφ = _get_n_points_in_arc_φ(tm, n_arc)
    n_arcθ = _get_n_points_in_arc_θ(tm, n_arc)

    θ = range(θMin, θMax, length = n_arcθ + 1)
    φ = range(φMin, φMax, length = n_arcφ + 1)

    [_transform_into_global_coordinate_system(CartesianPoint{T}((tm.r_torus + tm.r_tube*cos(θ))*cos(φ), (tm.r_torus + tm.r_tube*cos(θ))*sin(φ), tm.r_tube*sin(θ)), tm) for θ in θ for φ in φ]
end

function sample(tm::TorusMantle{T}, spacing::T)::Vector{CartesianPoint{T}} where {T}
    φMin, φMax = get_φ_limits(tm)
    θMin, θMax = get_θ_limits(tm)
    Δφ = abs(φMax - φMin)
    Δθ = abs(θMax - θMin)
    full2πφ = isnothing(tm.φ)
	full2πθ = isnothing(tm.θ)

	[_transform_into_global_coordinate_system(CartesianPoint{T}((tm.r_torus + tm.r_tube*cos(θ))*cos(φ), (tm.r_torus + tm.r_tube*cos(θ))*sin(φ), tm.r_tube*sin(θ)), tm) for θ in range(θMin, θMax - full2πθ*spacing/tm.r_tube, length = max(2,1+Int(ceil(Δθ*tm.r_tube/spacing)))) for φ in range(φMin, φMax - full2πφ*spacing/(tm.r_torus + tm.r_tube*cos(θ)), length = max(2,1+Int(ceil(Δφ*(tm.r_torus + tm.r_tube*cos(θ))/spacing))))]
end

function connections(tm::TorusMantle, n_arc::Int64)::Vector{Vector{Int64}}
    n_arcφ = _get_n_points_in_arc_φ(tm, n_arc)
    n_arcθ = _get_n_points_in_arc_θ(tm, n_arc)
    [[i+(n_arcφ+1)*j,i+1+(n_arcφ+1)*j,i+1+(n_arcφ+1)*(j+1),i+(n_arcφ+1)*(j+1)] for j in 0:n_arcθ-1 for i in 1:n_arcφ]
end

function _get_hor_lines_idx(tm::TorusMantle{T}, n_arcθ::Int64)::Vector{Int64} where {T}
    θMin, θMax = get_θ_limits(tm)
    if mod(θMax-θMin, T(2π)) == 0
        l_0 = Int(floor(abs(θMin/(θMax-θMin))*n_arcθ)) + 1
        l_π = Int(floor(abs((θMin-π)/(θMax-θMin))*n_arcθ)) + 1
        [l_0,l_π]
    else
        [1, n_arcθ+1]
    end
end

function connections(tm::TorusMantle{T}, n_arc::Int64, n_vert_lines::Int64)::Vector{Vector{Int64}} where {T}
    n_arcφ = _get_n_points_in_arc_φ(tm, n_arc)
    n_arcθ = _get_n_points_in_arc_θ(tm, n_arc)
    vcircs = [[i*(n_arcφ+1)+c, (i+1)*(n_arcφ+1)+c] for c in _get_vert_lines_range(tm,n_arcφ,n_vert_lines) for i in 0:n_arcθ-1]
    hcircs =  [[i + (c-1)*(n_arcφ+1), i + (c-1)*(n_arcφ+1) + 1] for c in _get_hor_lines_idx(tm, n_arcθ) for i in 1:n_arcφ]
    append!(vcircs, hcircs)
end

get_label_name(::TorusMantle) = "Torus Mantle"

const FullTorusMantle{T,D} = TorusMantle{T,Nothing,Nothing,D}

function lines(tm::FullTorusMantle{T}) where {T}
    top_circ_origin = CartesianPoint{T}(zero(T), zero(T),  tm.r_tube)
    top_circ_origin = _transform_into_global_coordinate_system(top_circ_origin, tm)
    top_circ = Circle{T}(r = tm.r_torus, origin = top_circ_origin, rotation = tm.rotation)
    bot_circ_origin = CartesianPoint{T}(zero(T), zero(T), -tm.r_tube)
    bot_circ_origin = _transform_into_global_coordinate_system(bot_circ_origin, tm)
    bot_circ = Circle{T}(r = tm.r_torus, origin = bot_circ_origin, rotation = tm.rotation)
    inner_circ = Circle{T}(r = tm.r_torus - tm.r_tube, origin = tm.origin, rotation = tm.rotation)
    outer_circ = Circle{T}(r = tm.r_torus + tm.r_tube, origin = tm.origin, rotation = tm.rotation)
    tube_circ_1_origin = CartesianPoint{T}(tm.r_torus, zero(T), zero(T))
    tube_circ_1_origin = _transform_into_global_coordinate_system(tube_circ_1_origin, tm)
    tube_circ_1 = Circle{T}(r = tm.r_tube, origin = tube_circ_1_origin, rotation = tm.rotation * RotX(T(π)/2))
    tube_circ_2_origin = CartesianPoint{T}(-tm.r_torus, zero(T), zero(T))
    tube_circ_2_origin = _transform_into_global_coordinate_system(tube_circ_2_origin, tm)
    tube_circ_2 = Circle{T}(r = tm.r_tube, origin = tube_circ_2_origin, rotation = tm.rotation * RotX(T(π)/2))
    (bot_circ, outer_circ, top_circ, inner_circ, tube_circ_1, tube_circ_2)
end

extremum(tm::TorusMantle{T}) where {T} = tm.r_torus + tm.r_tube

"""
    intersection(tm::TorusMantle{T}, l::Line{T}) where {T}

Calculates the intersections of a `Line` with a `TorusMantle`.

## Arguments
* `cm::TorusMantle{T}`: The `TorusMantle`.
* `l::Line{T}`: The `Line`.

!!! note
    The function will always return 4 CartesianPoint's.
    If the line just touches the mantle, the points will be the same.
    If the line does not touch the mantle at all, the points will have NaN's as there coordinates.
"""
function intersection(tm::TorusMantle{T}, l::Line{T}) where {T}
    obj_l = _transform_into_object_coordinate_system(l, tm) # direction is not normalized

    L1 = obj_l.origin.x
    L2 = obj_l.origin.y
    L3 = obj_l.origin.z
    D1 = obj_l.direction.x
    D2 = obj_l.direction.y
    D3 = obj_l.direction.z

    R = tm.r_torus
    r = tm.r_tube

    A = L1^2 + L2^2 + L3^2 + R^2 - r^2
    B = 2*(L1*D1 + L2*D2 + L3*D3)
    C = D1^2 + D2^2 + D3^2

    a = (2*B*C) / C^2
    b = (2*A*C + B^2 - 4*R^2*(D1^2 + D2^2)) / C^2
    c = (2*A*B - 8*R^2*(L1*D1 + L2*D2)) / C^2
    d = (A^2 - 4*R^2*(L1^2 + L2^2)) / C^2

    # Solve: `solve (sqrt((L1 + λ*D1)^2 + (L2 + λ*D2)^2)-R)^2 + (L3 + λ*D3)^2 = r^2 for λ`

	# λ1, λ2, λ3, λ4 = roots_of_4th_order_polynomial(a, b, c, d) # That does not work for all combinations of a, b, c, d...
	# fallback to Polynomials.jl, which is slower... We should improve `roots_of_4th_order_polynomial`...
    return broadcast(λ -> _transform_into_global_coordinate_system(obj_l.origin + λ * obj_l.direction, tm),
        real.(Polynomials.roots(Polynomial((d, c, b, a, one(T))))))
end

# """
#     roots_of_4th_order_polynomial(a::T, b::T, c::T, d::T, e::T)
#
# Calculate the 4 (possible) roots of `x^4 + ax^3 + bx^2 + cx + d = 0`
# """
# function roots_of_4th_order_polynomial(a::T, b::T, c::T, d::T) where {T}
# 	#=
# 	using Polynomials
# 	A, a, b, c, d = 1.0, 0.2, -5.0, -0.1, 2.5
# 	p = Polynomial(SVector{5,T}([d, c, b, a, A]))
# 	xs = -2.5:0.01:2.5; plot(xs, map(x -> p(x), xs))
# 	roots(p)
#
# 	rs = CSG.roots_of_4th_order_polynomial(a, b, c, d)
# 	vline!([rs...])
# 	@btime CSG.roots_of_4th_order_polynomial($a, $b, $c, $d)
# 	=#
#
# 	# There are some issues for certain combinations of a, b, c, d, ...
#
# 	term_1_1_1_1 = 2*b^3 - 9*a*b*c + 27*c^2 + 27*a^2*d - 72*b*d
# 	term_1_1_1_2 = b^2 - 3*a*c + 12*d
#
# 	comp_term = -4 * term_1_1_1_2^3 + term_1_1_1_1^2
#
# 	term_1_1_1_3 = term_1_1_1_1 + sqrt(Complex(comp_term))
# 	term_1_1_1 = 3*(term_1_1_1_3^(1/3))
#
# 	term_1_1 = 2^(1/3) * term_1_1_1_2 / term_1_1_1
# 	term_1_2 = (term_1_1_1_3 / 54)^(1/3)
# 	term_1 = a^2/4 - 2b/3 + term_1_1 + term_1_2
#
# 	if term_1 == 0
# 		term_6 = sqrt(Complex(b^2 - 4d))
# 		λ1 = sqrt(-term_6 - b)/sqrt(2)
# 		λ2 = -λ1
# 		λ3 = sqrt( term_6 - b)/sqrt(2)
# 		λ4 = -λ3
# 	else
# 		term_2_1 = (-a^3 + 4a*b - 8c) / (4*sqrt(term_1))
# 		term_2a = a^2/2 - 4b/3 - term_1_2 - term_1_2 - term_2_1
# 		term_2b = a^2/2 - 4b/3 - term_1_1 - term_1_2 + term_2_1
#
# 		term3  = sqrt(term_1)/2
# 		term4a = sqrt(term_2a)/2
# 		term4b = sqrt(term_2b)/2
# 		term5 = -a/4
#
# 		λ1 = term5 - term3 - term4a
# 		λ2 = term5 - term3 + term4a
# 		λ3 = term5 + term3 - term4b
# 		λ4 = term5 + term3 + term4b
# 	end
# 	real(λ1), real(λ2), real(λ3), real(λ4)
# end

# function distance_to_surface(pt::AbstractCoordinatePoint{T}, t::TorusMantle{T, Nothing})::T where {T}
#     pcy = CylindricalPoint(pt)
#     return distance_to_line(PlanarPoint{T}(pcy.r,pcy.z), Arc(t))
# end

# function distance_to_surface(pt::AbstractCoordinatePoint{T}, t::TorusMantle{T, <:AbstractInterval})::T where {T}
#     pcy = CylindricalPoint(pt)
#     if _in_φ(pt, t.φ)
#         return distance_to_line(PlanarPoint{T}(pcy.r,pcy.z), Arc(t))
#     else
#         φMin::T, φMax::T, _ = get_φ_limits(t)
#         Δφ = pcy.φ - _φNear(pcy.φ, φMin, φMax)
#         d, r_on_plane = pcy.r .* sincos(Δφ)
#         return hypot(d, distance_to_line(PlanarPoint{T}(r_on_plane, pcy.z), Arc(t)))
#     end
# end