text:
using DataDeps
using DataDepsGenerators
using Test
using ReferenceTests

@testset "UCI Air Quality" begin
    registration_code = generate(UCI(), "Air+Quality")

    @testset "Integration Test" begin
        eval(Meta.parse(registration_code))
        @test length(collect(readdir(datadep"Air Quality Data Set (UCI ML Repository)"))) > 0
    end

    @test_reference "references/UCI Air+Quality.txt" registration_code
end

@test_reference "references/UCI auto mpg.txt" generate(UCI(), "https://archive.ics.uci.edu/ml/datasets/auto+mpg")
@test_reference "references/UCI banking marketting.txt" generate(UCI(), "bank+marketing")
@test_reference "references/UCI BHP.txt" generate(UCI(), "Burst+Header+Packet+%28BHP%29+flooding+attack+on+Optical+Burst+Switching+%28OBS%29+Network", "Burst Header Packet (UCI)")

@test_reference "references/UCI Adult.txt" generate(UCI(), "Adult")

@testset "ForestFires" begin
    registration_block = generate(UCI(), "Forest+Fires")

    @test occursin("A Data Mining Approach to Predict Forest Fires using Meteorological Data", registration_block) # must get the citation info

    # All the above tests in this section (and more) are covered by the final check below
    # But more specific tests let us know which bits are broken
    @test_reference "references/UCI Forest Fires.txt" registration_block
end


text:
#### indexes

export DomainIdx, subdomains, alldomains, eachdomain, allcombinations

function get_idx(gm, en, ud)
    i_gm, i_en, i_ud = 0, 0, (0,)
    ud > 0 && (i_ud = Tuple([i for i in 1:ud]))
    en > 0 && (i_en = max(i_ud...) + 1)
    gm > 0 && en > 0 && (i_gm = i_en + 1)
    gm > 0 && !(en > 0) && (i_gm = max(i_ud...) + 1)
    return i_gm, i_en, i_ud
end

""" Constructor for DomainSet """
function DomainIdx(domain_set::DomainSet)
    i_gm, i_en, i_ud = get_idx(n_subdomains(domain_set)...)
    return DomainIdx{i_gm, i_en, i_ud}()
end

gridmean(::DomainIdx{GM, EN, UD}) where {GM, EN, UD} = GM
environment(::DomainIdx{GM, EN, UD}) where {GM, EN, UD} = EN
updraft(::DomainIdx{GM, EN, UD}) where {GM, EN, UD} = UD

has_gridmean(::DomainIdx{GM, EN, UD}) where {GM, EN, UD} = !(GM == 0)
has_environment(::DomainIdx{GM, EN, UD}) where {GM, EN, UD} = !(EN == 0)
has_updraft(::DomainIdx{GM, EN, UD}) where {GM, EN, UD} = !(UD == (0,))

# TODO: write a state-vector wrapper for the exported functions
"""
    subdomains(idx::DomainIdx)

Tuple of sub-domain indexes, including environment and updrafts
"""
subdomains(idx::DomainIdx) = (environment(idx), updraft(idx)...)

"""
    alldomains(idx::DomainIdx)

Tuple of all sub-domain indexes, including grid-mean, environment and updrafts
"""
alldomains(idx::DomainIdx) = (gridmean(idx), environment(idx), updraft(idx)...)

"""
    eachdomain(idx::DomainIdx)

Tuple of indexes including
 - grid-mean
 - environment
 - Tuple of updraft indexes
"""
eachdomain(idx::DomainIdx) = (gridmean(idx), environment(idx), updraft(idx))

"""
    allcombinations(idx::DomainIdx)

Tuple of indexes including
 - grid-mean
 - environment
 - Tuple of updraft indexes
 - Tuple of all sub-domains (environment and updrafts)
 - Tuple of all domains
"""
allcombinations(idx::DomainIdx) = (
    gridmean(idx),
    environment(idx),
    updraft(idx),
    subdomains(idx),
    alldomains(idx),
)

function DomainIdx(
    domain_set::DomainSet,
    domain_subset::DomainSubSet{BGM, BEN, BUD},
) where {BGM, BEN, BUD}
    i_gm, i_en, i_ud = get_idx(n_subdomains(domain_set, domain_subset)...)
    return DomainIdx{i_gm, i_en, i_ud}()
end

@inline get_i_state_vec(vm, a_map::AbstractArray, ϕ::Symbol, i) =
    vm[ϕ][a_map[i]]
@inline get_i_state_vec(domain_decomp::DomainDecomposition, ϕ::Symbol, i) =
    get_i_state_vec(domain_decomp.var_mapper, domain_decomp.a_map[ϕ], ϕ, i)
@inline get_i_var(a_map::AbstractArray, i) = a_map[i]
@inline get_i_var(domain_decomp::DomainDecomposition, ϕ, i) =
    get_i_var(domain_decomp.a_map[ϕ], i)

function var_suffix(vm, idx::DomainIdx, idx_ss::DomainIdx, ϕ::Symbol, i)
    if i == gridmean(idx)
        return "_gm"
    elseif i == environment(idx)
        return "_en"
    elseif i in updraft(idx)
        return "_ud_" * string(i)
    else
        throw(BoundsError(vm[ϕ], i))
    end
end

var_string(vm, idx::DomainIdx, idx_ss::DomainIdx, ϕ::Symbol, i) =
    string(ϕ) * var_suffix(vm, idx, idx_ss, ϕ, i)

function var_suffix(domain_decomp::DomainDecomposition, ϕ::Symbol, i = 0)
    i == 0 && (i = gridmean(domain_decomp))
    return var_suffix(
        domain_decomp.var_mapper,
        domain_decomp.idx,
        domain_decomp.idx_ss_per_var[ϕ],
        ϕ,
        i,
    )
end

function var_string(domain_decomp::DomainDecomposition, ϕ::Symbol, i = 0)
    i == 0 && (i = gridmean(domain_decomp))
    return var_string(
        domain_decomp.var_mapper,
        domain_decomp.idx,
        domain_decomp.idx_ss_per_var[ϕ],
        ϕ,
        i,
    )
end

ext: jl


text:
using TensorLabXD
using Documenter

makedocs(;
    modules=[TensorLabXD],
    authors="PhysicsCodesLab",
    sitename="TensorLabXD.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        mathengine = MathJax()
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => ["man/intro.md", "man/tutorial.md", "man/categories.md",
                    "man/spaces.md", "man/sectors.md", "man/tensors.md"],
        "Library" => [],
        "Index" => ["index/index.md"]
    ],
)

deploydocs(;
    repo="github.com/PhysicsCodesLab/TensorLabXD.jl.git",
    devbranch="master",
)


text:
const JsUndefined = Nothing
const JsNull = Nothing

const JsBoolean = Bool
const JsNumber = Float64
const JsString = String
const JsBigInt = BigInt
const JsSymbol = Symbol

macro wrap_get_env(name)
    func_name = Symbol(:get_, name)
    napi_name = Symbol(:napi_, func_name)
    @eval $func_name() = @napi_call $napi_name()::NapiValue
end

get_global(context=current_context()) = if isnothing(context)
    @napi_call napi_get_global()::NapiValue
else
    napi_value(context)
end
@wrap_get_env undefined
@wrap_get_env null

macro wrap_create_napi_value(name, T)
    create_func = Symbol(:napi_create_, name)
    @eval napi_value(v::$T) = @napi_call $create_func(v::$T)::NapiValue
end
@wrap_create_napi_value int32 Int32
@wrap_create_napi_value uint32 UInt32
@wrap_create_napi_value int64 Int64
@wrap_create_napi_value double Float64
# @wrap_create_napi_value bigint_int64 Int64
@wrap_create_napi_value bigint_uint64 UInt64

napi_value(::Nothing) = get_undefined()
napi_value(v::Bool) = @napi_call napi_get_boolean(v::Bool)::NapiValue
napi_value(v::Number) = napi_value(Float64(v))
napi_value(s::AbstractString) = @napi_call napi_create_string_utf8(s::Cstring, length(s)::Csize_t)::NapiValue
napi_value(s::Char) = napi_value(string(s))
napi_value(s::Symbol) = @napi_call napi_create_symbol(napi_value(string(s))::NapiValue)::NapiValue
napi_value(s::Cstring) = napi_value(unsafe_string(s))
function napi_value(v::BigInt)
    sign_bit = v < 0 ? 1 : 0
    @napi_call napi_create_bigint_words(sign_bit::Cint, v.size::Csize_t, v.d::Ptr{UInt64})::NapiValue
end

macro wrap_get_value(name, T)
    napi_func = Symbol(:napi_get_value_, name)
    @eval value(::Type{$T}, v::NapiValue) = @napi_call $napi_func(v::NapiValue)::$T
end
@wrap_get_value double Float64
@wrap_get_value bool Bool
@wrap_get_value int32 Int32
@wrap_get_value uint32 UInt32

value(::Type{T}, v::NapiValue) where T <: Number = T(value(Float64, v))
value(::Type{Int64}, v::NapiValue; is_bigint = false) = if is_bigint
    ret = Ref{Int64}()
    _lossless = @napi_call napi_get_value_bigint_int64(v::NapiValue, ret::Ptr{Int64})::Bool
    ret[]
else
    @napi_call napi_get_value_int64(v::NapiValue)::Int64
end
function value(::Type{UInt64}, v::NapiValue)
    ret = Ref{UInt64}()
    _lossless = @napi_call napi_get_value_bigint_uint64(v::NapiValue, ret::Ptr{UInt64})::Bool
    ret[]
end
function value(::Type{BigInt}, v::NapiValue)
    word_count = Ref{Csize_t}()
    @napi_call napi_get_value_bigint_words(v::NapiValue, C_NULL::Ptr{Cint}, word_count::Ptr{Csize_t}, C_NULL::Ptr{UInt64})
    sign_bit = Ref{Cint}()
    words = Vector{UInt64}(undef, word_count[])
    @napi_call napi_get_value_bigint_words(v::NapiValue, sign_bit::Ptr{Cint}, word_count::Ptr{Csize_t}, words::Ptr{UInt64})
    x = reverse(reinterpret(UInt8, words))
    hex = bytes2hex(x)
    parse(BigInt, hex, base=16) * (sign_bit[] & 1 == 0 ? 1 : -1)
end
function value(::Type{String}, v::NapiValue; is_string = false)
    if is_string
        len = @napi_call napi_get_value_string_utf8(v::NapiValue, C_NULL::Ptr{Cchar}, 0::Csize_t)::Csize_t
        buf = zeros(UInt8, len + 1)
        resize!(buf, len)
        @napi_call napi_get_value_string_utf8(v::NapiValue, buf::Ptr{Cchar}, (len+1)::Csize_t)::Csize_t
        String(buf)
    else
        @with_scope begin
            s = to_string(v)
            value(String, s; is_string = true)
        end
    end
end
value(::Type{Symbol}, v::NapiValue) = Symbol(v.description)


text:
function InertialProvider(arg0::Frame)
    return InertialProvider((Frame,), arg0)
end

function InertialProvider(arg0::Rotation)
    return InertialProvider((Rotation,), arg0)
end

function InertialProvider(arg0::Rotation, arg1::Frame)
    return InertialProvider((Rotation, Frame), arg0, arg1)
end

