text:
# # DIRS
# const HEERDEN_RAW_DATA_DIR =
#     joinpath(RAW_DATA_DIR, "Heerden2013")
# const HEERDEN_PROCESSED_DATA_DIR =
#     joinpath(PROCESSED_DATA_DIR, basename(HEERDEN_RAW_DATA_DIR))
# const HEERDEN_FIGURES_DIR =
#     joinpath(FIGURES_DATA_DIR, basename(HEERDEN_RAW_DATA_DIR))

# function _create_dirs()
#     for dir in [HEERDEN_PROCESSED_DATA_DIR, HEERDEN_FIGURES_DIR]
#         try; mkpath(dir); catch err; @warn("Error creating dir", dir, err); end
#     end
# end

# # FILES
# const HEERDEN_CONT_CUL_DATA_ORIG_FILE =
#     joinpath(HEERDEN_RAW_DATA_DIR, "heerden2013___cont_cult_data.tsv")
# const HEERDEN_CONT_CUL_DATA_CONV_FILE =
#     joinpath(HEERDEN_PROCESSED_DATA_DIR, basename(HEERDEN_CONT_CUL_DATA_ORIG_FILE))
# const HEERDEN_GENE_MODS_FILE =
#     joinpath(HEERDEN_RAW_DATA_DIR, "heerden2013___data", "heerden2013___gene_modifications.tsv")


text:
"""
    MOKernel

An abstract type for multi-output kernels.
"""
abstract type MOKernel <: Kernel end

ext: jl


text:
# Setting up Bayesian inference, plotting distributions and
# Bayesian linear regression
using Plots
using Distributions

μ = 5.0
σ = 4.0
n = 1000
lim = 10
x = range(0, stop = lim, length=1000)

f(μ, σ, x) = μ .+ (σ * rand(Normal(0,1),lim))

y = f(μ, σ, x)

function likelihood!(y, n, sd)
    yy = repeat(y, outer=(1,n))
    ran = range(0, stop = length(y), length = n)
    xx = repeat(ran, outer=(1, length(y)))'
    likelihood = exp.(-0.5 * sum((yy .- xx ).^2,dims=1) / σ^2)
    norm_likelihood = likelihood / sum(likelihood)
    return norm_likelihood
end

likelihood = likelihood!(y, n, σ)

ρ = exp.((-0.5 .* (x .- 0 ).^2) / 5^2)
ρ = ρ/sum(ρ)

plot(x, likelihood[1,:])
plot!(x, ρ)
scatter!(y, zeros(n))

#= Linear regression
Use linear regression to estimate the posterior mean and std
y = B*x + e
e = N(0,s^2)

y = N(B^T * x, s^2 * I)
=#

beta1 = 1 / σ^2
beta2 = 1 / 5^2

post_σ = 1 / √(lim * beta1 + beta2)
post_μ = (lim * beta1 * mean(y) + beta2 * 0) / (lim * beta1 + beta2)

posterior = pdf(Normal(post_μ, post_σ), x)
posterior = posterior/sum(posterior)

plot!(x, posterior)


text:
const FACTORY = () -> begin
    model = Cbc.Optimizer()
    MOI.set(model, MOI.Silent(), true)
    return model
end
const OPTIMIZER = FACTORY()
const MODEL = MOI.Utilities.CachingOptimizer(
    MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    OPTIMIZER,
)

# Feasible region: {x in R^3 :
#     x[1] + x[2] + x[3] <= 1
#     2x[1] - x[3] >= 0
#     0 <= x <= 1
#     x[3] in {0,1}
function _repopulate_model!(model)
    MOI.empty!(model)
    _x = [
        MOI.add_constrained_variable(model, MOI.Interval(0.0, 1.0)) for
        i in 1:3
    ]
    x = [MOI.SingleVariable(_x[i][1]) for i in 1:3]
    MOI.add_constraint(
        model,
        1.0 * x[1] + 1.0 * x[2] + 1.0 * x[3],
        MOI.LessThan(1.0),
    )
    MOI.add_constraint(model, 2.0 * x[1] - 1.0 * x[3], MOI.GreaterThan(0.0))
    MOI.add_constraint(model, x[3], MOI.ZeroOne())
    return x
end

@testset "maximum_activity" begin
    @testset "IntervalArithmetic" begin
        x = _repopulate_model!(MODEL)

        @test maximum_activity(
            DisjunctiveConstraints.IntervalArithmetic(),
            MODEL,
            1.0 * x[2] + 3.0 * x[3],
        ) ≈ 4.0
        @test minimum_activity(
            DisjunctiveConstraints.IntervalArithmetic(),
            MODEL,
            1.0 * x[2] + 3.0 * x[3],
        ) ≈ 0.0
    end

    @testset "LinearProgrammingRelaxation" begin
        x = _repopulate_model!(MODEL)

        @test maximum_activity(
            DisjunctiveConstraints.LinearProgrammingRelaxation(FACTORY),
            MODEL,
            1.0 * x[2] + 3.0 * x[3],
        ) ≈ 2.0
        @test minimum_activity(
            DisjunctiveConstraints.LinearProgrammingRelaxation(FACTORY),
            MODEL,
            1.0 * x[2] + 3.0 * x[3],
        ) ≈ 0.0
    end

    @testset "FullFormulation" begin
        x = _repopulate_model!(MODEL)

        @test maximum_activity(
            DisjunctiveConstraints.FullFormulation(FACTORY),
            MODEL,
            1.0 * x[2] + 3.0 * x[3],
        ) ≈ 1.0
        @test minimum_activity(
            DisjunctiveConstraints.FullFormulation(FACTORY),
            MODEL,
            1.0 * x[2] + 3.0 * x[3],
        ) ≈ 0.0
    end
end


text:
# EXAMPLE: scatter
# here we have 3 ranks

using MPI

MPI.Init()
comm = MPI.COMM_WORLD
rank = MPI.Comm_rank(comm)

if rank == 0
    data = [7.0,8.0,9.0,10.0,11.0,12.0]
else
    data = Float64[]
end

data = MPI.Scatter(data, 2, 0, comm) #each rank has 2 elements
@show rank, data

MPI.Finalize()

## output:
# (rank, data) = (1, [9.0, 10.0])
# (rank, data) = (2, [11.0, 12.0])
# (rank, data) = (0, [7.0, 8.0])


text:
@testset "TransformerAnalysis" begin

trace_edge_gap = 0.16e-3
trace_trace_gap = 0.13e-3
outer = copper_weight_to_meters(2.0)
inner = copper_weight_to_meters(1.0)
dielectric = (1.6e-3-2*outer-2*inner)/3
stackup = Stackup([copper,fr408,copper,fr408,copper,fr408,copper],
                  [outer,dielectric,inner,dielectric,inner,dielectric,outer])
pcb = PCB_Specification(trace_edge_gap,
                        trace_trace_gap,
                        stackup)
core_geometry = core_geometry_dict["e14_set"]
primaryturnsperlayer=4
secondaryturnsperlayer=1
isprimaryeries = true
issecondaryseries = false
w1 = windings(pcb, core_geometry, primaryturnsperlayer, secondaryturnsperlayer,
             prisec"S-P-P-S", isprimaryeries, issecondaryseries)
w2 = windings(pcb, core_geometry, primaryturnsperlayer, secondaryturnsperlayer,
            prisec"S-P-S-P", isprimaryeries, issecondaryseries)
w3 = windings(pcb, core_geometry, primaryturnsperlayer, secondaryturnsperlayer,
            prisec"S-P-S-P", false, true)
ferrite = ferrite_dict["4f1"]
t1 = Transformer(ferrite, w1)
t2 = Transformer(ferrite, w2)
t3 = Transformer(ferrite, w3)
@test PlanerTransformer.layer_resistance_tuple(windings(t1)) ==
  (0.002730986989493496, 0.07483432580022921, 0.07483432580022921, 0.002730986989493496)
@test PlanerTransformer.layer_resistance_tuple(windings(t2)) ==
  (0.002730986989493496, 0.07483432580022921, 0.005461973978986992, 0.037417162900114605)
@test PlanerTransformer.layer_resistance_tuple(windings(t3)) ==
  (0.002730986989493496, 0.07483432580022921, 0.005461973978986992, 0.037417162900114605)

frequency = 10e6
v_in = 15.0
i_out = -7.0
c1 = PushPull(v_in, i_out, frequency)
@test_throws ArgumentError PushPull(-v_in, i_out, frequency)
@test_throws ArgumentError PushPull(v_in, -i_out, frequency)
@test_throws ArgumentError PushPull(v_in, -i_out, -frequency)
duty = .25
c2 = Forward(v_in, i_out, frequency, duty)
@test_throws ArgumentError Forward(-v_in, i_out, frequency, duty)
@test_throws ArgumentError Forward(v_in, -i_out, frequency, duty)
@test_throws ArgumentError Forward(v_in, -i_out, -frequency, duty)

a_t1c1 = transformer_power_analysis(t1,c1)
@test capacitance(a_t1c1) ≈ 1.7949033519440568e-11
@test PlanerTransformer.pcb_thickness(a_t1c1) ≈ 0.0016000000000000003
@test flux_density(a_t1c1) ≈ 0.0031908618912421674
@test voltage(a_t1c1) == (15.0, 1.8223502434073622)
@test current(a_t1c1) == (0.8758650011589507, -7.0)
@test ac_winding_resistance(a_t1c1) == (0.2219529543698069, 0.004049950501870691)
@test dc_winding_resistance(a_t1c1) == (0.14966865160045842, 0.001365493494746748)
@test equilivent_resistance(a_t1c1) == (0.2219529543698069, 0.004049950501870691)
@test v_core(a_t1c1) ≈ 1.850699896920457
@test core_specific_power(a_t1c1) ≈ 42689.53481883259
@test core_total_power(a_t1c1) ≈ 0.012806860445649776
@test winding_power(a_t1c1) == (0.17026887849541192, 0.19844757459166384)
@test total_power(a_t1c1) ≈ 0.38152331353272556
@test r_core(a_t1c1) ≈ 267.44182331000746
@test input_power(a_t1c1) ≈ 13.13797501738426
@test output_power(a_t1c1) ≈ -12.756451703851535
@test efficiency(a_t1c1) ≈ 0.9709602649549957
@test leakage_inductance(a_t1c1) ≈ 1.4901268128799744e-7
@test abs(PlanerTransformer.power_error(a_t1c1))<1e-10

a_t1c2 = transformer_power_analysis(t1,c2)
@test capacitance(a_t1c2) ≈ 1.7949033519440568e-11
@test flux_density(a_t1c2) ≈ 0.0015971423866341192
@test voltage(a_t1c2) == (15.0, 1.8290333147449505)
@test current(a_t1c2) == (0.875598425879627, -7.0)
@test ac_winding_resistance(a_t1c2) == (0.2219529543698069, 0.004049950501870691)
@test dc_winding_resistance(a_t1c2) == (0.14966865160045842, 0.001365493494746748)
@test equilivent_resistance(a_t1c2) == (0.20388187867746976, 0.003378836250089705)
@test v_core(a_t1c2) ≈ 1.8526851684955783
@test core_specific_power(a_t1c2) ≈ 7391.29834419259
@test core_total_power(a_t1c2) ≈ 0.0022173895032577768
@test winding_power(a_t1c2) == (0.0390776626780815, 0.04139074406359889)
@test total_power(a_t1c2) ≈ 0.08268579624493817
@test r_core(a_t1c2) ≈ 386.99136174774026
@test input_power(a_t1c2) ≈ 3.2834940970486013
@test output_power(a_t1c2) ≈ -3.200808300803663
@test efficiency(a_t1c2) ≈ 0.9748177417711026
@test leakage_inductance(a_t1c2) ≈ 1.4901268128799744e-7
@test abs(PlanerTransformer.power_error(a_t1c2))<1e-10

a_t2c1 = transformer_power_analysis(t2,c1)
@test capacitance(a_t2c1) ≈ 2.6923550279160852e-11
@test flux_density(a_t2c1) ≈ 0.0031976841149601127
@test voltage(a_t2c1) == (15.0, 1.8318214383199423)
@test current(a_t2c1) == (0.8758678396008, -7.0)
@test ac_winding_resistance(a_t2c1) == (0.18581080298513264, 0.00326219262241758)
@test dc_winding_resistance(a_t2c1) == (0.11225148870034382, 0.001820657992995664)
@test equilivent_resistance(a_t2c1) == (0.18581080298513264, 0.00326219262241758)
@test v_core(a_t2c1) ≈ 1.8546567866768653
@test core_specific_power(a_t2c1) ≈ 42921.189476547006
@test core_total_power(a_t2c1) ≈ 0.012876356842964102
@test winding_power(a_t2c1) == (0.14254373043097793, 0.15984743849846142)
@test total_power(a_t2c1) ≈ 0.31526752577240347
@test r_core(a_t2c1) ≈  267.1370356007262
@test input_power(a_t2c1) ≈ 13.138017594011998
@test output_power(a_t2c1) ≈ -12.822750068239596
@test efficiency(a_t2c1) ≈ 0.9760034172951562
@test leakage_inductance(a_t2c1) ≈ 9.934178752533162e-8
@test abs(PlanerTransformer.power_error(a_t2c1))<1e-10

a_t2c2 = transformer_power_analysis(t2,c2)
@test capacitance(a_t2c2) ≈ 2.6923550279160852e-11
@test flux_density(a_t2c2) ≈ 0.0016005825568032011
@test voltage(a_t2c2) == (15.0, 1.8363631031362788)
@test current(a_t2c2) == (0.8756004041874144, -7.0)
@test ac_winding_resistance(a_t2c2) == (0.18581080298513264, 0.00326219262241758)
@test dc_winding_resistance(a_t2c2) == (0.11225148870034382, 0.001820657992995664)
@test equilivent_resistance(a_t2c2) == (0.16742097441393544, 0.002901808965062101)
@test v_core(a_t2c2) ≈ 1.8566757658917135
@test core_specific_power(a_t2c2) ≈ 7431.706030080913
@test core_total_power(a_t2c2) ≈ 0.0022295118090242738
@test winding_power(a_t2c2) == (0.032089413583281066, 0.03554715982201074)
@test total_power(a_t2c2) ≈ 0.06986608521431607
@test r_core(a_t2c2) ≈ 386.54705546931336
@test input_power(a_t2c2) ≈ 3.283501515702804
@test output_power(a_t2c2) ≈ -3.213635430488488
@test efficiency(a_t2c2) ≈ 0.9787220791949713
@test leakage_inductance(a_t2c2) ≈ 9.934178752533162e-8
@test abs(PlanerTransformer.power_error(a_t2c2))<1e-10

a_t3c1 = transformer_power_analysis(t3,c1)
@test capacitance(a_t3c1) ≈ 2.6923550279160852e-11
@test flux_density(a_t3c1) ≈ 0.006397992260147801
@test voltage(a_t3c1) == (15.0, 7.326737896892351)
@test current(a_t3c1) == (3.5050291347466387, -7.0)
@test ac_winding_resistance(a_t3c1) == (0.044695193801413506, 0.013561874982728374)
@test dc_winding_resistance(a_t3c1) == (0.024944775266743067, 0.00819296096848049)
@test equilivent_resistance(a_t3c1) == (0.044695193801413506, 0.013561874982728374)
@test v_core(a_t3c1) ≈ 3.7108355108857247
@test core_specific_power(a_t3c1) ≈ 248830.55742474739
@test core_total_power(a_t3c1) ≈ 0.0746491672274242
@test winding_power(a_t3c1) == (0.5490907015720118, 0.6645318741536903)
@test total_power(a_t3c1) ≈  1.2882717429531265
@test r_core(a_t3c1) ≈  184.46689628697774
@test input_power(a_t3c1) ≈ 52.57543702119958
@test output_power(a_t3c1) ≈ -51.287165278246455
@test efficiency(a_t3c1) ≈ 0.9754966992964097
@test leakage_inductance(a_t3c1) ≈ 2.4835446881332905e-8
@test abs(PlanerTransformer.power_error(a_t3c1))<1e-10

a_t3c2 = transformer_power_analysis(t3,c2)
@test capacitance(a_t3c2) ≈ 2.6923550279160852e-11
@test flux_density(a_t3c2) ≈ 0.0032027392435911095
@test voltage(a_t3c2) == (15.0, 7.34481751977721)
@test current(a_t3c2) == (3.503479779660588, -7.0)
@test ac_winding_resistance(a_t3c2) == (0.044695193801413506, 0.013561874982728374)
@test dc_winding_resistance(a_t3c2) == (0.024944775266743067, 0.00819296096848049)
@test equilivent_resistance(a_t3c2) == (0.039757589167745896, 0.012219646479166403)
@test v_core(a_t3c2) ≈ 3.7151775225656873
@test core_specific_power(a_t3c2) ≈ 43093.33059499487
@test core_total_power(a_t3c2) ≈ 0.01292799917849846
@test winding_power(a_t3c2) == (0.12199984556880146, 0.14969066936978845)
@test total_power(a_t3c2) ≈ 0.2846185141170884
@test r_core(a_t3c2) ≈ 266.9118367352115
@test input_power(a_t3c2) ≈ 13.138049173727206
@test output_power(a_t3c2) ≈ -12.853430659610117
@test efficiency(a_t3c2) ≈ 0.9783363184021069
@test leakage_inductance(a_t3c2) ≈ 2.4835446881332905e-8
@test abs(PlanerTransformer.power_error(a_t3c2))<1e-10
end # testset


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule PThreadPool_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("PThreadPool")
JLLWrappers.@generate_main_file("PThreadPool", UUID("a74eb8f4-56fd-5a14-80e9-60c06ef4dc0e"))
end  # module PThreadPool_jll


text:
using QuantumAnnealingTools, CUDA, DifferentialEquations
import QuantumAnnealingTools.⊗

function make_connect_tuples(n::Int64)
    """
    Make a list of all possible 2-body interactions for an
    [n] qubit system excluding duplicates, i.e. (i, j) where
    i < j and not both (i, j) and (j, i).

    Ex: [1, 2, 3] --> [(1, 2), (1, 3), (2, 3)]
    """
    connections = Array{Any,1}([])
    for prod in Iterators.product(1:n, 1:n)
        if (prod[1] < prod[2])
            push!(connections, prod)
        end
    end
    return connections
end

function make_spin_glass(n::Int64)
    """
    Create an [n] qubit spin-glass with couplings sampled
    (on the fly) from N(0, 1) distribution.
    """
    connections = make_connect_tuples(n)
    coeffs = randn(length(connections))
    return two_local_term(coeffs, connections, n)
end

function anneal_spin_glass(n::Int64, tf::T where T<:Number)
    """
    Performs a standard quantum anneal on an [n] qubit spin-glass for [tf]
    nano-seconds. Obtains spin-glass couplings from [load_dir].
    """
    # make Hamiltonians
    init_H = float_to_32_cpu(standard_driver(n))
    final_H = float_to_32_cpu(make_spin_glass(n))
    # get initial state
    init_state = ⊗([PauliVec[1][2] for _ in 1:n]...)
    H = DenseHamiltonian([(s)->1-s, (s)->s], [init_H, final_H], unit=:ħ)
    annealing = Annealing(H, init_state)
    sol = solve_schrodinger(annealing, tf)

    return sol
end

function anneal_spin_glass_gpu(n::Int64, tf::T where T<:Number)
    """
    Performs a standard quantum anneal on an [n] qubit spin-glass for [tf]
    nano-seconds. Obtains spin-glass couplings from [load_dir].
    """
    # make Hamiltonians
    init_H = float_to_32_gpu(cu(standard_driver(n)))
    final_H = float_to_32_gpu(cu(make_spin_glass(n)))
    # get initial state
    init_state = cu(⊗([PauliVec[1][2] for _ in 1:n]...))
    H = DenseHamiltonian([(s)->1-s, (s)->s], [init_H, final_H], unit=:ħ)
    annealing = Annealing(H, init_state)
    sol = solve_schrodinger_gpu(annealing, tf)

    return sol
end

function float_to_32_cpu(n::Array{Complex{Float64}})::Array{Complex{Float32}}
    """
    Convert to Float 64 arrays to Float 32 for improved runtime efficiency
    """
    return convert(Array{Complex{Float32}},n)
end

function float_to_32_gpu(n::CuArray{Complex{Float64}})::CuArray{Complex{Float32}}
    """
    Convert to Float 64 CUDA arrays to Float 32 for improved runtime efficiency
    """
    return convert(CuArray{Complex{Float32}},n)
end

@time anneal_spin_glass(3, 10);
@time anneal_spin_glass_gpu(3, 10);

max_stars_repo_name: naezzell/accelqat

max_issues_repo_name: naezzell/accelqat

max_forks_repo_name: naezzell/accelqat


text:
using Distributions, ProgressMeter, IterTools, LinearAlgebra
import StatsFuns.logistic

struct Data
  T::Array{Float64,1} # observation times
  δ::Array{Int64,1} # censoring indicator, 1 if exact observation
  Z::Array{Array{Float64,1},1} # covariate in dimension 1
end

struct Model
    data::Data
    obs::Array{ Float64, 1} # observed times, sorted and without repetitions
    last_inc::Float64 # to be added to the last observation to set overall time window in model
    incs::Array{ Float64, 1} # array of increments between observations
    E::Array{ Array{ Int64, 1}, 1} # array of indexes of exact observation corresponding to obs
    C::Array{ Array{ Int64, 1}, 1} # array of indexes of censored observation corresponding to obs
    f::Array{Function,1} # array of regression function in covariates corresponding to obs[i] CoRM dimension 1
    Hᵉ::Function # frequencies of exact observations in dimension 1
    Hᶜ::Function # frequencies of censored observations in dimension 1
    e_o::Array{ Int64, 1} # array of indexes of the exact observations in obs
    e_o_f::Function
end

# function to initialize model struct
@inline function Model(data::Data, last_inc::Float64, f::Array{Function,1})
      obs = sort( unique( data.T ) )
      incs = push!( diff( [ 0; obs ] ), last_inc)
      E = Array{Int64,1}[]
      C = Array{Int64,1}[]
      for i in 1:length( obs )
        push!( E, intersect( findall( data.T .== obs[i] ), findall( data.δ .== 1 ) ) )
        push!( C, intersect( findall( data.T .== obs[i] ), findall( data.δ .== 0 ) ) )
      end
      @inbounds hᵉ(k::Int64, β::Array{Float64,1}) = [ mapreduce( j -> f[k]( data.Z[j], β), +, E[i], init = 0.0 ) for i in 1:length(obs)]
      @inbounds hᶜ(k::Int64, β::Array{Float64,1}) = [ mapreduce( j -> f[k]( data.Z[j], β), +, C[i], init = 0.0 ) for i in 1:length(obs)]
      @inbounds Hᵉ(k::Int64, β::Array{Float64,1}) = [ cumsum( hᵉ(k,β)[end:-1:1] )[end:-1:1]; 0]
      @inbounds Hᶜ(k::Int64, β::Array{Float64,1}) = [ cumsum( hᶜ(k,β)[end:-1:1] )[end:-1:1]; 0]
      @inbounds Hᶜ(k::Int64, β::Array{Float64,1}) = [ cumsum( hᶜ(k,β)[end:-1:1] )[end:-1:1]; 0]
      e_o = findall( E .!= [ [] ] )
      e_o_f(s) = intersect( findall( obs .<= s ), e_o)
      return Model( data, obs, last_inc, incs, E, C, f, Hᵉ, Hᶜ, e_o, e_o_f)
end

# Random walk Metropolis-Hastings algorithm in the unit interval at log scale
@inline function metr_hast_chain_unif_a_b_interv_logscale_justend( lf::Function, l::Int64, a::Float64, b::Float64, σ::Float64)
    log_π_σ = lf(σ)
    for i in 1:l
        σ_prop = rand( Uniform(a,b) )
        log_π_σ_prop = lf(σ_prop)
        u = rand(Uniform())
        if log(u) < min( log_π_σ_prop - log_π_σ , 0. ) # symmetric q(x,y) AS IT IS gaussian random walk
            σ = σ_prop
            log_π_σ =log_π_σ_prop
        end
    end
    return σ
end

mi_logistic(x::Real,a::Real,b::Real) = (b + a*exp(-x))*inv(exp(-x) + one(x))
# Random walk Metropolis-Hastings algorithm at log scale
@inline function metr_hast_chain_randwalk_a_b_interv_logscale_justend( lf::Function, l::Int64, var::Float64, a::Float64, b::Float64, σ::Float64)
    log_π_σ = lf(σ)
    for i in 1:l
        σ_prop = mi_logistic( rand( Normal( log(σ-a) - log(b-σ), var) ), a, b)
        log_π_σ_prop = lf(σ_prop)
        log_q_given_σ = - log( σ_prop - a ) - log( b - σ_prop ) # + log(b-a) cancels out
        log_q_given_σ_prop = - log( σ - a ) - log( b - σ ) # + log(b-a) cancels out
        u = rand(Uniform())
        if log(u) < min( log_π_σ_prop - log_π_σ + log_q_given_σ_prop - log_q_given_σ, 0. ) # symmetric q(x,y) AS IT IS gaussian random walk
            σ = σ_prop
            log_π_σ = log_π_σ_prop
        end
    end
    return σ
end

# Random walk Metropolis-Hastings algorithm at log scale
@inline function metr_hast_chain_gauss_randwalk_logscale_justend( lf::Function, l::Int64, var::Float64, σ::Float64=0.0)
    log_π_σ = lf(σ)
    for i in 1:l
        σ_prop = rand(Normal(σ, var))
        log_π_σ_prop = lf(σ_prop)
        u = rand(Uniform())
        if log(u) < min( log_π_σ_prop - log_π_σ , 0. ) # symmetric q(x,y) AS IT IS gaussian random walk
            σ = σ_prop
            log_π_σ = log_π_σ_prop
        end
    end
    return σ
end

mi_LogNormal( m::Float64, v::Float64) = LogNormal( log(m) - 0.5*log( v/m^2 + 1. ), ( log( v/m^2 + 1. ) )^0.5 )

@inline function clayton_corm_lognorm_2dim_sampler(  δ::Float64, σ²::Float64)
    loc_u = rand( Uniform() )
    if loc_u < 1.0/2.0
        return [ rand( mi_LogNormal( 1.0*(1.0-δ) + δ, σ²)), rand( mi_LogNormal( 0.01*(1.0-δ) + δ, σ²)) ]
    else
        return [ rand( mi_LogNormal( 0.01*(1.0-δ) + δ, σ²)), rand( mi_LogNormal( 1.0*(1.0-δ) + δ, σ²)) ]
    end
end

# Uniform (0,1) Metropolis-Hastings algorithm at log scale
@inline function metr_hast_chain_unif_LogNorm_delta_logscale_justend( lf::Function,
        δ::Float64, V::Array{Array{Float64,1},1}, LogNorm_v::Float64, l_mc::Int64, l::Int64)
    log_π_V = lf(V)
    for i in 1:l
        δ_prop = rand( Uniform() )
        V_prop = [ clayton_corm_lognorm_2dim_sampler( δ_prop, LogNorm_v) for j in 1:l_mc ]
        log_π_V_prop = lf(V_prop)
        u = rand(Uniform())
        if log(u) < min( log_π_V_prop - log_π_V , 0.) # symmetric q(x,y) AS IT IS UNIFORM
            δ  = δ_prop
            V = V_prop
            log_π_V = log_π_V_prop
        end
    end
    return δ, V
end

mi_logistic(x::Real,a::Real,b::Real) = (b + a*exp(-x))*inv(exp(-x) + one(x))

# Uniform (0,1) Metropolis-Hastings algorithm at log scale
@inline function metr_hast_chain_randwalk_LogNorm_delta_logscale_justend( lf::Function,
        δ::Float64, V::Array{Array{Float64,1},1}, LogNorm_v::Float64, l_mc::Int64, l::Int64,
        var::Float64)
    log_π_V = lf(V)
    for i in 1:l
        δ_prop =  logistic( rand( Normal( log(δ) - log(1.0-δ), var) ) )
        V_prop = [ clayton_corm_lognorm_2dim_sampler( δ_prop, LogNorm_v) for j in 1:l_mc ]
        log_π_V_prop = lf(V_prop)
        log_q_given_δ = - log( δ_prop ) - log( 1.0 - δ_prop )
        log_q_given_δ_prop = - log( δ ) - log( 1.0 - δ )
        u = rand(Uniform())
        if log(u) < min( log_π_V_prop - log_π_V + log_q_given_δ_prop - log_q_given_δ, 0.)
            δ  = δ_prop
            V = V_prop
            log_π_V = log_π_V_prop
        end
    end
    return δ, V
end

function ψ_gamma( λ::Float64, α::Float64)
    return log(1.0 + λ/α)
end

function ψ( λ::Array{Float64,1}, α::Float64, V::Array{Array{Float64,1},1})
    return mean( [ ψ_gamma( sum( λ.*V[i] ), α) for i in 1:length(V)] )
end

function ψ_dif( λ::Array{Float64,1}, q::Array{Float64,1}, α::Float64, V::Array{Array{Float64,1},1})
    return mean( [ ψ_gamma( sum( (λ+q).*V[i]), α) - ψ_gamma(  sum( λ.*V[i]), α)  for i in 1:length(V) ] )
end

#function to calculate the product similar to the binomial theorem Π_{i ∈ I}(1-aᵢ) that I have in my likelihoods
# allowing for multiplicities 1 and 2
@inline function simil_binom_prod( q::Array{Float64}, i::Int64, β::Array{Float64,1}, α::Float64, v::Array{Float64,1}, model::Model)
    if length( model.E[i]) ==  1
        return log(  ( 1.0 + dot( v./α, [ model.f[j]( model.data.Z[model.E[i][1]], β) for j in 1:length(q) ] + q ) ) / ( 1.0 +
                            dot( v./α, q ) ) )
    elseif length( model.E[i]) ==  2
        q₁ = dot( v./α, [ model.f[j]( model.data.Z[model.E[i][1]], β) for j in 1:length(q) ] )
        q₂ = dot( v./α, [ model.f[j]( model.data.Z[model.E[i][2]], β) for j in 1:length(q) ] )
        q₀ = dot( v./α, q )
        return log( 1.0 + q₁*q₂ / ( 1.0 + q₁ + q₂ + 2.0*q₀ + q₁*q₀ + q₂*q₀ + q₀^2.0) )
    else
        println( length( model.E[i]) )
    end
end

#function to calculate the product similar to the binomial theorem Π_{i ∈ I}(1-aᵢ) that I have in my likelihoods
# allowing for multiplicity 1, length(model.E[i])=1 for all i
# @inline function simil_binom_prod( q::Array{Float64}, i::Int64, β::Array{Float64,1}, α::Float64, v::Array{Float64,1}, model::Model)
#     return log(  ( 1.0 + dot( v./α, [ model.f[j]( model.data.Z[model.E[i][1]], β) for j in 1:length(q) ] + q ) ) / ( 1.0 +
#                             dot( v./α, q ) ) )
# end

@inline function simil_binom_prod( q::Array{Float64}, i::Int64, β::Array{Float64,1}, α::Float64, v::Array{Float64,1}, model::Model)
    return log( 1.0 + ( 1.0/(α + dot(v,q)) ) * dot( v, [ model.f[j]( model.data.Z[model.E[i][1]], β) for j in 1:length(q) ] ) )
end

#function to calculate the product similar to the binomial theorem Π_{i ∈ I}(1-aᵢ) that I have in my likelihoods
# @inline function simil_binom_prod( q::Array{Float64}, i::Int64, β::Array{Float64,1}, α::Float64, v::Array{Float64,1}, model::Model)
#     return mapreduce( S -> (-1)^length(S) * ( ψ_gamma( sum( v.*( q + [ sum( [ model.f[i]( model.data.Z[l], β) for l in S]) for i in 1:length(q) ] ) ), α) - ψ_gamma( sum( v.*q ), α) ),
#               +, collect( subsets(model.E[i]) )[2:end], init = 0.0 )
# end
function e_contrib_func( q::Array{Float64,1}, i::Int64, β::Array{Float64,1}, α::Float64, V::Array{Array{Float64,1},1}, model::Model)
    r = mean( [ simil_binom_prod( q, i, β, α, V[j], model) for j in 1:length(V)])
    if r < 0.0
        println( "r<0!!!")
        global loc_q = q
        global loc_i = i
        global loc_β = β
        global loc_α = α
        global loc_V = V
    end
    return r
end

@inline function model_loglikelihood( β::Array{Float64,1}, A::Float64, α::Float64, V::Array{Array{Float64,1},1}, model::Model)
  obs = [0.; model.obs]
  Hᶜ = [ model.Hᶜ(i,β) for i in 1:length(model.f) ]
  Hᵉ = [ model.Hᵉ(i,β) for i in 1:length(model.f) ]
  return  - A * mapreduce( j -> ( obs[j+1] - obs[j] ) * ψ( [ Hᶜ[i][j] + Hᵉ[i][j] for i in 1:length(model.f) ], α, V), +, 1:length(model.obs), init = 0.0 ) +
          mapreduce( j -> log( e_contrib_func( [ Hᶜ[i][j] + Hᵉ[i][j+1] for i in 1:length(model.f) ], j, β, α, V, model) ), +, model.e_o, init = 0.0 ) +
          length(model.e_o)*log(A)
end

@inline function MCMC_comprisk_uniform( n::Int64, m::Int64, l_mc::Int64, A::Float64, α::Float64, model::Model,
    LogNorm_v::Float64, ch_init::Array{Float64},
    delta_init::Float64, a_interv::Float64, b_interv::Float64, β_run_map::Array{Float64}=Float64[], δ_run_map::Float64=0.0)
    chain_δ = zeros(Float64, n+1)
    chain_β = [ zeros(Float64, 1) for i in 1:(n+1) ]
    M = -Inf
    if ch_init != 1.0 .+ zeros(Float64,1)
        chain_β[1] = ch_init
        chain_δ[1] = delta_init
    else
        chain_β[1] = [ rand( Uniform(a_interv,b_interv) ) ]
        chain_δ[1] = rand( Uniform())
    end
    δ_map = δ_run_map
    β_map = β_run_map
    if δ_run_map == 0.0
        δ_map = chain_δ[1]
        β_map = chain_β[1]
    end
    loc_V = [ clayton_corm_lognorm_2dim_sampler( chain_δ[1], LogNorm_v) for j in 1:l_mc ]
    println("MCMC run has started.")
    prog = Progress( n, dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50)
    for i in 1:n
        @inbounds lf_pseud_marg(V) = model_loglikelihood( chain_β[i], A, α, V, model)
        @inbounds chain_δ[i+1], loc_V =  metr_hast_chain_unif_LogNorm_delta_logscale_justend( lf_pseud_marg, chain_δ[i], loc_V, LogNorm_v, l_mc, m)
        @inbounds lf_β_nocox(b) =  model_loglikelihood( [ b ], A, α, loc_V, model) # + logpdf( Beta(1.7,1.7), 0.5*b ) # logpdf( TruncatedNormal( 1.0, 1.0, 0.0, 2.0), b)
        @inbounds chain_β[i+1][1] =  metr_hast_chain_unif_a_b_interv_logscale_justend( lf_β_nocox, m, a_interv, b_interv, chain_β[i][1])
        M_temp = model_loglikelihood( chain_β[i+1], A, α, loc_V, model) # + logpdf( Beta(1.7,1.7), 0.5*chain_β[i+1][1])
        if M < M_temp
            @inbounds β_map = chain_β[i+1]
            @inbounds δ_map = chain_δ[i+1]
            M = M_temp
        end
        next!(prog)
    end
    return chain_δ, chain_β, δ_map, β_map
end

@inline function MCMC_comprisk_randwalk( n::Int64, m::Int64, l_mc::Int64, A::Float64, α::Float64, model::Model,
    LogNorm_v::Float64, δ_var::Float64, β_var::Float64, ch_init::Array{Float64}, delta_init::Float64, a_interv::Float64, b_interv::Float64, β_run_map::Array{Float64}=Float64[], δ_run_map::Float64=0.0)
    chain_δ = zeros(Float64, n+1)
    chain_β = [ zeros(Float64, 1) for i in 1:(n+1) ]
    M = -Inf
    if ch_init != 1.0 .+ zeros(Float64,1)
        chain_β[1] = ch_init
        chain_δ[1] = delta_init
    else
        chain_β[1] = [ rand( Uniform(a_interv,b_interv) ) ]
        chain_δ[1] = rand( Uniform())
    end
    δ_map = δ_run_map
    β_map = β_run_map
    if δ_run_map == 0.0
        δ_map = chain_δ[1]
        β_map = chain_β[1]
    end
    loc_V = [ clayton_corm_lognorm_2dim_sampler( chain_δ[1], LogNorm_v) for j in 1:l_mc ]
    println("MCMC run has started.")
    prog = Progress( n, dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50)
    for i in 1:n
        @inbounds lf_pseud_marg(V) = model_loglikelihood( chain_β[i], A, α, V, model)
        @inbounds chain_δ[i+1], loc_V = metr_hast_chain_randwalk_LogNorm_delta_logscale_justend( lf_pseud_marg, chain_δ[i], loc_V, LogNorm_v, l_mc, m, δ_var)
        @inbounds lf_β_nocox(b) =  model_loglikelihood( [ b ], A, α, loc_V, model) + logpdf( Beta(1.7,1.7), 0.5*b ) # logpdf( TruncatedNormal( 1.0, 1.0, 0.0, 2.0), b)
        @inbounds chain_β[i+1][1] =  metr_hast_chain_randwalk_a_b_interv_logscale_justend( lf_β_nocox, m, β_var, a_interv, b_interv, chain_β[i][1])
        M_temp = model_loglikelihood( chain_β[i+1], A, α, loc_V, model) # + logpdf( Beta(1.7,1.7), 0.5*chain_β[i+1][1])
        if M < M_temp
            @inbounds β_map = chain_β[i+1]
            @inbounds δ_map = chain_δ[i+1]
            M = M_temp
        end
        next!(prog)
    end
    return chain_δ, chain_β, δ_map, β_map
end

@inline function surv( t::Float64, z_new::Array{Float64,1}, β::Array{Float64,1}, A::Float64, α::Float64, V::Array{Array{Float64,1},1}, model::Model)
    obs = [0.; model.obs; Inf]
    Hᶜ = [ model.Hᶜ(i,β) for i in 1:length(model.f) ]
    Hᵉ = [ model.Hᵉ(i,β) for i in 1:length(model.f) ]
    H = [ [ model.Hᶜ(i,β)[j] + model.Hᵉ(i,β)[j] for i in 1:length(model.f) ] for j in 1:(length(model.obs) + 1) ]
    H_shift = [ [ model.Hᶜ(i,β)[j] + model.Hᵉ(i,β)[j+1] for i in 1:length(model.f) ] for j in 1:length(model.obs) ]
    f_new = [ model.f[i]( z_new, β) for i in 1:length(model.f) ]
    return ℯ^( - A * mapreduce( j -> ( min( t, obs[j+1]) - obs[j] ) * ( obs[j] <= t ) *
                                 ψ_dif( H[j], f_new, α, V), +, 1:(length(model.obs) + 1), init = 0.0 ) ) *
          mapreduce( j -> e_contrib_func( f_new + H_shift[j], j, β, α, V, model)/
                          e_contrib_func(  H_shift[j], j, β, α, V, model), *, model.e_o_f(t), init = 1.0 )
end

@inline function surv_vec_delta_prev( z_new::Array{Float64,1}, β::Array{Float64,1}, A::Float64, α::Float64, δ::Float64, LogNorm_v::Float64, l_mc::Int64, t_f::Float64, mesh::Float64, model::Model)
    V = [ clayton_corm_lognorm_2dim_sampler( δ, LogNorm_v) for j in 1:l_mc ]
    return @showprogress [ surv( j, z_new, β, A, α, V, model)  for j in collect(0.:1.0/mesh:t_f) ]
end

function mi_find(v::Array{Float64,1},e::Float64)
    index = 1
    while v[index] <= e
        index += 1
    end
    return index-1
end

@inline function surv_vec_delta( z_new::Array{Float64,1}, β::Array{Float64,1}, A::Float64, α::Float64, δ::Float64, LogNorm_v::Float64, l_mc::Int64, t_f::Float64, mesh::Float64, model::Model)
    t_vec = collect(0.:1.0/mesh:t_f)
    s_vec = zeros(Float64,length(t_vec))
    s_vec[1] = 1.0
    V = [ clayton_corm_lognorm_2dim_sampler( δ, LogNorm_v) for j in 1:l_mc ]
    obs = [model.obs; Inf]
    H = [ [ model.Hᶜ(i,β)[j] + model.Hᵉ(i,β)[j] for i in 1:length(model.f) ] for j in 1:(length(model.obs) + 1) ]
    H_shift = [ [ model.Hᶜ(i,β)[j] + model.Hᵉ(i,β)[j+1] for i in 1:length(model.f) ] for j in 1:length(model.obs) ]
    f_new = [ model.f[i]( z_new, β) for i in 1:length(model.f) ]
    obs_ind_run = 1
    exp_part_run = 1.0
    prod_part_run = 1.0
    for j in 2:length(t_vec)
        obs_ind_t = mi_find( obs[(obs_ind_run+1):end], t_vec[j]) + obs_ind_run
        if obs_ind_t == obs_ind_run
            exp_part_run = exp_part_run*exp( - A * ( t_vec[j] - t_vec[j-1] ) * ψ_dif( H[obs_ind_run], f_new, α, V))
        else
            exp_part_run = exp_part_run*exp( - A * ( obs[obs_ind_run+1] - t_vec[j-1] ) * ψ_dif( H[obs_ind_run], f_new, α, V))
            for i in (obs_ind_run+1):(obs_ind_t-1)
                exp_part_run = exp_part_run*exp( - A * (obs[i+1]-obs[i]) * ψ_dif( H[i], f_new, α, V))
                prod_part_run = prod_part_run*e_contrib_func( f_new + H_shift[i], i, β, α, V, model)/e_contrib_func( H_shift[i], i, β, α, V, model)
            end
            obs_ind_run = obs_ind_t
            exp_part_run = exp_part_run*exp( - A * ( t_vec[j] - obs[obs_ind_run] ) * ψ_dif( H[obs_ind_run], f_new, α, V))
            prod_part_run = prod_part_run*e_contrib_func( f_new + H_shift[obs_ind_run], obs_ind_run, β, α, V, model)/e_contrib_func( H_shift[obs_ind_run], obs_ind_run, β, α, V, model)
        end
        s_vec[j] = exp_part_run * prod_part_run
    end
    return s_vec
end

# Kaplan-Meier estiamtor function
function KMest(tiempo, muerte)
    t_KMest = sort( unique( tiempo ) )
    n_KMest = [count( i->(i >= j), tiempo) for j in t_KMest]
    d_KMest = [sum( muerte[ findall(j->j==i,tiempo ) ] ) for i in t_KMest]
    s_KMest = cumprod( 1 .- d_KMest./n_KMest )
    return t_KMest, s_KMest
end

# Kaplan-Meier estiamtor function
function make_plot_vec( t, s)
    if t[1] != 0.0
        t = [ 0.0; t ]
        s = [ 1.0; s ]
    end
    t_plot = [ t[1] ]
    s_plot = [ s[1] ]
    for i in 2:( length(t))
        push!( t_plot,  t[i], t[i] )
        push!( s_plot,  s[i-1], s[i] )
    end
    return t_plot, s_plot
end

function posterior( x::Array{Float64}, A::Float64, α::Float64, model::Model, LogNorm_v::Float64=0.1, l_mc::Int64=1000)
    V = [ clayton_corm_lognorm_2dim_sampler( x[1], LogNorm_v) for j in 1:l_mc ]
    return model_loglikelihood( x[2:2], A, α, V, model)
end

function cred_band_mat(S::Array{Float64,2},n::Int64)
    for j in 1:n
        b = zeros(Float64, size(S)[2])
        for i in 1:251
         b[i] = maximum(S[:,i])-minimum(S[:,i])
        end
        b_max_index = sortperm(b)[end]
        s_out_index = findmax( abs.( S[:,b_max_index] .- mean(S[b_max_index]) ) )[2]
        S = S[1:end .!= s_out_index, 1:end]
    end
    return S
end


text:
using Neuroimaging, DataDeps, Test
# Tests to ensure datasets are available and correct

# Basic biosemi data file
data_path = joinpath(datadep"BioSemiTestFiles", "Newtest17-2048.bdf")
s = read_SSR(data_path)
@test samplingrate(s) == 2048u"Hz"
@test length(channelnames(s)) == 16
@test length(s.triggers["Index"]) == length(s.triggers["Code"])
@test length(s.triggers["Code"]) == length(s.triggers["Duration"])

# Basic biosemi data file
data_path = joinpath(
    datadep"ExampleSSR",
    "Neuroimaging.jl-example-data-master",
    "neuroimaingSSR.bdf",
)
s = read_SSR(data_path)
@test samplingrate(s) == 8192u"Hz"
@test length(channelnames(s)) == 7
@test length(s.triggers["Index"]) == length(s.triggers["Code"])
@test length(s.triggers["Code"]) == length(s.triggers["Duration"])


text:
@inline function ExaTronKernels.ddot(n::Int,dx,incx::Int,
                      dy,incy::Int,
                      tx)
    # Currently, all threads compute the same dot product,
    # hence, no sync_threads() is needed.
    # For very small n, we may want to gauge how much gains
    # we could get by run it in parallel.

    v = 0.0
    @inbounds for i=1:n
        v += dx[i]*dy[i]
    end
    @synchronize
    return v
end


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

mutable struct CustomDnsConfigPropertiesFormat <: SwaggerModel
    fqdn::Any # spec type: Union{ Nothing, String } # spec name: fqdn
    ipAddresses::Any # spec type: Union{ Nothing, Vector{String} } # spec name: ipAddresses

    function CustomDnsConfigPropertiesFormat(;fqdn=nothing, ipAddresses=nothing)
        o = new()
        validate_property(CustomDnsConfigPropertiesFormat, Symbol("fqdn"), fqdn)
        setfield!(o, Symbol("fqdn"), fqdn)
        validate_property(CustomDnsConfigPropertiesFormat, Symbol("ipAddresses"), ipAddresses)
        setfield!(o, Symbol("ipAddresses"), ipAddresses)
        o
    end
end # type CustomDnsConfigPropertiesFormat

const _property_map_CustomDnsConfigPropertiesFormat = Dict{Symbol,Symbol}(Symbol("fqdn")=>Symbol("fqdn"), Symbol("ipAddresses")=>Symbol("ipAddresses"))
const _property_types_CustomDnsConfigPropertiesFormat = Dict{Symbol,String}(Symbol("fqdn")=>"String", Symbol("ipAddresses")=>"Vector{String}")
Base.propertynames(::Type{ CustomDnsConfigPropertiesFormat }) = collect(keys(_property_map_CustomDnsConfigPropertiesFormat))
Swagger.property_type(::Type{ CustomDnsConfigPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_CustomDnsConfigPropertiesFormat[name]))}
Swagger.field_name(::Type{ CustomDnsConfigPropertiesFormat }, property_name::Symbol) =  _property_map_CustomDnsConfigPropertiesFormat[property_name]

function check_required(o::CustomDnsConfigPropertiesFormat)
    true
end

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


text:
function longest_common_prefix(prefixa, prefixb)
    maxplength = min(sizeof(prefixa), sizeof(prefixb))
    maxplength == 0 && return ""
    idx = findfirst(i->(prefixa[i] != prefixb[i]), 1:maxplength)
    idx = idx === nothing ? maxplength : idx - 1
    prefixa[1:idx]
end

function skip_to_nl(str, idxend)
    while (idxend < sizeof(str)) && str[idxend] != '\n'
        idxend = nextind(str, idxend)
    end
    idxend > sizeof(str) ? prevind(str, idxend) : idxend
end

tostr(buf::IOBuffer) = _unescape_string(String(take!(buf)))

"""
parse_string_or_cmd(ps)

When trying to make an `INSTANCE` from a string token we must check for
interpolating opoerators.
"""
function parse_string_or_cmd(ps::ParseState, prefixed = false)
    sfullspan = ps.nt.startbyte - ps.t.startbyte
    sspan = 1 + ps.t.endbyte - ps.t.startbyte

    istrip = (kindof(ps.t) == Tokens.TRIPLE_STRING) || (kindof(ps.t) == Tokens.TRIPLE_CMD)
    iscmd = kindof(ps.t) == Tokens.CMD || kindof(ps.t) == Tokens.TRIPLE_CMD

    lcp = nothing
    exprs_to_adjust = []
    function adjust_lcp(expr, last = false)
        if isliteral(expr)
            push!(exprs_to_adjust, expr)
            str = valof(expr)
            (isempty(str) || (lcp !== nothing && isempty(lcp))) && return
            (last && str[end] == '\n') && return (lcp = "")
            idxstart, idxend = 2, 1
            while nextind(str, idxend) - 1 < sizeof(str) && (lcp === nothing || !isempty(lcp))
                idxend = skip_to_nl(str, idxend)
                idxstart = nextind(str, idxend)
                while nextind(str, idxend) - 1 < sizeof(str)
                    c = str[nextind(str, idxend)]
                    if c == ' ' || c == '\t'
                        idxend += 1
                    elseif c == '\n'
                        # All whitespace lines in the middle are ignored
                        idxend += 1
                        idxstart = idxend + 1
                    else
                        prefix = str[idxstart:idxend]
                        lcp = lcp === nothing ? prefix : longest_common_prefix(lcp, prefix)
                        break
                    end
                end
            end
            if idxstart != nextind(str, idxend)
                prefix = str[idxstart:idxend]
                lcp = lcp === nothing ? prefix : longest_common_prefix(lcp, prefix)
            end
        end
    end

    # there are interpolations in the string
    if prefixed != false || iscmd
        t_str = val(ps.t, ps)
        _val = istrip ? t_str[4:prevind(t_str, sizeof(t_str), 3)] : t_str[2:prevind(t_str, sizeof(t_str))]
        if iscmd
            _val = replace(_val, "\\\\" => "\\")
            _val = replace(_val, "\\`" => "`")
        else
            if endswith(_val, "\\\\")
                _val = _val[1:end - 1]
            end
            _val = replace(_val, "\\\"" => "\"")
        end
        expr = mLITERAL(sfullspan, sspan, _val, kindof(ps.t))
        if istrip
            adjust_lcp(expr)
            ret = EXPR(StringH, EXPR[expr], sfullspan, sspan)
        else
            return expr
        end
    else
        ret = EXPR(StringH, EXPR[], sfullspan, sspan)
        input = IOBuffer(val(ps.t, ps))
        startbytes = istrip ? 3 : 1
        seek(input, startbytes)
        b = IOBuffer()
        while !eof(input)
            c = read(input, Char)
            if c == '\\'
                write(b, c)
                write(b, read(input, Char))
            elseif c == '$'
                lspan = position(b)
                str = tostr(b)
                ex = mLITERAL(lspan + startbytes, lspan + startbytes, str, Tokens.STRING)
                push!(ret, ex)
                istrip && adjust_lcp(ex)
                startbytes = 0
                op = mOPERATOR(1, 1, Tokens.EX_OR, false)
                if peekchar(input) == '('
                    lparen = mPUNCTUATION(Tokens.LPAREN, 1, 1)
                    rparen = mPUNCTUATION(Tokens.RPAREN, 1, 1)
                    skip(input, 1)
                    ps1 = ParseState(input)

                    if kindof(ps1.nt) == Tokens.RPAREN
                        call = mUnaryOpCall(op, EXPR(InvisBrackets, EXPR[lparen, rparen]))
                        push!(ret, call)
                        skip(input, 1)
                    else
                        interp = @closer ps1 paren parse_expression(ps1)
                        call = mUnaryOpCall(op, EXPR(InvisBrackets, EXPR[lparen, interp, rparen]))
                        push!(ret, call)
                        seek(input, ps1.nt.startbyte + 1)
                    end
                    # Compared to flisp/JuliaParser, we have an extra lookahead token,
                    # so we need to back up one here
                elseif Tokenize.Lexers.iswhitespace(peekchar(input)) || peekchar(input) === '#'
                    push!(ret, mErrorToken(op, StringInterpolationWithTrailingWhitespace))
                else
                    pos = position(input)
                    ps1 = ParseState(input)
                    next(ps1)
                    if kindof(ps1.t) == Tokens.WHITESPACE
                        error("Unexpecte whitespace after \$ in String")
                    else
                        t = INSTANCE(ps1)
                    end
                    # Attribute trailing whitespace to the string
                    t = adjustspan(t)
                    call = mUnaryOpCall(op, t)
                    push!(ret, call)
                    seek(input, pos + t.fullspan)
                end
            else
                write(b, c)
            end
        end

        # handle last String section
        lspan = position(b)
        if b.size == 0
            ex = mErrorToken(Unknown)
        else
            str = tostr(b)
            if istrip
                str = str[1:prevind(str, prevind(str, sizeof(str), 2))]
                # only mark non-interpolated triple strings
                ex = mLITERAL(lspan + ps.nt.startbyte - ps.t.endbyte - 1 + startbytes, lspan + startbytes, str, length(ret) == 0 ? Tokens.TRIPLE_STRING : Tokens.STRING)
                adjust_lcp(ex, true)
            else
                str =  str[1:prevind(str, sizeof(str))]
                ex = mLITERAL(lspan + ps.nt.startbyte - ps.t.endbyte - 1 + startbytes, lspan + startbytes, str, Tokens.STRING)
            end
        end
        push!(ret, ex)

    end

    single_string_T = (Tokens.STRING, kindof(ps.t))
    if istrip
        if lcp !== nothing && !isempty(lcp)
            for expr in exprs_to_adjust
                for (i, a) in enumerate(ret.args)
                    if expr == a
                        ret.args[i].val = replace(valof(expr), "\n$lcp" => "\n")
                        break
                    end
                end
            end
        end
        # Drop leading newline
        if isliteral(ret.args[1]) && kindof(ret.args[1]) in single_string_T &&
                !isempty(valof(ret.args[1])) && valof(ret.args[1])[1] == '\n'
            ret.args[1] = dropleadlingnewline(ret.args[1])
        end
    end

    if (length(ret.args) == 1 && isliteral(ret.args[1]) && kindof(ret.args[1]) in single_string_T)
        ret = ret.args[1]
    end
    update_span!(ret)

    return ret
end

function adjustspan(x::EXPR)
    x.fullspan = x.span
    return x
end

dropleadlingnewline(x) = mLITERAL(x.fullspan, x.span, valof(x)[2:end], kindof(x))


text:
# http://sunlightlabs.github.io/datacommons
global OPEN_STATES_API = URI("http://openstates.org/")

# -------

function metadata_overview(auth::String, options...)
    args = Dict()

    sunlight_get(auth, OPEN_STATES_API, "/api/v1/metadata/", args; options...)
end

#get_metadata_overview(auth = "", options...) = get_metadata_overview(auth, options...)

# -------

function state_metadata(auth::String, state, options...)
    args = Dict()

    sunlight_get(auth, OPEN_STATES_API, "/api/v1/metadata/$state/", args; options...)
end

#get_state_metadata(state_str; auth = "", options...) = get_state_metadata(auth, state_str; options...)

# -------

function bill_search(auth::String; state=nothing, chamber=nothing, bill_id=nothing, bill_id__in=nothing,
                    q=nothing, search_window=nothing, updated_since=nothing, sponsor_id=nothing,
                    subject=nothing, bill_type=nothing, sort=nothing, page=nothing, per_page=nothing, options...)
    args = Dict()
    state != nothing && (args["state"] = state)
    chamber != nothing && (args["chamber"] = chamber)
    bill_id != nothing && (args["bill_id"] = bill_id)
    bill_id__in != nothing && (args["bill_id__in"] = bill_id__in)
    q != nothing && (args["q"] = q)
    search_window != nothing && (args["search_window"] = search_window)
    updated_since != nothing && (args["updated_since"] = updated_since)
    sponsor_id != nothing && (args["sponsor_id"] = sponsor_id)
    subject != nothing && (args["subject"] = subject)
    bill_type != nothing && (args["type"] = bill_type)
    sort != nothing && (args["sort"] = sort)
    page != nothing && (args["page"] = page)
    per_page != nothing && (args["per_page"] = per_page)

    sunlight_get(auth, OPEN_STATES_API, "/api/v1/bills/", args; options...)
end

#bill_search(; auth = "", options...) = bill_search(auth, entity_id; options...)

# -------

function bill_detail(auth::String; open_states_id=nothing, state=nothing, bill_id=nothing, cycle = nothing, options...)
    args = Dict()

    if open_states_id !=nothing
        sunlight_get(auth, OPEN_STATES_API, "/api/v1/bills/$open_states_id/", args; options...)
    else
        sunlight_get(auth, OPEN_STATES_API, "/api/v1/bills/$state/$cycle/$bill_id/", args; options...)
    end
end

#bill_search(entity_id; auth = "", options...) = bill_search(auth, entity_id; options...)

# -------

function legislator_search(auth::String; state=nothing, first_name=nothing, last_name=nothing, chamber=nothing,
                    active=nothing, term=nothing, district=nothing, party=nothing, options...)
    args = Dict()
    state != nothing && (args["state"] = state)
    first_name != nothing && (args["first_name"] = first_name)
    last_name != nothing && (args["last_name"] = last_name)
    chamber != nothing && (args["chamber"] = chamber)
    active != nothing && (args["active"] = active)
    term != nothing && (args["term"] = term)
    district != nothing && (args["district"] = district)
    party != nothing && (args["party"] = party)

    sunlight_get(auth, OPEN_STATES_API, "/api/v1/legislators/", args; options...)
end

# -------

function legislator_detail(auth::String, leg_id=nothing, options...)
    args = Dict()

    sunlight_get(auth, OPEN_STATES_API, "/api/v1/legislators/$leg_id/", args; options...)
end

# -------

function geo_lookup(auth::String; lat=nothing, long=nothing, options...)
    args = Dict()
    lat != nothing && (args["lat"] = lat)
    long != nothing && (args["long"] = long)

    sunlight_get(auth, OPEN_STATES_API, "/api/v1/legislators/geo/", args; options...)
end

# -------

function committee_search(auth::String; committee=nothing, subcommittee=nothing, chamber = nothing, state=nothing, options...)
    args = Dict()
    committee != nothing && (args["committee"] = committee)
    subcommittee != nothing && (args["subcommittee"] = subcommittee)
    chamber != nothing && (args["chamber"] = chamber)
    state != nothing && (args["state"] = state)

    sunlight_get(auth, OPEN_STATES_API, "/api/v1/committees/", args; options...)
end

# -------

function committee_detail(auth::String, committee_id=nothing, options...)
    args = Dict()

    sunlight_get(auth, OPEN_STATES_API, "/api/v1/committees/$committee_id/", args; options...)
end

# -------

function event_search(auth::String; state=nothing, event_type=nothing, format=nothing, options...)
    args = Dict()
    state != nothing && (args["state"] = state)
    event_type != nothing && (args["type"] = event_type)
    format != nothing && (args["format"] = format)

    sunlight_get(auth, OPEN_STATES_API, "/api/v1/events/", args; options...)
end

# -------

function event_detail(auth::String, event_id=nothing, options...)
    args = Dict()

    sunlight_get(auth, OPEN_STATES_API, "/api/v1/events/$event_id/", args; options...)
end

# -------

function district_search(auth::String; state=nothing, chamber=nothing, options...)
    args = Dict()
    if chamber != nothing
        sunlight_get(auth, OPEN_STATES_API, "/api/v1/districts/$state/$chamber/", args; options...)
    else
        sunlight_get(auth, OPEN_STATES_API, "/api/v1/districts/$state/", args; options...)
    end
end

# -------

function district_boundary_lookup(auth::String; boundary_id=nothing, options...)
    args = Dict()

    sunlight_get(auth, OPEN_STATES_API, "/api/v1/districts/boundary/$boundary_id/", args; options...)
end


text:
# TODO:
#   - Parsing factorials

#  https://github.com/sympy/sympy/blob/master/sympy/parsing/mathematica.py
"""
    mathematica2julia(s, additional_translations=[])

Converts a Mathematica expression string into a Julia string
"""
function mathematica2julia(s, additional_translations=[])
    # TODO: do creation of translations automatically??
    # TODO: test additional inputs

    checkinput(s)                       # Checking bracket pairs

    # Preprocessing
    s = applyrule(s, "implicitmult")    # uncover '*' hiding behind a whitespace
    s = applyrule(s, "whitespace")      # remove whitespace(s)
    s = applyrule(s, "add*_1")          # add omitted '*' character
    s = applyrule(s, "add*_2")          # add omitted '*' character

    s = replace(s, "!"=>"")

    try
        s = convertfunction(s)              # translate function
    catch
        @warn "Parse error on:\n$s"
        s = ""
    end

    s = applyrule(s, "Pi")              # Replace Pi with π
    s = applyrule(s, "I")               # Replace I with im

    return Meta.parse(s)
    # return s # Return a string or expression??
end

################################################################################
##                                                                            ##
#                     Creating a Translation Table                             #
#                                                                              #
################################################################################

Base.@kwdef struct JuliaTranslation
    fn::String = ""
    args::Vector{String} = []
    pat::Regex
end

"""
    JuliaTranslation(p::Pair)

Takes a Pair type which maps a Mathematica function to a Julia function and
creates a translation object to aid with function mapping.

Sample input: Log[x,y] => log(y,x)

# Fields:
 * `fn` Julia function name (eg, `log(y,x)`)
 * `args` arguments of the Julia function (eg, `["x", "y"]``)
 * `pat` pattern to match the arguments (eg, r"(?:\A|(?<=[^a-zA-Z]))(?:(x|y))(?=[^a-zA-Z])"))
"""
function JuliaTranslation(p::Pair)
    # s,j = p
    m = match(matfuncregex, p[1])
    args, _ = get_args(p[1], m)

    # convert '*x' to '\\*x' for regex
    re_args = replace.(args, "*"=>"\\*")
    xyz = "(?:(" * join(re_args, "|") * "))"

    # Template for argument pattern
    argpatregex = "(?:\\A|(?<=[^a-zA-Z]))$xyz(?=[^a-zA-Z])"
    pat = Regex(argpatregex)
    JuliaTranslation(p[2], args, pat)
end

# TODO: docs
"""
createtranslations(d::Dict)

Creates a dictionary where the keys are tuples of Mathematica functions and
their argument numbers (eg, "("Log", 2)"). The values of the dictionary are
JuliaTranslation types.
"""
function createtranslations(d::Dict)
    translations = Dict()
    for pair in d
        m = match(matfuncregex, pair[1])    # Match of Mathematica function
        jt = JuliaTranslation(pair)         # Create JuliaTranslation from pair
        key = (m.match, length(jt.args))    # Create a key from func and argnum
        push!(translations, key => jt)      # Push to dict
    end
    translations
end

################################################################################
##                                                                            ##
#                          Conversion Functions                                #
#                                                                              #
################################################################################

"""
    checkinput(s)

Checks to make sure that brackets are matched
"""
function checkinput(s)
    for bracket in [("[", "]"), ("{", "}"), ("(", ")")]
        if count(bracket[1], s) != count(bracket[2], s)
            @error "'$s' function form is invalid."
        end
    end
end

# TODO: Make this a mutating function
"""
    applyrule(s, rulename)

Applies a rule by name from a dictionary called rules
"""
function applyrule(s, rulename)
    pat, sub = rules[rulename]
    replace(s, pat => sub)
end

"""
    convertfunction(s)

Parse Mathematica function to SymPy one
"""
function convertfunction(s)
    scanned = ""                # converted string
    cur = 0                     # position cursor
    while true
        m = match(matfuncregex, s)

        if isnothing(m)
            # append the rest of string
            scanned *= s
            break
        end

        fm = m.match    # get Mathematica function name
        bgn = m.offset  # the start position of fm function

        # get arguments, and the end position of fm function
        args, stop = get_args(s, m)

        # convert Mathematica function to SymPy one
        s = convert_one_function(s, fm, args, bgn, stop)

        cur = bgn               # update cursor
        scanned *= s[1:cur-1]   # append converted part
        s = s[cur:end]          # shrink s
    end

    return scanned
end

# Create a new enumerate function which takes a variable start number
Base.enumerate(iter, start=1) = zip(Iterators.countfrom(start), iter)

# TODO: simplify
# TODO: write tests
"""
    get_args(s, m)

Get arguments of a Mathematica function
"""
function get_args(s, m)
    anc = m.offset + length(m.match) + 1 # pointing the first letter of arguments
    square, curly = [], []              # stack for brakets
    args = String[]
    # args = ""

    # current cursor
    cur = anc
    j = 0
    for (i, c) in enumerate(s[anc:end], anc)
        # extract one argument
        if c == ',' && isempty(square) #&& (!isempty(curly))
            push!(args, s[cur:i-1])       # add an argument
            cur = i + 1                   # move cursor
        end

        # # handle list or matrix (for future usage)
        if c == '{'
            append!(curly, c)
        elseif c == '}'
            pop!(curly)
        end

        # # seek corresponding ']' with skipping irrevant ones
        if c == '['
            append!(square, c)
        elseif c == ']'
            if !isempty(square)
                pop!(square)
            else
                push!(args, s[cur:i-1])
                break
            end
        end
        j = i
    end

    # the next position to "]" bracket (the function end)
    func_end = j + 1

    args, func_end
end

"""
    convert_one_function(s, fm, args, bgn, stop)

Inputs:
 * `s`     The entire string to be matched
 * `fm`    The Mathematica function
 * `args`  A list of arguments
 * `bgn`   The location in the string, s, where the function, fm, begins
 * `stop`  The location in the string, s, where the function, fm, ends

 Converts the first functions it see into Julia function
"""
function convert_one_function(s, fm, args, bgn, stop)
    # Function is in translations with the number of arguments
    if (fm, length(args)) in keys(translations)
        key = (fm, length(args))

        # x, y,... model arguments
        x_args = translations[key].args

        # make CORRESPONDENCES between model arguments and actual ones
        d = Dict(k => v for (k, v) in zip(x_args, args))

    # with variable-length argument
    # TODO: Implements variable-length arguments, perhaps with splatting
    elseif (fm, '*') in keys(translations)
        key = (fm, '*')

        # x, y,..*args (model arguments)
        x_args = translations[key].args

        # make CORRESPONDENCES between model arguments and actual ones
        d = Dict()
        for (i, x) in enumerate(x_args)
            if x[0] == '*'
                d[x] = join(args[i:end], ",")
                break
            end
            d[x] = args[i]
        end
    else    # The function is not in translations
        @error "$fm with $(length(args)) arguments is not in translations pairs."
    end

    template = translations[key].fn    # template string of converted function
    pat = translations[key].pat        # regex pattern for x_args

    scanned = ""
    cur = 0
    while true
        m = match(pat, template)

        if isnothing(m)
            scanned *= template
            break
        end

        x = m.match         # get model argument
        xbgn = m.offset - 1 # get a start position of the model argument

        # add the corresponding actual argument
        scanned *= template[1:xbgn] * d[x]

        # update cursor to the end of the model argument
        cur = m.offset + length(m.match)

        # shrink template
        template = template[cur:end]
    end

    # update to swapped string
    return s[1:bgn-1] * scanned * s[stop+1:end]
end

################################################################################
##                                                                            ##
#                          Default Dictionaries and Rules                      #
#                                                                              #
################################################################################
replacements = Dict(
    " " => "",
    "{" => "[",
    "}" => "]",
)

# Default dictionary of rules
rules = Dict(
    # a single whitespace to "*"
     "implicitmult" => (r"""
            (?<=[a-zA-Z\d])     # a letter or a number
            \                   # a whitespace
            (?=[a-zA-Z\d])      # a letter or a number
            """x,
            "*"),

    # Remove whitespace
    "whitespace" => (" ", ""),

    # add omitted "*" character
    "add*_1" => (r"""
            (?<=[])\d])         # ], ) or a number
                                # ""
            (?=[(a-zA-Z])       # ( or a single letter
            """x,
            "*"),

    # add omitted "*" character (variable letter preceding)
    "add*_2" => (r"""
            (?<=[a-zA-Z])       # a letter
            \(                  # ( as a character
            (?=.)               # any characters
            """x,
            "*("),

    # convert "Pi" to "π"
    "Pi" => (r"""
    (?:
    \A|(?<=[^a-zA-Z])
    )
    Pi                  # "Pi" is 3.14159... in Mathematica
    (?=[^a-zA-Z])
    """x,
    "π"),

    # convert "I" to "im"
    "I" => (r"""
            (?:
            \A|(?<=[^a-zA-Z])
            )
            I                  # "Pi" is 3.14159... in Mathematica
            (?=[^a-zA-Z])
            """x,
            "im"),
    )

matfuncregex = r"""
(?:
\A|(?<=[^a-zA-Z])   # at the top or a non-letter
)
[A-Z][a-zA-Z\d]*    # Function
(?=\[)              # [ as a character
"""x

# list or matrix pattern (for future usage)
matargsregex = r"\{.*\}"

#  TODO: perhaps manually write many more correspondences
correspondences = Dict(
    "Sqrt[x]" => "sqrt(x)",
    "Exp[x]" => "exp(x)",
    "Log[x]" => "log(x)",
    "Log[x,y]" => "log(y,x)",
    "Log2[x]" => "log(x,2)",
    "Log10[x]" => "log(x,10)",
    "Mod[x,y]" => "mod(x,y)",
    "Max[*x]" => "Max(*x)",
    "Min[*x]" => "Min(*x)",
    "Pochhammer[x,y]" =>"rf(x,y)",
    "ArcTan[x,y]" =>"atan2(y,x)",
    "ExpIntegralEi[x]" => "Ei(x)",
    "SinIntegral[x]" => "Si(x)",
    "CosIntegral[x]" => "Ci(x)",
    "AiryAi[x]" => "airyai(x)",
    "AiryAiPrime[x]" => "airyaiprime(x)",
    "AiryBi[x]"  =>"airybi(x)",
    "AiryBiPrime[x]"  =>"airybiprime(x)",
    "LogIntegral[x]" =>"li(x)",
    "PrimePi[x]" => "primepi(x)",
    "Prime[x]" => "prime(x)",
    "PrimeQ[x]" => "isprime(x)",
);

# Adds trig functions to the correspondences programmatically
triprefixes = ["", "Arc"]
trifunctions = ["Sin", "Cos", "Tan", "Cot", "Sec", "Csc"]
trisuffixes = ["", "h"]
for p in triprefixes, f in trifunctions, s in trisuffixes
    fm = "$p$f$s[x]"
    fs = length(p) > 0 ? "a$(lowercase(f))$s(x)" : "$(lowercase(f))$s(x)"
    push!(correspondences, fm => fs)
end

################################################################################
##                                                                            ##
#                                   Testing                                    #
#                                                                              #
################################################################################

using Test

translations = createtranslations(correspondences)

@testset "Replacement Rules" begin
    @test applyrule("2 x", "implicitmult") == "2*x"
    @test applyrule("(1 + 2)", "whitespace") == "(1+2)"
    @test applyrule("(1+2)x", "add*_1") == "(1+2)*x"
    @test applyrule("x(1+2)", "add*_2") == "x*(1+2)"
    @test applyrule("Sin[Pi]", "Pi")  == "Sin[π]"

    s = "Sqrt[x] + 2 x + 3 (1 + 3) + (2 + 4) d"
    s = applyrule(s, "implicitmult")    # uncover '*' hiding behind a whitespace
    s = applyrule(s, "whitespace")      # remove whitespace(s)
    s = applyrule(s, "add*_1")          # add omitted '*' character
    s = applyrule(s, "add*_2")          # add omitted '*' character
    @test s == "Sqrt[x]+2*x+3*(1+3)+(2+4)*d"
end;

# @testset "Translation" begin
#     @test mathematica2julia("Sqrt[x]") == "sqrt(x)"
#     @test mathematica2julia("ArcTan[x,y]") == "atan2(y,x)"
#     @test mathematica2julia("Sin[Cos[x]]") == "sin(cos(x))"
#     @test mathematica2julia("Sin[Cos[Pi]]") == "sin(cos(π))"
#     @test mathematica2julia("Sqrt[x] + 2 x + 3 (1 + 3) + (2 + 4) d") == "sqrt(x)+2*x+3*(1+3)+(2+4)*d"
#     @test mathematica2julia("2*I*Sqrt[2]") == "2*im*sqrt(2)"
# end;
@testset "Translation" begin
    @test mathematica2julia("Sqrt[x]") == :(sqrt(x))
    @test mathematica2julia("ArcTan[x,y]") == :(atan2(y,x))
    @test mathematica2julia("Sin[Cos[x]]") == :(sin(cos(x)))
    @test mathematica2julia("Sin[Cos[Pi]]") == :(sin(cos(π)))
    @test mathematica2julia("Sqrt[x] + 2 x + 3 (1 + 3) + (2 + 4) d") == :(sqrt(x)+2*x+3*(1+3)+(2+4)*d)
    @test mathematica2julia("2*I*Sqrt[2]") == :(2*im*sqrt(2))
end;

################################################################################
##                                                                            ##
#                              Benchmarking                                    #
#                                                                              #
################################################################################

using BenchmarkTools
using SymPy
using DelimitedFiles

const sympy_parsing_mathematica = SymPy.PyCall.pyimport("sympy.parsing.mathematica")
m2j(s::AbstractString, substitutions::Pair{<:AbstractString,<:AbstractString}...) =
SymPy.walk_expression(sympy_parsing_mathematica."mathematica"(s, Dict(substitutions...)))

m2j_subs = [r[1] => r[2] for r in  eachrow(readdlm("func_subs.csv", ';', String))]

function m2j(s::AbstractString; verbose=false)
	try
		return m2j(s, m2j_subs...)
	catch
		verbose && @warn "Couldn't parse:\n$(s)\n"
		return "Parse Error"
	end
end

##
@btime m2j("Csc[x + Sin[x]]") # 768.100 μs (964 allocations: 50.62 KiB)
@btime mathematica2julia("Csc[x + Sin[x]]") # 12.300 μs (131 allocations: 6.73 KiB)

@btime m2j("(2*I*Sqrt[2]*Log[Pi/4 - (I*x)/2, 2] * Sqrt[Sinh[x]])/Sqrt[I*Sinh[x]]") # 2.592 ms (3137 allocations: 175.34 KiB)
@btime mathematica2julia("(2*I*Sqrt[2]*Log[Pi/4 - (I*x)/2, 2] * Sqrt[Sinh[x]])/Sqrt[I*Sinh[x]]") # 28.300 μs (380 allocations: 20.55 KiB)

## Running on many functions
using JSON3

# Load a list of Expression Pairs
m2j_subs = Dict(r[1] => r[2] for r in  eachrow(readdlm("func_subs.csv", ';', String)))
translations = createtranslations(m2j_subs)

fname = joinpath("src", "rubitests.json")
tests_str = read(fname, String);
tests = JSON3.read(tests_str);

ints = [t.integrand for t in tests];    # Extract all the integrands
@show length(ints)

num = 1000
sample_exprs = rand(ints, num)
@time mathematica2julia.(sample_exprs);
@time m2j.(sample_exprs);

print("")

################################################################################
##                                                                            ##
#                              Experiments                                     #
#                                                                              #
################################################################################

m2j_subs = Dict(r[1] => r[2] for r in  eachrow(readdlm("func_subs.csv", ';', String)))
translations = createtranslations(m2j_subs)

@time jints =  mathematica2julia.(ints);    # 15x faster
# 12.475728 seconds (9.21 M allocations: 486.123 MiB, 2.47% gc time)

##
@time jints2 =  m2j.(ints);
# 196.246567 seconds (196.29 M allocations: 10.434 GiB, 3.02% gc time, 0.25% compilation time)
##

s = "(a + b*x)!^n*PolyGamma[0, 1 + a + b*x]"
mathematica2julia(s)
s = "Sec[a + b*x]^2*F[c, d, Tan[a + b*x], r, s]"
mathematica2julia(s)

##
["""("$(k[1])", $(k[2]))""" for k in keys(translations)] |> sort .|> println;

translations[("F", 5)]


text:
struct System
    _f::Function
    _raw::Vector{Vector{ComponentPoolIndex}}
    _scheme::Tuple
end
struct Kosmo
    _component_pool::ComponentPool #// Raw components entries
    _system_pool::Vector{System}
end
Kosmo() = Kosmo(ComponentPool(), Vector())

ComponentSet(k::Kosmo, v) = ComponentSet(k._component_pool, v)

function run_systems(k::Kosmo)
    for sys in k._system_pool
        for r in sys._raw
            # These are ordered
            args = [k._component_pool[x] for x in r]
            sys._f(args...)
        end
    end
end

function add!(k::Kosmo, cmp::Component)
    c = k._component_pool
    TT = typeof(cmp)
    # update system pool
    for sys in k._system_pool
        if TT in sys._scheme
            scheme = filter(x-> x!=T, sys._scheme)
            t = create_pair_combinations(create_combinations(scheme, k._component_pool), [cmp])
        end
        append!(sys._raw, collect.(t))
    end
    if TT in keys(c._raw)
        push!(c._raw[TT], cmp)
    else
        c._raw[TT] = [cmp]
    end
    return ComponentPoolIndex(TT, length(c._raw[TT]))
end


text:
#=
Define region of interest alignment
=#

struct ROI_Align
    crop_height::Int
    crop_width::Int
    extrapolation_value::Float64
    transform_fpcoor::Bool
end
#=
ROA_Align(crop_height, crop_width, extrapolation_value = 0.0f0, transform_fpcoor = true) = ROI_Align(crop_height, crop_width, extrapolation_value, transform_fpcoor)

function (c::ROA_Align)(feature_map, boxes, box_ind)
    x1 = boxes[:, 1]
    y1 = boxes[:, 2]
    x2 = boxes[:, 3]
    y2 = boxes[:, 4]

    img_height, img_width = Float.(size(feature_map)[1:2])

    if c.transform_fpcoor
        spacing_w = (x2 .- x1) ./ img_width
        spacing_h = (y2 .- y1) ./ img_height

        nx0 = (x1 .+ (spacing_w ./ 2.f0) .- 0.5f0) ./ (img_width .- 1.0f0)
        ny0 = (y1 .+ (spacing_h ./ 2.f0) .- 0.5f0) ./ (img_height .- 1.0f0)
        nw = spacing_w .* (c.crop_width .- 1.0f0) ./ (img_width .- 1.0f0)
        nh = spacing_h .* (c.crop_height .- 1.0f0) ./ (img_height .- 1.0f0)
        boxes = cat(ny0, nx0, ny0 .+ nh, nx0 .+ nw, dims = 2)

    else
        x1 = x1 ./ (img_width .- 1.0f0)
        x2 = x2 ./ (img_width .- 1.0f0)
        y1 = y1 ./ (img_height .- 1.0f0)
        y2 = y2 ./ (img_height .- 1.0f0)
        boxes = cat(y1, x1, y2, x2, dims = 2)
    end

    crop_and_resize(feature_map, boxes, box_ind, c.crop_height, c.crop_width)
end
=#


text:
@testset "Static Dense Modified Gram Schmidt" begin
	T = Float64
	n = 16
	A = rand(n,n)
	a = Vector{SVector{n,T}}([SVector{n,T}(A[:,i]) for i=1:n])
	gs_solver = ContactImplicitMPC.SDMGSSolver(n)
	ContactImplicitMPC.factorize!(gs_solver, A)
	ContactImplicitMPC.factorize!(gs_solver, a)
	ContactImplicitMPC.factorize!(gs_solver)

	b = rand(SVector{n,T})
	ContactImplicitMPC.qr_solve!(gs_solver, b)
	@test norm(A*gs_solver.xs - b, Inf) < 1e-10
	# @benchmark A\b

	n = 43
	m = 33
	A = rand(n,n)
	gs_solver = ContactImplicitMPC.SDMGSSolver(A)
	ContactImplicitMPC.factorize!(gs_solver, A)
	X = rand(n,m)
	B = rand(n,m)
	ContactImplicitMPC.qr_matrix_solve!(gs_solver, X, B)
	@test norm(A * X - B) < 1e-10
end


text:
export GridSpace

struct Region{D}
    mini::NTuple{D,Int}
    maxi::NTuple{D,Int}
end

"""
    Hood{D}
Internal struct for efficiently finding neighboring positions to a given position.
It contains pre-initialized neighbor cartesian indices and delimiters of when the
neighboring indices would exceed the size of the underlying array.
"""
struct Hood{D}
    whole::Region{D} # allowed values (only useful for non periodic)
    βs::Vector{CartesianIndex{D}} # neighborhood cartesian indices
end

# type P stands for Periodic and is a boolean
struct GridSpace{D,P,W} <: DiscreteSpace
    s::Array{Vector{Int},D}
    metric::Symbol
    hoods::Dict{Float64,Hood{D}}
    hoods_tuple::Dict{NTuple{D,Float64},Hood{D}}
    pathfinder::W
end

"""
    GridSpace(d::NTuple{D, Int}; periodic = true, metric = :chebyshev, pathfinder = nothing)
Create a `GridSpace` that has size given by the tuple `d`, having `D ≥ 1` dimensions.
Optionally decide whether the space will be periodic and what will be the distance metric
used, which decides the behavior of e.g. [`nearby_ids`](@ref).
The position type for this space is `NTuple{D, Int}`, use [`GridAgent`](@ref) for convenience.
In our examples we typically use `Dims{D}` instead of `NTuple{D, Int}` (they are equivalent).
Valid positions have indices in the range `1:d[i]` for the `i`th dimension.

`:chebyshev` metric means that the `r`-neighborhood of a position are all
positions within the hypercube having side length of `2*floor(r)` and being centered in
the origin position.

`:euclidean` metric means that the `r`-neighborhood of a position are all positions whose
cartesian indices have Euclidean distance `≤ r` from the cartesian index of the given
position.

`pathfinder`: Optionally provide an instance of [`Pathfinding.Pathfinder`](@ref) to enable
pathfinding capabilities based on the A* algorithm, see [Path-finding](@ref) in the docs.

An example using `GridSpace` is the [Forest fire](@ref).
"""
function GridSpace(
    d::NTuple{D,Int};
    periodic::Bool = true,
    metric::Symbol = :chebyshev,
    pathfinder::W = nothing,
    moore = nothing,
) where {D,W}
    s = Array{Vector{Int},D}(undef, d)
    if moore ≠ nothing
        @warn "Keyword `moore` is deprecated, use `metric` instead."
        metric = moore == true ? :chebyshev : :euclidean
    end
    for i in eachindex(s)
        s[i] = Int[]
    end

    # TODO: This is bad design. `AStar` should not be mentioned here,
    # nor any `Pathfinding` business. This file should be "pure".
    astar = pathfinder === nothing ? nothing : Pathfinding.AStar(d, periodic, pathfinder)

    return GridSpace{D,periodic,typeof(astar)}(
        s,
        metric,
        Dict{Float64,Hood{D}}(),
        Dict{NTuple{D,Float64},Hood{D}}(),
        astar,
    )
end

#######################################################################################
# %% Implementation of space API
#######################################################################################
function random_position(model::ABM{<:GridSpace})
    Tuple(rand(model.rng, CartesianIndices(model.space.s)))
end

function add_agent_to_space!(a::A, model::ABM{<:GridSpace,A}) where {A<:AbstractAgent}
    push!(model.space.s[a.pos...], a.id)
    return a
end

function remove_agent_from_space!(a::A, model::ABM{<:GridSpace,A}) where {A<:AbstractAgent}
    prev = model.space.s[a.pos...]
    ai = findfirst(i -> i == a.id, prev)
    deleteat!(prev, ai)
    return a
end

function move_agent!(a::A, pos::Tuple, model::ABM{<:GridSpace,A}) where {A<:AbstractAgent}
    remove_agent_from_space!(a, model)
    a.pos = pos
    add_agent_to_space!(a, model)
end

##########################################################################################
# %% Structures for collecting neighbors on a grid
##########################################################################################
# Most of the source code in this section comes from TimeseriesPrediction.jl, specifically
# the file github.com/JuliaDynamics/TimeseriesPrediction.jl/src/spatiotemporalembedding.jl
# It creates a performant envinroment where the cartesian indices of a given neighborhood
# of given radious and metric type are stored and re-used for each search.

Base.length(r::Region{D}) where {D} = prod(r.maxi .- r.mini .+ 1)

function Base.in(idx, r::Region{D}) where {D}
    @inbounds for φ in 1:D
        r.mini[φ] <= idx[φ] <= r.maxi[φ] || return false
    end
    return true
end

# This function initializes the standard cartesian indices that needs to be added to some
# index `α` to obtain its neighborhood
function initialize_neighborhood!(space::GridSpace{D}, r::Real) where {D}
    d = size(space.s)
    r0 = floor(Int, r)
    if space.metric == :euclidean
        # hypercube of indices
        hypercube = CartesianIndices((repeat([(-r0):r0], D)...,))
        # select subset of hc which is in Hypersphere
        βs = [β for β ∈ hypercube if LinearAlgebra.norm(β.I) ≤ r]
    elseif space.metric == :chebyshev
        βs = vec([CartesianIndex(a) for a in Iterators.product([(-r0):r0 for φ in 1:D]...)])
    else
        error("Unknown metric type")
    end
    whole = Region(map(one, d), d)
    hood = Hood{D}(whole, βs)
    space.hoods[float(r)] = hood
    return hood
end

function initialize_neighborhood!(space::GridSpace{D}, r::NTuple{D,Real}) where {D}
    @assert space.metric == :chebyshev "Cannot use tuple based neighbor search with the Euclidean metric."
    d = size(space.s)
    r0 = (floor(Int, i) for i in r)
    βs = vec([CartesianIndex(a) for a in Iterators.product([(-ri):ri for ri in r0]...)])
    whole = Region(map(one, d), d)
    hood = Hood{D}(whole, βs)
    push!(space.hoods_tuple, float.(r) => hood)
    return hood
end

"""
    grid_space_neighborhood(α::CartesianIndex, space::GridSpace, r::Real)

Return an iterator over all positions within distance `r` from `α` according to the `space`.

The only reason this function is not equivalent with `nearby_positions` is because
`nearby_positions` skips the current position `α`, while `α` is always included in the
returned iterator (because the `0` index is always included in `βs`).

This function re-uses the source code of TimeseriesPrediction.jl, along with the
helper struct `Hood` and generates neighboring cartesian indices on the fly,
reducing the amount of computations necessary (i.e. we don't "find" new indices,
we only add a pre-determined amount of indices to `α`).
"""
function grid_space_neighborhood(α::CartesianIndex, space::GridSpace, r::Real)
    hood = if haskey(space.hoods, r)
        space.hoods[r]
    else
        initialize_neighborhood!(space, r)
    end
    _grid_space_neighborhood(α, space, hood)
end

"""
    grid_space_neighborhood(α::CartesianIndex, space::GridSpace, r::NTuple)

Return an iterator over all positions within distances of the tuple `r`, from `α`
according to the `space`. `r` must have as many elements as `space` has dimensions.
For example, with a `GridSpace((10, 10, 10))` : `r = (1, 7, 9)`.
"""
function grid_space_neighborhood(
    α::CartesianIndex,
    space::GridSpace{D},
    r::NTuple{D,Real},
) where {D}
    hood = if haskey(space.hoods_tuple, r)
        space.hoods_tuple[r]
    else
        initialize_neighborhood!(space, r)
    end
    _grid_space_neighborhood(α, space, hood)
end

function _grid_space_neighborhood(
    α::CartesianIndex,
    space::GridSpace{D,true},
    hood,
) where {D}
    return ((mod1.(Tuple(α + β), hood.whole.maxi)) for β in hood.βs)
end

function _grid_space_neighborhood(
    α::CartesianIndex,
    space::GridSpace{D,false},
    hood,
) where {D}
    return Iterators.filter(x -> x ∈ hood.whole, (Tuple(α + β) for β in hood.βs))
end

grid_space_neighborhood(α, model::ABM, r) = grid_space_neighborhood(α, model.space, r)

##########################################################################################
# %% Extend neighbors API for spaces
##########################################################################################
function nearby_ids(pos::ValidPos, model::ABM{<:GridSpace}, r = 1)
    nn = grid_space_neighborhood(CartesianIndex(pos), model, r)
    s = model.space.s
    Iterators.flatten((s[i...] for i in nn))
end

"""
    nearby_ids(pos, model::ABM{<:GridSpace}, r::Vector{Tuple{Int,UnitRange{Int}}})

Return an iterable of ids over specified dimensions of `space` with fine grained control
of distances from `pos` using each value of `r` via the (dimension, range) pattern.

**Note:** Only available for use with non-periodic chebyshev grids.

Example, with a `GridSpace((100, 100, 10))`: `r = [(1, -1:1), (3, 1:2)]` searches
dimension 1 one step either side of the current position (as well as the current
position) and the third dimension searches two positions above current.

For a complete tutorial on how to use this method, see [Battle Royale](@ref).
"""
function nearby_ids(
    pos::ValidPos,
    model::ABM{<:GridSpace},
    r::Vector{Tuple{Int,UnitRange{Int}}},
)
    dims = first.(r)
    vidx = []
    for d in 1:ndims(model.space.s)
        idx = findall(dim -> dim == d, dims)
        dim_range = isempty(idx) ? Colon() :
            bound_range(pos[d] .+ last(r[only(idx)]), d, model.space)
        push!(vidx, dim_range)
    end
    s = view(model.space.s, vidx...)
    Iterators.flatten(s)
end

function bound_range(unbound, d, space::GridSpace{D,false}) where {D}
    return range(max(unbound.start, 1), stop = min(unbound.stop, size(space)[d]))
end

function nearby_positions(pos::ValidPos, model::ABM{<:GridSpace}, r = 1)
    nn = grid_space_neighborhood(CartesianIndex(pos), model, r)
    Iterators.filter(!isequal(pos), nn)
end

#######################################################################################
# %% Further discrete space functions
#######################################################################################
function positions(model::ABM{<:GridSpace})
    x = CartesianIndices(model.space.s)
    return (Tuple(y) for y in x)
end

function ids_in_position(pos::ValidPos, model::ABM{<:GridSpace})
    return model.space.s[pos...]
end

###################################################################
# %% pretty printing
###################################################################
Base.size(space::GridSpace) = size(space.s)

function Base.show(io::IO, space::GridSpace{D,P}) where {D,P}
    pathfinder = isnothing(space.pathfinder) ? "" : ", pathfinder=$(space.pathfinder)"
    s = "GridSpace with size $(size(space)), metric=$(space.metric), periodic=$(P)$pathfinder"
    print(io, s)
end


text:
using JondeaUtils
import Test: @test, @test_throws
import SpecialFunctions: hankelh1, besselj
import Colors: RGB
import OffsetArrays: OffsetArray
import StaticArrays: SVector

@test flatten([[1,2],[3,4]]) == [1,2,3,4]
@test flatten([[1.0,2.0],[3.0,4.0]]) == [1.0,2.0,3.0,4.0]

@test flatten([1 3; 2 4]) == [1,2,3,4]
@test flatten([1.0 3.0; 2.0 4.0]) == [1.0,2.0,3.0,4.0]

@test real_cols_to_complex([1 2;3 4]) == reshape([1+2im, 3+4im],2,1)

@test rel_abs_error(1.0+0.0im,1.0+0.1im) ≈ 0.1
@test max_rel_abs_error([1.0im,2.0im],[1.1im,2.4im]) ≈ 0.2

@test offsetmap(x->x^2, -2:2) == OffsetArray([4,1,0,1,4],-2:2)

finite_diff(f,x,ε) = (f(x+ε) - f(x))/ε
@test diffhankelh1(0,2.0+1.0im) ≈ finite_diff(z->hankelh1(0,z),2.0+1.0im,1.0e-8im) atol=1.0e-6
@test diffhankelh1(1,2.0+1.0im) ≈ finite_diff(z->hankelh1(1,z),2.0+1.0im,1.0e-8im) atol=1.0e-6
@test diffbesselj(0,2.0+1.0im) ≈ finite_diff(z->besselj(0,z),2.0+1.0im,1.0e-8im) atol=1.0e-6
@test diffbesselj(1,2.0+1.0im) ≈ finite_diff(z->besselj(1,z),2.0+1.0im,1.0e-8im) atol=1.0e-6

@test fresnels(1.0)::Float64 ≈ 0.4382591473903548
@test fresnelc(1.0)::Float64 ≈ 0.7798934003768228

@test fresnels(1.0f0)::Float32 ≈ 0.4382591f0
@test fresnelc(1.0f0)::Float32 ≈ 0.7798934f0

@test all(reds(3) .≈ [RGB(1.0,0.4,0.4) RGB(0.8,0.3,0.3) RGB(0.6,0.2,0.2)])
@test all(greens(3) .≈  [RGB(0.4,1.0,0.4) RGB(0.3,0.8,0.3) RGB(0.2,0.6,0.2)])
@test all(blues(3) .≈ [RGB(0.4,0.4,1.0) RGB(0.3,0.3,0.8) RGB(0.2,0.2,0.6)])

@test reds(1) == reshape([RGB(0.8,0.3,0.3)],1,1)
@test greens(1) == reshape([RGB(0.3,0.8,0.3)],1,1)
@test blues(1) == reshape([RGB(0.3,0.3,0.8)],1,1)

# Test some known values
@test angle([1,0]) ≈ 0
@test angle([1.0,1.0]) ≈ π/4
@test angle([0,1]) ≈ π/2
@test angle([-1,-1]) ≈ -3π/4

# Test using SVector type
@test angle(SVector(1,0)) ≈ 0
@test angle(SVector(1,1)) ≈ π/4
@test angle(SVector(0,1)) ≈ π/2
@test angle(SVector(-1,-1)) ≈ -3π/4

# Type stability
@test typeof(angle([0.9, -1.5]))     == Float64
@test typeof(angle([0.9f0, -1.5f0])) == Float32

# Test errors
@test_throws Exception angle([1])
@test_throws Exception angle([1,2,3])
@test_throws Exception angle(SVector(1))
@test_throws Exception angle(SVector(3))

# Branch cut at π should be same as built-in function
@test all(     sqrt.(exp.(im.*range(-π,stop=π,length=21)), π)
           .== sqrt.(exp.(im.*range(-π,stop=π,length=21))))

# Hand picked examples where specific branches should be picked
@test sqrt.(-1.0+1.0im, π/2) == -sqrt(-1.0+1.0im)
@test sqrt.( 1.0+1.0im, π/2) ==  sqrt( 1.0+1.0im)

# Values along real line should always be the same
@test all(     sqrt.(range(0,stop=100,length=5), range(-π,stop=π,length=7)')
           .== sqrt.(range(0,stop=100,length=5)))

# Test floor function
@test floor(11, 10) ≈ 10
@test floor(15, 10) ≈ 10
@test floor(19, 10) ≈ 10

# Test round function
@test round(11, 10) ≈ 10
@test round(15, 10) ≈ 20
@test round(19, 10) ≈ 20

# Test ceil function
@test ceil(11, 10) ≈ 20
@test ceil(15, 10) ≈ 20
@test ceil(19, 10) ≈ 20


text:
module UsefulModes
export link_start!

using DataStructures
using ReplMaker
using REPL.LineEdit
# a bug of julia? workaround:
modes = Dict{Symbol, Any}()
fs = Function[]
register(f) = push!(fs, f)

include("astmode.jl")
include("selectormode.jl")

register() do
    enter_mode!(modes[:selector])
end

link_start!() = begin
    for each in fs
        each()
    end
    println("available modes:")
    for each in keys(modes)
        print(each, " ")
    end
end
end # module


text:
#=
Given an array of strictly the characters 'R', 'G', and 'B', segregate the values of the array so that all the Rs come first, the Gs come second, and the Bs come last. You can only swap elements of the array.

Do this in linear time and in-place.

For example, given the array ['G', 'B', 'R', 'R', 'B', 'R', 'G'], it should become ['R', 'R', 'R', 'G', 'G', 'B', 'B'].
=#

function segregate_rgb(s::Array{Char, 1})
    for char_index in eachindex(s)
        if s[char_index] == 'R'
            if !isnothing(findfirst(isequal('G'), s)) && findfirst(isequal('G'), s) < char_index
                s = swap(s, findfirst(isequal('G'), s), char_index)
            end
        elseif s[char_index] == 'G'
            if !isnothing(findlast(isequal('R'), s)) && findlast(isequal('R'), s) > char_index
                s = swap(s, findlast(isequal('R'), s), char_index)
            elseif !isnothing(findfirst(isequal('B'), s)) && findfirst(isequal('B'), s) < char_index
                    s = swap(s, findfirst(isequal('B'), s), char_index)
            end
        elseif s[char_index] == 'B'
            if !isnothing(findlast(isequal('G'), s)) && findlast(isequal('G'), s) > char_index
                s = swap(s, findlast(isequal('G'), s), char_index)
            end
        end
    end
    return s
end

function swap(s::Array{Char, 1}, x::Int, y::Int)
    out_string = copy(s)
    out_string[x] = s[y]
    out_string[y] = s[x]
    return out_string
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using StaticArrays
using Dates

abstract type AbstractSpaceTimePoint end

struct SpacePoint{S} <: AbstractSpaceTimePoint
    x::S
    SpacePoint(x::S) where {S <: SVector} = new{S}(x)
end
SpacePoint(x::AbstractVector) = SpacePoint(SVector(x...))

Base.size(sp::SpacePoint) = size(sp.x)
Base.getindex(sp::SpacePoint, i::Int) = getindex(sp.x, i)
Base.setindex!(sp::SpacePoint, v, i::Int) = setindex!(sp.x, v, i)
location(sp::SpacePoint) = sp.x
Base.show(io::IO, sp::SpacePoint) = print(io, "SpacePoint: ", location(sp).data)

struct SpaceTimePoint{S, T} <: AbstractSpaceTimePoint
    x::S
    t::T
    function SpaceTimePoint(x::S, t::T) where {S<:SVector, T<:Dates.AbstractTime}
        return new{S, T}(x, t)
    end
end
function SpaceTimePoint(x::AbstractVector, t::T) where {S<:SVector, T<:Dates.AbstractTime}
    return SpaceTimePoint(SVector(x...), t)
end

Base.size(stp::SpaceTimePoint) = size(stp.x) + 1
Base.getindex(stp::SpaceTimePoint, i::Int) = getindex(stp.x, i)
Base.setindex!(stp::SpaceTimePoint, v, i::Int) = setindex!(stp.x, v, i)
location(stp::SpaceTimePoint) = stp.x

function Base.show(io::IO, ::MIME"text/plain", stp::SpaceTimePoint)
    print(io, "SpaceTimePoint: ", location(stp).data, ", ", stp.t)
end
Base.show(io::IO, stp::SpaceTimePoint) = print(io, "SpaceTimePoint: ", location(stp).data,", ",  stp.t)

# function Base.show(io::IO, ::MIME"text/plain", sp::SpacePoint{S}) where {S}
#     print(sp)
# end

struct Echogram{T,N,S} <: AbstractArray{T,N}
    data::Array{T,N}
    z::AbstractVector{T}
    s::AbstractVector{S}

    function Echogram(data::Array{T, N}, z::AbstractVector{T},
            s::AbstractVector{S}) where {T, N, S<:AbstractSpaceTimePoint}
        sz = size(data)
        @assert length(z) == sz[1]
        if length(sz) > 1
            @assert length(s) == sz[2]
        else
            @assert length(s) == 1
        end
        return new{T, N, S}(data, z, s)
    end
end
const Ping{T,S} = Echogram{T,1,S}
Ping(data::Vector{T}, z::Vector{T}, s::Vector{S}) where {T, N, S<:AbstractSpaceTimePoint} =
    Echogram(data, z, s)

Base.size(eg::Echogram) = size(eg.data)
Base.getindex(eg::Echogram, i) = getindex(eg.data, i)
Base.getindex(eg::Echogram{T,N}, inds...) where {T,N} = getindex(eg.data, inds...)
# Base.getindex(eg::Echogram{T,N}, I::Vararg{Int, N}) where {T,N} = getindex(eg.data, I)
Base.setindex!(eg::Echogram, i) = getindex(eg.data, i)
Base.setindex!(eg::Echogram{T,N}, v, I::Vararg{Int, N}) where {T,N} = setindex!(eg.data, v, I)

depths(eg::Echogram) = eg.z
spacetime(eg::Echogram, i) = eg.s[i:i]
spacetime(eg::Echogram) = eg.s
spacetime(png::Ping) = png.s[1]
getping(eg::Echogram, i::Integer) = Echogram(eg[:, i], depths(eg), spacetime(eg, i))

function Base.show(io::IO, ::MIME"text/plain", eg::Echogram{T,N,S}) where {T,N,S}
    println(io, "Echogram{$T,$N}")
    println(io, "Depth: $(minimum(eg.z)):$(maximum(eg.z))")
    println(io, "Locations: $(first(eg.s))...$(last(eg.s))")
    show(io, MIME"text/plain"(), eg.data)
end

max_stars_repo_name: ElOceanografo/ScatteringLayers

max_issues_repo_name: ElOceanografo/ScatteringLayers

max_forks_repo_name: ElOceanografo/ScatteringLayers


text:
# both src or dst may be a Symbol that should be treated as Node
# only possible for function we own
src_dst_symbol_as_node = (:DirectedEdge, :UndirectedEdge)
for f in src_dst_symbol_as_node
    @eval $f(src::Symbol, dst::Symbol) = $f(Node(src), Node(dst))
end

# either dst or src Symbol should be treated as Node (both is possibly type piracy)
dst_symbol_as_node = (:*, :^, src_dst_symbol_as_node...)
for f in dst_symbol_as_node
    @eval $f(src, dst::Symbol) = $f(src, Node(dst))
end

# src may be a Symbol that should be treated as Node
src_symbol_as_node = (:ModifyingNode, dst_symbol_as_node...)
for f in src_symbol_as_node
    @eval $f(src::Symbol, dst) = $f(Node(src), dst)
end

import Base.convert
convert(::Type{T}, x) where {T <: AbstractNode} = T(x)
convert(::Type{T}, x::Symbol) where {T <: AbstractNode} = SimpleNode(x)
# conversion of Nodes does not change anything
convert(::Type{T}, x::T) where {T <: AbstractNode} = x

import Base.promote_rule
# Symbols get promoted
promote_rule(::Type{T}, ::Type{Symbol}) where {T <: AbstractNode} = AbstractNode
promote_rule(::Type{T}, ::Type{Symbol}) where {T <: Node} = Node


text:
import Polynomials
using Polynomials

function polymult(z::AbstractVector{T}, l::AbstractVector{Int}) where {T<:Number}
#
# construct monic polynomial from its roots
# use sorting and a special order of processing to increase accuray.
# special implementation of poly()
# returns Real polynomal, if all conjugate roots with equal multiplicities.
#
#   INPUT:
#            z vector of distinct roots of polynomial
#            l vector of corresponding multiplicities
#   OUTPUT:
#            ff Poly representing the product (x - z[k]) ^ l[k]
#
    # println("polymult($(z), $(l)")
    # z = PolyZeros(z)
    products(z, l)
end

# Calculate the coefficient operator of the pejorative manifold
# PolyRoots is used to leverage calculation of (x-z[j]) - products
#
function polyG(z::AbstractVector{<:Number}, mult::AbstractVector{<:Integer})
    polyG(products(z, mult))
end
function polyG(p::AbstractArray{<:Number})
    p[2:end]
end

# The derivative of G at point z
function polyDG(zz::AbstractVector{S}, ll) where S<:Number
    m = length(zz)
    s = products(zz, ll .- 1)
    n0 = length(s)
    n = sum(ll)
    Df = zeros(S, n, m)

    for j = 1:m
        Df[1:n0,j] .= s
        ix = n0
        for k = 1:m
            if k != j && ll[j] > 0
                ix = product1!(ix, view(Df,1:n,j), 1, -zz[k])
            end
        end
        Df[1:n,j] .*= -ll[j]
    end
    isreal(Df) ? real(Df) : Df
end

function polyDG(zz::AbstractVector{S}, ll::AbstractVector{<:Integer}, p::AbstractArray{T}, zs::AbstractVector{<:Integer}) where {S<:Number,T<:Number}

    w = [1/max(abs(pp), 1) for pp in p]
    m = length(zz)
    n = sum(ll)
    Df = zeros(T, n, m)
    for j = 1:m
        zsj = zs[j]
        zsj <= 0 || zsj > j || continue
        q, r = adiv(p, -zz[j], w)
        q .*= -ll[j]
        if zsj <= 0
            Df[:,j] .= q
        else
            if zsj > j
                Df[:,j] .= real(q)
                Df[:,zsj] .= imag(q)
            end
        end
    end
    Df
end

function realtocoeff(z::AbstractVector{<:Real}, S::Type{<:Complex}, zs::AbstractVector)

    r = similar(z, S)
    all( zs .<= 0 ) && return r
    for j = 1:length(z)
        zsj = zs[j]
        if zsj <= 0
            r[j] = z[j]
        elseif zsj > j
            r[j] = z[j]/2 - im*z[zsj]/2
            r[zsj] = conj(r[j])
        end
    end
    r
end
function realtocoeff(z::AbstractVector, ::Type, zs::AbstractVector)
    z
end

"""

Calculate polynomial product of `p`and `q` given as vectors.
"""
product

function product!(r::AbstractVector{U}, p::AbstractVector{S}, q::AbstractVector{T}) where {U<:Number,S<:Number,T<:Number}

  np = length(p)
  nq = length(q)
  nr = length(r)
  nr >= np + nq - 1 || throw(ArgumentError("size of result vector must be exactly (length(a) + length(b) - 1 but $np + $nq - 1 != $nr"))

  nq > np && return product(q, p)
  U == promote_type(S, T, U) || throw(ArgumentError("result type not compatible with input"))

  if nq == 0
    return S[]
  end
  E = one(U)
  copy!(r, p)
  fill!(view(r, np+1:nr), zero(U))
  q1 = U(q[1])
  if q1 != E
    for j = 1:np
      r[j] *= q1
    end
  end
  for k = 2:nq
    q1 = U(q[k])
    for j = 1:np
      jk = j + k - 1
      r[jk] = muladd(q1, p[j], r[jk])
    end
  end
  r
end

function product(p::AbstractVector{S}, q::AbstractVector{T}) where {S<:Number,T<:Number}
  U = promote_type(S,T)
  product!(zeros(U, length(p) + length(q) - 1), p, q)
end

"""
    products(z::Vector{T}, ll::Vector{Int})

Calculate coefficients of polynomial given by its roots with multiplicities.
The z should be ordered by increasing absolute value for better accuracy.
"""
function products(z::AbstractVector{T}, ll::AbstractVector{Int}) where T <:Number
    m = length(z)
    m == length(ll) || throw(ArgumentError("root vector and muliplicities have differnet length"))
    all(ll .>= 0) || throw(ArgumentError("multiplicities must not be negative"))
    n = sum(ll)
    S, = zerostructure(z, ll)
    r = zeros(S, n+1)
    r[1] = 1
    Z = zero(S)
    ix = 1
    jj = sortperm(z, by=abs2, rev=false)
    z = view(z, jj)
    ll =  ll[jj]
    while sum(ll) > 0
        p = minpos(ll)
        jm = findall(x-> x >= p, ll)
        s = singprod(z[jm], S)
        s = power(s, p)
        ix = convolute!(ix, r, s)
        ll[jm] .-= p
    end
    r
end

"""
    singprod(z::Vector, S::Type))

Return coefficients of polynomial `prod((x - zz) for zz in z)`

"""
function singprod(z::AbstractVector{T}, ::Type{S}) where {T<:Number,S<:Number}
    m = length(z)
    r = Vector{S}(undef, m+1)
    r[1] = 1
    ix = 1
    if S<:Real && T<:Complex
        for k = 1:m
            zk = z[k]
            rr = -real(zk)
            ri = imag(zk)
            if iszero(ri)
                ix = product1!(ix, r, 1, rr)
            elseif ri > 0
                ix = product2!(ix, r, 1, 2*rr, abs2(zk))
            end
        end
    else
        for k = 1:m
            ix = product1!(ix, r, 1, -z[k])
        end
    end
    r
end

"""
    minpos(v::Vector{<:Integer})

Return the minimal positive integer in v
"""
function minpos(x::AbstractVector{T}) where T<:Integer
    y = typemax(T)
    minimum(x .+ y) - y
end

"""
    convolute!(ix, r, p)

If `r` is a vector representing a polynomial in `r[1:ix]` and p is a vector
containing another polynomial, modify `r` to contain the product.
Assume `ix >= 0` and `length(r) >= ix + length(p) - 1`.
"""
function convolute!(ix::Integer, r::AbstractVector, p::AbstractVector)
    ix <= 0 && return ix
    n = length(p) - 1
    ri = r[ix]
    for j = 0:n
        r[ix+j] = ri * p[j+1]
    end
    for i = ix-1:-1:1
        ri = r[i]
        for j = n:-1:1
            r[i+j] += ri * p[j+1]
        end
        r[i] = ri * p[1]
    end
    ix + n
end

"""
    convolute(p, q)

Return the convolution (polynomial product) of vectors p and q.
The coefficient order does not matter. No checks for zero coefficients.
"""
function convolute(p::AbstractVector{T}, q::AbstractVector{S}) where {T,S}
    R = promote_type(T,S)
    m = length(p)
    n = length(q)
    (n == 0 || m == 0) && return R[]
    r = similar(p, R, n + m - 1)
    copyto!(r, p)
    convolute!(m, r, q)
    r
end

"""
    product1!(ix, r, ll, p1)

multiply by (x + p1)^ll. ix is the current used size of r.
"""
function product1!(ix::Int, r::AbstractVector{T}, ll::Int, zk::T) where T<:Number
    for j = 1:ll
        r[ix+1] = r[ix] * zk
        for i = ix:-1:2
            r[i] = muladd(r[i-1], zk, r[i])
        end
        ix += 1
    end
    ix
end

"""
    product2!(ix, r, ll, p1, p2)

multiply by (x^2 + p1*x + p2)^ll. ix is the current used size of r.
"""
function product2!(ix::Int, r::AbstractVector{T}, ll::Int, p1::T, p2::T) where T<:Number
    for j = 1:ll
        r[ix+1] = r[ix+2] = 0
        for i = ix:-1:1
            rr = muladd(r[i+1], p1, r[i+2])
            r[i+2] = muladd(r[i], p2, rr)
        end
        r[2] = muladd(r[1], p1, r[2])
        ix += 2
    end
    ix
end

"""
    adiv(c::Vector, a::Number, weight::Vector) -> b::Vector

Divide Polynomial `x^n * c[1] + ... + c[n+1]` by `x + a` giving
`x^(n-1) * b[1] + ... + b[n]` under the condition that the division remainder is zero.

`weight` is a positive weight vector with the same size as `c` used to solve the
least squares problem `minimize norm(w .* (A * b - c)) where `A` is the `n+1 x n`
Cauchy matrix describing the multiplication with `x + a`.
The results tend to be much more accurate than that of the usual long division,
especially in the case of hughe coefficients, which are damped by `weight = 1 ./ abs.(c)`.
"""
function adiv(c::AbstractVector{<:Number}, a::S, w::AbstractVector{<:Real}) where S<:Number
    T = real(S)
    n = length(c) - 1
    aa2 = abs2(a)
    aa = sqrt(aa2)
    ac = conj(a)
    r = Vector{T}(undef, n)
    d = Vector{T}(undef, n)
    b = Vector{S}(undef, n)
    b[1] = 1
    b[n] = c[n+1] / a
    w1 = w[2]
    u1 = (c[2] - a) * w1
    for i = 2:n-2
        w2 = w[i+1]
        u2 = c[i+1] * w2
        z = d[i] = hypot(w1, w2*aa)
        w1z = w1 / z
        w2z = w2 / z
        r[i] = w2z * w2
        b[i] = w1z * u1 + w2z * u2 * ac
        u1 = -w2z * u1 * a + w1z * u2
        w1 *= w2z
    end
        w2 = w[n]
        u2 = (c[n] - b[n]) * w2
        z = d[n-1] = hypot(w1, w2*aa)
        w1z = w1 / z
        w2z = w2 / z
        r[n-1] = w2z * w2
        b[n-1] = w1z * u1 + w2z * u2 * ac
        u1 = -w2z * u1 * a + w1z * u2
    bi = b[n-1] = b[n-1] / d[n-1]
    for i = n-2:-1:2
        bi = b[i] = (b[i] - bi * r[i] * ac) / d[i]
    end
    b, u1
end

function power(p::AbstractArray, n::Integer)
    (Poly(p)^n).a
end

function Polynomials.roots(pc::AbstractVector{T}) where T<:Number
# redirect to implementation in Polynomials
	p = Poly(reverse(pc))
	roots(p)
end

function deriv(p::AbstractVector{T}) where T<:Number

    n = length(p) - 1
    map(k -> (p[k] * (n-k+1)), 1:n)
end