text:
struct TrialFESpace{B} <: SingleFieldFESpace
  space::SingleFieldFESpace
  dirichlet_values::AbstractVector
  cell_basis::CellField
  constraint_style::Val{B}

  function TrialFESpace(dirichlet_values::AbstractVector,space::SingleFieldFESpace)
    cell_basis = trialize_cell_basis(get_cell_basis(space))
    cs = constraint_style(space)
    B = get_val_parameter(cs)
    new{B}(space,dirichlet_values,cell_basis,cs)
  end
end

"""
"""
function TrialFESpace(space::SingleFieldFESpace)
  dirichlet_values = get_dirichlet_values(space)
  TrialFESpace(dirichlet_values,space)
end

"""
"""
function TrialFESpace(space::SingleFieldFESpace,objects)
  dirichlet_values = compute_dirichlet_values_for_tags(space,objects)
  TrialFESpace(dirichlet_values,space)
end

"""
"""
function TrialFESpace!(dir_values::AbstractVector,space::SingleFieldFESpace,objects)
  dir_values = compute_dirichlet_values_for_tags!(dir_values,space,objects)
  TrialFESpace(dir_values,space)
end

"""
"""
function TrialFESpace!(space::TrialFESpace,objects)
  dir_values = get_dirichlet_values(space)
  dir_values = compute_dirichlet_values_for_tags!(dir_values,space,objects)
  space
end

function TrialFESpace(space::TrialFESpace)
  space
end

function HomogeneousTrialFESpace(U::FESpace)
  dirichlet_values = zero_dirichlet_values(U)
  TrialFESpace(dirichlet_values,U)
end

function HomogeneousTrialFESpace!(dirichlet_values::AbstractVector,U::FESpace)
  fill!(dirichlet_values,zero(eltype(dirichlet_values)))
  TrialFESpace(dirichlet_values,U)
end

# Genuine functions

get_dirichlet_values(f::TrialFESpace) = f.dirichlet_values

get_cell_basis(f::TrialFESpace) = f.cell_basis

# Delegated functions

get_cell_axes(t::TrialFESpace)= get_cell_axes(t.space)

get_cell_axes_with_constraints(t::TrialFESpace)= get_cell_axes_with_constraints(t.space)

CellData.CellField(t::TrialFESpace,cell_vals) = CellField(t.space,cell_vals)

constraint_style(::Type{<:TrialFESpace{B}}) where B = Val{B}()

get_cell_dof_basis(f::TrialFESpace) = get_cell_dof_basis(f.space)

num_free_dofs(f::TrialFESpace) = num_free_dofs(f.space)

zero_free_values(f::TrialFESpace) = zero_free_values(f.space)

get_cell_dofs(f::TrialFESpace) = get_cell_dofs(f.space)

num_dirichlet_dofs(f::TrialFESpace) = num_dirichlet_dofs(f.space)

zero_dirichlet_values(f::TrialFESpace) = zero_dirichlet_values(f.space)

num_dirichlet_tags(f::TrialFESpace) = num_dirichlet_tags(f.space)

get_dirichlet_dof_tag(f::TrialFESpace) = get_dirichlet_dof_tag(f.space)

scatter_free_and_dirichlet_values(f::TrialFESpace,fv,dv) = scatter_free_and_dirichlet_values(f.space,fv,dv)

gather_free_and_dirichlet_values(f::TrialFESpace,cv) = gather_free_and_dirichlet_values(f.space,cv)

gather_free_and_dirichlet_values!(fv,dv,f::TrialFESpace,cv) = gather_free_and_dirichlet_values!(fv,dv,f.space,cv)

gather_dirichlet_values(f::TrialFESpace,cv) = gather_dirichlet_values(f.space,cv)

gather_dirichlet_values!(dv,f::TrialFESpace,cv) = gather_dirichlet_values!(dv,f.space,cv)

gather_free_values(f::TrialFESpace,cv) = gather_free_values(f.space,cv)

gather_free_values!(fv,f::TrialFESpace,cv) = gather_free_values!(fv,f.space,cv)

get_cell_isconstrained(f::TrialFESpace) = get_cell_isconstrained(f.space)

get_cell_constraints(f::TrialFESpace) = get_cell_constraints(f.space)


text:
include("abstract_agent.jl")
include("agent.jl")
include("dyna_agent.jl")

ext: jl


text:
# This file was generated, do not modify it. # hide
using MLJ, RDatasets, Random

data = dataset("datasets", "USArrests")
names(data)


text:
"Load data from file."
function loaddata(file::String)
    return readlines(file)
end

"Check password according to sled rental policy."
function checksled(low::Int, up::Int, char::Char, password::AbstractString)
    return count(string(char), password) ∈ low:up
end

"Check password according to toboggan rental policy."
function checktoboggan(low::Int, up::Int, char::Char, password::AbstractString)
    return (password[low] == char) ⊻ (password[up] == char)
end

"Count valid passwords according to given policy."
function countvalid(data::Array, rule::Function)
    n_count = 0
    for line in data
        low, up, char, password = split(replace(line, r"[-:]" => " "))
        n_count += rule(parse(Int, low), parse(Int, up), char[1], password)
    end
    return n_count
end

function main()
    data = loaddata("../assets/data/02.txt")

    # part one
    n_valid = countvalid(data, checksled)
    println("The product of the two numbers is ", n_valid, ".")

    # part two
    n_valid = countvalid(data, checktoboggan)
    println("The product of the two numbers is ", n_valid, ".")
end

main()


text:
using Automata
using Test

function allequal(arr1, arr2)
  return sum(arr1 .!= arr2) == 0
end

function tests()
  @testset "Elementary Tests" begin
    @test allequal(automaton(ElementaryRule(30), 10)[end, :], Bool[1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1])
    @test allequal(automaton(ElementaryRule(90), 10)[end, :], Bool[1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1])
    @test allequal(automaton(ElementaryRule(110), 10)[end, :], Bool[1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0])
  end
  @testset "Elementary Reversible Tests" begin
    @test allequal(automaton(ElementaryReversibleRule(30), 10)[end, :], Bool[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 0])
    @test allequal(automaton(ElementaryReversibleRule(111), 10)[end, :], Bool[0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 1, 0])
    @test allequal(automaton(ElementaryReversibleRule(137), 10)[end, :], Bool[0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0])
  end
  @testset "Totalistic Tests" begin
    @test allequal(automaton(TotalisticRule(177, 3), 10)[end, :], UInt8[0x02, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x02, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x00, 0x02])
    @test allequal(automaton(TotalisticRule(777, 3), 10)[end, :], UInt8[0x01, 0x01, 0x00, 0x00, 0x01, 0x02, 0x02, 0x02, 0x02, 0x01, 0x02, 0x02, 0x02, 0x02, 0x01, 0x00, 0x00, 0x01, 0x01])
    @test allequal(automaton(TotalisticRule(2530, 4), 10)[end, :], UInt8[0x00, 0x03, 0x02, 0x00, 0x02, 0x00, 0x00, 0x00, 0x01, 0x01, 0x01, 0x00, 0x00, 0x00, 0x02, 0x00, 0x02, 0x03, 0x00])
  end
end

tests()


text:
using Trajectories
# For plotting: discretize of the circular dynamics

"""
    discretize(x::Vector, Flow::Union{ZigZag1d, Boomerang1d}, dt)

Transform the output of the algorithm (a skeleton of points) to
a trajectory. Simple 1-d version.
"""
function discretize(x::Vector, Flow::Union{ZigZag1d, Boomerang1d}, dt0)
    k = 1
    _, ξ, θ = x[k]
    τ = x[k+1][1]
    clock = 0.0
    Ω = trajectory([clock => ξ])
    dt = dt0
    while k < length(x)-1
        while clock + dt <= τ
            if θ == 0.0
                clock += dt
            else
                clock, ξ, θ = move_forward(dt, clock, ξ, θ, Flow)
            end
            push!(Ω, clock => ξ)
            dt = dt0
        end
        Δt = τ - clock
        dt = dt - Δt
        if θ == 0.0
            clock = τ
        else
            clock, ξ, θ = move_forward(τ - clock, clock, ξ, θ, Flow)
        end
        k += 1
        _, ξ, θ = x[k]
        τ = x[k+1][1]
    end
    push!(Ω, clock => ξ)
    Ω
end

## THIS FUNCTION IS DEPRECATED. See src/trace.jl
"""
    discretize(x::Vector, Flow::Union{BouncyParticle, Boomerang}, dt)

Transform the output of the algorithm (a skeleton of points) to
a trajectory. multi-dimensional version.

Old version that would not work with the sticky Boomerang sampler not centered in 0
"""
function sdiscretize(x::Vector, Flow::Union{BouncyParticle, Boomerang}, dt0)
    k = 1
    _, ξ, θ, f = x[k]
    τ = x[k+1][1]
    clock = 0.0
    Ω = trajectory([clock => deepcopy(ξ)])
    dt = dt0
    while k < length(x)-1
        while clock + dt <= τ
            clock, ξ, θ = smove_forward!(dt, clock, ξ, θ, f, Flow)
            push!(Ω, clock => deepcopy(ξ))
            dt = dt0
        end
        Δt = τ - clock
        dt = dt - Δt
        clock, ξ, θ = smove_forward!(τ - clock, clock, ξ, θ, f, Flow)
        k += 1
        _, ξ, θ, f = x[k]
        τ = x[k+1][1]
    end
    Ω
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

function compile(libname, tarball_url, hash; prefix=BinaryProvider.global_prefix, verbose=false)
    # download to tarball_path
    tarball_path = joinpath(prefix, "downloads", libname)
    download_verify(tarball_url, hash, tarball_path; force=true, verbose=verbose)

    # unpack into source_path
    tarball_dir = joinpath(prefix, "downloads", dirname(first(list_tarball_files(tarball_path))))
    source_path = joinpath(prefix, "downloads", "src")
    verbose && @info("Unpacking $tarball_path into $source_path")
    rm(tarball_dir, force=true, recursive=true)
    rm(source_path, force=true, recursive=true)
    unpack(tarball_path, dirname(tarball_dir); verbose=verbose)
    mv(tarball_dir, source_path)

    # install libpng
    build_dir = joinpath(source_path, "build")
    mkdir(build_dir)
    verbose && @info("Compiling in $build_dir...")
    cd(build_dir) do
        # build in parallel hdf5 mode if mpi is installed
        mkpath(libdir(prefix))
        run(`$cmake_executable .. -DCMAKE_INSTALL_PREFIX=$(prefix.path)`)
        run(`make -j $(Sys.CPU_THREADS + 1)`)
        run(`make install -j $(Sys.CPU_THREADS + 1)`)
    end

    # remove old files if force=true
    manifest_path = manifest_from_url(tarball_url, prefix=prefix)
    isfile(manifest_path) && rm(manifest_path, force=true, recursive=true)

    # First, get list of files that are contained within the tarball
    file_list = list_tarball_files(tarball_path, verbose=true)

    # Save installation manifest
    mkpath(dirname(manifest_path))
    open(manifest_path, "w") do f
        for (i, path) in enumerate(file_list)
            file_list[i] = replace(file_list[i], dirname(first(list_tarball_files(tarball_path))) => "downloads/src" )
        end
        write(f, join(file_list, "\n"))
    end
end

max_stars_repo_name: FugroRoames/LibPNG

max_issues_repo_name: FugroRoames/LibPNG

max_forks_repo_name: FugroRoames/LibPNG


text:
module Day3

using DelimitedFiles

data_as_matrix(strings, bitcount) = begin
    m = zeros(Int, (length(strings), bitcount))
    for (row, string) in enumerate(strings)
        for bit in 1:bitcount
            b = parse(Int, string[bit])
            m[row, bit] = b
        end
    end
    m
end

bit_ratings(m) = begin
    hist = sum(m, dims=1) |> vec
    @assert size(hist) == (size(m, 2),)

    avgs = hist ./ size(m, 1)
    most_common = Int.(round.(avgs, RoundNearestTiesUp))
    least_common = 1 .- most_common
    (most_common, least_common)
end

rates_product(strings, bitcount) = begin
    m = data_as_matrix(strings, bitcount)

    γ_bits, ϵ_bits = bit_ratings(m)

    γ = parse(Int, join(γ_bits), base=2)
    ϵ = parse(Int, join(ϵ_bits), base=2)

    γ * ϵ
end

life_support_rating(strings, bitcount) = begin
    m = data_as_matrix(strings, bitcount)
    oxygen_candidates = co2_candidates = m

    for i in 1:bitcount
        if size(oxygen_candidates, 1) > 1
            most, _ = bit_ratings(oxygen_candidates)
            oxygen_mask = view(oxygen_candidates, :, i) .== most[i]
            oxygen_candidates = oxygen_candidates[oxygen_mask, :]
        end

        if size(co2_candidates, 1) > 1
            _, least = bit_ratings(co2_candidates)
            co2_mask = view(co2_candidates, :, i) .== least[i]
            co2_candidates = co2_candidates[co2_mask, :]
        end
    end

    oxygen_rating = parse(Int, join(oxygen_candidates[1, :]), base=2)
    co2_rating = parse(Int, join(co2_candidates[1, :]), base=2)

    oxygen_rating * co2_rating
end

test = ["00100"
        "11110"
        "10110"
        "10111"
        "10101"
        "01111"
        "00111"
        "11100"
        "10000"
        "11001"
        "00010"
        "01010"];

@assert rates_product(test, 5) == 198

@assert life_support_rating(test, 5) == 230

input = readdlm("resources/day3_input.txt", '\t', String, '\n') |> vec;

display(rates_product(input, 12))
display(life_support_rating(input, 12))

end


text:
# test_decimals.jl

@testset "Test decimal support" begin
    con = DBInterface.connect(DuckDB.DB)

    results = DBInterface.execute(
        con,
        "SELECT 42.3::DECIMAL(4,1) a, 4923.3::DECIMAL(9,1) b, 421.423::DECIMAL(18,3) c, 129481294.3392::DECIMAL(38,4) d"
    )

    # convert to DataFrame
    df = DataFrame(results)
    @test names(df) == ["a", "b", "c", "d"]
    @test size(df, 1) == 1
    @test df.a == [42.3]
    @test df.b == [4923.3]
    @test df.c == [421.423]
    @test df.d == [129481294.3392]

    DBInterface.close!(con)
end

# test returning decimals in a table function
function my_bind_function(info::DuckDB.BindInfo)
    DuckDB.add_result_column(info, "a", FixedDecimal{Int16, 0})
    DuckDB.add_result_column(info, "b", FixedDecimal{Int32, 1})
    DuckDB.add_result_column(info, "c", FixedDecimal{Int64, 2})
    DuckDB.add_result_column(info, "d", FixedDecimal{Int128, 3})
    return missing
end

mutable struct MyInitStruct
    pos::Int64

    function MyInitStruct()
        return new(0)
    end
end

function my_init_function(info::DuckDB.InitInfo)
    return MyInitStruct()
end

function my_main_function(info::DuckDB.FunctionInfo, output::DuckDB.DataChunk)
    init_info = DuckDB.get_init_info(info, MyInitStruct)

    a_array = DuckDB.get_array(output, 1, Int16)
    b_array = DuckDB.get_array(output, 2, Int32)
    c_array = DuckDB.get_array(output, 3, Int64)
    d_array = DuckDB.get_array(output, 4, Int128)
    count = 0
    multiplier = 1
    for i in 1:(DuckDB.VECTOR_SIZE)
        if init_info.pos >= 3
            break
        end
        a_array[count + 1] = 42 * multiplier
        b_array[count + 1] = 42 * multiplier
        c_array[count + 1] = 42 * multiplier
        d_array[count + 1] = 42 * multiplier
        count += 1
        init_info.pos += 1
        multiplier *= 10
    end

    DuckDB.set_size(output, count)
    return
end

@testset "Test returning decimals from a table functions" begin
    con = DBInterface.connect(DuckDB.DB)

    arguments::Vector{DataType} = Vector()
    DuckDB.create_table_function(con, "my_function", arguments, my_bind_function, my_init_function, my_main_function)
    GC.gc()

    # 3 elements
    results = DBInterface.execute(con, "SELECT * FROM my_function()")
    GC.gc()

    df = DataFrame(results)
    @test names(df) == ["a", "b", "c", "d"]
    @test size(df, 1) == 3
    @test df.a == [42, 420, 4200]
    @test df.b == [4.2, 42, 420]
    @test df.c == [0.42, 4.2, 42]
    @test df.d == [0.042, 0.42, 4.2]
end

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
# Monte Carlo solution with finite difference estimators.

"""
    get_sample_fd(sys::System{S,M}, pseudosps::NTuple{N,SamplingParameters{S,M,P}}, sp::SamplingParameters{S_,M,P})

Compute a sample for `sys` using `pseudosps` and `sp`.
"""
function get_sample_fd(sys::System{S,M}, pseudosps::NTuple{N,SamplingParameters{S,M,P}}, sp::SamplingParameters{S_,M,P}) where {S,S_,M,P,N}
    length(pseudosps) >= 1 || throw(DomainError(length(pseudosps), "At least one set of parameters."))

    # Choose a surface.
    s_ = sample(1:S_, sp.weights)

    # Sample coordinates.
    qs = zeros(P, M)
    for m in 1:M
        qs[:, m] .= rand(sp.mvns[m, s_])
    end

    # Calculate the numerator and denominator matrices.
    nums = [Matrix{Float64}(I, S, S) for _ in 1:N]
    denom = Matrix{Float64}(I, S_, S_)

    for i1 in 1:P
        # Next index in cyclic path.
        i2 = i1%P+1

        scaling = nothing

        # Numerator.
        for (j, pseudosp) in enumerate(pseudosps)
            FM_num, scaling = sampling_matrix_free_particle(pseudosp, qs[i1, :], qs[i2, :]; scaling=scaling)
            _, IM = sampling_matrix_interaction(sys, pseudosp, qs[i1, :])
            nums[j] *= IM * FM_num
        end

        # Denominator.
        FM_denom, _ = sampling_matrix_free_particle(sp, qs[i1, :], qs[i2, :]; scaling=scaling)
        denom *= FM_denom
    end

    traces_num = [tr(num) for num in nums]
    trace_denom = tr(denom)

    [trace_num/trace_denom for trace_num in traces_num]
end

"""
Monte Carlo solution for an arbitrary system using finite difference
estimators.
"""
struct SamplingFiniteDifference <: Sampling
    "Partition function."
    Z::Float64
    "Partition function standard error."
    Z_err::Float64
    "Energy."
    E::Float64
    "Energy standard error."
    E_err::Float64
    "Heat capacity."
    Cv::Float64
    "Heat capacity standard error."
    Cv_err::Float64
end

"""
    SamplingFiniteDifference(sys::System, beta::Float64, dbeta::Float64, P::Int, num_samples::Int; sampling_sys::Maybe{System}=nothing, progress_output::IO=stderr)

Calculate the solution for `sys` at `beta` with `P` links and `num_samples`
random samples, using finite difference step `dbeta`.

If `sampling_sys` is provided, it is used for sampling. Otherwise, sampling
defaults to the simplified diagonal subsystem of `sys`.

The progress meter is written to `progress_output`.
"""
function SamplingFiniteDifference(sys::System, beta::Float64, dbeta::Float64, P::Int, num_samples::Int; sampling_sys::Maybe{System}=nothing, progress_output::IO=stderr)
    sys_diag = diag(sys)
    sys_diag_simple = simplify(sys_diag)
    pseudosp = SamplingParameters(sys_diag_simple, beta, P)
    pseudosp_m = SamplingParameters(sys_diag_simple, beta-dbeta, P)
    pseudosp_p = SamplingParameters(sys_diag_simple, beta+dbeta, P)

    sampling_sys === nothing && (sampling_sys = sys_diag_simple)
    sp = SamplingParameters(sampling_sys, beta, P)

    samples = zeros(Float64, 3, num_samples)
    meter = Progress(num_samples, output=progress_output)
    for n in ProgressWrapper(1:num_samples, meter)
        samples[:, n] .= get_sample_fd(sys, (pseudosp, pseudosp_m, pseudosp_p), sp)
    end

    simple = Analytical(sampling_sys, beta)
    simple_m = Analytical(sampling_sys, beta-dbeta)
    simple_p = Analytical(sampling_sys, beta+dbeta)
    Zrat_m = simple.Z / simple_m.Z
    Zrat_p = simple.Z / simple_p.Z
    normalization = simple.Z

    f_E(samples, samples_m, samples_p) =
        simple.E .+
        1.0/(2dbeta) * (Zrat_m * samples_m .- Zrat_p * samples_p) ./ samples
    f_Cv(samples, samples_m, samples_p) =
        simple.Cv .+
        (1.0/dbeta^2 * (Zrat_m * samples_m .+ Zrat_p * samples_p) ./ samples .-
         2.0/dbeta^2 .-
         (E - simple.E)^2) * beta^2

    Z = mean(samples[1, :]) * normalization
    Z_err = std(samples[1, :]) / sqrt(num_samples) * normalization
    E, E_err = jackknife(f_E, [samples[n, :] for n in 1:size(samples, 1)]...)
    Cv, Cv_err = jackknife(f_Cv, [samples[n, :] for n in 1:size(samples, 1)]...)

    SamplingFiniteDifference(Z, Z_err, E, E_err, Cv, Cv_err)
end


text:
module TableViewExample

using Blink, CSV, DataFrames, WebIO, CSSUtil, Interact
using TableView.TableViewAdvanced

println("Starting [$(splitext(basename(@__FILE__))[1])]")

#=======================================
load the data for this example
=======================================#
current_file_folder = dirname(@__FILE__)
csv = CSV.File(joinpath(current_file_folder, "Data", "Data01.csv"))
df = DataFrame(csv)
# convert all the column to type Any to make it easier to process any type of imput and replace all missing values with spaces
for col in names(df)
    df[!,Symbol(col)] = Vector{Any}(df[!,Symbol(col)])
    replace!(df[!, col], missing => "");
end

#=======================================
define table options
=======================================#
table_options = Dict{Symbol, Any}(
    :columnDefs => Dict(
        :ObjectID => Dict(
            :headerName => "Object ID",
            :headerTooltip => "A Unique name to identify this row",
        ),
        :Comment => Dict(
            :headerName => "Comment",
        ),
        :ObjectType => Dict(
            :headerName => "Object Type",
        ),
        :Position => Dict(
            :headerName => "Position(x,y,z,θx,θy,θz)",
        ),
        :RelativeTo => Dict(
            :headerName => "Relative To",
        ),
        :Material => Dict(
            :headerName => "Material",
        ),
        :Parameters => Dict(
            :headerName => "Parameters (radium, thinkness, conic, etc)",
        ),

    ),
    :defaultColDef => Dict(
        :editable => false,
        :sortable => true,
        :resizable => true,
    ),
)

table = TableViewAdvanced.Table(df, table_options)

#=======================================
Create Blink Window
=======================================#
window_defaults = Blink.@d(
    :title => "Grid Example 01",
    :width=>1200,
    :height=>800,
)
win = Window(window_defaults)

# set the window as the table's "user data"  so we can use it inside grid events to run JS code
set_table_user_data(table, win)

#=======================================
create some of the ui widgets
=======================================#
toggle_dev_tools_button = button("Toggle Dev Tools")
on(toggle_dev_tools_button) do val
    Blink.AtomShell.tools(win)
end

show_in_repl_button = button("Show Data in REPL")
on(show_in_repl_button) do val
    @info df
end

ui_widgets = hbox(
    toggle_dev_tools_button,
    Node(:span, attributes=Dict(:style => "width: 10px;")),     # add a small space between the two buttons
    show_in_repl_button,
)

#=======================================
create some of the ui widgets
=======================================#

ui = Node(:dom,
    JS_loader_node(joinpath(dirname(@__FILE__), "examples.css")),                       # load the examples CSS file - not used in this example
    Node(:h1, "Grid Example 01"),
    Node(:h2, "In this example we create a read-only grid. We used the defaultColDef part of the options to indicate that :editable=false for all columns, unless a specific column override this option)."),
    Node(:h2, "We also show in this example how to change the header/caption of the presented columns in the grid."),
    Node(:hr),
    ui_widgets,
    Node(:hr),
    Node(:div, table_scope(table)),
)

body!(win, ui, async=false)
# Blink.AtomShell.opentools(win)

println("Ending [$(splitext(basename(@__FILE__))[1])]")
end # module TableViewExample


text:
using Revise
using AOM
using BrowseTables
##
plot_results(RequiemFusion;remove = ["Black Dahlia", "uncleb"])
##
res = update_results(RequiemFusion)
lw = lastweek(res)
open_html_table(lw)
pub = adjust_to_show(lw)
open_html_table(pub)
##
using DataFrames
filt=filter(row -> row.name != "Mazsola",lw)
##
range(colorant"red", stop=colorant"green", length=15)
p = AOM.plot_results(filt)
savefig(p,"/Users/dariosarra/Documents/House/AOM/Stats/03-May-2020.png")
##
plot_results(BoethiasChampions)
res = update_results(BoethiasChampions)
lw = lastweek(res)
pub = adjust_to_show(lw)
open_html_table(pub)
p = AOM.plot_results(lw)
##
using Plots
pow = dropmissing(res)
gd = groupby(pow,:date)
p = combine(gd,:total_power=>sum=>:clan_power)
scatter(p[:,:date],p[:,:clan_power])


text:
module Nektar
using Mesh2d

type NektarBndry <: Mesh2d.Bndry
    tag::Char
    elem::Int
    edge::Int
    params::Array{Float64,1}
    funs::Array{String, 1}
    NektarBndry(bt, el, ed, p, funs=Array(String,0)) = new(bt, el, ed, p, funs)
end

function section(flines, header)

    nl = length(flines)
    for i in 1:nl
        if length(search(flines[i], header)) > 0
            return i
        end
    end

    return -1
end

function element(rea, idx)
    if length(search(uppercase(rea[1]), "QUAD")) > 0
        el = Mesh2d.Quad()
    else
        el = Mesh2d.Tri()
    end
    el.id = idx
    el.verts[:,1] = [parsefloat(x) for x = split(strip(rea[2]))][1:nverts(el)]
    el.verts[:,2] = [parsefloat(x) for x = split(strip(rea[3]))][1:nverts(el)]
    return el
end

function mesh(rea)

    isec = section(rea, "MESH DATA") + 1
    nelem = parseint(split(strip(rea[isec]), " ")[1])
    elems = Array(Mesh2d.Element2d, nelem)
    isec += 1
    for e = 1:nelem
        elems[e] = element(rea[isec:(isec+2)], e)
        isec += 3
    end

    return elems

end

function bndry_sec(rea, isec, nterms=2)

    s = split(strip(rea[isec]))
    tag = strip(s[1])[1]
    nums = [parsefloat(strip(ss)) for ss in s[2:end]]
    el = convert(Int, round(nums[1]))
    ed = convert(Int, round(nums[2]))
    p = nums[3:end]
    if islower(tag)
        funs = rea[(isec+1):(isec+nterms)]
        isec += nterms + 1
    else
        funs = Array(String, 0)
        isec += 1
    end
    return NektarBndry(tag, el, ed, p, funs), isec
end

function bndry_conds(rea, elems)
    isec = section(rea, "BOUNDARY CONDITIONS")+2
    bndrylst = Array(Mesh2d.Bndry, 0)
    for el in elems
        ne = Mesh2d.nedges(el)
        for edge = 1:ne
            b, isec = bndry_sec(rea, isec)
            if b.tag == 'E' || b.tag == 'P'
                el.neigh[edge,1] = convert(Int, round(b.params[1]))
                el.neigh[edge,2] = convert(Int, round(b.params[2]))
            else
                push!(bndrylst, b)
                el.hasbc = true
            end
        end
        if ne == 3
            isec += 1
        end
    end
    nb = length(bndrylst)
    for i = 1:nb
        el = bndrylst[i].elem
        ed = bndrylst[i].edge
        elems[el].neigh[ed,2] = -i
    end
    return bndrylst
end

type NekCurve <: Mesh2d.Curve2d
    cname::String
    line::String
end

function curves(rea, elems)

    isec = section(rea, "CURVED SIDE DATA") + 1
    num = convert(Int, round(parsefloat(split(strip(rea[isec]))[1])))

    crv = Array(Mesh2d.Curve2d, num)
    if num == 0
        return crv
    end
    crvdict = Dict{String,Int}()
    isec += 1
    for i = 1:num
        ctype = lowercase(strip(rea[isec]))
        isec += 1
        if ctype == "circle"
            s = split(strip(rea[isec]))
            xc = parsefloat(s[1])
            yc = parsefloat(s[2])
            rad = parsefloat(s[3])
            inside = true
            if rad < 0
                rad = -rad
                inside = false
            end
            tag = strip(s[4])
            crvdict[tag] = i
            crv[i] = Mesh2d.Circle(xc, yc, rad, inside)
        elseif ctype=="str"
            crv[i] = Segment()
            tag = strip(rea[isec])
            crvdict[tag] = i
        elseif ctype=="fil"
            crv[i] = NekCurve(ctype, rea[isec])
            tag = strip(split(strip(rea[isec]))[2])
            crvdict[tag] = i
        else
            crv[i] = NekCurve(ctype, rea[isec])
        end
    end
    isec += 1
    num = convert(Int, round(parsefloat(split(strip(rea[isec]))[1])))
    isec += 1
    for i = 1:num
        s = split(strip(rea[isec]))
        iside = parseint(strip(s[1]))
        iel = parseint(strip(s[2]))
        tag = strip(s[3])
        elems[iel].curve[iside] = crvdict[tag]
        elems[iel].curved = true
        isec += 1
    end

    return crv

end

function nektar(filename)
    fd = open(filename, "r")
    rea = readlines(fd)
    close(fd)
    elems = mesh(rea)
    crv = curves(rea, elems)
    bcs = bndry_conds(rea, elems)

    return Mesh2d.Mesh(elems, crv, bcs, rea)
end

end


text:
if CUDA.functional()
    _cuda_log = CUDA.log
else
    _cuda_log = log
end

"""
    Fixed

Wrap a thing to indicate that it must not be learned.

# Fields
- `x`: Thing that must not be learned.
"""
struct Fixed
    x
end

# Make `Fixed` compatible with Flux's mechanisms of moving things to and from the GPU.
CUDA.cu(x::Fixed) = Fixed(CUDA.cu(x.x))
Flux.adapt_storage(T::Type{<:Array}, x::Fixed) = Fixed(Flux.adapt(T, x.x))

"""
    unwrap(x::Fixed)

Unwrap a thing wrapped by `Fixed`.

# Arguments
- `x::Fixed`: Wrapped thing.

# Returns
- Unwrapped thing.
"""
unwrap(x::Fixed) = x.x
unwrap(x) = x

"""
    untrack(model)

Untrack a model in Tracker.

# Arguments
- `model`: Model to untrack.

# Returns
- Untracked model.
"""
untrack(model) = fmap(x -> Tracker.data(x), model)

"""
    track(model)

Track a model with Tracker.

Any `AbstractArray{<:Number}` is tracked. Everything else is not tracked. This is similar
to how `Flux.params` operates.

# Arguments
- `model`: Model to track.

# Returns
- Tracked model.
"""
track(model) = fmap(_track, model)
_track(x) = x
_track(x::AbstractArray{<:Number}) = Tracker.param(x)

"""
    ceil_odd(x::T) where T<:Real

Ceil a number to the nearest odd integer.

# Arguments
- `x::T`: Number to ceil.

# Returns
- `Integer`: Nearest odd integer equal to or above `x`.
"""
ceil_odd(x::T) where T<:Real = Integer(ceil((x - 1) / 2) * 2 + 1)

"""
    insert_dim(x::T; pos::Integer) where T<:AA

# Arguments
- `x::T`: Array to insert dimension into.

# Keywords
- `pos::Integer`: Position of the new dimension.

# Returns
- `T`: `x` with an extra dimension at position `pos`.
"""
function insert_dim(x::T; pos::Integer) where T<:AA
    return reshape(x, size(x)[1:pos - 1]..., 1, size(x)[pos:end]...)
end

"""
    rbf(dists²::AA)

# Arguments
- `dists²::AA`: Squared distances.

# Returns
- `AA`: RBF kernel evaluated at squared distances `dists²`.
"""
rbf(dists²::AA) = exp.(-0.5f0 .* dists²)

"""
    compute_dists²(x::AA, y::AA)

Compute batched pairwise squared distances between tensors `x` and `y`. The batch
dimensions are the third to last.

# Arguments
- `x::AA`: Elements that correspond to the rows in the matrices of pairwise distances.
- `y::AA`: Elements that correspond to the columns in the matrices of pairwise distances.

# Returns:
- `AA`: Pairwise distances between and `x` and `y`.
"""
compute_dists²(x::AA, y::AA) = compute_dists²(x, y, Val(size(x, 2)))

compute_dists²(x::AA, y::AA, ::Val{1}) = (x .- batched_transpose(y)).^2

function compute_dists²(x::AA, y::AA, d::Val)
    y = batched_transpose(y)
    return sum(x.^2, dims=2) .+ sum(y.^2, dims=1) .- 2 .* batched_mul(x, y)
end

"""
    gaussian_logpdf(x::AA, μ::AA, σ::AA)

One-dimensional Gaussian log-pdf.

# Arguments
- `x::AA`: Values to evaluate log-pdf at.
- `μ::AA`: Means.
- `σ::AA`: Standard deviations.

# Returns
- `AA`: Log-pdfs at `x`.
"""
function gaussian_logpdf(x::AA, μ::AA, σ::AA)
    # We unroll the computation to avoid loop fusion, which causes GPU issues.
    logconst = 1.837877f0
    logdet = _cuda_log.(σ)
    logdet = 2 .* logdet
    z = (x .- μ) ./ σ
    quad = z .* z
    sum = logconst .+ logdet .+ quad
    return -sum ./ 2
end

"""
    gaussian_logpdf(x::AV, μ::AV, σ::AA)

Multi-dimensional Gaussian log-pdf.

# Arguments
- `x::AV`: Value to evaluate log-pdf at.
- `μ::AV`: Mean.
- `Σ::AM`: Covariance matrix.

# Returns
- `Real`: Log-pdf at `x`.
"""
gaussian_logpdf(x::AV, μ::AV, Σ::AM) =
    Tracker.track(gaussian_logpdf, x, μ, Σ)

function _gaussian_logpdf(x, μ, Σ)
    n = length(x)

    U = cholesky(Σ).U  # Upper triangular
    L = U'             # Lower triangular
    z = L \ (x .- μ)
    logconst = 1.837877f0
    # Taking the diagonal of L = U' causes indexing on GPU, which is why we equivalently
    # take the diagonal of U.
    logpdf = -(n * logconst + 2sum(_cuda_log.(diag(U))) + dot(z, z)) / 2

    return logpdf, n, L, U, z
end

gaussian_logpdf(x::CuOrVector, μ::CuOrVector, Σ::CuOrMatrix) =
    first(_gaussian_logpdf(x, μ, Σ))

@Tracker.grad function gaussian_logpdf(x, μ, Σ)
    logpdf, n, L, U, z = _gaussian_logpdf(Tracker.data.((x, μ, Σ))...)
    return logpdf, function (ȳ)
        u = U \ z
        eye = gpu(Matrix{Float32}(I, n, n))
        return ȳ .* -u, ȳ .* u, ȳ .* (u .* u' .- U \ (L \ eye)) ./ 2
    end
end

"""
    kl(μ₁::AA, σ₁::AA, μ₂::AA, σ₂::AA)

Kullback--Leibler divergence between one-dimensional Gaussian distributions.

# Arguments
- `μ₁::AA`: Mean of `p`.
- `σ₁::AA`: Standard deviation of `p`.
- `μ₂::AA`: Mean of `q`.
- `σ₂::AA`: Standard deviation of `q`.

# Returns
- `AA`: `KL(p, q)`.
"""
function kl(μ₁::AA, σ₁::AA, μ₂::AA, σ₂::AA)
    # We unroll the computation to avoid loop fusion, which causes GPU issues.
    logdet = _cuda_log.(σ₂ ./ σ₁)
    logdet = 2 .* logdet
    z = μ₁ .- μ₂
    σ₁², σ₂² = σ₁.^2, σ₂.^2
    quad = (σ₁² .+ z .* z) ./ σ₂²
    sum = logdet .+ quad .- 1
    return sum ./ 2
end

"""
    diagonal(x::AV)

Turn a vector `x` into a diagonal matrix.

# Arguments
- `x::AV`: Vector.

# Returns
- `AM`: Matrix with `x` on the diagonal.
"""
diagonal(x::AV) = Tracker.track(diagonal, x)

diagonal(x::Array{T, 1}) where T<:Real = convert(Array, Diagonal(x))

@Tracker.grad function diagonal(x)
    return diagonal(Tracker.data(x)), ȳ -> (diag(ȳ),)
end

"""
    batched_transpose(x)

Batch transpose tensor `x` where dimensions `1:2` are the matrix dimensions and dimension
`3:end` are the batch dimensions.

# Arguments
- `x`: Tensor to transpose.

# Returns
- Transpose of `x`.
"""
batched_transpose(x::AA) = Tracker.track(batched_transpose, x)

batched_transpose(x::CuOrArray) =
    permutedims(x, (2, 1, 3:ndims(x)...))

@Tracker.grad function batched_transpose(x)
    return batched_transpose(Tracker.data(x)), ȳ -> (batched_transpose(ȳ),)
end

"""
    batched_mul(x, y)

Batch matrix-multiply tensors `x` and `y` where dimensions `1:2` are the matrix
dimensions and dimension `3:end` are the batch dimensions.

# Args
- `x`: Left matrix in product.
- `y`: Right matrix in product.

# Returns
- Matrix product of `x` and `y`.
"""
batched_mul(x::AA, y::AA) = Tracker.track(batched_mul, x, y)

function _batched_mul(x, y)
    x, back = to_rank(3, x)
    y, _ = to_rank(3, y)
    return back(Flux.batched_mul(x, y)), x, y
end

batched_mul(x::CuOrArray, y::CuOrArray) = first(_batched_mul(x, y))

@Tracker.grad function batched_mul(x, y)
    z, x, y = _batched_mul(Tracker.data.((x, y))...)
    return z, function (ȳ)
        ȳ, back = to_rank(3, ȳ)
        return (
            back(Flux.batched_mul(ȳ, batched_transpose(y))),
            back(Flux.batched_mul(batched_transpose(x), ȳ))
        )
    end
end

"""
    logsumexp(x::AA; dims)

Safe log-sum-exp reduction of array `x` along dimensions `dims`.

# Arguments
- `x::AA`: Array to apply reductions to.
- `dims`: Dimensions along which reduction is applied.

# Returns
- `AA`: Log-sum-exp reduction of `x` along dimensions `dims`.
"""
function logsumexp(x::AA; dims=:)
    # Only do work if there is work to be done.
    !_must_work(x, dims) && (return x)
    u = maximum(Tracker.data(x), dims=dims)  # Do not track the maximum!
    # We unroll the computation to avoid loop fusion, which causes GPU issues.
    z = sum(exp.(x .- u), dims=dims)
    z = _cuda_log.(z)
    return u .+ z
end

function _must_work(x, dims)
    if dims == Colon()
        # We operate on all dimensions.
        return length(x) > 1
    else
        # We operate on a subset of the dimensions.
        return any([size(x, d) > 1 for d in dims])
    end
end

"""
    softmax(x::AA; dims)

Safe softmax array `x` along dimensions `dims`.

# Arguments
- `x::AA`: Array to apply softmax to.
- `dims`: Dimensions along which the softmax is applied.

# Returns
- `AA`: Softmax of `x` along dimensions `dims`.
"""
function softmax(x::AA; dims=:)
    u = maximum(Tracker.data(x), dims=dims)  # Do not track the maximum!
    x = exp.(x .- u)
    return x ./ sum(x, dims=dims)
end

"""
    softplus(x::AA)

Safe softplus.

# Arguments
- `x::AA`: Array to apply softplus to.

# Returns
- `AA`: Softplus applied to every element in `x`.
"""
function softplus(x::AA)
    # We unroll the computation to avoid loop fusion, which causes GPU issues.
    z = 1 .+ exp.(-abs.(x))
    z = _cuda_log.(z)
    return z .+ max.(x, 0)
end

"""
    repeat_cat(xs::MaybeAA...; dims::Integer)

Repeat the tensors `xs` appropriately many times to concatenate them along dimension `dims`.

# Arguments
- `xs::MaybeAA...`: Tensors to concatenate.
- `dims::Integer`: Dimensions to concatenate along.

# Returns
- `AA`: Concatenation of `xs`.
"""
repeat_cat(xs::AA...; dims::Integer) = cat(_repeat(xs...; dims=dims)..., dims=dims)
repeat_cat(xs::MaybeAA...; dims::Integer) =
    repeat_cat(_filter_nothings(xs)...; dims=dims)
repeat_cat(x::AA; dims::Integer) = x

"""
    repeat_merge(xs::MaybeAA...; dims::Integer)

Repeat the tensors `xs` appropriately many times to merge them along dimension `dims`.

# Arguments
- `xs::MaybeAA...`: Tensors to merge.
- `dims::Integer`: Dimensions to merge along.

# Returns
- `AA`: Merge of `xs`.
"""
function repeat_merge(xs::AA...; dims::Integer)
    xs = _repeat(xs...; dims=dims)

    # We can only merge if all inputs are equal.
    !all([xs[1] ≈ x for x in xs[2:end]]) && error("Not all inputs are approximately equal.")

    # Since all inputs are equal, their merge is simply any, so take the first one.
    return xs[1]
end
repeat_merge(xs::MaybeAA...; dims::Integer) =
    repeat_merge(_filter_nothings(xs)...; dims=dims)
repeat_merge(x::AA; dims::Integer) = x

_filter_nothings(xs) = filter(x -> !isnothing(x), collect(xs))

function _repeat(xs::MaybeAA...; dims::Integer)
    # Determine the maximum rank.
    max_rank = maximum(ndims.(xs))

    # Get the sizes of the inputs, extending to the maximum rank.
    sizes = [(size(x)..., ntuple(_ -> 1, max_rank - ndims(x))...) for x in xs]

    # Determine the maximum size of each dimension.
    max_size = maximum.(zip(sizes...))

    # Determine the repetitions for each input.
    reps = [div.(max_size, s) for s in sizes]

    # Do not repeat along the concatenation dimension.
    for i in 1:length(reps)
        reps[i][dims] = 1
    end

    # Repeat every element appropriately many times and return.
    return [repeat_gpu(x, r...) for (x, r) in zip(xs, reps)]
end

"""
    repeat_gpu(x::AA, reps::Integer...)

Repeat that is compatible with the GPU. This is not a full substitute of `repeat`: it can
only repeat along dimensions of size one.

# Arguments
- `x::AA`: Tensor to repeat.
- `reps::Integer...`: Repetitions.

# Returns
- `AA`: `x` repeated `reps` times along every dimension.
"""
function repeat_gpu(x::AA, reps::Integer...)
    # Only do work if there is work to be done.
    all(reps .== 1) && ndims(x) >= length(reps) && (return x)
    all(reps .== 1) && (return reshape(
        x, size(x)..., ntuple(_ -> 1, length(reps) - ndims(x))...
    ))
    return x .* ones_gpu(Float32, reps...)
end

"""
    expand_gpu(x::AV)

Expand a vector to a three-tensor and move it to the GPU.

# Arguments
- `x::AV`: Vector to expand.

# Returns
- `AA`: `x` as three-tensor and on the GPU.
"""
expand_gpu(x::AV) = reshape(x, :, 1, 1) |> gpu

"""
    slice_at(x::AA, dim::Integer, slice)

Slice `x` at dimension `dim` with slice `slice`.

# Arguments
- `x::AA`: Tensor to slice.
- `dim::Integer`: Dimension to slice.
- `slice`: Slice to get.

# Returns
- `AA`: `x[..., slice, ...]` with `slice` at position `dim`.
"""
slice_at(x::AA, dim::Integer, slice) =
    getindex(x, ntuple(i -> i == dim ? slice : Colon(), ndims(x))...)

"""
    split(x, dim::Integer)
    split(x)

Split `x` into two.

# Arguments
- `x`: Thing to split into two.
- `dim::Integer`: Dimension to split along.

# Returns
- `Tuple`: Two halves of `x`.
"""
function split(x, dim::Integer)
    mod(size(x, dim), 2) == 0 || error("Size of dimension $dim must be even.")
    i = div(size(x, dim), 2)  # Determine index at which to split.
    return slice_at(x, dim, 1:i), slice_at(x, dim, i + 1:size(x, dim))
end

function split(x)
    mod(length(x), 2) == 0 || error("Length of input must be even")
    i = div(length(x), 2)  # Determine index at which to split.
    return x[1:i], x[i + 1:end]
end

"""
    split_μ_σ(channels)

Split a three-tensor into means and standard deviations on dimension two.

# Arguments
- `channels`: Three-tensor to split into means and standard deviations on dimension two.

# Returns
- `Tuple{AA, AA}`: Tuple containing means and standard deviations.
"""
function split_μ_σ(channels)
    μ, transformed_σ = split(channels, 2)
    return μ, softplus(transformed_σ)
end

"""
    with_dummy(f, x)
    with_dummy(f)

Insert dimension two to `x` before applying `f` and remove dimension two afterwards.

# Arguments
- `f`: Function to apply.
- `x`: Input to `f`.

# Returns
- `f(x)`.
"""
with_dummy(f, x) = dropdims(f(insert_dim(x, pos=2)), dims=2)

"""
    to_rank(rank, x::AA)

Transform `x` into a tensor of rank `rank` by compressing the dimensions `rank + 1:end`.

# Arguments
- `rank::Integer`: Desired rank.
- `x::AA`: Tensor to compress.

# Returns
- `Tuple`: Tuple containing `x` as a `rank`-tensor and a function to transform back to
    the original dimensions.
"""
function to_rank(rank::Integer, x::AA)
    # If `x` is already rank `rank`, there is nothing to be done.
    ndims(x) == rank && (return x, identity)

    # If `x` is a tensor that broadcasts to anything, do nothing.
    size(x) == (1,) && (return x, identity)

    # Reshape `x` into a `rank`-tensor.
    size_x = size(x)
    return reshape(x, size_x[1:rank - 1]..., prod(size_x[rank:end])), function (y)
        size(y) == (1,) && (return y)
        return reshape(y, size(y)[1:rank - 1]..., size_x[rank:end]...)
    end
end

"""
    second(x)

Get the second element of `x`. This function complements `first`.

# Arguments
- `x`: Object to get second element of.

# Returns
- Second element of `x`.
"""
second(x) = x[2]

"""
    data_eltype(x)

Get the element type of the data of `x`. This removes any tracking from Tracker.

# Arguments
- `x`: Container to get element type of.

# Returns
- Type of the elements of `x` with all tracking removed.
"""
data_eltype(x) = eltype(Tracker.data(x))


text:
using Pkg
Pkg.add("Distributions")
Pkg.add("Distributed")
Pkg.add("Genie")

max_stars_repo_name: JulienPascal/MyFirstDockerApp

max_issues_repo_name: JulienPascal/MyFirstDockerApp

max_forks_repo_name: JulienPascal/MyFirstDockerApp


text:
function abs(obj::RealFieldElement)
    return jcall(obj, "abs", Object, ())
end

function round(obj::RealFieldElement)
    return jcall(obj, "round", jlong, ())
end


text:
module Preconditioners

using JuLIP: AbstractAtoms, Preconditioner, JVecs, JVecsF, Dofs, maxdist,
               constraint, bonds, cutoff, positions
using JuLIP.Potentials: PairPotential, AnalyticPotential
using JuLIP.Constraints: project!, FixedCell
using JuLIP.ASE: chemical_symbols, rnn
using PyAMG: RugeStubenSolver

import JuLIP: update!
import Base: A_ldiv_B!, A_mul_B!

# ================ AMGPrecon =====================
# this is wrapping some machinery around PyAMG

# TODO: allow direct solver as an alternative;
#       in this case we should to a cholesky factorisation
#       on `force_update!`

"""
`AMGPrecon{T}`: a preconditioner using AMG as the main solver

`AMGPrecon` stores a field `p` which is used to determine the preconditioner
matrix via
```julia
matrix(P, at::AbstractAtoms)
```
If `p` should be updated in response to an update of `at` then
the type can overload `JuLIP.Preconditioners.update_inner!`.

## Constructor:
```
AMGPrecon(p::Any, at::AbstractAtoms; updatedist=0.3, tol=1e-7)
```
## TODO:
* provide mechanism for automatically determining updatedist
from the nearest-neighbour distance in `at`
* should also have automatic updates every 10 or so iterations
"""
type AMGPrecon{T} <: Preconditioner
   p::T
   amg::RugeStubenSolver
   oldX::JVecsF
   updatedist::Float64
   tol::Float64
   updatefreq::Int
   skippedupdates::Int
end

function AMGPrecon(p, at::AbstractAtoms; updatedist=0.3, tol=1e-7, updatefreq=10)
   # make sure we don't use this in a context it is not intended for!
   @assert isa(constraint(at), FixedCell)
   P = AMGPrecon(p, RugeStubenSolver(speye(2)), copy(positions(at)),
                     updatedist, tol, updatefreq, 0)
   return force_update!(P, at)
end

A_ldiv_B!(out::Dofs, P::AMGPrecon, x::Dofs) = A_ldiv_B!(out, P.amg, x)
A_mul_B!(out::Dofs, P::AMGPrecon, f::Dofs) = A_mul_B!(out, P.amg, f)

need_update(P::AMGPrecon, at::AbstractAtoms) =
   (P.skippedupdates > P.updatefreq) ||
   (maxdist(positions(at), P.oldX) >= P.updatedist)

update!(P::AMGPrecon, at::AbstractAtoms) =
   need_update(P, at) ? force_update!(P, at) : (P.skippedupdates += 1; P)

function force_update!(P::AMGPrecon, at::AbstractAtoms)
   # perform updates of the potential p (if needed; usually not)
   P.p = update_inner!(P.p, at)
   # construct the preconditioner matrix ...
   A = project!( constraint(at), matrix(P.p, at) )
   # and the AMG solver
   P.amg = RugeStubenSolver(A, tol=P.tol)
   # remember the atom positions
   copy!(P.oldX, positions(at))
   # and remember that we just did a full update
   P.skippedupdates = 0
   return P
end

# ============== some tools to construct preconditioners ======================

function estimate_rnn(at::AbstractAtoms)
   sym = unique(chemical_symbols(at))
   return minimum([rnn(s) for s in sym])
end

# ============== Implementation of some concrete types ======================

# default implementation of update_inner!;
# most of the time no inner update is required

"update `AMGPrecon.p`, must return updated `p`"
update_inner!(p::Any, at::AbstractAtoms) = p

"return the three linear indices associated with an atom index"
atind2lininds(i::Integer) = (i-1) * 3 + [1;2;3]

"""
build the preconditioner matrix associated with a pair potential;
this is related to but not even close to the same as the hessian matrix!
"""
function matrix(p::PairPotential, at::AbstractAtoms)
   I = Int[]; J = Int[]; Z = Float64[]
   for (i, j, r, _, _) in bonds(at, cutoff(p))
      # the next 2 lines add an identity block for each atom
      # TODO: should experiment with other matrices, e.g., R ⊗ R
      ii = atind2lininds(i); jj = atind2lininds(j)
      z = p(r)
      for (a, b) in zip(ii, jj)
         append!(I, [a; a; b; b])
         append!(J, [a; b; a; b])
         append!(Z, [z; -z; -z; z])
      end
   end
   N = 3*length(at)
   return sparse(I, J, Z, N, N) + 0.001 * speye(N)
end

"""
A variant of the `Exp` preconditioner; see

### Constructor: `Exp(at::AbstractAtoms)`

Keyword arguments:

* `A=3.0`: stiffness of potential
* `r0=nothing`: if `nothing`, then it is an estimate nn distance
* `cutoff_mult`: cut-off multiplier
* `tol, updatefrew`: AMG parameters

### Reference

      D. Packwood, J. Kermode, L. Mones, N. Bernstein, J. Woolley, N. I. M. Gould,
      C. Ortner, and G. Csanyi. A universal preconditioner for simulating condensed
      phase materials. J. Chem. Phys., 144, 2016.
"""
function Exp(at::AbstractAtoms;
             A=3.0, r0=nothing, cutoff_mult=2.2, tol=1e-7, updatefreq=10)
   if r0 == nothing
      r0 = estimate_rnn(at)
   end
   rcut = r0 * cutoff_mult
   exp_shit = exp( - A*(rcut/r0 - 1.0) )
   pot = AnalyticPotential( :(exp( - $A * (r/$r0 - 1.0)) - $exp_shit),
                            cutoff = rcut )
   return AMGPrecon(pot, at, updatedist=0.2 * r0, tol=tol, updatefreq=updatefreq)
end

end # end module Preconditioners


text:
#DECLARE PHYSICAL CONSTANTS AND FUNCTIONS

const lightSpeed = 299_792_458.0
const μ₀ = 1.25663706212e-6
const Z₀ = μ₀*lightSpeed

#amplitude-scale due to divergence
function A(t₀::Float64)::Float64
   return sqrt(Z₀)/(sqrt(4π)*lightSpeed*t₀)
end

#---------------------------------------------
#DECLARE SPECIAL FUNCTIONS

function NaNnormalize(𝛏₀::Vector{Float64})::Vector{Float64}
   return ifelse( norm(𝛏₀)==0, zero(𝛏₀), normalize(𝛏₀) )
end

function angleBetween(𝛏₀::Vector{Float64},𝛏₁::Vector{Float64})::Float64
   return acos(dot(NaNnormalize(𝛏₁),NaNnormalize(𝛏₀)))
end

function distBetween(𝛏₀::Vector{Float64}, 𝛏₁::Vector{Float64})::Float64
   return norm(𝛏₀ - 𝛏₁)
end

function TXₜ2RXₜ(𝛏₀::Vector{Float64}, t₀::Float64, 𝐩ₛ::Function)::Float64
   return t₀ + distBetween(𝛏₀, 𝐩ₛ(t₀))/lightSpeed
end

function RXₜ2TXₜ(𝛏₀::Vector{Float64}, t₀::Float64, 𝐩ₛ::Function)::Float64
   return find_zero(t -> (TXₜ2RXₜ(𝛏₀, t, 𝐩ₛ) - t₀), t₀ + distBetween(𝛏₀, 𝐩ₛ(t₀)))
end


text:
## Definition of different types of links

# Link is used to define natural encoding of networks by users
# i.e. (tail,head,[capacity])
const Link = Tuple{Int,Int,Float64}
const SimpleLink = Tuple{Int,Int} # Unweighted links

function weightLink(l::SimpleLink,c::Float64)
  Link(l[1],l[2],c)
end
function weightLink(l::SimpleLink)
  weightLink(l,1.)
end

function zero_to_one!(links::Vector{Link})
  for (id,a) in enumerate(links)
    links[id] = (a[1] + 1, a[2] + 1, a[3])
  end
  return links
end
function zero_to_one!(links::Vector{SimpleLink})
  for (id,a) in enumerate(links)
    links[id] = (a[1] + 1, a[2] + 1)
  end
  return links
end

function addReverseLinks(links::Vector{Link})
  arcs = Vector{Link}()
  sizehint!(arcs, 2*length(links))
  n = 0
  for (id,a) in enumerate(links)
    push!(arcs,a)
    push!(arcs,(a[2],a[1],0.))
    n = max(n,a[1],a[2])
  end
  return arcs, n
end
function orientLinks(links::Vector{Link})
  arcs = Vector{Link}()
  sizehint!(arcs, 2*length(links))
  n = 0
  for (id,a) in enumerate(links)
    push!(arcs,a)
    push!(arcs,a[2],a[1],a[3])
    n = max(n,a[1],a[2])
  end
  return arcs, n
end

# Arc are oriented links. The network used a mixed adjacency list, i.e. nodes
# adjacency with storage of reverse arc (sym).
type Arc
  sym::Int
  head::Int
  cap::Float64
end
function Arc(h::Int,c::Float64)
  return Arc(0,h,c)
end
function ghostArc(h::Int)
  return Arc(0,h,0.)
end
function simpleArc(h::Int)
  return Arc(0,h,1.)
end


text:
module TestCore
include("preamble.jl")
using Transducers: has, reform

@testset "has" begin
    @test has(Count(), Count)
    @test has(Map(identity) |> Count(), Count)
    @test has(Count() |> Map(identity), Count)
end

@testset "reform" begin
    f = Completing(Transducers.push!!)  # bottom reducing function for eduction
    @testset for xf in [
            Count(),
            Map(identity) |> Count(),
            TeeZip(Map(identity) |> Count()),
            ]
        rf = eduction(xf, 1:1).rf
        @test rf === reform(rf, f)
        @test rf !== reform(rf, +)
    end
end

@testset "OnInit" begin
    @testset "mapfoldl" begin
        @testset for xs in iterator_variants(1:3)
            @test mapfoldl(
                Zip(Map(identity), Map(string)),
                push!,
                xs,
                init = OnInit(() -> [])
            ) == [(1, "1"), (2, "2"), (3, "3")]
        end
    end
end

end  # module


text:
using MortalityTables
using Test
using JSON
using CSV
using Pkg.Artifacts

# This is the path to the Artifacts.toml we will manipulate
MT_artifact_toml = joinpath(pkgdir(MortalityTables), "Artifacts.toml")
soa_tbl_dir = artifact_hash("mort.soa.org", MT_artifact_toml) |> artifact_path

include("CSV.jl")
include("basic.jl")
include("XTbML.jl")
include("parameterized_models.jl")
include("distribution.jl")
include("life_expectancy.jl")

# load tables to be used in subsequent tests
@test isa(MortalityTables.table(1),MortalityTables.MortalityTable)

include("projection_scale.jl")
include("MortalityTables.jl")
include("get_SOA_tables.jl")


text:
using .Networks

using Dates
using JuMP
using Statistics
using Printf

@with_kw mutable struct EnergyMarketAtFO <: AbstractMarket
    configs = EnergyMarketConfigs()
end

function aggregate_scenarios(context::AbstractContext, ech::Dates.DateTime)
    agg_scenario_name = aggregate_scenario_name(context, ech)

    agg_uncertainties_at_ech = UncertaintiesAtEch()
    uncertainties_at_ech = get_uncertainties(context, ech)
    for (injection_name, _) in uncertainties_at_ech
        for (ts, by_scenario_injections) in get_uncertainties(uncertainties_at_ech, injection_name)
            values = []
            for (_, val) in by_scenario_injections
                push!(values, val)
            end
            value_l = mean(values)
            PSCOPF.add_uncertainty!(agg_uncertainties_at_ech, injection_name, ts, agg_scenario_name, value_l)
        end
    end

    return agg_scenario_name, agg_uncertainties_at_ech
end

function run(runnable::EnergyMarketAtFO,
            ech::Dates.DateTime, firmness, TS::Vector{Dates.DateTime},
            context::AbstractContext)
    fo_start_time = TS[1] - get_fo_length(get_management_mode(context))
    if fo_start_time != ech
        msg = @sprintf("invalid step at ech=%s : EnergyMarketAtFO needs to be launched at FO start (ie %s)", ech, fo_start_time)
        throw( error(msg) )
    end

    problem_name_l = @sprintf("energy_market_at_FO_%s", ech)

    tso_actions = filter_tso_actions(get_tso_actions(context))

    gratis_starts = Set{Tuple{String,Dates.DateTime}}()
    if runnable.configs.CONSIDER_GRATIS_STARTS
        gratis_starts = init_gratis_start(context, runnable.configs.REF_SCHEDULE_TYPE)
    end

    agg_scenario_name, agg_uncertainties = aggregate_scenarios(context, ech)

    @assert(length(get_scenarios(agg_uncertainties)) == 1)
    @assert check_uncertainties(Uncertainties(ech=>agg_uncertainties), get_network(context))
    @assert check_uncertainties_contain_ts(Uncertainties(ech=>agg_uncertainties), get_target_timepoints(context))

    runnable.configs.out_path = context.out_dir
    runnable.configs.problem_name = problem_name_l

    return energy_market(get_network(context),
                        TS,
                        get_generators_initial_state(context),
                        [agg_scenario_name],
                        agg_uncertainties,
                        firmness,
                        get_market_schedule(context), #this uses original scenario names
                        get_tso_schedule(context),
                        tso_actions,
                        gratis_starts,
                        runnable.configs
                        )
end

function update_market_schedule!(context::AbstractContext, ech,
                                result::EnergyMarketModel,
                                firmness,
                                runnable::EnergyMarketAtFO)
    market_schedule = get_market_schedule(context)

    market_schedule.decider_type = DeciderType(runnable)
    market_schedule.decision_time = ech

    for ((gen_id, ts, s), p_injected_var) in result.pilotable_model.p_injected
        set_prod_definitive_value!(market_schedule, gen_id, ts, value(p_injected_var))
        if get_power_level_firmness(firmness, gen_id, ts) == DECIDED
            @assert( ismissing(get_prod_value(market_schedule, gen_id, ts))
                    || value(p_injected_var) ≈ get_prod_value(market_schedule, gen_id, ts) )
            set_prod_definitive_value!(market_schedule, gen_id, ts, value(p_injected_var))
        end
    end

    for ((gen_id, ts, _), b_on_var) in result.pilotable_model.b_on
        gen_state_value = parse(GeneratorState, value(b_on_var))
        set_commitment_definitive_value!(market_schedule, gen_id, ts, gen_state_value)
    end

    # TODO : may need to adapt cause result handles only one scenario
    # Capping
    update_schedule_capping!(market_schedule, context, ech, result.limitable_model, runnable.configs.CONSIDER_TSOACTIONS_LIMITATIONS)

    # Limitables levels : needs to be after capping update
    #FIXME : avoid re-aggregating scenarios
    agg_scenario_name, agg_uncertainties = aggregate_scenarios(context, ech)
    for gen in Networks.get_generators_of_type(get_network(context), Networks.LIMITABLE)
        gen_id = Networks.get_id(gen)
        for ts in get_target_timepoints(context)
            for s in get_scenarios(context)
                capped_l = safeget_capping(market_schedule, gen_id, ts, s)
                injected_l = get_uncertainties(agg_uncertainties,gen_id,ts,agg_scenario_name) - capped_l
                set_prod_value!(market_schedule, gen_id, ts, s, injected_l)
            end
        end
    end

    # loss_of_load (load-shedding)
    update_schedule_loss_of_load!(market_schedule, context, ech, result.lol_model)

    return market_schedule
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
""""""
function create_table_type_1(dbname, tablename)
    dbconn = SQLite.DB(dbname)
    SQLite.execute!(dbconn, """
        CREATE TABLE IF NOT EXISTS $tablename (
            id                  INTEGER PRIMARY KEY,
            link                TEXT NOT NULL,
            author              TEXT NOT NULL,
            author_website      TEXT NOT NULL,
            author_country      TEXT NOT NULL,
            author_city         TEXT NOT NULL,
            title               TEXT NOT NULL,
            subtitle            TEXT NOT NULL,
            content             TEXT NOT NULL,
            topic               TEXT NOT NULL,
            other               TEXT,
            method              TEXT NOT NULL,
            outcome             TEXT NOT NULL,
            importance          TEXT NOT NULL,
            target_group        TEXT NOT NULL,
            uniqueness          TEXT NOT NULL,
            on_shortlist        INTEGER,
            asessement_decision INTEGER,
            assessment_content  TEXT,
            community_winner    INTEGER,
            created_at          TEXT NOT NULL,
            modified_at         TEXT NOT NULL,
            pos_ratings         INTEGER NOT NULL,
            neg_ratings         INTEGER NOT NULL,
            num_comments        INTEGER NOT NULL,
            comments            TEXT
        );
    """)
end

""""""
function create_table_type_2(dbname, tablename)
    dbconn = SQLite.DB(dbname)
    SQLite.execute!(dbconn, """
        CREATE TABLE IF NOT EXISTS $tablename (
            id                  INTEGER PRIMARY KEY,
            ref                 TEXT NOT NULL,
            link                TEXT NOT NULL,
            title               TEXT NOT NULL,
            content             TEXT NOT NULL,
            pos_ratings         INTEGER NOT NULL,
            neg_ratings         INTEGER NOT NULL,
            category            TEXT NOT NULL,
            features            TEXT,
            num_comments        INTEGER NOT NULL,
            long                REAL NOT NULL,
            lat                 REAL NOT NULL,
            location            TEXT,
            mod_response        TEXT,
            assessment_content  TEXT,
            note                TEXT,
            author              TEXT NOT NULL,
            created_at          TEXT NOT NULL
        );
    """)
end

""""""
function create_table_type_3(dbname, tablename)
    dbconn = SQLite.DB(dbname)
    SQLite.execute!(dbconn, """
        CREATE TABLE IF NOT EXISTS $tablename (
            id          INTEGER PRIMARY KEY,
            number      INTEGER NOT NULL,
            content     TEXT NOT NULL,
            created_at  TEXT NOT NULL,
            author      TEXT NOT NULL,
            link        TEXT NOT NULL,
            pos_ratings INTEGER NOT NULL,
            neg_ratings INTEGER NOT NULL,
            reply_to    INTEGER
        );
    """)
end

""""""
function load(file, dbname, tablename)
    df = DataFrame(XLSX.readtable(file, "Sheet1")...)
    insertcols!(df, 1, id=1:nrow(df))
    db = SQLite.DB(dbname)
    SQLite.load!(df, db, tablename)
end

""""""
function advocate_europe(fileroot)
    dbname = "liqd_advocate_europe.sqlite"
    tablename = "contribution"
    create_table_type_1(dbname, tablename)
    load(fileroot * "/advocate_europe.xlsx", dbname, tablename)
end

""""""
function blankenburger_sueden(fileroot)
    dbname = "liqd_blankenburger_sueden.sqlite"
    tablename = "comment"
    create_table_type_3(dbname, "$(tablename)_a")
    create_table_type_3(dbname, "$(tablename)_b")
    create_table_type_3(dbname, "$(tablename)_c")
    create_table_type_3(dbname, "$(tablename)_lob_und_kritik")
    create_table_type_3(dbname, "$(tablename)_uebergreifende_themen")
    load("$fileroot/blankenburger_sueden_alternative_a.xlsx", dbname, "$(tablename)_a")
    load("$fileroot/blankenburger_sueden_alternative_b.xlsx", dbname, "$(tablename)_b")
    load("$fileroot/blankenburger_sueden_alternative_c.xlsx", dbname, "$(tablename)_c")
    load("$fileroot/blankenburger_sueden_lob_und_kritik.xlsx", dbname, "$(tablename)_lob_und_kritik")
    load("$fileroot/blankenburger_sueden_uebergreifende_themen.xlsx", dbname, "$(tablename)_uebergreifende_themen")
end

""""""
function laermorte(fileroot)
    dbname = "liqd_laermorte_melden.sqlite"
    create_table_type_2(dbname, "contribution")
    create_table_type_3(dbname, "comment")
    load("$fileroot/laermorte_melden.xlsx", dbname, "contribution")
    load("$fileroot/laermorte_melden_comments.xlsx", dbname, "comment")
end

""""""
function mauerpark(fileroot)
    dbname = "liqd_mauerpark.sqlite"
    create_table_type_2(dbname, "contribution")
    create_table_type_3(dbname, "comment")
    load("$fileroot/mauerpark.xlsx", dbname, "contribution")
    load("$fileroot/mauerpark_comments.xlsx", dbname, "comment")
end

""""""
function liqd(fileroot)
    advocate_europe(fileroot)
    blankenburger_sueden(fileroot)
    laermorte(fileroot)
    mauerpark(fileroot)
end


text:
using StatsBase

"""
    connectlocalizations!(connectID::Vector{Int64},
        clusterdata::Vector{Matrix{Float32}},
        params::ParamStruct, nframes::Int64)

Connect localizations in `clusterdata` by solving a linear assignment problem.

# Description
Connect localizations in `clusterdata` by giving "connected" localizations the
same integer value for their field `smd.connectID`.  Associations are made by
solving a linear assignment problem which designates the ideal connections
between localizations.
"""
function connectlocalizations!(connectID::Vector{Int64},
    clusterdata::Vector{Matrix{Float32}}, params::ParamStruct, nframes::Int64)

    # Loop through the clusters and solve the LAP.
    maxconnectID = maximum(connectID)
    idstocheck = findall(size.(clusterdata, 1) .> 1)
    for nn in idstocheck
        # Construct the cost matrix.
        costmatrix = create_costmatrix(clusterdata, params, nn, nframes)

        # Solve the LAP.
        assignment = solveLAP(costmatrix)

        # Update `connectID` to reflect the LAP solution.
        _, maxconnectID = linkclusters!(connectID, maxconnectID,
            Int64.(clusterdata[nn][:, 8]), assignment[1])
    end

    # Ensure `connectID` is contains the set of integers 1:nclusters.
    compress_connectID!(connectID)

    return connectID
end

"""
    connectID = connectlocalizations(connectID::Vector{Int64},
        clusterdata::Vector{Matrix{Float32}},
        params::ParamStruct, nframes::Int64)

Connect localizations in `clusterdata` by solving a linear assignment problem.

# Description
Connect localizations in `clusterdata` by giving "connected" localizations the
same integer value for their field `smd.connectID`.  Associations are made by
solving a linear assignment problem which designates the ideal connections
between localizations.
"""
function connectlocalizations(connectID::Vector{Int64},
    clusterdata::Vector{Matrix{Float32}}, params::ParamStruct, nframes::Int64)

    return connectlocalizations!(deepcopy(connectID), clusterdata, params, nframes)
end


text:
using ColorSchemes
include("geom_example.jl")
include("simplicialsurface.jl")

#= notes for the audiomp4:
starting beats:
0.05s : starts
[0.05s, 0.2s, 0.4s, 0.6s, 0.8s (big), 1.1s, 1.3s, 1.5s, 1.7s]
1.85s : beats end
2.1s-2.3s : RUN
2.5s : the drop starts
every 0.2s roughly a beat
first drum style: [3.3s, 5.1s, 6.9ss, 8.6s, 10.2s]
second drum snare? [2.7s, 3s, 3.2s, 4.7s, 6.1s, 6.4s, 6.6s, 8.2s, 9.5s, 9.8s, 10s]
big waaamk [2.5s, 6s, 9.6s]
10.4s : the drop ends
=#

function makescene(scene)
    #mesh = Makie.mesh!(scene, bunny, color=p1, show_axis=false, colormap=ColorSchemes.linear_bmw_5_95_c86_n256)
    #angle = -pi / 4 - pi / 10
    #quat = Makie.Quaternion((0., 0., sin(angle / 2), cos(angle / 2)))
    #Makie.rotate!(scene, quat)
    Makie.update_cam!(scene, 0.8 * Makie.Vec3f0(-1.4153199, 3.1283882, 0.1344321), Makie.Vec3f0(0., 0., -0.1))
    Makie.update!(scene)
    # 1.721643, 2.5823247, 0.6416913
    #Float32[-1.9936941, 3.4563544, 0.2008853]
    # realclose Float32[-1.7143533, 2.6653938, 0.2666605]
    # adjusted -1.4153199, 3.1283882, 0.1344321
    return scene
end

function bunnydance(scene)
    bunny = loadbunny()
    transvec = [Point3{Float32}(0., 0., 0.)]
    bunnyverts = coordinates(bunny)
    bunnyverts .+= transvec
    bs = constructsimplicialsurface(bunny)
    p1 = 3.0 .* bunnyeyepot2(0.)
    p2 = p1
    #p2 = solvepoissonproblem(bs, p1)

    #trans = Point3{Float32}(100., -100., 100.)
    #Makie.translate!(scene, trans...)
    #Makie.update!(scene)
    framerate = 60
    starttime = 0.5
    endtime = 10.5
    time = starttime:1/framerate:endtime

    numverts = length(vertices(bs))
    velocity = zeros(Float64, numverts)

    beatfreq = 2.0 * π * (1.0 / 1.6)
    beatfreqsmooth = 2.0 * π * 5.0
    forcingvalue = [cos(beatfreq * (t - 2.7)) for t in time]
    smoothingsteps = [cos(beatfreqsmooth * (t - 2.7)) for t in time]
    for i in 1:length(smoothingsteps)
        if smoothingsteps[i] < 0.0
            smoothingsteps[i] .* 0.1
        end
    end
    j = 0
    Makie.record(scene, "videos/generated/bunnydance.mp4", time; framerate = framerate) do i
        if mod(j, 10)  == 0
            @show minimum(p2)
            @show maximum(p2)
        end
        j += 1
        forcing = 1.0 .* forcingvalue[j].* bunnyeyepot2(0.)
        if i > 2.7
            p2[:], velocity[:] = solvewavestep(bs, p2, velocity, forcing, 100.0, 0.6)
            #@show forcing[1311]
            #@show norm(p2)
            #p2[:] = solveheatstep(bs, p2, 0.02)
            #p2[:], velocity[:] = solvewavestep(bs, p2, velocity, zeros(Float64, numverts), 1000.0, 0.1)
            newverts = solvecurvatureflowstep(bs, 0.00017 * smoothingsteps[j])
            updatemeshandsurface!(bunny, bs, newverts)
        end
        scene.plots[1][:mesh][] = bunny
        scene.plots[1][:color][] = p2
    end
    return scene
    # and then ran:
    # ffmpeg -i generated/bunnydance.mp4 -i audio/umbrella_run_sound_2.wav -map 0:v -map 1:a -c:v copy -shortest mixed/bunnydance_2.mp4
end

function bunnyeyepot2(t::Float64)
    eyepot = Array(sparsevec(
        # eyeleft, eyeright, max size, booty, ear left, ear right
        [1311, 1015, 1430, 680, 568, 747],
         [Float64(1 * cos(t)), Float64(1 * cos(t)), Float64(0), Float64(0), Float64(0), Float64(0)]))
    return eyepot
end

function generatesoundcurves(framerate=10)
    time = 0:1/framerate:10.5

    startingpulses = [0.1, 0.2, 0.4, 0.6, 0.8, 1.1, 1.3, 1.5, 1.7]
    minstart, maxstart = minimum(startingpulses), maximum(startingpulses)

    run = [2.1]

    mainbeat = collect(2.5:0.2:10.5)
    bigwaamk = [2.5, 6, 9.6]

    drums = [3.3, 5.1, 6.9, 8.6, 10.2]
    snaredrums = [2.7, 3, 3.2, 4.7, 6.1, 6.4, 6.6, 8.2, 9.5, 9.8, 10]

    heatstepcurve = zeros(Float64, length(time))
    for t in time
        @show t
    end

    return collect(time), heatstepcurve
    #return time, startingpulses, run, mainbeat, bigwaamk, drums, snaredrums, heatstepcurve
end

# bunny = loadbunny(); bs = constructsimplicialsurface(bunny); newverts = solvecurvatureflowstep(bs, 0.001); updatemeshandsurface!(bunny, bs, newverts); Makie.mesh(bunny, color=bunnyeyepot(0.), show_axis=false)

max_stars_repo_name: zoemcc/geometry

max_issues_repo_name: zoemcc/geometry

max_forks_repo_name: zoemcc/geometry


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

mutable struct ExpressRouteCircuitStats <: SwaggerModel
    primarybytesIn::Any # spec type: Union{ Nothing, Int64 } # spec name: primarybytesIn
    primarybytesOut::Any # spec type: Union{ Nothing, Int64 } # spec name: primarybytesOut
    secondarybytesIn::Any # spec type: Union{ Nothing, Int64 } # spec name: secondarybytesIn
    secondarybytesOut::Any # spec type: Union{ Nothing, Int64 } # spec name: secondarybytesOut

    function ExpressRouteCircuitStats(;primarybytesIn=nothing, primarybytesOut=nothing, secondarybytesIn=nothing, secondarybytesOut=nothing)
        o = new()
        validate_property(ExpressRouteCircuitStats, Symbol("primarybytesIn"), primarybytesIn)
        setfield!(o, Symbol("primarybytesIn"), primarybytesIn)
        validate_property(ExpressRouteCircuitStats, Symbol("primarybytesOut"), primarybytesOut)
        setfield!(o, Symbol("primarybytesOut"), primarybytesOut)
        validate_property(ExpressRouteCircuitStats, Symbol("secondarybytesIn"), secondarybytesIn)
        setfield!(o, Symbol("secondarybytesIn"), secondarybytesIn)
        validate_property(ExpressRouteCircuitStats, Symbol("secondarybytesOut"), secondarybytesOut)
        setfield!(o, Symbol("secondarybytesOut"), secondarybytesOut)
        o
    end
end # type ExpressRouteCircuitStats

const _property_map_ExpressRouteCircuitStats = Dict{Symbol,Symbol}(Symbol("primarybytesIn")=>Symbol("primarybytesIn"), Symbol("primarybytesOut")=>Symbol("primarybytesOut"), Symbol("secondarybytesIn")=>Symbol("secondarybytesIn"), Symbol("secondarybytesOut")=>Symbol("secondarybytesOut"))
const _property_types_ExpressRouteCircuitStats = Dict{Symbol,String}(Symbol("primarybytesIn")=>"Int64", Symbol("primarybytesOut")=>"Int64", Symbol("secondarybytesIn")=>"Int64", Symbol("secondarybytesOut")=>"Int64")
Base.propertynames(::Type{ ExpressRouteCircuitStats }) = collect(keys(_property_map_ExpressRouteCircuitStats))
Swagger.property_type(::Type{ ExpressRouteCircuitStats }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ExpressRouteCircuitStats[name]))}
Swagger.field_name(::Type{ ExpressRouteCircuitStats }, property_name::Symbol) =  _property_map_ExpressRouteCircuitStats[property_name]

function check_required(o::ExpressRouteCircuitStats)
    true
end

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


text:
## Relabel vector elements to 1:n ##

function relabel(label::Vector{T}) where {T}
    label_ = similar(label, Int)
    dict = Dict{T,Int}()
    for (i, l) = enumerate(label)
        if !haskey(dict, l)
            dict[l] = length(dict) + 1
        end
        label_[i] = dict[l]
    end
    return label_
end


text:
function cosineknn_test()
    A = load_matrix_network("bfs_example")
    S = cosineknn(A,2)

    A = speye(Int64,4)
    A[4,1] = 1
    CKNN = cosineknn(A,2)
    if ((CKNN[4,1] *2)^2 - 2) > 1e-7
        error("cosine knn failed")
    end
    return true
end


text:
# Closed-form formulas for the inner products between different kinds of components

# step x step
function GMStepStep(m::Int, M::Int, N::Int, μm::Float64, σm::Float64, μM::Float64, σM::Float64)::Float64

  return  (m*μm*μM-(M-m)*(1-μm)*μM+(N-M)*(1-μm)*(1-μM)) :: Float64

end
function GMStepStep(i::Int, j::Int, d, IT)::Float64
  	if i>j
    	return GMStepStep(j, i, IT, d.μ[STEP][j], d.σ[STEP][j], d.μ[STEP][i], d.σ[STEP][i]) :: Float64
    else
    	return GMStepStep(i, j, IT, d.μ[STEP][i], d.σ[STEP][i], d.μ[STEP][j], d.σ[STEP][j]) :: Float64
    end
    #return 0.0 :: Float64
end

# step x slope
function GMStepSlope(
    t::Int, l::Int, N::Int, μSTEP::Float64, σSTEP::Float64,
    μSLOPE::Float64, σSLOPE::Float64
    )::Float64

  if l > t
    return (
      t*μSTEP*μSLOPE - (l-t)*(1-μSTEP)*μSLOPE - (N+1)*l-μSLOPE*(N+1)
      + (N+1)*l*μSTEP + μSLOPE*(N+1)*μSTEP + (1/2)*(N+1)^2-(1/2)*N
      - (1/2)*μSTEP*(N+1)^2 + (1/2)*μSTEP*(N+1) + (l+1)*l+μSLOPE*(l+1)
      - (l+1)*l*μSTEP - μSLOPE*(l+1)*μSTEP - (1/2)*(l+1)^2+(1/2)*l
      + (1/2)*μSTEP*(l+1)^2 - (1/2)*μSTEP*(l+1)
      ) :: Float64
  else
    return (
      - (l+1)*l*μSTEP+(1/2)*t - (1/2)*N+(1/2)*(N+1)^2 + μSLOPE*(N+1)*μSTEP
      - μSLOPE*(l+1)*μSTEP - μSLOPE*(N+1) - (1/2)*μSTEP*(N+1)^2
      + (1/2)*μSTEP*(N+1) + (1/2)*μSTEP*(l+1)^2 - (1/2)*μSTEP*(l+1)
      + μSLOPE*(t+1) + l*μSTEP*μSLOPE - (1/2)*(t+1)^2 + (N+1)*l*μSTEP
      - (N+1)*l+(t+1)*l
      ) :: Float64
  end

  #return 0.0 :: Float64
end
function GMStepSlope(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMStepSlope(i, j, IT, d.μ[STEP][i], d.σ[STEP][i], d.μ[SLOPE][j], d.σ[SLOPE][j]) :: Float64
end
function GMSlopeStep(i::Int, j::Int, d, IT)::Float64
    return GMStepSlope(j, i, IT, d.μ[STEP][j], d.σ[STEP][j], d.μ[SLOPE][i], d.σ[SLOPE][i]) :: Float64
end

# step x spike
function GMStepSpike(
    i::Int, j::Int, N::Int, μSTEP::Float64, σSTEP::Float64,
    μSPIKE::Float64, σSPIKE::Float64
    )::Float64

  if j > i
    return (
      i*μSPIKE*μSTEP - (N-i)*(1-μSTEP)*μSPIKE + (1-μSTEP)*μSPIKE + (1-μSPIKE)*(1-μSTEP)
      )::Float64
  else
    return (
      i*μSPIKE*μSTEP - μSTEP*μSPIKE - μSTEP*(1-μSPIKE) - (N-i)*(1-μSTEP)*μSPIKE
      )::Float64
  end

  return 0.0 :: Float64
end
function GMStepSpike(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMStepSpike(i, j, IT, d.μ[STEP][i], d.σ[STEP][i], d.μ[SPIKE][j], d.σ[SPIKE][j])::Float64
end
function GMSpikeStep(i::Int, j::Int, d, IT)::Float64
    return GMStepSpike(j, i, IT, d.μ[STEP][j], d.σ[STEP][j], d.μ[SPIKE][i], d.σ[SPIKE][i])::Float64
end

# step x sine
function GMStepSin(
    t::Int, s::Int, N::Int, μSTEP::Float64, σSTEP::Float64,
    μSIN::Float64, σSIN::Float64,ω::Float64
    )::Float64

  return (
    μSTEP*μSIN*(t+1)-(1/2)*μSTEP*sin(ω)*cos((t+1)*ω)/(cos(ω)-1)
    +(1/2)*μSTEP*sin((t+1)*ω)-μSTEP*μSIN+(1/2)*μSTEP*sin(ω)*cos(ω)/(cos(ω)-1)
    -(1/2)*μSTEP*sin(ω)+(μSTEP*μSIN-μSIN)*(N+1)-(1/2)*(-1+μSTEP)*sin(ω)*cos((N+1)*ω)/(cos(ω)-1)
    +(-1/2+(1/2)*μSTEP)*sin((N+1)*ω)-(μSTEP*μSIN-μSIN)*(t+1)
    +(1/2)*(-1+μSTEP)*sin(ω)*cos((t+1)*ω)/(cos(ω)-1)-(-1/2+(1/2)*μSTEP)*sin((t+1)*ω)
    )::Float64

  #return GM::Float64
end
function GMStepSin(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMStepSin(i, j, IT, d.μ[STEP][i], d.σ[STEP][i], d.μ[SIN][j], d.σ[SIN][j],d.fs[j]) :: Float64
end
function GMSinStep(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMStepSin(j, i, IT, d.μ[STEP][j], d.σ[STEP][j], d.μ[SIN][i], d.σ[SIN][i],d.fs[i]) :: Float64
end

# step x cosine
function GMStepCos(
    t::Int, c::Int, N::Int, μSTEP::Float64, σSTEP::Float64,
    μCOS::Float64, σCOS::Float64,ω::Float64
    )::Float64

    return (
    μSTEP*μCOS*(t+1)+(1/2)*μSTEP*cos((t+1)*ω)
    +(1/2)*μSTEP*sin(ω)*sin((t+1)*ω)/(cos(ω)-1)
    -μSTEP*μCOS-(1/2)*μSTEP*cos(ω)-(1/2)*μSTEP*sin(ω)^2/(cos(ω)-1)
    +(μSTEP*μCOS-μCOS)*(N+1)+(-1/2+(1/2)*μSTEP)*cos((N+1)*ω)
    -(1/2)*(1-μSTEP)*sin(ω)*sin((N+1)*ω)/(cos(ω)-1)
    -(μSTEP*μCOS-μCOS)*(t+1)-(-1/2+(1/2)*μSTEP)*cos((t+1)*ω)
    +(1/2)*(1-μSTEP)*sin(ω)*sin((t+1)*ω)/(cos(ω)-1)
    )::Float64
end
function GMStepCos(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMStepCos(i, j, IT, d.μ[STEP][i], d.σ[STEP][i], d.μ[COS][j], d.σ[COS][j],d.fc[j]) :: Float64
end
function GMCosStep(i::Int, j::Int, d, IT)::Float64
    return GMStepCos(j, i, IT, d.μ[STEP][j], d.σ[STEP][j], d.μ[COS][i], d.σ[COS][i],d.fc[i]) :: Float64
end

# slope x slope
function GMSlopeSlope(m::Int, M::Int, N::Int, μm::Float64, σm::Float64, μM::Float64, σM::Float64)::Float64

  return (
    (1/6)*N - (1/6)*M + m*μM*μM - μM*(m+1)*μM - μM*(m+1)*m
    - (1/2)*(N+1)^2 + (1/3)*(N+1)^3 + (1/2)*(M+1)^2 - (1/3)*(M+1)^3
    + (1/2)*μM*(m+1)^2 - (1/2)*μM*(m+1) + (1/2)*m*(M+1)^2
    - (1/2)*m*(M+1) + (1/2)μM*(M+1)^2 - (1/2)*μM*(M+1)
    + (1/2)*M*(M+1)^2 - (1/2)*M*(M+1) - (1/2)*m*(N+1)^2 + (1/2)*m*(N+1)
    - (1/2)*μM*(N+1)^2 + (1/2)*μM*(N+1) - (1/2)*M*(N+1)^2
    + (1/2)*M*(N+1) - (1/2)*μM*(N+1)^2 - (M+1)*m*M
    - (M+1)*μM*M + (N+1)*m*M + μM*(N+1)*m + (N+1)*μM*M
    + μM*(N+1)*μM + (1/2)*μM*(N+1)
    ) :: Float64
end
function GMSlopeSlope(i::Int, j::Int, d::dataCD, IT::Int)::Float64
  	if j>i
    	return GMSlopeSlope(i, j, IT, d.μ[SLOPE][i], d.σ[SLOPE][i], d.μ[SLOPE][j], d.σ[SLOPE][j]) :: Float64
	else
		return GMSlopeSlope(j, i, IT, d.μ[SLOPE][j], d.σ[SLOPE][j], d.μ[SLOPE][i], d.σ[SLOPE][i]) :: Float64
	end
	return 0.0::Float64
end

# slope x spike
function GMSlopeSpike(
    l::Int, p::Int, N::Int, μSLOPE::Float64, σSLOPE::Float64,
                            μSPIKE::Float64, σSPIKE::Float64,
    )::Float64

  if p > l
    return (
      l*μSLOPE*μSPIKE + μSPIKE*(N+1)*l + μSPIKE*(N+1)*μSLOPE
      - (1/2)*μSPIKE*(N+1)^2 + (1/2)*μSPIKE*(N+1) - μSPIKE*(l+1)*l
      - μSPIKE*(l+1)*μSLOPE + (1/2)*μSPIKE*(l+1)^2
      - (1/2)*μSPIKE*(l+1) + (p-l-μSLOPE)*μSPIKE
      + (p-l-μSLOPE)*(1-μSPIKE)
      )::Float64
  else
    return (
      l*μSLOPE*μSPIKE + μSPIKE*(N+1)*l + μSPIKE*(N+1)*μSLOPE
      - (1/2)*μSPIKE*(N+1)^2 + (1/2)*μSPIKE*(N+1) - μSPIKE*(l+1)*l
      - μSPIKE*(l+1)*μSLOPE + (1/2)*μSPIKE*(l+1)^2
      - (1/2)*μSPIKE*(l+1) - μSLOPE*μSPIKE - μSLOPE*(1-μSPIKE)
      )::Float64
  end

  return 0.0::Float64
end
function GMSlopeSpike(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMSlopeSpike(i, j, IT, d.μ[SLOPE][i], d.σ[SLOPE][i], d.μ[SPIKE][j], d.σ[SPIKE][j]) :: Float64
end
function GMSpikeSlope(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMSlopeSpike(j, i, IT, d.μ[SLOPE][j], d.σ[SLOPE][j], d.μ[SPIKE][i], d.σ[SPIKE][i]) :: Float64
end

# slope x sine
function GMSlopeSin(
    l::Int, s::Int, N::Int, μSLOPE::Float64, σSLOPE::Float64,
    μSIN::Float64, σSIN::Float64,ω::Float64
    )::Float64

  return (
    μSLOPE*μSIN*(l+1)-(1/2)*μSLOPE*sin(ω)*cos((l+1)*ω)/(cos(ω)-1)
    +(1/2)*μSLOPE*sin((l+1)*ω)-μSLOPE*μSIN+(1/2)*μSLOPE*sin(ω)*cos(ω)/(cos(ω)-1)
    -(1/2)*μSLOPE*sin(ω)+(l*μSIN+(1/2)*μSIN+μSLOPE*μSIN)*(N+1)
    -(1/2)*μSIN*(N+1)^2+(1/2)*sin(ω)*(N+1)*cos((N+1)*ω)/(cos(ω)-1)-(1/2*(N+1))*sin((N+1)*ω)
    -(1/2)*(μSLOPE+l)*sin(ω)*cos((N+1)*ω)/(cos(ω)-1)
    +(1/2)*(l*cos(ω)+μSLOPE*cos(ω)-l-μSLOPE-1)*sin((N+1)*ω)/(cos(ω)-1)
    -(l*μSIN+(1/2)*μSIN+μSLOPE*μSIN)*(l+1)
    +(1/2)*μSIN*(l+1)^2-(1/2)*sin(ω)*(l+1)*cos((l+1)*ω)/(cos(ω)-1)
    +(1/2*(l+1))*sin((l+1)*ω)
    +(1/2)*(μSLOPE+l)*sin(ω)*cos((l+1)*ω)/(cos(ω)-1)
    -(1/2)*(l*cos(ω)+μSLOPE*cos(ω)-l-μSLOPE-1)*sin((l+1)*ω)/(cos(ω)-1)
    )::Float64
end
function GMSlopeSin(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMSlopeSin(i, j, IT, d.μ[SLOPE][i], d.σ[SLOPE][i], d.μ[SIN][j], d.σ[SIN][j], d.fs[j]) :: Float64
end
function GMSinSlope(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMSlopeSin(j, i, IT, d.μ[SLOPE][j], d.σ[SLOPE][j], d.μ[SIN][i], d.σ[SIN][i], d.fs[i]) :: Float64
end

# slope x cosine
@inline function GMSlopeCos(
    l::Int, c::Int, N::Int, μSLOPE::Float64, σSLOPE::Float64,
    μCOS::Float64, σCOS::Float64, ω::Float64
    )::Float64

  return (
    μSLOPE*μCOS*(l+1)+(1/2)*μSLOPE*cos((l+1)*ω)
    +(1/2)*μSLOPE*sin(ω)*sin((l+1)*ω)/(cos(ω)-1)
    -μSLOPE*μCOS-(1/2)*μSLOPE*cos(ω)-(1/2)*μSLOPE*sin(ω)^2/(cos(ω)-1)
    +(μSLOPE*μCOS+l*μCOS+(1/2)*μCOS)*(N+1)
    -(1/2)*μCOS*(N+1)^2-(1/2*(N+1))*cos((N+1)*ω)
    -(1/2)*sin(ω)*(N+1)*sin((N+1)*ω)/(cos(ω)-1)
    +(1/2)*(l*cos(ω)+μSLOPE*cos(ω)-l-μSLOPE-1)*cos((N+1)*ω)/(cos(ω)-1)
    +(1/2)*(l+μSLOPE)*sin(ω)*sin((N+1)*ω)/(cos(ω)-1)
    -(μSLOPE*μCOS+l*μCOS+(1/2)*μCOS)*(l+1)
    +(1/2)*μCOS*(l+1)^2+(1/2*(l+1))*cos((l+1)*ω)
    +(1/2)*sin(ω)*(l+1)*sin((l+1)*ω)/(cos(ω)-1)
    -(1/2)*(l*cos(ω)+μSLOPE*cos(ω)-l-μSLOPE-1)*cos((l+1)*ω)/(cos(ω)-1)
    -(1/2)*(l+μSLOPE)*sin(ω)*sin((l+1)*ω)/(cos(ω)-1)
    )::Float64
end
function GMSlopeCos(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMSlopeCos(i, j, IT, d.μ[SLOPE][i], d.σ[SLOPE][i], d.μ[COS][j], d.σ[COS][j], d.fc[j])::Float64
end
function GMCosSlope(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMSlopeCos(j, i, IT, d.μ[SLOPE][j], d.σ[SLOPE][j], d.μ[COS][i], d.σ[COS][i], d.fc[i])::Float64
end

# spike x spike
@inline function GMSpikeSpike(i::Int, j::Int, N::Int, μSPIKEi::Float64, σSPIKEi::Float64)::Float64

  if i != j
    return (
      (N-2)*μSPIKEi^2 - 2*μSPIKEi*(1-μSPIKEi)
    )::Float64
  else
    return float(N)::Float64
  end

  return 0.0::Float64
end
function GMSpikeSpike(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMSpikeSpike(i, j, IT, d.μ[SPIKE][i], d.σ[SPIKE][i])::Float64
end

# spike x sine
@inline function GMSpikeSin(
    p::Int, s::Int, N::Int, μSPIKE::Float64, σSPIKE::Float64,
    μSIN::Float64, σSIN::Float64, ω::Float64
    )::Float64

   return (
    μSPIKE*μSIN*(N+1)-(1/2)*μSPIKE*sin(ω)*cos((N+1)*ω)/(cos(ω)-1)
    +(1/2)*μSPIKE*sin((N+1)*ω)-μSPIKE*μSIN+(1/2)*μSPIKE*sin(ω)*cos(f[s])/(cos(f[s])-1)
    -(1/2)*μSPIKE*sin(ω)+μSPIKE*(sin(p*ω)-μSIN)+(1-μSPIKE)*(sin(p*ω)-μSIN)
    ) :: Float64
end
function GMSpikeSin(i::Int, j::Int, d::dataCD, IT::Int)
    return GMSpikeSin(i, j, IT, d.μ[SPIKE][i], d.σ[SPIKE][i], d.μ[SIN][j], d.σ[SIN][j], d.fs[j]) :: Float64
end
function GMSinSpike(i::Int, j::Int, d::dataCD, IT::Int)
    return GMSpikeSin(j, i, IT, d.μ[SPIKE][j], d.σ[SPIKE][j], d.μ[SIN][i], d.σ[SIN][i], d.fs[i]) :: Float64
end

# spike x cosine
@inline function GMSpikeCos(
    p::Int, c::Int, N::Int, μSPIKE::Float64, σSPIKE::Float64,
    μCOS::Float64, σCOS::Float64, ω::Float64
    )

  return (
    μSPIKE*μCOS*(N+1)+(1/2)*μSPIKE*cos((N+1)*ω)
    +(1/2)*μSPIKE*sin(ω)*sin((N+1)*ω)/(cos(ω)-1)
    -μSPIKE*μCOS-(1/2)*μSPIKE*cos(ω)-(1/2)*μSPIKE*sin(ω)^2/(cos(ω)-1)
    +μSPIKE*(cos(p*ω)-μCOS)+(1-μSPIKE)*(cos(p*ω)-μCOS)
    )::Float64
end
function GMSpikeCos(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMSpikeCos(i, j, IT, d.μ[SPIKE][i], d.σ[SPIKE][i], d.μ[COS][j], d.σ[COS][j],d.fc[j]) :: Float64
end
function GMCosSpike(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMSpikeCos(j, i, IT, d.μ[SPIKE][j], d.σ[SPIKE][j], d.μ[COS][i], d.σ[COS][i],d.fc[i]) :: Float64
end

# sine x sine
@inline function GMSinSin(i::Int, j::Int, N::Int, μSINi::Float64, σSINi::Float64, μSINj::Float64, σSINj::Float64, ωi::Float64, ωj::Float64)::Float64

  #ATENTION HERE
  #there still may be a ptoblem from cos(ωi) = cos(ωj) with i != j
  if i == j
    return float(N)::Float64
  elseif abs( cos(ωj)-cos(ωi) ) <= 1e-8
    #cancelation that need to happen for finiteness
    return (
      μSINj*μSINi*(N+1)
      -(1/2)*sin((N+1)*ωj)*sin((N+1)*ωi)-(1/2)*μSINj*sin(ωi)*cos((N+1)*ωi)/(-1+cos(ωi))
      -(1/2)*μSINi*sin(ωj)*cos((N+1)*ωj)/(cos(ωj)-1)+(1/2)*μSINj*sin((N+1)*ωi)
      +(1/2)*sin((N+1)*ωj)*μSINi-μSINj*μSINi
      +(1/2)*sin(ωi)*sin(ωj)
      +(1/2)*μSINj*sin(ωi)*cos(ωi)/(-1+cos(ωi))+(1/2)*μSINi*sin(ωj)*cos(ωj)/(cos(ωj)-1)
      -(1/2)*μSINj*sin(ωi)-(1/2)*μSINi*sin(ωj)
      )::Float64
  else
    return (
      μSINj*μSINi*(N+1)+(1/2)*sin(ωj)*cos((N+1)*ωj)*sin((N+1)*ωi)/(cos(ωj)-cos(ωi))
      -(1/2)*sin(ωi)*sin((N+1)*ωj)*cos((N+1)*ωi)/(cos(ωj)-cos(ωi))
      -(1/2)*sin((N+1)*ωj)*sin((N+1)*ωi)-(1/2)*μSINj*sin(ωi)*cos((N+1)*ωi)/(-1+cos(ωi))
      -(1/2)*μSINi*sin(ωj)*cos((N+1)*ωj)/(cos(ωj)-1)+(1/2)*μSINj*sin((N+1)*ωi)
      +(1/2)*sin((N+1)*ωj)*μSINi-μSINj*μSINi-(1/2)*sin(ωj)*cos(ωj)*sin(ωi)/(cos(ωj)-cos(ωi))
      +(1/2)*sin(ωi)*sin(ωj)*cos(ωi)/(cos(ωj)-cos(ωi))+(1/2)*sin(ωi)*sin(ωj)
      +(1/2)*μSINj*sin(ωi)*cos(ωi)/(-1+cos(ωi))+(1/2)*μSINi*sin(ωj)*cos(ωj)/(cos(ωj)-1)
      -(1/2)*μSINj*sin(ωi)-(1/2)*μSINi*sin(ωj)
      )::Float64
  end

  return 0.0::Float64
end
function GMSinSin(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMSinSin(i, j, IT, d.μ[SIN][i], d.σ[SIN][i], d.μ[SIN][j], d.σ[SIN][j], d.fs[i], d.fs[j]) :: Float64
end

# sine x cosine
@inline function GMSinCos(
    s::Int, c::Int, N::Int, μSIN::Float64, σSIN::Float64,
    μCOS::Float64, σCOS::Float64, ωs::Float64, ωc::Float64
    )::Float64

  if abs(cos(ωc)-cos(ωs)) >= 1e-5

    return (
      μSIN*μCOS*(N+1)+(1/2)*sin(ωs)*cos((N+1)*ωs)*cos((N+1)*ωc)/(cos(ωs)-cos(ωc))
      -(1/2)*sin((N+1)*ωs)*cos((N+1)*ωc)+(1/2)*sin(ωc)*sin((N+1)*ωs)*sin((N+1)*ωc)/(cos(ωs)-cos(ωc))
      +(1/2)*μSIN*cos((N+1)*ωc)-(1/2)*μCOS*sin(ωs)*cos((N+1)*ωs)/(cos(ωs)-1)
      +(1/2)*μSIN*sin(ωc)*sin((N+1)*ωc)/(cos(ωc)-1)+(1/2)*sin((N+1)*ωs)*μCOS-μSIN*μCOS
      -(1/2)*sin(ωs)*cos(ωs)*cos(ωc)/(cos(ωs)-cos(ωc))+(1/2)*sin(ωs)*cos(ωc)
      -(1/2)*sin(ωc)^2*sin(ωs)/(cos(ωs)-cos(ωc))-(1/2)*μSIN*cos(ωc)+(1/2)*μCOS*sin(ωs)*cos(ωs)/(cos(ωs)-1)
      -(1/2)*μSIN*sin(ωc)^2/(cos(ωc)-1)-(1/2)*sin(ωs)*μCOS
      )::Float64

  else
    return (
      μSIN*μCOS*(N+1)-(1/2)*cos(ωs)*cos((N+1)*ωs)^2/sin(ωs)-(1/2)*sin((N+1)*ωs)*cos((N+1)*ωs)
      +(1/2)*(cos(ωs)*μCOS+μSIN*sin(ωs)+μCOS)*cos((N+1)*ωs)/sin(ωs)
      +(1/2)*(-cos(ωs)*μSIN+μCOS*sin(ωs)-μSIN)*sin((N+1)*ωs)/sin(ωs)-μSIN*μCOS
      +(1/2)*cos(ωs)^3/sin(ωs)+(1/2)*cos(ωs)*sin(ωs)-(1/2)*(cos(ωs)*μCOS+μSIN*sin(ωs)+μCOS)*cos(ωs)/sin(ωs)
      +(1/2)*cos(ωs)*μSIN-(1/2)*μCOS*sin(ωs)+(1/2)*μSIN
      )::Float64

  end

  return 0.0::Float64
end
function GMSinCos(i::Int, j::Int, d::dataCD, IT::Int)
    return GMSinCos(i, j, IT, d.μ[SIN][i], d.σ[SIN][i], d.μ[COS][j], d.σ[COS][j],d.fs[i],d.fc[j]) :: Float64
end
function GMCosSin(i::Int, j::Int, d::dataCD, IT::Int)
    return GMSinCos(j, i, IT, d.μ[SIN][j], d.σ[SIN][j], d.μ[COS][i], d.σ[COS][i],d.fs[j],d.fc[i]) :: Float64
end

# cosine x cosine
@inline function GMCosCos(i::Int, j::Int, N::Int, μCOSi::Float64, σCOSi::Float64, μCOSj::Float64, σCOSj::Float64, ωi::Float64, ωj::Float64)::Float64

  if i == j
    return float(N)::Float64
  elseif abs(cos(ωi)-cos(ωj) ) >= 1e-8
    return (
      μCOSj*μCOSi*(N+1)-(1/2)*cos((N+1)*ωj)*cos((N+1)*ωi)
      +(1/2)*sin(ωi)*cos((N+1)*ωj)*sin((N+1)*ωi)/(cos(ωj)-cos(ωi))
      -(1/2)*sin(ωj)*sin((N+1)*ωj)*cos((N+1)*ωi)/(cos(ωj)-cos(ωi))
      +(1/2)*μCOSj*cos((N+1)*ωi)+(1/2)*cos((N+1)*ωj)*μCOSi
      +(1/2)*μCOSj*sin(ωi)*sin((N+1)*ωi)/(-1+cos(ωi))
      +(1/2)*μCOSi*sin(ωj)*sin((N+1)*ωj)/(cos(ωj)-1)-μCOSj*μCOSi+(1/2)*cos(ωj)*cos(ωi)
      -(1/2)*sin(ωi)^2*cos(ωj)/(cos(ωj)-cos(ωi))
      +(1/2)*sin(ωj)^2*cos(ωi)/(cos(ωj)-cos(ωi))-(1/2)*μCOSj*cos(ωi)
      -(1/2)*cos(ωj)*μCOSi-(1/2)*μCOSj*sin(ωi)^2/(-1+cos(ωi))-(1/2)*μCOSi*sin(ωj)^2/(cos(ωj)-1)
      )::Float64
  else
    return (
      μCOSj*μCOSi*(N+1)-(1/2)*cos((N+1)*ωj)*cos((N+1)*ωi)
      +(1/2)*μCOSj*cos((N+1)*ωi)+(1/2)*cos((N+1)*ωj)*μCOSi
      +(1/2)*μCOSj*sin(ωi)*sin((N+1)*ωi)/(-1+cos(ωi))
      +(1/2)*μCOSi*sin(ωj)*sin((N+1)*ωj)/(cos(ωj)-1)
      -μCOSj*μCOSi+(1/2)*cos(ωj)*cos(ωi)
      -(1/2)*μCOSj*cos(ωi)-(1/2)*cos(ωj)*μCOSi
      -(1/2)*μCOSj*sin(ωi)^2/(-1+cos(ωi))
      -(1/2)*μCOSi*sin(ωj)^2/(cos(ωj)-1)
      )::Float64
  end
  return 0.0::Float64
end
function GMCosCos(i::Int, j::Int, d::dataCD, IT::Int)::Float64
    return GMCosCos(i, j, IT, d.μ[COS][i], d.σ[COS][i],d.μ[COS][j], d.σ[COS][j],d.fc[i],d.fc[j])::Float64
end

# GM = Matrix{Function}(TOTALCOMPONENTS, TOTALCOMPONENTS)
# GM[STEP, STEP]   = GMStepStep
# GM[STEP, SPIKE]  = GMStepSpike
# GM[SPIKE, STEP]  = GMSpikeStep
# GM[STEP, SLOPE]  = GMStepSlope
# GM[SLOPE, STEP]  = GMSlopeStep
# GM[STEP, SIN]    = GMStepSin
# GM[SIN, STEP]    = GMSinStep
# GM[STEP, COS]    = GMStepCos
# GM[COS, STEP]    = GMCosStep
# GM[SPIKE, SPIKE] = GMSpikeSpike
# GM[SPIKE, SLOPE] = GMSpikeSlope
# GM[SLOPE, SPIKE] = GMSlopeSpike
# GM[SPIKE, SIN]   = GMSpikeSin
# GM[SIN, SPIKE]   = GMSinSpike
# GM[SPIKE, COS]   = GMSpikeCos
# GM[COS, SPIKE]   = GMCosSpike
# GM[SLOPE, SLOPE] = GMSlopeSlope
# GM[SLOPE, SIN]   = GMSlopeSin
# GM[SIN, SLOPE]   = GMSinSlope
# GM[SLOPE, COS]   = GMSlopeCos
# GM[COS, SLOPE]   = GMCosSlope
# GM[SIN, SIN]     = GMSinSin
# GM[SIN, COS]     = GMSinCos
# GM[COS, SIN]     = GMCosSin
# GM[COS, COS]     = GMCosCos

function GM2(c1,c2,i::Int, j::Int, d::dataCD, IT::Int)::Float64
  if c1 == STEP && c2 == STEP
    return GMStepStep(i::Int, j::Int, d, IT)::Float64
  elseif c1 == STEP && c2 == SPIKE
    return GMStepSpike(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SPIKE && c2 == STEP
    return GMSpikeStep(i::Int, j::Int, d, IT)::Float64
  elseif c1 == STEP && c2 == SLOPE
    return GMStepSlope(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SLOPE && c2 == STEP
    return GMSlopeStep(i::Int, j::Int, d, IT)::Float64
  elseif c1 == STEP && c2 == SIN
    return GMStepSin(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SIN && c2 == STEP
    return GMSinStep(i::Int, j::Int, d, IT)::Float64
  elseif c1 == STEP && c2 == COS
    return GMStepCos(i::Int, j::Int, d, IT)::Float64
  elseif c1 == COS && c2 == STEP
    return GMCosStep(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SPIKE && c2 == SPIKE
    return GMSpikeSpike(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SPIKE && c2 == SLOPE
    return GMSpikeSlope(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SLOPE && c2 == SPIKE
    return GMSlopeSpike(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SPIKE && c2 == SIN
    return GMSpikeSin(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SIN && c2 == SPIKE
    return GMSinSpike(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SPIKE && c2 == COS
    return GMSpikeCos(i::Int, j::Int, d, IT)::Float64
  elseif c1 == COS && c2 == SPIKE
    return GMCosSpike(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SLOPE && c2 == SLOPE
    return GMSlopeSlope(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SLOPE && c2 == SIN
    return GMSlopeSin(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SIN && c2 == SLOPE
    return GMSinSlope(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SLOPE && c2 == COS
    return GMSlopeCos(i::Int, j::Int, d, IT)::Float64
  elseif c1 == COS && c2 == SLOPE
    return GMCosSlope(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SIN && c2 == SIN
    return GMSinSin(i::Int, j::Int, d, IT)::Float64
  elseif c1 == SIN && c2 == COS
    return GMSinCos(i::Int, j::Int, d, IT)::Float64
  elseif c1 == COS && c2 == SIN
    return GMCosSin(i::Int, j::Int, d, IT)::Float64
  elseif c1 == COS && c2 == COS
    return GMCosCos(i::Int, j::Int, d, IT)::Float64
  else
    return 0.0::Float64
  end
end

max_stars_repo_name: joaquimg/ComponentFilterCD

max_issues_repo_name: joaquimg/ComponentFilterCD

max_forks_repo_name: joaquimg/ComponentFilterCD


text:
export getAveragedTemperature, getTemperature

struct FOTemp <: TemperatureSensor
  sd::SerialDevice
end

Base.close(t::FOTemp) = close(t.sd.sp)

function FOTemp(params::Dict)
  # Here we could put more parameters into the TOML file
  foTemp = FOTemp(params["connection"])
  return foTemp
end

"""
`fotemp(portAdress::AbstractString)`

Initialize Fotemp fiber optical temperature sensor. For an overview
over the high level API call `methodswith(SerialDevice{FOTemp})`.

Low level functions such as `getModelName(ft)` will return a function hash and
a whitespace, e.g. "#40 " leading the actual answer, which would be the model
name of the device in this case.
"""
function FOTemp(portAdress::AbstractString)
	pause_ms::Int=100
	timeout_ms::Int=100
	delim_read::String="\r\n"
	delim_write::String="\r"
	baudrate::Integer = 57600
	ndatabits::Integer=8
	parity::SPParity=SP_PARITY_NONE
	nstopbits::Integer=1
	rts::SPrts=SP_RTS_OFF
	cts::SPcts=SP_CTS_IGNORE
	dtr::SPdtr=SP_DTR_OFF
	dsr::SPdsr=SP_DSR_IGNORE
	xonxoff::SPXonXoff=SP_XONXOFF_DISABLED

	sp = SerialPort(portAdress)
	open(sp)
	set_speed(sp, baudrate)
	set_frame(sp,ndatabits=ndatabits,parity=parity,nstopbits=nstopbits)
	set_flow_control(sp,rts=rts,cts=cts,dtr=dtr,dsr=dsr,xonxoff=xonxoff)

	write(sp, "?40$delim_write")
	sleep(pause_ms/1000)
	answer = readuntil(sp, '\n', timeout_ms)
	flush(sp)
	if (answer=="#40 46 54 20 43 4F 4D 50 32\r\n")
		@info "Connection to fibre optical termometer established."
		return FOTemp( SerialDevice(sp,pause_ms,timeout_ms,delim_read,delim_write) )
	elseif (answer=="*FF\r\n")
		@error "Failed to establish connection. Try again."
		return nothing
	else
		@error "Connected to the wrong Device!"
		return nothing
	end
end

"""
Returns the averaged temperature of a channel.
"""
function getAveragedTemperature(ft::FOTemp,channel::Char)
	temp = query(ft.sd, "?01 $channel")
	acq = readuntil(ft.sd.sp, '\n', ft.sd.timeout_ms)
        return parse(Float64,split(temp," ")[2]) / 10
end

"""
Returns the averaged temperature of all channel.
"""
function getAveragedTemperature(ft::FOTemp)
	temp = query(ft.sd, "?02")
	acq = readuntil(ft.sd.sp, '\n', ft.sd.timeout_ms)
        return parse(Float64,split(temp," ")[2]) / 10
end

"""
Returns the current temperature of a channel.
"""
function getTemperature(ft::FOTemp, channel::Int)
	temp = query(ft.sd, "?03 $channel")
	#acq = readuntil(ft.sd.sp, '\n', ft.sd.timeout_ms)
    return parse(Float64, split(temp," ")[3]) / 10
end

"""
Returns the current temperature of all channel.
"""
function getTemperature(ft::FOTemp)
	temp = query(ft.sd, "?04")
	#acq = readuntil(ft.sd.sp, '\n', ft.sd.timeout_ms)
    A = split(temp," ")[2:end]

	T = Float64[]
	for t in A
	 t_ = tryparse(Float64, t)
	 if t_ != nothing
    	 push!(T, t_)
	 end
	end
    return T ./ 10
end

"""
Returns the number of channels the device has.
"""
function numChannels(ft::FOTemp)
	return parse(Int, split(query(ft.sd, "?0F"), " ")[2])
end

"""
Returns the the active channels. The active and inactive channels are the
bits of the ASCII encoded hexadecimal bytes returned after the whitespace
following the function hash. I.e. `"#10 1E"` represents that channels 2, 3, 4,
and 5 are on and the remaining channels 1,6,7, and 8 are switched off.

"""
function getActiveChannels(ft::FOTemp)
	return query(ft.sd, "?10")
end

"""
Sets the active channels. The parameter `channel` encodes the switched on and
off channels. The bits of the ASCII encoded hexadecimal bytes are the channels,
with bit 0 representing channel 1 to bit 7 for channel 8.

For example `channel="1E"` will switch on channels 2, 3, 4, and 5 and switch
off the remaining channels 1,6,7, and 8.
"""
function setActiveChannel(ft::FOTemp, channel::String)
	send(ft.sd, ":10 $channel")
	return nothing
end

"""
Returns the index of the active channel.
"""
function getMeasuringChannels(ft::FOTemp)
	return query(ft.sd, "?12")
end

"""
Returns the model name of the device.
"""
function getModelName(ft::FOTemp)
	return query(ft.sd, "?40")
end

"""
Returns the model name of the device.
"""
function getSerialNumber(ft::FOTemp)
	return query(ft.sd, "?41")
end

"""
Returns the firmware version.
"""
function getFirmwareVersion(ft::FOTemp)
	return query(ft.sd, "?42")
end

"""
Returns the averaging count of a given `channel`.
"""
function getTemperatureAveraging(ft::FOTemp,channel::Char)
	return query(ft.sd, "?53 $channel")
end

"""
Sets the moving averaging count of a given `channel` to 0<`numAverages`<21.
"""
function setTemperatureAveraging(ft::FOTemp,channel::Char,numAverages::Integer)
	(numAverages<0 || numAverages>20) && throw(BoundsError)
	send(ft.sd, ":53 $channel $numAverages")
	return nothing
end


text:
using Oceananigans

include("ConvergenceTests/ConvergenceTests.jl")

using  .ConvergenceTests.ForcedFlowFixedSlip: setup_and_run_xy, setup_and_run_xz

# Run 4 simulations:
Nx = [16, 32, 64, 128]
stop_time = 0.01
h = π / maximum(Nx)

Δt = 0.01 .* h^2 # Initial guess

# Modify so simulations stop on a dime
stop_iteration = round(Int, stop_time / Δt)
Δt = stop_time / stop_iteration

# Run (x, y) and (x, z) simulations
for N in Nx
    setup_and_run_xy(Nx=N, Δt=Δt, stop_iteration=stop_iteration)
    setup_and_run_xz(Nx=N, Δt=Δt, stop_iteration=stop_iteration)
end


text:
using ForwardDiff, Roots, Plots, LinearAlgebra, Optim

function backtracking_line_search(f, ∇f, x, d, α; p=0.5, β=0.1)
    y, g = f(x), ∇f(x)
    aux = 0
    while f(x + α*d) > y + β*α*(g'⋅d)
        α *= p
        aux = α
    end
    α
    return aux
end

f(x) = ((x[1]-1.0)^2) + (0.1*(x[2]-4)^2) + sin(5*pi*x[1])+sin(5*pi*x[2])
grad_f(x) = ForwardDiff.gradient(f,x)

# Ponto inicial
xk = [-1.0,1.0]

#direcao de busca
dk = -grad_f(xk)/norm(grad_f(xk))
beta = 0.1
alpha_f = 1.0

phi(alpha) = f(xk+alpha*dk)
d_phi(alpha) = grad_f(xk+alpha*dk)'dk
l_alpha(alpha) = phi(0)+alpha*beta*d_phi(0)

#condicao de Armijo
armj = []
for alpha = 0:0.01:alpha_f
    if phi(alpha) <= l_alpha(alpha)
        push!(armj,alpha)
    end
end

#Armijo com Backtracking
armjbt = backtracking_line_search(f,grad_f,xk,dk,alpha_f)

plot(phi,0,alpha_f,lw=2,lab="phi")
plot!(l_alpha,0,alpha_f,lw=2,lab="L_phi",c=:black)
scatter!(armj,phi.(armj),lab="Armijo",ms=3,c=:red)
scatter!([armjbt],[phi(armjbt)],lab="Armijo - BT",c=:yellow)
plot!(legend=:outerbottomright)
plot!(title="beta = 0.1",xlabel = "alpha", ylabel = "phi(alpha), l_alpha(alpha")


text:
vis = Visualizer()
open(vis)

# Parameters
origin = Origin{Float64}()
pbody = Sphere(0.1, 1.0)
cbody = Sphere(0.1, 1.0)
joint0to1 = JointConstraint(Floating(origin, pbody))
joint1to2 = JointConstraint(FixedOrientation(pbody, cbody))
bodies = [pbody, cbody]
joints = [joint0to1, joint1to2]
contact_p = contact_constraint(pbody, [0.0, 0.0, 1.0],
    friction_coefficient=0.5,
    contact_origin=[0.0, 0.0, 0.0],
    contact_radius=0.1,
    contact_type=:nonlinear)
contact_c = contact_constraint(cbody, [0.0, 0.0, 1.0],
    friction_coefficient=0.5,
    contact_origin=[0.0, 0.0, 0.0],
    contact_radius=0.1,
    contact_type=:nonlinear)
contacts = [contact_p, contact_c]
mech = Mechanism(origin, bodies, joints, contacts,
            gravity=-9.81,
            timestep=0.1)
mech.bodies[1].state.x2 = [0.0, 0.25, 0.5]
mech.bodies[2].state.x2 = [0.0, -0.25, 1.0]

@show adjacency_matrix(mech.joints, mech.bodies, mech.contacts)

storage = simulate!(mech, 1.0,
    verbose=true,
    record=true)

visualize(mech, storage, vis=vis)


text:
export getSensMatVec

function getSensMatVec(v::Vector,m::Vector,pFor::FWIparam)

    # extract pointers
    M    			= pFor.Mesh
    omega 			= pFor.omega
    gamma 			= pFor.gamma
    Q     			= pFor.Sources
    P     			= pFor.Receivers
    Ainv   			= pFor.ForwardSolver
    nsrc 			= size(Q,2)
    nrec 			= size(P,2)
	n 				= prod(M.n.+1);
	batchSize 		= pFor.forwardSolveBatchSize;
	select  		= pFor.sourceSelection;

	if length(select) > 0
		nsrc = length(select);
	end

	if batchSize > nsrc
		batchSize = nsrc;
	end

	numBatches 	= ceil(Int64,nsrc/batchSize);

	# derivative of mass matrix

	## ALL AT ONCE CODE
	H = spzeros(ComplexF64,n,n);
	if isa(Ainv,ShiftedLaplacianMultigridSolver)
		H = GetHelmholtzOperator(M,m,omega, gamma, true,useSommerfeldBC);
		Ainv.helmParam = HelmholtzParam(M,gamma,m,omega,true,useSommerfeldBC);
		H = H + GetHelmholtzShiftOP(m, omega,Ainv.shift[1]);
		H = sparse(H');
		# H is actually shifted laplacian now...
	elseif isa(Ainv,JuliaSolver)
		H = GetHelmholtzOperator(M,m,omega, gamma, true,useSommerfeldBC);
	end

	Jv = zeros(ComplexF64,nrec,nsrc);
	t = ((1.0.-1im*vec(gamma./omega)).*v);
	# t = ((1.0.+1im*vec(gamma)).*v);
	for k_batch = 1:numBatches
		batchIdxs = (k_batch-1)*batchSize + 1 : min(k_batch*batchSize,nsrc);
		if pFor.useFilesForFields
			filename = getFieldsFileName(omega);
			file     = matopen(filename);
			U        = read(file,string("Ubatch_",k_batch));
			close(file);
		else
			U = pFor.Fields[:,batchIdxs];
		end
		U = t.*U;
		U,Ainv = solveLinearSystem(H,U,Ainv,0);
		Jv[:,batchIdxs] = omega^2*(P'*U);
	end

	# if isa(Ainv,ShiftedLaplacianMultigridSolver)
		# Experimental: J does not necessarily needs to clear the setup.
		# clear!(Ainv.MG);
	# end

    return vec(Jv)
end


text:
# obstacle parameters

# obstacle representation for optimization problem
nOb =  3 			# number of obstacles
vOb = [3 3 2]		# number of vertices of each obstacle, vector of dimenion nOb
vObMPC = vOb-1		# adjustment for optimizaton problem

# for optimization problem
lOb = [   [ [-20;5], [-1.3;5], [-1.3;-5]]  ,
 	  [ [1.3;-5] , [1.3;5] , [20;5] ] ,
	  [ [20;11], [-20;11]]	]		#vetices given in CLOCK-WISE direction

# for plotting
lObPlot = [   [ [-20;5], [-1.3;5], [-1.3;-5], [-20;-5], [-20;5] ]  ,
      [ [1.3;5], [20;5], [20;-5], [1.3;-5], [1.3;5] ] ,
      [ [-20;15], [20;15], [20;11], [-20,11], [-20;15] ]        ]       #vetices given in CLOCK-WISE direction
nObPlot =  3        # number of obstacles
vObPlot = [4 4 4]   # number of vertices of each obstacle, vector of dimenion nOb

# build obstacles for Hybrid A* algorithm
ox = Float64[]
oy = Float64[]

# obstacle 1
for i = -12:0.1:-1.3
    push!(ox, Float64(i))
    push!(oy, 5.0)
end
for i in -2:5
    push!(ox, -1.3)
    push!(oy, Float64(i))
end

# obstacle 2
for i in -2:5
    push!(ox, 1.3)
    push!(oy, Float64(i))
end
for i = 1.3:0.1:12
    push!(ox, Float64(i))
    push!(oy, 5.0)
end

# obstacle 3
for i = -12:12
    push!(ox, Float64(i))
    push!(oy, 11.0)
end

#           [x_lower, x_upper, -y_lower,   y_upper  ]
XYbounds =  [ -15   , 15      ,  1      ,  10       ]

max_stars_repo_name: treyfortmuller/barc

max_issues_repo_name: treyfortmuller/barc