function get_attitude(obj::InertialProvider, arg0::FieldPVCoordinatesProvider, arg1::FieldAbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", FieldAttitude, (FieldPVCoordinatesProvider, FieldAbsoluteDate, Frame), arg0, arg1, arg2)
end

function get_attitude(obj::InertialProvider, arg0::PVCoordinatesProvider, arg1::AbsoluteDate, arg2::Frame)
    return jcall(obj, "getAttitude", Attitude, (PVCoordinatesProvider, AbsoluteDate, Frame), arg0, arg1, arg2)
end


text:
@testset "svg" begin
    @testset "draw" begin
        ASSETS_DIR = joinpath(dirname(@__FILE__), "..", "..", "assets")
        demomol = sdftomol(open(joinpath(ASSETS_DIR, "test", "demo.mol")))
        drawsvg(demomol, 200, 200)
        # dest = open(joinpath(ASSETS_DIR, "image", "demo.svg"), "w")
        # write(dest, drawsvg(demomol, 200, 200))
    end
end


text:
module PropagationWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbsolutePVCoordinates = @jimport org.orekit.utils.AbsolutePVCoordinates
const AbstractPropagator = @jimport org.orekit.propagation.AbstractPropagator
const AdditionalStateProvider = @jimport org.orekit.propagation.AdditionalStateProvider
const Attitude = @jimport org.orekit.attitudes.Attitude
const AttitudeProvider = @jimport org.orekit.attitudes.AttitudeProvider
const BoundedPropagator = @jimport org.orekit.propagation.BoundedPropagator
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const Enum = @jimport java.lang.Enum
const EventDetector = @jimport org.orekit.propagation.events.EventDetector
const Field = @jimport org.hipparchus.Field
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldAbsolutePVCoordinates = @jimport org.orekit.utils.FieldAbsolutePVCoordinates
const FieldAbstractPropagator = @jimport org.orekit.propagation.FieldAbstractPropagator
const FieldAdditionalStateProvider = @jimport org.orekit.propagation.FieldAdditionalStateProvider
const FieldAttitude = @jimport org.orekit.attitudes.FieldAttitude
const FieldBoundedPropagator = @jimport org.orekit.propagation.FieldBoundedPropagator
const FieldEventDetector = @jimport org.orekit.propagation.events.FieldEventDetector
const FieldOrbit = @jimport org.orekit.orbits.FieldOrbit
const FieldOrekitFixedStepHandler = @jimport org.orekit.propagation.sampling.FieldOrekitFixedStepHandler
const FieldOrekitStepHandler = @jimport org.orekit.propagation.sampling.FieldOrekitStepHandler
const FieldPVCoordinatesProvider = @jimport org.orekit.utils.FieldPVCoordinatesProvider
const FieldPropagator = @jimport org.orekit.propagation.FieldPropagator
const FieldSpacecraftState = @jimport org.orekit.propagation.FieldSpacecraftState
const FieldTimeInterpolable = @jimport org.orekit.time.FieldTimeInterpolable
const FieldTransform = @jimport org.orekit.frames.FieldTransform
const Frame = @jimport org.orekit.frames.Frame
const Frames = @jimport org.orekit.frames.Frames
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Map = @jimport java.util.Map
const MultiSatStepHandler = @jimport org.orekit.propagation.sampling.MultiSatStepHandler
const Object = @jimport java.lang.Object
const Optional = @jimport java.util.Optional
const Orbit = @jimport org.orekit.orbits.Orbit
const OrekitFixedStepHandler = @jimport org.orekit.propagation.sampling.OrekitFixedStepHandler
const OrekitStepHandler = @jimport org.orekit.propagation.sampling.OrekitStepHandler
const PVCoordinatesProvider = @jimport org.orekit.utils.PVCoordinatesProvider
const PropagationType = @jimport org.orekit.propagation.PropagationType
const Propagator = @jimport org.orekit.propagation.Propagator
const PropagatorsParallelizer = @jimport org.orekit.propagation.PropagatorsParallelizer
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const Stream = @jimport java.util.stream.Stream
const TimeInterpolable = @jimport org.orekit.time.TimeInterpolable
const TimeStampedFieldPVCoordinates = @jimport org.orekit.utils.TimeStampedFieldPVCoordinates
const TimeStampedPVCoordinates = @jimport org.orekit.utils.TimeStampedPVCoordinates
const Transform = @jimport org.orekit.frames.Transform

export AbstractPropagator
export AdditionalStateProvider
export BoundedPropagator
export FieldAbstractPropagator
export FieldAdditionalStateProvider
export FieldBoundedPropagator
export FieldPropagator
export FieldSpacecraftState
export PropagationType
export Propagator
export PropagatorsParallelizer
export SpacecraftState
export add_additional_state
export add_additional_state_provider
export add_event_detector
export clear_events_detectors
export compare_to
export describe_constable
export ensure_compatible_additional_states
export equals
export get_a
export get_abs_pva
export get_additional_state
export get_additional_state_providers
export get_additional_states
export get_attitude
export get_attitude_provider
export get_class
export get_date
export get_declaring_class
export get_default_law
export get_e
export get_equinoctial_ex
export get_equinoctial_ey
export get_events_detectors
export get_field
export get_frame
export get_generated_ephemeris
export get_hx
export get_hy
export get_i
export get_initial_state
export get_keplerian_mean_motion
export get_keplerian_period
export get_le
export get_lm
export get_lv
export get_managed_additional_states
export get_mass
export get_max_date
export get_min_date
export get_mode
export get_mu
export get_name
export get_orbit
export get_propagators
export get_pv_coordinates
export has_additional_state
export hash_code
export interpolate
export is_additional_state_managed
export is_orbit_defined
export name
export notify
export notify_all
export ordinal
export propagate
export reset_initial_state
export set_attitude_provider
export set_ephemeris_mode
export set_master_mode
export set_slave_mode
export shifted_by
export to_spacecraft_state
export to_string
export to_transform
export value_of
export values
export wait

include("AnalyticalWrapper/AnalyticalWrapper.jl")
include("ConversionWrapper/ConversionWrapper.jl")
include("EventsWrapper/EventsWrapper.jl")
include("IntegrationWrapper/IntegrationWrapper.jl")
include("NumericalWrapper/NumericalWrapper.jl")
include("SamplingWrapper/SamplingWrapper.jl")
include("SemianalyticalWrapper/SemianalyticalWrapper.jl")
include("abstract_propagator.jl")
include("additional_state_provider.jl")
include("bounded_propagator.jl")
include("field_abstract_propagator.jl")
include("field_additional_state_provider.jl")
include("field_bounded_propagator.jl")
include("field_propagator.jl")
include("field_spacecraft_state.jl")
include("propagation_type.jl")
include("propagator.jl")
include("propagators_parallelizer.jl")
include("spacecraft_state.jl")

end


text:
using Latexify, DataFrames
import CSV, Plots

include("../src/shared_variables.jl")

function standard_plot_setup()
  Plots.plot(
    xlabel = "iterations",
    ylabel = "normalized duality gap",
    yaxis = :log,
    linestyle = :dash,
    legend = :topright,
  )
end

function get_residual(data::DataFrame)
  #return data.kkt_error_average_iterate
  return data.average_normalized_gap
  #return data.primal_delta_norms .+ data.dual_delta_norms
end

function select_data_before_iteration_limit(
  data::DataFrame,
  iteration_limit::Int64,
)
  indicies_before_iteration_limit = findall(data.iteration .<= iteration_limit)
  return data[indicies_before_iteration_limit, :]
end

function generic_restart_plot(;
  data::DataFrame,
  color::Symbol,
  markershape::Symbol,
  label::String,
  linestyle::Symbol,
  y_min::Float64,
  iteration_limit::Int64,
)
  data = select_data_before_iteration_limit(data, iteration_limit)
  residual = get_residual(data)
  # Plot the residuals
  Plots.plot!(
    data.iteration,
    residual,
    label = "",
    ylims = (y_min, maximum(residual)),
    color = color,
    linestyle = linestyle,
  )
  # Plot the restart points
  restart_indicies = findall(data.restart_occurred)
  Plots.plot!(
    data.iteration[restart_indicies],
    residual[restart_indicies],
    color = color,
    linealpha = 0.0,
    markershape = markershape,
    label = "",
    markerstrokewidth = 0,
    markerstrokecolor = :auto,
    markersize = 5.0,
  )
  # Plot the last active set change
  no_change_indicies = findall(data.number_of_active_set_changes .> 0)
  if length(no_change_indicies) > 0
    last_active_set_change_index = no_change_indicies[end]
    Plots.plot!(
      [data.iteration[last_active_set_change_index]],
      [residual[last_active_set_change_index]],
      color = color,
      markershape = :star5,
      linealpha = 0.0,
      label = "",
      markerstrokecolor = :auto,
      markersize = 7.0,
    )
  end
  # Artificial plot for the legend.
  restart_plt = Plots.plot!(
    data.iteration[restart_indicies][1:1],
    residual[restart_indicies][1:1],
    color = color,
    markershape = markershape,
    label = label,
    markerstrokecolor = :auto,
    linestyle = linestyle,
  )
  return restart_plt
end

function plot_no_restart_and_adaptive(
  results_directory::String,
  y_min::Float64,
  iteration_limit::Int64,
)
  no_restarts_df =
    CSV.read(joinpath(results_directory, "no_restarts.csv"), DataFrame)
  no_restarts_df =
    select_data_before_iteration_limit(no_restarts_df, iteration_limit)
  residual = no_restarts_df.current_normalized_gap
  Plots.plot!(
    no_restarts_df.iteration,
    residual,
    color = :black,
    label = "No restarts",
  )

  restarts_plt = generic_restart_plot(
    data = CSV.read(
      joinpath(results_directory, "adaptive_restarts.csv"),
      DataFrame,
    ),
    color = :blue,
    markershape = :circle,
    label = "Adaptive restarts",
    linestyle = :solid,
    y_min = y_min,
    iteration_limit = iteration_limit,
  )
  return restarts_plt
end

function plot_no_restart_adaptive_and_fixed_frequency_results(
  results_directory::String,
  restart_lengths::Vector{Int64},
  y_min::Float64,
  iteration_limit::Int64,
)
  # Figure out which restart lengths to plot
  df_restart_performance = DataFrame(
    restart_length = Int64[],
    first_approx_opt_index = Float64[],
    final_function_value = Float64[],
  )
  for i in 1:length(restart_lengths)
    restart_length = restart_lengths[i]
    fixed_frequency_df = CSV.read(
      joinpath(results_directory, "restart_length$(restart_length).csv"),
      DataFrame,
    )
    residuals = get_residual(fixed_frequency_df)
    approx_opt_indicies = findall(residuals .< y_min)
    if length(approx_opt_indicies) > 0
      first_approx_opt_index = approx_opt_indicies[1]
    else
      first_approx_opt_index = Inf
    end
    final_function_value = residuals[end]
    append!(
      df_restart_performance,
      DataFrame(
        restart_length = restart_length,
        first_approx_opt_index = first_approx_opt_index,
        final_function_value = final_function_value,
      ),
    )
  end
  sort!(
    df_restart_performance,
    [:first_approx_opt_index, :final_function_value],
    rev = (false, false),
  )
  # Pick three best restart lengths
  subset_of_restart_lengths = df_restart_performance.restart_length[1:3]
  sort!(subset_of_restart_lengths)

  standard_plot_setup()

  colors = [:red, :green, :purple, :orange, :pink]
  markers = [
    :dtriangle,
    :rect,
    :diamond,
    :hexagon,
    :cross,
    :xcross,
    :utriangle,
    :rtriangle,
    :ltriangle,
    :pentagon,
    :heptagon,
    :octagon,
    :vline,
    :hline,
  ]
  for i in 1:length(subset_of_restart_lengths)
    restart_length = subset_of_restart_lengths[i]
    fixed_frequency_df = CSV.read(
      joinpath(results_directory, "restart_length$(restart_length).csv"),
      DataFrame,
    )
    generic_restart_plot(;
      data = fixed_frequency_df,
      color = colors[i],
      markershape = markers[i],
      label = "Restart length = $restart_length",
      linestyle = :dot,
      y_min = y_min,
      iteration_limit = iteration_limit,
    )
  end

  restarts_plt =
    plot_no_restart_and_adaptive(results_directory, y_min, iteration_limit)
  return restarts_plt
end

function plot_dynamic_adaptive_and_no_restarts(
  results_directory::String,
  restart_lengths::Vector{Int64},
  y_min::Float64,
  iteration_limit::Int64,
)
  standard_plot_setup()
  plot_no_restart_and_adaptive(results_directory, y_min, iteration_limit)

  dynamic_adaptive_restarts_df = CSV.read(
    joinpath(results_directory, "dynamic_adaptive_restarts.csv"),
    DataFrame,
  )
  residual = get_residual(dynamic_adaptive_restarts_df)
  flexible_restarts_plt = Plots.plot!(
    dynamic_adaptive_restarts_df.iteration,
    residual,
    label = "Flexible restarts",
    ylims = (y_min, maximum(residual)),
    color = :blue,
    linestyle = :dot,
  )
  return flexible_restarts_plt
end

function first_iteration_to_hit_tolerance(
  df::DataFrame,
  target_tolerance::Float64,
)
  indicies_below_tolerance = findall(
    min.(df.kkt_error_average_iterate, df.kkt_error_current_iterate) .<=
    target_tolerance,
  )
  if length(indicies_below_tolerance) > 0
    return df.iteration[indicies_below_tolerance[1]]
  else
    return Inf
  end
end

function create_dictionary_of_iterations_to_hit_tolerance(
  problem_name::String,
  results_directory::String,
  restart_lengths::Vector{Int64},
  target_tolerance::Float64,
)

  dictionary_hits = Dict()
  dictionary_hits["problem_name"] = problem_name
  dictionary_hits["no_restarts"] = first_iteration_to_hit_tolerance(
    CSV.read(joinpath(results_directory, "no_restarts.csv"), DataFrame),
    target_tolerance,
  )
  dictionary_hits["adaptive_restarts"] = first_iteration_to_hit_tolerance(
    CSV.read(joinpath(results_directory, "adaptive_restarts.csv"), DataFrame),
    target_tolerance,
  )
  dictionary_hits["flexible_restarts"] = first_iteration_to_hit_tolerance(
    CSV.read(
      joinpath(results_directory, "dynamic_adaptive_restarts.csv"),
      DataFrame,
    ),
    target_tolerance,
  )
  best_fixed_frequency_iterations = Inf
  for restart_length in restart_lengths
    fixed_frequency_df = CSV.read(
      joinpath(results_directory, "restart_length$(restart_length).csv"),
      DataFrame,
    )
    fixed_frequency_iterations =
      first_iteration_to_hit_tolerance(fixed_frequency_df, target_tolerance)
    if fixed_frequency_iterations < best_fixed_frequency_iterations
      best_fixed_frequency_iterations = fixed_frequency_iterations
    end
  end
  dictionary_hits["best_fixed_frequency"] = best_fixed_frequency_iterations

  return dictionary_hits
end

results_directory = ARGS[1]
@assert length(ARGS) == 1

function main()
  y_min = 1e-7
  plot_iteration_limit_dict = Dict(
    "qap10" => ITERATION_LIMIT,
    "qap15" => ITERATION_LIMIT,
    # manually make nug08-3rd shorter so the plot is interesting
    "nug08-3rd" => 10000,
    "nug20" => ITERATION_LIMIT,
  )

  for problem_name in ALL_PROBLEM_NAMES
    subdirectory = joinpath(results_directory, problem_name)
    restarts_plt = plot_no_restart_adaptive_and_fixed_frequency_results(
      subdirectory,
      RESTART_LENGTHS_DICT[problem_name],
      y_min,
      ITERATION_LIMIT,
    )
    Plots.savefig(
      restarts_plt,
      joinpath(subdirectory, "$(problem_name)-adaptive-residuals.pdf"),
    )
    flexible_plt = plot_dynamic_adaptive_and_no_restarts(
      subdirectory,
      RESTART_LENGTHS_DICT[problem_name],
      y_min,
      ITERATION_LIMIT,
    )
    Plots.savefig(
      flexible_plt,
      joinpath(subdirectory, "$(problem_name)-flexible-residuals.pdf"),
    )
  end

  ####################
  # table of results #
  ####################
  target_tolerance = 1e-6
  df_hit_tolerance = DataFrames.DataFrame(
    problem_name = String[],
    no_restarts = Float64[],
    best_fixed_frequency = Float64[],
    adaptive_restarts = Float64[],
    flexible_restarts = Float64[],
  )
  for problem_name in ALL_PROBLEM_NAMES
    subdirectory = joinpath(results_directory, problem_name)
    push!(
      df_hit_tolerance,
      create_dictionary_of_iterations_to_hit_tolerance(
        problem_name,
        subdirectory,
        RESTART_LENGTHS_DICT[problem_name],
        target_tolerance,
      ),
    )
  end

  latex_table_file = open(joinpath(results_directory, "latex_table.txt"), "w")
  write(latex_table_file, latexify(df_hit_tolerance, env = :table))
  close(latex_table_file)
end

main()

max_stars_repo_name: pedersor/google-research

max_issues_repo_name: pedersor/google-research

max_forks_repo_name: pedersor/google-research


text:
mutable struct TileStates
    quads::Array{Tuple{Integer, Integer}, 2}
    chars::Array{Char, 2}
    rands::Array{UInt8, 2}

    TileStates() = new(Matrix{Tuple{Integer, Integer}}(undef, 0, 0), Matrix{Char}(undef, 0, 0), Matrix{Integer}(undef, 0, 0))
end

mutable struct DrawableRoom
    map::Map
    room::Room

    fgTileStates::TileStates
    bgTileStates::TileStates

    rendering::Union{Layer, Nothing}
    layers::Array{Layer, 1}

    fillColor::Union{Nothing, colorTupleType}
end

Base.size(state::TileStates) = size(state.rands)

getTileStateSeed(name::String, package::String="", fg::Bool=false) = foldl((a, b) -> a + Int(b) * 128, collect(package * (fg ? "fg" : "bg") * name), init=0)
getTileStateSeed(room::Room, package::String="", fg::Bool=false) = getTileStateSeed(room.name, package, fg)

function updateTileStates!(room::String, package::String, states::TileStates, width::Integer, height::Integer, fg::Bool=false)
    seed = getTileStateSeed(room, package, fg)
    rng = MersenneTwister(seed)

    states.quads = fill((-1, -1), (height, width))
    states.chars = fill('0', (height, width))
    states.rands = rand(rng, 0:255, (height, width))
end

updateTileStates!(room::Room, package::String, states::TileStates, width::Integer, height::Integer, fg::Bool=false) = updateTileStates!(room.name, package, states, width, height, fg)

function getDrawingLayers()
    return Layer[
        Layer("bgParallax"),
        Layer("bgTiles", clearOnReset=false),
        Layer("bgDecals"),
        Layer("entities"),
        Layer("fgTiles", clearOnReset=false),
        Layer("fgParallax"),
        Layer("fgDecals"),
        Layer("triggers"),

        Layer("tools", selectable=false),
        Layer("all", dummy=true)
    ]
end

function roomVisible(camera::Camera, width::Integer, height::Integer, room::Room)
    actuallX = camera.x / camera.scale
    actuallY = camera.y / camera.scale

    actuallWidth = width / camera.scale
    actuallHeight = height / camera.scale

    cameraRect = Rectangle(actuallX, actuallY, actuallWidth, actuallHeight)
    roomRect = Rectangle(Int.(room.position)..., Int.(room.size)...)

    return checkCollision(cameraRect, roomRect)
end

redrawRenderingLayer(renderingLayer::Layer, layers::Array{Layer, 1}) = renderingLayer.redraw || any(map(layer -> layer.redraw, layers))

function getRoomBackgroundColor(room::Room)
    if 0 <= room.color < length(colors.background_room_color_coded_fill)
        return colors.background_room_color_coded_fill[room.color + 1]

    else
        return colors.background_room_fill
    end
end

function DrawableRoom(map::Map, room::Room)
    return DrawableRoom(
        map,
        room,

        TileStates(),
        TileStates(),

        Layer("rendering"),
        getDrawingLayers(),

        nothing
    )
end

function destroy(dr::DrawableRoom)
    ctx = getSurfaceContext(dr.rendering.surface)
    Cairo.destroy(ctx)

    for layer in dr.layers
        ctx = getSurfaceContext(layer.surface)
        Cairo.destroy(ctx)
    end
end


text:
import Random

sample_snapshot_id() = sample_snapshot_ids(1) |> first
sample_snapshot_ids(num_samples) = rand(get_snapshot_ids(), num_samples)

function sample_daughters(snapshot, organism_to_replace::Organism, genotype, num_samples)
    daughters = []

    while length(daughters) < num_samples
        append!(daughters, get_daughters(snapshot, organism_to_replace, genotype))
    end

    return daughters[1:num_samples]
end

function sample_fitness_values(snapshot, organism_to_replace::Organism, genotype, num_samples)
    [get_fitness(snapshot, organism_to_replace, genotype) for _ in 1:num_samples]
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
export g₀, rocket_thrust, mass_flow
# These are equations related to rocketry from Tsiolkovsky rocket equation
#  https://en.wikipedia.org/wiki/Tsiolkovsky_rocket_equation

const global g₀ = 9.80665 # Acceleration of gravity on earth in m/s²

"""
    rocket_thrust(Isp, mass_flow)

Calculate the thrust generated by a rocket engine with specific impulse of `Isp` and
mass flow of `mass_flow`. Specific impulse says something about efficiency of engine:
how much force you get out of one unit of propellant.
"""
rocket_thrust(Isp::Number, mass_flow::Number) = g₀ * Isp * mass_flow

"""
    mass_flow(thrust, Isp)
Get mass flow in Kg/s of propellant for an engine with given `thrust` and specific impulse `Isp`.

## Example

```@example
# Calculating mass flow for a Falcon 9 rocket.

engine_thrust = 0.845e6            # Newton
Isp           = 282                # Specific impulse in seconds
thrust  = engine_thrust * 9        # There are 9 merlin engines on Falcon 9
mflow1  = mass_flow(thrust, Isp)   # Total mass flow in rocket. Kg/s

# Should get similar result from using burn time and total propellant
burn_time = 162          # Number of seconds first stage engines will burn
propellant_mass = 395700 # About 400 tons
mflow2 = propellant_mass / burn_time
```
"""
mass_flow(thrust::Number, Isp::Number) = thrust / (Isp * g₀)


text:
using DataStructures

"""
    process_input(sentences, mincount, window)

Preprocessing function for input to be used in Word2Vec model.  Takes in an array of
tokenized sentences, an integer for the minimum word frequency, and the context
window size (one sided).

Returns a frequency table to use in downsampling, a vocabulary hash, a vocabulary
array, a PriorityQueue for constructing a HuffmanTree, and the data for Word2Vec
training.
"""
function process_input(sentences::Array{Array{T, 1}, 1}, mincount::Integer, window::Integer, model::Symbol) where T <: AbstractString
    pq, freq_table, vocab_hash, vocab = _dicts(_unigrams(sentences), mincount)
    data = _data(sentences, window, vocab_hash, freq_table, model)
    return vocab_hash, vocab, pq, data
end

function _unigrams(sentences::Array{Array{T, 1}, 1}) where T <: AbstractString
    unigrams = Dict{String, Integer}()
    for sentence in sentences
        for word in sentence
            in(word, keys(unigrams)) ? unigrams[word] += 1 : unigrams[word] = 1
        end
    end
    unigrams
end

# creates PriorityQueue for HuffmanTree, as well as a frequency table for downsampling
# a vocabulary hash dict, and vocabulary array
function _dicts(unigrams::Dict, mincount::Integer)
    vocab_hash = Dict{String, Int}()
    vocab = Array{String, 1}()
    count_table = _dropmin(unigrams, mincount)
    for (i,j) in enumerate(count_table)
        (k,l) = j
        vocab_hash[k] = i
        push!(vocab, k)
    end
    pq = PriorityQueue(count_table)
    vocab_size = sum(values(count_table))
    freq_table = Dict(vocab_hash[key]=>value/vocab_size for (key, value) in count_table)
    return pq, freq_table, vocab_hash, vocab
end

# drop out unigrams that don't meet the minimum count threshold and return a
# count table
function _dropmin(unigrams, mincount)
    count_table = Dict{String, Int}()
    for (i, j) in unigrams
        if j ≥ mincount
            count_table[i] = j
        else
            continue
        end
    end
    return count_table
end

# takes in array of sentences, returns a sateteful IterTools.Partition for each
# tokenized sentence in the doc
function _data(sentences::Array{Array{T, 1}, 1}, window::Integer, vocab_hash::Dict, freq_table::Dict, model::Symbol) where T <: AbstractString
    data = word2vecdata(length(keys(vocab_hash)),freq_table,  model)
    for sentence in sentences
        if isempty(sentence)
            continue
        else
            dropmin_sentence = Int[]
            for word in sentence
                if in(word, keys(vocab_hash))
                    push!(dropmin_sentence, vocab_hash[word])
                else
                    continue
                end
            end
            if isempty(dropmin_sentence) || length(dropmin_sentence) == 1
                continue
            else
                push!(data, Context(window, dropmin_sentence))
            end
        end
    end
    return data
end


text:
function rob(nums::Array{T,1}) where T
    table = Dict{T,UInt64}()
    function inner(nums, start, len, table)
        try
            return table[start]
        catch end

        if start == len
            table[start] = nums[start]
            return nums[start]
        elseif start == len + 1
            table[start] = 0
            return 0
        end

        a = max(
            maximum([inner(nums, x, len, table) for x in start+1:len]),
            nums[start] + inner(nums, start+2, len, table)
        )

        table[start] = a
        a
    end
    inner(nums, 1, length(nums), table)
end

using Test
@test rob([183, 219, 57, 193, 94, 233, 202, 154, 65, 240, 97, 234, 100, 249, 186, 66, 90, 238, 168, 128, 177, 235, 50, 81, 185, 165, 217, 207, 88, 80, 112, 78, 135, 62, 228, 247, 211]) == UInt64(3365)
@test rob([9, 297, 196, 336, 435, 2, 343, 159, 146, 359, 45, 470, 265, 131, 17, 271, 74, 242, 448, 402, 55, 423, 414, 240, 430, 135, 322, 468, 422, 351, 441, 463, 30, 399, 132, 439, 463, 260, 399, 32, 374, 383, 276, 166, 104, 315, 314, 445, 458, 422, 104, 251, 382, 230, 484, 127, 355, 332, 317, 362, 257, 493, 474, 401, 40, 93, 433, 464, 136, 342, 98, 159, 223, 110, 89, 47, 53, 179, 219, 314, 486, 301, 307, 453, 37, 366, 334, 355, 26, 484, 124, 408, 346, 133, 420, 280, 124, 210, 358, 140]) == UInt64(15632)

# use this to test time
# @time rob([9, 297, 196, 336, 435, 2, 343, 159, 146, 359, 45, 470, 265, 131, 17, 271, 74, 242, 448, 402, 55, 423, 414, 240, 430, 135, 322, 468, 422, 351, 441, 463, 30, 399, 132, 439, 463, 260, 399, 32, 374, 383, 276, 166, 104, 315, 314, 445, 458, 422, 104, 251, 382, 230, 484, 127, 355, 332, 317, 362, 257, 493, 474, 401, 40, 93, 433, 464, 136, 342, 98, 159, 223, 110, 89, 47, 53, 179, 219, 314, 486, 301, 307, 453, 37, 366, 334, 355, 26, 484, 124, 408, 346, 133, 420, 280, 124, 210, 358, 140, 9, 297, 196, 336, 435, 2, 343, 159, 146, 359, 45, 470, 265, 131, 17, 271, 74, 242, 448, 402, 55, 423, 414, 240, 430, 135, 322, 468, 422, 351, 441, 463, 30, 399, 132, 439, 463, 260, 399, 32, 374, 383, 276, 166, 104, 315, 314, 445, 458, 422, 104, 251, 382, 230, 484, 127, 355, 332, 317, 362, 257, 493, 474, 401, 40, 93, 433, 464, 136, 342, 98, 159, 223, 110, 89, 47, 53, 179, 219, 314, 486, 301, 307, 453, 37, 366, 334, 355, 26, 484, 124, 408, 346, 133, 420, 280, 124, 210, 358, 140, 9, 297, 196, 336, 435, 2, 343, 159, 146, 359, 45, 470, 265, 131, 17, 271, 74, 242, 448, 402, 55, 423, 414, 240, 430, 135, 322, 468, 422, 351, 441, 463, 30, 399, 132, 439, 463, 260, 399, 32, 374, 383, 276, 166, 104, 315, 314, 445, 458, 422, 104, 251, 382, 230, 484, 127, 355, 332, 317, 362, 257, 493, 474, 401, 40, 93, 433, 464, 136, 342, 98, 159, 223, 110, 89, 47, 53, 179, 219, 314, 486, 301, 307, 453, 37, 366, 334, 355, 26, 484, 124, 408, 346, 133, 420, 280, 124, 210, 358, 140])

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
using ArrayViews
import ArrayViews: view

typealias CVecView{T} ContiguousView{T,1,Vector{T}}

immutable SparseVectorView{Tv,Ti<:Integer} <: AbstractSparseVector{Tv,Ti}
    n::Int                  # the number of elements
    nzind::CVecView{Ti}     # the indices of nonzeros
    nzval::CVecView{Tv}     # the values of nonzeros

    function SparseVectorView(n::Integer, nzind::CVecView{Ti}, nzval::CVecView{Tv})
        n >= 0 || throw(ArgumentError("The number of elements must be non-negative."))
        length(nzind) == length(nzval) ||
            throw(DimensionMismatch("The lengths of nzind and nzval are inconsistent."))
        new(convert(Int, n), nzind, nzval)
    end
end

### Construction

SparseVectorView{Tv,Ti}(n::Integer, nzind::CVecView{Ti}, nzval::CVecView{Tv}) =
    SparseVectorView{Tv,Ti}(n, nzind, nzval)

view(x::AbstractSparseVector) =
    SparseVectorView(length(x), view(nonzeroinds(x)), view(nonzeros(x)))

### Basic properties

length(x::SparseVectorView) = x.n
size(x::SparseVectorView) = (x.n,)
nnz(x::SparseVectorView) = length(x.nzval)
countnz(x::SparseVectorView) = countnz(x.nzval)
nonzeros(x::SparseVectorView) = x.nzval
nonzeroinds(x::SparseVectorView) = x.nzind


text:
module SurveyDocumentation
include("SamplingProcedures.jl")
Include("SurveyAnalysis.jl")
include("SurveyDesign.jl")

end


text:
### MittagLefflerE

## FIXME: organize this more like polylog is organized
## in general, we need a way to combine effective pattern matching with efficient dispatch.
## i.e. we want fewer `isa(...)`s.
## Maybe one of the Julia matching packages will do this for us.
## PatternDispatch.jl ?
## Match.jl appears to not write functions taking advantage of multiple dispatch
## Stefan mentioned something about ameliorating ambiguity hell in v0.6.
##
## for alpha == 0 , must have abs(z) < 1
## Mma does not return conditions, I think.
## We could optionally return conditions.

@sjdoc MittagLefflerE """
    MittagLeffler(α,β,z)

    MittagLeffler(α,z)

represents the Mittag-Leffler function.
"""

@mkapprule MittagLefflerE :nargs => 2:3

@doap MittagLefflerE(α,z) = mittagleffler(mx,α,z)

@doap function MittagLefflerE(α,β,z)
    β == 1  && return mittagleffler(mx,α,z)
    mittagleffler(mx,α,β,z)
end

function mittagleffler(mx,α,z)
    α == 1//2 && return Exp(mpow(z,2)) * Erfc(mminus(z))
    α == -1//2 && return Exp(mpow(z,2)) * Erfc(z)
    α == 0 && return mpow(mminus(1,z),-1)
    α == 1 && return Exp(z)
    α == 2 && return Cosh(Sqrt(z))
    α == 3 && return (1//3)* (Exp(mpow(z,(1//3))) + 2*Exp(-mpow(z,1//3)/2) * Cos(Sqrt(3)/2 * mpow(z,1//3)))
    α == 4 && return (1//2)* ( Cosh(mpow(z,1//4)) + Cos(mpow(z,1//4)))
    if (isa(α,Number) && isa(z, Number)) &&
        (isa(α,AbstractFloat) || isa(z, FloatRC))
        return MittLeff.mittleff(float(α),float(z))
    end
    z == 0 && return mpow(Gamma(β),-1)
    mx
end

function mittagleffler(mx,α,β,z)
    if (isa(α,AbstractFloat) || isa(β,AbstractFloat) || isa(z, FloatRC)) &&
        (isa(α,Number) && isa(β,Number) && isa(z, Number))
       return MittLeff.mittleff(float(α),float(β),float(z))
    end
    z == 0 && return mpow(Gamma(β),-1)
    if α == 1
        β == 2 && return (mmul(mplus(Exp(z),-1),mpow(z,-1)))
        return (mpow(:E,z)*mpow(z,mminus(1,β)) * mxpr(:GammaRegularized, mminus(β,1),0,z))
    end
    α == 2 && β == 2 && return mmul(Sinh(mpow(z,1//2)), mpow(mpow(z,1//2),-1))
    mx
end


text:
"""
    struct Polygon{N,T} <: AbstractPlanarSurfacePrimitive{T}

Surface primitive describing a planar polygon, e.g. the base of a [`RegularPrism`](@ref).

## Parametric types
* `T`: Precision type.
* `N`: Number of vertices of the `Polygon`.

## Fields
* `points::SVector{N, CartesianPoint{T}}`: Vertices of the polygon in the order with which they are connected.
"""
struct Polygon{N,T} <: AbstractPlanarSurfacePrimitive{T}
    points::SVector{N, CartesianPoint{T}}
end
const Triangle{T} = Polygon{3, T}
const Quadrangle{T} = Polygon{4, T}

normal(p::Polygon) = normalize((p.points[2] - p.points[1]) × (p.points[3] - p.points[1]))

vertices(p::Polygon) = p.points
vertices(p::Polygon, n::Int64) = vertices(p)
connections(p::Polygon{N}) where {N} = [collect(1:N)]
connections(p::Polygon, ::Int64) = connections(p)
connections(p::Polygon{N}, ::Int64, ::Int64) where {N} = [[i%N+1, (i+1)%N+1] for i in 0:N-1]
get_label_name(::Polygon) = "Polygon"

extreme_points(p::Polygon) = p.points

function edges(p::Triangle{T}) where {T}
    vs = vertices(p)
    return (
        Edge(vs[1], vs[2]),
        Edge(vs[2], vs[3]),
        Edge(vs[3], vs[1])
    )
end
function edges(p::Quadrangle{T}) where {T}
    vs = vertices(p)
    return (
        Edge(vs[1], vs[2]),
        Edge(vs[2], vs[3]),
        Edge(vs[3], vs[4]),
        Edge(vs[4], vs[1])
    )
end

flip(p::Polygon) = Polygon(reverse(p.points))

lines(p::Polygon) = edges(p)

Plane(p::Polygon{N, T}) where {N, T} = Plane{T}(p.points[1], (p.points[2] - p.points[1]) × (p.points[3] - p.points[1]))

function _get_rot_for_rotation_on_xy_plane(p::Polygon{<:Any, T}) where {T}
    n = normal(p)
    rot_angle_rad = CartesianVector{T}(0,0,1) ⋅ n
    return if abs(rot_angle_rad) == 1
        RotationVec(zero(T), zero(T), zero(T))
    else
        rot_angle = -acos(rot_angle_rad)
        rot_axis = rot_angle * normalize(CartesianVector{T}(0,0,1) × n)
        RotationVec(rot_axis[1], rot_axis[2], rot_axis[3])
    end
end

function _rotate_on_xy_plane(p::Polygon{<:Any, T}) where {T}
    rot = _get_rot_for_rotation_on_xy_plane(p)
    map(pt -> rot * pt, p.points)
end

function in(pt::CartesianPoint{T}, p::Polygon{N,T}) where {N,T}
    plane = Plane(p)
    if (pt - origin(plane)) ⋅ normal(plane) == 0
        rot = _get_rot_for_rotation_on_xy_plane(p)
        vs = vertices(p)
        pts2d = SVector{N+1, SVector{2, T}}(
            view(rot * vs[i%N + 1], 1:2) for i in 0:N
        )
        # PolygonOps.inpolygon -> in = 1, on = -1, out = 0)
        b = PolygonOps.inpolygon(view(rot * pt, 1:2), pts2d) != 0
    end
    return b
end

function _above_or_below_polygon(pt::AbstractCoordinatePoint, p::Quadrangle{T}) where {T}
    rot = _get_rot_for_rotation_on_xy_plane(p)
    vs = vertices(p)
    pts2d = SVector{5, SVector{2, T}}(
        view(rot * vs[1], 1:2),
        view(rot * vs[2], 1:2),
        view(rot * vs[3], 1:2),
        view(rot * vs[4], 1:2),
        view(rot * vs[1], 1:2),
    )
    # PolygonOps.inpolygon -> in = 1, on = -1, out = 0)
    return PolygonOps.inpolygon(view(rot * pt, 1:2), pts2d) != 0
end

function _transform_into_global_coordinate_system(poly::Polygon{N, T}, p::AbstractPrimitive{T}) where {N, T}
    Polygon{N,T}(broadcast(pt -> _transform_into_global_coordinate_system(pt, p), poly.points))
end

function distance(pt::CartesianPoint, p::Polygon)
    return if _above_or_below_polygon(pt, p)
        distance(pt, Plane(p))
    else
        es = edges(p)
        ds = map(e -> distance(pt, e), es)
        minimum(ds)
    end
end


text:
using Test
@testset "ResourcePools" begin
    include("pooled_resource.jl")
    include("resource_pool.jl")
    include("pooled_abstract_array.jl")
    include("pooled_dense_array.jl")
end


text:
include("../src/lap_algs/helpers.jl")

k, l = 1, 1

A = rand(123, 123, 123)

var1 = A[l, :, :]
time1 = @benchmark A[l, :, :] = rand(123,123)
time1 = mean(time1.times)

B = rand(12345, 4, 4)

var2 = A[:, l, :]
time2 = @benchmark A[:, l, :] = rand(123,123)
time2 = mean(time2.times)

var3 = A[:, :, l]
time3 = @benchmark A[:, :, l] = rand(123,123)
time3 = mean(time3.times)

##
A = rand(2, 2, 65536);
@benchmark A[k, l, :] .= rand(65536)

A = rand(65536, 2, 2);
@benchmark A[:, k, l] .= rand(65536)

##

N = 5
C = rand(N, N, 12345)
d = rand(N, 12345)

time4 = @benchmark multi_mat_div_gem(C, d)
time4 = mean(time4.times)

time5 = @benchmark multi_mat_div(C, d)
time5 = mean(time5.times)

time5/time4

out4 = multi_mat_div_gem(C, d)
out5 = multi_mat_div(C, d)

maximum(abs.(out4-out5)) < 10e-7

## test single lap with and without gem

img, imgw, flow = gen_init();
fhs = 23
window = [47, 47]

showflow(flow)

# no gem
oldflow = LAP_julia.lap(img, imgw, fhs, window)
time6 = @benchmark LAP_julia.lap(img, imgw, fhs, window)
time6 = mean(time6.times)

# gem
newflow = LAP_julia.lap(img, imgw, fhs, window)
time7 = @benchmark LAP_julia.lap(img, imgw, fhs, window)
time7 = mean(time7.times)

time6/time7

## multi mat div 2    with views and without

@benchmark (@views multi_mat_div_gem($C, $d))
@benchmark multi_mat_div_gem($C, $d)

#--> views are better here

## try cumsum of matrix
img = rand(2565, 2565)

function cumsum2d(A)
    result = cumsum(A, dims = 1)
    cumsum!(result, result, dims = 2)
    return result
end

function cumsum2d!(result, A)
    cumsum!(result, A, dims = 1)
    cumsum!(result, result, dims = 2)
    return result
end

window_sum2
S(x,y) = a(x,y) + S(x-1,y) + S(x,y-1) - S(x-1,y-1)

include("../src/lap_algs/helpers.jl")

imgg = copy(img)
window_sum!(imgg, img, (256, 256), (47,47))

imgg

##

small = [1 2 3; 4 5 6; 7 8 9]

@inline function get_sum_around_ind(ind, B, radius) # square window
    hw = radius
    x, y = ind[1], ind[2]
    return B[x + hw, y + hw] - B[x + hw, y - hw - 1] - B[x - hw - 1, y + hw] + B[x - hw - 1, y - hw - 1]
end

function window_sum2!(result, pixels, img_size, window, timer)
    @timeit timer "borders" begin
        w = Int64((window-1)/2)
        summed = padarray(pixels, Pad(:symmetric, (w, w),
                                           (w, w)))
    end
    @timeit timer "cumsum" begin
        summed = cumsum2d!(summed, summed)
        B = padarray(summed, Fill(0, (1, 1), (1, 1)))
    end

    @timeit timer "for loop" begin
        @simd for ind in CartesianIndices(pixels)
            result[ind] = get_sum_around_ind(ind, B, w)
        end
    end
    return result
end

function window_sum3!(result, pixels, img_size, window, timer)
    @timeit timer "borders" begin
        w = Int64((window-1)/2)
        summed = padarray(pixels, Pad(:symmetric, (w, w),
                                           (w, w)))
    end
    @timeit timer "cumsum" begin
        summed = cumsum2d!(summed, summed)
        B = padarray(summed, Fill(0, (1, 1), (0, 0)))
    end
    c = size(pixels, 1)

    @timeit timer "calculations" begin
        result .= view(B, (1+w):(w+c), (1+w):(w+c)) .-
                  view(B, (1+w):(w+c), (-w):(c-w-1)) .-
                  view(B, (-w):(c-w-1), (1+w):(w+c)) .+
                  view(B, (-w):(c-w-1), (-w):(c-w-1))
    end
end

A = rand(234, 234)
result3 = Float64.(copy(A))
result2 = Float64.(copy(A))
result1 = Float64.(copy(A))

# timed
timer = TimerOutput();

reset_timer!(timer)
@timeit timer "whole" window_sum3!(result3, A, size(A), 21, timer)
print_timer(timer)

reset_timer!(timer)
@timeit timer "whole" window_sum2!(result2, A, size(A), 21, timer)
print_timer(timer)

result3
result2

A = small
# bench
@benchmark window_sum3!(result3, A, size(A), 9, timer)
@benchmark window_sum3!(result3, A, size(A), 9)
@benchmark window_sum2!(result2, A, size(A), 31, timer)
@benchmark window_sum!(result1, A, size(A), (9,9))

maximum(abs.(result1 - result2))

## time classic alg improved
TimerOutputs.enable_debug_timings(LAP_julia.registration)

img, imgw, flow = gen_init()

fhs = 23
window = [47, 47]

fhs = 3
window = [7, 7]
newflow = LAP_julia.lap(img, imgw, fhs, window)


text:
using JuMP.MathOptInterface: AbstractOptimizer, OptimizerWithAttributes
using JuMP: Model, @variable, @objective, @constraint, optimize!,
            termination_status, objective_value, value, primal_status

# solver status
function is_lp_optimal(status)
    return status == JuMP.MathOptInterface.OPTIMAL
end

function is_lp_infeasible(status; strict::Bool=false)
    if status == JuMP.MathOptInterface.INFEASIBLE
        return true
    end
    if strict
        return false
    end
    return status == JuMP.MathOptInterface.INFEASIBLE_OR_UNBOUNDED
end

function is_lp_unbounded(status; strict::Bool=false)
    if status == JuMP.MathOptInterface.DUAL_INFEASIBLE
        return true
    end
    if strict
        return false
    end
    return status == JuMP.MathOptInterface.INFEASIBLE_OR_UNBOUNDED
end

function has_lp_infeasibility_ray(model)
    return primal_status(model) == JuMP.MathOptInterface.INFEASIBILITY_CERTIFICATE
end

# solve a linear program (in the old MathProgBase interface)
function linprog(c, A, sense::Char, b, l::Number, u::Number, solver)
    n = length(c)
    m = length(b)
    return linprog(c, A, fill(sense, m), b, fill(l, n), fill(u, n), solver)
end

function linprog(c, A, sense, b, l, u, solver)
    n = length(c)
    model = Model(solver)
    @variable(model, l[i] <= x[i=1:n] <= u[i])
    @objective(model, Min, c' * x)
    eq_rows, ge_rows, le_rows = sense .== '=', sense .== '>', sense .== '<'
    @constraint(model, A[eq_rows, :] * x .== b[eq_rows])
    @constraint(model, A[ge_rows, :] * x .>= b[ge_rows])
    @constraint(model, A[le_rows, :] * x .<= b[le_rows])
    optimize!(model)
    return (
        status = termination_status(model),
        objval = objective_value(model),
        sol = value.(x),
        model = model
    )
end


text:
using Documenter, Amincoacid_to_RNA_Translation

makedocs(
    modules = [Amincoacid_to_RNA_Translation],
    format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Kathatinka",
    sitename = "Amincoacid_to_RNA_Translation.jl",
    pages = Any["index.md"]
    # strict = true,
    # clean = true,
    # checkdocs = :exports,
)

deploydocs(
    repo = "github.com/Kathatinka/Amincoacid_to_RNA_Translation.jl.git",
    push_preview = true
)


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

# Tests for deprecated functionality.
#
# These can't be run with --depwarn=error, so currently require special
# treatment when run inside the test system.

using Test
using Logging

using Base: remove_linenums!

module DeprecationTests # to test @deprecate
    f() = true

    # test the Symbol path of @deprecate
    @deprecate f1 f
    @deprecate f2 f false # test that f2 is not exported

    # test the Expr path of @deprecate
    @deprecate f3() f()
    @deprecate f4() f() false # test that f4 is not exported
    @deprecate f5(x::T) where T f()

    # test deprecation of a constructor
    struct A{T} end
    @deprecate A{T}(x::S) where {T, S} f()

    # test that @deprecate_moved can be overridden by an import
    Base.@deprecate_moved foo1234 "Foo"
    Base.@deprecate_moved bar "Bar" false
end # module
module Foo1234
    export foo1234
    foo1234(x) = x+1
end

# issue #21972
struct T21972
    @noinline function T21972()
        Base.depwarn("something", :T21972)
        new()
    end
end

@testset "@deprecate" begin
    using .DeprecationTests
    using .Foo1234
    @test foo1234(3) == 4
    @test_throws ErrorException DeprecationTests.bar(3)

    # 22845
    ex = :(module M22845; import ..DeprecationTests: bar;
                          bar(x::Number) = x + 3; end)
    @test_warn "importing deprecated binding" eval(ex)
    @test @test_nowarn(DeprecationTests.bar(4)) == 7

    # enable when issue #22043 is fixed
    # @test @test_warn "f1 is deprecated, use f instead." f1()
    # @test @test_nowarn f1()

    # @test_throws UndefVarError f2() # not exported
    # @test @test_warn "f2 is deprecated, use f instead." DeprecationTests.f2()
    # @test @test_nowarn DeprecationTests.f2()

    # @test @test_warn "f3() is deprecated, use f() instead." f3()
    # @test @test_nowarn f3()

    # @test_throws UndefVarError f4() # not exported
    # @test @test_warn "f4() is deprecated, use f() instead." DeprecationTests.f4()
    # @test @test_nowarn DeprecationTests.f4()

    # @test @test_warn "f5(x::T) where T is deprecated, use f() instead." f5(1)
    # @test @test_nowarn f5(1)

    # @test @test_warn "A{T}(x::S) where {T, S} is deprecated, use f() instead." A{Int}(1.)
    # @test @test_nowarn A{Int}(1.)

    # issue #21972
    @noinline function f21972()
        T21972()
    end
    @test_deprecated "something" f21972()
end

f24658() = depwarn24658()

depwarn24658() = Base.firstcaller(backtrace(), :_func_not_found_)

@testset "firstcaller" begin
    # issue #24658
    @test eval(:(if true; f24658(); end)) == (Ptr{Cvoid}(0),StackTraces.UNKNOWN)
end

# issue #25130
f25130() = Base.depwarn("f25130 message", :f25130)
# The following test is for the depwarn behavior of expressions evaluated at
# top-level, so we can't use the usual `collect_test_logs()` / `with_logger()`
testlogger = Test.TestLogger()
prev_logger = global_logger(testlogger)
# Each call at top level should be distinct. This won't be true if they're
# attributed to internal C frames (including generic dispatch machinery)
f25130()
f25130()
testlogs = testlogger.logs
@test length(testlogs) == 2
@test testlogs[1].id != testlogs[2].id
@test testlogs[1].kwargs[:caller].func == Symbol("top-level scope")
@test all(l.message == "f25130 message" for l in testlogs)
global_logger(prev_logger)

#-------------------------------------------------------------------------------
# BEGIN 0.7 deprecations

@testset "parser syntax deprecations" begin
    # #15524
    # @test (@test_deprecated Meta.parse("for a=b f() end")) == :(for a=b; f() end)
    @test_broken length(Test.collect_test_logs(()->Meta.parse("for a=b f() end"))[1]) > 0
end

# END 0.7 deprecations

@testset "tuple indexed by float deprecation" begin
    @test_deprecated getindex((1,), 1.0) === 1
    @test_deprecated getindex((1,2), 2.0) === 2
    @test_throws Exception getindex((), 1.0)
    @test_throws Exception getindex((1,2), 0.0)
    @test_throws Exception getindex((1,2), -1.0)
end

max_stars_repo_name: xgdgsc/julia

max_issues_repo_name: xgdgsc/julia

max_forks_repo_name: xgdgsc/julia


text:
import Test: @test, @test_throws, @testset
import Logging: @info

testdir = dirname(@__FILE__)

import Luna: set_fftw_mode
set_fftw_mode(:estimate)

@testset "All" begin

@testset "Maths" begin
    @info("================= test_maths.jl")
    include(joinpath(testdir, "test_maths.jl"))
end

@testset "PhysData" begin
    @info("================= test_physdata.jl")
    include(joinpath(testdir, "test_physdata.jl"))
end

@testset "Capillary" begin
    @info("================= test_capillary.jl")
    include(joinpath(testdir, "test_capillary.jl"))
end

@testset "Rectangular Modes" begin
    @info("================= test_rect_modes.jl")
    include(joinpath(testdir, "test_rect_modes.jl"))
end

@testset "ODE Solver" begin
    @info("================= test_rk45.jl")
    include(joinpath(testdir, "test_rk45.jl"))
end

@testset "Ionisation" begin
    @info("================= test_ionisation.jl")
    include(joinpath(testdir, "test_ionisation.jl"))
end

@testset "Output" begin
    @info("================= test_output.jl")
    include(joinpath(testdir, "test_output.jl"))
end

@testset "Multimode" begin
    @info("================= test_multimode.jl")
    include(joinpath(testdir, "test_multimode.jl"))
end

@testset "Polarisation" begin
    @info("================= test_polarisation.jl")
    include(joinpath(testdir, "test_polarisation.jl"))
    @info("================= test_polarisation_field.jl")
    include(joinpath(testdir, "test_polarisation_field.jl"))
    @info("================= test_polarisation_env.jl")
    include(joinpath(testdir, "test_polarisation_env.jl"))
end

@testset "Tools" begin
    @info("================= test_tools.jl")
    include(joinpath(testdir, "test_tools.jl"))
end

@testset "Utils" begin
    @info("================= test_utils.jl")
    include(joinpath(testdir, "test_utils.jl"))
end

@testset "Gradients" begin
    @info("================= test_gradient.jl")
    include(joinpath(testdir, "test_gradient.jl"))
end

@testset "Tapers" begin
    @info("================= test_tapers.jl")
    include(joinpath(testdir, "test_tapers.jl"))
end

@testset "Scans" begin
    @info("================= test_scans.jl")
    include(joinpath(testdir, "test_scans.jl"))
end

@testset "Raman" begin
    @info("================= test_raman.jl")
    include(joinpath(testdir, "test_raman.jl"))
end

@testset "Kerr" begin
    @info("================= test_kerr.jl")
    include(joinpath(testdir, "test_kerr.jl"))
end

@testset "LinearOps" begin
    @info("================= test_linops.jl")
    include(joinpath(testdir, "test_linops.jl"))
end

@testset "Modes" begin
    @info("================= test_modes.jl")
    include(joinpath(testdir, "test_modes.jl"))
end

@testset "Radial Propagation" begin
    @info("================= test_radial.jl")
    include(joinpath(testdir, "test_radial.jl"))
end

@testset "Full 3D Propagation" begin
    @info("================= test_full_freespace.jl")
    include(joinpath(testdir, "test_full_freespace.jl"))
end

@testset "Antiresonant modes" begin
    @info("================= test_antiresonant.jl")
    include(joinpath(testdir, "test_antiresonant.jl"))
end

@testset "Fields" begin
    @info("================= test_fields.jl")
    include(joinpath(testdir, "test_fields.jl"))
end

@testset "Processing" begin
    @info("================= test_processing.jl")
    include(joinpath(testdir, "test_processing.jl"))
end

@testset "Vector plasma" begin
    @info("================= test_vectorplasma.jl")
    include(joinpath(testdir, "test_vectorplasma.jl"))
end

@testset "Statistics" begin
    @info("================= test_stats.jl")
    include(joinpath(testdir, "test_stats.jl"))
end

@testset "Gas mixtures" begin
    @info("================= test_mixtures.jl")
    include(joinpath(testdir, "test_mixtures.jl"))
end

@testset "Interface" begin
    @info("================= test_interface.jl")
    include(joinpath(testdir, "test_interface.jl"))
end

@testset "Linear propagation" begin
    @info("================= test_linearprop.jl")
    include(joinpath(testdir, "test_linearprop.jl"))
end

end


text:
"""
Deletes existing decision tree.
"""
function clear_tree!(p::PISPlanner)
    p.tree = nothing
end

"""
Constructs a PISTree and choose an action.
"""
POMDPs.action(p::PISPlanner, s) = first(action_info(p, s))

"""
Constructs a PISTree and choose the best action.
"""
function POMDPModelTools.action_info(p::PISPlanner, s; tree_in_info=false, β=0.0, γ=1.0, schedule=0.1)
    local a::Any  # actiontype(p.mdp)
    info = Dict{Symbol, Any}()
    try
        if isterminal(p.mdp, s)
            error("MCTS cannot handle terminal states: s = $s")
        end

        tree = p.tree
        if !p.solver.keep_tree || tree == nothing
            # TODO(kykim): Temporary hack to get around the RMDP case.
            # Should ideally be actiontype(p.mdp).
            tree = PISTree{statetype(p.mdp),Any}()
            p.tree = tree
        end
        snode = insert_state_node!(tree, s)

        timer = p.solver.timer
        start_s = timer()
        timeout_s = start_s + p.solver.max_time
        n_iterations = p.solver.n_iterations
        p.solver.show_progress ? progress = Progress(n_iterations) : nothing

        sim_channel = Channel{Task}(min(10_000, n_iterations)) do channel
            # TODO(kykim): Try pulling out the add_sample! and update! calls out
            # from the simulate() method.
            for n in 1:n_iterations
                put!(channel, Threads.@spawn simulate(p, snode, p.solver.depth, 0.0, 0.0,
                                                      β, γ, schedule, timeout_s))
            end
        end

        nquery = 0
        for sim_task in sim_channel
            if timer() > timeout_s
                p.solver.show_progress ? finish!(progress) : nothing
                @show Timeout reached
                break
            end

            try
                fetch(sim_task)  # Throws a TaskFailedException if failed.
                nquery += 1
                p.solver.show_progress ? next!(progress) : nothing
            catch err
                throw(err.task.exception)  # Throw the underlying exception.
            end
        end

        p.reset_callback(p.mdp, s)  # Optional: Leave the MDP in the current state.
        info[:search_time_s] = (timer() - start_s)
        info[:tree_queries] = nquery
        if p.solver.tree_in_info || tree_in_info
            info[:tree] = tree
        end

        sanode = best_sanode(tree, snode)
        a = sanode.a_label
    catch ex
        a = next_action(p.next_action, p.mdp, s, snode)
        info[:exception] = ex
    end

    return a, info
end

"""
Adds a sample to the tree.
"""
function add_sample!(tree::PISTree, cost::Float64, weight::Float64)
    Base.@lock tree.costs_weights_lock begin
        push!(tree.costs, cost)
        push!(tree.weights, weight)
    end
    return cost, weight
end

"""
Returns the reward for one iteration of MCTS.
"""
function simulate(dpw::PISPlanner, snode::PISStateNode, d::Int,
                  cost::Float64=0.0, weight::Float64=0.0,
                  β::Float64=0.0, γ::Float64=1.0, schedule::Float64=0.1,
                  timeout_s::Float64=0.0)
    sol = dpw.solver
    tree = dpw.tree
    s = snode.s_label
    action_distrib_fn = sol.nominal_distrib_fn
    action_distrib = action_distrib_fn(dpw.mdp, s)

    Base.@lock tree.costs_weights_lock begin
        n_samples = length(tree.costs)
        for i in tree.cdf_est.last_i+1:n_samples
            # Depending on how fast/slow this is, making copies of costs and weights first then
            # immediately releasing the lock can be better. But, of course, as we collect more
            # samples, copying can increasingly be slow.
            ImportanceWeightedRiskMetrics.update!(tree.cdf_est, tree.costs[i], exp(tree.weights[i]))
        end
    end

    if isterminal(dpw.mdp, s)
        return add_sample!(tree, cost, weight)
    elseif d == 0
        out_cost, out_weight = estimate_value(dpw.mdp, s, d, cost, weight,
                                              sol.cost_reduction,
                                              action_distrib_fn)
        return add_sample!(tree, out_cost, out_weight)
    end

    # Action progressive widening.
    if sol.enable_action_pw
        if n_children(snode) <= sol.k_action * total_n(snode)^sol.alpha_action
            a = next_action(dpw.next_action, dpw.mdp, s, snode)
            insert_action_node!(tree, snode, a,
                                init_N(sol.init_N, dpw.mdp, s, a),
                                init_Q(sol.init_Q, dpw.mdp, s, a))
        end
    elseif n_children(snode) == 0
        for a in support(actions(dpw.mdp, s))
            insert_action_node!(tree, snode, a,
                                init_N(sol.init_N, dpw.mdp, s, a),
                                init_Q(sol.init_Q, dpw.mdp, s, a))
        end
    end

    sanode, q_logprob = sample_sanode_UCB(dpw, snode, action_distrib, β, γ, schedule)
    a = sanode.a_label
    w_node = compute_weight(q_logprob, a, action_distrib)

    # State progressive widening.
    spnode = nothing
    new_node = false
    if ((sol.enable_state_pw && n_a_children(sanode) <= sol.k_state * n(sanode)^sol.alpha_state) ||
        n_a_children(sanode) == 0)
        sp, r = @gen(:sp, :r)(dpw.mdp, s, a, dpw.rng)
        Base.@lock tree.state_nodes_lock begin
            if haskey(tree.state_nodes, sp)
                spnode = tree.state_nodes[sp]
            else
                spnode = insert_state_node!(tree, sp)
                new_node = true
            end
        end
        Base.@lock sanode.a_lock begin
            push!(sanode.transitions, (sp, r))
            if !(sp in sanode.unique_transitions)
                push!(sanode.unique_transitions, sp)
                sanode.n_a_children += 1
            end
        end
    else
        sp, r = rand(dpw.rng, sanode.transitions)
        spnode = Base.@lock tree.state_nodes_lock begin; tree.state_nodes[sp]; end
    end

    new_weight = weight + w_node
    new_cost = update_cost(cost, r, sol.cost_reduction)
    if new_node
        out_cost, out_weight = estimate_value(dpw.mdp, sp, d - 1, new_cost,
                                              new_weight, sol.cost_reduction,
                                              action_distrib_fn)
        add_sample!(tree, out_cost, out_weight)
        q = discount(dpw.mdp) * out_cost
    else
        out_cost, out_weight = simulate(
            dpw, spnode, d - 1, new_cost, new_weight, β, γ, schedule, timeout_s)
        q = discount(dpw.mdp) * out_cost
    end

    # Backpropagate and update node values.
    Base.@lock sanode.a_lock begin
        Base.@lock snode.s_lock begin
            snode.total_n += 1
            delete!(snode.a_selected, sanode.a_label)
            sanode.n += 1
            sanode.q += (q - sanode.q) / sanode.n
            ImportanceWeightedRiskMetrics.update!(
                sanode.c_cdf_est, q, exp(out_weight))
        end
    end

    return out_cost, out_weight
end

"""
Computes importance weights i.e., p/q.
"""
function compute_weight(q_logprob, a, distribution)
    if distribution == nothing
        w = -q_logprob
    else
        w = logpdf(distribution, a) - q_logprob
    end
    return w
end

"""
Returns a sampled action and the corresponding log probability.
"""
function select_action(sanodes, action_probs)
    sanode_idx = sample(1:length(sanodes), Weights(action_probs))
    sanode = sanodes[sanode_idx]
    q_logprob = log(action_probs[sanode_idx])
    return sanode, q_logprob
end

"""
Computes the probabilities with which to sample from the actions.
"""
function compute_action_probs(a_selection::Symbol, ucb_scores,
                              est_α_probs, est_α_costs, nominal_probs,
                              β, γ)
    a_selection == :ucb && return ucb_probs(ucb_scores)
    a_selection == :ucb_softmax && return ucb_softmax_probs(ucb_scores)
    a_selection == :expected_cost && return expected_cost_probs()
    a_selection == :adaptive &&
        return adaptive_probs(est_α_probs, est_α_costs, nominal_probs, β, γ)
    # TODO(kykim): Boltzmann exploration type strategy.
    throw("Unsupported action selection strategy: $(a_selection)")
end

"""
Returns an adaptive α to use for an action.
"""
function action_α(a_n, schedule, α)
    w_annealed = 1.0 / (1.0 + schedule * a_n)
    return w_annealed + (1 - w_annealed) * α
end

function sample_sanode_UCB(dpw::PISPlanner, snode::PISStateNode, action_distrib,
                           β::Float64, γ::Float64, schedule::Float64)
    tree = dpw.tree
    sol = dpw.solver
    c = sol.exploration_constant
    virtual_loss = sol.virtual_loss
    α = sol.α

    ucb_scores = []
    sanodes = []
    ltn = log(total_n(snode))
    est_α_costs = []
    est_α_probs = []
    nominal_probs = []

    # Instead of creating a copy, we can presumably just compute the no. of children and read only
    # that many elements of children(snode). But, this is also technically not safe as the vector
    # can grow in size and be resized around the same time when we access the element.
    sa_children = Base.@lock snode.s_lock begin; deepcopy(children(snode)) end
    for a_label in sa_children
        Base.@lock tree.state_action_nodes_lock begin
            state_action_key = sanode_key(snode.s_label, a_label)
            if !haskey(tree.state_action_nodes, state_action_key)
                continue
            end
            sanode = tree.state_action_nodes[state_action_key]
        end
        a_n = n(sanode)
        a_q = q(sanode)
        if (ltn <= 0 && a_n == 0) || c == 0.0
            UCB = a_q
        else
            UCB = a_q + c * sqrt(ltn / a_n)
        end

        # If applicable, apply virtual loss to the score.
        vloss = Base.@lock snode.s_lock begin; (a_label in snode.a_selected ? virtual_loss : 0.0) end;
        UCB -= vloss

        @assert !isnan(UCB) "UCB was NaN (q=$a_q, c=$c, ltn=$ltn, n=$a_n)"
        @assert !isequal(UCB, -Inf) "UCB was -Inf (q=$a_q, c=$c, ltn=$ltn, n=$a_n)"

        push!(ucb_scores, UCB)
        push!(sanodes, sanode)

        a_α = action_α(a_n, schedule, α)
        est_quantile = ImportanceWeightedRiskMetrics.quantile(tree.cdf_est, a_α)
        c_tail = ImportanceWeightedRiskMetrics.tail_cost(sanode.c_cdf_est, est_quantile)
        c_cdf = ImportanceWeightedRiskMetrics.cdf(sanode.c_cdf_est, est_quantile)
        push!(est_α_costs, c_tail)
        push!(est_α_probs, 1.0 - c_cdf)
        push!(nominal_probs, pdf(action_distrib, a_label))
    end

    action_probs = compute_action_probs(sol.action_selection, ucb_scores,
                                        est_α_probs, est_α_costs, nominal_probs,
                                        β, γ)
    sanode, q_logprob = select_action(sanodes, action_probs)
    Base.@lock snode.s_lock begin; push!(snode.a_selected, sanode.a_label); end
    return sanode, q_logprob
end

"""
Returns the best action.
Some publications say to choose action that has been visited the most
e.g., Continuous Upper Confidence Trees by Couëtoux et al.
"""
function best_sanode(tree::PISTree, snode::PISStateNode)
    best_Q = -Inf
    best_sa = nothing
    for a_label in children(snode)
        state_action_key = sanode_key(snode.s_label, a_label)
        sanode = tree.state_action_nodes[state_action_key]
        if q(sanode) > best_Q
            best_Q = q(sanode)
            best_sa = sanode
        end
    end
    return best_sa
end


text:
function OceanTidesReader(arg0::JString)
    return OceanTidesReader((JString,), arg0)
end

function load_data(obj::DataLoader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end


text:
include("../src/PowerFlowNetworks.jl")
using .PowerFlowNetworks
using SQLite
using Dates
using DataFrames
using ArgParse

ArgParse.parse_item(::Type{Vector{String}}, x::AbstractString) = split(x, ",")

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--dbpath", "-p"
            help = "Path of the sqlite database file."
            default = "data/PowerFlowNetworks.sqlite"
        "--indirs_rawgo"
            help = "Directory/ies containing the RAWGO networks."
            arg_type = Vector{String}
            default = readdir("data/RAWGO"; join=true)
        "--indirs_matpowerm"
            help = "Directory/ies containing the MATPOWERM networks."
            arg_type = Vector{String}
            default = ["data/MATPOWERM"]
    end
    return parse_args(s)
end

function clean_dirs(dirs; verbose=true)
    cleaned_dirs = [d for d in dirs if isdir(d)]
    if length(cleaned_dirs) < length(dirs)
        @warn "$(setdiff(dirs, cleaned_dirs)) are not valid directory/ies (skipped)."
    end
    return cleaned_dirs
end

function load_rawgo_instance!(db::SQLite.DB, path::AbstractString)
    date = Dates.now()
    name = splitdir(path)[end]
    scenarios = readdir(path)
    @info "Loading $(name)"
    for scenario in scenarios
        !isdir(joinpath(path, scenario)) && continue
        scenario_nb = parse(Int, scenario[10:end])
        source_path = abspath(joinpath(path, scenario, "case.raw"))
        @info "  $(scenario_nb)"
        try
            load_instance_in_db!(db, name, scenario_nb, source_path, "RAWGO", date)
        catch e
            if isa(e, SQLiteException) && e.msg == "UNIQUE constraint failed: instances.name, instances.scenario"
                @warn "$(name) scenario $(scenario_nb) already in the database. Loading ignored."
            else
                rethrow()
            end
        end
    end
end

function load_rawgo_instances!(db::SQLite.DB, dirs)
    for dir in dirs
        for instance_dir in readdir(dir; join=true)
            !isdir(instance_dir) && continue
            load_rawgo_instance!(db, instance_dir)
        end
    end
end

function load_matpower_mat_instance!(db::SQLite.DB, path::AbstractString)
    date = Dates.now()
    name = splitext(basename(path))[end - 1]
    scenario_nb = 0
    source_path = abspath(path)
    source_type = "MATPOWERM"
    @info "Loading $(name)"
    try
        load_instance_in_db!(db, name, scenario_nb, source_path, source_type, date)
    catch e
        if isa(e, SQLiteException) && e.msg == "UNIQUE constraint failed: instances.name, instances.scenario"
            @warn "$(name) scenario $(scenario_nb) already in the database. Loading ignored."
        else
            rethrow()
        end
    end
end

function load_matpower_mat_instances!(db::SQLite.DB, dirs)
    for dir in dirs
        files = readdir(dir; join=true)
        for f in files
            !startswith(basename(f), "case") && continue
            load_matpower_mat_instance!(db, f)
        end
    end
end

function main()
    parsed_args = parse_commandline()
    indirs_rawgo = clean_dirs(parsed_args["indirs_rawgo"])
    indirs_matpowerm = clean_dirs(parsed_args["indirs_matpowerm"])
    dbpath = parsed_args["dbpath"]
    isdb = isfile(dbpath)
    db = SQLite.DB(dbpath)
    tables = ["instances", "decompositions", "mergers", "combinations", "solve_results"]
    if !isdb || any(map(x -> !(x in SQLite.tables(db)[:name]), tables))
        db = setup_db(dbpath)
    end
    load_matpower_mat_instances!(db, indirs_matpowerm)
    load_rawgo_instances!(db, indirs_rawgo)
end

main()

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using POMDPs # This defines the interface for mdps
using POMDPPolicies # Used just for the "FunctionPolicy" type for rendering
using POMDPGym # This is where the environments are defined (including rendering code)
using Crux # This the deep RL library with the various solvers
using Flux # This is the ML library for defining models

# Define the lavaworld mdp
mdp = LavaWorldMDP(lava = [GWPos(2,5), GWPos(3,5), GWPos(4,5), GWPos(5,5)], goal=GWPos(7,5))
S = state_space(mdp) # this is the state space of the mdp
as = actions(mdp) # Actions of the mdp
idim = prod(Crux.dim(S)) # input dimension

# See the lavaworld rendered
render(mdp)

# Define the Q network and solve with DQN
Q() = DiscreteNetwork(Chain(flatten, Dense(idim, 64, relu), Dense(64, 64, relu), Dense(64, length(as))), as)
solver = DQN(π=Q(), S=S, N=10000)
π_exp = solve(solver, mdp)

# Record some expert trajectories
Neps = 10
B_expert = ExperienceBuffer(episodes!(Sampler(mdp, π_exp), Neps=Neps))
mean_return = sum(B_expert[:r]) / Neps

## Lets see how our expert does
# Plot the policy
render(mdp, policy=π_exp, return_compose=true)

# Plot the value function (The factor of 10 is a quirk of an external render function that assumes values are between -10 and 10)
render(mdp, color = s->10.0*maximum(value(π_exp, convert_s(Vector{Float32}, s, mdp))))

# Plot the occupancy of different grid points in the expert data
occupancy = POMDPGym.gen_occupancy(B_expert, mdp)
render(mdp, color = s->200f0*occupancy[s])

## Now lets solve using imitation learning

# Start with a simple BC baseline
𝒮_bc = BC(π=Q(), S=S, 𝒟_demo=B_expert, opt=(;epochs=10000), log=(;period=100,))
π_bc = solve(𝒮_bc, mdp)
render(mdp, policy=π_bc, return_compose=true)

# then use off-policy GAIL
D() = ContinuousNetwork(Chain(flatten, DenseSN(idim + 4, 64, relu), DenseSN(64, 64, relu), DenseSN(64, 2)))
𝒮_gail = OffPolicyGAIL(D=D(),
                       𝒟_demo=B_expert,
                       solver=DQN,
                       π=Q(),
                       S=S,
                       N=10000,
                       )
π_gail = solve(𝒮_gail, mdp)

render(mdp, policy=π_gail, return_compose=true)

# Now lets see how it does with a different lava configuration:
mdp2 = LavaWorldMDP(lava = [GWPos(2,3), GWPos(3,3), GWPos(4,3), GWPos(5,3)], goal=GWPos(7,5))
render(mdp2, policy=π_gail, return_compose=true)


text:
#module Streaming

import Base: *, +

using Interpolations
#using FastGaussQuadrature
using LinearAlgebra
using SpecialFunctions
using ForwardDiff
using DiffRules

import ForwardDiff:Dual,value,partials,derivative,extract_derivative

"""
    Params(ϵ,Re)

Set the parameters for the streaming solution. The problem is scaled so
that the radius of the cylinder is unity. Reynolds number is defined as Re = ΩR²/ν,
where Ω is the angular frequency of the oscillatory motion, e.g. sin(Ωt).
"""
struct Params
    ϵ :: Float64
    Re :: Float64
    γ² :: ComplexF64
    γ :: ComplexF64
    λ :: ComplexF64
    λ² :: ComplexF64
    H₀ :: ComplexF64
    C :: ComplexF64
  end

function Params(ϵ::Number,Re::Number)
        γ² = im*Re
        γ = exp(im*π/4)*√Re
        λ = √2*γ
        λ² = 2*γ²
        H₀ = hankelh1(0,γ)
        C = hankelh1(2,γ)/H₀
        Params(ϵ,Re,γ²,γ,λ,λ²,H₀,C)
    end

function Base.show(io::IO, p::Params) where {N}
        println(io, "Streaming flow parameters with Re = $(p.Re), ϵ = $(p.ϵ)")
end

"""
    ComplexFunc(f)

Provides a wrapper for a function expected to return complex values, for use in
dispatch in automatic differentiation with `ForwardDiff`.
"""
struct ComplexFunc{FT}
    fcn::FT
end

(f::ComplexFunc)(x) = f.fcn(x)

macro create_dual(f,nu,c,d,hankel)
    SF = :SpecialFunctions
    H = :($hankel)
    _f = Symbol("_",f)
    defs = quote
                $(_f)(x)  = $SF.$H($nu,$c*x)/$d
                $f = ComplexFunc($(_f))
           end
    return esc(defs)
end

macro extend_H(hankel)
    SF = :SpecialFunctions
    H = :($hankel)
    _, dH = DiffRules.diffrule(SF,H,:(ν),:(v))
    Hdef = quote
              function $SF.$H(ν,d::Complex{<:Dual{T}}) where {T}
                    dr, di = reim(d)
                    v = value(dr)+im*value(di)
                    return Dual{T}(real($SF.$H(ν,v)),real($dH)*partials(dr)-imag($dH)*partials(di)) +
                        im*Dual{T}(imag($SF.$H(ν,v)),imag($dH)*partials(dr)+real($dH)*partials(di))
              end
            end
    return esc(Hdef)
end

# Derivatives #

@inline function ForwardDiff.derivative(f::ComplexFunc{F}, x::R) where {F,R<:Real}
    T = typeof(ForwardDiff.Tag(f.fcn, R))
    return ForwardDiff.extract_derivative(T, real(f.fcn(Dual{T}(x, one(x))))) +
        im*ForwardDiff.extract_derivative(T, imag(f.fcn(Dual{T}(x, one(x)))))
end

"""
    D²(f::ComplexFunc,K::Int) -> ComplexFunc

Return the radial part of the Laplacian at `r` of a function having cylindrical form

\$ f(r)e^{iKθ} \$
"""
@inline function D²(f::ComplexFunc,K::Int)
      rdf = ComplexFunc(r -> r*derivative(f,r))
      drdf = ComplexFunc(r -> derivative(rdf,r))
      return ComplexFunc(r -> drdf(r)/r - K^2*f(r)/r^2)
end

"""
    curl(f::ComplexFunc,K::Int) -> Tuple{ComplexFunc,ComplexFunc}

Return the radial part of the cylindrical components of the curl at `r` of a
function having cylindrical form

\$ f(r)e^{iKθ} \$
"""
@inline function curl(f::ComplexFunc,K::Int)
  return ComplexFunc(r -> K*f(r)/r),ComplexFunc(r -> -derivative(f,r))
end

# Integrals #
"""
    Integral(f,a<:Real,b<:Real[,fake_infinity=1e5],[,length=10000])

Set up an integral

\$ \\int_a^r f(\\tau)d\\tau \$

for a < r < b
If `b` is set to `Inf`, then this sets up the integral

\$ \\int_r^\\infty f(\\tau)d\\tau \$

where a < r < ∞
The optional argument `length` sets the size of the interpolation table, and
`fake_infinity` sets the maximum value at which to truncate the limit at
infinity.
"""
struct Integral{FT,V,IL}
    fcn::FT
    a::Float64
    b::Float64
    fake_infinity::Float64
    table::Interpolations.AbstractInterpolation
end

function Integral(fcn,a,b;fake_infinity=1e6,length=10000)
    V = typeof(fcn(1))
    v = zeros(V,length)
    τ = range(0,stop=1,length=length)
    if b == Inf
        # Integration variable τ = c/x - c/fake_infinity, where c is
        # such that τ goes from 1 to 0 when x goes from a to fake_infinity.
        lim_diff = 1/a-1/fake_infinity
        c = 1/lim_diff
        xi = c/(0+c/fake_infinity)
        fcn_old = fcn(xi)*xi^2/c
        for (i,τi) in enumerate(τ)
          if i > 1
            xi = c/(τi+c/fake_infinity)
            fcn_new = fcn(xi)*xi^2/c
            v[i] += v[i-1] + 0.5*step(τ)*(fcn_new+fcn_old)
            fcn_old = fcn_new
          end
        end
        table = LinearInterpolation(τ, v, extrapolation_bc = Flat())
        Integral{typeof(fcn),V,true}(fcn,a,b,fake_infinity,table)
    else
        # Integration variable τ = c(x - a), where c is such
        # that τ goes from 0 to 1 when x goes from a to b
        lim_diff = b-a
        c = 1/lim_diff
        fcn_old = fcn(a)
        for (i,τi) in enumerate(τ)
          if i > 1
            xi = a + τi*lim_diff
            fcn_new = fcn(xi)
            v[i] += v[i-1] + 0.5*step(τ)*lim_diff*(fcn_new+fcn_old)
            fcn_old = fcn_new
          end
        end
        table = LinearInterpolation(τ, v, extrapolation_bc = Flat())
        Integral{typeof(fcn),V,false}(fcn,a,b,fake_infinity,table)
    end
end

function Base.show(io::IO, I::Integral{FT,V,IL}) where {FT,V,IL}
    if IL
      println(io, "Integral from $(I.a) < r to $(I.b)")
    else
      println(io, "Integral from $(I.a) to r < $(I.b)")
    end
end

ComplexIntegral(fcn::ComplexFunc,a,b;ops...) = ComplexFunc(Integral(fcn,a,b;ops...))

ComplexIntegral(fcn::Function,a,b;ops...) = ComplexFunc(Integral(ComplexFunc(fcn),a,b;ops...))

function (I::Integral{FT,V,false})(r) where {FT,V<:Number}
  return r > I.a ? I.table((r-I.a)/(I.b-I.a)) : zero(V)
end

function (I::Integral{FT,V,true})(r) where {FT,V<:Number}
  return I.table((1/r-1/I.fake_infinity)/(1/I.a-1/I.fake_infinity))
end

#(I::Integral{FT,V,true} where {FT,V<:Number})(d::Dual{T}) where {T} = Dual{T}(I(value(d)),-I.fcn(value(d))*partials(d))

#(I::Integral{FT,V,false} where {FT,V<:Number})(d::Dual{T}) where {T} = Dual{T}(I(value(d)),I.fcn(value(d))*partials(d))

# Now define how they behave with Dual types. The following defines the automatic differentiation:
function (I::Integral{FT,V,true} where {FT<:ComplexFunc,V<:Number})(d::Dual{T}) where {T}
          return Dual{T}(real(I(value(d))),-real(I.fcn(value(d)))*partials(d)) +
       im*Dual{T}(imag(I(value(d))),-imag(I.fcn(value(d)))*partials(d))
end

(I::Integral{FT,V,false} where {FT<:ComplexFunc,V<:Number})(d::Dual{T}) where {T} =
         Dual{T}(real(I(value(d))),real(I.fcn(value(d)))*partials(d)) +
      im*Dual{T}(imag(I(value(d))),imag(I.fcn(value(d)))*partials(d))

# Set up functions #

@extend_H(hankelh1)
@extend_H(hankelh2)

@create_dual(X,0,p.γ,p.H₀,hankelh1)
@create_dual(Y,1,p.γ,p.H₀,hankelh1)
@create_dual(Z,2,p.γ,p.H₀,hankelh1)

@create_dual(H11,1,p.λ,1,hankelh1)
@create_dual(H12,1,p.λ,1,hankelh2)
@create_dual(H21,2,p.λ,1,hankelh1)
@create_dual(H22,2,p.λ,1,hankelh2)

abstract type Order end
abstract type First <: Order end
abstract type Second <: Order end

struct ComplexAmplitude{FT,K}
    f :: FT
end

ComplexAmplitude(f,K) = ComplexAmplitude{typeof(f),K == 1 ? First : Second}(f)

function (A::ComplexAmplitude{FT,Second})(x,y) where {FT}
    r2 = x^2+y^2
    r = sqrt(r2)
    return 2*real(A.f(r))*x*y/r2
end

#### Construct the first and second order solutions

#end


text:
import RowEchelon

#=
Homework 2
Section 1.4
Problem 3
=#

# Create the first Matrix
A = [1 2 ;
     -3 1 ;
     1 6]

# Create the second matrix
B = [-2 ; 3]

# Multiply the A and B
matrixProduct = A*B

# Print the results
println()
println("Original Matrix, A:")
println(A)
println("Original Matrix, B")
println(B)
println("Product Matrix, A*B:")
println(matrixProduct)
println()


text:
function main()
    println(devnull, 1)
    println(devnull, 2) # COV_EXCL_LINE
    println(devnull, 3)
    # COV_EXCL_START
    println(devnull, 4)
    # COV_EXCL_STOP
    println(devnull, 5)
end


text:
#=##############################################################################
# DESCRIPTION
    Tools for aerodynamic calculations.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

"""
    Calculates the aerodynamic force at each element of a VLM system using it's
current `Gamma` solution, the velocity induced by the particle field `pfield`,
and the kinematic velocity calculated from `prev_vlm_system`. It saves the force
as the field `vlm_system.sol["Ftot"]`
"""
function calc_aerodynamicforce(vlm_system::Union{vlm.Wing, vlm.WingSystem},
                                prev_vlm_system, pfield, Vinf, dt, rho; t=0.0,
                                per_unit_span=false,
                                Vout=nothing, lenout=nothing,
                                lencrit=-1)

    m = vlm.get_m(vlm_system)    # Number of horseshoes

    # Nodes of every horseshoe
    Ap = _get_Xs(vlm_system, "Ap")
    A = _get_Xs(vlm_system, "A")
    B = _get_Xs(vlm_system, "B")
    Bp = _get_Xs(vlm_system, "Bp")

    # Midpoints of bound vortices
    ApA = (Ap .+ A)/2
    AB = (A .+ B)/2
    BBp = (B .+ Bp)/2

    # Evaluate VPM on each midpoint
    Xs = vcat(ApA, AB, BBp)

    ## NOTE: Instead of calling the VPM, we use what was calculated
    ## by `solve()`, which includes Rotor-on-VLM velocities
    # Vvpm = Vvpm_on_Xs(pfield, Xs; dt=dt)
    Vvpm = vcat(vlm_system.sol["Vvpm_ApA"], vlm_system.sol["Vvpm_AB"], vlm_system.sol["Vvpm_BBp"])

    # Evaluate VLM on each midpoint
    Vvlm = vlm.Vind.(Ref(vlm_system), Xs; t=t, ign_col=true, ign_infvortex=true)

    # Evaluate Vinf on each midpoint
    Vinfs = Vinf.(Xs, t)

    # Evaluate kinematic velocity on each node
    Vtran = _Vkinematic(vlm_system, prev_vlm_system, dt; t=t,
                                                targetX=["Ap", "A", "B", "Bp"])

    Ftot = [zeros(3) for i in 1:m]

    for i in 1:m                 # Iterate over horseshoes
        for j in 1:3             # Iterate over bound vortices (1==ApA, 2==AB, 3==BBp)

            # Bound vortex' length vector
            if j==1
                l = (A[i]-Ap[i])
            elseif j==2
                l = (B[i]-A[i])
            else
                l = (Bp[i]-B[i])
            end

            # Bound vortex' midpoint
            # X = Xs[i + m*(j-1)]

            # Kinematic velocity at midpoint
            Vtranmid = (Vtran[i + m*(j-1)] + Vtran[i + m*j])/2

            # Effective velocity at midpoint
            V = Vvpm[i + m*(j-1)] + Vvlm[i + m*(j-1)] + Vinfs[i + m*(j-1)] + Vtranmid
            if Vout != nothing
                push!(Vout, [Vvpm[i + m*(j-1)], Vvlm[i + m*(j-1)], Vinfs[i + m*(j-1)], Vtranmid])
            end

            # Circulation
            Gamma = vlm_system.sol["Gamma"][i]

            # Normalization factor
            if per_unit_span
                # len = norm(cross(V/norm(V),l))
                unitV = V/sqrt(V[1]*V[1]+V[2]*V[2]+V[3]*V[3])
                projl = zeros(3)
                projl[1] = unitV[2]*l[3] - unitV[3]*l[2]
                projl[2] = unitV[3]*l[1] - unitV[1]*l[3]
                projl[3] = unitV[1]*l[2] - unitV[2]*l[1]
                len = sqrt(projl[1]*projl[1]+projl[2]*projl[2]+projl[3]*projl[3])
            else
                len = 1
            end
            if lenout != nothing
                push!(lenout, per_unit_span==false || len>lencrit ? len : -10*lencrit)
            end

            # Kutta–Joukowski theorem: F = rho * V × vecGamma
            if per_unit_span==false || len > lencrit # NOTE: Use lencrit to avoid dividing by zero
                Ftot[i][1] += rho * Gamma * (V[2]*l[3] - V[3]*l[2]) / len
                Ftot[i][2] += rho * Gamma * (V[3]*l[1] - V[1]*l[3]) / len
                Ftot[i][3] += rho * Gamma * (V[1]*l[2] - V[2]*l[1]) / len
            end

        end
    end

    vlm._addsolution(vlm_system, per_unit_span ? "ftot" : "Ftot", Ftot; t=t)

    return Ftot
end

"""
Receives an array of vectors `Fs` and decomposes them into components in the
`ihat`, `jhat`, `khat` unitary orthogonal directions (`khat` is deduced).
"""
function decompose(Fs, ihat, jhat)

    if abs(ihat[1]*ihat[1]+ihat[2]*ihat[2]+ihat[3]*ihat[3] - 1) >= 1e-12
        error("Got non-unitary vector ihat=$ihat")
    elseif abs(jhat[1]*jhat[1]+jhat[2]*jhat[2]+jhat[3]*jhat[3] - 1) >= 1e-12
        error("Got non-unitary vector jhat=$jhat")
    end

    khat = zeros(3)
    khat[1] = ihat[2]*jhat[3] - ihat[3]*jhat[2]
    khat[2] = ihat[3]*jhat[1] - ihat[1]*jhat[3]
    khat[3] = ihat[1]*jhat[2] - ihat[2]*jhat[1]

    Fis = [zeros(3) for i in 1:length(Fs)]
    Fjs = [zeros(3) for i in 1:length(Fs)]
    Fks = [zeros(3) for i in 1:length(Fs)]

    for (n, F) in enumerate(Fs)
        Fis[n] .= (F[1]*ihat[1]+F[2]*ihat[2]+F[3]*ihat[3]) * ihat
        Fjs[n] .= (F[1]*jhat[1]+F[2]*jhat[2]+F[3]*jhat[3]) * jhat
        Fks[n] .= (F[1]*khat[1]+F[2]*khat[2]+F[3]*khat[3]) * khat
    end

    return Fis, Fjs, Fks
end

"""
    `remove_particles_lowstrength(crit_Gamma2, every_nsteps)`

Returns an extra_runtime_function that every `step` steps removes all
particles that have a squared-magnitude Gamma smaller than `crit_Gamma2`.
"""
function remove_particles_lowstrength(crit_Gamma2::Real, step::Int)

    function wake_treatment(sim, PFIELD, T, DT, args...; optargs...)
        if sim.nt%step==0

            for i in vpm.get_np(PFIELD):-1:1
                P = vpm.get_particle(PFIELD, i)

                if P.Gamma[1]*P.Gamma[1] + P.Gamma[2]*P.Gamma[2] + P.Gamma[3]*P.Gamma[3] < crit_Gamma2
                    vpm.remove_particle(PFIELD, i)
                end
            end

        end

        return false
    end

    return wake_treatment
end

"""
    `remove_particles_lowstrength(crit_Gamma2, every_nsteps)`

Returns an extra_runtime_function that every `step` steps removes all
particles that are outside of a box of minimum and maximum vertices `Pmin`
and `Pmax`.
"""
function remove_particles_box(Pmin, Pmax, step::Int)

    function wake_treatment(sim, PFIELD, T, DT, args...; optargs...)
        if sim.nt%step==0

            for i in vpm.get_np(PFIELD):-1:1
                P = vpm.get_particle(PFIELD, i)

                if (  (P.X[1] < Pmin[1] || P.X[1] > Pmax[1])
                        ||
                      (P.X[2] < Pmin[2] || P.X[2] > Pmax[2])
                        ||
                      (P.X[3] < Pmin[3] || P.X[3] > Pmax[3])
                   )
                    vpm.remove_particle(PFIELD, i)
                end
            end

        end

        return false
    end

    return wake_treatment
end


text:
function material_scattering_coefficients(scat_field::ScatteringCoefficientsField;
        numdensity = (x1, s1) -> number_density(s1),
        rtol::AbstractFloat = 1e-2,
        maxevals::Int = Int(5e4)
    )

    lm2n = lm_to_spherical_harmonic_index
    rθφ2xyz = radial_to_cartesian_coordinates

    v = regular_basis_function(scat_field.medium, scat_field.ω)

    particle_radius = maximum(outer_radius.(scat_field.material.species))
    R = outer_radius(scat_field.material.shape)

    function kernel(rθφ)
        x = rθφ2xyz(rθφ)

        vs = conj.(v(scat_field.basis_order + scat_field.basis_field_order, x))
        fs = scat_field.coefficient_field(x)

        (rθφ[1]^2 * sin(rθφ[2]) * numdensity(x,scat_field.material.species[1])) .* [
            sum(
                gaunt_coefficient(l,m,dl,dm,l1,m-dm) * vs[lm2n(l1,m-dm)] * fs[lm2n(dl,dm)]
            for dl = 0:scat_field.basis_order for dm = -dl:dl for l1 in max(abs(m-dm),abs(dl-l)):(dl+l))
        for l = 0:scat_field.basis_field_order for m = -l:l]
    end

    (v,err) = hcubature(kernel, SVector(0.0,0.0,-π), SVector(R - particle_radius,π,π)
        ; rtol=rtol, maxevals = maxevals
    )

    println("Integration error of $err")

    return v
end

function material_scattering_coefficients(scat_field::ScatteringCoefficientsField{Sphere{T,3},Acoustic{T,3},RadialSymmetry{3}};
        numdensity = (x1, s1) -> number_density(s1),
        rtol::AbstractFloat = 1e-4,
        maxevals::Int = Int(1e5),
    ) where T

    k = scat_field.ω / scat_field.medium.c
    R = outer_radius(scat_field.material.shape)

    lm2n = lm_to_spherical_harmonic_index
    ls = 0:scat_field.basis_order

    rθφ2xyz = radial_to_cartesian_coordinates

    particle_radius = maximum(outer_radius.(scat_field.material.species))

    function kernel(r)
        x = rθφ2xyz([r,0.0,0.0])

        js = sbesselj.(ls, k*r)

        fs = scat_field.coefficient_field(x)

        (4.0 * π * r^2 * numdensity([r,0,0],scat_field.material.species[1])) .* sum(
            (-1.0) .^ ls .* sqrt.(2.0 .* ls .+ 1.0) .* fs[lm2n.(ls,0)] .* js
        )
    end

    (v,err) = hquadrature(kernel, 0.0, R - particle_radius; rtol=rtol, maxevals = maxevals)

    println("Integration error of $err")

    return [v]

end

max_stars_repo_name: arturgower/EffectiveWaves


text:
#####
# AbstractBody
#####

abstract type AbstractBody{T} end

get_shape(body::AbstractBody) = body.shape

get_position(body::AbstractBody) = body.position
set_position!(body::AbstractBody{T}, position::SA.SVector{2, T}) where {T} = body.position = position

get_position_change(body::AbstractBody) = body.position_change
set_position_change!(body::AbstractBody{T}, position_change::SA.SVector{2, T}) where {T} = body.position_change = position_change

add_position_change!(body::AbstractBody{T}, position_change::SA.SVector{2, T}) where {T} = set_position_change!(body, get_position_change(body) .+ position_change)

function apply_position_change!(body)
    position = get_position(body)
    position_change = get_position_change(body)
    set_position!(body, position .+ position_change)
    set_position_change!(body, zero(position_change))
    return body
end

get_velocity(body::AbstractBody) = body.velocity
set_velocity!(body::AbstractBody{T}, velocity::SA.SVector{2, T}) where {T} = body.velocity = velocity

get_velocity_change(body::AbstractBody) = body.velocity_change
set_velocity_change!(body::AbstractBody{T}, velocity_change::SA.SVector{2, T}) where {T} = body.velocity_change = velocity_change

add_velocity_change!(body::AbstractBody{T}, velocity_change::SA.SVector{2, T}) where {T} = set_velocity_change!(body, get_velocity_change(body) .+ velocity_change)

function apply_velocity_change!(body)
    velocity = get_velocity(body)
    velocity_change = get_velocity_change(body)
    set_velocity!(body, velocity .+ velocity_change)
    set_velocity_change!(body, zero(velocity_change))
    return body
end

get_inv_mass(body::AbstractBody) = body.inv_mass


text:
function poly_derivative!(
    df::T1,
    f::T1,
    pdm::T2,
) where {T1<:AbstractVector,T2<:AbstractMatrix}
    @assert length(f) == size(pdm, 1)
    for i in eachindex(df)
        df[i] = dot(f, pdm[i, :])
    end
end


text:
@testset "research.jl" begin

    data_surv = CSO.dataset("cso_raw")
    data_trudeau = CSO.dataset("Montreal-Trudeau_15min_precipitations")
    data = innerjoin(data_surv, data_trudeau, on=:Date)
    dropmissing!(data, :Duration)

    topredict = "Duration"

    train, test = train_test_year(data, 2020)

    result = CR_linreg(topredict, [:d15min, :d45min, :d1h, :d6h, :d24h], train, test)
    varnames = result[1, :Variables]
    RMSD = result[1, :RMSD]

    #The RMSD is well derived
    @test RMSD == linreg("Duration", varnames, train, test)

    #The df is well sorted
    @test RMSD == minimum(result.RMSD)

    #The best set of variables found stepwise corresponds to the best set of the CR_linreg
    @test countmap(varnames) == countmap(stepRMSD_linreg(topredict, [:d45min, :d15min, :d1h, :d6h, :d24h], [], train, test, false))

end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO